Fourth order Runge-Kutta (RK4) instead of Euler integrator for 3D physics (Quaternion)

Questions
1

My port of the rigid body simulator from “Physics for Game Developers” (boat / plane / vehicle in fluid) starts to oscillate when a certain value gets very close to zero (steady state) and blows up. A smaller dt helps a little - I assume the Euler integrator to be the problem. So I’d like to try a different integrator. Unfortunately, the book does not show other implementations. “RK4” seems to be popular. There are some general implementations but to my embarrassment I lack the math and physics skills to apply them to the 3D simulation. So I need help. :confused:

I’ve found Integration Basics | Gaffer On Games that shows how to implement the algorithm for scalars (1D) and a spring-damper system example.
I had hope that I’d be able to implement the linear motion to get started (vectors only).

However, a key function looks like this:

    float acceleration( const State & state, double t )
    {
        const float k = 15.0f;
        const float b = 0.1f;
        return -k * state.x - b * state.v;
    }

Note it uses only the state (x:position, v:velocity) to determine the acceleration, but not t.
The calling function passed t+dt from its own parameters but does not explain where t originates. I assume it’s just there for completeness’s sake and would most likely be zero, since the simulation uses relative time dt.

In my simulation, acceleration is determined by constant mass and some forces (in addition to thrust) that depend on orientation and speed (drag, lift).

I don’t understand how I could calculate this for future steps, while I’m in the middle of calculating it for the next step.
I currently simply use the acceleration vector obtained previously, but I’m pretty sure that’s wrong.

Here is the JSFiddle:

I thought my initial attempt at the linear motion RK4 was somewhat working so I could post it, but it’s not doing anything plausible. Once it will, I’ll provide a fiddle.

Thinking ahead, the next step would be to implement it for the angular motion (inertia Matrix, rotation Quaternion) - I have no idea how :cry:.

Any help appreciated!

1 Like
2

No progress yet, but I’d like to share Physics in 3D | Gaffer On Games where the concepts introduced in the first link are extended to 3D, although only little code is provided.

1 Like
3

Here is the JSFiddle using my Vector3 extension of the scalar RK4 code from Glenn Fiedler / Gaffer On Games:

It does something, but it can’t be right.
As mentioned before, I don’t understand how to get acceleration(state, t+dt ).
So currently, it uses the acceleration calculated in the current step (thrust, lift+drag from speed) for all of the RK4 evaluations…
With the Euler integrator, acceleration gets calculated at every step, based on the previously integrated state.
So at best, the acceleration computed in last step could be used for t and the newly computed could be used for t+dt?
What am I missing?

1 Like
4

In “Physics in 3D”, it is suggested to “drive the simulation” using forces instead of acceleration.
Another method is suggested for ease of use with rotations.
Using Euler, it’s easy to get it working:

    // EULER 2
    // integrate force to momentum
    this.vMomentum.add(this.vForces.clone().multiplyScalar(dt));  // dp/dt = F
    // convert momentum to velocity by dividing it by mass
    this.vVelocity = this.vMomentum.divideScalar(this.fMass); // v = p/m
    // integrate velocity to get position
    this.vPosition.add(this.vVelocity.clone().multiplyScalar(dt));  // dx/dt = v

That means that the RK4 code needs to change, and instead of acceleration() it will probably require force() - but the problem remains (function of time).

1 Like
5

Runge-Kutta is not terribly complicated, but you do kind of need to know what a time integration function is actually doing numerically. Essentially you have a known state (list of variables describing your system), and you have a rate of change for those variables (derivatives of the variables with respect to time). The rate of change of position is velocity, the rate of change of velocity is acceleration, etc. etc.

For many simple cases with small time steps, Forward Euler works just fine:

x[n+1] = x[n] + v[n] * delta_t // Linear extrapolation of position
v[n+1] = v[n] + a[n] * delta_t // Linear extrapolation of velocity

where n is the index of the current state. Assuming constant dt, the elapsed simulation time would be t = n * dt. Since v = dx/dt is how much x is currently changing with respect to time, we multiply it by the current time step (dt) to get how much x will change over the interval of dt. NOTE: This assumes that v is constant over this entire interval. The same holds true for updating the v using a = dv/dt. But if a is constant over this interval, then v will be changing over this interval, thus invalidating our previous assumption that v is constant over the interval.

A slightly more accurate forward Euler integration may look like this:

x[n+1] = x[n] + dt * (v[n] + 0.5 * dt * a[n]) // Quadratic extrapolation of position
v[n+1] = v[n] + dt * a[n] // Linear extrapolation of velocity

Most of the time, this will be sufficient to get the job done, but you have to accept the fact that this is likely only a rough approximation, since your acceleration is probably changing over the interval as well. For example, if I were doing orbital dynamics simulations with this integrator, the acceleration due to gravitational forces is a function of the relative positions of all of the bodies in the system. Since these positions are constantly changing in time, then the acceleration of each body is also constantly changing. In these situations, the above approximations will gradually diverge from the exact solution as more and more errors accumulate in the positions as a result of the constant acceleration assumption.

The Runga-Kutta family of time integrators improves the situation somewhat by sampling the acceleration at multiple points over the interval and incorporating that information into the updates of the position and velocity. For the simplest RK2 algorithm (second order accurate):

x[n,0] = x[n] // position at the beginning of the time interval
v[n,0] = v[n] // velocity at the beginning of the interval
// NOTE: F is a forcing function computed from the current state variables
a[n,0] = F(x[n,0],v[n,0]) / m // acceleration at the beginning of the interval 

// Update the state variables using linear extrapolation to t = t[n] + alpha*dt
// for some 0 < alpha <= 1.
alpha = 1.0 // NOTE: This may vary for different RK2 implementations
x[n,1] = x[n,0] + alpha * dt * v[n] // updated position
v[n,1] = v[n,0] + alpha * dt * a[n] // updated velocity
a[n,1] = F(x[n,1],v[n,1]) / m // Acceleration recomputed using the updated states in X[n,1]

// Refine the update using a combination of the previous states and accelerations
x[n+1] = x[n,0] + 0.5*(v[n,0] + v[n,1]) * dt // Refined position update using an average of the velocities at the ends of the interval
v[n+1] = v[n,0] + 0.5*(a[n,0] + a[n,1]) * dt // Refined velocity update using an average of the accelerations at the ends of the interval

Higher order Runga-Kutta methods operate on the same basic principle, but they use more samples of the state variables and their derivatives across the interval to obtain better and better estimates of their (potentially non-linear) behaviors.

BTW: The code in your first JSFiddle is likely blowing up due to a divide by zero error. Whenever you use division or vector normalization (dividing a vector by its length), you need to be sure that the denominator does not go to zero. If it does, you will usually get a NaN, but you may also get undefined or otherwise erratic behavior.

2 Likes
6

Hi and thanks for your help! I’ll probably give it another try in the next few weeks.