我试图在Python中使用solve_ivp重现ode 45求解器的一些结果。尽管所有参数、初始条件、步长以及“atol”和“rtol”（分别为1 e-6和1 e-3）相同，但得到的解不同。两个解都收敛到周期解，但类型不同。由于solve_ivp使用相同的rk 4（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中的ode 45和Python中的solve_ivp得到类似的结果。我如何才能准确地重现Python中的ode 45的结果？出现差异的原因是什么？