用于简单重力模拟的 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 的文档中，无论哪种方式，不正确的运动都可能是数值不稳定性:

note: For new code, use scipy.integrate.solve_ivp to solve a differential equation.

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))。