Runge Kutta Fourth Order (RK4) Method Python Program

This program implements Runge Kutta fourth order method for solving ordinary differential equation in Python programming language.

Output of this Python program is solution for dy/dx = x + y with initial condition y = 1 for x = 0 i.e. y(0) = 1 and we are trying to evaluate this differential equation at y = 1 using RK4 method ( Here y = 1 i.e. y(1) = ? is our calculation point)

Python Source Code: RK4 Method

In this Python program x0 & y0 represents initial condition. xn is calculation point on which value of yn corresponding to xn is to be calculated using Runge Kutta method. step represents number of finite step before reaching to xn.


# RK-4 method python program

# function to be solved
def f(x,y):
    return x+y

# or
# f = lambda x: x+y

# RK-4 method
def rk4(x0,y0,xn,n):
    
    # Calculating step size
    h = (xn-x0)/n
    
    print('\n--------SOLUTION--------')
    print('-------------------------')    
    print('x0\ty0\tyn')
    print('-------------------------')
    for i in range(n):
        k1 = h * (f(x0, y0))
        k2 = h * (f((x0+h/2), (y0+k1/2)))
        k3 = h * (f((x0+h/2), (y0+k2/2)))
        k4 = h * (f((x0+h), (y0+k3)))
        k = (k1+2*k2+2*k3+k4)/6
        yn = y0 + k
        print('%.4f\t%.4f\t%.4f'% (x0,y0,yn) )
        print('-------------------------')
        y0 = yn
        x0 = x0+h
    
    print('\nAt x=%.4f, y=%.4f' %(xn,yn))

# Inputs
print('Enter initial conditions:')
x0 = float(input('x0 = '))
y0 = float(input('y0 = '))

print('Enter calculation point: ')
xn = float(input('xn = '))

print('Enter number of steps:')
step = int(input('Number of steps = '))

# RK4 method call
rk4(x0,y0,xn,step)

Output

Output of above implementation to solve ordinary differential equation by RK4 is:

Enter initial conditions:
x0 = 0
y0 = 1
Enter calculation point: 
xn = 1
Enter number of steps:
Number of steps = 2

--------SOLUTION--------
-------------------------
x0	y0	yn
-------------------------
0.0000	1.0000	1.7969
-------------------------
0.5000	1.7969	3.4347
-------------------------

At x=1.0000, y=3.4347