Click here to Skip to main content
15,867,999 members
Articles / Programming Languages / Python
Tip/Trick

Fourth Order Runge-Kutta Method in Python

Rate me:
Please Sign up or sign in to vote.
5.00/5 (11 votes)
29 Jul 2014CPOL2 min read 131.9K   797   12   8
Implementation of the fourth order Runge-Kutta method in Python for solving n-dimensional ordinary differential equations

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):

Python
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)$
Python
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):

Python
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], k1 is represented by [a1, b1, c1], k2 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:

Python
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:

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

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


Written By
United Kingdom United Kingdom
This member has not yet provided a Biography. Assume it's interesting and varied, and probably something to do with programming.

Comments and Discussions

 
QuestionMr tyebillion Pin
Member 1475786227-Feb-20 21:19
Member 1475786227-Feb-20 21:19 
Questionthanks Pin
Member 128379927-Nov-16 16:45
Member 128379927-Nov-16 16:45 
AnswerRe: thanks Pin
tyebillion26-Jan-17 3:45
tyebillion26-Jan-17 3:45 
Questionhelp for line 9 in method 1 Pin
Member 1269267018-Aug-16 20:18
Member 1269267018-Aug-16 20:18 
AnswerRe: help for line 9 in method 1 Pin
tyebillion19-Aug-16 10:44
tyebillion19-Aug-16 10:44 
GeneralMy vote of 5 Pin
Member 1131046914-Dec-14 14:58
Member 1131046914-Dec-14 14:58 
GeneralMy vote of 5 Pin
CatchExAs8-Jul-14 13:08
professionalCatchExAs8-Jul-14 13:08 
GeneralRe: My vote of 5 Pin
tyebillion9-Jul-14 1:45
tyebillion9-Jul-14 1:45 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.