[Haskell-cafe] representing differencial equations in haskell

[Paul L]

Here is a minimal answer using Yampa-like Signal Function and Arrow
notation. You have to load this using "ghci -farrows".

> import Control.Arrow

The differential equation you gave are indeed:

  i = integral (cos i) + i0
  c = integral (alpha * (i - c)) + c0

where i0 and c0 are the initial constant values for i and c. Turn it
into a Haskell definition using Arrow Notation:

> cSF :: SF () Double
> cSF = proc () -> do
>      rec
>        i <- integral i0 -< cos i
>        c <- integral c0 -< alpha * (i - c)
>      returnA -< c

Notice that it matches the math definition almost line to line. And
this is indeed the solution. (unfoldSF cSF) will give you an infinite
list of the sequence sampled at interval dt. Below are the
stream-based definitions for the signal function SF using Arrow.

> dt = 0.01
> alpha = 1
> i0 = 0
> c0 = 0

The signal function is implemented as a stream function/transformer.

> newtype SF a b = SF { runSF :: [a] -> [b] }

It's made instances of Arrow and ArrowLoop.

> instance Arrow SF where
>    arr f = SF g
>       where g (x:xs) = f x : g xs
>    first (SF f) = SF g
>       where g l = let (x, y) = unzip l
>                    in zip (f x) y
>    (SF f) >>> (SF g) = SF (g . f)

> instance ArrowLoop SF where
>   loop (SF f) = SF g
>       where g x = let (y, z) = unzip' (f (zip x z))
>                    in y

Note that for ArrowLoop we must use a modified unzip that works on
streams (instead of eager pattern matching).

> unzip' ~((x, y):l) = let (xs, ys) = unzip' l
>                       in (x:xs, y:ys)

The integral function implements simple Euler's rule.

> integral :: Double -> SF Double Double
> integral i = SF (g i)
>       where g i (x:xs) = i : g (i + dt * x) xs

Unfolding the SF into a list by given it an input stream.

> unfoldSF :: SF () b -> [b]
> unfoldSF (SF f) =  f inp
>   where inp = ():inp

Hope it helps to serve an example of FRP/Yampa and Arrow. The full
Yampa implementation is continuation based rather than stream, because
it's difficult to have event switches using streams, but that's
another story.

Regards,
Paul Liu

On 9/25/07, Thomas Girod <girodt at gmail.com> wrote:
> Hi there.
>
> Let's say I have mathematical model composed of several differential
> equations, such as :
>
> di/dt = cos(i)
> dc/dt = alpha * (i(t) - c(t))
>
> (sorry my maths are really bad, but I hope you get the point)
>
> I would like to approximate the evolution of such a system iteratively. How
> would you do that in haskell ?
>
> cheers,
>
> Thomas
>
> _______________________________________________
> Haskell-Cafe mailing list
> Haskell-Cafe at haskell.org
> http://www.haskell.org/mailman/listinfo/haskell-cafe
>
>