As a realtime 3D graphics developer, I need to wage many battles. I fight with artists over polygon counts, with graphics card manufacturers over incomplete or incorrect drivers, and with some producers’ tendencies to continuously expand feature lists. However, some of the greatest battles I have fought have been with myself. I fight to bring back the knowledge I have long since forgotten. I fight my desire to play the latest action game when more pressing needs are at hand (deadlines, the semblance of a social life).
In this I document one of the less glamorous battles — the battle of the physics simulator. It’s not going to be fun. It’s going to be a bit bloody. However, if I ever hope to achieve a realistic and interesting physics simulation, it’s a battle that must be fought. So, my brave warriors, join me. Sharpen your pencils, stock your firstaid kit with plenty of aspirin, drag out the calculus book, and fire up the coffeepot. Let’s get started.
I hope you all had a chance to play around with the soft body dynamics simulator from my article titled Collision Response: Bouncy, Trouncy Fun. The demo highlighted an interesting problem — the need for stability. While creating my dynamics simulation, I waged a constant battle for stability. However, in order to wage the war effectively, I need to understand the roots of the instability in the system. Last month, I implied that the problem resulted from my use of a simple Euler integrator. But I didn’t really explain why that caused the problem. Let me fix that right now.
Integrators and You
Many game programmers never realize that when they create the physics model for their game, they are using differential equations. One of my first programs on the Apple II was a spaceship flying around the screen. My "physics" loop looked like this:
ShipPosition = ShipPosition + ShipVelocity;
ShipVelocity = ShipVelocity + ShipAcceleration;
Look familiar to anyone? It’s a pretty simple physics model, but it turns out that even here I was integrating. If you look at the Euler integrator from last month, I had
Position = Position + (DeltaTime * Velocity);
Velocity = Velocity + (DeltaTime * Force * OneOverMass);
Now for my simple physics model, DeltaTime = 1 and Mass = 1. Guess what? I was integrating with Euler’s method and didn’t even know it. If I had made this Apple II physics model any more complex, this integrator could have blown up on me. These sorts of problems can be difficult to track down, so it’s important to understand the causes.
When Things Go Wrong
The reason that the Euler integrator can blow up is that it’s an approximation. I’m trying to solve a differential equation by using an iterative numerical method. The approximation can differ from the true value and cause error. When this error gets too large, the simulation can fail. A concrete example may help to explain. Last month, I added a viscous drag force to the simulation to add stability. The formula for this force was
(Eq. 1)
In this formula, kd represents the coefficient of drag that is multiplied by the velocity of the particle. This coefficient determines how fast the velocity of the object is dragged down to zero. This is a very simple differential equation. In fact, it’s simple enough to be satisfied for v directly by the formula. I can use this exact solution to check the accuracy of my numerical integrator:
(Eq. 2)
Euler’s method is used to approximate the integral curve of Equation 2 with a series of line segments along this path. Each step along this path is taken every time, interval h, via the formula
(Eq. 3)
In all cases, the viscous drag force should approach zero. However, the size of the step h and coefficient of drag kd determine how well the approximation performs. Take a look at Figure 1.

Figure
1. A decent approximation.

With the given step size and drag coefficient, Euler’s method may not be a great approximation, but it gives the desired result. The velocity converges on zero. But take a look at the relationship between the step size and drag coefficient in Equation 3.
Ifthen the approximation step will overshoot zero, as you can see in Figure 2.

Figure
2. This looks a lot worse.

By increasing the step size, I was trying to get a system that converged to zero more quickly — but I got something entirely different. Things really start to get bad when the drag coefficient increases more, as in Figure 3. As each step is taken, not only does the approximation oscillate across zero, but it also actually diverges from zero, and eventually explodes the system. This is exactly what was happening in the spring demonstration from last month, when the box blew up.

Figure
03. Kaboom!

How Can I Prevent Explosions?
If you find a situation where your simulator blows up, there’s an easy way to see if this kind of numerical instability is the cause. Reduce the step size. If you reduce the size of the step and the simulation works, then this numerical instability is the problem.
The easy solution is always to take small steps. However, realize that each step requires quite a few calculations. The simulation will run faster if it can take fairly large step sizes. Unfortunately, when you get lots of objects interacting, these instability problems appear even more. So, just when things start to get interesting, you need to reduce the step size and slow things down.
I’d rather create an integrator that would allow me to take large step sizes without sacrificing stability. To do this, I need to look at the origins of Euler’s method.