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.
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.