Templated numerical integrators in C++

I find myself writ­ing numer­i­cal inte­gra­tors for physics solvers often enough1 that I wished I didn’t have to keep rewrit­ing them every time I used a differ­ent linear alge­bra library or changed the dimen­sions what I was solv­ing.

This is exactly what C++ templates are great for—writing code that works with multi­ple types with simi­lar inter­faces. So I whipped up Euler, midpoint, and Runge-Kutta order 4 (RK4) inte­gra­tors with a generic C++ template inter­face.

Now I know templates are scary and all, but you don’t need to know them to use my code. Here’s a simple exam­ple of a physics simu­la­tor loop for solv­ing 1D Newton­ian physics with constant accel­er­a­tion:

#include "Integrators.h"

// acceleration evaluation function - just return a constant
float constAccel(float pos, float vel, float t) {
    return 1.f;
}

void simLoop() {
    IntegratorRK4 integrator;

    // 1D particle states (positions and velocities) & buffers for next state
    std::vector<float> poss, vels, nextPoss, nextVels;

    // system time state and timestep size
    float t = 0.f, dt = 1.f / 100;
    while (true) {
        // integrate accelerations to velocities & velocities to positions
        for (size_t i = 0; i < poss.size(); i++) {
            integrator.evaluate(poss[i], vels[i], constAccel, t, dt, nextPoss[i], nextVels[i]);
        }
        // take a timestep
        t += dt;
        // advance state to the next state
        poss.swap(nextPoss);
        vels.swap(nextVels);
    }
}

It’s pretty straight­for­ward: every physics iter­a­tion, make the follow­ing call for every parti­cle i: integrator.evaluate(poss[i], vels[i], constAccel, t, dt, nextPoss[i], nextVels[i]). Thanks to type infer­ence, the actual templates are hidden from you. But what are those argu­ments?

  • poss[i] and vels[i] are the particle’s state (posi­tion and veloc­ity)
  • constAccel is a func­tion that returns the accel­er­a­tion for that parti­cle (in this case a constant)
  • t and dt are the system time and timestep
  • nextPoss[i] and nextVels[i] are the outputs to the call—they get writ­ten with the next state of the parti­cle

Of course, those argu­ments don’t have to be limited to those in that exam­ple. For exam­ple, we can do:

template<typename Vector, typename Scalar>
struct DragAndConstAccel {
    const Vector accel;
    const Scalar drag;

    DragAndConstAccel(const Vector &accel, double drag) :
            accel(accel), drag(drag) {
    }

    Vector operator()(const Vector &pos, const Vector &vel, Scalar t) {
        return accel - vel * drag;
    }
};

void simLoop() {
    IntegratorRK4 integrator;
    std::vector<std::complex<double> > poss, vels, nextPoss, nextVels;
    double t = 0.0, dt = 1.0 / 100;
    while (true) {
        for (size_t i = 0; i < poss.size(); i++) {
            integrator.evaluate(poss[i],
                    vels[i],
                    DragAndConstAccel<std::complex<double>, double>(1.0, 0.25),
                    t,
                    dt,
                    nextPoss[i],
                    nextVels[i]);
        }
        t += dt;
        poss.swap(nextPoss);
        vels.swap(nextVels);
    }
}

Here I made our system posi­tions and veloc­i­ties be std::complexs (why would you ever do this?), the time and timestep be double, and replaced the constant accel­er­a­tion func­tion with a neat func­tor that can main­tain state. In this case, the accel­er­a­tion func­tor applies a constant accel­er­a­tion with a drag accel­er­a­tion that’s linear to the veloc­ity.

But look how the integrator.evaluate() call and the simu­la­tion loop look almost the same as before: that’s the beauty of templates2. You can change around your vector & time types and still keep the same loop code.

You wanna see some­thing crazy? Let’s stick in Eigen Vector3f and make our accel­er­a­tion func­tion be a C++11 lambda:

void simLoop() {
    using namespace Eigen;
    IntegratorRK4 integrator;
    std::vector<Vector3f> poss, vels, nextPoss, nextVels;
    double t = 0.0, dt = 1.0 / 100;
    while (true) {
        for (size_t i = 0; i < poss.size(); i++) {
            integrator.evaluate(poss[i], vels[i], [](const Vector3f &pos, const Vector3f &vel, double t) {
                return Vector3f(1.f, 0.f, -1.f) - vel * 0.25f;
            }, t, dt, nextPoss[i], nextVels[i]);
        }
        t += dt;
        poss.swap(nextPoss);
        vels.swap(nextVels);
    }
}

Neat, innit? Anyways, all the inte­gra­tors are avail­able in a single header file under the MIT/X11 license. If you plan to use it, it would be nice to give me an email shout.

You may have noticed that the inte­gra­tors don’t handle mass at all—that’s why you pass it an accel­er­a­tion eval­u­a­tion func­tion rather than a force eval­u­a­tion func­tion, as it has no knowl­edge of mass. Extend­ing the inte­gra­tors to use mass is a triv­ial exer­cise left up to the reader. 🙂

Also, weath­ered C++ coders may notice that I spec­ify the inte­gra­tors’ inter­face using a CRTP base class. There’s no compelling soft­ware engi­neer­ing reason for this; I actu­ally find CRTP static poly­mor­phism to be a mostly useless construct, as you can’t swap around your derived classes with­out recom­pil­ing. But it does make for a good inter­face spec.

  1. Oh dear god it hurt me to write that… []
  2. Yes, yes, I know that “the beauty of templates” is a contra­dic­tion in terms []