Limegarden.net Personal site of Wouter Lindenhof

Runge-Kutta 4th order

Finally I got RK4 working. I never had highly sophisticated math and more then half of the books you come across start easy, telling you how to program (the part I do know) and then suddenly they shift in high gear and expect you to suddenly understand all the math that is throw at you.
It seems that most of these physics books are written for physics students and not for programmers. But also wikipedia assumes that I'm a math guru, which I'm not. In fact I like things that don't require me to think, you don't want to know how an engine exactly works if you only need to use it.

So I decided to post my own C++ implementation of RK4, which you are free to steal and use, but you can't hold me liable. I would like to hear if it worked for you or not, but don't feel obligated. My own version is a bit advancer than the one you will see, but that is only because I need to have drag, impulse, etc. I'm not going to explain as I simply don't understand every step, it is more stable than euler (read: it takes longer to break).

If you are really intersted in how RK4 works I suggest that you contact your math teacher and keep asking until you know the exact ins and outs of it.

struct Particle
{
	Vec3f Position;	// Vec3f is a vector (3d position) which should initilizes with all zero
	Vec3f Velocity;
	Vec3f Acceleration;
};

void EulerIntegration(Particle p)
{
	p.Position += p.Velocity;
	p.Velocity += p.Acceleration;
}

// We will use this function later, it's just to make it c++ correct
// DT is DeltaTime, just use 1.0f if you don't know what I mean with that.
void RK4Eval(const Particle p, float DT, const Vec3f d[2], Vec3f* output);

// DT can be 1.0f if you don't know what I mean with that
void RungeKuttaIntegrator(float DT, Particle& p)
{
	Vec3f k0[2]; Vec3f k1[2]; Vec3f k2[2]; Vec3f k3[2]; Vec3f k4[2]; // remember all should be initiliazed zero
	RK4Eval(p,    0.0f, k0, k1);
	RK4Eval(p, DT*0.5f, k1, k2);
	RK4Eval(p, DT*0.5f, k2, k3);
	RK4Eval(p, DT     , k3, k4);

	const Vec3f dxdt = (1.0f/6.0f)*(k1[0] + 2.0f*(k2[0]+k3[0]) + k4[0] );
	const Vec3f dvdt = (1.0f/6.0f)*(k1[1] + 2.0f*(k2[1]+k3[1]) + k4[1] );

	p.Position = p.Position + dxdt * DT;
	p.Velocity = p.Velocity + dvdt * DT;
}

inline void RK4Eval(const Particle p, float DT, const Vec3f d[2], Vec3f* output)
{
	const int x = 0; const int v = 1;	// otherwise I keep thinking arrays
	Vec3f state[2];
	state[x] = p->Position + d[x]*DT;
	state[v] = p->Velocity + d[v]*DT;

	output[x] = state[v];
	output[v] = p->Acceleration;
}

Reading material

• Game Physics - A book recommend by my school, which I don't recommend, as it's a lot of math while the code is simple.

• Game Physics Engine Development - A book I do recommend, it's quite easy going and it brings you to a goal namely a physics engine. So far no complains.
• gaffer.org - It told me the step I was missing (the derivatives)