haskelllambdanumerical-methodsverlet-integration

Verlet Integration in Haskell


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)2x 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 xn+1 = xn + dt*vn+1/2*a(xn,vn,n), while the Verlet method described in the Wikipedia goes xn+1 = 2*xn - xn-1+a(xn,vn,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.


Solution

  • The faithful Verlet sequence has xn depending on the previous two values of x -- xn-1 and xn-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.