求解 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 的情况下计算特定时间 t 中每个隔间的数量， k2、k3 等...

对于一阶微分方程组，我可以像上面那样用矩阵法求解吗？我在 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 的符号方法...

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

Solving a system of first and second order differential equations in Python

python

matrix

exponentiation

differential-equations

odeint