Simulating Newtonian Particles

CS 493/693 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 Newtonian 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 abuse notation and call this timestep "dt".

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

   change in P = dt*V
   change in V = dt*A

Again, this usually works, but 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 somewhat wrong!

Acceleration is related to the net force and mass according to:
    F = m A
which is more useful to rearrange as:
    A = F/m

For gravitational forces:
    F = G m₁ m₂ / r²

G is the gravitational constant. m₁ is the mass of the particle. m₂ is the mass of the Earth. r is the distance between the centers of the objects.

In computing the acceleration of the particle, the mass of the particle m₁ just cancels out, leaving you with a constant acceleration:
   A = G m₂ / r² = 9.8 m/s² near Earth's surface

Wind Resistance

One beautiful thing about simulations is that it's usually easy to add new terms. For example, we can add wind resistance by just applying a new force in the direction opposite the particle's velocity:
F = -k V

You can download the spreadsheet simulator we built in class. One thing to notice is how very high wind resistance at large dt values causes the simulation to "overcorrect": it pushes back so hard against the existing velocity, that it actually reverses the velocity's sign. Even larger wind resistance constants, unlike actual hyperviscous situations like the Queensland Pitch Drop Experiment, actually add energy to the system, with the velocity reversing sign and increasing magnitude with each step. 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.

For example, dt=0.1 and wind-drag=-21 exceeds 100km/s within ten seconds. This is nonphysical, but surprisingly common, especially when you're pushing the simulation timestep upward, to improve performance.