解决Python中的一阶和二阶微分方程组

Question

我需要解决以下微分方程组：

$\frac{dx_1}{dt} = -k_1x_1+k_2x_2-(K_R)x_1y_1$
$\frac{dx_2}{dt} =  k_1x_1-k_2x_2-k_3x_2-(K_R)x_2y_2$
$\frac{dx_3}{dt} =  k_3x_3$
$\frac{dy_1}{dt} = -k_1y_1+k_2y_2-(K_R)x_1y_1$
$\frac{dy_2}{dt} =  k_1y_1-k_2y_2-k_3y_2-(K_R)x_2y_2$
$\frac{dy_3}{dt} =  k_3y_3$
$\frac{dz_1}{dt} = -k_1z_1+k_2z_2+(K_R)x_1y_1$
$\frac{dz_2}{dt} =  k_1z_1-k_2z_2-k_3z_2+(K_R)x_2y_2$
$\frac{dz_3}{dt} =  k_3z_3$

t = 0时的初始条件是x2 =1。而在t = 1时，将化合物y引入y2隔室中，y2 =10。KR的值是1e-3。

---

我已经解决了使用矩阵求幂的简单得多的系统，并且想知道是否可以使用类似的方法来解决上述系统。

我有一个隔离模型系统X，其简化版本如下：

然后是微分方程组：

我可以使用以下矩阵方法来求解此方程组。

首先，我编写速率矩阵[R]。从[R]可以得到一个新的矩阵[A]，方法是先用每个行元素之和的负数替换[R]的每个对角元素，然后对其进行转置：

我可以通过执行以下操作来计算每个隔室中的量：

在python中：

RMatrix = model_matrix.as_matrix()
row, col = np.diag_indices_from(RMatrix)
RMatrix[row, col] = -(RMatrix.sum(axis=1)-RMatrix[row,col])
AMatrix = RMatrix.T

def content(t):
    cont = np.dot(linalg.expm(t*AMatrix), x0))

此方法对我来说效果很好。

---

[上面的模型（原始问题）比系统X复杂一些。在该模型中，系统X和Y的1和2隔室中的反应物结合在一起，得到系统Z中的产物。

X + Y-> Z，反应常数为KR。

，相应的微分方程组将是：

我正在努力解决这种微分方程组（一阶和二阶）的方法，以便在给定初始条件KR和传输速率k1，k2，k3的情况下，在特定时间t计算每个隔室中的量等等...

对于一阶微分方程组，我是否可以使用上述矩阵方法求解？我在Python中还有哪些其他选择？

提前感谢！

Answer 1

嗯，正如评论中指出的那样，您的（更复杂的）ODE是非线性的。因此，矩阵指数方法将不再起作用。

通常，有两种解决ODE的一般方法。首先，您可以尝试找到一种符号解决方案。在大多数情况下，您会根据有根据的猜测采取某种方法。有几种类型的ODE已知符号解决方案。

但是，绝大多数ODE并非如此。因此，我们通常会提出一个数值解决方案，本质上是基于右侧对ODE进行数值积分。

结果不是显式函数，而是某个点上函数值的近似值。在python中，您可以使用scipy这样解决ODE。基于您的右手边（除非出现我的任何错误），这看起来像这样：

import numpy as np

import scipy.integrate 

k_1 = 1
k_2 = 1
k_3 = 1
K_R = 1

def eval_f(v, t):
    [x, y, z] = np.split(v, [3, 6])

    return np.array([-k_1*x[0] +k_2*x[1] - (K_R)*x[0]*y[0],
                     k_1*x[0] - k_2*x[1] - k_3*x[1] - (K_R)*x[1]*y[1],
                     k_3*x[2],
                     - k_1*y[0] + k_2*y[1] - (K_R)*x[0]*y[0],
                     k_1*y[0] - k_2*y[1] - k_3*y[1] - (K_R)*x[1]*y[1],
                     k_3*y[2],
                     - k_1*z[0] + k_2*z[1] + (K_R)*x[0]*y[0],
                     k_1*z[0] - k_2*z[1] - k_3*z[1] + (K_R)*x[1]*y[1],
                     k_3*z[2]])

initial = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])

t = np.linspace(0, 1, 10)

values = scipy.integrate.odeint(eval_f, initial, t)

# variable x[0]
print(values[:,0])

这将产生x1的以下值：

[1.         0.70643591 0.49587121 0.35045691 0.25034256 0.1809533
 0.13237994 0.09800056 0.07338967 0.05557138]

基于网格点

[0.         0.11111111 0.22222222 0.33333333 0.44444444 0.55555556
 0.66666667 0.77777778 0.88888889 1.        ]

[如果您想查看函数的行为，积分器可能就足够了。否则，我建议您在教科书中阅读有关ODE的符号方法...