Topic: SDOF Spring-Mass-Damper System
Subject: Numerical Methods & Mechanical
Vibration
Tools: C++, CodeBlocks, GitHub Copilot, LabPlot,
Python, VS Code, and QCAD
- 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 4th-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 2nd-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}$$
Subject to initial conditions
$$y(0) = y_0 \quad \text{and} \quad y'(0) = v_0 \tag{2}$$
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
- Underdamped: \(0 < \zeta < 1\)
- Overdamped: \(\zeta > 1\)
- 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 <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
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 matplotlib.pyplot as plt
def runge_kutta_4_system(funcs, y0, t0, t_end, h):
"""
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:
# 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
- D.I.Inman. Engineering Vibration. 2nd Ed. Prentice Hall International, Inc., New Jersey 07458. 2001.
- GitHub Copilot. AI Tools. Accessed on May 8, 2025 within the VS Code environment.
