In Class: 2nd Order Runge-Kutta Orbit

The main utility of functions is that we can reuse them in a code. Let’s look at this by extending our Euler integration of an orbit to do 2nd-order Runge-Kutta integration.

The Euler method was based on a first-order difference approximation to the derivative. But a centered-difference approximation to \(df/dt\) is second order accurate at the midpoint in time, so we can try to update our system in the form:

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

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

Then the updates are:

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

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

This is locally third-order accurate (but globally second-order accurate), but we don’t know how to compute the state at the half-time.

To find the \(n+1/2\) state, we first use Euler’s method to predict the state at the midpoint in time. We then use this provisional state to evaluate the accelerations at the midpoint and use those to update the state fully through \(\tau\).

The two step process appears as:

\[{\bf r}^\star = {\bf r}^n + (\tau / 2) {\bf v}^n\]

\[{\bf v}^\star = {\bf v}^n + (\tau / 2) {\bf a}^n\]

then we use this for the full update:

\[{\bf r}^{n+1} = {\bf r}^n + \tau \, {\bf v}^\star\]

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

Graphically this looks like the following:

First we take a half step and we evaluate the slope at the midpoint:

Fig. 8 The predictor part of RK2 integration.

Fig. 9 The corrector / final update of RK2 integration.

Notice how the final step (the red line) is parallel to the slope we computed at \(t^{n+1/2}\). Also note that the solution at \(t^{n+1}\) is much closer to the analytic solution than in the figure from Euler’s method.

To implement this, it would be nice to have a function that updates the state based on the derivatives, like:

OrbitState update_state(const OrbitState& state, const double dt, const OrbitState& state_derivs);

This way we can easily do the update both times as needed via a simple function call.

Tip

Later we’ll learn how to use operator overloading to do something like:

OrbitState state{};
double dt{};
OrbitState state_derivs{};

...

auto state_new = state + dt * state_derivs;

Here’s a solution:

solution

Listing 34 orbit_rk2_example.cpp

#include <iostream>
#include <iomanip>
#include <vector>
#include <cmath>

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

OrbitState rhs(const OrbitState& state);
OrbitState update_state(const OrbitState& state, const double dt, const OrbitState& state_derivs);
void write_history(const std::vector<OrbitState>& history);
std::vector<OrbitState> integrate(const double a, const double tmax, const double dt_in);

const double GM = 4.0 * M_PI * M_PI;   // G * Mass in AU, year, solar mass units

OrbitState rhs(const OrbitState& state) {

    OrbitState dodt{};

    // dx/dt = vx; dy/dt = vy

    dodt.x = state.vx;
    dodt.y = state.vy;

    // d(vx)/dt = - GMx/r**3; d(vy)/dt = - GMy/r**3

    double r = std::sqrt(state.x * state.x + state.y * state.y);

    dodt.vx = - GM * state.x / std::pow(r, 3);
    dodt.vy = - GM * state.y / std::pow(r, 3);

    return dodt;

}

OrbitState update_state(const OrbitState& state, const double dt, const OrbitState& state_derivs) {

    OrbitState state_new{};

    state_new.t = state.t + dt;
    state_new.x = state.x + dt * state_derivs.x;
    state_new.y = state.y + dt * state_derivs.y;
    state_new.vx = state.vx + dt * state_derivs.vx;
    state_new.vy = state.vy + dt * state_derivs.vy;

    return state_new;
}

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

    for (auto o : history) {
        std::cout << std::setw(14) << o.t
                  << std::setw(14) << o.x
                  << std::setw(14) << o.y
                  << std::setw(14) << o.vx
                  << std::setw(14) << o.vy << std::endl;

    }

}

std::vector<OrbitState> integrate(const double a, const double tmax, const double dt_in) {

    // how the history of the orbit

    std::vector<OrbitState> orbit_history{};

    // set initial conditions
    OrbitState state{};

    // assume circular orbit on the x-axis, counter-clockwise orbit

    state.t = 0.0;
    state.x = a;
    state.y = 0.0;
    state.vx = 0.0;
    state.vy = std::sqrt(GM / a);

    orbit_history.push_back(state);

    double dt = dt_in;

    // integration loop
    while (state.t < tmax) {

        if (state.t + dt > tmax) {
            dt = tmax - state.t;
        }

        // get the derivatives
        auto state_derivs = rhs(state);

        // get the derivatives at the midpoint in time
        auto state_star = update_state(state, 0.5 * dt, state_derivs);
        state_derivs = rhs(state_star);

        // update the state
        state = update_state(state, dt, state_derivs);

        orbit_history.push_back(state);
    }

    return orbit_history;

}

int main() {

    double tmax = 1.0;
    double dt = 0.05;
    double a = 1.0;      // 1 AU

    auto orbit_history = integrate(a, tmax, dt);
    write_history(orbit_history);

}