在python中用ODEINT求解微分方程

Question

我有一个耦合微分方程组，我已经在 Excel 中使用 Euler 求解了该方程组。现在我想在 python 中使用 ODE 求解器使其更精确。 但是，我的代码中一定有一个错误，因为曲线看起来与 Excel 中的不同。我不希望曲线最终达到 1 和 0。

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

# define reactor
def reactor(x,z):
    n_a = x[0]
    n_b = x[1]
    n_c = x[2]

    dn_adz = A * (-1) * B * (n_a/(n_a + n_b + n_c)) / (1 + C * (n_c/(n_a + n_b + n_c)))
    dn_bdz = A * (1) * B * (n_a/(n_a + n_b + n_c)) / (1 + C * (n_c/(n_a + n_b + n_c)))
    dn_cdz = A * (1) * B * (n_a/(n_a + n_b + n_c)) / (1 + C * (n_c/(n_a + n_b + n_c)))
    dxdz = [dn_adz,dn_bdz,dn_cdz]
    return dxdz

# initial conditions
n_a0 = 0.5775
n_b0 = 0.0
n_c0 = 0.0
x0 = [n_a0, n_b0, n_c0]

# parameters
A = 0.12
B = 3.1e-9
C = 4.02e15

# number of steps
n = 100

# z step interval (m)
z = np.linspace(0,0.0274,n)

# solve ODEs
x = odeint(reactor,x0,z)


# Plot the results
plt.plot(z,x[:,0],'b-')
plt.plot(z,x[:,1],'r--')
plt.plot(z,x[:,2],'k:')
plt.show()

初始条件保持不变且不随步骤变化是否存在问题？ 是不是应该像 Excel with Euler 一样，下一步使用珍贵步骤的条件/值？

Answer 1

从右侧的结构中，您可以得到状态变量n_a+n_b=n_a0+n_b0 和n_a+n_c=n_a0+n_c0 的常量组合。这意味着动态减少到n_a 的一维动态。

根据第一个方程，n_a 的导数对正的n_a 是负的，因此解正朝着n_a=0 下降。根据动力学常数，n_b 收敛到 n_a0+n_b0 和 n_c 收敛到 n_a0+n_c0。

尚不清楚如何在某些组件中收敛到 1，因为初始条件不支持这一点。除此之外，所描述的odeint 结果符合这种定性行为。