Example: Planetary Orbit

Let’s consider an example of integrating the orbit of Earth around the Sun (we’ll consider the mass of the Earth to be tiny compared to the Sun so we can neglect the Sun’s motion around the center of mass).

Kepler’s third law tells us:

\[4\pi^2 a^3 = G M_\star P^2\]

where \(a\) is the semi-major axis and \(P\) is the period of the orbit. If we work in units of solar masses, astronomical units, and years, then this becomes:

\[a^3 = P^2\]

which implies that \(G M_\star = 4 \pi^2\).

The evolution of the planet is given by

\[\dot{\bf r} = {\bf v}\]
\[\dot{\bf v} = - \frac{G M_\star {\bf r}}{r^3}\]

where \({\bf r} = (x, y)\) and \({\bf v} = (v_x, v_y)\).

Euler’s Method

The simplest method to advance the solution begins with a first-order accurate difference approximation to the derivative:

\[\frac{{\bf r}^{n+1} - {\bf r}^n}{\tau} = {\bf v}^n + \mathcal{O}(\tau)\]
\[\frac{{\bf v}^{n+1} - {\bf v}^n}{\tau} = {\bf a}({\bf r}^n,{\bf v}^n) + \mathcal{O}(\tau)\]

or

\[{\bf r}^{n+1} = {\bf r}^n + \tau {\bf v}^n + \mathcal{O}(\tau^2)\]
\[{\bf v}^{n+1} = {\bf v}^n + \tau {\bf a}({\bf r}^n,{\bf v}^n) + \mathcal{O}(\tau^2)\]

This shows that each time we take a step, we make an error (called truncation error) of \(\mathcal{O}(\tau^2)\). Over the course of the orbit we will take \(N = T / \tau\) steps, so the global error will be \(\mathcal{O}(\tau)\).

This method is first-order accurate in time.

Graphically, Euler’s method looks like:

_images/rk2_Euler.png

Fig. 9 Euler’s method: we evaluate the slope at the start of the timestep and then follow this through to the new time.

Implementation

We’ll store the state of the solution at a single instance in time in a struct:

struct OrbitState {
    double t;
    double x;
    double y;
    double vx;
    double vy;
};

and we’ll save the solution at each step in a vector:

std::vector<OrbitState> orbit_history{};

This corresponds to the “array-of-structs” storage scheme we discussed earlier.

To make our code flexible, we’ll write the following functions:

  • A righthand side function, with the forward declaration:

    OrbitState rhs(const OrbitState& state);
    

    we can evaluate this as needed by passing in an OrbitState and returning the derivative of that state with respect to time.

  • A function to write out the history:

    void write_history(const std::vector<OrbitState>& history);
    

    We can adapt this as desired to write to the screen or a file.

  • An integration function:

    std::vector<OrbitState> integrate(const double a, const double tmax, const double dt);
    

    This manages the time-stepping, implementing Euler’s method and storing the solution in our history vector each step. At the end of integration it will return the history.

try it…

Implement the solver by filling in the details of the functions described above.

Tip

We can plot the solution using gnuplot.

First, when you run, redirect the output to a file:

./orbit_example > orbit.dat

Now run gnuplot:

gnuplot

A simple plot in gnuplot is made using the plot command, specifying the file with the data, the columns to use for x and y axes, and any styling (like w l for with-lines):

plot 'orbit.dat' using 2:3 w  l

You can save the plot as:

set term png
set output "orbit.png"
replot

This will make something like:

_images/orbit.png

For more info, see: http://www.gnuplot.info/