Templated numerical integrators in C++

2012/3/11

I find myself writing numerical integrators for physics solvers often enough^* that I wished I didn’t have to keep rewriting them every time I used a different linear algebra library or changed the dimensions what I was solving.

This is exactly what C++ templates are great for—writing code that works with multiple types with similar interfaces. So I whipped up Euler, midpoint, and Runge-Kutta order 4 (RK4) integrators with a generic C++ template interface.

Now I know templates are scary and all, but you don’t need to know them to use my code. Here’s a simple example of a physics simulator loop for solving 1D Newtonian physics with constant acceleration:

#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 straightforward: every physics iteration, make the following call for every particle i: integrator.evaluate(poss[i], vels[i], constAccel, t, dt, nextPoss[i], nextVels[i]). Thanks to type inference, the actual templates are hidden from you. But what are those arguments?

poss[i] and vels[i] are the particle’s state (position and velocity)
constAccel is a function that returns the acceleration for that particle (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 written with the next state of the particle

Of course, those arguments don’t have to be limited to those in that example. For example, 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 positions and velocities be std::complex<double>s (why would you ever do this?), the time and timestep be double, and replaced the constant acceleration function with a neat functor that can maintain state. In this case, the acceleration functor applies a constant acceleration with a drag acceleration that’s linear to the velocity.

But look how the integrator.evaluate() call and the simulation loop look almost the same as before: that’s the beauty of templates^†. You can change around your vector & time types and still keep the same loop code.

You wanna see something wild? Let’s stick in Eigen Vector3f and make our acceleration function 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 integrators are available 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.

Integrators.h

You may have noticed that the integrators don’t handle mass at all—that’s why you pass it an acceleration evaluation function rather than a force evaluation function, as it has no knowledge of mass. Extending the integrators to use mass is a trivial exercise left up to the reader. 🙂

Also, weathered C++ coders may notice that I specify the integrators’ interface using a CRTP base class. There’s no compelling software engineering reason for this; I actually find CRTP static polymorphism to be a mostly useless construct, as you can’t swap around your derived classes without recompiling. But it does make for a good interface spec.

Templated numerical integrators in C++

Every Post by Year