In [45]:
{-# LANGUAGE TypeFamilies #-}
import Graphics.Rendering.Chart.Easy
import Data.Default.Class
import Control.Lens
import Control.Arrow
import Data.VectorSpace
In [13]:
type  = Double
type Time = Double

s :: Time ->
s t = fromIntegral (round t) - t
In [26]:
signalPlot :: String -> [(Time,)] -> EC (Layout  ) ()
signalPlot capt sig = do
layout_x_axis . laxis_title .= "t"
plot (line capt [takeWhile ((<3).fst) sig])
In [27]:
toRenderable $signalPlot "s(t)" [(t,s t) | t <- [0,0.01..]] In [42]: rk₄ :: (VectorSpace y, Scalar y ~ ) => (Time -> y -> y) -- the function r_t(t,r) -> Time -- time-step for the solver -> [(Time, y)] -- sequence of solution snapshots [r(t)] rk₄ f h = go 0 zeroV where go t y = (t,y) : go (t+h) (y ^+^ (h/6)*^(k₁ ^+^ 2*^k₂ ^+^ 2*^k₃ ^+^ k₄)) where k₁ = f t y k₂ = f (t + h/2) (y ^+^ (h/2)*^k₁) k₃ = f (t + h/2) (y ^+^ (h/2)*^k₂) k₄ = f (t + h) (y ^+^ h*^k₃) In [43]: toRenderable$ forM [0 .. 8] $\η -> signalPlot ("η = "++show η)$ rk₄ (\t r -> η * (s t - r)) 0.01
In [63]:
toRenderable $forM [0, 5 .. 30]$ \η -> signalPlot ("η = "++show η) $second fst <$> rk₄ (\t (r,u) -> (u, η^2 * (s t - r) - η*u)) 0.005
In [ ]: