Python中使用sympy计算雅可比矩阵后调用numpy求逆报错问题

错误原因 #

• Sympy生成的雅可比矩阵是符号类型，Numpy的inv方法不支持该类型输入，需要提前转换为数值类型的Numpy数组
• 代码存在多处笔误：
  • 幂运算逻辑错误：3.10**7、3.10^7为书写错误，其中^是Python的异或运算符，不是幂运算，对应3乘10的7次方应写为3e7
  • 定义雅可比的表达式字符串被强制换行，触发语法错误
  • 后续用来存储逆矩阵的变量名j和雅可比计算函数重名，会导致后续调用异常
  • Numpy导入别名py不符合通用习惯，容易和其他库混淆，建议改为np

修正后的可运行代码 #

import numpy as np
from numpy.linalg import inv
from sympy import Matrix, symbols
import warnings

def f1(y1, y2, y3, y1_old, dt):
    return y1_old + (-0.04*y1 + 1e4*y2*y3)*dt

def f2(y1, y2, y3, y2_old, dt):
    return y2_old + (0.04*y1 - 1e4*y2*y3 - 3e7*(y2**2))*dt

def f3(y1, y2, y3, y3_old, dt):
    return y3_old + (3e7*(y2**2))*dt

def calc_jac(y1,y2,y3):
    y1_sym, y2_sym, y3_sym = symbols('y1 y2 y3')
    expr_list = [
        -0.04*y1_sym + 1e4*y2_sym*y3_sym,
        0.04*y1_sym - 1e4*y2_sym*y3_sym - 3e7*(y2_sym**2),
        3e7*(y2_sym**2)
    ]
    # 计算符号雅可比矩阵
    jac_sym = Matrix(expr_list).jacobian([y1_sym, y2_sym, y3_sym])
    # 代入数值转换为numpy浮点数组
    jac_np = np.array(jac_sym.subs({
        y1_sym: y1,
        y2_sym: y2,
        y3_sym: y3
    }), dtype=np.float64)
    return jac_np

warnings.filterwarnings("ignore", category=DeprecationWarning)

y_old = np.zeros((3,1))
y_old[0] = 1

# 隐式变量初始猜测值
y_guess = 2*np.ones((3,1))

# 初始新值设定
y_new = np.ones((3,1))

F = np.copy(y_new)

start_time = 0
end_time = 10
dt = 0.01
nt = np.arange(start_time,end_time,dt)

error = 9e9
tol = 1e-10
alpha = 0.8

jac = calc_jac(y_guess[0],y_guess[1],y_guess[2])

F[0] = f1(y_guess[0],y_guess[1],y_guess[2], y_old[0], dt)
F[1] = f2(y_guess[0],y_guess[1],y_guess[2], y_old[1], dt)
F[2] = f3(y_guess[0],y_guess[1],y_guess[2], y_old[2], dt)

jac_inv = inv(jac)

额外优化建议 #

如果不需要保留符号计算逻辑，完全可以直接写数值形式的雅可比矩阵，不用调用Sympy，运行效率会高很多：

def calc_jac_numeric(y1, y2, y3):
    return np.array([
        [-0.04, 1e4*y3, 1e4*y2],
        [0.04, -1e4*y3 - 6e7*y2, -1e4*y2],
        [0, 6e7*y2, 0]
    ], dtype=np.float64)

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