The previous post introduced 2-DOF systems—adding one mass doubled our natural frequencies and revealed the concept of mode shapes. But real-world structures rarely stop at two degrees of freedom.

Consider: a finite element model of an aircraft wing might have millions of degrees of freedom. An automotive suspension system, thousands. A building frame, hundreds. Writing individual equations for each mass quickly becomes intractable—we need the power of matrix methods.

This post, based on Chapter 6 of Rao’s Mechanical Vibrations, generalizes our analysis to Multi-Degree-of-Freedom (MDOF) Systems. We will cover:

• Matrix formulation of the equations of motion
• Three approaches for deriving system matrices (Newton, influence coefficients, Lagrange)
• The eigenvalue problem that yields natural frequencies and mode shapes
• Modal analysis—the elegant technique that decouples $n$ coupled equations into $n$ independent oscillators

The General Equation of Motion

When analyzing systems with $n$ masses, writing out $n$ separate scalar equations becomes impractical. For a 100-DOF system, we would need 100 coupled differential equations! Instead, we organize variables into vectors and system properties into matrices.

The equation of motion for an $n$-DOF damped system is:

$$ [\mathbf{M}]\ddot{\mathbf{x}} + [\mathbf{C}]\dot{\mathbf{x}} + [\mathbf{K}]\mathbf{x} = \mathbf{F} $$

This single matrix equation encapsulates all $n$ equations of motion:

• Inertial forces: $[\mathbf{M}]\ddot{\mathbf{x}}$ (mass matrix × acceleration)
• Damping forces: $[\mathbf{C}]\dot{\mathbf{x}}$ (velocity-dependent dissipation)
• Restoring forces: $[\mathbf{K}]\mathbf{x}$ (elastic spring forces)
• External forces: $\mathbf{F}$ (applied loads)

The positive definiteness of $[\mathbf{M}]$ ensures all natural frequencies are real. The symmetry of $[\mathbf{K}]$ is guaranteed by Maxwell’s reciprocity theorem. For damping, Rayleigh (proportional) damping $[\mathbf{C}] = \alpha[\mathbf{M}] + \beta[\mathbf{K}]$ is commonly assumed, as it preserves the structure of mode shapes.

Building the System Matrices: Three Approaches

The matrix equation is elegant, but where do the actual numbers come from? How do we populate $[\mathbf{M}]$, $[\mathbf{C}]$, and $[\mathbf{K}]$ for a real system? Rao presents three complementary methods, each suited to different problem types:

Newton’s Second Law (Direct Method)

The most intuitive approach: draw a free-body diagram for each mass, identify all forces (spring forces, damping forces, external forces), and apply $\Sigma F = ma$.

For a chain of masses connected by springs, this gives equations like:

$$ m_i \ddot{x}_i = -k_i(x_i - x_{i-1}) - k_{i+1}(x_i - x_{i+1}) + F_i $$

Pros: Intuitive, directly visualizable
Cons: Tedious for complex systems; difficult when rotational DOFs or constraints are involved

Influence Coefficients (Structural Method)

This method is particularly powerful for continuous structures (beams, plates) where we discretize at specific points.

The matrices are inverses: $[\mathbf{K}] = [\mathbf{A}]^{-1}$. Maxwell’s reciprocity theorem guarantees symmetry: $k_{ij} = k_{ji}$, meaning the force at $i$ due to displacement at $j$ equals the force at $j$ due to the same displacement at $i$.

Lagrange’s Equations (Energy Method)

For complex systems, energy methods avoid the vector complications of force analysis. We define:

• $T$ = Kinetic energy: $T = \frac{1}{2}\dot{\mathbf{x}}^T[\mathbf{M}]\dot{\mathbf{x}}$
• $V$ = Potential energy: $V = \frac{1}{2}\mathbf{x}^T[\mathbf{K}]\mathbf{x}$
• $R$ = Rayleigh dissipation: $R = \frac{1}{2}\dot{\mathbf{x}}^T[\mathbf{C}]\dot{\mathbf{x}}$

Then apply: $$ \frac{d}{dt}\left(\frac{\partial T}{\partial \dot{x}_i}\right) - \frac{\partial T}{\partial x_i} + \frac{\partial R}{\partial \dot{x}_i} + \frac{\partial V}{\partial x_i} = F_i $$

The Eigenvalue Problem: Natural Frequencies and Mode Shapes

With the system matrices in hand, we arrive at the central question of vibration analysis: At what frequencies will this structure naturally vibrate, and what patterns will it assume?

The answer lies in the eigenvalue problem—a cornerstone of linear algebra that maps directly onto the physics of oscillating systems.

From Differential Equation to Algebraic Problem

For free vibration (no damping, no external force), the equation of motion simplifies to:

$$ [\mathbf{M}]\ddot{\mathbf{x}} + [\mathbf{K}]\mathbf{x} = \mathbf{0} $$

We seek solutions where all parts of the system oscillate at the same frequency—synchronous harmonic motion:

$$ \mathbf{x}(t) = \mathbf{X} \cos(\omega t) $$

where $\mathbf{X}$ is the amplitude vector (shape) and $\omega$ is the circular frequency. Substituting this assumed solution (noting that $\ddot{\mathbf{x}} = -\omega^2 \mathbf{X} \cos(\omega t)$) yields the Generalized Eigenvalue Problem:

$$ \boxed{\left[ [\mathbf{K}] - \omega^2[\mathbf{M}] \right] \mathbf{X} = \mathbf{0}} $$

The Characteristic Equation

For a non-trivial solution (where $\mathbf{X} \neq \mathbf{0}$, meaning the system actually moves), the coefficient matrix must be singular. This requires:

$$ \det \left| [\mathbf{K}] - \omega^2[\mathbf{M}] \right| = 0 $$

Expanding this determinant produces the Characteristic Equation—an $n$-th degree polynomial in $\omega^2$:

$$ a_n(\omega^2)^n + a_{n-1}(\omega^2)^{n-1} + \cdots + a_1(\omega^2) + a_0 = 0 $$

Eigenvalues and Eigenvectors

The $n$ roots of the characteristic equation give us $n$ natural frequencies: $\omega_1 \leq \omega_2 \leq \cdots \leq \omega_n$. For each frequency $\omega_i$, back-substituting into the eigenvalue equation yields a corresponding eigenvector $\mathbf{X}^{(i)}$—the mode shape describing how the masses move relative to each other at that frequency.

Modal Analysis: Decoupling Through Coordinate Transformation

We now have $n$ natural frequencies and $n$ mode shapes—but the original equations of motion remain coupled. Solving $n$ simultaneous differential equations is computationally expensive and obscures physical understanding.

Modal analysis provides an elegant solution: by transforming to a special coordinate system aligned with the mode shapes, the coupled system becomes $n$ independent single-degree-of-freedom oscillators—each solvable with the techniques we mastered in Parts 2–4.

The Key Insight: Orthogonality of Mode Shapes

The mathematical foundation of modal analysis rests on a remarkable property: mode shapes are orthogonal with respect to both the mass and stiffness matrices. For any two distinct modes $i$ and $j$ (where $i \neq j$):

$$ \mathbf{X}^{(i)T}[\mathbf{M}]\mathbf{X}^{(j)} = 0 \quad \text{(mass-orthogonality)} $$

$$ \mathbf{X}^{(i)T}[\mathbf{K}]\mathbf{X}^{(j)} = 0 \quad \text{(stiffness-orthogonality)} $$

This means different modes are energetically independent—they don’t exchange energy with each other during vibration!

The Modal Transformation

We exploit orthogonality by constructing the Modal Matrix $[\mathbf{\Phi}]$, whose columns are the eigenvectors (mode shapes):

$$ [\mathbf{\Phi}] = \begin{bmatrix} \mathbf{X}^{(1)} & \mathbf{X}^{(2)} & \cdots & \mathbf{X}^{(n)} \end{bmatrix} $$

Introducing the principal (modal) coordinates $\mathbf{q}(t)$ through the transformation:

$$ \mathbf{x}(t) = [\mathbf{\Phi}]\mathbf{q}(t) $$

Substituting into the equation of motion and pre-multiplying by $[\mathbf{\Phi}]^T$ yields:

$$ [\mathbf{\Phi}]^T[\mathbf{M}][\mathbf{\Phi}] \ddot{\mathbf{q}} + [\mathbf{\Phi}]^T[\mathbf{K}][\mathbf{\Phi}] \mathbf{q} = [\mathbf{\Phi}]^T\mathbf{F} $$

Thanks to the orthogonality properties, the matrix products $[\mathbf{\Phi}]^T[\mathbf{M}][\mathbf{\Phi}]$ and $[\mathbf{\Phi}]^T[\mathbf{K}][\mathbf{\Phi}]$ become diagonal matrices (all off-diagonal terms are zero). This decouples the system completely into:

$$ [\tilde{\mathbf{M}}] \ddot{\mathbf{q}} + [\tilde{\mathbf{K}}] \mathbf{q} = \tilde{\mathbf{F}} $$

The Decoupled Equations

Each modal coordinate $q_i(t)$ satisfies an independent SDOF equation:

$$ \tilde{m}_i \ddot{q}_i + \tilde{k}_i q_i = \tilde{Q}_i(t) $$

where the modal mass, modal stiffness, and modal force are:

The natural frequency of each mode remains: $\omega_i = \sqrt{\tilde{k}_i / \tilde{m}_i}$

Computational Example: Eigenvalue Solver in Julia

Let’s put theory into practice. Julia’s LinearAlgebra standard library provides efficient routines for the generalized eigenvalue problem—perfect for vibration analysis.

We’ll solve for the natural frequencies and mode shapes of a 3-DOF spring-mass chain (based on Example 6.11 in Rao). The system consists of three equal masses connected by springs to a fixed wall:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44

using LinearAlgebra

# System parameters
m = 1.0  # Mass (kg)
k = 1.0  # Stiffness (N/m)

# Define Mass and Stiffness Matrices
M = m * I(3)  # Diagonal mass matrix

K = k * [2  -1   0;
        -1   2  -1;
         0  -1   1]  # Tridiagonal stiffness matrix

# Solve the Generalized Eigenvalue Problem: [K]X = λ[M]X
eigenresult = eigen(K, M)
eigenvalues = eigenresult.values      # λᵢ = ωᵢ²
eigenvectors = eigenresult.vectors    # Mode shapes

# Extract and sort natural frequencies
ω = sqrt.(real.(eigenvalues))
sorted_idx = sortperm(ω)
ω = ω[sorted_idx]
Φ = eigenvectors[:, sorted_idx]

# Mass-normalize mode shapes: Xᵀ[M]X = 1
for i in 1:size(Φ, 2)
    modal_mass = Φ[:, i]' * M * Φ[:, i]
    Φ[:, i] ./= sqrt(modal_mass)
end

# Display results
println("Natural Frequencies (rad/s):")
println("  ω = ", round.(ω, digits=4))

println("\nMode Shape Matrix [Φ] (columns = modes):")
display(round.(Φ, digits=4))

# Verify orthogonality
println("\nOrthogonality Check:")
println("Φᵀ[M]Φ (should be Identity):")
display(round.(Φ' * M * Φ, digits=4))

println("\nΦᵀ[K]Φ (should be diag(ω²)):")
display(round.(Φ' * K * Φ, digits=4))

Expected Output:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

Natural Frequencies (rad/s):
  ω = [0.445, 1.247, 1.8019]

Mode Shape Matrix [Φ] (columns = modes):
3×3 Matrix{Float64}:
 0.328  -0.737  -0.591
 0.591  -0.328   0.737
 0.737   0.591  -0.328

Orthogonality Check:
Φᵀ[M]Φ (should be Identity):
3×3 Matrix{Float64}:
  1.0  -0.0  0.0
 -0.0   1.0  0.0
  0.0   0.0  1.0

Φᵀ[K]Φ (should be diag(ω²)):
3×3 Matrix{Float64}:
  0.1981  -0.0     0.0
 -0.0      1.555  -0.0
  0.0      0.0     3.247

Interpreting the Results:

• Natural Frequencies: The three values $\omega_1 < \omega_2 < \omega_3$ represent the resonant speeds of the system
• Mode Shapes: Each column of $[\Phi]$ shows how masses move relative to each other:
  • Mode 1 ($\omega_1 = 0.445$ rad/s): All masses move in phase (fundamental mode)
  • Mode 2 ($\omega_2 = 1.247$ rad/s): First mass opposes the others
  • Mode 3 ($\omega_3 = 1.802$ rad/s): Alternating pattern with two nodal regions
• Orthogonality Verification: $\Phi^T M \Phi = I$ confirms mass-normalization; $\Phi^T K \Phi = \text{diag}(\omega_i^2)$ confirms stiffness decoupling

Summary

This post marked our transition from simple oscillators to the real-world complexity of Multi-Degree-of-Freedom systems. The key concepts we covered form the foundation of modern structural dynamics:

The elegance of modal analysis lies in its physical interpretation: every complex vibration is simply a superposition of independent modes, each oscillating at its own natural frequency. This insight—combined with modal truncation—makes it possible to analyze structures with thousands of degrees of freedom by focusing on just the first few dominant modes.

What’s Next?

For massive systems like a 100-story building or an aircraft fuselage, even setting up the eigenvalue problem can be computationally prohibitive. In Part 7, we will explore Numerical Methods for Determination of Natural Frequencies—including Dunkerley’s Formula, Rayleigh’s Method, and Holzer’s Method—that allow engineers to estimate frequencies without explicitly solving the full characteristic polynomial.

References

Rao, S. S. (2018). Mechanical Vibrations (6th ed.). Pearson.