Getting my head around RK4

Started by
6 comments, last by JimmyDeemo 13 years ago
So i want to add some physics to my little demo game, and after reading Glenn Fiedler's blog here, i decided to try and implement an RK4 integrator.

So i generally understand how RK4 work, how you get the derivatives from the current object state, applying forces at different points in the time step. Then you move on the simulation based on some kind of weighted average of those 4 derivatives.

My problem is understanding how external influences are applied to an object in the system, as Glenn Fiedler's example seems to be dependant on time rather on an object moving under the influence of external forces.

Let me give you an example, here is a function and description from the source code in Glenn's article titled 'Physics in 3D':

/// Evaluate derivative values for the physics state at future time t+dt
/// using the specified set of derivatives to advance dt seconds from the
/// specified physics state.

static Derivative evaluate(State state, float t, 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(state, t+dt, output.force, output.torque);
return output;
}

/// Calculate force and torque for physics state at time t.
/// Due to the way that the RK4 integrator works we need to calculate
/// force implicitly from state rather than explictly applying forces
/// to the rigid body once per update. This is because the RK4 achieves
/// its accuracy by detecting curvature in derivative values over the
/// timestep so we need our force values to supply the curvature.

static void forces(const State &state, float t, Vector &force, Vector &torque)
{
// attract towards origin

force = -10 * state.position;

// sine force to add some randomness to the motion

force.x += 10 * Mathematics::sin(t*0.9f + 0.5f);
force.y += 11 * Mathematics::sin(t*0.5f + 0.4f);
force.z += 12 * Mathematics::sin(t*0.7f + 0.9f);

// sine torque to get some spinning action

torque.x = 1.0f * Mathematics::sin(t*0.9f + 0.5f);
torque.y = 1.1f * Mathematics::sin(t*0.5f + 0.4f);
torque.z = 1.2f * Mathematics::sin(t*0.7f + 0.9f);

// damping torque so we dont spin too fast

torque -= 0.2f * state.angularVelocity;
}



From what i can gather 't is the time since the start of the simulation. In this function he uses the time to give a cube object motion. So firstly i want to ask, is this an integral part of the RK4 implementation? Does the force calculation of a derivative need to be derived from a total time value?

What i want to know is how do you use this model in a simple example such as gravity and say a ships thrusters acting against it?

My thought is that forces acting on an object would be stored as a variable in the State class, so you would add gravity as a force acting down each update, and the thrusters added to this store acting in the opposite direction on each update the thrusters are active. The forces function would then be modified to be interested in the time step rather than the total time. Giving you something that looks like this:

static void forces(const State &state, float dt, Vector &force, Vector &torque)
{
//Force of derivative is equal to the external forces over the time step.
force = state.totalExternalForce * dt;

//...omited the rotation for clarity...
}



Am i on the right lines here? Have i misunderstood something fundamental? Thanks for reading my gibberings, any help understanding this would be much appreciated.
Advertisement
Hi Jimmy,

I just wanted to check - is there any particular reason you're using RK4? Do you have a lot of pendulums or springs in your physics world? The reason i ask is that you can go a very very long way with simple euler integration, despite its inherent rubbishness... Lots of games (even physics ones) go into the shops with nothing more than simple euler integration :)

And its a hell of a lot easier to work with than RK4...

Cheers, Paul.

I just wanted to check - is there any particular reason you're using RK4?


In all honesty, its simply as a learning experience, i don't intend my little game to be particularly physics heavy at all. The code i am writing is mainly for portfolio/demo purposes and i though i would push the boat out a bit. :D In addition i am creating my own little library of stuff so i figured i'd write it as reusable.

After reading the articles on Glenn's blog it seemed that implementing an RK4 integrator might not be much more work anyway. I have done [color=#1C2837][size=2]euler before, but am i so far off the mark with this that i am best to stick with it?
[color="#1C2837"]

[color="#1C2837"]I also know it might be a good idea to simply use an existing framework, but as i said its for learning and Glenn's blog provided what seemed like a good explanation. In particular the rotation aspect of it, as this was something that confused me even in euler integration.


Anyone have any idea about this? Perhaps i am not explaining myself very well? Let me know if more info would help.
Hello,

You don't need to make your derivative function explicitly dependent on time. You have the option to, but you don't have to. Typically, games do not.

Your state however depends on time, in the sense that it evolves with time. But as long as you have your state for the wanted time, you don't need to input the actual time into your derivative function.


So, how implement RK4 in practice?

Wikipedia has a decent article on this.

Lets start with definitions:

y is your state. This is the complete state, including both positions and velocities. A point moving in three dimensions would have the state (px, py, pz, vx, vy, vz). This can be written more concisely as (p, v)

y' is the time derivative of the state. So y' = (p, v)' = (v, a). Since a = 1/m * f, the forces of the system needs to be calculated to evaluate y'.

f(t, y) give you y' for a given state at a given time. As I said before, you don't have to, and in this case shouldn't, have a direct dependence on time in this function. Just calculating the forces on the system given by the current state (y) is sufficient.

Anyhow, lets calculate the k-values of RK4:

k1 = f(t_n; y_n)

Here you need to evaluate the forces of the system at state y_n, and save the forces divided by masses in k1 and the current velocities.

k2 = f(t_n + 0.5h; y_n + 0.5*h*k1)

Here you first need to create a temporary state, y_n + 0.5*h*k1. Evaluate the forces for this temporary state to obtain k2.

k3 = f(t_n + 0.5h; y_n + 0.5*h*k2)

Here you need another temporary state, y_n + 0.5*h*k2. Evaluate forces.

k4 = f(t_n + h; y_n + h*k3)

I guess you can figure this out :)


Finally, y_{n+1} is given by:

y_{n+1} = y_n + h / 6 * (k1 + 2k2 + 2k3 + k4)


When it comes to objects interacting with each other, like collision handling etc, you need to evaluate this every time you evaluate the force for a given state. I.e., in your temporary state needed when calculating k3 two objects might collide with each other while not colliding with each others in the temporary state needed to calculate k4. Collision forces should be generate in the first case and not the other. (given that a force based collision avoidance scheme is used)


Thats it basically, given that you have the ability to create these temporary states and calculate forces from them, implementing RK4 is not much more difficult than Euler integration (but lots more code :) ). It is however way more computationally expensive to use, which is probably the reason it's not commonly used in commercial games. But it is a fun exercise to implement!


Hope that clear things up for you.
//Emil
Emil Jonssonvild
Hi there,

Imho the RK4 algorithm has gotten more attention than it deserves. Despite that it seems to preserve energy very well, it is in fact not a symplectic algorithm, meaning that energy will change over time (decrease towards zero in this case). In stead, take a look at David Wysongs symplectic 4th, 6th, and 8th order algorithms. They give better result with less cpu-cycles and are easier to implement than RK4.

http://read.pudn.com/downloads72/sourcecode/math/261769/symplectic.cpp__.htm

Cheers,
Mike
I simulated a cannonball shot with both Euler and RK4. The high start velocity leads to a very high air resistance (rapid deacceleration), and with Euler
you need very small timesteps to be accurate. Just as an example, 4 (or 5) RK4 steps were more accurate than 100 Euler steps.
RK4 is basically 4 Euler type of steps, so if we assume it's 4 times as cpu extensive, then atleast in that case, RK4 was more effective.
Hi all,

I've not actually had chance to come back to this since my last post, but will do soon.

@cannonicus Thanks for you detailed post. I'm not sure that you understood where my confusion comes from. I am quite happy with how the integrator works, and its ins an outs. I am confused about applying forces to an object so i can actually make use of it.

Perhaps it might be best to deal with this by example. I have an object that i want gravity to act upon. When using Euler integration i would apply this force each update, solve and then integrate to give me the new position of the object.

Where does this application of force fit within the RK4 integrator? I have to put it all in one function right? The equivalent of [color="#1C2837"]Glenn Fiedler's static void forces() function. Is it simply a case of storing external influences on an object and then solving them in each evaluation step?

This topic is closed to new replies.

Advertisement