Using 4th order Runge Kutta to solve the 2nd order differential equation of a damped oscillator

I think there may be some confusion over the external forcing term and the Runge Kutta derivative helper function F. The F in RK4 returns the derivative dX/dt of the system of first order differential equations X. The forcing term in a damped oscillator is unfortunately also called F but it is a function of t.

One of your issues is that the arity (number of parameters) of your RungeKutta() function and your call to that function do not match: you tried to do RungeKutta(F, y0, t, func), but the RungeKutta() function only takes arguments (f, y0, x) in that order.

In other words, the f parameter in your current RungeKutta() function should encapsulate the forcing function F(t).

You can do this with helpers:

# A constant function in your case, but this can be any function of `t`
def applied_force(t):
    # Note, you did not provide a value for `omega0`
    return omega0 ** 2

def rk_derivative_factory(c, k, F):
    return lambda X, t: np.array([X[1], -c * X[1] - k * X[0] + F(t)])

The rk_derivative_factory() is a function which takes a damping coefficient, a spring constant, and a forcing function F(t), and returns a function which takes a system X and a time step t as arguments (because this is what is demanded of you by the implementation of RungeKutta()).

Then,

omega0 = 0.234
m, c, k = 1, 0.25, 2
oscillator = HarmonicOscillator(m, c, k)

f = rk_derivative_factory(oscillator, applied_force)
x_osc = oscillator.solve(f, [1, 0])

Where solve() is defined like so:

def solve(self, func, y0):
    t = np.linspace(0,10 * self.period(), 1000)
    return RungeKutta(f, y0, t)

As an aside, I strongly recommend being more consistent about your variable names. You named the initial state of your oscillator x0, and were passing it to RungeKutta() as the argument for the parameter y0, and the x parameter of RungeKutta() represents time... Gets pretty confusing.

Full solution

Lastly, your implementation of RK4 wasn't quite correct, so I've fixed that and made some other slight improvements.

Note that one thing you might want to consider is making the HarmonicOscillator.solve() function take a solver. Then you could play around with different integrators.

import numpy as np


def RungeKutta(f, y0, x):
    y = np.zeros((len(x), len(y0)))
    y[0, :] = np.array(y0)
    h = x[1] - x[0]
    for i in range(0, len(x) - 1):
        # Many slight changes below
        k1 = np.array(f(y[i, :], x[i]))
        k2 = np.array(f(y[i, :] + h * k1 / 2, x[i] + h / 2))
        k3 = np.array(f(y[i, :] + h * k2 / 2, x[i] + h / 2))
        k4 = np.array(f(y[i, :] + h * k3, x[i] + h))
        y[i + 1, :] = y[i, :] + (h / 6) * (k1 + 2 * k2 + 2 * k3 + k4)
    return y


# A constant function in your case, but this can be any function of `t`
def applied_force(t):
    # Note, you did not provide a value for `omega0`
    return omega0 ** 2


def rk_derivative_factory(osc, F):
    return lambda X, t: np.array([X[1], (F(t) - osc.c * X[1] - osc.k * X[0]) / osc.m])


class HarmonicOscillator:
    def __init__(self, m, c, k):
        if (type(m) in (int, float)) and (m > 0):
            self.m = m 
        else:
            raise ValueError("Parameter 'm' must be a positive number")
        if (type(c) in (int, float)) and (c > 0):
            self.c = c
        else:
            raise ValueError("Parameter 'c' must be a positive number")
        if (type(k) in (int, float)) and (k > 0):
            self.k = k 
        else:
            raise ValueError("Parameter 'k' must be a positive number")
        
        self.T = 2 * np.pi * (self.m / self.k)**(0.5)

    def period(self):
        return self.T

    def solve(self, func, y0):
        t = np.linspace(0, 10 * self.period(), 1000)
        return RungeKutta(func, y0, t)

Demo:

import plotly.graph_objects as go

omega0 = 0.234
m, c, k = 1, 0.25, 2
oscillator = HarmonicOscillator(m, c, k)

f = rk_derivative_factory(oscillator, applied_force)
x_osc = oscillator.solve(f, [1, 0])

x, dx = x_osc.T
t = np.linspace(0, 10 * oscillator.period(), 1000)

fig = go.Figure(go.Scatter(x=t, y=x, name="x(t)"))
fig.add_trace(go.Scatter(x=t, y=dx, name="x'(t)"))

Output: