Added Stormer Verlet and Velocity Verlet (algorithm-archivists#398)

Author: jiegillet

contents/verlet_integration/code/haskell/verlet.hs

@@ -1,36 +1,54 @@
 -- submitted by Jie
-type Position     = [Double]
-type Speed        = [Double]
-type Time         = Double
-type Particle     = (Position, Speed, Acceleration, Time)
-type Acceleration = [Double]
-
-verletStep :: (Particle -> Acceleration)
-              -> Time
-              -> Particle
-              -> Particle
-              -> Particle
-verletStep acc dt (xOld, _, aOld, _) (x, v, a, t) = (x', v', a', t+dt)
+type Time = Double
+
+type Position = Double
+
+type Speed = Double
+
+type Acceleration = Double
+
+type Particle = (Position, Speed, Acceleration, Time)
+
+type Model = Particle -> Acceleration
+
+type Method = Model -> Time -> Particle -> Particle -> Particle
+
+verlet :: Method
+verlet acc dt (xOld, _, _, _) (x, _, a, t) = (x', v', a', t + dt)
   where
-  x' = zipWith3 (\xOld x a -> 2*x - xOld + a*dt^2 ) xOld x a
-  v' = zipWith3 (\v a aOld -> v + 0.5*(aOld + a)*dt) v a aOld
-  a' = acc (x', v', [], t+dt)
-
-trajectory :: (Particle -> Acceleration)
-              -> Time
-              -> Particle
-              -> [Particle]
-trajectory acc dt p0@(x, v, a, t0) = t
+    x' = 2 * x - xOld + a * dt ^ 2
+    v' = 0
+    a' = acc (x', v', a, t + dt)
+
+stormerVerlet :: Method
+stormerVerlet acc dt (xOld, _, _, _) (x, _, a, t) = (x', v', a', t + dt)
+  where
+    x' = 2 * x - xOld + a * dt ^ 2
+    v' = (x' - x) / dt
+    a' = acc (x', v', a, t + dt)
+
+velocityVerlet :: Method
+velocityVerlet acc dt (xOld, _, aOld, _) (x, v, a, t) = (x', v', a', t + dt)
   where
-  t  = p0 : p1 : zipWith (verletStep acc dt) t (tail t)
-  p1 = (x', v', acc (x', v', [], t0+dt), t0+dt)
-  x' = zipWith3 (\x v a -> x + v*dt + 0.5*a*dt^2 ) x v a
-  v' = zipWith (\v a -> v + a*dt) v a
+    x' = 2 * x - xOld + a * dt ^ 2
+    v' = v + 0.5 * (aOld + a) * dt
+    a' = acc (x', v', a, t + dt)
 
-freeFall :: Particle
-freeFall = last $ takeWhile (\([x],_,_,_) -> x > 0) $ trajectory acc dt p0
+trajectory :: Method -> Model -> Time -> Particle -> [Particle]
+trajectory method acc dt p0@(x, v, a, t0) = traj
   where
-  p0    = ([5], [0], [-10], 0)
-  dt    = 0.001
-  acc _ = [-10]
+    traj = p0 : p1 : zipWith (method acc dt) traj (tail traj)
+    p1 = (x', v', acc (x', v', a, t0 + dt), t0 + dt)
+    x' = x + v * dt + 0.5 * a * dt ^ 2
+    v' = v + a * dt
 
+main :: IO ()
+main = do
+  let p0 = (5, 0, -10, 0)
+      dt = 0.001
+      freefall _ = -10
+      aboveGround (x, _, _, _) = x > 0
+      integrate m = last $ takeWhile aboveGround $ trajectory m freefall dt p0
+  print $ integrate verlet
+  print $ integrate stormerVerlet
+  print $ integrate velocityVerlet

contents/verlet_integration/verlet_integration.md

@@ -41,8 +41,7 @@ Here is what it looks like in code:
 {% sample lang="py" %}
 [import:1-9, lang:"python"](code/python/verlet.py)
 {% sample lang="hs" %}
-Unfortunately, this has not yet been implemented in haskell, so here's Julia code:
-[import:1-13, lang:"julia"](code/julia/verlet.jl)
+[import:14-21, lang:"haskell"](code/haskell/verlet.hs)
 {% sample lang="scratch" %}
 Unfortunately, this has not yet been implemented in scratch, so here's Julia code:
 [import:1-13, lang:"julia"](code/julia/verlet.jl)
@@ -88,8 +87,7 @@ However, the error for this is $$\mathcal{O}(\Delta t)$$, which is quite poor, b
 {% sample lang="py" %}
 [import:11-21, lang:"python"](code/python/verlet.py)
 {% sample lang="hs" %}
-Unfortunately, this has not yet been implemented in scratch, so here's Julia code:
-[import:15-31, lang:"julia"](code/julia/verlet.jl)
+[import:23-28, lang:"haskell"](code/haskell/verlet.hs)
 {% sample lang="scratch" %}
 Unfortunately, this has not yet been implemented in scratch, so here's Julia code:
 [import:15-31, lang:"julia"](code/julia/verlet.jl)
@@ -149,8 +147,7 @@ Here is the velocity Verlet method in code:
 {% sample lang="py" %}
 [import:23-32, lang:"python"](code/python/verlet.py)
 {% sample lang="hs" %}
-Unfortunately, this has not yet been implemented in haskell, so here's Julia code:
-[import:33-45, lang:"julia"](code/julia/verlet.jl)
+[import:30-35, lang:"haskell"](code/haskell/verlet.hs)
 {% sample lang="scratch" %}
 Unfortunately, this has not yet been implemented in scratch, so here's Julia code:
 [import:33-45, lang:"julia"](code/julia/verlet.jl)