diff --git a/contents/verlet_integration/code/haskell/verlet.hs b/contents/verlet_integration/code/haskell/verlet.hs index 78a1580c7..df82004e6 100644 --- a/contents/verlet_integration/code/haskell/verlet.hs +++ b/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 diff --git a/contents/verlet_integration/verlet_integration.md b/contents/verlet_integration/verlet_integration.md index 881a87eaf..734d3cf05 100644 --- a/contents/verlet_integration/verlet_integration.md +++ b/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)