Higher Order ODE Methods: The Verlet Scheme

In the last post, we learned how to use the Euler-Cromer method to simulate a nonlinear system like the simple pendulum. Recall that in deriving the Euler-Cromer method, we took the Taylor expansion of $\theta(t)$ and discarded terms of order $t^2$ or higher. Perhaps if we kept these higher order terms, we might be able to derive a method that would give us better accuracy. Recall that we want to expand $\theta(t)$ about $t=t_0$ , since $\theta(t)$ is our unknown and we are trying to find an approximate value for it using our known value $\theta(t_0)$ . Doing this, we get:

$\theta(t) = \theta(t_0) + \dot{\theta}(t_0)(t-t_0) + \frac{1}{2}\ddot{\theta}(t_0)(t-t_0)^2 + \frac{1}{6}\dddot{\theta}(t_0)(t-t_0)^3 + O(t^4)$

This is the same result as before except this time I have explicitly included some higher order terms. Doing that was easy enough, but now it’s a bit unclear how to proceed. Recall that before we were able to solve for $\dot{\theta}(t_0)$ to obtain an approximation for the first derivative, which we then substituted into our system of ODEs. However, we can’t do that here, since we have too many unknowns – each extra time derivative of $\theta(t_0)$ is a new unknown, and we add one with every extra higher-order term! We can resolve this with an interesting trick. First, let’s do a change of variables to more explicitly address that we are moving “forward” in time. Let’s introduce the change of variables:

$t\equiv t+\tau \\[10pt] t_0\equiv t$

Hence $t$ is now our “current” time, the time at which we know the displacement, and $t+\tau$ is our “future” time, the time at which we want to find the displacement. With this change of variables, our Taylor expansion becomes:

$\theta(t+\tau) = \theta(t) + \dot{\theta}(t)\tau + \frac{1}{2}\ddot{\theta}(t)\tau^2 + \frac{1}{6}\dddot{\theta}(t)\tau^3 + O(\tau^4)$

Ignore that our system is deterministic for a second and imagine that we want to find the state of our system at a “past” time, given information about the “current” time. For this inverse problem, we just do the change of variables:

$t\equiv t-\tau \\[10pt] t_0\equiv t$

With this change of variables, we get the Taylor expansion:

$\theta(t-\tau) = \theta(t) + \dot{\theta}(t)(-\tau) + \frac{1}{2}\ddot{\theta}(t)(-\tau)^2 + \frac{1}{6}\dddot{\theta}(t)(-\tau)^3 + O(\tau^4)$

which becomes:

$\theta(t-\tau) = \theta(t) - \dot{\theta}(t)\tau + \frac{1}{2}\ddot{\theta}(t)\tau^2 - \frac{1}{6}\dddot{\theta}(t)\tau^3 + O(\tau^4)$

Adding these two expressions together, we see that the odd terms cancel out, giving us:

$\theta(t+\tau)+\theta(t-\tau) = 2\theta(t) + \ddot{\theta}(t)\tau^2 + O(\tau^4)$

Since our goal is to step forward in time, we move all the “current” and “past” values to the left hand side of the equation to obtain our final update equation:

$\theta(t+\tau)= 2\theta(t) -\theta(t-\tau) + \ddot{\theta}(t)\tau^2 + O(\tau^4)$

Note that by incorporating the “past” value of the displacement, we have improved our first-order method to a third-order method. Recall that the term $\dot{\theta}(t)$ is the acceleration and is a function of $\theta(t)$ , so it is known anytime the “current” displacement is known. The method we have derived is called the Verlet method, and it is a handy numerical integrator for more complicated, nonlinear systems. Note that the velocity, $\dot{\theta}(t)$ is completely ignored in our update equation, which is fine for us since we really only care about the displacement in our pendulum system. It will be useful to record the velocity anyway for energy calculations, so we revert to the typical Euler method to calculate it:

$\dot{\theta}(t+\tau) = \frac{\theta(t+\tau)-\theta(t-\tau)}{2\tau} + O(\tau^2)$

Note that we take two Euler steps instead of one, since the past position $\theta(t-\tau)$ is known anyway. Note that this step also must come after the position calculation, since it is not until then that $\theta(t+\tau)$ is known. Translating these update equations into code, we have:

theta_new = 2*theta_curr - theta_old + accel_curr*tau^2
omega_new = (theta_new - theta_old)/2*tau

Note that we have used $\verb|accel_curr|$ to represent $\ddot{\theta}(t)$ . To implement this update scheme, we will have to be careful with how we update the values in our loop. The core update loop, written in Python, looks like this:

for i in range(0,numSteps):
    # make the "current" the "past"
    theta_old = theta_curr
    # make the "new" the "current"
    theta_curr = theta_new
    omega_curr = omega_new
    # compute the acceleration using the "current"
    accel_curr = -(g/l)*sin(theta_curr)
    # step forward in time
    theta_new = 2*theta_curr - theta_old + accel_curr*tau**2
    omega_new = (theta_new - theta_old)/(2*tau)
    # record the values
    time_vec[i] = tau*i
    theta_vec[i] = theta_new
    omega_vec[i] = omega_new
    KE_vec[i] = (1/2)*m*l**2*omega_new**2
    PE_vec[i] = m*g*l*(1-cos(theta_new))

However, we run into a problem when starting our loop for the first time. Since we are only given the initial conditions at one time $t=0$ , how can we start the Verlet method when we don’t have knowledge of conditions at the time $t=0-\tau$ ? This is a fundamental limitation of the Verlet method, as in order to start the scheme we will have to find some way to obtain the “initial conditions” at two times. This will introduce some error, but in the long run it doesn’t really matter since we are only introducing this error during one step, the first step. One standard way to “jump-start” the scheme is to use basic kinematics to find $\theta_1$ , the displacement at $t=\tau$ , given the known initial conditions given at $t=0$ , denoted by $\theta_0$ , $\dot{\theta}_0$ , and $\ddot{\theta}_0$ :

$\theta_1 = \theta_0 + \dot{\theta}_0\tau + \frac{1}{2}\ddot{\theta}_0\tau^2$

We implement this in the code by doing the following:

theta_curr = theta0
omega_curr = omega0
theta_new = theta_curr + omega_curr*tau - (1/2)*(g/l)*sin(theta_curr)*tau**2
omega_new = (theta_new - theta_curr)/tau

Running the simple pendulum simulation with the Verlet method results in the outputs:

It doesn’t look that much different from our results using the Euler-Cromer method! How much benefit does extending our truncation to the third order really give us? One explicit way of visualizing this is noting that the Verlet method is exact for systems with a constant acceleration. After all, fully expanded, the update equation looks like:

$\theta(t+\tau)= 2\theta(t) -\theta(t-\tau) + \ddot{\theta}(t)\tau^2 + \theta^{(4)}(t)\tau^4 + \theta^{(6)}(t)\tau^6 + ...$

But if the acceleration is constant, then all time derivatives of $\theta(t)$ of order three and higher are just zero. Thus, when we truncate the $O(\tau^4)$ terms and higher, we are really truncating nothing! This can be shown by considering the case of uniform gravitational acceleration, where the acceleration is indeed constant and is simply given by $\vec{a}=\vec{g}$ . For simplicity, we will only consider the one-dimensional case. The Euler-Cromer update equations are: