Steady state

steadystate.hs

{-# LANGUAGE RecordWildCards #-}
import Data.Ratio
import Linear.Matrix
import Linear.V3
import Linear.Epsilon
import Control.Monad.Fix

type State = V3 Double

data System = System {
    x0 :: State,
    a :: M33 Double
}

(===) :: State -> State -> Bool
a === b = nearZero $ a - b

steadyState :: System -> State
steadyState System{..} = fix t x0 where
    t rec xk = let xk1 = a !* xk in
        if xk === xk1 then
            xk
        else
            rec xk1

dtmc = System {
    x0 = V3 1 0 0,
    a = V3 (V3 0.5 0.2 0.3)
           (V3 0.3 0.8 0.3)
           (V3 0.2 0.0 0.4)
}