C++ and Python for Numerical Solution to Spring-Mass-Damper Model

Topic: SDOF Spring-Mass-Damper System
Subject: Numerical Methods & Mechanical Vibration
Tools: C++, CodeBlocks, GitHub Copilot, LabPlot, Python, VS Code, and QCAD

By: Gani Comia, May 2025

• Spring-Mass-Damper Model and Tools for Numerical Solution

Numerical solutions, such as Euler's method or Runge-Kutta method, can be used to determine the free response of a damped single-degree-of-freedom system. These frequently employed numerical techniques are utilized to solve vibration problems for which closed-form analytical solutions are not readily obtainable. Runge-Kutta simulations provide a higher degree of accuracy and quality compared to those generated by the Euler method.

C++ and Python programming language will be used in this article to solve numerically the response from a Spring-Mass-Damper Model using the Runge-Kutta method of 4^th-Order. One may refer to an article in this blog entitled “Simple Harmonic Motion Solution using Runge-Kutta RK4” for the discussion of RK4 numerical method as applied to 2^nd-Order ODE.

As observed, free oscillations in most systems eventually decay and the motion ceases entirely. This phenomenon can be modelled by adding a damper to the spring-mass model of the mechanical systems. The equation of motion of the spring-mass-damper model is represented by

$$m \, y''(t) + c \, y'(t) + k \, y(t) = 0  \tag{1}$$

Where:

\(m\) – mass, \(kg\)

\(c\) – damping coefficient, \(kg/s\)

\(k\) – spring constant, \(N/m\)

\(y’’(t)\) – acceleration as a function of time, \(m/s^2\)

\(y’(t)\) – velocity as a function of time, \(m/s\)

\(y(t)\) – displacement as a function of time, \(m\)

Subject to initial conditions

$$y(0) = y_0 \quad \text{and} \quad y'(0) = v_0  \tag{2}$$

Where:

\(y_0\) – initial position, \(m\), at \(t = 0 \; s\)

\(v_0\) – initial velocity, \(m/s\), at \(t = 0 \, s\)

Equation (1) is also referred to as a damped single-degree-of-freedom (SDOF) system. In analyzing the damped SDOF system, the non-dimensional number, \(\zeta\), called the damping ratio dictates the type of its response. Damping ratio is defined in Equation (3)

$$\zeta = \frac{c}{2 \, \sqrt{k \, m}}  \tag{3}$$

Damping ratio defines the types of motion of the response as follows

1. Underdamped: \(0 < \zeta < 1\)
2. Overdamped: \(\zeta > 1\)
3. Critically Damped:  \(\zeta = 1\)

• Application

Figure 1 shows a sample application where the response in terms of displacement and velocity as a function of time are to be solved.

Figure 1. Schematic of a SDOF Spring-Mass-Damper Vibration Model.

The given sample of damped system can be represented by the Equation (4).

$$3 \, y''(t) + y'(t) + 2 \, y(t) = 0  \tag{4}$$

With initial conditions as follows

$$y(0) = 0 \quad \text{and} \quad y'(0) = 0.25  \tag{5}$$

• Numerical Solution Using C++, Codeblocks, Github Copilot, and LabPlot.

Figure 2 is the plot of the numerical solution or the response on the system in terms of displacement and velocity. The code is a modified solution through the assistance of Github Copilot (AI). Codeblocks is used as the Integrated Development Environment (IDE) for the C++. The calculation results from C++ are exported as .csv file. LabPlot is the application used to import the .csv file to make the plot as shown in Figure 2.

Figure 2. Plot of Numerical Solution to Damped System using C++ and LabPlot.

• C++ Code

#include <iostream>
#include <fstream>
#include <vector>
#include <cmath>
#include <iomanip>

// Define the system of ODEs
// y1' = y2
// y2' = f(t, y1, y2) = -(2/3)*y1-(1/3)*y2

double f1(double t, double y1, double y2) {
    return y2; // y1' = y2
}
double f2(double t, double y1, double y2) {
    // Example: y'' = -(2/3)*y1-(1/3)*y2 (Simple Harmonic Oscillator)
    return -(2.0/3.0)*y1-(1.0/3.0)*y2; // y2' = -(2/3)*y1-(1/3)*y2
}
// Runge-Kutta 4th Order Method
void rungeKutta4(double t0, double y10, double y20, double tEnd, double h, const std::string& filename) {
    double t = t0;
    double y1 = y10;
    double y2 = y20;
    // Open the CSV file for writing
    std::ofstream file(filename);
    if (!file.is_open()) {
        std::cerr << "Error: Could not open file " << filename << " for writing.\n";
        return;
    }
    // Write the header to the CSV file
    file << "t,y1,y2\n";
    // Write the initial conditions
    file << t << "," << y1 << "," << y2 << "\n";
    while (t < tEnd) {
        // Compute RK4 coefficients for y1
        double k1_y1 = h * f1(t, y1, y2);
        double k1_y2 = h * f2(t, y1, y2);
        double k2_y1 = h * f1(t + h / 2.0, y1 + k1_y1 / 2.0, y2 + k1_y2 / 2.0);
        double k2_y2 = h * f2(t + h / 2.0, y1 + k1_y1 / 2.0, y2 + k1_y2 / 2.0);
        double k3_y1 = h * f1(t + h / 2.0, y1 + k2_y1 / 2.0, y2 + k2_y2 / 2.0);
        double k3_y2 = h * f2(t + h / 2.0, y1 + k2_y1 / 2.0, y2 + k2_y2 / 2.0);
        double k4_y1 = h * f1(t + h, y1 + k3_y1, y2 + k3_y2);
        double k4_y2 = h * f2(t + h, y1 + k3_y1, y2 + k3_y2);
        // Update y1 and y2
        y1 += (k1_y1 + 2 * k2_y1 + 2 * k3_y1 + k4_y1) / 6.0;
        y2 += (k1_y2 + 2 * k2_y2 + 2 * k3_y2 + k4_y2) / 6.0;
        // Update time
        t += h;
        // Write the results to the CSV file
        file << t << "," << y1 << "," << y2 << "\n";
    }
    // Close the file
    file.close();
    std::cout << "Results saved to " << filename << "\n";
}

int main() {
    // Initial conditions
    double t0 = 0.0;    // Initial time
    double y10 = 0.0;   // Initial value of y (y1)
    double y20 = 0.25;   // Initial value of y' (y2)
    double tEnd = 20.0; // End time
    double h = 0.1;     // Step size
    // Output file name
    std::string filename = "rk4_2nd_ode_output.csv";
    // Solve the ODE and save results to the CSV file
    rungeKutta4(t0, y10, y20, tEnd, h, filename);
    return 0;
}

• Numerical Solution Using Python, Github Copilot, and VS Code.

Figure 3 is the plot for the response of the same damped system also in terms of displacement and velocity. Github Copilot (Artificial Intelligence) is used to generate the Python script with minor modification. VSCode is used as IDE to run the script and plot the results.

Figure 3. Plot of Numerical Solution to Damped System using Python Script.

• Python Script

import numpy as np
import matplotlib.pyplot as plt

def runge_kutta_4_system(funcs, y0, t0, t_end, h):
    """

    Solve a system of first-order ODEs using the 4th-order Runge-Kutta method (RK4).
    Parameters:
        funcs (list): List of functions representing the system of ODEs.
        y0 (list): Initial values for the system [y1(0), y2(0), ...].
        t0 (float): Initial time.
        t_end (float): End time.
        h (float): Step size.
    Returns:
        tuple: Arrays of t and y values for each variable.
    """
    t_values = np.arange(t0, t_end + h, h)
    y_values = np.zeros((len(t_values), len(y0)))
    y_values[0] = y0
    for i in range(1, len(t_values)):
        t = t_values[i - 1]
        y = y_values[i - 1]
        k1 = [h * f(t, *y) for f in funcs]
        k2 = [h * f(t + h / 2, *(y + np.array(k1) / 2)) for f in funcs]
        k3 = [h * f(t + h / 2, *(y + np.array(k2) / 2)) for f in funcs]
        k4 = [h * f(t + h, *(y + np.array(k3))) for f in funcs]
        y_values[i] = y + (np.array(k1) + 2 * np.array(k2) + 2 * np.array(k3) + np.array(k4)) / 6
    return t_values, y_values

# Example application
if __name__ == "__main__":
    # Define the system of ODEs for 3y'' + y' + 2y = 0
    # Let y1 = y and y2 = -(2.0/3.0)*y1-(1.0/3.0)*y2, then:

    # y1' = y2
    # y2' = -(2.0/3.0)*y1-(1.0/3.0)*y2
    def f1(t, y1, y2):
        return y2
    def f2(t, y1, y2):
        return -(2.0/3.0)*y1-(1.0/3.0)*y2

    # Initial conditions
    y0 = [0, 0.25]  # y(0) = 0, y'(0) = 0.25
    t0 = 0       # Initial time
    t_end = 20   # End time
    h = 0.1      # Step size
    # Solve the system
    t_values, y_values = runge_kutta_4_system([f1, f2], y0, t0, t_end, h)
    # Extract y1 (the solution to the original ODE)
    y1_values = y_values[:, 0]
    y2_values = y_values[:, 1]

    # Plot the solution
    plt.figure(figsize=(8, 6))
    plt.plot(t_values, y1_values, linewidth=3, label="y(t)", color="blue")
    plt.plot(t_values, y2_values, linestyle='--', linewidth=2, label="y'(t)", color="red")
    plt.title("Using Python to Solve 2nd-Order ODE with RK4")
    plt.xlabel("Time, t")
    plt.ylabel("Displacement, y(t) & Velocity, y'(t)")
    plt.grid(True)
    plt.ylim(-0.3, 0.3)
    plt.legend()
    plt.show()

The damping ratio for the given application can be solved as shown

$$\zeta = \frac{c}{2 \, \sqrt{k \, m}} = \frac{1}{2 \sqrt{(2)(3)}} = 0.204$$

With the damping ratio, \(\zeta = 0.204\), the motion of the response of the system is considered underdamped.

C++ and Python programming languages using Codeblocks and VS Code as IDEs, and LabPlot for plotting application are some powerful engineering tools to analyze a system of mechanical vibration. AI through the use of Github Copilot provided the applicable codes for the solution.

Feel free to comment for inquiry, clarification, correction or suggestion for improvement. Drop your email to make a request to the author.

References

1. D.I.Inman. Engineering Vibration. 2^nd Ed. Prentice Hall International, Inc., New Jersey 07458. 2001.
2. GitHub Copilot. AI Tools. Accessed on May 8, 2025 within the VS Code environment. 

Solving an Undamped System Under Harmonic Force Using Laplace Transform

Topic: Harmonically Excited Vibration
Subject: Mechanical Vibration
Tool: Euler Math Toolbox, QCAD, & Scilab

By: Gani Comia, Mar. 2025

• Undamped System Under Harmonic Force - Equation of Motion

The undamped system under harmonic force is a dynamic system subjected to external force or excitation. Consider the harmonic excitations of the form

$$F(t) = F_0 \cos(\omega t) \tag{1}$$

Where:
\(F(t)\) – excitation of forcing function
\(F_0\) – constant static force
\(\omega\) – excitation angular frequency
\(t\) - time

If the force \(F(t)\) acts on an undamped spring-mass system, the equation of motion is

$$m\,y'' + k\,y = F(t) \tag{2}$$

or

$$m\,y'' + k\,y = F_0 \cos(\omega t) \tag{3}$$

Where:
\(m\) – mass
\(k\) – spring constant
\(y’’\) – acceleration of mass
\(y\) – displacement from rest position

Under a harmonic excitation, the response of the system will also be harmonic. If the frequency of excitation \(\omega\) coincides with the natural frequency \(\omega_n\) of the system, the response of the system will be very large. The condition, \(\omega = \omega_n\), is what we called resonance. Resonance is to be avoided to prevent system’s failure as in this case mechanical damage.

$$\omega = \omega_n \quad \left(\omega_n = \sqrt{k / m}\right) \tag{4} $$

Where:
\(\omega_n\) – natural angular frequency

• Laplace Transform of Second-Order Linear Differential Equation

Laplace transform is one of the methods to solve for the response for the undamped system under harmonic force. Consider a general second-order linear differential equation

$$a\,\frac{d^2 y}{d t^2} + b\, \frac{dy}{dt} + c\,y = u(t) \quad \left(t \geq 0 \right) \tag{5}$$

Subject to initial conditions

$$t = 0, \quad y(0) = y_0, \quad \frac{dy}{dt}(0) = \nu_0 \tag{6}$$

Such differential equations may model the dynamics of some system for which the variable \(y(t)\) determines the response of the system to the forcing or excitation term \(u(t)\). The terms system input and system output are also frequently used for \(u(t)\) and \(y(t)\) respectively.

• Laplace and Inverse Laplace Transform of an Undamped System Under Harmonic Force

Taking the Laplace transform of each term of Equation (3)

$$m\, \mathcal{L} \left\{ \frac{d^2 y}{d t^2} \right\} + k\, \mathcal{L} \{y\} = F_0 \, \mathcal{L} \{ \cos(\omega t)\} \tag{7}$$

With initial conditions

$$t = 0, \quad y(0) = \frac{F_0}{k - m\,\omega^2}, \quad \frac{dy}{dt}(0) = 0 \tag{8}$$

Leads to

$$m \, \left( s^2 \, Y(s) - \frac{s\,F_0}{k - m \, \omega^2} \right) + k\,Y(s) = \frac{s \, F_0}{\omega^2 + s^2} \tag{9}$$

By rearranging and solving for \(Y(s)\) gives

$$Y(s) = - \frac{s\,F_0}{m\,\omega^4 + \omega^2 \, (m\,s^2 - k) - k\,s^2} \tag{10}$$

The inverse Laplace transform solves the system output \(y(t)\) for \(\omega > 0\)

$$\mathcal{L^{-1}} \{Y(s)\} = \mathcal{L^{-1}} \left\{ - \frac{s\,F_0}{m\,\omega^4 + \omega^2 \, (m\,s^2 - k) - k\,s^2} \right\} \quad (\omega > 0) \tag{11}$$

The steady-state response \(y(t)\) is

$$y(t) = \frac{F_0 \, \cos(\omega \, t)}{m \, \omega^2 - k} \tag{12}$$

The response of the system can be identified to be of three types.

Case 1. When \(0 < {\omega}/{\omega_n} < 1\), the harmonic response of the system \(y(t)\) is said to be lagging the force.

Case 2. When \({\omega}/{\omega_n} > 1\), the response of the system to a harmonic force of very high frequency is close to zero. The response of the system is said to be leading the force.

Case 3. When \({\omega}/{\omega_n} = 1\), the response is called resonance and the amplitude of \(y(t)\) becomes infinite. This phenomenon can result on large displacements and stresses.

Consider a sample application illustrated in Figure 1, let us analyze the response \(y(t)\) in terms of excitation \(F(t)\) based on the Case 1 and 2 responses. Case 3 will not be possible to plot. For this let us assume an excitation frequency, \(\omega\), of less than and more than the natural frequency, \(\omega_n\).

Figure 1. Application of an Undamped System under Harmonic Force

Below is the Scilab script for plotting both the system input and output for the illustrated problem.

• Scilab Script

// Copyright (C) 2025 - Gani Comia
// Date of creation: 22 Mar 2025
// Undamped System under Harmonic Force using Laplace Transform
clear;clc;
// (1) solution y(t) from m y''(t) + k y(t) = F0 cos(wt)
function y=position(t)
    y = -(F0.*cos(omega.*t))./((m.*omega.^2)-k)
endfunction

// (2) excitation of forcing function
function F=force(t)
    F = F0.*cos(omega.*t)
endfunction

// (3) primary parameters from m y''(t) + k y(t) = F0 cos(wt)
m = 100/9.8;                        // kg, mass (where W = 100 N)
k = 2000;                           // N/m, spring constant
F0 = 25;                            // N, constant static force

// (4) secondary parameters
omega_n = sqrt(k./m);               // Hz, natural angular frequency
mprintf("omega_n (Hz) = %3.1f\n", omega_n)

clf;
// (5) figure properties
fig = gcf()
fig.figure_size = [700,800];

// (6) for excitation or forcing function
omega = omega_n
omega_t = linspace(0,2*%pi,100);
t = omega_t./omega;

// (6-1) plot of excitation force
subplot(3,1,1)
F = force(t)
plot(t,F,"m-","linewidth",4)
title("$\Large\text{Excitation Function}\quad F(t)=F_0\cos(\omega t)$")
ylabel("$\Large \mathbf{F(t)}$")
xlabel("$\Large \text{time,}\; t$")
legend("$\LARGE\omega = \omega_n$",4,with_box=%F)
xgrid(3,0)
// (6-2) axes properties
ax = gca()
ax.x_location = 'origin' 
ax.data_bounds = [0 -30; 0.5 30]

// (7) for Case 1. 0 < omega < omega_n
omega = 0.95*omega_n
omega_t = linspace(0,2*%pi,100);
t = omega_t./omega;

// (7-1) plot of a sampe case 1
subplot(3,1,2)
y = position(t)
plot(t,y,"b-","linewidth",4)
title("$\Large\text{Displacement at}\quad \mathbf{\omega < \omega_n}$")
ylabel("$\Large \mathbf{y(t)}$")
xlabel("$\Large \text{time,}\; t$")
legend("$\LARGE \omega = \text{0.95} \; \omega_n$",4,with_box=%F)
xgrid(3,0)
ax = gca()
ax.x_location = 'origin'
ax.data_bounds = [0 -0.2; 0.5 0.2]

// (8) for Case 2. omega > omega_n
omega = 1.05*omega_n
omega_t = linspace(0,2*%pi,100);
t = omega_t./omega;

// (8-1) plot of a sampe case 2
subplot(3,1,3)
y = position(t)
plot(t,y,"r-","linewidth",4)
title("$\Large\text{Displacement at}\quad \mathbf{\omega > \omega_n}$")
ylabel("$\Large \mathbf{y(t)}$")
xlabel("$\Large \text{time,}\; t$")
legend("$\LARGE \omega = \text{1.05} \; \omega_n$",1,with_box=%F)
xstring(0.31,-0.2,"https://gani-mech-toolbox.blogspot.com")
xgrid(3,0)
ax = gca()
ax.x_location = 'origin'
ax.data_bounds = [0 -0.2; 0.5 0.2] 

• Scilab Output (Figure 2)

Figure 2. Forcing Function and Responses to Undamped System with Harmonic Excitation

The harmonic steady-state response \(y(t)\) for case 1 is said to be lagging with the forcing function \(F(t)\). As for the case 2 the response \(y(t)\) is said to be leading with external harmonic force \(F(t)\).

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

References

1. Singiresu S. Rao. Mechanical Vibrations. 2^nd Ed. Addison-Wesley Publishing Company. 1990.
2. G. James, D. Burley, P. Dyke, J. Searl, N. Steele, J. Wright. Advanced Modern Engineering Mathematics. Addison-Wesley Publishing Company, Inc. 1993.

Dynamic Amplitude Response of an Undamped System under Harmonic Force

Topic: Harmonically Excited Vibration

Subject: Mechanical Vibration

Tool: Scilab & QCAD

By: Gani Comia, Mar. 2025

The undamped system under harmonic force is a dynamic system subjected to external force or excitation. This excitation is called the forcing or excitation function. The excitation function is usually time-dependent and may be harmonic or nonharmonic. 

This article presents a dynamic response of a single degree of freedom system under the harmonic motions of the form

$$F(t) = F_0\,\cos(\omega\,t + \phi) \tag{1}$$

Where:
\(F(t)\) – excitation of forcing function
\(F_0\) – constant static force
\(\omega\) – excitation angular frequency
\(t\) - time
\(\phi\) – phase angle of harmonic excitation

If the force \(F(t)\) acts on an undamped spring-mass system with \(\phi = 0\), the equation of motion is

$$m\,{y''} + k\,y = F(t) \tag{2}$$

Or

$$m\,{y''} + k\,y = F_0 \, \cos(\omega\,t) \tag{3}$$

Where:
\(m\) – mass
\(k\) – spring constant
\({y''}\) – acceleration of mass
\(y\) – displacement from rest position

The dynamic amplitude of the undamped spring-mass system is a vibration response which refers to the maximum displacement under dynamic conditions subjected to an oscillatory force or motion. This can be expressed as

$$Y = \large \frac{\delta_{st}}{1 - \left(\frac{\omega}{\omega_n}\right)^2} \tag{4}$$

Where:
\(Y\) – dynamic amplitude
\(\delta_{st}\) – static deflection
\(\omega\) – excitation angular frequency
\(\omega_n\) – natural angular frequency

The static deflection and natural angular frequency are defined as follows

$$\delta_{st} = \frac{F_0}{k}\;, \quad \omega_n = \sqrt{\frac{k}{m}} \tag{5}$$

The response of the system can be identified to be of three types.

Case 1. When \(\mathbf{ 0 < {\omega}/{\omega_n} < 1 }\), the harmonic response of the system \(y(t)\) is said to be lagging the external force.

Case 2. When \(\mathbf{ \omega / \omega_n > 1 }\), the response of the system to a harmonic force of very high frequency is close to zero. The response of the system is said to be leading the force.

Case 3. When \(\mathbf{ \omega / \omega_n = 1 }\), the response is called resonance and the amplitude of \(y(t)\) becomes infinite. This phenomenon can result on large displacements and stresses.

• Application

Consider a sample application illustrated in Figure 1, let us analyze its dynamic amplitude response in terms of excitation frequency based on the three types of responses.

Figure 1. Undamped System under Harmonic Force – Sample Problem

Below is the calculation of dynamic amplitude for the illustrated application with the use of Scilab script.

• Scilab Script

// Copyright (C) 2025 - Gani Comia
// Date of creation: 20 Mar 2025
// Dynamic Amplitude from Undamped System under Harmonic Force
clear;clc;
// (1) dynamic amplitude function
function Y = dynAmp(w)
    Y = delta./(1 - (w./wn).^2)
endfunction

// (2) primary parameters
m = 100/9.8                     // kg, mass
k = 2000                        // N/m, spring stiffness or constant
F0 = 25                         // N, static force (constant)
delta = F0./k                   // m, static amplitude or deflection

// (3) secondary parameters
wn = sqrt(k./m)                 // Hz, natural angular frequency
mprintf("Natural frequency, w_n = %3.1f",wn)

// (4) frequency domain
w_l = linspace(0,wn-0.05,100)    // left side of natural frequency
w_r = linspace(wn+0.05,60,100)   // right side of natural frequency

// (5) calculation of dynamic amplitude
Y_l = dynAmp(w_l)
Y_r = dynAmp(w_r)

clf;
// (6) figure properties
fig = gcf()
fig.figure_size = [700,700]

// (7) plot of results
plot(w_l,Y_l,"b-","linewidth",3.5)
plot(w_r,Y_r,"m-","linewidth",3.5)
plot([wn wn],[-1,1],"r--","linewidth",2)
note1 = "Vibration Response of an Undamped System w/ Harmonic Force"
title(note1,"fontsize",3.75)
ylabel("$\Large\text{Dynamic Amplitude (m),}\;\;\mathbf{Y(\omega)}$")
xlabel("$\Large\text{Mechanical Frequency (Hz),}\;\;\mathbf{\omega}$")
xgrid(3,0)

// (8) additional plot information
note2 = "$\LARGE m\ddot{y}(t) + k y(t) = F_0\,cos(\omega\,t)$"
xstring(34,0.425,note2)
note3 = "$\LARGE Y(\omega\;<\;\omega_n)$"
note4 = "$\LARGE Y(\omega\;>\;\omega_n)$"
note5 = "$\LARGE Y(\omega\;=\;\omega_n)$"
legend(note3,note4,note5,with_box=%F)
note6 = "$\Large \text{Resonance}$"
xstring(16,0.4,note6,-90)
xstring(13.5,-0.8,note6,-90)

// (9) axes properties
ax = gca()
ax.x_location = 'origin' 
ax.data_bounds = [0 -1; 60 1];

• Scilab Output (Figure 2)

Figure 2. Vibration Response to Undamped System with Harmonic Excitation

Figure 2 shows the dynamic amplitude in three types of responses to harmonic excitation. For the illustrated problem, a suitable forcing frequency \(\omega\) can be chosen based on the requirements. Resonance frequency shows an infinite displacement.

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

Reference

1. Singiresu S. Rao. Mechanical Vibrations. 2^nd Ed. Addison-Wesley Publishing Company. 1990.

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.

Mechanical Vibration Measurement using Accelerometer and Scilab

Topic: Vibration Detection and Measurement;

Subject: Mechanical Vibration;

Tool: Phyphox and Scilab;

By: Gani Comia, Nov. 2024;

A vibrating object is said to be moving and accelerating.  Accelerometer sensors, a device that measures the acceleration of an object, are used to detect the mechanical vibration. The Phyphox app, an application that can be installed in an Android or IOS phone, can access your phone’s sensor to measure acceleration. It has the capability to store and download the acceleration data in the .csv file used for numerical analysis.

Here is the Scilab script for visualization of the example downloaded data from the Phyphox app.  The app’s acceleration data is organized into five columns: time, accel x-axis, accel y-axis, accel z-axis, and accel abs.  An example plot of the time and accel abs (absolute acceleration) and its script are presented.

Scilab Script

// Plotting data from accelerometer using Phyphox & Scilab
// by: Gani Comia, Aug. 2024
clear;clc;

// extracting data
clear importdata;

function [data]=importdata(filename)
    data = csvRead(filename, ",", ".", "double")
endfunction

[A] = importdata("accelerometerData.csv");

// assigning variable name
time = A(:,1);
accel_x = A(:,2);
accel_y = A(:,3);
accel_z = A(:,4);
accel_abs = A(:,5);

maxTime = max(time);
maxAccel = max(accel_abs);
disp(maxTime,maxAccel)

// plotting acceleration data
clf;
f=gcf();
f.figure_size=[700,600];
plot(time,accel_abs,"red")
title("Vibration Analysis using Abs Acceleration","fontsize",4)
xlabel("$time (sec)$","fontsize",4)
ylabel("$acceleration\;(m/s^2)$","fontsize",4)

// displaying information
format("v",6)
xstring(50,maxAccel,["Max Accel:", string(maxAccel)])
xstring(maxTime,0.15,["Max Time:", string(maxTime)],-90)

Executing the script will give you the plot for visualization of the acceleration data.

Fig. 1. Vibration Analysis using Absolute Acceleration.

The Scilab simulation tool facilitates visualization of numerical data, as in this case the acceleration of the vibrating object. Further numerical analysis to determine the velocity, displacement, and frequency of vibration can be done as well using the tool.

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