One of the most widely used numerical algorithms for solving differential equation is the 4th order Runge-Kutta method. This post shows how the Runge-Kutta method can be written naturally as a fold over the set of points where the solution is needed. These points do not need to be evenly spaced.

Given a differential equation of the form

with initial condition *y*(*t*_{0}) = *y*_{0}, the 4th order Runge-Kutta method advances the solution by an amount of time *h* by

where

The Haskell code for implementing the accumulator function for `foldl`

looks very much like the mathematical description above. This is a nice feature of the `where`

syntax.

rk (t, y) t' = (t', y + h*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0) where h = t' - t k1 = f t y k2 = f (t + 0.5*h) (y + 0.5*h*k1) k3 = f (t + 0.5*h) (y + 0.5*h*k2) k4 = f (t + 1.0*h) (y + 1.0*h*k3)

Suppose we want to solve the differential equation *y* ‘ = (*t*^{2} – *y*^{2}) sin(*y*) with initial condition *y*(0) = -1, and we want to approximate the solution at [0.01, 0.03, 0.04, 0.06]. We would implement the right side of the equation as

f t y = (t**2 - y**2)*sin(y)

and fold the function `rk`

over our time steps with

foldl rk (0, -1) [0.01, 0.03, 0.04, 0.06]

This returns (0.06, -0.9527). The first part, 0.06, is no surprise since obviously we asked for the solution up to 0.06. The second part, -0.9527, is the part we’re more interested in.

If you want to see the solution at all the specified points and not just the last one, replace `foldl`

with `scanl`

.

scanl rk (0, -1) [0.01, 0.03, 0.04, 0.06]

This returns

[(0.0, -1.0), (0.01, -0.9917), (0.03, -0.9756), (0.042, -0.9678), (0.06, -0.9527)]

As pointed out in the previous post, writing algorithms as folds separates the core of the algorithm from data access. This would allow us, for example, to change `rk`

independently, such as using a different order Runge-Kutta method. (Which hardly anyone does. Fourth order is a kind of sweet spot.) Or we could swap out Runge-Kutta for a different ODE solver entirely just by passing a different function into the fold.

Next see how you could use a fold to compute mean, variance, skewness, and kurtosis in one pass through a stream of data using a fold.

Thanks! I’m sure you’re overflowing with future blog ideas, but…

Would you write about (cheap and/or easy) ways to get control over energy, or other conserved quantities?

If I understand correctly, Runge-Kutta usually adds friction or damping of some kind, but Euler on the other hand usually explodes? Sometimes that’s fine – e.g. if I’m already explicitly damping and turning a knob regarding how much, maybe the energy leak from using Runge-Kutta is hidden.

Do I need to change coordinates to make the conserved quantity explicit?

I skimmed a book one time on numerical DE methods that preserve energy. Maybe by Hairer? It’s not something I know about.

You can use symplectic integration schemes. Preserves conserverd quantities at some costs.