# Runge-Kutta Implemented as a Fold in Haskell

### 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.

· Performance Zone · Tutorial
Save
2.74K Views

one of the most widely used numerical algorithms for solving differential equations 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.

Topics:
algorithm, performance, optimization

Published at DZone with permission of John Cook, DZone MVB.

Opinions expressed by DZone contributors are their own.