Scipy solve_ivp 没有给出预期的解决方案

如何解决Scipy solve_ivp 没有给出预期的解决方案

我需要解决一个二阶 ODE，我已经将它解耦为两个一阶 ODE。我已经尝试使用 solve_ivp 解决它，但它似乎没有提供我期望的解决方案。我提供了下面的代码。

import numpy as np
import matplotlib.pyplot as plt
from scipy.misc import derivative
from scipy.integrate import solve_ivp
from matplotlib import rc
rc('text',usetex = True)

V0 = 2*10**-10
A = 0.130383
f = 0.129576
phi_i_USR2 = 6.1
phi_Ni_USR2 = 1.2

def V_USR2(phi):
    return V0*(np.tanh(phi/np.sqrt(6)) + A*np.sin(1/f*np.tanh(phi/np.sqrt(6))))**2

def V_phi_USR2(phi):
    return derivative(V_USR2,phi)

N = np.linspace(0,66,1000)
def USR2(N,X):
    phi,g = X
    return [g,(g**2/2 - 3)*(V_phi_USR2(phi)/V_USR2(phi) + g)]

X0 = [phi_i_USR2,phi_Ni_USR2]
sol = solve_ivp(USR2,(0,66),X0,method = 'LSODA',t_eval = N)
phi_USR2 = sol.y[0]
phi_N_USR2 = sol.y[1]
N_USR2 = sol.t
plt.plot(phi_USR2,phi_N_USR2)
plt.xlabel("$\phi$")
plt.ylabel("$\phi'$")
plt.title("Phase plot for USR2")
plt.show()

solve_ivp 给我以下情节：

问题是原点附近应该有一个振荡，solve_ivp 没有很好地捕捉到它。然而，对于相同的方程和初始条件，Mathematica 给了我我想要的：

我也想在 Python 中使用相同的绘图。我在 solve_ivp 中尝试了各种方法，例如 RK45、LSODA、Radau 和 BDF，但它们都显示出相同的问题（LSODA 尝试描绘振荡但失败了，但其他方法甚至没有越过该点振荡开始的地方）。如果有人能阐明这个问题，那就太好了。提前致谢。