ODE45 solver implementation in Python

Question

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

Answer 1

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