在python中求解微分方程组

Posted 2023-02-18

【中文标题】在python中求解微分方程组

【英文标题】：Solve system of differential equation in python

【发布时间】：2021-05-27 00:21:22

【问题描述】：

我正在尝试在 python 中求解一个微分方程组。 我有一个由两个方程组成的系统，其中我有两个变量 A 和 B。 初始条件是 A0=1e17 和 B0=0，它们同时变化。 我使用 ODEINT 编写了以下代码：

import numpy as np
from scipy.integrate import odeint

def dmdt(m,t):
    A, B = m

    dAdt = A-B
    dBdt = (A-B)*A

    return [dAdt, dBdt]

# Create time domain
t = np.linspace(0, 100, 1)

# Initial condition
A0=1e17
B0=0

m0=[A0, B0]

solution = odeint(dmdt, m0, t)

显然我得到了一个不同于预期的输出，但我不明白这个错误。 有人能帮我吗？ 谢谢

【参考方案1】：

来自A*A'-B'=0 一个结论

B = 0.5*(A^2 - A0^2)

插入第一个给出的方程

A' = A - 0.5*A^2 + 0.5*A0^2
   = 0.5*(A0^2+1 - (A-1)^2)

这意味着A动态有两个固定点大约在A0+1和-A0+1，在那个区间内增长，上固定点是稳定的。但是，在标准浮点数中，1e17 和 1e17+1 之间没有区别。如果你想看到差异，你必须单独编码。

另请注意，atol 和 rtol 介于 1e-6 和 1e-9 之间的范围内的标准误差容限与最初所述的问题规模完全不兼容，这也突出了重新调整和移位的必要性将问题转化为更可观的值范围。

将A = A0+u 与|u| 设置为1..10 的预期比例，然后给出

B = 0.5*u*(2*A0+u)

u' = A0+u - 0.5*u*(2*A0+u) = (1-u)*A0 - 0.5*u^2

现在建议将时间尺度减少A0，设置t=s/A0。另外，B = A0*v。将直接参数化插入原系统得到

du/ds = dA/dt / A0 = (A0+u-A0*v)/A0            = 1 + u/A0 - v

dv/ds = dB/dt / A0^2 = (A0+u-A0*v)*(A0+u)/A0^2 =  (1+u/A0-v)*(1+u/A0)

u(0)=v(0)=0

现在在浮点数和u 的预期范围内，我们得到1+u/A0 == 1，所以实际上u'(s)=v'(s)=1-v 给出了

u(s)=v(s)=1-exp(-s)`,

A(t) = A0 + 1-exp(-A0*t)  +  very small corrections
B(t) = A0*(1-exp(-A0*t))  +  very small corrections

s,u,v 中的系统应该可以被任何求解器在默认容差下很好地计算。