{-# LANGUAGE RankNTypes, TypeFamilies, FlexibleContexts #-}
import Data.VectorSpace
import Data.AffineSpace
import Data.Manifold.Types
type ODESolver
= forall x . ( AffineSpace x, VectorSpace (Diff x)
, Scalar (Diff x) ~ ℝ )
=> (x -> Diff x) -- ^ Function to integrate
-> ℝ -- ^ Time step
-> x -- ^ Start state
-> [x] -- ^ 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=-18.277723235168448, rBound=22.059872841157212, bBound=-22.031233945119673, tBound=27.82502109487232, xResolution=640, yResolution=480}