As seen in the previous article in this series, everything is a continuous multi-DoF SMD system. However, practical engineering models are usually built by isolating the prominent masses of a structure and identifying the springs and dampers that connect them. This approach — replacing a continuous system with a finite number of discrete masses, springs, and dampers — is called a lumped SMD model.
A Simple Horizontal 2DoF SMD System
To demonstrate the power of simulating SMD systems, consider two masses sliding horizontally on a frictionless surface, each constrained by springs and dampers. Because there are two masses each free to move in one direction, the system has two degrees of freedom — a 2DoF system.
A horizontally sliding 2-mass (2DoF) SMD system constrained on both sides. The masses are assumed to slide without friction.
The Model
Applying Newton's second law \(F = ma\) to each mass and summing all forces gives the equations of motion. For mass 1, three forces act: the wall spring and damper on the left, and the shared spring and damper coupling it to mass 2. Taking rightward displacement as positive:
\[m_1\ddot{x}_1 = -k_1 x_1 - c_1\dot{x}_1 - k_2(x_1 - x_2) - c_2(\dot{x}_1 - \dot{x}_2)\]
Collecting terms by kinematic quantity:
\[m_1\ddot{x}_1 + (c_1 + c_2)\dot{x}_1 - c_2\dot{x}_2 + (k_1 + k_2)x_1 - k_2 x_2 = 0\]
For mass 2, the shared coupling element pulls it toward mass 1 while the right-side spring and damper restore it toward the fixed wall:
\[m_2\ddot{x}_2 = k_2(x_1 - x_2) + c_2(\dot{x}_1 - \dot{x}_2) - k_3 x_2 - c_3\dot{x}_2\]
Collecting terms:
\[m_2\ddot{x}_2 - c_2\dot{x}_1 + (c_2 + c_3)\dot{x}_2 - k_2 x_1 + (k_2 + k_3)x_2 = 0\]
These two equations take the same compact matrix form as the single-DoF system from the previous article:
\[\mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} + \mathbf{K}\mathbf{x} = \mathbf{F}\]
where the mass, damping, and stiffness matrices are
\[\mathbf{M} = \begin{bmatrix} m_1 & 0 \\ 0 & m_2 \end{bmatrix}, \quad \mathbf{C} = \begin{bmatrix} c_1+c_2 & -c_2 \\ -c_2 & c_2+c_3 \end{bmatrix}, \quad \mathbf{K} = \begin{bmatrix} k_1+k_2 & -k_2 \\ -k_2 & k_2+k_3 \end{bmatrix}\]
Note
The off-diagonal terms in \(\mathbf{C}\) and \(\mathbf{K}\) are negative and equal to the shared element values. They encode the coupling: a positive displacement of mass 1 produces a restoring force on mass 1 and an attractive force on mass 2.
In Python, the matrices are constructed directly from the system parameters:
Python
m1 = 1.0
m2 = 2.0
k1 = 100.0
k2 = 50.0
k3 = 100.0
c1 = 2.0
c2 = 1.0
c3 = 2.0
M = np.array([
[m1, 0 ],
[0, m2]
])
C = np.array([
[c1 + c2, -c2 ],
[-c2, c2 + c3 ]
])
K = np.array([
[k1 + k2, -k2 ],
[-k2, k2 + k3 ]
])
Simulating Perturbation
The first scenario is a free response: mass 1 is pulled slightly aside, held while everything is still, then released. The initial conditions set mass 1 displaced by 0.1 units while everything else starts from rest:
Python
y0 = np.array([0.1, 0.0, 0.0, 0.0])
The state vector \(\mathbf{y} = [x_1,\, x_2,\, \dot{x}_1,\, \dot{x}_2]\) is integrated over 10 seconds using SciPy's solve_ivp:
Python
t_eval = np.linspace(0, 10, 1000)
sol = solve_ivp(
fun=lambda t, y: eom(t, y, M, C, K),
t_span=(0, 10),
y0=y0,
t_eval=t_eval
)
The equation-of-motion function rearranges the matrix form to isolate the acceleration vector:
Python
import numpy as np
from scipy.integrate import solve_ivp
def eom(t, y, M, C, K):
x = y[0:2]
v = y[2:4]
f = np.array([0.0, 0.0]) # no external forcing
a = np.linalg.solve(M, f - C @ v - K @ x)
return np.hstack((v, a))
The key line evaluates
\[\ddot{\mathbf{x}} = \mathbf{M}^{-1}(\mathbf{F} - \mathbf{C}\dot{\mathbf{x}} - \mathbf{K}\mathbf{x})\]
by solving the linear system \(\mathbf{M}\ddot{\mathbf{x}} = \mathbf{F} - \mathbf{C}\dot{\mathbf{x}} - \mathbf{K}\mathbf{x}\) directly — passing \(\mathbf{M}\) separately avoids the numerical cost of explicitly inverting it.
Displacement of both masses over time with mass 1 starting at 0.1 and mass 2 at rest. The coupling transfers energy between the masses while damping dissipates it until both settle back to equilibrium.
The graph confirms the initial conditions: mass 1 begins at 0.1 and mass 2 at zero, both with zero velocity — shown by the horizontal tangent of each curve at t = 0. The shared spring and damper couple the two masses so that energy displaced into mass 1 is gradually transferred to mass 2. Both masses oscillate at a blend of the two natural frequencies of the system, and damping steadily drains the energy until both return to rest.
Forced Vibration
When a steady oscillating force is applied to the system, it never fully settles — even with damping present, the input keeps feeding energy in. To isolate this behaviour the system now starts from rest, with zero displacement and zero velocity:
Python
y0 = np.array([0.0, 0.0, 0.0, 0.0])
A sinusoidal force is applied to mass 1 alone by replacing the zero forcing vector inside eom with a time-dependent one:
Python
f = np.array([10.0 * np.sin(5.0 * t), 0.0])
This drives mass 1 with an amplitude of 10 at an angular frequency of 5 rad/s — about \(5/2\pi \approx 0.80\) Hz. Note that the forcing now depends on \(t\), so the time argument already passed to eom is finally put to use.
Displacement of both masses with a 5 rad/s (≈0.80 Hz) sinusoidal force applied to mass 1, starting from rest.
The response shows mass 1 reacting to the force first while mass 2 lags behind, dragged along through the shared spring and damper. After the initial transient dies away, both masses settle into a steady-state oscillation at the driving frequency, moving together in the same direction. The natural question is what happens as the driving frequency changes — and that is where the natural frequencies provide great insight.
Mathematical Vibration Modeling
As explained in the article on the complex exponential, \(e^{st}\) is the natural ansatz for an equation like this: differentiating it only ever multiplies it by \(s\), so substituting \(\mathbf{x} = \boldsymbol{\phi}\,e^{st}\) turns a differential equation in \(t\) into plain algebra in \(s\). That article works out why this holds on a mathematical level; the demo below shows what the same building block actually does in motion.
The solid curve is the real part of \(e^{st}\), \(e^{\sigma t}\cos\omega t\); the dashed curves are its envelope \(\pm e^{\sigma t}\). Drag \(\sigma\) negative and the curve decays like a damped vibration dying out; push it positive and the motion runs away; set it to zero and \(\omega\) alone produces a pure, undying oscillation.
Natural Frequencies and Mode Shapes
Every undamped SMD system has a set of special frequencies at which it naturally prefers to vibrate, each with an associated shape of motion. To find them, set damping and forcing aside and start from the equation of motion:
\[\mathbf{M}\ddot{\mathbf{x}} + \mathbf{K}\mathbf{x} = \mathbf{0}\]
As covered in the article on eigenvalues and eigenvectors, the eigen equation
\[\mathbf{A}\mathbf{v} = \lambda\mathbf{v}\]
singles out the directions a matrix leaves unrotated, scaling each by its eigenvalue \(\lambda\).
With no damping or forcing, the \(s\) in \(\mathbf{x} = \boldsymbol{\phi}\,e^{st}\) must be purely imaginary, \(s = i\omega\), so \(\mathbf{x} = \boldsymbol{\phi}\,e^{i\omega t}\). Differentiating twice brings down a factor of \((i\omega)^2 = -\omega^2\), so \(\ddot{\mathbf{x}} = -\omega^2\boldsymbol{\phi}\,e^{i\omega t}\), and substituting both into the equation of motion gives
\[-\omega^2\mathbf{M}\boldsymbol{\phi}\,e^{i\omega t} + \mathbf{K}\boldsymbol{\phi}\,e^{i\omega t} = \mathbf{0}\]
The exponential is never zero, so it cancels from both terms, leaving the generalized eigenvalue problem on the mass and stiffness matrices:
\[\mathbf{K}\boldsymbol{\phi} = \omega^2\,\mathbf{M}\boldsymbol{\phi}\]
Here \(\omega\) is a natural frequency and \(\boldsymbol{\phi}\) its mode shape — the relative amplitudes of the masses when the system oscillates purely at that frequency. SciPy solves this directly with eigh, which is built for symmetric matrix pairs like \(\mathbf{K}\) and \(\mathbf{M}\):
Python
from scipy.linalg import eigh
import numpy as np
omega_squared, Phi = eigh(K, M)
omega = np.sqrt(omega_squared)
freq_hz = omega / (2 * np.pi)
print("Natural frequencies (rad/s):", omega)
print("Natural frequencies (Hz):", freq_hz)
print("Mode shapes (columns):\n", Phi)
For the parameter values used above, the two natural frequencies are
and the corresponding mode shapes, normalized so that mass 1 has unit amplitude, are
Because the system has two degrees of freedom, the eigenvalue problem returns exactly two solutions — two natural frequencies, each paired with its own mode shape. In general an N-DoF system has N natural frequencies and N normal modes, one pair per degree of freedom.
The first mode is in-phase — both masses move in the same direction, with mass 2 swinging further. The second mode is out-of-phase — the masses move in opposite directions. Any free motion of the system, like the perturbation simulated earlier, is just a blend of these two modes.
Resonance at the Natural Frequencies
The natural frequencies are more than abstract eigenvalues — they are the frequencies at which the system is most responsive. Driving mass 1 with a sinusoidal force tuned to a natural frequency feeds energy in step with the system's preferred motion, so the amplitude builds over several cycles before the damping caps it at a large steady-state swing. The two figures below apply the same forcing as before, but now tuned to each natural frequency in turn.
Response to a sinusoidal force on mass 1 at the first natural frequency (ω₁ ≈ ), starting from rest. The masses swing in phase, tracing the first mode shape, and the amplitude grows well beyond the off-resonance forced response.
Response to a sinusoidal force on mass 1 at the second natural frequency (ω₂ ≈ ), starting from rest. The masses move in opposition, tracing the second mode shape.
Each driving frequency excites its own mode: at \(\omega_1\) the masses move together, and at \(\omega_2\) they move in opposition — exactly the mode shapes found above. This selective amplification is why natural frequencies matter so much in practice, where an unintended resonance can drive a structure to failure.
Sweeping the driving frequency across a whole range and recording the steady-state amplitude of each mass traces out the system's frequency response function (FRF). For a harmonic force \(\mathbf{F}e^{i\omega t}\) the steady-state amplitude follows directly from the dynamic stiffness matrix:
\[\mathbf{X}(\omega) = \left(\mathbf{K} - \omega^2\mathbf{M} + i\omega\mathbf{C}\right)^{-1}\mathbf{F}\]
Plotting \(|\mathbf{X}(\omega)|\) for each mass reveals two sharp peaks — one at every natural frequency — separated by a dip where the masses barely respond at all.
Frequency response function: steady-state amplitude of each mass against driving frequency, with the two natural frequencies marked. Both masses peak at the natural frequencies and fall away in between.
The amplitude tells only half the story. The phase records how far each mass lags the driving force: near zero at low frequencies, then dropping by roughly 180° as the sweep passes through each resonance — the signature of a mass crossing over from moving with the force to moving against it.
Phase of each mass's steady-state response relative to the driving force. Every resonance adds about 180° of lag, so mass 2 — which responds strongly to both modes — accumulates close to 360° across the sweep.