In [45]:
{-# LANGUAGE TypeFamilies #-}
import Graphics.Rendering.Chart.Easy
import Data.Default.Class
import Control.Lens
import Control.Monad
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 [ ]: