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::complex`s (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 []