python - 尝试同时求解微分方程

Question

我正在尝试构建化学 react 器设计代码，该代码能够解决 react 器的压降、转化率和温度问题。所有这些参数都有微分方程，因此我尝试在函数内定义它们，以便能够使用 ODEINT 对它们进行积分。然而，我构建的函数似乎有一个错误，我无法弄清楚是什么阻碍了我集成它。

我遇到的错误:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-32-63c706e84be5> in <module>
     1 X0=[0.05,1200,2]
----> 2 y=func(X0,0)

<ipython-input-27-6cfd4fef5ee2> in func(x, W)
     8     kp=np.exp(((42311)/(R*T))-11.24)
     9     deltah=-42471-1.563*(T-1260)+0.00136*(T**2 -1260**2)- 2,459*10e-7*(T**3-1260**3)
---> 10     ra=k*np.sqrt((1-X)/X)*((0.2-0.11*X)/(1-0.055*X)*(P/P0)-(x/(kp*(1-x)))**2)
    11     summ = 57.23+0.014*T-1.94*10e-6*T**2
    12     dcp=-1.5625+2.72*10e-3*T-7.38*10e-7*T**2

TypeError: unsupported operand type(s) for -: 'int' and 'list'

这是完整的代码:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

fi = 0.45
gas_density = 0.054  #density
Pres0 = 2  #pressure
visc = 0.09
U = 10
Ac = 0.0422
T0 = 1400  #and 1200 also
gc = 4.17 * 10**8
bed_density = 33.8
Ta = 1264.6  #wall temp
fa0 = 0.188

def func(x,W):
    X = x[0]
    T = x[1]
    P = x[2]
    P0 = 2
    R = 0.7302413
    k = np.exp((-176008 / T) - (110.1 * np.log(T) + 912.8))
    kp = np.exp(((42311) / (R * T)) - 11.24)
    deltah = -42471 - 1.563 * (T - 1260) + 0.00136 * (T**2 - 1260**2) - 2,459 * 10e-7 * (T**3 - 1260**3)
    ra = k * np.sqrt((1 - X) / X) * ((0.2 - 0.11 * X) / (1 - 0.055 * X) * (P / P0) - (x / (kp * (1 - x)))**2)
    summ = 57.23 + 0.014 * T - 1.94 * 10e-6 * T**2
    dcp = -1.5625 + 2.72 * 10e-3 * T - 7.38 * 10e-7 * T**2
    dxdw = 5.31 * k * np.sqrt((1 - X) / X) * ((0.2 - 0.11 * X) / (1 - 0.055 * X) * (P / P0) - (x / (kp * (1 - x)))**2)
    dpdw = (((-1.12 * 10**-8) * (1 - 0.55 * X) * T) / P) * (5500 * visc + 2288)
    dtdw = (5.11 * (Ta - T) + (-ra) * deltah) / (fa0 * (summ + x * dcp))
    return [dxdw, dpdw, dtdw]

X0 = [0.05, 1200, 2]
y = func(X0, 0)

提前致谢

Answer 1

内线

ra=k*np.sqrt((1-X)/X)*((0.2-0.11*X)/(1-0.055*X)*(P/P0)-(x/(kp*(1-x)))**2)

您可能想使用 ...X/(kp*(1-X))... 而不是 ...x/(kp*(1-x) )...(即使用上面的X)，下面的x是列表类型。

如果您想在某处使用某些列表变量l作为多个值，请将其转换为numpy数组la = np.array(l)并使用la 在 numpy 矢量化表达式中。