Stable Simulation: Avoiding Overshoot

CS 493 Lecture, Dr. Lawlor

Almost all interesting simulations involve some sort of closed-loop feedback, where a change in A causes a change in B which in turn affects A again.

For example, consider the velocity of an object including drag forces from air resistance.
    F = -k * V
    A = F/m  (from F = mA)
    dV / dt = A

Here F is the net force on the object, V is the object's velocity, k is the object's wind resistance (newtons of resistance per meter per second of velocity), A is the object's acceleration, and m is the mass.

If drag is the only force involved, and discretizing to first order, we get:
    dV = (-k*dt/m) * V
    V = V + (-k*dt/m) * V = V + speed * V
where we've defined speed = (-k*dt/m).

That is, the change in V itself varies with V.  If that speed term (-k*dt/m) is positive, then V will feed back on itself leading to exponential growth.  Luckily, speed is always negative because real drag constants, timesteps, and masses are all positive.  In the theoretical non-discretized equations, a negative speed term causes V to drop exponentially, though possibly with a managably small exponent.  A small negative speed term causes the discrete V to also drop exponentially.

But weird stuff starts happening if the speed term gets negative enough:
Anytime speed<-2.0, V is going to spiral out of control fairly quickly.  This is similar to an inexperienced driver skidding on ice: they steer too hard, overcorrecting the car's orientation, and end up making the skid worse instead of better.  The effect is that the simulation "blows up", where the points have vanished off the screen, or hit floating-point infinity or not-a-number.  Not good.

Note that speed = -k*dt/m, so speed is going to get out of control if:
Many simulations have a stability limitation on the timestep.  If your steps get bigger than this, you can expect the simulation to explode.

Stability in Spring Systems

Try out this spring system stability spreadsheet--notice that:
Generally, if the spring forces can push the mass to complete a cycle in less than pi timesteps, the mass starts overshooting and the simulation blows up.  Here, a simple harmonic oscillator has period 2*pi*sqrt(m/k).  I've seen programs that compute their own timestep by figuring out the minimum stable period at startup.

Note also that if you update P and V in the wrong order, the system is unconditionally unstable.  You can see this by updating cell B5 to "=B4+A$2*C4" (using the *old* velocity, instead of the new updated velocity), and choosing "fill down".  Note it now always blows up.  Using the *new* velocity (C5) makes this "Leapfrog integration"; using the old velocity (C4) makes it the less accurate and unstable Euler.

The Speed of Sound

The other big limit on the timestep is the speed of sound.  If I take a steel bar, and break it into 100 small springs, the only way a force exerted on one end will eventually move the other end is via each of the 100 springs.  This must take at least 100 timesteps, since each spring computes its forces based on the *old* position of each neighbor.

The problem is the speed of sound in solid materials is very fast; for example, the speed of sound in steel is about 6 km/sec (just under orbital velocity!).  Thus, sound crosses a 1m steel bar in 1/6,000 of a second.  To simulate 100 smaller springs, we need 600,000 physics updates per second, which is usually way too many.

One saving grace for computer graphics is that it's not easy to see a slower-than-correct speed of sound. And in 3D, using larger (hence faster-reacting) tets inside an object can help reduce this effect.  But the net result is that accurate physics typically takes at least several physics steps per animation frame.