I wrote the simple python code for evaluating matrix exponential using truncated Taylor series given by the formula, where M is the matrix, t is time and k is the order of truncation, which isn't converging to correct values even for large orders. I wrote the same code in Mathematica, which is converging correctly. I am not able to find the error in the python code due to which it is producing incorrect result.
The Python code is-
def matrix_exponential_taylor(m, t, order):
result = np.eye(m.shape[0])
for i in range(1,order+1):
term=np.linalg.matrix_power(m*t,i)/np.math.factorial(i)
result+=term.astype(np.float64)
return result
A = np.array([[2, 1, 0, 1],
[1, 1, 1, 0],
[0, 1, 3, 1],
[1, 0, -1, 2]])
t = 3
order = 1000
result = matrix_exponential_taylor(A, t, order)
print(result)
output from Python code-
[[3594.20250036 1511.17187957 -1249.84655916 3085.96687505]
[5087.48182172 2074.0643393 769.10695263 4399.3840113 ]
[7966.0379292 4335.76396246 1637.40121583 7160.40223197]
[ -480.68704288 -1249.84655916 -3584.7820205 -1947.37752023]]
Mathematica Code-
MatrixExponentialTaylor[m_, t_, order_] :=
Module[{result, term}, result = IdentityMatrix[Length[m]];
For[i = 1, i < order, i++, term = MatrixPower[m*t, i]/Factorial[i];
result += term;];
result]
m = {{2, 1, 0, 1}, {1, 1, 1, 0}, {0, 1, 3, 1}, {1, 0, -1, 2}}
t = 3;
order = 1000;
result = N[MatrixExponentialTaylor[m, t, order]]
result // MatrixForm
output from Mathematica code-
{
{6110.33, 1480.1, -4384.45, 3441.73},
{10249., 4630.23, -942.716, 8307.71},
{17076.6, 7826.18, -1178.47, 14172.3},
{-5327.16, -4384.45, -5403.36, -6581.83}
}
You'll get faster and closer results by updating your term at every iteration rather than recomputing it from 0:
def matrix_exponential_taylor(m, t, order):
term = np.eye(m.shape[0])
result = np.eye(m.shape[0])
mt = m*t
for i in range(1,1+order):
term = term @ mt / i
result += term
return result