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 dif­fer­ent lin­ear alge­bra library or changed the dimen­sions what I was solv­ing.

This is exact­ly what C++ tem­plates are great for—writing code that works with mul­ti­ple types with sim­i­lar inter­faces. So I whipped up Euler, mid­point, and Runge-Kut­ta order 4 (RK4) inte­gra­tors with a gener­ic C++ tem­plate inter­face.

Now I know tem­plates are scary and all, but you don’t need to know them to use my code. Here’s a sim­ple exam­ple of a physics sim­u­la­tor loop for solv­ing 1D New­to­ni­an physics with con­stant 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 pret­ty straight­for­ward: every physics iter­a­tion, make the fol­low­ing call for every par­ti­cle i: integrator.evaluate(poss[i], vels[i], constAccel, t, dt, nextPoss[i], nextVels[i]). Thanks to type infer­ence, the actu­al tem­plates are hid­den from you. But what are those argu­ments?

  • poss[i] and vels[i] are the particle’s state (posi­tion and veloc­i­ty)
  • constAccel is a func­tion that returns the accel­er­a­tion for that par­ti­cle (in this case a con­stant)
  • t and dt are the sys­tem time and timestep
  • nextPoss[i] and nextVels[i] are the out­puts to the call—they get writ­ten with the next state of the par­ti­cle

Of course, those argu­ments don’t have to be lim­it­ed 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 sys­tem posi­tions and veloc­i­ties be std::complexs (why would you ever do this?), the time and timestep be double, and replaced the con­stant 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 con­stant accel­er­a­tion with a drag accel­er­a­tion that’s lin­ear to the veloc­i­ty.

But look how the integrator.evaluate() call and the sim­u­la­tion loop look almost the same as before: that’s the beau­ty of tem­plates2. You can change around your vec­tor & time types and still keep the same loop code.

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

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? Any­ways, all the inte­gra­tors are avail­able in a sin­gle head­er 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 han­dle 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 read­er. 🙂

Also, weath­ered C++ coders may notice that I spec­i­fy the inte­gra­tors’ inter­face using a CRTP base class. There’s no com­pelling soft­ware engi­neer­ing rea­son for this; I actu­al­ly find CRTP sta­t­ic poly­mor­phism to be a most­ly use­less con­struct, as you can’t swap around your derived class­es 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 beau­ty of tem­plates” is a con­tra­dic­tion in terms []