## Introduction

Hi, I’m Glenn Fiedler and welcome to the first article in my series on Game Physics.

If you have ever wondered how the physics simulation in a computer game works then this series of articles will explain it for you. I assume you are proficient with C++ and have a basic grasp of physics and mathematics. Nothing else will be required if you pay attention and study the example source code.

A game physics simulation works by making many small predictions based on the laws of physics. These predictions are actually quite simple, and basically boil down to something like “given that the object is here, and is traveling this fast, in a short amount of time from now should be over there”. We perform these predictions using a mathematical technique called integration.

Exactly how to implement this integration is the subject of this article.

## Integrating the Equations of Motion

You should remember from high school or university physics that force equals mass times acceleration. We can switch this equation around to see that acceleration equals force divided by mass. This makes sense because the more an object weighs, the less acceleration it will receive from the same amount of force.

f= maa=f/m

Since acceleration is the rate of change in velocity over time, and acceleration is force divided by mass, we can say that force divided by mass is the rate of change in velocity:

dv/dt =a=F/m

Similarly, velocity is the rate of change in position over time:

dx/dt =v

This means that if we know the current position and velocity of an object, and the forces that will be applied to it, we can integrate to find its position and velocity at some point in the future.

## Numerical Integration

For those who have not formally studied differential equations at university, take heart for you are in just as good a position as those who have! This is because we are not going to solve the differential equations as you would normally do in first year mathematics, instead we are going to numerically integrate to find the solution.

Here is how numerical integration works. First, start at a certain initial position and velocity, then take a small step forward in time to find the velocity and position at the next time value. Do this repeatedly to go forward in time in small increments, each time taking the results of the previous integration as the starting point for the next.

For simplicity, lets assume we are stepping forward one second at a time. First we start at our initial velocity and position, then work out how much the velocity and position change over the course of one second. Then, we add this change to the velocity and position and repeat the process. Starting with the velocity values we now have at 1 second, we work out how much the values will now change over the next second to find the values at 2 seconds and so on for 3, 4, 5, 6… etc.

But how do we find the amount of change in velocity and position each step? The answer lies in the equations of motion presented in the paragraph above. We work out how much acceleration is being applied by dividing force by mass. Since acceleration is literally change in velocity per second, we realize that the change in velocity is acceleration * (time in seconds). Similarly, the change in position is velocity * (time in seconds).

We will call our current time ‘t’ and the change in time between each step ‘dt’ (delta time). Now we can present the equations of motion in a form that should make sense to anyone:

- acceleration = force divided by mass
- change in velocity = acceleration times delta time
- change in position = velocity times delta time

It makes sense because we know if you are traveling 60 kilometers per hour in a car, then in one hour we’ll be exactly 60 kilometers further down the road. Similarly, if a car is accelerating at 10 kilometers per hour per second from a standing start, then in 10 seconds it must be traveling at exactly 100 kilometers per hour.

Lets put this into code. Assuming that our initial position and velocity are zero, we can integrate the equations of motion over time simply by adding the change in position and velocity to the previous values each timestep:

float t = 0; float dt = 1; float velocity = 0; float position = 0; float force = 10; float mass = 1; while ( t <= 10 ) { position = position + velocity * dt; velocity = velocity + ( force / mass ) * dt; t = t + dt; }

Notice how the integrated velocity feeds into the integration of the position next frame. This uses the most recently integrated velocity value to update the position of the car over time so that at each second we know how far the car has travelled down the road and how fast it is traveling.

## Why Euler is not (always) so great

What we just did above is called Euler Integration. Specifically it is a type of integration called explicit euler.

To save you future embarrassment, I must point out now that Euler is pronounced “Oiler” not “yew-ler” as it is the last name of the Swiss mathematician Leonhard Euler who first discovered this technique.

Euler integration is the most basic of numerical integration techniques. It is only 100% accurate if the rate of change is constant over the timestep. For example it is accurate in integrating velocity over time in the example above because acceleration is a constant 10 meters per second per second. But when we integrate velocity to get the position in the example, it is not completely accurate. This is because velocity is not constant, it is changing over time due to the acceleration.

This means that the values for **x**(t) and **v**(t) that the integrator produces over time are different from the true values for the functions of position and velocity over time. In the example above the distance travelled down the road is somewhat inaccurate because of the integrator.

Just how inaccurate is the Euler integrator?

From basic physics we know:

s = ut + 0.5at^2

That is, given a constant acceleration a, and an initial velocity u, we can work out how much distance (s) is travelled over the time t.

In our case, u=0 because the car is at a standing start, so:

s = 0.5at^2 s = 0.5(10)(10)^2 s = 0.5(10)(100) s = 500 meters

If we use the Euler integrator above instead of solving exactly we get:

t=0: position = 0 velocity = 0 t=1: position = 0 velocity = 10 t=2: position = 10 velocity = 20 t=3: position = 30 velocity = 30 t=4: position = 60 velocity = 40 t=5: position = 100 velocity = 50 t=6: position = 150 velocity = 60 t=7: position = 210 velocity = 70 t=8: position = 280 velocity = 80 t=9: position = 360 velocity = 90 t=10: position = 450 velocity = 100

So when using the Euler integrator we have an inaccuracy of 50 meters for this simple calculation after only 10 seconds! The bad news is that this error keeps getting larger as time goes on.

One solution to this problem is to decrease the time step. For example we can reduce the time step to dt = 1/100th of second each Euler integration instead of one second. Doing this will get us a result much closer to the exact value.

However, no matter how much you reduce your timestep, the simple truth is that the error will always be there and that it will keep increasing over time. Given that this is an extremely simple simulation, imagine something as simple as adding a torque curve to the car, or adding gears. Suddenly the car’s acceleration is not a constant, it changes over time. Now there is error when integrating the velocity, and these errors are magnified even further when integrating the position.

## Many different integration methods exist

You can see that explicit euler has a large amount of error in the case above, but it’s a gross simplification to say that it is terrible or that it should never be used. Just be aware of its shortcomings. There are also several other Euler method variants with different properties and different regions of stability.

For example, semi-implicit euler (symplectic euler) is a method that changes only the point of time where velocity is sampled, eg. it adds acceleration to velocity first, and then updates position from that updated velocity, vs. the other way around. This gives the integrator different properties and it is more stable and energy conserving across a wider range of time steps, although it still has integration error order 1.

On a different note is implicit euler (backwards euler). This method is better for simulating stiff equations (eg. stiff springs) that break down and become unstable with other methods. It is more involved and requires solving a system of equations per-timestep vs. just substituting values and adding velocity * dt to position and so on. It feels quite different to implement vs. the explicit and semi-implicit euler methods.

An excellent option for greater accuracy and less memory when simulating a large number of particles is verlet integration. In this case, the integrator does not require storing velocity per-particle, as it derives it from the difference in the previous two frames position values. This makes collision detection and position fix-up (eg. pushing particles out of collision) really quick and easy to implement. Here is an excellent example.

In this article I’m going to introduce the “RK4″ integrator, eg. Runge Kutta order 4, so named for the two german mathematicians that discovered it: Carl Runge and Martin Kutta. This is not to suggest that this is automatically “the best” integrator for all applications (it’s not), or that you should always use this integrator over other methods listed above (you shouldn’t).

But it is one of the most accurate general purpose integration techniques available. Read on for the details of how it works.

## Implementing RK4 is not as hard as it looks

I could present to you RK4 in form of general mathematical equations which you could use to integrate some function using its derivative f(x,t) as is usually done, but seeing as my target audience for this article are programmers, not mathematicians, I will explain using code instead. From my point of view, processes like this are best described in code anyway.

So before we go any further I will define the state of an object as a struct in C++ so that we have both position and velocity values conveniently stored in one place.

struct State { float x; // position float v; // velocity };

We’ll also need a struct to store the derivatives of the state values so we can easily pass them around to functions. The derivatives we store are velocity and acceleration:

struct Derivative { float dx; // dx/dt = velocity float dv; // dv/dt = acceleration };

The final piece of the puzzle that we need to implement the RK4 integrator is a function that can advance the physics state ahead from t to t+dt using one set of derivatives, then once there, recalculate the derivatives using this new state. This routine is the heart of the RK4 integrator and when implemented in C++ it looks like this:

Derivative evaluate( const State &initial, float t, float dt, const Derivative &d ) { State state; state.x = initial.x + d.dx*dt; state.v = initial.v + d.dv*dt; Derivative output; output.dx = state.v; output.dv = acceleration(state, t+dt); return output; }

It is absolutely critical that you understand what this method is doing. First it takes the current state of the object (its position and velocity) and advances it ahead dt seconds using an Euler Integration step with the derivatives that were passed in (velocity and acceleration). Once this new position and velocity are calculated, it calculates new derivatives at this point in time using the integrated state. These derivatives will be different to the derivatives that were initially passed in to the method if the derivatives are not constant over the timestep.

In order to calculate the derivatives it copies the current state velocity into the derivative struct (this is the trick to doing position and velocity integration simultaneously) then it calls the acceleration function to calculate the acceleration for the current state at the time t+dt. The acceleration function is what drives the entire simulation and in the example source code for this article I define it as follows:

float acceleration( const State &state, float t ) { const float k = 10; const float b = 1; return -k * state.x - b*state.v; }

This method calculates a spring and damper force and returns it as the acceleration assuming unit mass. What you write here of course is completely simulation dependent, but it is crucial that you structure your simulation in such a way that it can calculate the acceleration or force derivatives completely from inside this method given the current state and time, otherwise your simulation cannot work with the RK4 integrator.

Finally we get to the integration routine itself which integrates the state ahead from t to t+dt using RK4:

void integrate( State &state, float t, float dt ) { Derivative a,b,c,d; a = evaluate( state, t, 0.0f, Derivative() ); b = evaluate( state, t, dt*0.5f, a ); c = evaluate( state, t, dt*0.5f, b ); d = evaluate( state, t, dt, c ); float dxdt = 1.0f / 6.0f * ( a.dx + 2.0f*(b.dx + c.dx) + d.dx ); 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; }

Notice the multiple calls to evaluate that make up this routine. RK4 samples derivatives four times to detect curvature instead of just once in Euler integration. The important thing to understand is how it actually does this sampling.

If you look at the code above it should be clear what is going on. Notice how it uses the previous derivative when calculating the next one. When calculating derivative b it steps ahead from t to t+dt*0.5 using derivative a, then to calculate c it steps ahead using derivative b, and finally d is calculated by stepping ahead with c. This feedback of the current derivative into the calculation of the next one is what gives the RK4 integrator its accuracy.

Once the four derivative samples have been evaluated, the best overall derivative is calculated using a weighted sum of the derivatives which comes from the Taylor Series expansion. This single value can then be used to advance the position and velocity over dt as we did before in the Euler integrator.

Notice that even when using a complicated integrator such as RK4, it all boils down into something = something + change in something * change in time. This is because differentiation and integration are fundamentally linear operations. For now we are just integrating single values, but rest assured that it all ends up like this when integrating vectors, or even quaternions and matrices for rotational dynamics.

## Conclusion

I have demonstrated how to implement an RK4 integrator for a basic physics simulation. At this point if you are serious about learning game physics you should study the example source code that comes with this article and play around with it.

Here are some ideas for study:

- Switch from integrating velocity directly from acceleration to integrating momentum from force instead (the derivative of momentum is force). You will need to add “mass” and “inverseMass” to the State struct and I recommend adding a method called “recalculate” which updates velocity = momentum * inverseMass whenever it is called. Every time you modify the momentum value you need to recalculate the velocity. You should also rename the acceleration method to “force”. Why do this? Later on when working with rotational dynamics you will need to work with angular momentum directly instead of angular velocity, so you might as well start now.
- Try modifying the integrate method to implement an Euler integrator. Compare the results of the simulation against the RK4 integrator. How much can you increase the spring constant k before the simulation explodes with Euler? How large can you make it with RK4?
- Extend position, velocity and force to 3D quantities using vectors. If you use your intuition you should easily be able to extend the RK4 integrator to do this.

**Next:** Fix Your Timestep!

If you enjoyed this article please consider making a small donation. __Donations encourage me to write more articles!__

I’m implementing a physics engine using multistep integration (5-step Adams-Bashforth-Moulton); I currently use symplectic Euler to compute starting values (A-B-M is not self-starting), and am considering switching to RK4.

I have two questions:

1) Assuming RK4 more accurately calculates starting values, is there likely to be a significant increase in overall stability or accuracy as a result?

2) When implementing RK4, is there an advantage to taking a symplectic approach when evaluating the intermediate values?

There is not likely to be a massive reason to use RK4 for your application. Most physics engines don’t need to use Euler for rigid bodies under constant acceleration symplectic Euler is good enough.

cheers

I cannot thank you enough for a real example in code. That alone helped cement the concept for me.

No problem! Glad you found it useful — cheers

Heyo, for the first example, Euler integration, doesn’t using average velocity for that time step give more accurate results than using either velocity from the start or end of the frame. This is assuming that acceleration is constant, but if velocity = 0 at the start of the frame and 1 at the end, then the integral for that has the same area under the curve as a steady velocity of 0.5. Since area under the velocity curve, when integrated maps directly to increase in position, this small improvement can reduce the error inherent in that form of modelling.

It does. This is equivalent to the ballistic equation of motion: s = ut + 1/2at^2

But this is only accurate for constant acceleration. As soon as you have an acceleration that is a function of distance between two points, eg. springs — it is no longer 100% accurate.

I have two questions: 1) dxdt is velocity, so why don’t we set state.v to it and ignore integrating for velocity? Or, better yet, integrate for velocity and use that result to get position? 2) In a model rocket, mass is a function of time. So, if I’m writing a rocket simulation, I can’t integrate Force, can I?

Yes, you can do that and avoid storing velocity explicitly. This is known as a verlet integrator. It’s really good for particle systems. More here: http://web.archive.org/web/20080410171619/http://www.teknikus.dk/tj/gdc2001.htm

Hello, Thanks for the article !; I have a question, how can i add friction in these formulas ?

https://en.wikipedia.org/wiki/Friction

Hi Glenn,

Thank s a lot for the wonderful article.

The function acceleration

float acceleration( const State &state, float t )

{

const float k = 10;

const float b = 1;

return -k * state.x – b*state.v;

}

is the dynamic model of a damped-spring.

Is this function a “dummy” function to see how RK4 works?

The reason of my question is. If I write a program on my own, where the model of a quadrotor is implemented,. should I put in that function the acceleration of my model?

I hope my question is clear

Yes it’s just a dummy function that you replace with your own model returning the force at a given object state and (optional) time.

A question about this part…

float acceleration( const State &state, float t )

{

const float k = 10;

const float b = 1;

return -k * state.x – b*state.v;

}

why is “float t” there? It doesn’t do anything !

Why “b*state.v” when b will always be 1 ?

float t is not used in this case, but in other cases it could be (but I would recommend a double). For example, maybe a force f=sin(t)?

Next, for the b this is the damping amount. You can change it from 1 to other values, higher values mean more drag.

cheers

There is something to be said about complexity for sake of complexity. This is generally true with the madness of religious calculus.

god bless you for the oiler tip. It would of been awkward.

When you demonstrated how Euler’s method produces an error of 50m, how come the position is 0 when t=1? As the velocity changes to 10 as t=1, shouldn’t the position increase as well as t=1?

Because initial velocity is zero at t=0, and then you add that zero to position(0) to get position(1) of zero. You could of course, add acceleration to velocity first, and then add that to position to get position(1) of 10 if you wish. Work through it and you’ll see it still has error. Interestingly, in the case of constant acceleration if you average the two velocities you’ll get a 100% accurate result, eg. trapezoidal rule, or closed form s = ut+1/2at^2. cheers

clarify this for me please

Derivative a = evaluate(state, t, 0.0f, Derivative());

Derivative b = evaluate(state, t, dt*0.5f, a);

Derivative c = evaluate(state, t, dt*0.5f, b);

Derivative d = evaluate(state, t, dt, c);

Looks like “a” derivative always will be zero. Not just because of undefined

members of struct “Derivative”, but because of “dt”(which is 0.0f by now) is using here:

state.x = initial.x + d.dx*dt;

state.v = initial.v + d.dv*dt;

so the “state” wont ever change(and considering “initial” members all 0 at start point – state.x and .v will be 0 as well). Then we use that “a”(output.dx & .dv = 0) for evaluating “b” – derivatives zero again, and so on.

Apologize in advance if its such silly question.

Nope. Derivative a is the derivative at the point in time t (t+0). Derivative b is the evaluation at t+dt/2 using derivative a. Derivative c is the derivative estimate at t+dt/2 using derivative b. d is the derivative at t+dt using the derivative c. The linear combination of these together provides the RK4 integration method.

Hi Glenn,

Question about this little bit:

“it is crucial that you structure your simulation in such a way that it can calculate the acceleration or force derivatives completely from inside this method given the current state and time, otherwise your simulation cannot work with the RK4 integrator.”

I have a function that calculates the force direction of user input. It basically takes the length and direction of the mouse movement and applies a constant.

I then use this in the acceleration function after the spring acceleration is computed by doing this:

accel.x += inputForce.x;

…

Is this incorrect? It seems to produce good results

If you are happy with the result then it is perfectly fine.

Actually I should’ve asked my question in your timestep article…

http://gafferongames.com/game-physics/integration-basics/#comment-241186

And after I checked the comments there I see you’ve already answered it June 2007:

> interpolation is between two known values, so its “perfect” – just lagged by up to one [physics] frame

Looking at the article’s sample code:

while (accumulator>=dt)

{

accumulator -= dt;

previous = current;

integrate(current, t, dt);

t += dt;

}

State state = interpolate(previous, current, accumulator/dt);

And interpolate function:

State interpolate(const State &previous, const State ¤t, float alpha)

{

State state;

state.x = current.x*alpha + previous.x*(1-alpha);

state.v = current.v*alpha + previous.v*(1-alpha);

return state;

}

You’re not really rendering an interpolation of future state, right? Instead what you’re doing is you render an interpolation between between previous state and current state (which is rounded down to the most recent physics update).

Suppose we are on physics update 10.3 (meaning 30% of the way from physics update 10 to 11). Your code interpolates between 30% previous + 70% current. Previous state is from update 9 and current state is from update 10, so you are interpolating the state between 9 and 10. 30%*9 + 70%*10 means your interpolate() gives you an approximation for state at update 9.3.

In other words, you are rendering one physics frame in the past. If we are on physics frame 10.3, you interpolaate between previous (9) and current (10) to get 9.3. Did I miss something? So you’re rendering one physics frame in the past. Correct?

I notice in this other article ( http://www.koonsolo.com/news/dewitters-gameloop/ ), he interpolates one frame into the future. He does this by using a constant speed between physics updates (ie between update 10.0 and 10.3). So we do a full physics update for 9 and 10 and 11. But for 10.3, we do current (10) + a low-computation-cost approximation. This low-computation-cost approximation is to just assume a constant velocity.

So the gafferongames.com article gives a more accurate interpolation but it’s one frame behind (eg frame 9.3). The koonsolo.com article uses a low-computation-cost approximation of the future frame (eg frame 10.3) such as assuming a constant velocity. This is only for sub-frames; once we get to frame 11 we’ll do the full update computation.

Any errors in my understanding?

Hey,

great article. I’m just wondering how you would handle a wrapping t value?

That’s a tricky one. I don’t know! The problem here is that you have usually different functions of t with different periods so there is no safe place to wrap

Hello,

Interesting article, many thanks.

However, when I change the dt in the main function of your sample code the results a very strange. For instance, when

dt = .1f it takes about 23.0 secs for the spring to stop

dt = .2f and dt = .3f….. also 23 secs which is very nice

dt = .4f….. 22.4 secs also nice

dt = .5f….. 20.0 secs

dt = .6f……15.6 seconds

So, it means that if I run your code on a slow computer the simulation gets out of hand.

I’d expect the spring to stop always at 23 seconds.

Am I wrong?

Very interesting!

Hello,

over here https://sites.google.com/site/mcasualsdazscripts3/physicsrk4_001

my Daz Script ( Qt script ) implementation running in Daz Studio (all versions)

i’ll probably make an HTML5 version if google-sites lets me embed it

There’s also 2 verlet based experiments

particles https://sites.google.com/site/mcasualsdazscripts/particle-physics

cloth sim https://sites.google.com/site/mcasualsdazscripts/cloth-physics-experiment-kit

Has anyone tried to implement this in java for android devices?

I am having lots of fun handling memory allocations.

I don’t know why but I lol’d at your post 😀

This is a really useful article. I’ve gone over it like five times though I’ve noticed that when it comes to RK4, pretty much everyone refers to this article and copies and pastes the code you have here. That’s fine, but I don’t want to do that, I want to understand why you made the decisions you did here. As someone else mentioned above, your implementation of RK4 doesn’t exactly match the traditional definition of RK4 on Wikipedia and elsewhere, which goes something like this:

k1 = v(t, x) * dt

k2 = v(t + 0.5*dt, x + 0.5*k1) * dt

k3 = v(t + 0.5*dt, x + 0.5*k2) * dt

k4 = v(t + dt, x + k3) * dt

x_next = x + ((k1 + 2*k2 + 2*k3 + k4) / 6)

In plain English, we want to calculate four offsets for x:

* one using the velocity calculated using the current values of t and x (k1)

* one using the velocity calculated from advancing t by a half-timestep and x by half the k1 offset (k2)

* one using the velocity calculated from advancing t by a half-timestep and x by half the k2 offset (k3)

* and finally one using the velocity calculated from advancing t by a full timestep and x by the full k3 offset (k4)

* then average those offsets to get a final offset we can apply to x.

The first difference in your code, as someone else has said, is that instead of simply calculating position and velocity, you’re doing position, velocity, and acceleration.

Secondly, your process goes rather more like this:

* at t = 0: get an initial state of acc, vel, and pos (k1)

* at t = 0.5: integrate acc from k1 to get vel, vel from k1 to get pos, and then re-calculate acc (k2)

* at t = 0.5: integrate acc from k2 to get vel, vel from k2 to get pos, and then re-calculate acc (k3)

* at t = 1: integrate acc from k3 to get vel, vel from k3 to get pos, and then re-calculate acc (k4)

* then average the acc’s and vel’s and integrate one final time to get final vel and pos

I can see how this is similar to RK4 in that we are sampling four points (0, 1/2, 1/2, 1) and averaging them together, but I’m trying to figure out how you took the RK4 concept and translated it into code. For instance, this line:

k2 = v(t + 0.5*dt, x0 + 0.5*k1) * dt

is translated as this:

k2 = evaluate(state, t, dt*0.5, k1)

which is the same as saying

x = x0 + k1.dx * 0.5*dt

v = v0 + k1.dv * 0.5*dt

k2.dx = v

k2.dv = acc(x, v, t + 0.5*dt)

The first thing I notice is that you’re calculating slope not offset, so we chop off the dt. Since k1 is now a slope we have to add dt after k1 when we advance x0:

k1 = (a slope, i.e., dx)

k2 = v(t + 0.5*dt, x0 + 0.5*(k1*dt))

If we flip the 0.5 and k1 we get something resembling Euler:

k1 = (a slope, i.e., dx)

k2 = v(t + 0.5*dt, x0 + k1 * 0.5*dt)

Second you also have velocity so we can change that to:

k1 = {dx, dv}

k2 = a(t + 0.5*dt, x0 + k1.dx * 0.5*dt, v0 + k1.dv * 0.5*dt)

And now that’s basically identical to what you’re doing.

Is this kind of the general idea or am I off base here?

I think so. It’s been so long since I wrote this article I’ve basically forgotten

According to Erin Catto’s latest GDC presentation, it is pretty much stupid to use springs as they are difficult to tune in practice. He explains an alternative based on the harmonic oscillator called soft constraints, which are much more stable

I agree

I’m nearly two years late on this, but I recently used this article to write an integrator for a simulation. I saw this comment, rewrote my code, realized what was happening, and unrewrote my code.

The first thing that the evaluate() function does is compute x + delta * k[n]. Then it computes k[n+1], which we pass into the next call. This is RK4, it just looks a little strange because the stepping of the input parameters occurs within the function call due to there being so many input parameters that we want to step. I was writing a function to compute this stepped forward state, and realized that I was copy-pasting the beginning of evaluate().

“If you use Euler then you are a bloody idiot”

I had to laugh!

Awesome article. Made me go back to my Mechanical engineering days.

What you think about multi-step methods ?

http://burtleburtle.net/bob/math/multistep.html

Thanks for the excellent description! I just finished implementing both Euler and RK4 integrator’s for my 2D physics engine. Euler outputs exactly as you show above and RK4 hits the mark perfect! I still want to implement Symplectic Euler and Midpoint methods for my engine and compare how accurate they are.

Good call with midpoint/symplectic euler RK4 is overkill for many things — cheers

Brilliant explanation, thanks!

Hi, thank you for great series of articles.

Just wanted to point out you’re perhaps unduly harsh on Euler. If you swap the position and velocity calculations you get the symplectic or semi-implicit Euler, which would actually be better over long time periods for Hamiltonian systems like springs (without friction). RK4 is dissipative and has quite bad energy and angular momentum conservation properties.

If you need 4th order or higher there are methods which are probably better than RK4 (which is itself about hundred years old and developed for calculations by hand). Cash & Karp, Dormand & Prince, Calvo, Sans-Serna are all mathematicians who have lent their names to ODE solvers that might be helpful. I’ve tested some of these with the Kepler problem and they all performed better than RK4, especially over long time periods.

Possibly this is overkill for a game though

Thanks again for the articles

Agreed 100%. Thanks for the info!

Hey, good job dude. It’s crazy how I found your article on the web and you just happen to sit 2 desks away from me at work! I was trying to figure out how to apply the math from Wikipedia and found your article which was really helpful for making a connection between the theory and a practical implementation. I was inspired to write my own article which shows how to derive the version you have from the Wikipedia math, and in doing so, I was able to show how to implement any of the RK methods. I provided an example for RK4/5 which is useful if you need adaptive time steps to ensure accuracy within a certain tolerance.

You can find my article here:

http://buttersblog.com/runge-kutta

Great article! Good job

great work

velocityX1 = velocityX;

velocityX2 = velocityX + CorrectedAccel * 9.80665 *dt*0.5f;

velocityX3 = velocityX2 + CorrectedAccel * 9.80665 *dt*0.5f;

velocityX4 = velocityX3 + CorrectedAccel * 9.80665 * dt;

velocityX5 = (velocityX1+2.0f*velocityX2+2.0f*velocityX3+velocityX4)/6;

positionX1 = positionX;

positionX2 = positionX + velocityX5 * dt*0.5f;

positionX3 = positionX2 + velocityX5 * dt*0.5f;

positionX4 = positionX3 + velocityX5 * dt;

positionX5 = (positionX1+2.0f*positionX2+2.0f*positionX3+positionX4)/6;

Eh Glenn is this correct?Acceleration in my case is output from the IMU,dt is 0.01

Hi again, i ve noticed smth strange when put RK4 into my program. As i told u i’m working on a fluid sim and seems that i have wall stuck artifact as soon as i integrated with RK4 whereas euler seems to leave room for mistakes and works fine..strange ah?

Hi, and sorry for the silly query

Hi, i am doing a fluid simulation myself and i was wondering what should i put as a t (time).

for example i have timestep (dt) set to 0.006 for example but what about t?

t would be the accumulated time since the simulation started, eg. 0.000, 0.006, 0.012, 0.018 etc…

you probably don’t need that for your fluid simulation though, unless it feeds into your differential equations

yeah i figure this out 5 mins after i posted. i was just new to the method and i didn’t paid so much attention.

Thnx

Can someone help me?I am currently working on a program for an imu.I need to find out how much the imu moved linearly using linear acceleration collected from the imu. I am considering using the rk4 integrator to integrate my linear acceleration twice. However i am not sure where to place the linear acceleration value in the rk4 integrator?Looking at your sample codes,output.dv = acceleration(state, t+dt); i would like to replace that acceleration with my linear acceleration. So how can it be done?

Hi,

I’m a recent Computer Games Technology graduate from the University of Abertay Dundee, Scotland. For my 4th year honours project I conducted an investigation concerning the computational efficiency and the aesthetic quality of using both RK4 and Implicit Euler numerical integration methods for simulation of simple 3D rigid body motion. I developed a small demo application in XNA 4.0, where the user could apply impulses to a box model (making use of a standard box inertia tensor) and evaluate the time taken to compute each time-step summation using either RK4 or Implicit Euler (the user can change the integration methods in real time). What I found was that for simple 3D rigid body motion, there is no real difference between using RK4 or Implicit Euler, as long as a correct time-stepping value (h) is used ([Eberly 2010], [DeLoura et al 2010] and [Bourg 2001] kind of point this out). You can have a look at a copy of my dissertation if you are interested:

http://vladetastojanovic.weebly.com/uploads/1/0/2/3/10232848/vladeta_stojanovic_dissertation_final.pdf

And finally, I would like to thank you for providing these excellent tutorials in such a simple and straightforward manner.

– Vlad

Fantastic info, thanks!

Wait, isn’t this just single paneled Simpson’s quadrature? Runge Kutta methods are usually used to describe implicit integration such that the integral itself is part of the function being integrated (think of PDEs with initial value conditions), but since we can model numerical integration as an ODE without the y parameter, wouldn’t it be more efficient to just not evaluate f(x+step/2) twice?

Also, one nice property of Simpson’s is that, as the 4 in RK4 suggests, it is of degree 4, which means that for panel integration, the absolute error is O(step^4) as long as that step is less than 4, and since the error is proportional to the bound on the 4th derivative, all polynomials of degree 3 or less are perfectly integrated, which is why it works so great for constant acceleration whereas Euler (degree 2, hence only integrates linear functions perfectly) is less than perfect

For a simple jump algorithm, would the acceleration function return a constant?

Yep. Acceleration due to gravity is typically modeled as constant for jumping / ballistic projectiles etc.

Also, do you have any examples of collision detection and response?

http://www.metanetsoftware.com/technique/tutorialA.html

“For example it is accurate in integrating velocity over time in the example above because acceleration is a constant 10 meters per second per second.”

“But when we integrate velocity to get the position in the example, it is not completely accurate”

If you can stomach my ignorance, these statements seem to cancel each other out. I know I’m just not getting it though.

We are integrating twice. First from acceleration to velocity, which is 100% accurate because the acceleration is constant. Then from velocity to position, which is not 100% accurate because the velocity is changing over time (and in fact, in “reality” is changing in between our sample points dt apart…)

cheers

In other words, the first integration is 100% accurate, but the second one (using the result of the first integration, which changes over time…) is not.

Hello and thank you for such an amazing article! I would like to know how to move from Acceleration and Velocity to Force and Momentum. I’ve figured the acceleration() (now called force()) method needs to have mass multiplied, and then state.v and derivative.dvdt becomes state.m and derivative.dmdt, but where does velocity appear? Isn’t it just about changing “state.x = state.x + dxdt*dt;” to state.x = state.x + dxdt*dt/mass;” in integrade()? And maybe “state.x = initial.x + d.dx*dt/mass;” inside evaluate()?

F = ma, therefore acceleration = F/m

p = mv, therefore momentum = velocity * mass

Basically, you integrate force directly to get linear momentum. Just don’t divide force by mass to get acceleration before integration. cheers

And yes, you go from momentum -> velocity, then add the velocity to the previous position. Velocity is still used to update positions. It’s just that momentum is a better quantity to use when dealing with ANGULAR motion, because angular momentum is conserved while angular velocity is not (except in special cases).

Basically, if you are just doing linear motion — there is no benefit to using momentum instead of velocity…

Thank you for a great article! I´m having trouble implementing force and momentum, though. I´ve got it almost working, but trying to rewrite my code I´ve rewritten everything into oblivion and now I´ve got… Nothing.

Anyway, what really seems to mess it up for me is how to integrate moment from force in the “output state”. I just can´t seem to get it right. Could you perhaps shed some light how to do this. If I only get momentum it should be downhill from there, right?

Thank you!

Linear momentum is mass * velocity. Force is the time derivative of linear momentum. So to take force and go to momentum you just integrate force over time. There is an example of this in the source code for this article. cheers

Apparently I got it right on my first try. I just didn´t realize I integrated momentum from force directly (which is what you wrote in your “physics for 3d” article) and then I spent 4 days rewriting my code to understand how momentum would fit in. But I got it right in the end and I even know what I´m doing now.

Thank´s again for a great series of articles! It has really helped me out a lot!

It’s a useful exercise to think it through and rewrite like that. Congrats!

Thanks a lot for a nice article, but I suppose that a brief piece mathematical explanation would still be useful here. Who doesn’t need it, could ignore it, but who understands a bit, would find it really helpful.

You’re right. I’ve been meaning to pass over this article and flesh it out for quite a while, for example I’ve rewritten the “Fix your Timestep” article and this one is next in the queue. cheers

In article:

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);

should be

Derivative a = evaluate(state, t, 0.0f, Derivative());

Derivative b = evaluate(state, t, dt*0.5f, a);

Derivative c = evaluate(state, t, dt*0.5f, b);

Derivative d = evaluate(state, t, dt, c);

since you’re already adding dt inside:

output.dv = acceleration(state, t+dt);

This can cause differences if your acceleration depends on time. I’ve tested Euler integration with 100’000 iterations per frame vs RK4 with 10 iterations per frame using alternating fields – this indeed needs a fix.

Besides that, thanks for great article, it helped a lot.

Yes, your fix looks correct. Thanks!

Works like a charm with a vector implementation. Thanks for making it and for the follow-ups!

Excellent article, helped me out no end.

Thanks a lot for article

Everywhere on web there’s only a lot of empty equations and theory with no easy-to-understand samples provided, thanks, now I fully understand RK4 😉

Hi Glenn,

This article is great, especially having some code to mull and reuse! Thank you.

(Sorry for this long post)

I’m doing a 3-space n-body gravitational sim and am starting to come to the unpleasant conclusion that with RK4 I’m going to have to compute the effects due to all n bodies on each other all at once inside evaluate() and acceleration(), because I can’t isolate the effects of each as I could with Euler or a modifed Euler. If I don’t update all the bodies in each evaluate() call, I’ll see the 2nd, 3rd, and 4th evaluate() calls start to accumulate excess error because the other body positions are stale. Your comment above for the damped spring:

> you integrate both bodies together, ideally you have a single method that returns an array of forces to apply to the points — you cant update each body independently with RK4

seems to support this gloomy conclusion. In a high-gradient field, i.e. when near another massive body, this will become quite important I suspect. I’ll have to fully compute updated state for each of the other n-1 bodies including inside each of the 4 evaluate() calls in integrate(), much as you describe above for the spring case. So I end up doing them all at once.

evaluate() and accelerate() look like becoming much more complex, and computationally quite expensive for large n. integrate() will have some loops, but will not be especially troublesome, I believe.

I think the changes needed, in rough outline, are:

– turn State and Derivative into arrays[n] of structs, each struct containing 6 components (x, y, z, vx, vy, vz), with n being the number of bodies, and have caller pass states[] into integrate() along with t and dt, also making available masses[] and G somehow, maybe as parameters. More on integrate() later.

– accelerate() needs to compute G*mj/r^2 (to get acceleration mi on body bi) for each pair of bodies with masses and return the results in outputs[n]

– evaluate() needs to update all the elements of states[] with the initial.x + d.dx*dt and initial.vx + d.dv*dt computations for each of the 3 axes, initialize outputs[].dx, .dy, and.dz, then invoke acceleration() on states[], which will copy results to outputs[] for evaluate().

– integrate() needs to make the 4 calls to evaluate(), starting with a Derivatives[] initialized to {0.0, 0.0, 0.0, 0.0, 0.0, 0.0} in each array element. Then it can serially (at last, no parallel computations!) compute the dxdt, dydt, dzdt, dvxdt, dvydt, and dvzdt values for each array element and update that element’s state. values.

Does this sound about right?

Yes

really fantastic

Fantastic explanation!

Thank you so much.

Greetings from Germany

A visual explanation of Euler and RK4 http://www.youtube.com/watch?v=LbKKzMag5Rc

Hi,

I hope I didn’t miss any comment already asking this question. I am wondering how you can use RK4 in a normal game. Acceleration will change over time. But it changes unpredictable because the user can change it by f.e. pressing a key. Will the RK4 still work? Do you have to save the last 3 states? Maybe I missunderstand something (perhaps because I’m not native english speaker).

Generally you do 4 RK4 steps for every “real” physics step each frame. This is the bit you are missing! This also makes RK4 pretty expensive in real world situations so you’d want to try a few other integrators first, I think. cheers

Hi

This is graet article. So far for my (usually simple) simulations I used Euler technique of integration thinking RK4 is hard to implement. I have just changed some of my code and it wasn’t so hard. The toughest part was to move forces calculation into one routine as I had just addForce method that was called form many parts of code. I want to be sure i did it right (I have many points in my sim linked with springs – every spring links TWO bodies so it’s force depends on the positions of both bodies). If you could read pseudocode below and tell if it’s right (I divided “evaluate” method to two methods):

1) Do not calculate new state as dt for first “evaluate” is 0 than state is the final state from previous iteration.

2) Calculate derivative A using states from step 1 (using the “accelerate” that iterates all springs attached to the point and adds constant gravity. Each spring contibutes half the force (the second half is for other body)).

3) Calculate new state for each point for t=T+dT/2 dt=dT and using derivative A from step 2

4) Calculate derivative B using states from step 3 with the same “accelerate”

5),6),7),8) etc. (giving C, and D derivatives)

9) Calculate weightes sum of A,B,C and D and calculate the final state

10) show result to screen and next iteration

Is that right?

Sounds about right but for particle simulations with lots of bodies you may find Verlet integration easier to implement. cheers

Ok, I’m not sure if I was clear. My apologizes, I’m not good to comunicate even in my first language. Let me try one more time: talking specifically about rigid body simulation, I’ve heard that some people run a full collision detection routine at middle steps of a RK4 integrator to find a new full set of forces to apply at these steps. Certainly this is much more slow than performing collision detection only one time at the beginning of integration. I would like to know if even ignoring these probable new forces at meddle steps, RK4 is more accurate than Euler. Then, I’ve tried to run your first example in the article using RK4:

http://www.inf.ufsc.br/~sms/testes_programacao/euler_rk4.cc

And what I’ve realized with this test is that RK4 was more accurate than Euler (at least for this case) since the final position is now 500 as it should be.

Then I ask: Is my test correct? Does it make sence?

I’d be glad with any answer.

Thank you!

Well, I’ve just realized that with what I said the integrator could be simplified to something like:

void integrate(State &state, float t, float dt) {

float acc = acceleration(state, t);

float adx = state.v;

float bdx = state.v + acc*dt*0.5f;

float cdx = state.v + acc*dt;

const float dxdt = 1.0f/6.0f * (adx + 4.0f*bdx + cdx);

state.x += dxdt * dt;

state.v += acc * dt;

}

what is that? a kind of RK3 or what? 😛

Thanks!

Is RK4 much more accurate than Euler even if my simulator isn’t able to calculate new values for acceleration at middle steps (t+dt/2)? Thanks for the article!

no, i don’t believe it will be. the sampling of derivatives in the internal timesteps is what gives RK4 the extra accuracy – cheers

I’m trying to convert this code from C++ to C (as well as making it usable for 3-Dimensional movement). My grasp of C++ is rather limited so I was just wondering if someone could clearly explain the purpose of “Derivative()” in the 3rd line of the “integrate” function so that I can properly convert it.

Thanks

derivative is a struct that stores the velocity and acceleration. it’s just two floating point values in a struct to make it easy to pass around. to convert this to 3-dimensional movement you’d want to make it so that position, velocity, acceleration and force are all vector quantities, eg. x,y,z — everything else will work unchanged. cheers

I have the same question as Anonymous but likely just don’t understand your response. In the “integrate” function, why is “Derivative()” being pass to the “evaluate” function to calculate “a”? I undestand that Derivative is a struct but it is seemingly being passed as a function. Is this a C++ specific thing?

That’s just the default constructor. It’s a C++ syntax that says, pass in a “default” version of the derivative object. In this case, I use that to pass in an initial derivative of zero.

I couldn’t help but notice that the “t” variable in evaluate() and accelerate() doesn’t actually do anything. I read back and it seems in 2006, I may have asked the same question (heh – I can’t remember that far back, so it may have very well been a different Alex).

The reason I state this is because I am trying to integrate your ideas in Fix Your Timestep with chipmunk physics, where the “t” variable for simulation doesn’t actually exist.

Take care,

-Alex

The t value does not need to be used, it’s there if you want to drive some time-based forces (eg. sin(t)) but they are not used in this demo (forces are distance and velocity based.)

Perfect, thanks for the extra explaination!

I would like to make a question about collision detection:

Suppose I’m using RK4 to simulate a projectile falling to the ground.

During the simulation, in a particular instant of time, the projectile will be a little above the level of the ground. One timestep later, it will be a little below the level of the ground. So, in this situation, the positions and velocities of the projectile a little before and a little after the collision are known.

To detect the collision, I can check if the projectile is below or at the level of the ground.

But, once the collision is detected, I will have to retrocede the simulation a little so that the projectile is put exactly at the level of the ground. I imagine one way of doing this:

By using the known positions before and after the collision, estimate where the projectile crosses the level of the ground. Then, put the projectile in this position.

The problem here is estimating what the velocity would be at the time of crossing.

I could also literally retrocede the simulation, but with a smaller timestep, which would have to be estimated somehow.

Any ideas?

Thank you in advance.

I’m not really sure of the best approach here. Perhaps you could look into “timewarp” by Brian Mirtich ( http://portal.acm.org/citation.cfm?id=344866 ) where he has some sort of rewinding to point of collision, you could also investigate verlet integration instead of RK4 which remove the need to explicitly patch up velocity post-adjusting position for collisions or subdividing the timestep on a rewind:

http://codeflow.org/entries/2010/sep/01/hard-constraints-easy-solutions/

cheers

Thanks so much for the article and allowing all the smart people here to post questions and ideas.

I am trying to figure out how to use RK4 to simulate the dynamic response of a mass being moved with a flexible linkage and returned by a spring. The real life case is the valve of an engine being moved by the cam, the pushrod is flexible. I have been trying to figure out how to apply RK4 to this for what seems like an eternity and just can’t figure it out.

The inputs are:

cam_displacement //an array of cam actuator input positions in time increment

mass_follower //mass of the first body directly actuated by the cam

mass_valve //mass of the body to be moved by the linkage

stiffness_return_spring //the spring rate of the return spring that acts on the mass_valve

stiffness_pushrod //the spring rate of the linkage between the cam actuator and

//the mass_valve

damping //resistance to motion relative to velocity

The desired return is the dynamic response of the mass_valve at the increments of the input array.

Does anyone know how to set that up using RK4?

First of all, this article is very informative and very useful.

I have a question about how to apply this method to a gravitational N-body simulation. The doubt concerns the fact that all the N bodies affect each other simultaneously. Would the following procedure work, or is there a better way?

Given all the N bodies’ initial positions and velocities in all dimensions:

1. Go through all the N bodies and, for each one, sum the accelerations produced by the other (N – 1) bodies in each dimension. Now, I have the net acceleration of each body at a specific moment.

2. Go through all the N bodies again and, for each one, apply the numerical integration by making use of the net acceleration previously calculated. Now, I’ve found the new position and velocity for each body.

Then, when going back to Step 1, the net acceleration would be recalculated for the new positions.

Thank you in advance.

I believe that option #2 is the most correct because the force between objects depends on their positions — you should advance the whole simulation state forward in one big integration multiple times for the RK4.

May be a good time to bring out the SSE intrinsics!

cheers

Thanks so much for a very informative article. It has been incredibly useful. I think I may have found a slight typo, but please correct me if I’m wrong. I think that:

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);

Should read:

Derivative a = evaluate(state, t, 0.0f, Derivative());

Derivative b = evaluate(state, t, dt*0.5f, a);

Derivative c = evaluate(state, t, dt*0.5f, b);

Derivative d = evaluate(state, t, dt, c);

Thanks again!

Euler can be perfectly fine if your ODE is not stiff. At the end of the day the maximum timestep you can integrate has to be the draw frequency of the interactive simulation. That means that if you’re using RK4 on a non-stiff ODE there reaches a point where the accuracy benefits of RK4 become imperceptible and you begin wasting CPU time.

In my simulations I use a linearizing control scheme for my characters which makes my ODE non-stiff. I can integrate at 100 Hz using forward Euler and achieve stability for most cases. I find that RK2 works best overall when the passive non-linear centrifugal forces kick in.

agreed

Nice article…though euler is not as crappy as u describe it…or have u ever worked in a perturbed environment using RK methods? No? Why? Cause it wont work Thats when u can only use euler….

Nice article. Anyway, I don’t quite get, why the function accelerate() has the parameter float t – it’s not used in the function body at all.

It’s not used but it could be for forces which are functions of time, eg. sinusoidal wobble for a rocket thrust

Hey Glenn,

Thanks for such a informative article. I was looking at the sample code of your next article where you implement a cube rotating spinning around a sin-cos curve and I observed that the way RK4 integrator works is by calculating forces implicitly. Wouldn’t be very tricky to simulate impulse based system like a pool game? In such case, an external input, such as force from other colliding ball or impact with the rail would have to be taken into account.

Yes everything here is force driven. In an impulsed based simulation you’d probably be much better off using a simpler integrator – really the only reason why a higher order integrator is useful here is because I’m attempting to do everything via forces and the integrator lets me have stiffer springs for a given timestep

cheers

Hey Glenn,

Thanks very much for clearing that up. Much appreciated.

That’s already been addressed, and the answer is no.

If you CAN integrate algebraically, would it be faster/valid to just do that? Or Newton-Raphson (my personal favourite numerical solution to a function).

It’s interesting, theoretically you could have a closed form solution which although being entirely accurate (up to floating point accuracy) has so many terms that it’s actually slower than integrating directly

Generally, folks use the closed form solution when they don’t want to worry about drift over time accumulating values in an integrator. It’s not always an option, for example consider that if you have player input affecting the solution, I’m not sure that you could actually get a closed form solution for that — the function of player input is unknown at compile time!

But yeah, for stuff like stiff springs, suspension etc. I’ve heard a lot of people used closed form here for specific cases, just not the “general” case if you get me

cheers

hi everybody, can 4th order runge kutta be used for 3 axis accelerometer? since the values are discrete and we do not know the slopes or step sizes. if not, then why are so many people using it?

Stumbling across this in my search to make an accurate physics simulation. I am quite sorry I haven’t read all the comments, too many to follow being 3yrs worth. There will be more than 1 physical object in my simulation so I am wondering how the method works with multiple objects.

For instance to I do each RK4 step for each object:

PhysicSystem::Update()

{

for EachObject Compute Derivative a = evaluate(state, t, 0.0f, Derivative());

for EachObject Compute Derivative b = evaluate(state, t+dt*0.5f, dt*0.5f, a);

for EachObject Compute Derivative c = evaluate(state, t+dt*0.5f, dt*0.5f, b);

for EachObject Compute Derivative d = evaluate(state, t+dt, dt, c);

for EachObject Compute Finalize

}

Where the ‘accelerate’ function within evaluate will add any impulse and constant forces on the object. (So a constant force like gravity would be added 4 times).

I think that the method presented on this site is a modified version of RK4. It is different from the one presented on Wikipedia. On wikipedia they integrate 1st order ODE’s. y’ = f(t,y). The method presented on this should be able to do y” = f(t,y,y’) which is different from the RK4 formulated on wikipedia. The method used here looks like two first order RK4 methods stacked on top of eachother. Anyway, you can numerically solve 2nd order ODE’s with euler (error = O(h^2) per step, strange idea) or taylor (error = O(h^3) per step) or by using other methods.

For example, by using a simple hack, one can make taylor (y(t) = y(0) + y'(0)t + y”(0)/2 t^2) (having 2nd order accuracy) third order. Just add a fourth term: y(t) = y(0) + y'(0)t + y”(0)/2 t^2 + y”'(0)/6 * t^3. The third derivative can be obtained by using the backward difference on the second derivative: y”'(0) = (y”(0) – y”(-h)) / h. The backward difference method is known to have an error of O(h). When using the formula mentioned earlier, you get:

(t is replaced by h)

y(h) = y(0) + y'(0)h + y”(0)/2 h^2 + (y”'(0) + O(h))/6 * h^3 = y(0) + y'(0)h + y”(0)/2 h^2 + y”'(0)/6 * h^3 + O(h^4)

For the integration of the first derivative use: y'(h) = y'(0) + y”(0)h + (y”'(0) + O(h))/2 h^2 = y'(0) + y”(0)h + y”'(0)/2 h^2 + O(h^3)

The program:

http://www.wouterbuddingh.nl/index.php?page=projecten/physics/physics.php

uses this method.

(Note that the above method does not require mid step evaluation.)

The method on this site should get an error of O(h^5). It is unclear to me weather this involves the accuracy of the 0th derivative of the first.

I should have another look at it in the future. I tested the method of this site with the equation y” = -y, and the results were near perfect.

Probably, a better test would be y” = -y’ – y. The case y” = -y leads to a sine. “sin(t)” has an even distribution of derivatives at t = 0. This might cause some errors to remain unnoticed. Also, methods exist to verify the accuracy of a numerical method. One of them is richardson extrapolation. I have a book in which it is explained and i do not consider it to be easy.

I took a course on approximating partial differential equations. Here richardson extrapolation was used to analyse the shape of the error.

Keep up the good work! ODE’s rock! I recently read some nice things on networking too. (I need to port my linux socket code to windows some day.)

Regards,

Wouter

hi i’m struggling with implemeting RK4 to my code. Arent you solving 2nd order ODE with rk4 ?

hence formula taken from your examle looks like that x” = f(t, x, x’)

float acceleration(const State &state, float t)

{

const float k = 10;

const float b = 1;

return -k * state.x – b*state.v;

}

can we sove it like that ? i heard that it’s not allowed to solve 2nd order ODE with RK but i may be wrong 😛

Anyways in my code i’m solving symilar sytuation its oscillator with damper on a conveyor belt.

My eqation of motion looks like this x” = FrictionForce + Acos(wt) – k*x – b*x’ assuming mass of 1 kg

friction force given by conveyor belt drives this system into motion.

Will it be correct to solve it your way ?

float acceleration(const State &state, float t)

{

const float k = 10;

const float b = 1;

const float w = 0.1;

const float A = 1;

return FrictionForce + A*cos(w*t) -k * state.x – b*state.v;

}

i am attempting to solve a 2nd order ODE with RK4 – i’m not a mathematician, but this seems to work, and seems roughly equivalent to how one does symplectic euler integration for 2nd order DEs

I tried posting before but seems it wasn’t considered.

It was a bit messy anyway.

I really want to solve this problem, so I’d like any honest intelligent feedback.

My solution using an Euler type step works and gets the exact values.

Using the Euler method the standard way we would:

1. Update old velocity from the acceleration, v = v + a * dt, where dt is change in time.

2. Add new velocity to the old position, s = s + v, where s is position.

This over estimates the true value because it doesn’t consider the average velocity for the time step taken.

We have to consider the initial velocity and final velocity to calculate precise distance traveled in that period of time.

Also, it only works stepping by 1. Otherwise for example:

• dt = 5, a = 10, thus v = 10*5 = 50,

• if we consider time, s = 0 + 50*5 = 250

• not considering time, s = 0 + 50 = 50 … both wrong.

Note: If I have missed something?? Please say.

Using the Euler method my way:

1.Update velocity from the acceleration, v = v + a * dt, where dt is change in time.

2.Update position from 0.5 * velocity, or from 0.5 * dt (if you like! ;-)), s = s + 0.5 * v * dt, where s is position.

Example:

a = 10.

v = a * dt, where dt is the time-step.

Position(x) = position(x-1) + velocity(x) (left column) Euler original

pos(x) = dt * 0.5 * velocity(x) (right column) Modified Euler

dt=0: position = 0 velocity = 0 pos = dt * (0.5 * velocity)

dt=1: position = 10 velocity = 10 pos = 1 * 5 = 5

dt=2: position = 30 velocity = 20 pos = 2 * 10 = 20

dt=3: position = 60 velocity = 30 pos = 3 * 15 = 45

dt=4: position = 100 velocity = 40 pos = 4 * 20 = 80

dt=5: position = 150 velocity = 50 pos = 5 * 25 = 125

dt=6: position = 210 velocity = 60 pos = 6 * 30 = 180

dt=7: position = 280 velocity = 70 pos = 7 * 35 = 245

dt=8: position = 360 velocity = 80 pos = 8 * 40 = 320

dt=9: position = 450 velocity = 90 pos = 9 * 45 = 405

dt=10: position = 550 velocity = 100 pos = 10 * 50 = 500

You should notice that the exact value (right column) is an average of the last 2 Euler estimates (left column).

e.g. dt = 6, pos = 180 = avg(210,150).

With the method I made, it solves the problem of changing velocity.

A position can be calculated with any time step because the average velocity is taken into consideration.

The inaccuracy of changing acceleration is not present because the position is updated from the velocity.

* Velocity, is first updated from, acceleration * dt. (v = v + a*dt).

This is the velocity for the current time step.

* Position change is then calculated using the average velocity. (s = s + 0.5*v*dt).

I don’t understand how acceleration can affect it.

Acceleration changes the Velocity by a constant amount for each unit of time.

look mate it’s real simple: make the acceleration a function of position, velocity or time – such that the amount of acceleration given over a period of time dt is not easily predictable, and is certainly not constant — now you see when RK4 is useful

This method is called “Euler-trapezoidal”, which is known to be a simple enhancement of straight Euler.

Derivative b = evaluate(state, t+0.5f, dt*0.5f, a);

Shouldn’t this be

Derivative b = evaluate(state, t+dt*0.5f, dt*0.5f, a);

?

And similarly for the next line?

yes, you are correct – sorry that’s an error in my edit to the article vs. the actual source code i have on my drive, note also the comment above, the code is out of date vs. the article

thanks

Fantastic article Glenn. It must be the best one I’ve read on the subject of RK4. I just noticed a small difference between the code on this page and the source code. In the integrate function on this page you add 0.5f to the value of time for the call to evaluate:

Derivative b = evaluate(state, t+0.5f, dt*0.5f, a);

whereas in your source code you don’t bother:

Derivative b = evaluate(state, t, dt*0.5f, a);

I’m assuming that the source code is the correct implementation.

Thanks.

Hey I’ve used this code and article to help me understand RK4, how it works, and why it’s superior to Euler.

Anyways, this question was asked before, how does RK4 perform in terms of stability? although I’m not that familiar with integration stability I have read it refers to error accumulation over the life of the sim. In terms of stability how does RK4 compare to semi-implicit euler? When does integrator stability become a factor one needs to care about.

I’ve been looking for answers but not many people talk about integrator stability. Thanks, and great article!

i’m not really qualified to say, i’m not a mathematician

plus the stability of an integrator depends very much on the function being integrated… RK4 may be a very a good general purpose integrator with high accuracy, there are cases where other integrators are better suited.

one thing to consider is the difference between semi-implicit euler which tend to have better stability than explicit euler, at the cost of some accuracy. i’m sure that other higher order integrators exist with similar properties tending towards stability over accuracy – eg. higher order symplectic integrators, but i haven’t used them.

http://en.wikipedia.org/wiki/Semi-implicit_Euler_method

this page may also be of some interest to you:

http://beltoforion.de/pendulum_revisited/pendulum_revisited_en.html

i was also digging about and found this page contains a useful discussion about integrator choice, hope it helps:

http://www.cs.toronto.edu/~wayne/research/thesis/depth/node7.html

forgive my ignorance, but isnt Euler supposed to give perfectly accurate results with constant acceleration?

s = ut + 1/2at^2 gives exact results, but simple euler integration like this: v += a * dt, s += v * dt does not

the reason being that the constant acceleration means that velocity is changing constantly, so the inaccuracy is incurred in the s += v * dt term

cheers

got it, thanks

thank you!

Thanks for writing this great article. I tested the code yesterday, and was thinking of using it for improving integration of measured data when I suddenly realized that this was not possible:

I was thinking of measuring acceleration and using RK4 to get better results than just calculating x=v0t+1/2at^2. In order to use RK4 or any other numerical integration you need to know how the acceleration is going to evolve between the starting value and the end value for a continuous function (f(t,x)). This is not the case when doing measurements of discrete data. Am I correct?

to be honest, i’m not sure – perhaps it would be better, but you would end up comparing an RK4 on a data sample where you use 4 samples, vs. a s = ut+1/2at^2 per-sample, i’m not sure which would win, perhaps you can try it on your data and see how it goes?

When I start to think about it you actually only need three points (a(0),a(1) and a(2)) with the same time period in between, not a continuous function. The algorithm seems to be working fine, but the displacement s seems to be drifting a bit compared to the ordinary function, that is even if I am high pass filtering it. I will look a bit more into it to compare the two methods.

Hi Glenn….

I am unable to compile your C++ code since I don’t code in C++ but I was interested to know what the resulting effect of this code should be. I have ported it to flash as best I can and wanted to make sure that I have the correct results. The data traced out are as follows… ( sorry for the large amount of splurge)

[ object ParticleState x=100, v=0]

[ object ParticleState x=95.20416666666667, v=-93.57916666666667]

[ object ParticleState x=81.88127307291667, v=-169.42547116319446]

[ object ParticleState x=62.09968928156894, v=-222.06922653612696]

[ object ParticleState x=38.34043852154505, v=-248.75047511755145]

[ object ParticleState x=13.223832821344754, v=-249.42161761040853]

[ object ParticleState x=-10.751027285679594, v=-226.49385796348292]

[ object ParticleState x=-31.43053241877831, v=-184.37576179993152]

[ object ParticleState x=-47.176906610963826, v=-128.86724709778446]

[ object ParticleState x=-56.973670392117775, v=-66.47994304603026]

[ object ParticleState x=-60.4624457865241, v=-3.755052134795548]

[ object ParticleState x=-57.91416230690944, v=53.35668147004745]

[ object ParticleState x=-50.14362181820809, v=100.00030063526725]

[ object ParticleState x=-38.380872488645416, v=132.7704915180671]

[ object ParticleState x=-24.11563785802634, v=149.89498870909898]

[ object ParticleState x=-8.932043928178269, v=151.2464396885771]

Do these correlate ( or at least pretty close ) with the expected values from the C++ version of the code. Many thanks and thanks for the tutorials..they’re great.

RipX

it looks about right but i can’t be sure, if i remember correctly it is meant to oscillate from right to left as a damped spring — try visualizing the output as a point on the screen and see what you get?

cheers

Thanks for your great articles! This one actually inspired me to write my own article, so I made some investigations regarding the relationship between global error and overall calculation time for several different integration schemes solving the differential equation of the damped harmonic oscilator and the magnetic pendulum. The results are on my home page:

http://beltoforion.de/pendulum_revisited/pendulum_revisited_en.html

Best Regards

nice work ingo, cheers

The picture shows a significant difference. However, i also put the RK4 implementation of this site in the test. (Copied the 1D source and changed the differential equation to a = -d.

RK4 handles the ODE flawlessly (compared to the modified RK2 method).

I just read a post about integrating objects simultaniously.

This is something i was worried about,, since this was not part of the examples.

If you don’t do this, RK4 loses it’s accuracy because the derivatives sampled will not have first order accuracy any more. I think this violation will cause RK4 to not be fourth order.

In my previous posts, I was very critical about this RK4 implementation.

It is hard to see that the implementation is indeed 4th order.

(The implementation is not identical to the RK4 spec on wikipedia.)

A lot of aspects seem questionable.

You should emphasize the accuracy of the 1D algorithm by showing a few pictures of RK4(state based) in comparison with other techniques.

(ODE integration plots, in combination with the precise solution).

The problem is that not many people are able to verify the workings of the integrator in an analytical way.

The article i wrote about the modified RK2 tries to do just that.

The problem is that such an article can contain errors.

A proof that RK4(state based) is fourth order accurate will be much more difficult. Testing the integrator is a reasonably good approach.

It seems like RK4 is a good total solution.

Ehw,, don’t know if the SDL program will compile on a Mac (it should). (It does compile on Linux.)

Anyway, running the program is not really relevant, i made a picture that illustrates the difference between a method that is totally second order (blue), and the tweaked third/second order method (cyan).

Again, the ODE is a = -d. The reference sine is omitted (to make it clearer).

http://www.wouterbuddingh.nl/projecten/physics/plot.png

This is what i got when entering the formula into a test program:

http://www.wouterbuddingh.nl/projecten/physics/RK2.zip

Note that the test is subjective.

Use the left mouse button to drag objects around.

The right mouse button draws a convex polygon on the background.

“Enter” makes this polygon an object (the moments of inertia should be correct).

Greetings

how about some graphs showing relative error over time for euler, RK2, modified Rk2, RK4? — really this stuff can all be measured, pick an ODE and lets see how it works…

A moment ago i uploaded an application that demonstrates the integrator.

The differential equation used is a = -d, (a sine).

I think this is a good choice.

The reference sine is red.

The green sine is integrated using Euler.

“The blue sine is integrated by using RK2″.

I don’t really know how to just use RK2 to integrate a physics state, that’s why i modified it. The position is integrated by using a second order Taylor polynomial.

The velocity is indeed RK2 (the midpoint method was used).

Both position and velocity have an error of O(h^3).

The cyan sine, is integrated by the method illustrated in the article i posted yesterday (based on RK2). Here the position has a third order accuracy and the velocity has a second order accuracy. Note that the integration of velocity for both the blue and the cyan sine is identical. The integration of position makes the only difference.

The program can be found here:

http://www.wouterbuddingh.nl/projecten/physics/integrators.zip

Use W,A,S,D (gaming controls), to modify the number of periods and step density.

Today i have written an article about using RK2 to integrate game physics.

The zero’th derivative (position) can be integrated with third order accuracy, and the first derivative can be integrated with second order accuracy.

The article is primarily based upon the article Gaffer handed to me today.

Basically RK2 just makes up an extra derivative.

A normal taylor series expansion would give second order accuracy for position, and first order accuracy for velocity (Euler).

This concept could be used in RK4 as well (don’t know how).

This way, position integration could be done with fifth order accuracy, velocity with 4rh order. RK4 would simply add tree derivatives.

The article can be found at:

http://www.wouterbuddingh.nl/projecten/physics/RK2.pdf

(It contains a disclamer :)).

nice! now lets see how this integrator performs in real world cases!

Yesterday, I have rewritten the article. Now it is just about RK2 without extra positional derivative extension.

Today, i made a draft of an RK2 integrator. (Based on the article that was supplied.)

The error of each step is of O(h^4). (3’rd degree accuracy)

Another thing that it seems to show is that it doesn’t affect the order of error if you use Euler steps instead of Taylor steps.

Probably, this also applies for RK4.

It seems like it’s better to just use Euler.

When this integrator is finshed (and checked for error’s), i think i’ll publish it.

Greetings

The article is great. It does a better explanation than some YouTube videos.

A simple error analysis of the midpoint method does not prove that using a Taylor gives a better result that is one order better than Euler. (When you are looking at the formula alone.) The order of error of the final step is not in Taylor form and therefore unknown. However if this error were O(h^2), the midpoint method would not be O(h^3) any more. This is where i get stuck.

Perhaps the article can lead to a solution.

I’d like to comment on the use of the word vector in my previous post. Here, i do not mean (x,y,z) but a vector of derivatives like (position, velocity).

yes, my suspicion is that it either does not work – or it does not change the actual order of the error, remember the order in RK4 comes from the fact that 4 derivative samples are taken and recombined to form a taylor’s series expansion with error O(5) — not from the actual prediction to determine the velocity between each step

Dear Gaffer,

Today, and yesterday i have been occupied a lot with numerical integration. In my previous post i mentioned that i haven’t looked at the RK4 implementation of this site for a long time. Today i checked the sourcecode of NetworkedPhysics.

The process is actually reasonably clear.

I especially looked at the evaluate function. This function does an Euler step according to a given derivative and then evaluates the forces at that given time/position/velocity.

This is a copy of that function:

static Derivative evaluate(const Input &input, const std::vector &planes, State state, float dt, const Derivative &derivative)

{

state.position += derivative.velocity * dt;

state.momentum += derivative.force * dt;

state.orientation += derivative.spin * dt;

state.angularMomentum += derivative.torque * dt;

state.recalculate();

Derivative output;

output.velocity = state.velocity;

output.spin = state.spin;

forces(input, planes, state, output.force, output.torque);

return output;

}

It does what it is supposed do do. An Euler step is done, later the derivative is queried. The velocity is copied from the original state. If you modify this, i think you could model Impulse as well. I think impulse is considered to be a secondary physical quantity, so it is quite understandable if this is omitted.

The thing that i want to address is that the Euler step gives a prediction. This prediction has an error of O(h^2). In the (wikipedia) secification of RK4, also Euler steps are used. This is because in the specification, the problem does not contain any unnecessary complexities.

I claim that by making the error of these steps, smaller, the resulting error of the enitire RK4 step will also be smaller.

The evaluation of a state is done with an Euler step. This can be improved because you do have more information about the curvature.

Instead of writing “state.position += derivative.velocity * dt;”, you can write

“state.position += derivative.velocity * dt + 0.5f * (derivative.force / mass) * dt^2″.

This increases the accuracy of position by one order.

A similar approach applies for orientation. (By using a weighted sum of two radial vectors. This sum needs to be converted to a quaternion.)

The momentum and angular momentum are just fine because there are no more derivatives left. (i.e. you cannot do better than Euler)

Sometimes Euler steps are taken to just get a quick computation. However, i always find it hard to make an integrator accurate. (Made a few.) I therefore think that all Euler steps should be replaced with Taylor steps if possible.

Greetings

Hey Wouter, although that logically makes sense from a programmer point of view, i’m not sure that the mathematics are correct – you should perhaps try to work through the derivation of RK4 from the original taylor’s series expansion, then you can explore the mathematics side and see if it works there

here is a good article that may help

cheers

another good article to look at:

http://pathfinder.scar.utoronto.ca/~dyer/csca57/book_P/node50.html

Glenn,

I’ll try to read the article. The thing you are suggesting (mathematically proving this), i tried to do this for RK2. The problem was that i got stuck somewhere.

I’m still planning to do this some day, maybe your article will help.

The principle of RK2 is the same, but less complex. In RK2, Euler steps are done. The objects that are integrated are numbers. Their derivatives are analogue to velocities.

The best way to integrate these objects without using any unknown derivatives, would be to use Euler.

In the case of a physics state, the best way to integrate without using any unkown derivatives would be to use a Taylor series expansion.

RK2/4 are based on real numbers, applying it on a vector might change a few things.

But you are right, in order to know this makes a difference, it should be proven by using the derivation of RK4 or 2 (that’s my goal). It might be that this does make a difference, but not a significant one (i.e. an order).

Greetings

Recently i have also been occupied with the integration problem.

I have not read all the articles on this page so forgive me if i’m posting something that has already been discussed a number of times.

In a previous post of mine, i promoted the usage of taylor polynomials instead of RK4. I think implementing RK4 is too complicated. I once copied part of the Networked Physics source. And wondered, what if i made a mistake.

I also think that implementing RK4 in a more complex setting leads to problems.

Before i explain my views, let me give an overview of some other methods, that can also do the integration thing.

Euler (O(h^2) error)

Taylor d(h) = d_0 + v_0 * h + 0.5 * a_0 * h^2 (O(h^3) error)

Taylor + RK2 (Probably O(h^4) error. Still trying to prove this.)

RK4 implementation of this site (O(h^5))

(I don’t think i really understand this, haven’t looked at it recently.)

Now, about the approximation methods. What puzzles me about this site, is that the inaccuracy of Euler is illustrated (100% right). The solution that is proposed is to use RK4 instead. The problem was that the integration was inaccurate, the solution is to use a complex numerical method that was designed to approximate solutions for ordinary differential equations to do the integration.

This method is appropriate if you want to model ODE’s in your system. If you just want correct integration, using your high school physics will give just the right results (a taylor series expansion). (I already mentioned this in my previous post.) Some people (like me) come to this site with the problem “my Euler integration sucks”. They are handed a bazooka to kill a mouse.

Coming back to my previous statement that RK4 leads to problems.

Suppose you have 5 objects. They are all connected with springs. (They form a complete graph i.e. every object is connected with every other object.)

The forces are evaluated at t = 0. The forces applied on the object you are trying to integrate depend on the positions of all the objects. (This is not yet a problem.)

By using the forces obtained, you do an Euler integration to the middle of the frame (t = 0.5). Now, a slight problem arises. To evaluate the forces at this time, the other objects also need to be “Eulered”. They must be in the same time.

The problem that i’m trying to illustrate is, that objects aren’t integrated independently. To handle this consistently you should integrate the entire system instead of each object separate. This problem can be solved but you need to think a lot. Using RK2 instead of RK4 makes this a lot easyer. By using RK2, the time between each evaluation always increases by 0.5. Usually it is preferable to just have an update function that is called successively than to have an update function that is parameter dependent. Defining game actions incrementally is preferable in my view.

Currently i am toying with an RK2 + taylor integrator. A double state is kept. Every force evaluation, the entire system is taken into account. Also all evaluations are called successively. This makes the “game programming interface” a lot easyer.

Greetings

I’m just trying to figure out whether it’s worth adding RK4 to my current project, replacing the current Euler method. I’ve got characters walking around whose legs are represented by springs. At the moment, I’m seeing quite a bit of jitter when the characters come to rest, because of the spring forces fighting gravity across the timestep (the character shoots up, gravity takes control, then it bounces up again).

Would an integral method such as RK4 allow me to let these springs come to rest without jittering?

first step is to consider using springs that are less stiff, this will reduce the jitter

at this point if you want even more accuracy, then you can try switching to RK4

Yes, I’ve played with the spring constant a little. Lowering it too much, while allowing th springs to come to a smoother equilibrium, means that the equilibrium is at about mid-shin-height, which means characters’ legs are bent when they’re standing still. Also, it gives a less desirable efect when moving and jumping. I’ve found a value of the spring constant that gives me good equilibrium extension and movement, but jitters.

I’m going to give RK4 a try – thanks for your thorough explanation!

excellent, now you should be able to get your springs significantly tighter using RK4 – you may find this enough, but if you still find it too saggy (after all springs are reactive, they only work to correct error that has already happened) then there are other techniques you can explore

there is an old article about this about “velocity verlet” based integration of particles for ragdolls you should look at, it’s by thomas jakobsen

http://www.teknikus.dk/tj/gdc2001.htm

and lots more resources on this page:

http://en.wikipedia.org/wiki/Ragdoll_physics

the key point is that when you want “stiff” ragdolls you probably need to use different methods than spring based physics

cheers

Found this looking for a replacement for my Euler integration.

Thanks so much for making this available–quite a good read

thanks, glad you enjoyed it!

Thank you!

just doing some very simple scripting for 3dsMAX and found your stuff extremly useful.

/Mattias

one of the ones using Euler…

OK i forgive you

You are a god send.

I have been working on my first advanced game over the past month and have put physics off to the very end because I was not sure how to handle time stepping.

Thanks to your articles I feel that things will go a bit smoother than planned.

I will definitely be frequenting this site in the future, all the articles here are great.

Awesome article, but just wanted to give you uber-kudos for the link to the study about incompetence. It rings so true of the ‘higher-ups’ at our studio it made me almost fall off my seat. I just hope I’m not equally incompetent and unaware 😐

Love the whole series, tho!

So, I was looking at the RK4 code and such, and though it was quite readable, I had to step through each step to understand what it was doing. In the end I realized that it seemed to just be a taylor polynomial for the position of the object, which I had originally thought about using before coming across this article. And after solving with variables it turned out that they were the exact some equations. Correct me if I’m wrong.

Basically, I realized that

n=infinity

x(dt) = Sigma ( (nth derivitive of x() at 0) / n! * dt^n )

n=0

I just learned about taylor polynomials in calculus, so I thought about it and it makes sense. What that taylor polynomial is saying is basically that x = initial x + initial velocity * dt + initial accelleration / 2! * dt^2 + initial jerk / 3! * dt^3 + …

So basically, you can take more terms to get more precise, and RK4 just goes up to the jerk. To make the changes for velocity you just take the derivative of the polynomial you used for position, and so on for acceleration. Though this may be slightly different as I am coding the integration for a physics simulation where the acceleration is independent of the object’s current state, (i.e. it depends on a list of forces acting on the object in stead of an equation that depends on the object’s current state). So the RK4 code can be reexpressed as such supposing you store acceleration in the state as ‘a’ and calculate the acceleration from a list of forces:

float newA = Acceleration(forces)

float jerk = newA – state.a

state.x = state.x + state.v * dt + state.a * 1.0/2.0 * dt^2 + jerk * 1.0/6.0 * dt^3

state.v = state.v + state.a * dt + jerk * 1.0/2.0 * dt^2

state.a = state.a + jerk * dt

I felt that these equations displayed this way were mathematically more expressive of what they actually do. Though if the acceleration depends on the object’s state, I think these equations should still work.

About reducing the Euler time step: you touched on this briefly but, in general is it possible to reduce the time step so much that Euler Integration would be as good or better than RK4 Integration was at the original time step?

(Or put a different way, ignoring any questions about processing speed, can you tell us more about what RK4 Integration is doing that is superior in nature to Euler Integration?)

Great article, thanks.

For a given timestep, you get more accurate integration using RK4, because RK4 has error order 5.

Yes you can decrease the timestep with euler to get as accurate or more accurate than RK4, but you would have to do such ridiculously small step sizes to reach this point that it becomes inefficient

It is better to use RK4, because you get better accuracy per-step, and don’t need to drop to really small step sizes when you have stiffer springs

cheers

here is a very good article and demo showing the accuracy differences between euler and RK4, take a look!

http://www.gamedev.net/community/forums/topic.asp?topic_id=380499

Thank you for this great article. I’ve learned a lot. However, I’m not sure how I would implement a generic rigid body

acceleration()(orforce()) method. I’m used to the explicitaddForce, which is quite simple to understand and implement. Without it, how can I “add force”, say, once a collision occur?so you may model spring like forces for collisions, making the collision forces a function of the particle position and velocity over time — take a look at the spring physics example, or the source code from the networked physics example, it does this

now i’m not saying this is the best way, but it’s a way, and a good way to get started – another name for it is “penalty method”. of course it has the downside that all collisions are reactive on springs so it’s hard to tune and objects have a bit of cushiness in their collisions, and penetrate into each other when resting

other methods include impulse methods, i think for these there is really no need to use RK4 unless you have some other hard to integrate DE, the reason being that i use RK4 in this article series to make the spring forces somewhat stable for collision response

cheers

I found your article very informative. However, I was wondering why you were so quick to dismiss Euler integration.

You can’t get something from nothing. Improving an algorithm in one way means worsening it in another.

From what it looks like, RK4 works best where the fifth derivative is negligible. While this is a reasonable assumption in many types of simulations, it by no means proves that RK4 is superior under all circumstances. In cases with high order derivatives, Euler integration actually gives less error per sample point.

Also, sorry if I accidentally posted this multiple times. I’ve been having a lot of trouble trying to post a comment.

1. the article says nothing about RK4’s load on computer resources. If I have, say a couple dozen objects falling down – will I have jerks?

2. For realistic computer modelling for, say, training purposes this is great. But what sort of Flash games (I was refered to here from Flash forum) requires such precision? Does it really matter if a calculated position of an object will be 500 pixels instead of 450?

RK4 is slightly more expensive, given that it requires 4 steps for each timestep, but conversely, it may allow you to take larger physics timesteps so it could end up being more efficient overall.

whether or not it will be too slow depends entirely on your own physics simulation. the cost in the example you give would be more likely due to the collision detection, the integration is usually a very tiny part of the CPU cost – if you have thousands of particles however, you should probably not use RK4

whether or not accuracy counts is also entirely simulation dependent. if you are simulating springs and need accuracy, then Rk4 provides more of it, – if you are just doing simple ballistic motion, then you can get away with something simpler, perhaps s = ut+1/2at^2

Here is a really good lecture series covering many different numerical integration methods, and their mathematical foundation:

http://www.damtp.cam.ac.uk/lab/people/sd/lectures/nummeth98/odes.htm

if anybody is interested, here is a paper describing some alternative approaches to RK4

http://citeseer.ist.psu.edu/cache/papers/cs/28838/http:zSzzSzwww.uni-tuebingen.dezSzunizSzopxzSzreportszSzhairer_182.pdf/hairer03geometric.pdf

To get the benefits of RK4 it is essential that you recalculate all forces at each step of RK4, otherwise your drag force is constant over the frame time – and you are really just doing an euler integration

cheers

I want to implement drag in my game simulation, which I think can be fairly approximated as Fd = v^2 * some_scalar.

Should I use the velocity derivative samples, or just the initial velocity? From a code design perspective, I was hoping I could calculate all the force vectors working on objects before the actual integration step.

Though I’m asking for what would be the most ‘correct’ approach, I know it doesn’t really impact the game simulation. In the event that the drag force becomes bigger than the opposing force, acceleration changes into retardation, until the decreased velocity makes the forces cancel each other out, i.e the acceleration becomes zero. Right?

the solution you presented is not exact, because your assumption of constant acceleration over the frame is incorrect.

why? well for any interesting simulation acceleration is a continuous function of other variables such as current velocity (for example, drag), and position (springs). since these quantities change over the frame, the force which is a function of them

alsochanges over the frame.in these cases, RK4 provides better accuracy – because is is able to detect and correct for these changes in acceleration over the timestep

if you don’t believe me, then code up a spring simulation using your method, then run with a few different timesteps, if as you say the timestep is irrelevant because your solution is exact, then of course you’ll get exactly the same result from your simulation no matter how large a timestep you choose. but when you don’t get the same result with different timesteps, come back here and read the article again, and see if you can now begin to understand what RK4 is actually doing

cheers

Why use RK4 for this problem when it can be solved by a simple integral.

I tried the approach of RK4 described on your site. It worked.

Within my (numerical) simulation, accelleration (or force) is constant during the current frame.

Forces get accumulated (into a variable) at the beginning of the frame.

So, accelleration is constant, which means the velocity is linear, which means that the distance is quadratic.

Using basic integration, you get something like this:

d_0 + v_0 * t + 0.5 * a_0 * t^2 (represents the distance at (t))

Where the _0 denotes that those are the values at the beginning of the frame.

At the end of the timestep, the velocity has changed, (v_1 = v_0 + a_0 * t). v_1 will be v_0 for the next timestep.

The term d_0 + v_0 * t + 0.5 *a_0 * t^2, is a Taylor (actually a Maclaurin series).

d_0 * t^0 / 0! + v * t^1 / 1! + a * t^2 / 2! (….)

You can simply make an array of derivatives (starting with der. 0 to define your motion as precise as you want).

This solution is exact, that means that the size of your timestep is not relevant. Your timestep should be small enough to get a descent user/envoirement interaction.

I thought of it, but it sounds memory costly as, since the program needs to keep all steps in memory, instead of just the result after one RK4 step.

Instead I’ll see what happens, if I save the position of the first mass and use that (and not its new position) for the second mass.

yeah that should work, provided only the distance between particles feeds in to the force calculation (this should be the case for simple F=-fx springs, but anything more complex, you’ll eventually end up just duplicating the state 4X)

cheers

why don’t you try evolving the whole system RK4 simultaneously?

eg. take four steps of the whole system to perform one RK4 step — only sample forces from the distance between points from the same step

cheers

I have now implemented the program for multiple masses and springs. It works, but there is problem (I think). Testing it with 2 masses connected with a spring shows that the masses don’t move the same relative to the CG point. This is because one mass is first integrated, which changes the length of the spring. Because of this, the force is slightly different when the other mass is integrated . This problem becomes bigger when the ratio k/b gets larger, in which case the system is displaced to one side, i.e. CG has moved to one side.

oh duh, I didn’t realize the article was about non-uniform acceleration until after I read the comments, because your euler example uses uniform acceleration with a constant of 10. I guess I should have paid closer attention to the title.

thanks for the articles, very informative.

well if the car physics has constant acceleration it would be OK, but if the acceleration was non-constant (eg. changing gears, or changing the weight on the accelerator), then RK4 will give you less error than euler

that’s really the point of this article, as per the best integrator for a simple constant acceleration situation, try this one: s = ut + 1/2at^2 it’s an exact solution

cheers

What’s wrong with setting dt = currentTime – startingTime, in the Euler example, where currentTime = time(), and startingTime is whatever time() is when the car starts moving? wouldn’t that be completely accurate and useable code?

Thanks for a nice tutorial. It was easy to implement for a single mass and spring, and now I want to expand it to multiple particles and springs.

In many other examples I’ve read (e.g. http://www.pixar.com/companyinfo/research/pbm2001/pdf/notesc.pdf), the state phase vector of a particle has a force entry, and the order is to first calculate the force and then to calculate the new positions, but in your code above the calculation of the force (acceleration) is part of the integration that calculates the positions. Any comment on this would be much appreciated.

Hi everyone,

I am implementing a spring mass system and I need to solve my differential equations. I tried the code posted here but my masses oscillate progressively. As I understood, I need to integrate all my objects (masses) at once. However I did not understand how it will be implemented. I tried various schemes but it did not work. I guess I need to recalculate all spring forces acting on a mass 4 times in RK4. Any help or tip will be appreciated…

Thanks

Zeki

yes, you are right about the re-evaluating forces at multiple points — this makes it hard to structure your simulation (you need an “evaluate” function that calculates forces on demand, instead of just an “applyForce” function)

for this reason, i would not recommend RK4 unless your simulation suits this model well – RK4 is easy once you understand it, but structuring your simulation in this way is non-traditional, and difficult unless you are doing everything with spring forces

cheers

Oh and yeah another thing, since the force function accumulates all forces acting on all objects for the particular time step, this would also include collision response forces right? It sounds like doing a full collision detection and response pass for all objects four times per frame might run the risk of getting very expensive, unless you can set the update rate four times lower when using RK4. Is this not true or is it at least very likely that I can set the update rate this much lower without numerical inaccuracies growing too high?

Thanks

Two things stood out to me in this article compared to other RK4 articles/explanations, and both are positive.

1. You skip the rigid body factors to focus on the integrator

2. Your code is actually readable as opposed to another implementation I’ve seen where the author (who shall not be named) seem to have accidentally pressed caps lock and written the entire code piece without looking up on the screen.

The thing that has always bothered me with RK4 (and other multi step methods) is this thing about needing

to be able to evaluate a force function for every object at any point in time. This means that I must be able to take into account everything from input all the way through AI, game logic, etc. Unless I’ve misunderstood things, this would require quite the special kind of game loop so if the game (and engine) hasn’t been designed for this, then you’re pretty much screwed. And I’m not sure I would even know how to design this kind of system robustly when introducing threads.

I would love it if you could tell me that I’ve got it all wrong and that this is easy as pie but until then I’m sticking with verlet integration.

Thanks for a pleasant read.

Regards

its to do with inertia tensors, if you have a proper inertia tensor setup, the angular velocity can actually change while the angular momentum is constant

the classic example of this is a spinning top

search for “conservation of angular momentum” and you may find more details on this

but yes, since you are simulating billiards they have a uniform tensor (or perhaps you just ignore its effect, same thing), so you could actually use angular velocity in this case with no problems

cheers

Thanks for a great article. I implemented your integration in an xna billiards game I made just to test drive RK4. It was both fun and I learned much.

You say:

Later on when working with rotational dynamics you will need to work with angular momentum directly instead of angular velocity, so you might as well start now

Could you elaborate on that please. I don’t understand the “need” to work in momentum when dealing with rotations. I think of angular velocity as simply a vector similar to lineal velocity but whose accelerations are tempered by an inertial tensor instead of mass.

andre, i recommend the right integrator for the right job. for cloth especially an RK4 integrator may be too expensive, a simpler integrator that gets good enough results is perfectly fine

the idea of me presenting RK4 here is just to demystify it – people tend to think its super hard to do, but its not!

right tool for the job, and all that

cheers

dm: that would actually work perfectly, and without error – given contant acceleration, but what about the case where acceleration is non-linear?

in this case RK4 works much better, because it can take into account change down to the fourth derivative, any error in RK4 is O(5) from the taylor’s series expansion

so in the general case RK4 is better than the integrator that you present, but if you have a linear acceleration only, perhaps a ballistic point under

constantgravity, then the integrator you presented works just finecheers

Thanks for the article. And that APA link too!

Glenn

What are your thoughts regarding the use of a Verlet integrator over RK4?

Described here:

http://www.teknikus.dk/tj/gdc2001.htm

I’m intrigued, if acceleration is applied during the timestep, then the velocity isn’t constant and so the position must be calculated with respect to the acceleration. So, why wouldn’t the following work?

while ( t <= 10 )

{

float accel ( force / mass );

position+= (velocity*dt) + ((accel*(dt*dt))/2);

velocity+= accel * dt;

t+= dt;

}

http://www.apa.org/journals/features/psp7761121.pdf

And I’m even more of a pratt, now as I’ve read the other comments and found youv’e already addressed this. Gah!

Apparently I’m a bloody idiot 😉 I agree with you about RK4 being vastly superior to Oily’s method, but am struggling to work out how to apply RK4 to quaternions for angular velocity..

Sorry for my stupid question above. I’ve just looked at the source code and found that you’ve overloaded the evaluate functions. I was looking at the wrong one!. D’oh.

Keep up the good work.

Hi,

Great article!. I’m having trouble understanding the evaluate function. Yes I see that you’re getting the current state vector by using the derivative values passed into the function, but what I don’t understand is why you’re setting the rate of change of velocity in the derivative to the current velocity.

Thanks.

Great job my friend.

trajectory isn’t useful either. I get links, but nothing that talks about it from a numerical integration standpoint. I’d like to write an article on it for wikipedia, but I’d need a definate reference on what to call it to proceed :/

@48 – trajectory may be a better search term – remembering that g is just a constant a

I’ve been trying to find out what the integration method for position: r(t) = r(0) + v(0) * t + a(0)/2 * t ^ 2 is called. You seem to call it a ballistic integrator. I’ve also heard it called a “parabolic” integrator by Erin Catto in the Bullet Physics forum. However googling for either terms doesn’t provide any useful references. The only reference I can find to it is on Wikipedia where it’s a member of the SUVAT equations (so I call it a SUVAT integrator for lack of a better term).

I’d greatly appreciate any links to references on what to call it. Makes me feel less stupid when I discuss it if I know what it’s called 😛

sure, thats just

s = ut + 1/2at^2

its a ballistic integrator, works perfectly for constant acceleration

but, game engines don’t often have constant acceleration all the time =)

…it did strike me when reading that if you were doing updates at each dt=1 you’d be more accurate to use… newpos = pos + vel + 0.5 * acceleration …which would also save on a bit of math!

Excellent indeed!

Thank you for publishing this (in fact thanks for the entire series)!

Brilliant explanation! Thank you.

well, passing is as reference avoids a copy – and the state structure can get quite large, so this is important

now, when we want to avoid the copy and dont want to change the value, we use a const reference to pass in

this is all standard c++ practice!

cheers

Hi,

Brilliant article (and others too). I have a question though. Why do you use references of variables as a parameters in a method calls if you do not intend to change their value?

Derivative evaluate(const State &initial // <-

float acceleration(const State &state // <-

void integrate(State &state // <-

In my opinion code gets really messy.. But what’s the bright side?

Cheers.

This article explains how:

http://www.gaffer.org/game-physics/physics-in-3d

Great page been looking 4 something like this. I was wondering if you could help. I need to RK4 quaternions, how would the code above need to be changed to achieve this? thanx

Euler, Improved Euler and Runge-Kutta (4th Order) are all taught as numerical methods at my university in the “Differential Equations” class.

I hated doing them by hand because I would always mess up so I wrote a C++ program for each method (yes, it was easier and it worked). I remember we had to compare the error values and such for the three methods and Euler’s method is only accurate if you choose a very small step size, whereas Runge-Kutta is much MUCH more accurate with a pretty large step size.

Great work on explaining, hope to see more on other math principles around the site.

May I just ask why the position has to be updated before the velocity?

Why not the other way?

If you do it the other way, what are the consequences?

Thank you.

sure, but what you described there is ballistic integration not euler integration – whats your point?

Well, the problem with updating our displacement in your original Euler example is that you didn’t do it right. Under acceleration change in displacement is the integral of change in velocity. An integrator such as what you have demonstrated approximates this integral, thus the name “integrator”. But there is no need to approxiamate, we can just include the integral in our calculations:

d’ = d + v*t + 0.5*a*t*t

v’ = v + a*t

we can quite easily see that this is “right” because the change in velocity is the exact derivative of the change in displacement with respect to change in time. If we take the derivative of change in velocity with respect to change in time, we get a’ = a, which again holds to the basic definitions of velocity being the first and acceleration being the second derivative of displacement.

I think the point is that the method you’ve chosen for integration of velocity to get position is inappropriate. Euler integration works well to take acceleration to velocity because acceleration is assumed to be constant over the interval. Velocity however is not constant; it varies linearly throughout the time step. Euler integration will not work well in this case. Instead we usually employ trapezoidal integration (which under the assumption of constant acceleration is equivalent to capn_midnight’s approach and, again under the assumption of constant acceleration is an exact approach). It is usually implemented as:

v’ = v + a*dt

p’ = p + dt*(v + v’)/2

Your article seems to be suggesting that one integration method can be used for all cases and while RK4 might work in a practical sense for the types of simulation you are doing, I don’t think this is generally true. You might eventually find that it doesn’t really work as well as you think it might be working. I can honestly say that in the 30 or so years I’ve been writing real-time simulations I have never seen anyone use RK4 integration (beyond an academic environment).

I like your presentation of the RK4 method and I agree it is not the “big deal” that most people seem to make of it. I also agree that it does not impose much additional overhead on the simulation; however it is not necessary in this context. I have seen missile simulations run accurately and correctly as low as 10 Hz using the appropriate methods described above. I do not like your presentation of Euler integration in an inappropriate context. It propagates a common mistake that I have seen and battled throughout the years and to this day causes me no end of grief.

Thanks for the great document Glenn!

i’ve been laughed at for saying yewler, live and learn

Very interesting and educational piece, just what I was looking for.

While you are right about the pronunciation of Euler being “Oiler” in his native tongue, I doubt you’ll ever hear people pronounce Einstein as “AynShtine” or Heisenberg as “HaysenBehrg”, so I would like to think that “Yewler” is a perfectly acceptable, anglicised way of saying his name (Just like we say “IneStine” or “Hisenburg” instead of their actual names).

That’s the only point in your write-up that I had any sort of issue with ,)

Thanks again!

I love your writing style! “If you use Euler then you are a bloody idiot” Love it 😉 Thanks for the great explanation!

Don’t disparage Verlet methods – RK4 is accurate enough for game simulations, but unfortunately does not conserve energy over long periods of time (which doesn’t really matter for game physics, but does if you want to do real physics). Google “symplectic integrators” for more details.

Awesome and very cool way of explaining things to programmers. I liked your approach. Thank you bro…

Thanks! Yes, you’re probably right. According to wikipedia the a,b,c,d are different slopes calculated during the interval, which in the case of discrete acceleration is not possible to get at.

Problem, from what I’ve read, is that using euler with an accelerometer is bound to give lots of errors very fast (ie http://instruct1.cit.cornell.edu/courses/ee476/FinalProjects/s2005/mouse%20webpage%20KM249_AK288/conclusion.htm). Dunno what other methods to use, I may have read something about using Simpsons Rule, whatever that is.

actually, i take this back – no i dont believe you’ll get any benefit with RK4, with the way i described it above, because fundamentally you need to be able to sample the derivative quantity in ways which your device does not support

note that in the calculation of derivatives a,b,c,d that you pass the previous derivative each time you evaluate the next, eg. evaluating derivative b, requires derivative a, c requires b etc. — this is meaningless when you just have a discrete function of a(t), there is no way for you to do this calculation

so i’d say, you probably need a different technique, which will be more like a filter, not RK4

failing that, you could just try euler integration of acceleration to find velocity, then euler integrate velocity to find position – and see what sort of results you get

cheers

just hooking up the current value into the acceleration() wont work

you’d be best to take the raw acceleration values to a buffer – then process them 4 at a time, using the middle ones as the extra samples

failing this, try doing it with a simple explicit euler integration first, just to get your head around it

cheers

sure, but you’ll get drift over time – try it and see

I’m currently dabbling with using an embedded 3-axis accelerometer. I’m trying to find a way to get the velocity and position from it by twice-integrating the acceleration values from the device.

Do you think it would be possible to use RK4 for this, ie just rewrite your acceleration() function so the it simply returns the current reading from the device?

(bridges): Midpoint or RK2 may be what you need is you are just using constant acceleration, Euler as from the article, fails at even constant acceleration.

(Glenn): Awesome article, thanks to your explanation of RK4, I was able to create a separate class structure for integrating ODEs(numerically) separate from the physics(although that is its primary purpose). The ODE solvers rely on operators and functions from a generic derivative class, and the physics just interprets its specific derivatives(velocity, force, and such). Thanks again Glenn.

yes. in the case of constant acceleration, there is not much point to using RK4, because in the 4 steps of the integration, you are just doing extra work to come up with the same constant dv

first of all, thanks for this article, its best explanation i´ve found so far.

but i have a problem understanding the RK4 integrator.

is it right that in your example the acceleration is changing over time?

i am working on a physics simulation and i have a constant acceleration.

so is dv the same in every derivative in case of constant acceleration?

there are special RK methods for ‘second order’ ODE’s.

like a = mF.I did both numerical analysis and theoretical physics and thereby know that it takes some effort to cram the two together, in fact i made a job out of it.

BTW , chapeau for your article

best regards

joop renes

you integrate both bodies together, ideally you have a single method that returns an array of forces to apply to the points — you cant update each body independently with RK4

cheers

I do have a question considering springs. If a spring is connected to two bodies and we are using rk4 to integrate, how would our acceleration/force method look like?

We are changing the state of only ONE body a a time, but the spring uses the position of both bodies to calculate the current force. Shouln’t we change the state of BOTH bodies in order to calculate the correct force at the specified point in time?

Two things.

First. Great article. clear, lucid, and technically accurate. Best explanation of Runge Kutta I’ve ever come across for programmers (not bad for mathematicians, either. Don’t get me wrong).

Second. As a teacher of technical topics (Software, hardware) I’m jealous of your style.

Thanks.

Take care,

Sam

Thankyou, a very clear explanation of RK4, you also saved me from looking silly as I have been saying “Yewler” instead of “Oiler”!

Regards Mike.

so is runge pronounced runj or run-ge?

I’ve only heard it pronounced “run-gah”.

yes, this parameter could be used to apply forces parameterized by time, eg. f = sin(t), however probably best to use doubles instead as float will lose precision significantly as the t value increases.

In the acceleration() function, you declare the parameter float t, but you never use it. Is this intended?

That’s because Hooke’s force, that is causing the spring to behave like a spring is not a function of time. The acceleration function uses float t if you somehow need time to calculate the acceleration. For example when you have a function of a = a(t). In case of a spring, the elastic force is calculated as a function of absolute deformation of a spring: F = -k*x. where k is value called “the Spring constant”.

(P.S. I hope I named all the physics stuff correct. My native language is russian and we have other names for these physical constants and stuff. )

Yes this is correct

Thank you for the valuable information.

Hey, thank you for providing actual code for RK4! I’m a little curious about RK4’s stability, because right now I’m using Euler @ 100Hz and it looks fine. If my simulation starts falling apart below 30Hz, would RK4 be stable down to 20Hz, or some sort of quantified value for the increased stability?

I was looking for Runge-Kutta Integration Techniques on the web and stumbled on this article. Thank you for this wonderful explanation.