Discretizing Differentials: Numerical Integration

CS 480 Lecture, Dr. Lawlor

Tons of physical phenomena can be described by a "partial differential equation", an equation involving derivatives. For example, velocity V is the derivative of position P with respect to time t:

    V = dP / dt

Similarly, acceleration A is the time derivative of velocity:

    A = dV / dt

Basic physics, of course. But this is all continuum mathematics: in a real computer simulation, we can't take infinitesimal steps dt, so we have to pick some finite discrete timestep. We'll call this timestep h (instead of, say, delta t).

It's easy but inaccurate to discretize the above equations by algebraically substituting the infinitesimal step dt for the finite timestep h:

    h * V = change in P
    h * A = change in V

Again, this usually works, but it's also usually somewhat wrong--that is, the error in this approximation is nonzero, sometimes enough to screw up your simulation. To see why this is wrong, ask yourself: should I first update P (using the old V), then update V, or vice versa (using the new V)? When dt was infinitesimal, this didn't matter; but for any nonzero h, you'll get a slightly different answer if you update P then V versus V then P. When you can't tell which of two things is the right answer, it's a good sign both are wrong!

Higher Order Integration

The problem with our discretization above is clearer if you substitute in the Taylor Series for P (again, see the MathWorld article, or Wikipedia):
    P (t+h) = P(t) + h*P'(t) + h²*P''(t)/2! + h³*P'''(t)/3! + ...

That is, the value of P at some future time t+h really depends not just on the current position P and the current velocity P', but also on the current acceleration P'', and all the higher derivatives.

To "zeroth order", we can ignore all the h-dependent terms:
    P(t+h) = P(t)
That is, we can assume objects don't move during our timestep. (Say, stationary objects!)

To "first order", we can ignore all the powers of h higher than one:
    P(t+h) = P(t) + h*P'(t) = P + h*V
That is, we're assuming the object moves in a straight line during our timestep. (Say, no gravity, no friction.)

To "second order", we can ignore all the h powers higher than two:
    P(t+h) = P(t) + h*P'(t) + h²*P''(t)/2! = P + h*V + h*h*A/2
That is, we're assuming the object moves in a parabolic arc during our timestep. (Say, constant gravity without friction.)

And of course you can trivially substitute in yet higher orders straight from the Taylor Series, assuming you know or can estimate the corresponding derivatives. There's even a cool trick where you keep one old timestep called Verlet Integration (read it!).

As it turns out, higher order integration can substantially improve your simulation's bottom line--you get bigger timesteps (hence a faster simulation), more accurate results, or both!

There is a flip side to higher order integration, though--you often have to keep more history data, compute more derivatives, or otherwise compute more flops. Also, low-order integration can be less problematic around discontiuities, which can cause ringing or other strangeness with a higher order integrator. I've seen research simulation codes that used up to fourth order integration, but most codes use second or third order integration, with a few basic codes (like most of mine!) as low as first order integration.