这对于将重力建模为二阶 ODE 是否正确? [英] Is this correct for modeling gravity as a second order ODE?

问题描述

这是我在这里的第一个问题，如果格式不正确，请见谅.

我想在 Python 中将牛顿万有引力定律建模为二阶微分方程，但结果图没有意义.作为参考，

这对我来说很有意义.你有一个质量从一个大质量中逸出，但起始距离和速度令人难以置信，所以 r(t) 应该在时间上几乎是线性的.然后我把299195800的速度降为0，结果

This is my first question on here, so apologies if the formatting is off.

I want to model Newton's Universal Law of Gravitation as a second-order differential equation in Python, but the resulting graph doesn’t make sense. For reference, here's the equation and [here's the result][2]. This is my code

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

I used this website as a base for the code I have virtually no experience with using Python as an ODE solver. What am I doing wrong? Thank you!

解决方案

To integrate a second order ode, you need to treat it like 2 first order odes. In the link you posted all the examples are second order, and they do this.

 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

In python this looks like

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

I also made some changes to your time list. What I got for the graph looks like so

which makes sense to me. You have a mass escaping away from a large mass but at an incredible starting distance and speed, so r(t) should pretty much be linear in time. Then I brought the speed of 299195800 down to 0, resulting in

这篇关于这对于将重力建模为二阶 ODE 是否正确?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！