计算和使用 Jacobian for solve_ivp 方法答案

【问题标题】：Computing and using Jacobian for solve_ivp method计算和使用 Jacobian for solve_ivp 方法
【发布时间】：2019-08-11 18:53:34
【问题描述】：

我正在尝试计算大型一阶 ODE 系统的雅可比行列式。我需要在我的 solve_ivp 代码中使用 Jacobian 矩阵。

这是我尝试计算雅可比的方法。

C1 = sp.symbols('C[0], C[1], C[2], C[3], C[4], C[5], C[6], C[7], C[8], C[9], C[10], C[11], C[12], C[13], C[14], C[15], C[16], C[17]')

X = sp.Matrix([[-1.0*final_kcm[0]*C1[0] - final_kcm[1]*C1[0]*C1[11] - final_kcm[2]*C1[0]*C1[12] - final_kcm[3]*C1[0]*C1[14] - final_kcm[4]*C1[0]*C1[15] - final_kcm[5]*C1[0]*C1[16],
       final_kcm[2]*C1[0]*C1[12] - final_kcm[8]*C1[1]*C1[11] + final_kcm[13]*C1[17]*C1[12],
       final_kcm[1]*C1[0]*C1[11] + final_kcm[6]*C1[12]*C1[11] + final_kcm[7]*C1[11]*C1[13] + final_kcm[8]*C1[11]*C1[1] + final_kcm[9]*C1[11]*C1[9] + final_kcm[10]*C1[11]*C1[10] +  final_kcm[12]*C1[11]*C1[17] + 2*final_kcm[20]*C1[11]*C1[8],
       2*final_kcm[20]*C1[12]*C1[8],
       final_kcm[15]*C1[15]*C1[17],
       final_kcm[7]*C1[11]*C1[13] - final_kcm[18]*C1[5]*C1[11],
       final_kcm[14]*C1[17]*C1[13],
       final_kcm[19]*C1[15]*C1[8],
       final_kcm[17]*C1[15] - 2*final_kcm[19]*C1[12]*C1[8] - final_kcm[20]*C1[12]*C1[8],
       final_kcm[3]*C1[0]*C1[14] - final_kcm[9]*C1[11]*C1[9],
       final_kcm[5]*C1[0]*C1[16] - final_kcm[10]*C1[11]*C1[10],
       final_kcm[0]*C1[0] - final_kcm[1]*C1[0]*C1[11] - final_kcm[6]*C1[12]*C1[11] - final_kcm[7]*C1[11]*C1[13] - final_kcm[8]*C1[11]*C1[1] - final_kcm[9]*C1[11]*C1[9] - final_kcm[10]*C1[11]*C1[10] - final_kcm[11]*C1[11]*C1[10] - final_kcm[12]*C1[11]*C1[17] + final_kcm[14]*C1[17]*C1[14] + final_kcm[15]*C1[15]*C1[17] + final_kcm[16]*C1[13] - final_kcm[18]*C1[5]*C1[11] + final_kcm[19]*C1[12]*C1[8] - 2*final_kcm[20]*C1[12]*C1[8],
       final_kcm[0]*C1[0] - final_kcm[2]*C1[0]*C1[12] - final_kcm[6]*C1[12]*C1[11] + final_kcm[8]*C1[11]*C1[1] - final_kcm[13]*C1[17]*C1[12],
       final_kcm[1]*C1[0]*C1[11] + final_kcm[2]*C1[0]*C1[12] + final_kcm[3]*C1[0]*C1[14] + final_kcm[4]*C1[0]*C1[15] + final_kcm[5]*C1[0]*C1[16] - final_kcm[7]*C1[11]*C1[13] - final_kcm[16]*C1[13],
       -1.0*final_kcm[3]*C1[0]*C1[14] + final_kcm[9]*C1[11]*C1[9] + final_kcm[11]*C1[11]*C1[10] - final_kcm[14]*C1[17]*C1[14],
       -1.0*final_kcm[4]*C1[0]*C1[15] + final_kcm[12]*C1[11]*C1[17] + final_kcm[13]*C1[17]*C1[12] - final_kcm[15]*C1[15]*C1[17] - final_kcm[17]*C1[15] - final_kcm[19]*C1[15]*C1[8],
       -1.0*final_kcm[5]*C1[0]*C1[16] + final_kcm[10]*C1[11]*C1[10] + final_kcm[18]*C1[5]*C1[11],
       final_kcm[4]*C1[0]*C1[15] + final_kcm[6]*C1[12]*C1[11] - final_kcm[11]*C1[11]*C1[10] - final_kcm[12]*C1[11]*C1[17] - final_kcm[13]*C1[17]*C1[12] - final_kcm[14]*C1[17]*C1[13] - final_kcm[15]*C1[15]*C1[17] + final_kcm[16]*C1[13]]])
Y = sp.Matrix(C1)
J =  X.jacobian(Y)

这会返回一个矩阵元素，如Matrix([[..]....])。

如何将其转换为我需要在solve_ivp 问题中使用它的内容？

这是我的solve_ivp 问题的代码。

Temperature = 500.0

Temp_K = Temperature + 273.15
r = 8.314459848 #J/K*mol
R = r/1000.0 #kJ/K*mol

Ea = [342,7,42,45,34,48,13,12,4,6,15,0,56,31,30,61,84,90,70,20,70]

k_0 =[5.9 * (10**15), 1.3 * (10**13), 1.0 * (10**12), 5.0 * (10**11), 1.2 * (10**13), 2.0 * (10**11), 1.0 * (10**13), 1.0 * (10**13), 1.7 * (10**13), 1.2 * (10**13), 1.7 * (10**13), 9.1 * (10**10), 1.2 * (10**14), 5.0 * (10**11), 2.0 * (10**10), 3.0 * (10**11), 2.1 * (10**14), 5.0 * (10**14), 2.0 * (10**13), 1.0 * (10**14), 1.6 * (10**14)]



fk1 = [mp.exp((-1.0 * x)/(R*Temp_K)) for x in Ea]
final_k = [fk1*k_0 for fk1,k_0 in zip(fk1,k_0)]
final_kcm = [mpf(x) for x in final_k]

def rhs(t, C):
return[(-1.0*final_kcm[0]*C[0]) - final_kcm[1]*C[0]*C[11] - final_kcm[2]*C[0]*C[12] - final_kcm[3]*C[0]*C[14] - final_kcm[4]*C[0]*C[15] - final_kcm[5]*C[0]*C[16],
       final_kcm[2]*C[0]*C[12] - final_kcm[8]*C[1]*C[11] + final_kcm[13]*C[17]*C[12],
       final_kcm[1]*C[0]*C[11] + final_kcm[6]*C[12]*C[11] + final_kcm[7]*C[11]*C[13] + final_kcm[8]*C[11]*C[1] + final_kcm[9]*C[11]*C[9] + final_kcm[10]*C[11]*C[10] +  final_kcm[12]*C[11]*C[17] + 2*final_kcm[20]*C[11]*C[8],
       2*final_kcm[20]*C[12]*C[8],
       final_kcm[15]*C[15]*C[17],
       final_kcm[7]*C[11]*C[13] - final_kcm[18]*C[5]*C[11],
       final_kcm[14]*C[17]*C[13],
       final_kcm[19]*C[15]*C[8],
       final_kcm[17]*C[15] - 2*final_kcm[19]*C[12]*C[8] - final_kcm[20]*C[12]*C[8],
       final_kcm[3]*C[0]*C[14] - final_kcm[9]*C[11]*C[9],
       final_kcm[5]*C[0]*C[16] - final_kcm[10]*C[11]*C[10],
       final_kcm[0]*C[0] - final_kcm[1]*C[0]*C[11] - final_kcm[6]*C[12]*C[11] - final_kcm[7]*C[11]*C[13] - final_kcm[8]*C[11]*C[1] - final_kcm[9]*C[11]*C[9] - final_kcm[10]*C[11]*C[10] - final_kcm[11]*C[11]*C[10] - final_kcm[12]*C[11]*C[17] + final_kcm[14]*C[17]*C[14] + final_kcm[15]*C[15]*C[17] + final_kcm[16]*C[13] - final_kcm[18]*C[5]*C[11] + final_kcm[19]*C[12]*C[8] - 2*final_kcm[20]*C[12]*C[8],
       final_kcm[0]*C[0] - final_kcm[2]*C[0]*C[12] - final_kcm[6]*C[12]*C[11] + final_kcm[8]*C[11]*C[1] - final_kcm[13]*C[17]*C[12],
       final_kcm[1]*C[0]*C[11] + final_kcm[2]*C[0]*C[12] + final_kcm[3]*C[0]*C[14] + final_kcm[4]*C[0]*C[15] + final_kcm[5]*C[0]*C[16] - final_kcm[7]*C[11]*C[13] - final_kcm[16]*C[13],
       (-1.0*final_kcm[3]*C[0]*C[14]) + final_kcm[9]*C[11]*C[9] + final_kcm[11]*C[11]*C[10] - final_kcm[14]*C[17]*C[14],
       (-1.0*final_kcm[4]*C[0]*C[15]) + final_kcm[12]*C[11]*C[17] + final_kcm[13]*C[17]*C[12] - final_kcm[15]*C[15]*C[17] - final_kcm[17]*C[15] - final_kcm[19]*C[15]*C[8],
       (-1.0*final_kcm[5]*C[0]*C[16]) + final_kcm[10]*C[11]*C[10] + final_kcm[18]*C[5]*C[11],
       final_kcm[4]*C[0]*C[15] + final_kcm[6]*C[12]*C[11] - final_kcm[11]*C[11]*C[10] - final_kcm[12]*C[11]*C[17] - final_kcm[13]*C[17]*C[12] - final_kcm[14]*C[17]*C[13] - final_kcm[15]*C[15]*C[17] + final_kcm[16]*C[13]]
res = solve_ivp(rhs, (0, 3), [0.9999999999929631, 1.8518729652310317e-12, 3.7964305690400015e-12, 4.3194364603720453e-32, 8.840289274135903e-29, 7.58855704948555e-17, 3.747886894978005e-21, 7.443392243056976e-25, 2.713933190549257e-19, 1.3888890038089627e-12, 3.778455148948492e-26, 1.3426298908781515e-12, 4.630881200895968e-13, 2.823935453386071e-12, 9.7200025431433e-13, 3.024125323625077e-19, 1.7428556044758935e-25, 2.314889117721353e-12],  method = 'Radau', jac = J)

当我运行代码时，结果出现以下错误。

TypeError: __array__() takes 1 positional argument but 2 were given

我该怎么办？

【问题讨论】：

由于您无论如何都在使用符号，您可能想看看this module of mine，它使用符号输入并且可以使用符号推导完全自动地为雅可比生成（精确）函数。解决 IVP 是可能的后端之一。它还可以为您提供相关的速度提升。
非常感谢，我会继续努力，希望它能奏效。

标签： python python-3.x matrix scipy ode

【解决方案1】：

编辑：我想我使用lambdify解决了它，如下所示：

C1 = sp.symbols('C[0], C[1], C[2], C[3], C[4], C[5], C[6], C[7], C[8], C[9], C[10], C[11], C[12], C[13], C[14], C[15], C[16], C[17]') 
t1 = sp.symbols('t1') 
ydot = rhs(t1, C1) 
J = sp.Matrix(ydot).jacobian(C1) 
J_func = sp.lambdify((t1,C1), J)

然后可以在solve_ivp方法中使用J_func，如下所示。

Y0 = [0.9999999999929631, 1.8518729652310317e-12, 3.7964305690400015e-12, 4.3194364603720453e-32, 8.840289274135903e-29, 7.58855704948555e-17, 3.747886894978005e-21, 7.443392243056976e-25, 2.713933190549257e-19, 1.3888890038089627e-12, 3.778455148948492e-26, 1.3426298908781515e-12, 4.630881200895968e-13, 2.823935453386071e-12, 9.7200025431433e-13, 3.024125323625077e-19, 1.7428556044758935e-25, 2.314889117721353e-12]

res = solve_ivp(rhs, (0, 30), Y0 , method = 'Radau', jac = J_func)

编辑 2：这是生成相同内容的更简洁的代码。

def jacob(x, C):
  ydot = rhs(x,C)
  J =  sp.Matrix(ydot).jacobian(C)
  J_func = sp.lambdify((x,C) , J)
  return J_func

【讨论】：

您似乎在使用有限差分来逼近雅可比行列式。这是不必要的，因为当您不提供雅可比行列时，Solve IVP 会自行执行此操作。
想详细了解如何计算雅可比？我似乎真的被卡住了，似乎没有办法计算它。我在数学交换中找到了这种方法并且它有效，所以我只是假设它是正确的。编辑：我在问题中使用的雅可比方法似乎也有效，但我无法让它与solve_ivp一起使用。你知道这是为什么吗？
你的方法不正确（AFAICT），只是画蛇添足。 Jacobian 是基于您右手边的偏导数的clearly defined thing。你可以在纸上计算它（考虑到你的问题的复杂性，我不建议这样做），用计算机象征性地计算它（这很有意义，因为你无论如何都在使用符号），或者使用finite differences 进行数字估计。你在这里做的是后者，但这也是 Solve IVP 自己做的事情。
好的，谢谢，当我看到 epsilon 值时，我有一个小小的倾向，但我并不完全确定任何事情。当我使用'J = X.jacobian（Y）'时，它似乎象征性地计算它，但我无法弄清楚如何让它工作。当我重新开始工作时，我会下载你的模块并发布结果，如果我能让它工作的话。感谢您的帮助！
编辑：似乎只是对矩阵进行了简单化处理，才允许使用它。我仍然不是 100% 确定一切都是正确的，但这似乎更合理。我在这个网站上找到了一个指南link