Is This Correct For Modeling Gravity As A Second Order ODE?

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, b

Solution 1:

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