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