计算和使用雅可比行列式的solve_ivp方法

问题描述 投票：0回答：1

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

这是我如何尝试计算雅可比。

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([[..]....])这样的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

我该怎么办？

python python-3.x matrix scipy ode
1个回答
0
投票

编辑：我想我用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
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.