Python使用scipy.odeint求解常微分方程组报错结果异常如何解决

问题原因分析 #

• 微分方程中存在1/(1-X)**2、1/(1+X)**2这类项，当X或Z的求解结果趋近±1时，项值会直接发散到无穷大，ODE出现奇点，求解器被迫将步长压缩到极小，最终触发算力超额的警告，求解结果自然不符合预期。
• 方程存在变量误用错误：dHdt中误用了X的导数I作为阻尼项变量，应该用Z对应的导数H，同理dOdt的公式也存在对应错误。
• 求解参数设置不合理：tmax=1000、步长0.001的设置带来了1e6个时间点，过大的求解跨度和过密的采样点进一步加重了求解负担。
• 存在math、scipy.integrate重复导入的冗余代码。

修复步骤 #

1. 修正方程错误与奇点问题 #

首先修正变量误用的问题，同时给分母添加极小的保护值eps，避免X/Z接近±1时出现除零发散；dOdt可以直接简化为X和Z的作用力项差值，避免重复抄写代码出错。

2. 调整求解参数 #

先缩小tmax做小范围验证，确认求解逻辑正确后再扩大求解跨度，同时给odeint设置更小的误差阈值提升求解稳定性。

3. 修复后完整代码 #

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

##############################
# 常数定义
alpha = 1      # p1
beta = 12      # p2
mu = 0.338     # p3
omega = 0.5    # p5
gamma = 0.25   # p6
Acoef = 0.5
tmax = 20      # 先小范围测试，验证正确后再调大
t = np.arange(0.0, tmax, 0.01) # 步长无需设置过小
eps = 1e-6     # 防止分母为0的保护值
################################
# 初始条件设置
X0 = 0
I0 = 0
Z0 = 0.06
H0 = 0.04
E0 = 0
O0 = 0.02

# 修正后的微分方程
def deriv(y,t):
    X, I, Z, H, E, O = y
    # X对应的作用力项，添加保护值避免发散
    term_X = -(alpha)*X - (beta)*(X**3) - mu*I \
            + gamma*(1/( (1-X+eps)**2 )) - gamma*(1/( (1+X+eps)**2 )) \
            + Acoef*(1/( (1-X+eps)**2 ))*(np.sin(omega*t))
    # Z对应的作用力项，添加保护值避免发散
    term_Z = -(alpha)*Z - (beta)*(Z**3) - mu*H \
            + gamma*(1/( (1-Z+eps)**2 )) - gamma*(1/( (1+Z+eps)**2 )) \
            + Acoef*(1/( (1-Z+eps)**2 ))*(np.sin(omega*t))
    dXdt = I
    dIdt = term_X
    dZdt = H
    dHdt = term_Z
    dEdt = O
    dOdt = term_X - term_Z
    return dXdt, dIdt, dZdt, dHdt, dEdt, dOdt

y0 = [X0, I0, Z0, H0, E0, O0]
# 添加误差阈值参数提升求解稳定性
ret = odeint(deriv, y0, t, rtol=1e-8, atol=1e-8)
X, I, Z, H, E, O = ret.T

# 绘制相图
plt.figure(figsize=(8,6))
plt.plot(E, O)
plt.xlabel('E')
plt.ylabel('O')
plt.xlim([-10, 10])
plt.grid(True)
plt.show()
e1 = E
e2 = O

4. 可选优化 #

如果后续需要求解到tmax=1000，可以改用scipy官方更推荐的solve_ivp接口，选择适配刚性方程的Radau或BDF求解器，求解效率和稳定性会更好。

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