## Introduction

The Python code presented here is for the fourth order Runge-Kutta method in `n`

-dimensions. The Runge-Kutta method is a mathematical algorithm used to solve systems of ordinary differential equations (ODEs). The general form of these equations is as follows:

$\Large\begin{aligned} \dot{x}&=f(t, x) \\ x(t_{0})&=x_{0} \end{aligned}$

Where x is either a scalar or vector. The fourth order Runge-Kutta method is given by:

$\Large\begin{aligned} x_{i+1}&=x_{i}+(k_{1}+2(k_{2}+k_{3})+k_{4})/6 \\ t_{i+1}&=t_{i}+h \end{aligned}$

where `h > 0`

is a step size parameter, `i=1, 2, 3,...`

and:

$\Large\begin{aligned} k_{1}&=f(t_{i}, x_{i})h \\ k_{2}&=f(t_{i}+\frac{h}{2}, x_{i}+\frac{k_{1}}{2})h \\ k_{3}&=f(t_{i}+\frac{h}{2}, x_{i}+\frac{k_{2}}{2})h \\ k_{4}&=f(t_{i}+h, x_{i}+k_{3})h \end{aligned}$

The Runge-Kutta method offers greater accuracy than the method of multiplying each function in the ODEs by a step size parameter and adding the results to the current values in `x`

.

## Implementation

It is common practise to eliminate `t`

with a suitable substitution such as:

$\Large c=\omega t$

Hence:

$\Large\dot{c}=\omega$

Presented here are two techniques for implementing the fourth order Runge-Kutta method. Here is a general Python function for the method in `n`

-dimensions implemented using arrays (technique 1):

def rKN(x, fx, n, hs):
k1 = []
k2 = []
k3 = []
k4 = []
xk = []
for i in range(n):
k1.append(fx[i](x)*hs)
for i in range(n):
xk.append(x[i] + k1[i]*0.5)
for i in range(n):
k2.append(fx[i](xk)*hs)
for i in range(n):
xk[i] = x[i] + k2[i]*0.5
for i in range(n):
k3.append(fx[i](xk)*hs)
for i in range(n):
xk[i] = x[i] + k3[i]
for i in range(n):
k4.append(fx[i](xk)*hs)
for i in range(n):
x[i] = x[i] + (k1[i] + 2*(k2[i] + k3[i]) + k4[i])/6
return x

Both `x`

and `fx`

are arrays, the latter is an array of functions, and `n`

is the number of dimensions. The function definitions correspond with the ODEs being solved. Here is the ODE and an example usage of the `rKN`

function for the forced Van der Pol oscillator:

$\Large\ddot{y}=\mu(1-y^2)\dot{y}-y+Asin(\omega t)$

Let:

$\Large x=(\dot{y},y,\omega t)$

def fa1(x):
return 0.9*(1 - x[1]*x[1])*x[0] - x[1] + math.sin(x[2])
def fb1(x):
return x[0]
def fc1(x):
return 0.5
def VDP1():
f = [fa1, fb1, fc1]
x = [1, 1, 0]
hs = 0.05
for i in range(20000):
x = rKN(x, f, 3, hs)

In the above functions `μ=0.9`

, `A=1`

and `ω=0.5`

. Here is a general Python function for the fourth order Runge-Kutta method in 3-dimensions (technique 2):

def rK3(a, b, c, fa, fb, fc, hs):
a1 = fa(a, b, c)*hs
b1 = fb(a, b, c)*hs
c1 = fc(a, b, c)*hs
ak = a + a1*0.5
bk = b + b1*0.5
ck = c + c1*0.5
a2 = fa(ak, bk, ck)*hs
b2 = fb(ak, bk, ck)*hs
c2 = fc(ak, bk, ck)*hs
ak = a + a2*0.5
bk = b + b2*0.5
ck = c + c2*0.5
a3 = fa(ak, bk, ck)*hs
b3 = fb(ak, bk, ck)*hs
c3 = fc(ak, bk, ck)*hs
ak = a + a3
bk = b + b3
ck = c + c3
a4 = fa(ak, bk, ck)*hs
b4 = fb(ak, bk, ck)*hs
c4 = fc(ak, bk, ck)*hs
a = a + (a1 + 2*(a2 + a3) + a4)/6
b = b + (b1 + 2*(b2 + b3) + b4)/6
c = c + (c1 + 2*(c2 + c3) + c4)/6
return a, b, c

This function performs the same calculations as `rKN`

, but specifically in 3-dimensions and with the loops unravelled. Numerics need to be passed to the parameters `a`

, `b`

and `c`

. Functions, each taking three numerics, and each returning numerics need to be passed to the parameters `fa`

, `fb`

and `fc`

. So `x`

is represented here by `[a, b, c]`

, k_{1} is represented by `[a1, b1, c1]`

, k_{2} by `[a2, b2, c2]`

and so on. The variables `ak`

, `bk`

and `ck`

are temporary variables used to optimize the calculations. Here is an example usage of the `rK3`

function, again for the forced Van der Pol oscillator:

def fa2(a, b, c):
return 0.9*(1 - b*b)*a - b + math.sin(c)
def fb2(a, b, c):
return a
def fc2(a, b, c):
return 0.5
def VDP2():
a, b, c, hs = 1, 1, 0, 0.05
for i in range(20000):
a, b, c = rK3(a, b, c, fa2, fb2, fc2, hs)

The program *TestRK.py* demonstrates that technique 2 is faster by a factor of about 3.

## Code Generation

The program *GenRK.py* can be used to generate the Python code (for technique 2) in any dimension up to `n=26`

. Here the function `rK3`

was created using the generator program *GenRK.py* with `n=3`

. The generator works using a number of `for`

loops in a function called `genRK`

. Below is the code for *GenRK.py*:

def getChr(i):
return chr(i + 97)
def genRK(n):
print("# fourth order Runge-Kutta method in " + str(n) + " dimensions")
u = ""
v = ""
f = ""
for i in range(n):
c = getChr(i)
if i != 0:
u += ", "
v += ", "
f += ", "
u += c
v += c + "k"
f += "f" + c
print("def rK" + str(n) + "(" + u + ", " + f + ", hs" + "):")
for i in range(n):
c = getChr(i)
print("\t" + c + "1 = f" + c + "(" + u + ")*hs")
for i in range(n):
c = getChr(i)
print("\t" + c + "k = " + c + " + " + c + "1*0.5")
for i in range(n):
c = getChr(i)
print("\t" + c + "2 = f" + c + "(" + v + ")*hs")
for i in range(n):
c = getChr(i)
print("\t" + c + "k = " + c + " + " + c + "2*0.5")
for i in range(n):
c = getChr(i)
print("\t" + c + "3 = f" + c + "(" + v + ")*hs")
for i in range(n):
c = getChr(i)
print("\t" + c + "k = " + c + " + " + c + "3")
for i in range(n):
c = getChr(i)
print("\t" + c + "4 = f" + c + "(" + v + ")*hs")
for i in range(n):
c = getChr(i)
print("\t" + c + " = " + c + " + (" + c + "1 + 2*(" + c + "2 + " + c + "3) + " + c + "4)/6")
print("\treturn " + u)

## Conclusion

Whilst the first technique is easier to implement, it is somewhat slower than the second. Technique 2 becomes harder to implement in higher dimensions, though it is relatively easy to generate the code.