Simulation of Physical Phenomena

CS 493 Lecture, Dr. Lawlor

In a class about simulations, you'd be correct in thinking we're going to design and analyze quite a few simulations.  These are going to differ in really fundamental ways, but there are many broad concepts that show up repeatedly in this field.

Units

Physical quantities have units, like meters or seconds or watts.  Most programming languages don't know or care about the units, so you can easily make unit errors like:
    float energy_in_joules = distance_in_feet * force_in_pounds;

Of course, foot * pound force gives you ft-lb of energy, not joules.  But all the variables are "float", so the compiler happily lets these errors creep in. 

The best way to avoid unit errors is to pick your units when you begin building the simulation, and stick with them throughout.  Metric units (meters, kilograms, seconds) seem to be the most common, although I've worked with everything from astrophysical simulations where "1.0" is the size of the universe, and molecular dynamics simulations where "1.0" is the size of an atom.

It is possible to go back and figure out the effective units you're working in after the fact, but it's less effort and more reliable to figure them out as you go.  Sometimes I'll work without units just to try to capture the overall "feel" of a phenomena before actually matching it to reality.

Differential Equations

The usual definition of velocity is the time derivative of position:
    V = dP/dt

This is extremely mathematically useful, since we can use all the tools of calculus on a differential definition like this.  But it's less computationally useful, since computers don't do infinitesimals.  There are a variety of ways to discretize this equation into actual numerical steps; but the simplest is to replace the infinitesimal derivative with a finite difference:
    V = delta_P / delta_t
which means
    V = (new_P - old_P) / delta_t

Given a set of positions, this equation lets us estimate the velocity between known positions. We can also solve for new_P, to get the Euler update:
    new_P = old_P + V * delta_t

There are many variations on this method with better properties, but this is the basic idea.  One big problem here is to best approximate the underlying partial differential equation, "delta_t" should be as small as possible.  Small enough, and roundoff begins to change your simulation results.  Also, making the timestep half as big needs twice as many steps for the same simulated time, so you quickly run into computational limitations: even a simple partial differential equation can consume arbitrary computational power and still not be perfect!

Conserved Quantities, like Energy

Energy is conserved mostly because:

  1. If you keep adding energy, everything explodes (the world "ends in fire").
  2. If you keep removing energy, everything freezes in place (the world "ends in ice").
  3. If you don't add or remove energy, everything keeps working reasonably.
Many simulations calculate and monitor the total energy of everything in the simulation, as a way to verify that things are working correctly.  This is a good way to detect even surprisingly subtle flaws in the physics or numerics.  Be aware that conserving energy is a goal, not an absolute requirement, and often minor conservation errors from things like roundoff are difficult to avoid. 

You can also lose or gain energy at the borders of the simulation, where things like the energy added to the wind as your object passes by, or the electrical energy input to the motor, aren't accounted for.  You can often fairly easily add these border energy terms back in, and sometimes they can be quite useful (e.g., noting that most of the drag energy loss in a boat is happening at the tail).

Linear and angular momentum are also conserved quantities.  Loss of momentum conservation tends to look weird, causing objects to fly off into the distance.

Boundary Conditions

Often the trickiest part of a simulation is not "knowing how to simulate", it's "knowing when to stop".  Since our computers are finite, we need to do something at the borders of the domain that we can afford to simulate.  That "something" can be very simple, but there are times when a simple boundary condition messes up the simulation in the interior.  Simulating waves or fluid flow is especially tricky, since it's easy to get wind resistance or wave reflection off the boundary.

One common trick is to make the boundaries periodic, essentially wrapping the simulation around at the boundaries.