# Discrete dynamical systems with Haskell

A discrete dynamical system is given by **x***ₖ*₊₁ = *A***x***ₖ* for *k* = 0, 1, 2, … with a matrix *A* and an initial state vector **x**₀. Often times it’s useful to observe the long-term behavior of a system by finding its steady state vector—that is, a vector **q** for which **q** = *A***q**. We can do this in Haskell by declaring a linear transformation **x** ↦ *A***x** and finding its fixed point given the initial state. In other words, we will repeatedly apply the linear transformation, starting with the initial state, until we converge to a single vector (conventionally, one would just solve the equation **q** = *A***q**). Conveniently enough, Haskell defines some data structures and functions for computing fixed points, and I’ll cheat a bit by using a linear algebra library, although implementing one from scratch might be left to a future post.

For simplicity, we’ll work with 3 × 3 matrices and 3-vectors. First, we must define a data structure for state vectors and systems.

```
{-# LANGUAGE RecordWildCards #-}
module Main where
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
}
```

Now for a brief interlude on fixed points. The module `Control.Monad.Fix`

defines a function `fix`

that infinitely applies a function with a starting parameter, like so: f(f(…f(x)…)).

```
fix :: (a -> a) -> a
fix f = let x = f x in x
```

We can “fix” recursive functions by redefining them with an auxiliary parameter `rec`

that is used for the recursive call instead of the function itself. For example, the factorial function is normally defined as the following.

`f n = if n == 0 then 1 else n * f $ n - 1`

Our “fixed” version is:

`f rec n = if n == 0 then 1 else n * rec $ n - 1`

This essentially converts a recursive function to a non-recursive one, because `fix`

will keep passing said function as `rec`

. As a result, `fix`

enables support for recursion in primitive languages like the lambda calculus. Furthermore, notice that we always need to have a base case/terminating condition, so that `fix`

doesn’t recur forever.

So how do we get from `fix`

to a system’s steady state vector? First, we need to define a terminating condition: we want our linear transformation `t`

to terminate when **q** = *A***q**. However, we can’t check for strict equality because (1) we’re dealing with floating-point numbers and (2) we’re checking for *convergence* with this method, and not actual equality. With that in mind, let’s define an operator `===`

that returns whether or not two vectors are close enough to each other to be considered equal.

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

Now we’re ready to define a function `steadyState`

that takes a system and returns its steady state vector. Thus, we will pass `fix`

a function `t`

that will return the result of computing *A***x** if it `===`

**x**, otherwise it will recur on *A***x**. Once we pass `fix`

the initial state vector as well, it will magically converge to the steady state of the system.

```
steadyState :: System -> State
steadyState System{..} = fix t x0 where
t rec xk = let xk1 = a !* xk in
if xk1 === xk then
xk1 -- We found it, so terminate
else
rec xk1 -- Keep going
```

Now that we’ve defined everything, let’s test it on a discrete-time Markov chain, with an initial probability vector and stochastic matrix.

```
dtmc = System {
x0 = V3 1 0 0,
a = V3 (V3 0.5 0.2 0.3) -- Row-major representation
(V3 0.3 0.8 0.3)
(V3 0.2 0.0 0.4)
}
```

Let’s load up all of this code in `ghci`

and run it.

```
$ ghci steadystate.hs
...
*Main> steadyState dtmc
V3 0.30000038146975283 0.5999988555908204 0.10000076293942693
```

That’s reasonably close to the vector (.3, .6, .1), which is the right answer according to my textbook. The source for this project is on github.