module Numeric.DDE (
integ
, integ'
, integRk4_2D
, integHeun2_2D
, Input (..)
, InputSnapshot (..)
, State (..)
, HistorySnapshot (..)
, Stepper1 (..)
, rk4
, heun2
) where
import qualified Data.Vector.Storable as V
import qualified Data.Vector.Storable.Mutable as VM
import System.IO.Unsafe ( unsafePerformIO )
import Numeric.DDE.Types
infixl 6 .+.
(.+.) :: V.Vector Double -> V.Vector Double -> V.Vector Double
(.+.) = V.zipWith (+)
infixl 7 *.
(*.) :: Double -> V.Vector Double -> V.Vector Double
(*.) c = V.map (* c)
rk4 :: Stepper1
rk4 = Stepper1 _rk4
where
_rk4 hStep rhs' (State !xy) !(!x_tau1, !x_tau1') !(!u1, !u1') = State xy_next
where
xy_next = xy .+. over6 *. (a .+. 2 *. b .+. 2 *. c .+. d)
over6 = recip 6
over2 = recip 2
! (State a') = rhs' (State xy, toHist1 x_tau1, inp1)
! a = hStep *. a'
! (State b') = rhs' (State $ xy .+. over2 *. a, toHist1 x_tau1_b, inp1_b)
! b = hStep *. b'
! (State c') = rhs' (State $ xy .+. over2 *. b, toHist1 x_tau1_c, inp1_c)
! c = hStep *. c'
! (State d') = rhs' (State $ xy .+. c, toHist1 x_tau1', inp1')
! d = hStep *. d'
! x_tau1_b = (x_tau1 + x_tau1') / 2
! x_tau1_c = x_tau1_b
! inp1 = Inp u1
! inp1_b = Inp $ (u1 + u1') / 2
! inp1_c = inp1_b
! inp1' = Inp u1'
toHist1 :: Double -> HistorySnapshot
toHist1 = Hist. V.singleton
heun2 :: Stepper1
heun2 = Stepper1 _heun2
where
_heun2 hStep rhs' (State !xy) !(!x_tau1, !x_tau1') !(!u1, !u1') = State xy_next
where
! (State f1) = rhs' (State xy, toHist1 x_tau1, Inp u1)
! xy' = xy .+. hStep *. f1
! (State f2) = rhs' (State xy', toHist1 x_tau1', Inp u1')
! xy_next = xy .+. (hStep / 2.0) *. (f1 .+. f2)
integ'
:: (State -> (Double, Double) -> (Double, Double) -> State)
-> Int
-> Int
-> Int
-> (State, V.Vector Double, Input)
-> (State, V.Vector Double)
integ' iter1 len1 krecord total (!xy0, !hist0, !(Input in1)) = a
where
a = unsafePerformIO $ do
! v <- VM.new (len1 + total)
copyHist v hist0
xy' <- go v len1 xy0
trace <- V.unsafeFreeze v
return (xy', V.slice (len1 + total krecord) krecord trace)
copyHist !v !hist = do
mapM_ (\i -> VM.unsafeWrite v i (hist V.! i)) [0..(V.length hist) 1]
go !v !i !xy
| i == len1 + total = do
return xy
| otherwise = do
x_tau1 <- VM.unsafeRead v (i len1)
x_tau1' <- VM.unsafeRead v (i len1 + 1)
let u1 = in1 V.! (i len1)
u1' = in1 V.! (i len1 + 1)
!xy' = iter1 xy (x_tau1, x_tau1') (u1, u1')
!(State xy'1) = xy'
!x' = xy'1 V.! 0
VM.unsafeWrite v i x'
go v (i + 1) xy'
integ
:: Stepper1
-> State
-> V.Vector Double
-> Int
-> Double
-> RHS
-> Input
-> (State, V.Vector Double)
integ (Stepper1 stp) state0 hist0 len1 hStep rhs' inp@(Input in1) = r
where
! totalIters = V.length in1 1
! iterator = stp hStep rhs'
! r = integ' iterator len1 totalIters totalIters (state0, hist0, inp)
integRk4_2D :: Int
-> Double
-> RHS
-> Input
-> (State, V.Vector Double)
integRk4_2D len1 = integ rk4 state0 hist0 len1
where
! state0 = State (V.replicate 2 0.0)
! hist0 = V.replicate len1 0
integHeun2_2D :: Int
-> Double
-> RHS
-> Input
-> (State, V.Vector Double)
integHeun2_2D len1 = integ heun2 state0 hist0 len1
where
! state0 = State (V.replicate 2 0.0)
! hist0 = V.replicate len1 0