Mismatch in result for matrix exponential calculation using Python and Mathematica

Question

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}
 }

Answer 1

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