Simulating Mechanical Systems

Mechanical systems form the backbone of countless engineering applications — from the suspension of an automobile absorbing road bumps to a precision machine tool vibrating under cutting forces. Simulating these systems numerically allows engineers to predict behavior, tune parameters, and avoid destructive resonance conditions long before physical prototypes are built. The spring-mass-damper system is the canonical model for studying oscillatory mechanical behavior. Its elegance lies in the fact that a single second-order ordinary differential equation (ODE) captures the interaction of three fundamental physical phenomena: inertia, energy dissipation, and elastic restoring force. Mastering the numerical simulation of this system — and understanding how its analytical solution compares to numerical integration — builds the intuition needed to tackle far more complex mechanical models.

The governing equation of motion for a spring-mass-damper system is:

m·ẍ + c·ẋ + k·x = F(t)

Here, m is the mass (kg), c is the viscous damping coefficient (N·s/m), k is the spring stiffness (N/m), x is the displacement from equilibrium (m), ẋ is velocity (m/s), ẍ is acceleration (m/s²), and F(t) is an externally applied force that may vary with time. Each term has clear physical meaning: m·ẍ is the inertial force resisting acceleration, c·ẋ is the damping force opposing velocity (think of a shock absorber), and k·x is the spring force pulling the mass back toward equilibrium. When these three forces balance the applied force, the motion is fully described.

Two dimensionless and dimensional quantities derived from m, c, and k characterize the system's behavior completely. The natural frequency is defined as:

ωₙ = √(k / m)    [radians per second]

This is the rate at which the undamped system would oscillate if displaced from rest and released. A stiffer spring or lighter mass increases ωₙ, meaning faster oscillations. A heavier mass or softer spring lowers it. The damping ratio is defined as:

ζ = c / (2·√(m·k))    [dimensionless]

The denominator 2·√(m·k) is called the critical damping coefficient, the smallest value of c that prevents oscillation. The damping ratio ζ compares the actual damping to this critical value. Together ωₙ and ζ collapse the three separate parameters m, c, k into a form that immediately reveals the qualitative character of the system's response.

Because numerical integration methods (Euler, Runge-Kutta, etc.) are designed to advance first-order systems, the single second-order ODE must be converted into two coupled first-order ODEs. The standard approach introduces a state vector:

y = [x, v]    where v = ẋ

Differentiating each component:

dy₁/dt = v
dy₂/dt = [F(t) - c·v - k·x] / m

The second equation comes directly from solving the equation of motion for acceleration ẍ. This pair of first-order equations is the form handed to any numerical integrator. A Python implementation of the derivative function for free vibration (F(t) = 0) looks like this:

def spring_mass_damper(t, y, m, c, k):
    x, v = y
    dxdt = v
    dvdt = (-c * v - k * x) / m
    return [dxdt, dvdt]

This function receives the current time t, the current state y = [x, v], and the system parameters, then returns the time derivatives of each state variable. The numerical integrator calls this function repeatedly to march the solution forward in time.

The three qualitatively distinct response types — underdamped, critically damped, and overdamped — arise depending on the value of ζ.

An underdamped system (ζ < 1) oscillates. The mass overshoots equilibrium, reverses direction, overshoots again on the other side, and gradually the amplitude decays toward zero. The oscillation does not occur at the natural frequency ωₙ but at the slightly lower damped natural frequency:

ωd = ωₙ · √(1 - ζ²)

For lightly damped systems (ζ ≪ 1), ωd ≈ ωₙ, and the difference is negligible. For moderate damping — say ζ = 0.5 — the damped frequency is noticeably reduced: ωd = ωₙ·√(0.75) ≈ 0.866·ωₙ. Engineers designing vibration isolators, vehicle suspensions, or seismic dampers almost always operate in the underdamped regime, accepting some oscillation in exchange for better isolation performance.

A critically damped system (ζ = 1) returns to equilibrium as quickly as possible without oscillating at all. This is the fastest non-oscillatory response achievable. Instrument pointer mechanisms (galvanometers, analog meters) and door closing mechanisms are often designed near critical damping to avoid the needle or door panel bouncing back and forth.

An overdamped system (ζ > 1) also returns to equilibrium without oscillation, but does so more slowly than a critically damped system. The excess damping fights the restoring force as well as velocity, creating a sluggish return. Very thick viscous fluids, or deliberately over-stiffened shock absorbers, produce overdamped behavior. In practice, overdamping is rarely desirable because performance is sacrificed without any benefit over critical damping.

To appreciate these distinctions quantitatively, consider a system with m = 1 kg and k = 100 N/m, giving ωₙ = 10 rad/s. Choosing ζ = 0.1, 1.0, and 2.0 requires c values of 2, 20, and 40 N·s/m respectively. Simulating all three from the same initial condition x(0) = 1 m, v(0) = 0 reveals the three response shapes starkly on a single plot.

The analytical solution for free vibration in the underdamped case provides the exact trajectory for comparison against numerical results. Starting from the equation of motion with F(t) = 0, the solution is:

x(t) = A · e^(-ζ·ωₙ·t) · cos(ωd·t + φ)

The constants A (amplitude) and φ (phase angle) are determined by initial conditions. For initial displacement x₀ and initial velocity v₀:

A = √( x₀² + ((v₀ + ζ·ωₙ·x₀) / ωd)² )
φ = arctan( -(v₀ + ζ·ωₙ·x₀) / (ωd·x₀) )

In NumPy this analytical solution is computed element-wise over a time array with no looping required:

import numpy as np
import matplotlib.pyplot as plt

m, c, k = 1.0, 2.0, 100.0
x0, v0  = 1.0, 0.0

wn  = np.sqrt(k / m)
zeta = c / (2 * np.sqrt(m * k))
wd  = wn * np.sqrt(1 - zeta**2)

A   = np.sqrt(x0**2 + ((v0 + zeta*wn*x0)/wd)**2)
phi = np.arctan2(-(v0 + zeta*wn*x0), wd*x0)

t   = np.linspace(0, 5, 1000)
x_analytical = A * np.exp(-zeta * wn * t) * np.cos(wd * t + phi)

plt.plot(t, x_analytical, label='Analytical')
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.legend()
plt.show()

Overlaying the numerical solution on the same axes immediately reveals any discrepancies. A coarse time step in the numerical integrator will cause the two curves to drift apart, providing a concrete visual demonstration of truncation error. The maximum absolute error between the two solutions is computed as:

max_error = np.max(np.abs(x_numerical - x_analytical))

This scalar metric quantifies how accurately the chosen time-stepping scheme reproduces the true physics. Halving the time step and observing how the maximum error decreases verifies the convergence order of the method — first order for Euler, fourth order for RK4.

Forced vibration extends the model by adding a harmonic driving force to the right-hand side:

m·ẍ + c·ẋ + k·x = F₀·sin(Ω·t)

Here F₀ is the force amplitude (N) and Ω is the forcing frequency (rad/s), which is an independently tunable parameter distinct from ωₙ. The system now has two frequencies competing: the natural frequency ωₙ at which it wants to oscillate, and the forcing frequency Ω at which energy is being pumped in.

Resonance is the condition Ω ≈ ωₙ. At resonance, energy is added to the system in perfect synchrony with its natural oscillation, and the response amplitude grows dramatically. In a lightly damped system the amplitude at resonance can be many times the static deflection F₀/k. The dynamic amplification factor (or frequency response function magnitude) at steady state is:

|H(Ω)| = 1 / √((1 - r²)² + (2·ζ·r)²)
where r = Ω / ωₙ

At r = 1 (resonance), |H| = 1/(2ζ). For ζ = 0.05 (lightly damped), |H| = 10 — the response amplitude is ten times the static deflection. For ζ = 0.5, |H| = 1.0 at resonance. This makes quantitative the intuition that light damping produces sharp, dangerous resonance peaks.

A frequency response curve is produced by sweeping Ω from near zero to well above ωₙ, running the simulation to steady state at each forcing frequency, measuring the peak amplitude, and plotting amplitude versus Ω. The resulting curve shows a peak near ωₙ whose height is controlled by ζ. The derivative function for forced vibration is a simple extension:

def forced_smd(t, y, m, c, k, F0, Omega):
    x, v = y
    force = F0 * np.sin(Omega * t)
    dxdt  = v
    dvdt  = (force - c * v - k * x) / m
    return [dxdt, dvdt]

The complete numerical implementation ties all these concepts together. The state vector [x, v] is initialized from initial conditions, and a time-stepping loop advances the solution using the chosen integrator. Using SciPy's solve_ivp with the RK45 method (adaptive 4th/5th-order Runge-Kutta) is the recommended approach for accuracy and convenience:

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

# System parameters
m, c, k = 1.0, 2.0, 100.0
F0, Omega = 5.0, 10.0        # forced vibration
x0, v0   = 0.0, 0.0          # initial conditions at rest

wn   = np.sqrt(k / m)
zeta = c / (2 * np.sqrt(m * k))
wd   = wn * np.sqrt(1 - zeta**2)

def eom(t, y):
    x, v = y
    dxdt = v
    dvdt = (F0 * np.sin(Omega * t) - c * v - k * x) / m
    return [dxdt, dvdt]

t_span = (0, 10)
t_eval = np.linspace(*t_span, 2000)
sol    = solve_ivp(eom, t_span, [x0, v0], t_eval=t_eval, method='RK45')

x_num = sol.y[0]   # displacement history
v_num = sol.y[1]   # velocity history

plt.figure(figsize=(10, 4))
plt.plot(sol.t, x_num, label='Numerical (RK45)')
plt.axhline(0, color='gray', linewidth=0.5)
plt.xlabel('Time (s)')
plt.ylabel('Displacement (m)')
plt.title('Forced Spring-Mass-Damper Response')
plt.legend()
plt.tight_layout()
plt.show()

When Omega is set equal to wn in the code above, the resonance condition is triggered. Starting from rest with light damping, the displacement envelope grows steadily — demonstrating visually why engineers must identify and detune resonant frequencies. Adding the analytical free-vibration solution as an overlay (when F0 = 0) immediately validates whether the numerical scheme is accurately tracking the true trajectory, and computing np.max(np.abs(sol.y[0] - x_analytical)) provides the numerical accuracy metric.

Together, these components — the equation of motion, the state-space reformulation, the damping ratio classification, the analytical solution, and the forced vibration model — form a complete framework for simulating and understanding fundamental mechanical dynamics. Every more complex system, whether a multi-degree-of-freedom structure, a rotating shaft, or a nonlinear pendulum, is built on exactly these foundations.