每次运行Python odeint微分方程代码输出图形不一致的问询

问题背景 #

在Anaconda Spyder中运行下述Python代码时，每次生成的输出图形均存在差异：所有运行的图形在x轴60微秒之前完全一致，但60微秒之后的图形每次都不相同。预期是每次运行都能得到与60微秒之前一致的重复性曲线。

原代码 #

import math
from scipy.integrate import odeint
import numpy as np
import matplotlib.pyplot as plt

#Input
c = 1500
rho= 1000
S= 0.07
mu =0.001
R0 =0.00001
PA=5
ω =100000
PO=1 
PGO=1
Pa=1
A= (ω*R0/c)
B=Pa*(10**5)/(rho*(ω**2)*(R0**2))
C=(2*S)/(rho*(ω**2)*(R0**3))
D=(4*mu)/(rho*ω*(R0**2))
          
#Model Definition
def f(X, T):
    β = X[0]
    z = X[1]
    dβdT = z
    dzdT = (1/(β-(A*β*z)+(A*D)))*((-1.5*z**2)+(0.5*A*z**3)+(B*PGO/(β**3))-(C/β)-(D*z/β)-(B*PO)+(B*PA*math.sin(T))+
                              (A*B*PGO*z/(β**3))-(A*C*z/β)-(A*D*(z**2)/β)-(A*B*PO*z)+(A*B*PA*z*math.sin(T))-
                              (3*A*B*PGO*z/(β**3))+(A*C*z/β)+(A*D*(z**2)/β)+(A*B*PA*β*math.cos(T)))
    return [dβdT, dzdT]

#Initial Boundary Condition
X0 = [1, 0]

#Time Interval
T = np.linspace(0, 14, 250)

sol = odeint(f, X0, T)
β = sol[:,0]
z = sol[:,1]

t= (T/ω)*(10**6)   #micro seconds
R= β*R0*(10**6)    #micro meters
V= z*R0*ω

plt.plot(t,R)
plt.xlabel('t')
plt.legend('R')

问题原因 #

• 该微分方程属于刚性系统，原代码中dzdT的表达式存在大量冗余计算，放大了数值舍入误差。
• 积分后期，分母项β-(A*β*z)+(A*D)可能接近零，导致dzdT计算值剧烈波动，引发数值发散。
• odeint默认的误差容忍度（rtol=1e-6, atol=1e-9）不足以约束该系统的误差积累，多次运行时误差走向不同，导致结果不一致。

解决措施 #

1. 简化ODE表达式，消除冗余计算 #

合并dzdT中的同类项，减少计算过程中的误差来源：

def f(X, T):
    β = X[0]
    z = X[1]
    dβdT = z
    # 合并冗余项，简化计算
    term1 = -1.5 * z**2 + 0.5 * A * z**3
    term2 = B * PGO / (β**3) - C/β - D*z/β - B*PO + B*PA*math.sin(T)
    term3 = -2 * A*B*PGO*z/(β**3) - A*B*PO*z + A*B*PA*z*math.sin(T) + A*B*PA*β*math.cos(T)
    denominator = β - A*β*z + A*D
    dzdT = (term1 + term2 + term3) / denominator
    return [dβdT, dzdT]

2. 严格约束积分误差 #

调用odeint时手动设置更小的误差容忍度，强制积分器更精确地跟踪解：

sol = odeint(f, X0, T, rtol=1e-10, atol=1e-12)

3. 更换为现代积分器solve_ivp #

odeint是老旧接口，改用scipy.integrate.solve_ivp并选择专门的刚性求解器（如Radau），适配刚性系统的求解需求：

from scipy.integrate import solve_ivp

# 注意solve_ivp的函数接口为f(t, y)，与odeint的f(y, t)相反
def f_new(T, X):
    β = X[0]
    z = X[1]
    dβdT = z
    term1 = -1.5 * z**2 + 0.5 * A * z**3
    term2 = B * PGO / (β**3) - C/β - D*z/β - B*PO + B*PA*math.sin(T)
    term3 = -2 * A*B*PGO*z/(β**3) - A*B*PO*z + A*B*PA*z*math.sin(T) + A*B*PA*β*math.cos(T)
    denominator = β - A*β*z + A*D
    dzdT = (term1 + term2 + term3) / denominator
    return [dβdT, dzdT]

# 使用Radau刚性求解器，设置严格误差阈值
sol = solve_ivp(f_new, [T[0], T[-1]], X0, t_eval=T, method='Radau', rtol=1e-10, atol=1e-12)
β = sol.y[0]
z = sol.y[1]

4. 监控分母项的稳定性 #

在积分过程中添加对分母项的监控，若其接近零，需重新审视模型的物理合理性，或添加约束避免β/z进入导致分母趋近于零的区域：

# 在f函数中添加监控（可选）
denominator = β - A*β*z + A*D
if abs(denominator) < 1e-10:
    print(f"Warning: Denominator near zero at T={T}, β={β}, z={z}")

效果验证 #

通过上述调整后，积分过程的数值稳定性会大幅提升，多次运行后60微秒之后的曲线将保持一致，符合预期的重复性要求。

内容的提问来源于stack exchange，提问作者user20838423