Matplotlib.odeint 为 acc = Gm/r**2 的简单公式给出了错误的解决方案

问题描述

我编写了一个程序来计算任何物体（在这种情况下是月球）在地球引力作用下的运动，在初始速度为零的情况下，月球应该在地球上沿直线来回摆动，直到它完全正确地停在地球上在地球顶部，所以绘制时间与距离的关系图我应该得到一个类似于阻尼信号的图，相反，我得到的是一个无限递减的图。

我试过了:-

1. 给月球赋予不同的初始速度，但得到了相同的结果，好像我用odeint求解微分方程的方式不对？不确定。对编码非常陌生。

2. 假设 1000 秒不足以让这种情况发生，所以我将时间增加到 1e+5,1e+10,1e+20，似乎 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)

我应该使用其他函数来求解这个微分方程吗？

3. 将质量减少到 10,20，将 G 和 r 减少到 10,10，并收到与情况 (2) 相同的警告

任何反馈都有帮助

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

G=6.67408e-11 #N-m2/kg2 #N-m2/kg2
m1 = 5.972e+24 # kg,mass of earth
m2 = 7.348e+22  # kg,mass of moon


def dvdt(S,t):
    #  v = dr./dt,so a = dv/dt
    r,v = S
    return [v,-G*m1 / r**2]

# initial values
r10 = 0 # position of earth in meters
r20 = 4e+8 # position of moon from earth in meters
v10 = 0 # veLocity of earth m/s
v20 = 0 # veLocity of moon relative to earth m/s
S0 = [r20,v20]
t = np.linspace(0,1000,100)

# solving the differential eqn
acc = odeint(dvdt,S0,t)
r,v = acc.T

# plotting
plt.plot(t,r)
plt.xlabel('Time')
plt.ylabel('distance between earth and moon')
plt.show()

解决方法

该场景假设两个天体“异相”，因此它们彼此之间的行为就像暗物质一样，没有弱电或强核力。

地球半径以下的有效重力由当前半径内的质量决定，外壳的影响加为零。

力矢量始终指向中心，对于一维运动，力的符号必须始终与坐标的符号相反。

总共是这样

R1 = 6.7e+6 # m

acc = -sign(r) * G*(m1*min(1,abs(r)/R1)**3) / abs(r)**2
    = -G*m1 * r/max(R1,abs(r))**3

如果实现这一点，就会得到预期的振荡图