Topic: Oscillatory Motion and Numerical Solution of ODE
Subject: Physics and Numerical Methods
Tool: Scilab
The 4th -order Runge-Kutta methods (RK4) is also used for numerical technique of approximating solution to the initial value problem (IVP) of 2nd- 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_2\;=\;y' \tag{3}$$
Then, the nonhomogeneous LDE can be re-written as a system
of 1st-order ODE in Equation (4) and (5):
$$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 2nd-order ODE to a system
of 1st-order to use the RK4 method in finding its approximate
solution.
The SDOF spring-mass model oscillation can be described by
the basic equation of SHM shown in Equation (6).
Where:
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.
