python - 这对于将重力建模为二阶 ODE 是否正确？

Question

这是我在这里的第一个问题，所以如果格式被关闭，我很抱歉。
我想将牛顿万有引力定律建模为 Python 中的二阶微分方程，但结果图没有意义。供引用，here's the equation和[这是结果][2]。这是我的代码

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


# dy/dt

def model(r, t):
g = 6.67408 * (10 ** -11)
m = 5.972 * 10 ** 24
M = 1.989 * 10 ** 30
return -m * r[1] + ((-g * M * m) / r ** 2)


r0 = [(1.495979 * 10 ** 16), 299195800]

t = np.linspace(-(2 * 10 ** 17), (2 * 10 ** 17))
r = odeint(model, r0, t)

plt.plot(t, r)
plt.xlabel('time')
plt.ylabel('r(t)')
plt.show()

我用过 this website作为代码的基础
我几乎没有使用 Python 作为 ODE 求解器的经验。我究竟做错了什么？谢谢!

Answer 1

要集成二阶颂歌，您需要将其视为 2 个一阶颂歌。在您发布的链接中，所有示例都是二阶的，他们这样做了。

 m d^2 r/ dt^2 = - g M m / r^2
r = u[0]
dr / dt = u[1]

(1) d/dt(u[0]) = u[1]
m * d/dt(u[1]) = -g M m / u[0]^2 =>
(2) d/dt(u[1]) = -g M / u[0]^2

在 python 中，这看起来像

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

def model(u, t):
    g = 6.67408 * (10 ** -11)
    M = 1.989 * 10 ** 30
    return (u[1], (-g * M ) / (u[0] ** 2))

r0 = [(1.495979 * 10 ** 16), 299195800]

t = np.linspace(0, 5 * (10 ** 15), 500000)
r_t = odeint(model, r0, t)
r_t = r_t[:,0]

plt.plot(t, r_t)
plt.xlabel('time')
plt.ylabel('r(t)')
plt.show()

我还对您的时间表进行了一些更改。我得到的图表看起来像这样
这对我来说很有意义。你有一个质量从一个大质量中逸出，但起始距离和速度令人难以置信，所以 r(t) 应该在时间上几乎是线性的。
然后我把299195800的速度降为0，结果