如何正确设置scipy集成模块'solve_ivp'中的'rtol'和'atol'，用于求解未知解析解的常微分方程组？

问：

我试图使用 solve_ivp 在 Python 中重现 ode45 求解器的一些结果。尽管所有参数、初始条件、步长以及“atol”和“rtol”（即 1e-6 和 1e-3）都相同，但我得到了不同的解决方案。这两种解决方案都收敛为周期性解决方案，但种类不同。由于 solve_ivp 使用与 ode45 相同的 rk4（5） 方法，因此最终结果中的这种差异并不完全不稳定。我们怎么知道哪一个是正确的解决方案？ 代码如下

import sys
import numpy as np
from scipy.integrate import solve_ivp
#from scipy import integrate
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): # beware of the order of the aruments
"""Return the first derivatives of y = theta1, z1, theta2, z2, z3."""
a, c, b, d, e = y

#c, s = np.cos(theta1-theta2), np.sin(theta1-theta2)
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.

y0 = np.array([3.15, -0.1, 3.13, 0.1, 0])


# Do the numerical integration of the equations of motion



sol = integrate.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

我期待 MATLAB 中的 ode45 和 Python 中的 solve_ivp 获得类似的结果。如何在python中准确重现ode45的结果？差异的原因是什么？