Verlet Integration in Haskell

Question

I am new to Haskell, and as an exercise I've been trying to implement some of the code (written in Mathematica) from Joel Franklin's book Computational Methods in Physics. I've written the following code to take a lambda expression (acceleration) as the first argument. In general, acceleration is of the form x'' = f(x',x,t), but not always of all three variables.

-- Implementation 2.1: The Verlet Method
verlet _  _  _  _  0 = [] 
verlet ac x0 v0 dt n = take n $ zip [0,dt..] $ verlet' ac x0 v0 0 dt
  where verlet' ac x0 v0 t0 dt = 
          let
            xt = x0 + dt*v0 + 0.5*(dt^2)*(ac x0 v0 t0)
            vt = v0 + dt*(ac x0 v0 t0)
          in
            xt:(verlet' ac xt vt (t0+dt) dt)

In ghci, I would run this code with the following command (the acceleration function a = -(2pi)²x comes from an exercise in the book):

verlet (\x v t -> -(2*pi)^2*x) 1 0 0.01 1000

My problem is that this isn't really the Verlet method - here x_n+1 = x_n + dt*v_n+1/2*a(x_n,v_n,n), while the Verlet method described in the Wikipedia goes x_n+1 = 2*x_n - x_n-1+a(x_n,v_n,n). How would I re-write this function to more faithfully represent the Verlet integration method?

Tangentially, is there a way to write this more elegantly and succinctly? Are there linear algebra libraries that would make this easier? I appreciate your advice.

Answer 1

The faithful Verlet sequence has x_n depending on the previous two values of x -- x_n-1 and x_n-2. A canonical example of such a sequence is the Fibonacci sequence which has this one-liner Haskell definition:

fibs :: [Int]
fibs = 0 : 1 : zipWith (+) fibs     (tail fibs)
                        -- f_(n-1)  f_n

This defines the Fibonacci sequence as an infinite (lazy) list. The self-reference to tail fibs gives you the previous term, and the reference to fibs gives you the term before that. The terms are then combined with (+) to yield the next term in the sequence.

You can take the same approach with your problem as follows:

type State = (Double, Double, Double)  -- (x, v, t) -- whatever you need here

step :: State -> State -> State
step s0 s1 = -- determine the next state based on the previous two states

verlet :: State -> State -> [State]
verlet s0 s1 = ss
  where ss = s0 : s1 : zipWith step ss (tail ss)

The data structure State holds whatever state variables you need - x, v, t, n, ... The function step is analogous to (+) in the Fibonacci case and computes the next state given the previous two. The verlet function determines the entire sequence of states given the initial two states.

Answer 2

Actually, if you read on you will find that both variants are represented on the wikipedia page.

---

Basic Verlet

The basic second order central difference quotient discretization for second order ODE x''(t)=a(x(t)) is

x_n+1 - 2*x_n + x_n-1 = a_n*dt^2

Note that there are not velocities in the iteration and also not in the acceleration function a(x). This is because Verlet integration is only then superior to other integration methods when the dynamical system is conservative, that is, -m*a(x) is the gradient of some potential function, and potential functions are static objects, they depend only on the position, not on time and not on velocity. Many frictionless mechanical systems fall into this category.

---

Velocity Verlet

Now set, using first order central difference quotients, the velocity at time t_n to

v_n*(2*dt) = x_n+1 - x_n-1

and add and subtract this to and from the first equation to obtain

-2*x_n + 2*x_n-1 = -2*v_n*dt + a_n*dt^2

2*x_n+1 - 2*x_n = 2*v_n*dt + a_n*dt^2

or

v_n = (x_n - x_n-1)/dt + 0.5*a_n*dt

x_n+1 = x_n + v_n*dt + 0.5*a_n*dt^2

This is one variant to write the velocity-Verlet algorithm.

---

(updated to have all state variables correspond to the same time before and after the iteration step)

Using the equations of the previous step from n-1 to n, one can replace x_n-1 by v_n-1 and a_n-1 in the velocity computation. Then

v_n = v_n-1 + 0.5*(a_n-1 + a_n)*dt

To avoid having two instances of any of the vectors x,v,a, one can arrange the update process so that everything is in place. Assume that at entry of the iteration step, the stored data corresponds to (t_n-1,x_n-1,v_n-1,a_n-1). Then the next state is computed as

v_n-0.5 = v_n-1 + 0.5*a_n-1*dt

x_n = x_n-1 + v_n-0.5*dt

Do collision detection with x_n and v_n-0.5

a_n = a(x_n)

v_n = v_n-0.5 + 0.5*a_n*dt

Do statistics with x_n and v_n

or as code

v += a*0.5*dt;
x += v*dt;
do_collisions(x,v);
a = eval_a(x);
v += a*0.5*dt;
do_statistics(x,v); 

Changing the order of these operations will destroy the Verlet scheme and significantly alter the results, rotation of the oparations is possible, but one has to be careful about the interpretation of the state after the iteration step.

The only initialization needed is the computation of a₀ = a( x₀ ).

---

Leapfrog Verlet

From the formulas of velocity Verlet one can see that for the updates of the position one does not need the velocities v_n, but only the half point velocities v_n+0.5. Then

a_n = a(x_n)

v_n+0.5 = v_n-0.5 + a_n*dt

x_n+1 = x_n + v_n+0.5*dt

or in code

a = eval_a(x);
v += a*dt;
x += v*dt;

Again, the order of these operations is fundamentally important, changes will lead to strange results for conservative systems.

---

(Update) However, one may rotate the execution sequence to

x += v*dt;
a = eval_a(x);
v += a*dt;

This corresponds to the iteration of the triples (t_n,x_n,v_n+0.5) as in

x_n = x_n-1 + v_n-0.5*dt

a_n = a(x_n)

v_n+0.5 = v_n-0.5 + a_n*dt

The initialization only needs to compute

v_0+0.5 = v₀ + 0.5*a( x₀ )*dt

(end update)

Any statistics that one computes using x_n and v_n-0.5 or v_n+0.5 will be off by an error proportional to dt, since the time indices do not match. Collisions can be detected with the position vector alone, but in the deflection the velocities need to be sensibly updated too.

Answer 3

Here is my solution after implementing the suggestion from user5402:

-- 1-Dimensional Verlet Method
type State = (,,) Double Double Double -- x, x', t

first :: State -> Double
first (x, _, _) = x

second :: State -> Double
second (_, x, _) = x

third :: State -> Double
third (_, _, x) = x

verlet1 :: (State -> Double) -> State -> Double -> Int -> [State]
verlet1 f s0 dt n = take n ss
  where
    ss = s0 : s1 : zipWith step ss (tail ss)
      where 
        s1 = (first s0 + dt*(second s0) + 0.5*dt^2*(f s0), 
              second s0 + dt*(f s0), third s0 + dt)
        step :: State -> State -> State
        step s0 s1 = (2*(first s1) - first s0 + dt^2*(f s1), 
                      second s1 + dt*(f s1), third s1 + dt)

I ran it in ghci using the following command:

verlet1 (\x -> -(2*pi)^2*(first x)) (1, 0, 0) 0.01 100

That seems to yield exactly what I was expecting - apparently sinusoidal motion! I have yet to plot x (if anyone has any suggestions how to do that in Haskell, they're welcome). Also, if you see any obvious refactorings, feel free to point them out. Thanks!