求解 python 中的微分方程组

Question

我正在尝试求解 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)

显然我获得了与预期不同的输出，但我不明白其中的错误。有人能帮我吗？谢谢

Answer 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 之间没有区别。要想看区别，就得单独编码了。

另请注意，1e-6 和 1e-9 之间范围内的标准误差容限 atol 和 rtol 与最初的问题规模完全不相容陈述，还强调需要重新调整并将问题转移到更可感知的值范围内。

将 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 中的系统应该可以在默认容差下由任何求解器很好地计算。

求解 python 中的微分方程组

Solve system of differential equation in python

precision

scipy

numerical-methods

differential-equations

odeint