Discretizing Differentials: Numerical Integration
CS 480 Lecture,
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
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
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) + h2*P''(t)/2! + h3*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) + h2*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.