2
\$\begingroup\$

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()
\$\endgroup\$
2
\$\begingroup\$

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

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.