使用odeint求解快速振荡微分方程报错结果异常的解决方法

• 抛出Excess work done警告的本质是odeint底层调用的LSODA求解器，在默认最大步数限制下无法完成指定区间的积分。你的目标函数存在快振荡特性，同时参数t、m、initial_value的取值跨了十到二十余个数量级，属于典型的多尺度刚性问题，默认求解配置下的步长策略无法追踪快速变化的振荡分量，步长收缩到下限仍达不到精度要求时就会终止计算，返回无效值。
• 代码本身存在实现缺陷：一是在两层循环内部重复定义微分方程函数dzdt，产生不必要的性能开销；二是混用mpmath库的mp.coth与numpy数组运算，mpmath为任意精度运算设计，对numpy数组的广播支持极差，会大幅拖慢计算速度，还可能触发类型不匹配的数值异常；三是while循环的索引逻辑存在漏洞，持续累加j最终会触发initial_value的索引越界错误。
• 积分节点设置不合理：直接使用跨26个数量级的logspace序列作为积分输出节点，要求求解器覆盖从极短瞬态到长时演化的全区间，默认仅500步的最大迭代步数完全无法覆盖这么大跨度的演化过程。

• 首先清理无效依赖、修正基础逻辑：删除未使用的导入项，将mp.coth替换为numpy原生实现1/np.tanh(x)，完全移除mpmath调用；把while循环改为直接遍历initial_value的索引，避免索引越界。
• 调整求解器配置适配刚性振荡问题：
  • 调大odeint的mxstep参数，对快振荡场景可设置为100000（即1e5），给求解器足够的迭代步数追踪振荡过程；
  • 若调整mxstep后仍有求解失败的区间，可更换为刚性问题专用求解器，比如scipy.integrate.solve_ivp提供的Radau或BDF方法，对多尺度刚性问题的稳定性远优于默认的LSODA；
  • 调试阶段可开启full_output=1参数，拿到求解失败的具体诊断信息，精准定位问题参数区间。
• 从根源降低求解难度：对微分方程做无量纲化处理，通过变量代换消去方程中跨数量级的常数项和参数，把方程转化为参数在O(1)量级的形式，可大幅降低求解器的步长压力；积分节点不需要全程等密度取log间隔，在阻尼小、振荡快的小T区间适当加密节点，在阻尼大、振荡衰减的大T区间适当放疏节点，平衡计算精度和速度。

import numpy as np
import scipy.integrate as spi
import time

# 提前定义常数项，避免循环内重复计算
COEFF_A = 3 * 1.4441e-6 * np.sqrt(0.69)
COEFF_B = 1.5 * np.sqrt(0.69) * 1e-6 * 1.4441

initial_value = np.logspace(24, 27, 100)
t = np.logspace(-20, 6, 100)
m = np.logspace(0, 6, 100)

start_time = time.perf_counter()
phi_m = {}
phi_m_prime = {}

# 替换有漏洞的while循环为直接遍历索引
for j in range(len(initial_value)):
    i = np.pi * 2.435 * initial_value[j]
    phi = []
    phi_prime = []
    for mx in m:
        # 传参固定当前m值，替换coth为numpy原生实现
        def dzdt(z, T, mx_val):
            coth_term = 1 / np.tanh(COEFF_B * mx_val * T)
            damp = COEFF_A * mx_val * coth_term
            return [z[1], -damp * z[1] - z[0]]
        
        z0 = [i, 0]
        ts = t / mx
        # 调大mxstep适配快振荡求解，开启full_output方便调试
        zs, info = spi.odeint(dzdt, z0, ts, args=(mx,), mxstep=100000, full_output=True)
        phi.append(zs[-1, 0])
        phi_prime.append(zs[-1, 1])
    phi_m[j] = phi
    phi_m_prime[j] = phi_prime

end_time = time.perf_counter()
print(end_time - start_time, "seconds")

注：如果修正后仍有部分参数区间计算失败，优先对该区间的方程做无量纲化，或者换用solve_ivp的刚性求解器，不要盲目继续调大mxstep，否则会让计算速度变得极慢。

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