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.
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.