python - 用于简单重力模拟的Python scipy.integrate.odeint失败

我试图写一个非常简单的重力模拟的质量轨道的起源。我用scipy.integrate.odeint来积分微分方程。
问题是我收到以下错误消息：

ODEintWarning: Excess work done on this call (perhaps wrong Dfun type). Run with full_output = 1 to get quantitative information.
  warnings.warn(warning_msg, ODEintWarning)

除此之外，很明显有些地方出了问题——方程没有正确地积分，运动也不正确。以下是初始条件下的运动图，应给出原点周围的圆周运动：

这是代码：

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

G=1
m=1
def f_grav(y, t):
    x1, x2, v1, v2 = y
    m = t
    dydt = [v1, v2, -x1*G*m/(x1**2+x2**2)**(3/2), -x2*G*m/(x1**2+x2**2)**(3/2)]
    return dydt

t = np.linspace(0, 100, 1001)
init = [0, 1, 1, 0]
ans = odeint(f_grav, init, t)
print(ans)

x = []
y = []
for i in range (100):
    x.append(ans[i][0])
    y.append(ans[i][1])
plt.plot(x, y)
plt.show()

注意，我以前使用过这个函数，并且为SHM微分方程编写几乎相同的代码可以获得正确的结果。更改t中的数字没有帮助。有人知道为什么这次失败会这么严重吗？

最佳答案

不正确的运动可能是数值不稳定性，不管怎样，从odeint的文档中可以看出：
注意：对于新代码，使用scipy.integrate.solve_ivp来求解微分方程。
solve_ivp实际上只取边界并决定点数，这样积分法对方程是稳定的。也可以选择集成方法。

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

G=1
m=1
def f_grav(t, y):
    x1, x2, v1, v2 = y
    m = t
    dydt = [v1, v2, -x1*G*m/(x1**2+x2**2)**(3/2), -x2*G*m/(x1**2+x2**2)**(3/2)]
    return dydt

domain = (0, 100)
init = [0, 1, 1, 0]
ans = solve_ivp(fun=f_grav, t_span=domain, y0=init)

plt.plot(ans['y'][0], ans['y'][1])
plt.show()

这样我就不会收到任何警告，而且模拟看起来更好（注意，函数的参数必须按(t, y)的顺序排列）。