python - 这对于将重力建模为二阶 ODE 是否正确?
问题描述
这是我在这里的第一个问题,如果格式关闭,请道歉。
我想将牛顿的万有引力定律建模为 Python 中的二阶微分方程,但结果图没有意义。作为参考,这里是等式和[这里是结果] [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()
我使用这个网站作为代码的基础,我几乎没有使用 Python 作为 ODE 求解器的经验。我究竟做错了什么?谢谢!
解决方案
要集成二阶颂歌,您需要将其视为 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,导致
推荐阅读
- c++ - 我们如何在 C++ 中初始化一个所有值都为 0 的向量
- r - 如何在 R 中创建包容性分箱函数?
- gitlab - 从 Gitlab 本地服务器克隆 Android Studio 项目
- python-3.x - 使用 xpath 使用 scrapy 转到下一页
- excel - 如果另一个单元格中的文本与“A”列文本匹配,如何对“A”列的行求和(excel)
- ios - 网址重定向不适用于 WKWebview ,而适用于 UIWebview
- shell - shell - 为(已验证的)提供的数字或数字范围迭代循环
- go - 我们应该停止这种 Ticker 吗?
- kubernetes - 为 kubernetes 静态 Pod 创建的服务 IP 的时间和地点
- r - R中带有MLR包的随机森林的Information.Gain错误