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 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 = 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 omega = y 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()