Last time I implemented a Runge-Kutta method was in December 2001. Back then, I implemented the classical 4th order Runge-Kutta (RK4) method in C. However, in the past decade, whenever I needed a numerical solution of an ordinary differential equation (ODE), I used MATLAB. Now, I would like to use Haskell instead.
Consider the following initial value problem (IVP)
where the vector field and the initial condition are given. Solving an instance of the IVP consists of finding a function that satisfies both the differential equation and the initial condition.
Let be the time-step, and let be the -th time instant, where . Suppose that we solve a given instance of the IVP and obtain a closed-form solution ; we then discretize the closed-form solution to obtain an infinite sequence where . However, not all instances of the IVP have a closed-form solution (or even an analytic solution), and even in the cases where a closed-form solution does exist, it can be quite hard to find such a solution. Therefore, instead of discretizing the closed-form solution, one is tempted to discretize the ODE itself to obtain an approximation of the sequence , which we will denote by . We will call the approximate sequence a numerical solution.
Consider now the following discretized initial value problem (DIVP)
where depends on what particular numerical method one wants to use. We will use the (classical) 4th order Runge-Kutta (RK4) method , which is given by
where the variables are given by
In a nutshell, one hopes that by making the time-step “small enough”, the numerical solution will be “close enough” to , the discretization of the actual solution . This hoping involves a fair amount of faith. Especially so if floating-point arithmetic is used, as underflow, overflow, and round-off error are not to be dismissed.
Haskell implementation of RK4
The following Haskell script implements the RK4 method:
-- define 4th order Runge-Kutta map (RK4) rk4 :: Floating a => (a -> a) -> a -> a -> a rk4 f h x = x + (1/6) * (k1 + 2*k2 + 2*k3 + k4) where k1 = h * f (x) k2 = h * f (x + 0.5*k1) k3 = h * f (x + 0.5*k2) k4 = h * f (x + k3)
Let us now test this rather succinct implementation.
Consider the following IVP where , , and
The closed-form solution of this IVP is . We load the script above into GHCi and compute the numerical solution of the IVP using the higher-order function iterate:
*Main> -- discretize closed-form solution *Main> let xs = [ exp (-t) | t <- [0,0.1..]] *Main> -- define function f *Main> let f x = -x *Main> -- define time-step h *Main> let h = 1 / 100000 *Main> -- define initial condition *Main> let x0 = 1.0 *Main> -- compute numerical solution *Main> let xs_tilde = iterate (rk4 f h) x0 *Main> -- decimate numerical solution *Main> let xs_dec = [ xs_tilde !! k | k <- [0,1..], rem k 10000 == 0]
Note that the time-step is . Hence, the numerical solution will have samples in between and , which probably is too fine a discretization. Thus, we decimate the numerical solution by a factor of , so that we are left with only eleven samples in between and , with a uniform spacing of between consecutive samples. We can now compare the discretization of the closed-form solution with the heavily decimated numerical solution, as follows:
*Main> -- print discretized closed-form solution *Main> take 11 xs [1.0,0.9048374180359595,0.8187307530779818, 0.7408182206817179,0.6703200460356392, 0.6065306597126333,0.5488116360940264, 0.49658530379140947,0.44932896411722156, 0.4065696597405991,0.36787944117144233] *Main> -- print decimated numerical solution *Main> take 11 xs_dec [1.0,0.9048374180359563,0.8187307530779842, 0.7408182206817214,0.6703200460356442, 0.6065306597126376,0.5488116360940266, 0.49658530379141175,0.44932896411722334, 0.40656965974059905,0.367879441171439] *Main> -- compute difference of two solutions *Main> let es = zipWith (-) xs_dec xs *Main> -- print difference *Main> take 11 es [0.0,-3.219646771412954e-15,2.3314683517128287e-15, 3.552713678800501e-15,4.9960036108132044e-15, 4.3298697960381105e-15,2.220446049250313e-16, 2.275957200481571e-15,1.7763568394002505e-15, -5.551115123125783e-17,-3.3306690738754696e-15] *Main> -- compute energy of the difference *Main> sum $ take 11 $ map (^2) es 9.161263559461204e-29
Note how small the values of difference sequence are, of the order of . The implementation of the RK4 method seems to be working. Who needs MATLAB now? ;-)
 Richard W. Hamming, Numerical Methods for Scientists and Engineers, Dover Publications, 1973.