Game Physics

In my first post I talked about two different ways of integrating the laws of motion, Euler integration and Verlet integration. There I mentioned that Verlet integration, although it looks very similar to Euler integration at a first glance, is actually better. I thought I’d take the time to show just how much better with two simple examples.

An object falling under the influence of gravity

Our first example is probably one of the first laws of physics you learn in school, how an object moves when falling under the influence of gravity. To simplify things, we will set the value of the gravitational acceleration to 10 instead of the actual value of about 9.82:

\(a = 10\ \text{m/s}^2\).

Say we drop an object from the height of 500 meters. We start from rest, so the initial velocity at time \(t=0\ \text{s}\) is \(v_0 = 0\ \text{m/s}\) and the initial position is \(y_0=500\ \text{m}\). We can quite easily find the exact path of the object, since the acceleration is constant. An object with constant acceleration moves according to the formula

\(y(t) = y_0 + v_0t + \frac{\displaystyle 1}{\displaystyle 2}a_0t^2\).

In our case we get

\((1)\qquad\qquad y(t) = 500 + 0\cdot t + \frac{1}{2}\cdot (-10)\cdot t^2 = 500 – 5\ t^2,\)

which is the actual path the object takes. Let’s see what we get when using Euler integration to compute the path of the object. For simplicity, we set the time \(\Delta t\) between the steps to 10. So, let’s where the object will be after 10 seconds:

time	position	velocity
0	500	0
1	500	10
2	490	20
3	470	30
4	440	40
5	400	50
6	350	60
7	290	70
8	220	80
9	140	90
10	50	100

So, after 10 seconds the object will be 50 meters above ground, if it’s position is computed using Euler integration with time steps of 1 second. Let’s compare this with the actual position the object will have:

\(y(10) = 500 – 0.5\cdot 10\cdot 10^2 = 500 – 0.5\cdot 1000 = 500 – 500 = 0\).

After 10 seconds the object would actually hit the ground, not be 50 meters above ground. How does Verlet integration measure up? Well, as I mentioned in my previous post, Verlet integration is exact as long as the acceleration is constant. You can have a look here for a nice derivation of Verlet integration by Jonathan Dummer where this fact is easy to see.

One important thing to note when using Verlet integration is that you need to be careful when setting up the initial conditions. With Euler integration, it’s fairly straightforward. In the above example, the object starts at position \(y_0 = 500\ \text{m}\) and since we’re dropping the object it’s initial velocity is \(v_0 = 0\ \text{m/s}^2\). With Verlet integration, we need to know the current and previous position of the object. The current position we know, it’s \(y_0 = 500\ \text{m}\). But what about the previous position \(y_{-1}\)? For Verlet integration to be accurate, we need to set \(y_{-1}\) to be the physically accurate position, assuming the current physically accurate situation, which is that the object is being accelerated downwards at 10 m/s². So we use formula (1) above, and get

\(y_{-1} = y(-1) = 500 – 5\cdot(-1)^2 = 500 – 5 = 495\).

In the graph below I have plotted the physically correct path (black), the path computed using Euler integration (blue), and the path computed using Verlet integration (red).

Harmonic oscillation

So what about a situation where the acceleration is not constant, and hence Verlet integration is not exact. How much better is it than Euler integration then? Let’s look at one of the simplest such situations, undamped harmonic oscillation. Here, the force on the object is proportional to its position:

\( (2)\qquad\qquad F = -ky\),

where \(k\) is a constant that determines the ‘stiffness’. For simplicity, lets say the object has mass 1 (ignoring units for now), and the constant is

\( \displaystyle k = \frac{\pi^2}{4}\).

Since \(F = ma\) and \(m = 1\), equation (2) becomes

\( (3)\qquad\qquad \displaystyle a(t) = -\frac{\pi^2}{4}y(t)\)

Assuming that the initial conditions are \(y_0 = 1\) and \(v_0 = 0\), we get the exact solution

\(\displaystyle y(t) = \cos\bigl(\frac{\pi}{2}t\bigr).\)

Below is a plot, as above, where the black line is the exact path, the blue line is the Euler integration approximation (initial values \(y_0 = 1\), \(v_0 = 0\)), and the red line is the Verlet integration approximation (initial values \(y_0 = 1\), \(y_{-1} = \cos(\frac{\pi}{2}\cdot (-0.1))\)), using \(\Delta t = 0.1\ \text{s}\), running it for just over one cycle (5 seconds).

As you can see, the Verlet approximation is very close to the true path, while the Euler approximation is quite far off and very unstable. My conclusion is that you should stay away from Euler integration if possible. Note also that there are many other integration methods (Runge-Kutta for example) that also have very good properties. However, Verlet integration has one very good property which is the ease of which you can apply positional constraints. Hopefully I’ll have the time to write up something about that soon.

Symplectic Euler

With Verlet integration, you keep track of two positions, instead of position and velocity. If you find this strange, and you really want position and velocity, there is a version of Euler integration that performs very similar to Verlet integration, called Symplectic Euler integration. Here, you first update the velocity, and then use this new velocity to update the position, i.e.

\(v_1 = v_0 + a_0\Delta t\),
\(p_1 = p_0 + v_1\Delta t\).

To see how similar this is to Verlet, suppose we have set up the initial velocity \(v_0\) based on the current and past positions \(p_0\) and \(p_{-1}\),

\(v_0 = (p_0 – p_{-1})/\Delta t\).

Then we see that

\(p_1 = p_0 + \Delta t v_1 = p_0 + (v_0 + \Delta t a_0) \Delta t = p_0 + v_0 \Delta t + a_0\Delta t^2\)
\( = p_0 + p_0 – p_{-1} + a_0\Delta t^2 = 2p_0 – p_{-1} + a_0\Delta t^2,\)

which is exactly the same expression as for updating the positions in Verlet integration. Note that there is a difference between Verlet and Symplectic Euler, and that is in the initial conditions.

Resources 

A Simple Time-Corrected Verlet Integration Method by Jonathan Dummer

So, I thought I’d try playing around with game physics. A long time ago, when I was about 14, I started programming. I really wanted to become a game programmer, and started learning C and C++, writing DOS programs on my 486. Writing directly to video memory (0A000h) in VGA mode 13h, optimizing the code with little assembler snippets. My highpoint was a demo for a friends BBS written completely in assembler, with a lovely lens effect and flames:

Soon though, my interest shifted more towards physics, and later mathematics, and I ended up with a PhD in mathematics earned by investigating wonderful properties of infinite dimensional representations of semi-simple Lie algebras. Now, many years down the line, my interest has move back towards programming, and I’ve decided to try and understand at least the basics of game physics modeling.

So, I started looking at simple integration methods. I was well acquainted with Euler and Runge-Kutta integration, but soon read about Verlet integration. Whereas in Euler integration you keep track of the position and velocity of a particle, in Verlet integration you keep track of the current and previous position of a particle. In that respect, they are equivalent, since you can get a velocity from the current and previous position, and vice versa. Where they differ are in the dynamics, how the new data are computed based on the current acceleration.

In Euler integration, you do

\(p_1 = p_0 + \Delta t v_0\),
\(v_1 = v_0 + \Delta t a_0\),

where \(v_0\) and \(p_0\) are the previous position and velocity, \(v_1\) and \(p_1\) are the new position and velocity, \(\Delta t\) is the time difference between frames, and \(a_0\) is the acceleration. Basically, at each time point you extrapolate what the position and velocity will be in the future based on the current velocity, position and acceleration.

In Verlet integration, you do the following instead:

\(p_1 = 2p_0 – p_{-1} + \Delta t^2a_0\),

where \(p_0\) is the current position, \(p_{-1}\) is the previous position, and \(p_1\) is the next position. If we solve the above equation for \(a_0\), we get

\(\displaystyle a_0 = \frac{p_1 – 2p_0 + x_{-1}}{\Delta t^2} = \frac{(p_1 – p_0) – (p_0 – p_{-1})}{\Delta t^2} = \frac{\frac{p_1 – p_0}{\Delta t} – \frac{p_0 – p_{-1}}{\Delta t}}{\Delta t}\).

Getting creative with our notation, we can write this as

\(\displaystyle a_0 = \frac{v_{\scriptstyle 1/2} – v_{-1/2}}{\Delta t}\),

i.e. we effectively use the velocity between our timepoints to compute the new positon \(p_1\). For various reasons, this gives a much better approximation than Euler integration. Actually, if the acceleration is constant, Verlet integration gives exact results.

So, to try this out I wrote a little test using OpenGL of a ball on a string affected by gravity, using Verlet integration. The string is modelled as a sequence of stiff springs, so in this example I also have constraints on the position of the endpoints of the springs, but I’ll leave that for the next point.