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

Question

我正在尝试编写一个围绕原点旋转的质量的非常简单的重力模拟。我使用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中的数字无济于事。有没有人知道为什么这可能会如此糟糕？

Answer 1

从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)顺序的参数）。