Вычисление и использование якобиана для метода solve_ivp - PullRequest
Новая
       4

Вычисление и использование якобиана для метода solve_ivp

0 голосов
/ 21 марта 2019

Я пытаюсь вычислить якобиан большой системы ОДУ первого порядка. Мне нужен якобиан для использования в моем коде 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

Что мне делать?

1 Ответ

0 голосов
/ 03 апреля 2019

РЕДАКТИРОВАТЬ: Я думаю, что я решил это с помощью 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)

J_func может затем использоваться в методе solve_ivp следующим образом. 

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