orbit.ipynb

Kernel: Python 3 (ipykernel)

Orbiting the Sun

Modeling and Simulation in Python

Copyright 2021 Allen Downey

License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

# install Pint if necessary

try:
    import pint
except ImportError:
    !pip install pint

# download modsim.py if necessary

from os.path import basename, exists

def download(url):
    filename = basename(url)
    if not exists(filename):
        from urllib.request import urlretrieve
        local, _ = urlretrieve(url, filename)
        print('Downloaded ' + local)
    
download('https://github.com/AllenDowney/ModSimPy/raw/master/modsim.py')

# import functions from modsim

from modsim import *

In a previous example, we modeled the interaction between the Earth and the Sun, simulating what would happen if the Earth stopped in its orbit and fell straight into the Sun.

Now let's extend the model to two dimensions and simulate one revolution of the Earth around the Sun, that is, one year.

At perihelion, the distance from the Earth to the Sun is 147.09 million km and its velocity is 30,290 m/s.

r_0 = 147.09e9     # initial distance m
v_0 = 30.29e3      # initial velocity m/s

Here are the other constants we'll need, all with about 4 significant digits.

G = 6.6743e-11     # gravitational constant N / kg**2 * m**2
m1 = 1.989e30      # mass of the Sun kg
m2 = 5.972e24      # mass of the Earth kg
t_end = 3.154e7    # one year in seconds

Exercise: Put the initial conditions in a State object with variables x, y, vx, and vy. Create a System object with variables init and t_end.

# Solution goes here

# Solution goes here

Exercise: Write a function called universal_gravitation that takes a State and a System and returns the gravitational force of the Sun on the Earth as a Vector.

Test your function with the initial conditions; the result should be a Vector with approximate components:

x   -3.66e+22
y   0

# Solution goes here

# Solution goes here

Exercise: Write a slope function that takes a timestamp, a State, and a System and computes the derivatives of the state variables.

Test your function with the initial conditions. The result should be a sequence of four values, approximately

0.0, -30290.0, -0.006, 0.0

# Solution goes here

# Solution goes here

Exercise: Use run_solve_ivp to run the simulation. Save the return values in variables called results and details.

# Solution goes here

You can use the following function to plot the results.

from matplotlib.pyplot import plot

def plot_trajectory(results):
    x = results.x / 1e9
    y = results.y / 1e9

    make_series(x, y).plot(label='orbit')
    plot(0, 0, 'yo')

    decorate(xlabel='x distance (million km)',
             ylabel='y distance (million km)')

plot_trajectory(results)

You will probably see that the earth does not end up back where it started, as we expect it to after one year. The following cells compute the error, which is the distance between the initial and final positions.

error = results.iloc[-1] - system.init
error

offset = Vector(error.x, error.y)
vector_mag(offset) / 1e9

The problem is that the algorithm used by run_solve_ivp does not work very well with systems like this. There are two ways we can improve it.

run_solve_ivp takes a keyword argument, rtol, that specifies the "relative tolerance", which determines the size of the time steps in the simulation. Lower values of rtol require smaller steps, which yield more accurate results. The default value of rtol is 1e-3.

Exercise: Try running the simulation again with smaller values, like 1e-4 or 1e-5, and see what effect it has on the magnitude of offset.

The other way to improve the results is to use a different algorithm. run_solve_ivp takes a keyword argument, method, that specifies which algorithm it should use. The default is RK45, which is a good, general-purpose algorithm, but not particularly good for this system. One of the other options is RK23, which is usually less accurate than RK45 (with the same step size), but for this system it turns out to be unreasonably good, for reasons I don't entirely understand. Yet another option is 'DOP853', which is particularly good when rtol is small.

Exercise: Run the simulation with one of these methods and see what effect it has on the results. To get a sense of how efficient the methods are, display details.nfev, which is the number of times run_solve_ivp called the slope function.

details.nfev

Animation

You can use the following draw function to animate the results, if you want to see what the orbit looks like (not in real time).

xlim = results.x.min(), results.x.max()
ylim = results.y.min(), results.y.max()

def draw_func(t, state):
    x, y, vx, vy = state
    plot(x, y, 'b.')
    plot(0, 0, 'yo')
    decorate(xlabel='x distance (million km)',
             ylabel='y distance (million km)',
             xlim=xlim,
             ylim=ylim)

# animate(results, draw_func)