snippetMinor
How can I make this Euler/RK4 implementation more elegant?
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:
And then we define the euler method itself:
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}
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 xAnd 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
The first hit is
Bingo. So we know that
And we fill in the blanks
A similar transformation can be applied to
I'm not sure if I'd actually suggest this, but you could also introduce another function,
Finally, I think 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 seeuntil :: (a -> Bool) -> (a -> a) -> a -> auntil p f yields the result of applying f until p holds.Bingo. So we know that
euler is going to look likeeuler :: 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 comprehensionpam :: [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.