Euler vs Verlet

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