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

Question

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} + 2*\mathbf{k_2} + 2*\mathbf{k_3} + \mathbf{k_4}) $$

And the code I came up with:

rk4 :: Fractional a => [([a] -> a)] -> ([a] -> Bool) -> a -> [a] -> [a]
rk4 fs target h xs 
    | target xs = xs                                                            
    | otherwise = let k1 = pam fs xs      
                      k2 = pam fs ((+.) xs (map (h / 2 *) k1))           
                      k3 = pam fs ((+.) xs (map (h / 2 *) k2))       
                      k4 = pam fs ((+.) xs (map (h *) k3))       
                      dx = map (h / 6 *) (foldr (+.) (repeat 0) [k1, k2, k2, k3, k3, k4])
                      xn = (+.) xs dx                           
                  in rk4 fs target h xn                         
                  where (+.) = zipWith (+)

---

A very simple set of functions to test this:

target :: (Fractional a, Ord a) => [a] -> Bool                                  
target [x,y] = x >= 1.4                                                         

f1 :: Fractional a => [a] -> a                                                  
f1 [x,y] = 1                                                                    

f2 :: Fractional a => [a] -> a                                                  
f2 [x,y] = x + 3 * y   

Usage example:

euler [f1,f2] target 0.0001 [0,0]

---

I'm particularly unhappy with the need to define +. for element-wise addition (surely that exists, but Google didn't help me), and the way in which I added the various k values to create dx.

I also feel like I'm declaring too many variables (if that's what they're called?) in the let statement.

Oh, and my brackets are breeding like rabbits.

Answer 1

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 ]