Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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.

like image 283
castle-bravo Avatar asked Mar 09 '14 05:03

castle-bravo


2 Answers

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.

like image 172
ErikR Avatar answered Sep 23 '22 08:09

ErikR


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

xn+1 - 2*xn + xn-1 = an*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 tn to

vn*(2*dt) = xn+1 - xn-1

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

-2*xn + 2*xn-1 = -2*vn*dt + an*dt^2

2*xn+1 - 2*xn = 2*vn*dt + an*dt^2

or

vn = (xn - xn-1)/dt + 0.5*an*dt

xn+1 = xn + vn*dt + 0.5*an*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 xn-1 by vn-1 and an-1 in the velocity computation. Then

vn = vn-1 + 0.5*(an-1 + an)*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 (tn-1,xn-1,vn-1,an-1). Then the next state is computed as

vn-0.5 = vn-1 + 0.5*an-1*dt

xn = xn-1 + vn-0.5*dt

Do collision detection with xn and vn-0.5

an = a(xn)

vn = vn-0.5 + 0.5*an*dt

Do statistics with xn and vn

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 a0 = a( x0 ).


Leapfrog Verlet

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

an = a(xn)

vn+0.5 = vn-0.5 + an*dt

xn+1 = xn + vn+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 (tn,xn,vn+0.5) as in

xn = xn-1 + vn-0.5*dt

an = a(xn)

vn+0.5 = vn-0.5 + an*dt

The initialization only needs to compute

v0+0.5 = v0 + 0.5*a( x0 )*dt

(end update)

Any statistics that one computes using xn and vn-0.5 or vn+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.

like image 28
Lutz Lehmann Avatar answered Sep 22 '22 08:09

Lutz Lehmann