使用odeint求解不同阶数多极ODE遇结果异常问题求助

我尝试用scipy.integrate.odeint求解由1个一阶ODE和2个二阶ODE组成的多极方程组。单独求解部分方程能得到合理结果，但合并所有方程同时求解时，odeint输出完全不符合预期。试过重新定义函数、对时间t做缩放（t→1e31*t）调整后得到部分结果，但仍需排查代码问题。

变量说明 #

• p[0] = y
• p[1] = dy/dt
• p[2] = q
• p[3] = u
• p[4] = du/dt

初始代码（单独求解部分方程） #

import numpy as np
import time
from scipy.integrate import odeint

g = 1e-25
k = 4.14 * (10**-5)
t = (10**31) * np.logspace(0.247237, 2.35443, 10**7)
t = t.tolist()
z0 = [2.435*1e26, 0, 1/3400]

def my_eqs1(p, t):
    return [
        p[1],
        -2*p[2]*1.4441*(10**-33)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)*p[1] - (r*p[2])**2*p[0],
        p[2]**2*1.4441*(10**-33)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)
    ]

start_time = time.perf_counter()
r = 1e-29
z1 = odeint(my_eqs1, z0, t)
end_time = time.perf_counter()
print(end_time - start_time, "seconds")

合并求解的代码（结果异常） #

def my_eqs2(p, t):
    return [
        p[1],
        -2*p[2]*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)*p[1] - ((1e31*r*p[2])**2)*p[0],
        p[2]**2*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69),
        1e-62*p[4],
        -((k**2)+(g*p[1]/k))*p[3]
    ]
z0 = [2.435*1e26, 0, 1/3400, 1/np.sqrt(2*k), 0.35927359523]

时间缩放后的调整代码 #

t = np.logspace(0.247237, 2.35443, 10**7)

def my_eqs3(p, t):
    return [
        p[1],
        -2*p[2]*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69)*p[1] - ((1e31*r*p[2])**2)*p[0],
        p[2]**2*1.4441*(10**-2)*np.sqrt(0.31*(p[2]**-4)*(p[2]+1/3400)+0.69),
        1e-62*p[4],
        -((k**2)+(g*p[1]/k))*p[3]
    ]

可能的问题点 #

1. 量级差异引发的数值稳定性：方程组中变量和系数的量级跨度极大（如1e26、1e-33、1e31），单独求解时变量量级相对协调，合并后不同方程的量级冲突会导致odeint的积分器精度丢失或发散。
2. 时间缩放的方程变换错误：当对时间做t→τ=1e31*t缩放时，原方程的导数项需同步调整（例如dy/dt = 1e31 * dy/dτ），但当前代码仅替换了时间数组，未修正导数的缩放关系，导致方程形式不符合变换后的动力学规律。
3. 系数的人为误改：对比my_eqs1与my_eqs2/my_eqs3，发现1.4441*(10**-33)被替换为1.4441*(10**-2)，量级突变直接改变了方程的动力学行为，这很可能是结果异常的核心原因。
4. 初始条件的量级匹配问题：合并后的初始条件中，p[3]、p[4]与其他变量（如p[0]=2.435e26）的量级差距过大，可能导致积分初期就出现数值爆炸。

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