Python与Mathematica矩阵指数泰勒展开计算结果不一致问题

问题背景 #

用截断泰勒级数实现矩阵指数计算，其中M为输入矩阵、t为时间参数、k为截断阶数。Python版本即使使用高阶截断（如order=1000），结果仍无法收敛至正确值，但相同逻辑的Mathematica代码收敛正常，需定位Python代码错误。

Python代码及输出 #

import numpy as np

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)

输出结果：

[[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代码及输出 #

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

输出结果：

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

错误分析与修复 #

核心问题1：高次矩阵幂计算的精度与效率问题 #

Python原代码每次循环直接计算np.linalg.matrix_power(m*t, i)，当i很大时，矩阵高次幂的计算会引入大量数值误差，同时效率极低。正确的泰勒级数递推应利用前一项结果：
$$ term_i = term_{i-1} \times (M \times t) / i $$
这种方式避免了重复计算高次幂，大幅提升精度与效率。

核心问题2：阶乘的浮点溢出问题 #

np.math.factorial(i)返回整数，当i>170时，170!的数值超过float64的最大可表示范围（~1e308），转成float64会变成inf，导致后续项全部为0，累加结果提前停止收敛。而Mathematica的符号计算会自动处理大数精度，不会出现此类溢出。

修复后的Python代码 #

import numpy as np

def matrix_exponential_taylor(m, t, order):
    n = m.shape[0]
    result = np.eye(n, dtype=np.float64)
    mt = m * t  # 提前计算M*t，避免重复计算
    term = np.eye(n, dtype=np.float64)  # 初始项对应i=0，i从1开始递推
    for i in range(1, order+1):
        term = term @ mt / i  # 递推计算当前项
        result += term
    return result

A = np.array([[2, 1, 0, 1],
               [1, 1, 1, 0],
               [0, 1, 3, 1],
               [1, 0, -1, 2]], dtype=np.float64)
t = 3
order = 1000
result = matrix_exponential_taylor(A, t, order)
print(result)

验证结果 #

修复后的代码输出与Mathematica结果一致：

[[ 6110.33206138  1480.09927727 -4384.44814731  3441.72876233]
 [10249.00381732  4630.22870622  -942.71606734  8307.70977726]
 [17076.59720829  7826.17943347 -1178.46821826 14172.29819923]
 [-5327.15948497 -4384.44814731 -5403.35830832 -6581.8297493 ]]

额外说明 #

Mathematica的循环条件是i < order（循环到order-1），而Python原代码是range(1, order+1)（循环到order），但这个差异在order=1000时对结果影响可忽略，核心问题还是递推方式与浮点溢出。

内容的提问来源于stack exchange，提问作者XZark