HiveBrain v1.2.0
Get Started
← Back to all entries
snippetMinor

How can I make this Euler/RK4 implementation more elegant?

Submitted by: @import:stackexchange-codereview··
0
Viewed 0 times
thisrk4canmakemoreeulereleganthowimplementation

Problem

One of the things that I'm doing to teach myself is converting some numerical methods from existing Python code (they seem to me to lend themselves to functional programming quite well).

I'd like to know how I can improve the general style of my code.
In particular, I feel like I'm overusing brackets, that I should be able to use more inbuilt functions (that I'm not aware of), and that I may be misusing the let statement.

Background

A quick refresher on the terminology I'm using (slightly different to wiki*):

Given a state vector \$\mathbf{x}\$, and a function vector \$\mathbf{f}(\mathbf{x})\$ (which is our system of first order differential equations), we want to iterate using step size \$h\$ until a given target \$t(\mathbf{x})\$ is reached.

*Different in that I don't make the useless distinction between dependent and independent variables that some texts, and wiki, use.

Euler Method

For this method, we use a very simple approximation:

$$
\mathbf{k_1} = \mathbf{f}(\mathbf{x}) \\
\mathbf{x}_{new} = \mathbf{x} + h*\mathbf{k_1}
$$

In the code, I use what is basically a 'reverse map' (ie, apply a list of functions to a single argument) that I found online:

pam :: [a -> b] -> a -> [b]
pam f x = map g f
    where g h = h x


And then we define the euler method itself:

euler :: Fractional a => [([a] -> a)] -> ([a] -> Bool) -> a -> [a] -> [a]
euler fs target h xs
    | target xs = xs           
    | otherwise = let k1 = pam fs xs
                      xn = (+.) xs (map (h *) k1)       
                  in euler fs target h xn  
                  where (+.) = zipWith (+)


RK4 (4th order Runge-Kutta)

For this method, the approximation is a bit more complex:

$$
\mathbf{k_1} = \mathbf{f}(\mathbf{x}) \\
\mathbf{k_2} = \mathbf{f}(\mathbf{x} + \frac{h}{2}\mathbf{k_1}) \\
\mathbf{k_3} = \mathbf{f}(\mathbf{x} + \frac{h}{2}\mathbf{k_2}) \\
\mathbf{k_4} = \mathbf{f}(\mathbf{x} + h*\mathbf{k_3}) \\
\mathbf{x_{new}} = \mathbf{x} + \frac{h}{6}(\mathbf{k_1}

Solution

I feel ... that I should be able to use more inbuilt functions (that
I'm not aware of)

Let me introduce you to your new best friend... Hoogle. Seriously, it's the best thing since sliced... monad analogies. Install a local copy as well.

We're producing a sequence of values by repeated application of a function, so let's search Hoogle for

(a -> a) -> a -> [a]


The first hit is iterate :: (a -> a) -> a -> [a], which looks like a good candidate. But a little further down the page, we see


until :: (a -> Bool) -> (a -> a) -> a -> a


until p f yields the result of applying f until p holds.

Bingo. So we know that euler is going to look like

euler :: Fractional a => [[a] -> a] -> ([a] -> Bool) -> a -> [a] -> [a]
euler fs target h = until target step
    where step x = ???


And we fill in the blanks

euler :: Fractional a => [[a] -> a] -> ([a] -> Bool) -> a -> [a] -> [a]
euler fs target h = until target step
    where step x = zipWith (+) x (map (h *) (pam fs x))


A similar transformation can be applied to rk4:

rk4step fs h x = x +. dx
    where k1 = pam fs x
          k2 = pam fs $ x +. map (h / 2 *) k1
          k3 = pam fs $ x +. map (h / 2 *) k2
          k4 = pam fs $ x +. map (h *) k3
          dx = map (h / 6 *) $ k1 +. k2 +. k2 +. k3 +. k3 +. k4
          (+.) = zipWith (+)

rk4 :: Fractional a => [[a] -> a] -> ([a] -> Bool) -> a -> [a] -> [a]
rk4 fs target h = until target (rk4step fs h)


I'm not sure if I'd actually suggest this, but you could also introduce another function, *.:

rk4step fs h x = x +. dx
    where k1 = pam fs x
          k2 = pam fs $ x +. ((h / 2) *. k1)
          k3 = pam fs $ x +. ((h / 2) *. k2)
          k4 = pam fs $ x +. (h *. k3)
          dx = (h / 6) *. (k1 +. (2 *. k2) +. (2 *. k3) +. k4)
          (+.) = zipWith (+)
          (*.) n = map (n *)


Finally, I think that pam would be more readable as a list comprehension

pam :: [a -> b] -> a -> [b]
pam fs x = [ f x | f <- fs ]

Code Snippets

(a -> a) -> a -> [a]
euler :: Fractional a => [[a] -> a] -> ([a] -> Bool) -> a -> [a] -> [a]
euler fs target h = until target step
    where step x = ???
euler :: Fractional a => [[a] -> a] -> ([a] -> Bool) -> a -> [a] -> [a]
euler fs target h = until target step
    where step x = zipWith (+) x (map (h *) (pam fs x))
rk4step fs h x = x +. dx
    where k1 = pam fs x
          k2 = pam fs $ x +. map (h / 2 *) k1
          k3 = pam fs $ x +. map (h / 2 *) k2
          k4 = pam fs $ x +. map (h *) k3
          dx = map (h / 6 *) $ k1 +. k2 +. k2 +. k3 +. k3 +. k4
          (+.) = zipWith (+)

rk4 :: Fractional a => [[a] -> a] -> ([a] -> Bool) -> a -> [a] -> [a]
rk4 fs target h = until target (rk4step fs h)
rk4step fs h x = x +. dx
    where k1 = pam fs x
          k2 = pam fs $ x +. ((h / 2) *. k1)
          k3 = pam fs $ x +. ((h / 2) *. k2)
          k4 = pam fs $ x +. (h *. k3)
          dx = (h / 6) *. (k1 +. (2 *. k2) +. (2 *. k3) +. k4)
          (+.) = zipWith (+)
          (*.) n = map (n *)

Context

StackExchange Code Review Q#60585, answer score: 3

Revisions (0)

No revisions yet.