ODE45 solver implementation in Python I have a project where I need ODE solver without dependencies to libraries like Scipy. I decide to implement ODE45. According to tutorials from internet and from what I remember from classes I implement it somehow.

The code below contains also example function (pendulum) with exactly same values as in the scipy ODE tutorial https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.integrate.odeint.html. The resulting plot from my code looks the same.

However, I am not really sure if I done it in the good way. Especially I would really appreciate any hints how to improve the functions ode45_step and ode45 and verification whether my implementation is correct. I do not need to check the input variables, it will be done somewhere else.

import numpy as np
import matplotlib.pylab as plt

def ode45_step(f, x, t, dt, *args):
    """
    One step of 4th Order Runge-Kutta method
    """
    k = dt
    k1 = k * f(t, x, *args)
    k2 = k * f(t + 0.5*k, x + 0.5*k1, *args)
    k3 = k * f(t + 0.5*k, x + 0.5*k2, *args)
    k4 = k * f(t + dt, x + k3, *args)
    return x + 1/6. * (k1 + k2 + k3 + k4)

def ode45(f, t, x0, *args):
    """
    4th Order Runge-Kutta method
    """
    n = len(t)
    x = np.zeros((n, len(x0)))
    x[0] = x0
    for i in range(n-1):
        dt = t[i+1] - t[i] 
        x[i+1] = ode45_step(f, x[i], t[i], dt, *args)
    return x

def f(t, y, b, c):
    """
    Pendulum example function.
    """
    theta = y[0]
    omega = y[1]
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return np.array(dydt)


b = 0.25
c = 5.0    

N = 101

x0 = np.array([np.pi - 0.1, 0.0])
t = np.linspace(0, 10, N)

x = ode45(f, t, x0, b, c)

plt.plot(x)
plt.show()
The following line is incorrect:

return x + 1/6. * (k1 + k2 + k3 + k4)


Should be:

return x + 1/6. * (k1 + 2*k2 + 2*k3 + k4)


Cross reference:
This code