Chris Hecker wrote an awesome series of articles on physics simulation in 1996, and they're available on his site.

I've been doing some sorts of simulation since I was very young (QBasic years). Let's say you're trying to simulate a ball with gravity. You keep track of the ball's x and y position, and x and y velocity. Every frame, you add the velocity to the position to get the new position. Then you add some constant to the y velocity to simulate gravity.

When the ball hits something (like the ground, or a wall), you have to bounce it. To do this you can just invert the component of the velocity according to the wall it bounced off (for ground/ceiling invert the y component, and for side walls, the x component).

But things are always tricky. Sometimes, you're not careful, and the ball goes right into the floor and instead of bouncing, it gets stuck, because it's always colliding every frame so you keep inverting the velocities and it never goes anywhere. Or you discover that the ball isn't bouncing quite perfectly right. Or you try to make a billiards simulation, but the pool balls end up interpenetrating too far and your balls bounce at really odd angles. These things always frustrated me and I never knew how to deal with them well until I read Chris's papers.

One thing I did figure out in order to make my simulation smoother was to add a time parameter to the simulation. Every frame, you take the amount of real-world time that passed since the last frame and use that as the value for delta-T. Then, instead of adding the velocity itself to the position, you multiply by the delta-T, so that no matter what interrupts your program or makes it slower (more objects on screen), the simulation makes objects move across the screen at the same real-world speed. But this also came in useful later.

Anyway, I wanted to make a better physics simulator than I had done as a kid, so I dug up the old Chris Hecker. The idea in the above paper series I always liked (and you can read it in part 3) is the idea of subdividing time in order to find the exact point that two objects begin colliding, so you can do the nice collision response. By "exact" here, I really mean "to the desired epsilon," which is not actually exact at all but it's a hell of a lot better than the simplistic method I learned as a kid.

The way time subdividing works is that you have three outputs from your collision detector: "No collision", "Interpenetrating" and "Colliding". When simulating, you advance the positions of the objects, then check for collisions. If there's no collision, you go onto the next frame; if they're interpenetrating, you went too far -- back up. Go back to the previous gamestate (where, presumably, you knew all the objects were in valid positions) and then set delta-T to be half of what it was before. Advance the position again -- this time things only move by half as much. Test for collision again in the same way.

When you actually see the result "Colliding," it means that you've subdivided time enough (moved the objects in tiny enough increments) for the collision detector to say that two objects are touching, but not penetrating. This is where the desired epsilon comes in. Your collision detector will return "colliding" when the objects are not quite touching but they're within the tiny epsilon of doing so. At this time, you apply collision response as if they were actually colliding, and expect that this estimate is close enough to the correct answer. And most of the time, it is.

Once you apply that collision response, you can continue simulating the gamestate normally. You have to finish out that simulation frame, so you step forward in time until you get back to the original step size. If the collision response worked (and nothing else is colliding), the collision tester will return "No collision" and you're done with the frame.

I decided to implement it in Haskell, and the time subdivision algorithm came out really beautifully. I picked Haskell because I figured it would be fun to write collision detection and response methods using a functional style: accept a GameState and return "Colliding" or "Not Colliding, here's your new GameState". And it would be easy to return to a previous gamestate. So that's what I did.

Here's the type signature of my collision detection/response function, as well as the game state it operates on:

data GameState = GameState {
          t :: Double                                       -- The current time
          ,objs :: [((Double, Double), (Double, Double))]   -- List of objects with (x,y) and (xvel, yvel)
          }

collides :: GameState -> Maybe GameState

The spec for this function is: If the state is interpenetrating, it returns Nothing. Otherwise, if the state is colliding, apply the appropriate collision response and return the new game state. Otherwise (no collision), just return the same game state.

It's also worth looking at the (very simple) advance function, which just advances the position of every object by its velocity times the delta-T:

advance :: Double -> GameState -> GameState
advance dt (GameState t objs) = GameState (t+dt) $ map (\ ((x,y), (xv, yv)) -> ((x+xv*dt, y+yv*dt), (xv,yv))) objs

Now, you can appreciate how easy and beautiful it was to write the simulation function, which accepts a gamestate and a delta-T, and (using the given 'collides' function), simulates delta-T worth of time using the physics:

simulate :: Double -> GameState -> GameState
simulate dt gs =
        let gs_dt = advance dt gs in
        case (collides gs_dt) of
                Just gs' -> gs'
                Nothing ->      -- Break it down.
                        let gs' = simulate (dt / 2) gs
                        in        simulate (dt / 2) gs'

As you can see, it cleanly, recursively implements the above time-subdivision algorithm. "Try to advance the whole time step. If you can't, advance it halfway, and then advance it halfway again." It was a bit mind-blowing the first time I tried to write it. "It can't possibly be that simple!" But it actually is that simple, and trying to analyze the runtime of this function can help explain its workings.

What are the cases? If no collision or colliding, the function returns immediately (Just gs' -> gs'). So that's constant time (pretending 'advance' and 'collides' run in constant time).

If interpenetrating, the function calls itself recursively twice. So it's tree-recursive. But (if there's only one collision to deal with), one of the recursive calls will always return immediately. Additionally, once dt is small enough, BOTH calls will return immediately.

Worst case: you keep calling simulate until dt is on the order of epsilon. Since you're dividing by two each time you recur, the number of calls to simulate to find one collision is about O(-log(epsilon)) -- if epsilon is a small number like .00001, -log(epsilon) is 5.

This algorithm handles gamestates with multiple collisions. Simulate will end up subdividing time until it resolves the first one, then it will do it for the second, and so on until they're all resolved. So the actual runtime is O(c * -log(epsilon)) where c is the number of different collisions to be resolved in advance gs.