{-# LANGUAGE RankNTypes, TypeFamilies #-}
import Data.VectorSpace
import Data.Manifold.Types
type ODESolver
= forall v . (VectorSpace v, Scalar v ~ ℝ)
=> (v -> v) -- ^ Function to integrate
-> ℝ -- ^ Time step
-> v -- ^ Start state
-> [v] -- ^ Sequence of calculated steps
euler :: ODESolver
euler f h y = y : euler f h (y ^+^ h*^f y)
rk₄ :: ODESolver
rk₄ f h y = y : rk₄ f h (y ^+^ (h/6)*^(k₁ ^+^ 2*^k₂ ^+^ 2*^k₃ ^+^ k₄))
where k₁ = f y
k₂ = f $ y ^+^ (h/2)*^k₁
k₃ = f $ y ^+^ (h/2)*^k₂
k₄ = f $ y ^+^ h*^k₃
import Graphics.Dynamic.Plot.R2
import Data.List
plotWindow
[ legendName name
$ plotLatest
[ lineSegPlot [(x,y) | (x,y,_) <- segment]
| segment <- take 400 <$> iterate (drop 10) traject ]
| (solver, name) <- [(euler, "Euler"), (rk₄,"RK₄")]
, let traject = solver (\(x,y,z) -> (σ*(y-x), x*(ρ-z), x*y - β*z))
0.002
(10,11,12)
where ρ=28; σ=10; β=8/3
]
GraphWindowSpecR2{lBound=-13.874009118754717, rBound=25.510352576781575, bBound=-19.831199555474484, tBound=35.94014689231086, xResolution=640, yResolution=480}