如何设置solve_ivp的rtol/atol？及与MATLAB ode45结果差异排查

我尝试用Scipy的solve_ivp复现MATLAB中ode45求解器的结果，参数、初始条件、步长以及rtol（1e-3）和atol（1e-6）都保持一致，但得到的解却不同。两种解都收敛到周期解，但周期特性存在差异。由于solve_ivp与ode45采用相同的RK4(5)方法，这种结果差异令人费解。我想知道哪一个解是正确的，以及如何在Python中精准复现ode45的结果，差异产生的原因是什么？

相关Python代码 #

import sys
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
from matplotlib.patches import Circle

# Pendulum rod lengths (m), bob masses (kg).
L1, L2, mu, a1 = 1, 1, 1/5, 1
m1, m2, B = 1, 1, 0.1
# The gravitational acceleration (m.s-2).
g = 9.81

# The forcing frequency,forcing amplitude
w, a_m = 10, 4.5
A = (a_m * w**2) / g
A1 = a_m / g

def deriv(t, y, mu, a1, B, w, A):
    """Return the first derivatives of y = theta1, z1, theta2, z2, z3."""
    a, c, b, d, e = y
    
    adot = c
    cdot = (-(1 - A*np.sin(e))*(((1+mu)*np.sin(a)) - (mu*np.cos(a-b)*np.sin(b))) 
            - ((mu/a1)*((d**2) + (a1*np.cos(a-b)*c**2))*np.sin(a-b)) 
            - (2*B*(1 + (np.sin(a-b))**2)*c) 
            - ((2*B*A/w)*(2*np.sin(a) - (np.cos(a-b)*np.sin(b)))*np.cos(e))) / (1 + mu*(np.sin(a-b))**2)
    
    bdot = d
    ddot = ((-a1*(1+mu)*(1 - A*np.sin(e))*(np.sin(b) - (np.cos(a-b)*np.sin(a)))) 
            + (((a1*(1+mu)*c**2) + (mu*np.cos(a-b)*d**2))*np.sin(a-b)) 
            - ((2*B/mu)*(((1 + mu*(np.sin(a-b))**2)*d) + (a1*(1-mu)*np.cos(a-b)*c))) 
            - ((2*B*a1*A/(w*mu))*(((1+mu)*np.sin(b)) - (2*mu*np.cos(a-b)*np.sin(a)))*np.cos(e))) / (1 + mu*(np.sin(a-b))**2)
    edot = w
    return adot, cdot, bdot, ddot, edot

# Initial conditions: theta1, dtheta1/dt, theta2, dtheta2/dt, z3.
y0 = np.array([3.15, -0.1, 3.13, 0.1, 0])

# Numerical integration of the equations of motion
sol = solve_ivp(deriv, [0, 40000], y0, args=(mu, a1, B, w, A), 
                method='RK45', t_eval=np.arange(0, 40000, 0.005), 
                dense_output=True, rtol=1e-3, atol=1e-6)

T = sol.t
Y = sol.y

差异原因分析 #

• 自适应步长实现细节差异：尽管两者都用RK4(5)方法，但MATLAB和Scipy在步长调整的具体逻辑（比如误差阈值的加权方式、步长缩放因子的计算）上存在细微差别。长时间积分（如40000单位时间）下，微小的步长差异会累积，导致周期特性偏离。
• 浮点数运算底层差异：MATLAB和Python依赖的数学库（BLAS/LAPACK）版本、编译器优化不同，会引入微小的舍入误差，在混沌或强非线性系统中，这种误差会被快速放大，最终导致解的差异。
• 导数函数实现偏差：需严格核对MATLAB和Python的导数代码，确保运算符优先级、括号位置、三角函数参数完全一致。比如复杂分式的分子分母括号是否匹配，避免因计算顺序不同导致结果偏差。
• t_eval插值影响：Python中指定t_eval会强制在固定时间点插值，而MATLABode45默认输出自适应步长的计算点，插值过程可能引入额外误差。

验证解正确性的方法 #

• 短区间对比：先积分短时间（如0到10），对比两种方法的结果。若初始阶段一致，说明是长时间误差累积导致；若初始就不同，需检查导数函数或参数设置。
• 高精度验证：将rtol设为1e-6、atol设为1e-9，重新计算。若此时两种方法结果趋近一致，说明原差异是精度设置不足导致。
• 解析解对照：若系统存在简单场景的解析解，用解析解验证两种求解器的正确性。

精准复现ode45结果的方案 #

• 匹配步长控制参数：设置solve_ivp的max_step参数与MATLABode45的默认MaxStep一致（MATLAB默认MaxStep为积分区间的1/10，即4000），同时可指定initial_step为MATLAB的初始步长。
• 移除t_eval插值：先不指定t_eval，让solve_ivp自动输出计算点，再与MATLAB的输出点对比。若一致，再插值到目标时间点。
• 严格对齐导数函数：逐行核对MATLAB和Python的导数计算代码，确保每一个运算、系数、括号完全匹配，消除代码层面的偏差。
• 使用相同的浮点数精度：在Python中设置np.float64类型，与MATLAB的默认精度对齐，避免因精度类型差异导致的计算偏差。

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