4

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!

baddieb
  • 41
  • 2
  • oopsie doopsie, the result didn’t post...Basically, the graph does not look right at all – baddieb Feb 28 '21 at 18:43
  • you need parenthesis around the `r ** 2` – NobbyNobbs Feb 28 '21 at 19:30
  • @NobbyNobbs I did that, but nothing's changed. Also, this error comes up "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)" – baddieb Feb 28 '21 at 20:19
  • 1
    You lost some links, equation and plot are no longer present. You need to check again the notation. It can be confusing, but in gravity you have both radius and position in the formula, you mixed them here in the wrong way. Where does the linear term come from? Also, check the difference of force and acceleration. – Lutz Lehmann Mar 01 '21 at 08:43

1 Answers1

3

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 plot

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

plot

EntropicFroggy
  • 143
  • 10