Runge-Kutta (RK4) integration for game physics

Gaffer on Games has a great article about using RK4 integration for better game physics. The implementation is straightforward, but the math behind it confuses me. I understand derivatives and integrals on a conceptual level, but haven't manipulated equations in a long while.

Here's the brunt of Gaffer's implementation:

void integrate(State &state, float t, float dt)
{
     Derivative a = evaluate(state, t, 0.0f, Derivative());
     Derivative b = evaluate(state, t+dt*0.5f, dt*0.5f, a);
     Derivative c = evaluate(state, t+dt*0.5f, dt*0.5f, b);
     Derivative d = evaluate(state, t+dt, dt, c);

     const float dxdt = 1.0f/6.0f * (a.dx + 2.0f*(b.dx + c.dx) + d.dx);
     const float dvdt = 1.0f/6.0f * (a.dv + 2.0f*(b.dv + c.dv) + d.dv)

     state.x = state.x + dxdt * dt;
     state.v = state.v + dvdt * dt;
}

Can anybody explain in simple terms how RK4 works? Specifically, why are we averaging the derivatives at 0.0f, 0.5f, 0.5f, and 1.0f? How is averaging derivatives up to the 4th order different from doing a simple euler integration with a smaller timestep?

---

After reading the accepted answer below, and several other articles, I have a grasp on how RK4 works. To answer my own questions:

Can anybody explain in simple terms how RK4 works?

> RK4 takes advantage of the fact that > we can get a much better approximation > of a function if we use its > higher-order derivatives rather than > just the first or second derivative. > That's why the Taylor series > converges much faster than Euler > approximations. (take a look at the > animation on the right side of that > page)

Specifically, why are we averaging the derivatives at 0.0f, 0.5f, 0.5f, and 1.0f?

> The Runge-Kutta method is an > approximation of a function that > samples derivatives of several points > within a timestep, unlike the Taylor > series which only samples derivatives > of a single point. After sampling > these derivatives we need to know how > to weigh each sample to get the > closest approximation possible. An > easy way to do this is to pick > constants that coincide with the > Taylor series, which is how the > constants of a Runge-Kutta equation > are determined. > > This article made it clearer for > me. Notice how (15) is the Taylor > series expansion while (17) is the > Runge-Kutta derivation.

How is averaging derivatives up to the 4th order different from doing a simple euler integration with a smaller timestep?

> Mathematically, it converges much > faster than doing many Euler > approximations. Of course, with enough > Euler approximations we can gain equal > accuracy to RK4, but the computational > power needed doesn't justify using > Euler.

This may be a bit oversimplified so far as actual math, but meant as an intuitive guide to Runge Kutta integration.

Given some quantity at some time t1, we want to know the quantity at another time t2. With a first-order differential equation, we can know the rate of change of that quantity at t1. There is nothing else we can know for sure; the rest is guessing.

Euler integration is the simplest way to guess: linearly extrapolate from t1 to t2, using the precisely known rate of change at t1. This usually gives a bad answer. If t2 is far from t1, this linear extrapolation will fail to match any curvature in the ideal answer. If we take many small steps from t1 to t2, we'll have the problem of subtraction of similar values. Roundoff errors will ruin the result.

So we refine our guess. One way is to go ahead and do this linear extrapolation anyway, then hoping it's not too far off from truth, use the differential equation to compute an estimate of the rate of change at t2. This, averaged with the (accurate) rate of change at t1, better represents the typical slope of the true answer between t1 and t2. We use this to make a fresh linear extrapolation from to t1 to t2. It's not obvious if we should take the simple average, or give more weight to the rate at t1, without doing the math to estimate errors, but there is a choice here. In any case, it's a better answer than Euler gives.

Perhaps better, make our initial linear extrapolation to a point in time midway between t1 and t2, and use the differential equation to compute the rate of change there. This gives roughly as good an answer as the average just described. Then use this for a linear extrapolation from t1 to t2, since our purpose it to find the quantity at t2. This is the midpoint algorithm.

You can imagine using the mid-point estimate of the rate of change to make another linear extrapolation of the quantity from t1 to the midpoint. With the differential equation we get an better estimate of the slope there. Using this, we end by extrapolating from t1 all the way to t2 where we want an answer. This is the Runge Kutta algorithm.

Could we do a third extrapolation to the midpoint? Sure, it's not illegal, but detailed analysis shows diminishing improvement, such that other sources of error dominate the final result.

Runge Kutta applies the differential equation to the intial point t1, twice to the midpoint, and once at the final point t2. The in-between points are a matter of choice. It is possible to use other points between t1 and t2 for making those improved estimates of the slope. For example, we could use t1, a point one third the way toward t2, another 2/3 the way toward t2, and at t2. The weights for the average of the four derivatives will be different. In practice this doesn't really help, but might have a place in testing since it ought to give the same answer but will provide a different set of round off errors.

RK4 in the simplest sense is making a approximation function that that is based on 4 derivatives and point for each time step: Your initial condition at starting point A, a first approximated slope B based on data point A at your time step/2 and the slope from A, a third approximation C , which has an correction value for the slope at B to reflect the shape changes of your function, and finally a final slope based on the corrected slope at point C.

So basically this method lets you calculate using a starting point, an averaged midpoint which has corrections built into both parts to adjust for the shape, and a doubly corrected endpoint. This makes the effective contribution from each data point 1/6 1/3 1/3 and 1/6, so most of your answer is based on your corrections for the shape of your function.

It turns out that the order of an RK approximation (Euler is considered an RK1) corresponds to how its accuracy scales with smaller time steps.

The relationship between RK1 approximations is linear, so for 10 times the precision you get roughly 10 times better convergence.

For RK4, 10 times the precision yields you about 10^4 times better convergence. So while your calcuation time increases linearly in RK4, it increases your accuracy polynomially.

As to your question why: I recall once writing a cloth simulator where the cloth was a series of springs interconnected at nodes. In the simulator, the force exerted by the spring is proportional to how far the spring is stretched. The force causes acceleration at the node, which causes velocity which moves the node which stretches the spring. There are two integrals (integrating acceleration to get velocity, and integrating velocity to get position) and if they are inaccurate, the errors snowball: Too much acceleration causes too much velocity which causes too much stretch which causes even more acceleration, making the whole system unstable.

It is difficult to explain without graphics, but I'll try: Say you have f(t), where f(0) = 10, f(1) = 20, and f(2) = 30.

A proper integration of f(t) over the interval 0 < t < 1 would give you the surface under the graph of f(t) over that interval.

The rectangle rule integration approximates that surface with a rectangle where the breadth is the delta in time and the length is the new value of f(t), so in the interval 0 < t < 1 , it will yield 20 * 1 = 20, and in the next interval 1

Now if you were to plot these points and draw a line through them you'll see that it is actually triangular, with a surface of 30 (units), and therefore the Euler integration is inadequate.

To get a more accurate estimation of the surface (integral) you can take smaller intervals of t, evaluating at for example f(0), f(0.5), f(1), f(1.5) and f(2).

If you're still following me, the RK4 method is then simply a way of estimating values of f(t) for t0 < t < t0+dt invented by people smarter than myself for getting accurate estimates of the integral.

(but as others have said, read the Wikipedia article for a more detailed explanation. RK4 is in the category of numerical integration)