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

Question

So im tasked with using the 4th order Runge Kutta Meathod to solve the 2nd order differential equation of a damped occilator.

my function for the runge-kutta meathod looks as such

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):
        k1=h*np.array(f(y[i,:],x[i]))
        k2=h*np.array(f(y[i,:]+k1/2,x[i]+h/2))
        k3=h*np.array(f(y[i,:]+k2/2,x[i]+h/2))
        k4=h*np.array(f(y[i,:]+k3,x[i]+h))
        y[i+1,:]=y[i,:]+k1/6+k2/3+k3/3+k4/6
    return y

the rungeKutta function works fine, and I have tested it with a list of example inputs so that doesnt seem to be the problem

im given question parameters and have to make a class to solve the problem

class harmonicOscillator:
    def __init__(self,m,c,k):
        if((m>0) and ((type(m) == int) or (type(m) ==  float))):
            self.m = m 
        else:
            raise ValueError
        if((c>0) and ((type(c) == int) or (type(c) ==  float))):
            self.c = c
        else:
            raise ValueError
        if((k>0) and ((type(k) == int) or (type(k) ==  float))):
            self.k = k
        else:
            raise ValueError
    
    def period(self):
        self.T = 2 * np.pi * (self.m / self.k)**(0.5)
        return(self.T)

    def solve(self, func, y0):
        m = self.m
        c = self.c
        k = self.k
        T = self.T
        t = np.linspace(0,10*T,1000)

but im unsure where to really progress. ive tried turning the 2nd order differential equation into a lambda function like such

F = lambda X,t: [X[1], (-c) * X[1] + (-k) * X[0] + func(t)]

and then passing that into my RungeKutta function

result = RungeKutta(F, y0, t, func)
return(result)

but im not really well versed in lambda functions and am clearly going wrong somewhere.

an example input that it should pass would be something like this

####### example inputs #######
m=1
c=0.5
k=2
a = harmonicOscillator(m,c,k)
a.period()
x0 = [0,0]
tHO,xHO= a.solve(lambda t: omega0**2,x0)

would really appreciate some help. the requirments for the questions are that I have to use the above rungeKutta function, but im just kind of lost at this point thanks.

Answer 1

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: