用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

FromA*A'-B'=0one conclusion

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