Simple Harmonic Motion Solution using Runge-Kutta RK4 Method

Topic: Oscillatory Motion and Numerical Solution of ODE

Subject: Physics and Numerical Methods

Tool: Scilab

By: Gani Comia, Jan. 2025

The 4^th -order Runge-Kutta methods (RK4) is also used for numerical technique of approximating solution to the initial value problem (IVP) of 2^nd- order ordinary differential equation (ODE). This article presents the solution to the basic equation of simple harmonic motion (SHM) using RK4 algorithm in Scilab.

SHM is a form of nonhomogeneous linear differential equation (LDE) of second-order as shown in Equation (1).

$$y''\;+\;p(t)\,y'\;+\;q(t)\,y\;=\;r(t) \tag{1}$$ 

In order to solve a given higher-order differential equation, it is required to be converted to a system of first-order differential equations as presented below.

Let us define the following dependent variables, \(y_1\) and \(y_2\):

$$y_1\;=\;y \tag{2}$$

$$y_2\;=\;y' \tag{3}$$

Then, the nonhomogeneous LDE can be re-written as a system of 1^st-order ODE in Equation (4) and (5):

$$y'_1\;=\;y_2 \tag{4}$$

$$y'_2\;=\;-p(t)\,y_2\;-\;q(t)\,y_1\;+\;r(t) \tag{5}$$

With that presented, let us define and illustrate a problem using the basic equation of SHM and convert its 2^nd-order ODE to a system of 1^st-order to use the RK4 method in finding its approximate solution.

Figure 1. Single-Degree of Freedom (SDOF) Spring-Mass Model Oscillation Problem

The SDOF spring-mass model oscillation can be described by the basic equation of SHM shown in Equation (6).

$$m\,\frac{d^2y}{dt^2}\,+\,k\,y\,=\,0, \quad y(t_0)=y_0\;\;\text{and}\;\;y'(t_0)=y'_0 \tag{6}$$

Where:

\(m\) – mass

\(k\) – spring constant

\(t\) – time

\(t_0\) - initial time

\(y_0\) – initial displacement

\(y’_0\) – initial velocity

The solution to the ODE of SHM can be approximated using the Runge-Kutta (RK4) Method. The RK4 will be using the same algorithm presented in the previous article “First-Order IVP ODE Solution using Runge-Kutta RK4 Method” with the link shown: https://gani-mech-toolbox.blogspot.com/2025/01/first-order-ivp-ode-solution-using.html. Figure 2 shows the solution to the problem illustrated in Figure 1.

Figure 2. Solution to Spring-Mass Model System with RK4 using Scilab Script

The Scilab script to solve the given problem is presented below. Each block of the code is commented comprehensively for better understanding of the program sequence by the reader.

Scilab Script:

clear;clc;
// single-degree of freedom spring-mass model oscillator system
// m*y'' + k*y = 0  or  y'' + (k/m)*y = 0
// nonhomogeneous linear differential equation (LDE) of second-order
// y'' + p(t)y' + q(t)y = g(t) 
// where: p(t)=0; q(t)=k/m; g(t)=0

// (1) 2nd-order LDE re-written in the system of 1st-order ODEs
function dydt=f(t, y)
    m = 1;                  // kg, mass of the object
    k = 20;                 // N/m, spring constant
    p = 0;                  // coefficient of y'
    q = k/m;                // coefficient of y
    g = 0;                  // forcing function
    dydt = [y(2); -p*y(2) - q*y(1) + g];
endfunction

As you may compare the \(rk4\) function in block (2) to the previous article, the final time, \(t_f\), is used as one of the input arguments instead of the number of intervals, \(n\). Both \(rk4\) functions can be used to come up with the same solution.

// (2) RK4 function or algorithm for 1st-order ODEs
function [t, y]=rk4(f, t0, y0, h, tf)
    n = round((tf - t0)./h);
    t(1) = t0;          // initial time, t0
    y(:, 1) = y0;       // y(t0)=y0 & y'(t0)=y0' in column vector

    for i = 1:n
        // approximate values of slopes of solution curves
        k1 = h * f(t(i), y(:, i));
        k2 = h * f(t(i) + h / 2, y(:, i) + k1 / 2);
        k3 = h * f(t(i) + h / 2, y(:, i) + k2 / 2);
        k4 = h * f(t(i) + h, y(:, i) + k3);
        // estimates of numerical solution
        y(:, i + 1) = y(:, i) + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
        t(i+1) = t(i) + h;
    end
endfunction

// (3) initial conditions and parameters
t0 = 0;             // sec, initial time
tf = 5;             // sec, final time
h = 0.01;           // sec, time step size
y0 = [0.05; 0];     // initial disp & velocity: y(0)=0.05, y'(0)=0

// (4) solution to the ODE by calling RK4 defined function
[t, y] = rk4(f, t0, y0, h, tf);

// (5) solution plot
clf;
plot(t, y(1, :),"b-","linewidth",4);       // displacement plot
plot(t, y(2, :),"r--","linewidth",1.8);    // velocity plot
xlabel("$\Large t \quad \text{(time)}$");
ylabel("$\Large y(t)\;\;\text{or}\;\;\frac{dy(t)}{dt}$");
note1 = "Spring-Mass Model System using Runge-Kutta RK4";
title(note1,"fontsize",3.5);
note2 = ["$\Large y(t)$","$\Large \frac{dy(t)}{dy}$"];
legend(note2, with_box = %F);
note3 = "https://gani-mech-toolbox.blogspot.com"
xstring(0.1,0.35,note3);
xgrid(3,1);

// (6) other plot information
disp(max(y(1,:)), min(y(1,:)))
disp(max(y(2,:)), min(y(2,:)))
ax = gca();
ax.data_bounds = [0 -0.4; 5 0.4];

// (7) natural frequency
format(5);
k = 20;                     // N/m, spring constant
m = 1;                      // kg, mass of object
omega = sqrt(k/m);          // Hz, natural frequency
disp(omega);
note4 = ["$\LARGE\mathbf{\omega_n}\;=\;\text{4.5 Hz}$"];
xstring(0.1,0.25,note4)

The Scilab script in this article can be one of your references to solve any engineering problems that can be described by the nonhomogeneous linear differential equation (LDE) of second-order.

Feel free to comment for inquiry, clarification, or suggestion for improvement. Drop your email to request the soft copy of the file.