Coursewide Review for Final Exam

CS 480 Lecture, Dr. Lawlor

Over the past months, we've looked at several very different ways to simulate a range of real-world phenomena. Our simulations have been very simple, but I claim they capture the sorts of simulations really used by pure scientists, applied engineers, and game and movie developers.

Specifically, simulating physical phenomena boils down to three (interrelated) steps:

Pick a data structure to represent your problem's data. Often, this is a set of discrete objects like particles or triangles, or a continuum field like a 2D or 3D texture.
Pick a mathematical model to represent your problem's behavior. This is usually a partial differential equation, at least in all the literature generated since the time of Newton and Leibnitz.
Figure out how to apply your mathematical model to your data structure. This usually means discretizing your partial differential equation in space (taking steps equal to your data structure's size) and time (taking small time steps, like a few milliseconds). You can usually evaluate spatial derivatives by looking at your neighbors, and thus figure out how your value changes with time (see below).

There's often also a step 4, one of "figure out how to keep it stable", "figure out how to make it run fast", or "figure out why it's not simulating what you want it to simulate".

Data Structures

Particles. Separate disconnected entities. Easy to store arbitrary data on each particle. Need a spatial search data structure like a tree or grid if the particles interact with anything spatially (e.g., other particles, grids, etc). Typical representation for critters, fireworks, dust, etc.
Unstructured meshes, like triangle or tetrahedron meshes. They're basically a set of particles with new force-generating entities between fixed neighboring points. You normally generate the meshes using a separate special program or library, not part of the main simulation. Typical representation for solid objects.
Grids, like 2D or 3D textures. Grids are divided into pixels or voxels that each represent a small part of the problem's space. Typical representation for flowing fields like liquids or gasses.

Mathematical Models

Newton's Equations

    F = mA
    dV/dt = A
    dP/dt = V
    F: net force (vector, usually computed from something interesting!)
    A, V, P: acceleration, velocity, and position (vectors)
    m: mass (scalar)

Blurring model: du / dt = k * d² u / dx²
Navier-Stokes
σ=Ee (stress is stiffness time strain)

Applying Models to Data Structures

Generally, the way you compute discrete values from a partial differential equation model is to substitute in the Taylor series. Here's the Taylor Series for a function P that changes with time (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' (dP/dt), but also on the current acceleration P'' (d²P/dt²), 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.) This mathematics is *exactly* the formula we use to compute the new P:
    P = P + dt*V

For derivatives in space, like dP/dx, you can expand the Taylor series in position:
    P(x+h) = P(x) + h * dP/dx (x) + h²*d²P/dx²(x)/2! + ...

Again, to zero'th order,
    P(x+h) = P(x)
which means dP/dx isn't involved (hmm... not very useful for computing dP/dx!).

To first order,
    P(x+h) = P(x) + h * dP/dx (x)
Unlike the derivatives in time, where we knew the derivative, but not the left hand side, now we know the left hand side (we know P everywhere), but we want to find the derivative. So we rearrange to the well-known one-sided "finite difference equation":
    (P(x+h) - P(x)) / h = dP/dx (x)

So, the value of our x derivative dP/dx is proportional to the difference between our right neighbor's value P(x+h), and our value P(x).

We can substitute in a negative value for h to get:
    (P(x) - P(x-h)) / h = dP/dx (x)
This is the difference between us and our left neighbor. If you average both of these estimates, you get:
    (P(x+h) - P(x-h)) / 2h = dP/dx
This is the same "right minus left" we've written perhaps fifteen times in class!

Similarly, the value of our y derivative is this same finite difference, but measured against our top and bottom neighbors. Again, we can actually compute any derivative in space or time using the Taylor series, and even (with a bit of practice) figure out what partial differential equation we're actually solving given a chunk of code!