求解二阶ODE时出现警告且解轨迹不完整的问题求助

问题现象 #

• 当参数g=0时，ODE求解正常，相轨迹完整
• 当g取非零值（如0.9）时，触发ODEintWarning警告，提示使用full_output=1获取定量信息
• 生成的相轨迹不完整，调整步长、mxstep参数、可视化限制均无法解决

问题根源 #

1. H的实数有效性问题：H的表达式包含平方根，当根号内的项为负时，会产生复数，导致数值求解崩溃。g≠0时，g*(phi_1**3)项在phi_1为负时会放大负向值，可能让根号内整体变为负数，使H无法得到有效实数，引发求解中断。
2. 数值求解稳定性问题：H出现异常后，微分方程右端项会出现NaN或无穷大，odeint无法继续积分，导致轨迹提前终止。

修复方案 #

1. 确保H的实数性：计算H前检查根号内数值，若为负则强制设为极小正值，避免复数产生
2. 拆分时间积分方向：原时间区间同时包含正负时间，可能加剧数值不稳定，拆分正向和反向时间分别积分后合并轨迹
3. 启用full_output排查细节：添加full_output=1参数，可获取每个初始条件下的求解错误详情，定位具体中断点

修改后的代码 #

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

m=0.5
g=0.9
# 拆分时间区间为正向和反向，分别积分以提升稳定性
time_forward = np.linspace(0, 10, 10000)
time_backward = np.linspace(0, -10, 10000)

def system(phi,t):
    phi_1=phi[0]
    phi_2=phi[1]
    dphi1_dt=phi_2
    # 计算根号内项，添加极小值避免负数
    sqrt_term = (8*np.pi/3)*(0.5*(phi_2**2)) + (0.5*(phi_1**2)*m**2) + (g*(phi_1**3))
    sqrt_term = max(sqrt_term, 1e-10)  # 确保根号内非负
    H=np.sqrt(sqrt_term)
    dphi2_dt=(-3*H*phi_2)-(phi_1*m**2)-(3*g*phi_1**2)
    return [dphi1_dt, dphi2_dt]

init=np.linspace(-0.2,0.2,5)
init_2=np.linspace(-0.2,0.2,5)

plt.figure(figsize=(8,6))
for i in init:
    for j in init_2:
        # 正向时间积分
        phi_forward = odeint(system, [i,j], time_forward, mxstep=100000)
        # 反向时间积分
        phi_backward = odeint(system, [i,j], time_backward, mxstep=100000)
        # 合并轨迹：反转反向轨迹并去掉重复初始点
        full_phi = np.vstack([phi_backward[::-1][1:], phi_forward])
        plt.plot(full_phi[:,0], full_phi[:,1])

plt.xlabel("$\phi$",fontsize=12)
plt.ylabel("$d \phi/dt$",fontsize=12)
plt.xticks(fontsize=12)
plt.yticks(fontsize=12)
plt.xlim(-0.2,0.2)
plt.ylim(-0.2,0.2)
plt.show()

额外说明 #

• sqrt_term = max(sqrt_term, 1e-10)确保H始终为实数，避免复数导致的求解中断
• 拆分时间区间分别积分再合并，能有效降低双向时间积分带来的数值不稳定
• 若需要进一步排查问题，可在调用odeint时添加full_output=1参数，通过返回的第二个元素查看详细求解日志

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