# Introduction

These are some very hasty notes on Runge-Kutta methods and IRK2 in particular. I make no apologies for missing lots of details. I may try and put these in a more digestible form but not today.

# Some Uncomprehensive Theory

In general, an implicit Runge-Kutta method is given by

where

and

Traditionally this is written as a Butcher tableau:

and even more succintly as

For a Gauß-Legendre method we set the values of the to the zeros of the shifted Legendre polynomials.

An explicit expression for the shifted Legendre polynomials is given by

The first few shifted Legendre polynomials are:

Setting

then the co-location method gives

For we have and thus and the Butcher tableau is

that is, the implicit RK2 method aka the mid-point method.

For we have and thus and and the Butcher tableau is

that is, the implicit RK4 method.

# Implicit RK2

Explicitly

where

Substituting in the values from the tableau, we have

where

and further inlining and substitution gives

which can be recognised as the mid-point method.

# Implementation

A common package for solving ODEs is gsl which Haskell interfaces via the hmatrix-gsl package. Some of gsl is coded with the conditional compilation flag *DEBUG* e.g. in msbdf.c but sadly not in the simpler methods maybe because they were made part of the package some years earlier. We can add our own of course; here’s a link for reference.

Let’s see how the Implicit Runge-Kutta Order 2 method does on the following system taken from the gsl documenation

which can be re-written as

but replacing *gsl_odeiv2_step_rk8pd* with *gsl_odeiv2_step_rk2imp*.

Here’s the first few steps

```
rk2imp_apply: t=0.00000e+00, h=1.00000e-06, y:1.00000e+00 0.00000e+00
-- evaluate jacobian
( 2, 2)[0,0]: 0
( 2, 2)[0,1]: 1
( 2, 2)[1,0]: -1
( 2, 2)[1,1]: -0
YZ:1.00000e+00 -5.00000e-07
-- error estimates: 0.00000e+00 8.35739e-20
rk2imp_apply: t=1.00000e-06, h=5.00000e-06, y:1.00000e+00 -1.00000e-06
-- evaluate jacobian
( 2, 2)[0,0]: 0
( 2, 2)[0,1]: 1
( 2, 2)[1,0]: -0.99997999999999998
( 2, 2)[1,1]: 1.00008890058234101e-11
YZ:1.00000e+00 -3.50000e-06
-- error estimates: 1.48030e-16 1.04162e-17
rk2imp_apply: t=6.00000e-06, h=2.50000e-05, y:1.00000e+00 -6.00000e-06
-- evaluate jacobian
( 2, 2)[0,0]: 0
( 2, 2)[0,1]: 1
( 2, 2)[1,0]: -0.999880000000002878
( 2, 2)[1,1]: 3.59998697518904009e-10
YZ:1.00000e+00 -1.85000e-05
-- error estimates: 1.48030e-16 1.30208e-15
rk2imp_apply: t=3.10000e-05, h=1.25000e-04, y:1.00000e+00 -3.10000e-05
-- evaluate jacobian
( 2, 2)[0,0]: 0
( 2, 2)[0,1]: 1
( 2, 2)[1,0]: -0.999380000000403723
( 2, 2)[1,1]: 9.6099972424212865e-09
YZ:1.00000e+00 -9.35000e-05
-- error estimates: 0.00000e+00 1.62760e-13
rk2imp_apply: t=1.56000e-04, h=6.25000e-04, y:1.00000e+00 -1.56000e-04
-- evaluate jacobian
( 2, 2)[0,0]: 0
( 2, 2)[0,1]: 1
( 2, 2)[1,0]: -0.996880000051409643
( 2, 2)[1,1]: 2.4335999548874554e-07
YZ:1.00000e+00 -4.68500e-04
-- error estimates: 1.55431e-14 2.03450e-11
```

Let’s see if we can reproduce this in a fast and loose way in Haskell

To make our code easier to write and read let’s lift some arithmetic operators to act on lists (we should really use the Naperian package).

```
> module RK2Imp where
```

```
> instance Num a => Num [a] where
> (+) = zipWith (+)
> (*) = zipWith (*)
```

The actual algorithm is almost trivial

```
> rk2Step :: Fractional a => a -> a -> [a] -> [a] -> [a]
> rk2Step h t y_n0 y_n1 = y_n0 + (repeat h) * f (t + 0.5 * h) ((repeat 0.5) * (y_n0 + y_n1))
```

The example Van der Pol oscillator with the same parameter and initial conditions as in the gsl example.

```
> f :: Fractional b => a -> [b] -> [b]
> f t [y0, y1] = [y1, -y0 - mu * y1 * (y0 * y0 - 1)]
> where
> mu = 10.0
```

```
> y_init :: [Double]
> y_init = [1.0, 0.0]
```

We can use co-recursion to find the fixed point of the Runge-Kutta equation arbitrarily choosing the 10th iteration!

```
> y_1s, y_2s, y_3s :: [[Double]]
> y_1s = y_init : map (rk2Step 1.0e-6 0.0 y_init) y_1s
```

```
> nC :: Int
> nC = 10
```

```
> y_1, y_2, y_3 :: [Double]
> y_1 = y_1s !! nC
```

This gives

```
ghci> y_1
[0.9999999999995,-9.999999999997499e-7]
```

which is not too far from the value given by gsl.

`y:1.00000e+00 -1.00000e-06 `

Getting the next time and step size from the gsl *DEBUG* information we can continue

```
> y_2s = y_1 : map (rk2Step 5.0e-6 1.0e-6 y_1) y_2s
>
> y_2 = y_2s !! nC
```

```
ghci> y_2
[0.999999999982,-5.999999999953503e-6]
```

gsl gives

`y:1.00000e+00 -6.00000e-06`

```
> y_3s = y_2 : map (rk2Step 2.5e-5 6.0e-6 y_2) y_3s
> y_3 = y_3s !! nC
```

One more time

```
ghci> y_3
[0.9999999995194999,-3.099999999372456e-5]
```

And gsl gives

`y:1.00000e+00 -3.10000e-05`

# Nix

The nix-shell which to run this is here and the C code can be built using something like

`gcc ode.c -lgsl -lgslcblas -o ode`

(in the nix shell).

The Haskell code can be used inside a ghci session.