In an undamped multi-degree-of-freedom system, the generalized displacement vector \( \mathbf{q}(t)\in\mathbb{R}^n \) satisfies a second-order linear ODE of the form

\[ \mathbf{M}\,\ddot{\mathbf{q}}(t) + \mathbf{K}\,\mathbf{q}(t)=\mathbf{0}, \quad \mathbf{M}\in\mathbb{R}^{n\times n},\; \mathbf{K}\in\mathbb{R}^{n\times n}. \]

The matrices \( \mathbf{M} \) (mass/inertia) and \( \mathbf{K} \) (stiffness) are typically symmetric. For a stable structure without rigid-body motions, one usually has: \( \mathbf{M}\succ 0 \) and \( \mathbf{K}\succ 0 \). (If rigid-body modes exist, \( \mathbf{K}\succeq 0 \) and some natural frequencies become zero.)

A normal mode is a special motion where every coordinate oscillates sinusoidally at the same frequency, with a fixed shape (relative amplitudes). We seek solutions of the form \( \mathbf{q}(t)=\boldsymbol{\phi}\cos(\omega t) \), leading to an eigenvalue problem for \( \omega \) and \( \boldsymbol{\phi} \).

flowchart TD
  A["Start: M qdd + K q = 0"] --> B["Assume harmonic: q(t) = phi * cos(w t)"]
  B --> C["Substitute -> (K - w^2 M) phi = 0"]
  C --> D["Nontrivial phi requires det(K - w^2 M) = 0"]
  D --> E["Solve for w_i^2 (eigenvalues)"]
  E --> F["Compute phi_i (mode shapes)"]
  F --> G["Normalize modes (e.g., phi^T M phi = 1)"]
  G --> H["Check orthogonality: phi_i^T M phi_j = 0 for i != j"]
        

Consider the undamped free-vibration equation: \( \mathbf{M}\ddot{\mathbf{q}}(t)+\mathbf{K}\mathbf{q}(t)=\mathbf{0} \). Assume a single-frequency trial solution:

\[ \mathbf{q}(t)=\boldsymbol{\phi}\cos(\omega t), \quad \ddot{\mathbf{q}}(t)=-\omega^2\boldsymbol{\phi}\cos(\omega t). \]

Substitution yields

\[ \left(\mathbf{K}-\omega^2\mathbf{M}\right)\boldsymbol{\phi}\cos(\omega t)=\mathbf{0} \;\;\Longrightarrow\;\; \left(\mathbf{K}-\omega^2\mathbf{M}\right)\boldsymbol{\phi}=\mathbf{0}. \]

A nontrivial mode shape \( \boldsymbol{\phi}\neq\mathbf{0} \) exists iff the matrix is singular:

\[ \det\!\left(\mathbf{K}-\omega^2\mathbf{M}\right)=0. \]

Let \( \lambda=\omega^2 \). Then we obtain the symmetric generalized eigenproblem:

\[ \mathbf{K}\boldsymbol{\phi}=\lambda\,\mathbf{M}\boldsymbol{\phi}. \]

The set \( \{\lambda_i,\boldsymbol{\phi}_i\}_{i=1}^n \) provides the squared natural frequencies \( \omega_i=\sqrt{\lambda_i} \) and corresponding normal modes.

Suppose \( \mathbf{M}\succ 0 \) and \( \mathbf{K}\succeq 0 \) are symmetric. For any nonzero vector \( \mathbf{x}\neq\mathbf{0} \), define the Rayleigh quotient:

\[ \mathcal{R}(\mathbf{x}) \;=\; \frac{\mathbf{x}^\mathsf{T}\mathbf{K}\mathbf{x}}{\mathbf{x}^\mathsf{T}\mathbf{M}\mathbf{x}}. \]

Because \( \mathbf{x}^\mathsf{T}\mathbf{M}\mathbf{x} > 0 \) for \( \mathbf{x}\neq\mathbf{0} \) and \( \mathbf{x}^\mathsf{T}\mathbf{K}\mathbf{x} \ge 0 \) when \( \mathbf{K}\succeq 0 \), we have \( \mathcal{R}(\mathbf{x}) \ge 0 \).

Key identity (eigenpairs). If \( \mathbf{K}\boldsymbol{\phi}=\lambda\mathbf{M}\boldsymbol{\phi} \) with \( \boldsymbol{\phi}\neq\mathbf{0} \), then:

\[ \lambda = \frac{\boldsymbol{\phi}^\mathsf{T}\mathbf{K}\boldsymbol{\phi}}{\boldsymbol{\phi}^\mathsf{T}\mathbf{M}\boldsymbol{\phi}} = \mathcal{R}(\boldsymbol{\phi}) \;\;\ge\; 0. \]

Proof. Left-multiply by \( \boldsymbol{\phi}^\mathsf{T} \):

\[ \boldsymbol{\phi}^\mathsf{T}\mathbf{K}\boldsymbol{\phi} = \lambda\,\boldsymbol{\phi}^\mathsf{T}\mathbf{M}\boldsymbol{\phi} \;\;\Longrightarrow\;\; \lambda=\frac{\boldsymbol{\phi}^\mathsf{T}\mathbf{K}\boldsymbol{\phi}}{\boldsymbol{\phi}^\mathsf{T}\mathbf{M}\boldsymbol{\phi}}. \]

Therefore \( \lambda=\omega^2 \) is real and nonnegative. If, additionally, \( \mathbf{K}\succ 0 \), then \( \lambda>0 \) and hence \( \omega>0 \).

Energy interpretation. For a displacement \( \mathbf{q} \), define the potential energy \( V=\tfrac{1}{2}\mathbf{q}^\mathsf{T}\mathbf{K}\mathbf{q} \) and kinetic energy \( T=\tfrac{1}{2}\dot{\mathbf{q}}^\mathsf{T}\mathbf{M}\dot{\mathbf{q}} \). For a mode \( \mathbf{q}(t)=\boldsymbol{\phi}\cos(\omega t) \), the quotient \( \boldsymbol{\phi}^\mathsf{T}\mathbf{K}\boldsymbol{\phi}/\boldsymbol{\phi}^\mathsf{T}\mathbf{M}\boldsymbol{\phi} \) is precisely \( \omega^2 \).

Let \( (\lambda_i,\boldsymbol{\phi}_i) \) and \( (\lambda_j,\boldsymbol{\phi}_j) \) be two eigenpairs of \( \mathbf{K}\boldsymbol{\phi}=\lambda\mathbf{M}\boldsymbol{\phi} \) with \( \lambda_i\neq\lambda_j \). Then the corresponding eigenvectors are orthogonal in the \( \mathbf{M} \)-inner product:

\[ \boldsymbol{\phi}_i^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_j = 0. \]

They are also orthogonal with respect to \( \mathbf{K} \):

\[ \boldsymbol{\phi}_i^\mathsf{T}\mathbf{K}\boldsymbol{\phi}_j = 0. \]

Proof (mass orthogonality). Start with

\[ \mathbf{K}\boldsymbol{\phi}_i=\lambda_i\mathbf{M}\boldsymbol{\phi}_i, \quad \mathbf{K}\boldsymbol{\phi}_j=\lambda_j\mathbf{M}\boldsymbol{\phi}_j. \]

Left-multiply the first equation by \( \boldsymbol{\phi}_j^\mathsf{T} \) and the second by \( \boldsymbol{\phi}_i^\mathsf{T} \):

\[ \boldsymbol{\phi}_j^\mathsf{T}\mathbf{K}\boldsymbol{\phi}_i=\lambda_i\,\boldsymbol{\phi}_j^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_i, \quad \boldsymbol{\phi}_i^\mathsf{T}\mathbf{K}\boldsymbol{\phi}_j=\lambda_j\,\boldsymbol{\phi}_i^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_j. \]

Since \( \mathbf{K}=\mathbf{K}^\mathsf{T} \), the left-hand sides are equal scalars: \( \boldsymbol{\phi}_j^\mathsf{T}\mathbf{K}\boldsymbol{\phi}_i=\boldsymbol{\phi}_i^\mathsf{T}\mathbf{K}\boldsymbol{\phi}_j \). Subtracting yields

\[ (\lambda_i-\lambda_j)\,\boldsymbol{\phi}_i^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_j = 0. \]

Because \( \lambda_i\neq\lambda_j \), it follows that \( \boldsymbol{\phi}_i^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_j=0 \). Stiffness orthogonality follows immediately by substituting the eigenrelation: \( \boldsymbol{\phi}_i^\mathsf{T}\mathbf{K}\boldsymbol{\phi}_j=\lambda_j\,\boldsymbol{\phi}_i^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_j=0 \).

Normalization conventions. A convenient choice is mass normalization:

\[ \boldsymbol{\phi}_i^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_i = 1 \quad \Longrightarrow \quad \boldsymbol{\phi}_i^\mathsf{T}\mathbf{K}\boldsymbol{\phi}_i=\lambda_i=\omega_i^2. \]

With the modal matrix \( \mathbf{\Phi}=[\boldsymbol{\phi}_1\;\cdots\;\boldsymbol{\phi}_n] \), mass-normalization implies

\[ \mathbf{\Phi}^\mathsf{T}\mathbf{M}\mathbf{\Phi}=\mathbf{I}, \quad \mathbf{\Phi}^\mathsf{T}\mathbf{K}\mathbf{\Phi}=\mathbf{\Omega}^2 \;\;=\;\;\operatorname{diag}(\omega_1^2,\dots,\omega_n^2). \]

Many numerical libraries solve a standard symmetric eigenproblem \( \mathbf{A}\mathbf{u}=\lambda\mathbf{u} \). For \( \mathbf{M}\succ 0 \), we can convert the generalized problem \( \mathbf{K}\boldsymbol{\phi}=\lambda\mathbf{M}\boldsymbol{\phi} \) into a standard one.

Let \( \mathbf{M}=\mathbf{L}\mathbf{L}^\mathsf{T} \) be the Cholesky factorization (lower triangular \( \mathbf{L} \) invertible). Define \( \mathbf{u}=\mathbf{L}^\mathsf{T}\boldsymbol{\phi} \). Then \( \boldsymbol{\phi}=\mathbf{L}^{-\mathsf{T}}\mathbf{u} \) and

\[ \mathbf{K}\mathbf{L}^{-\mathsf{T}}\mathbf{u} = \lambda\,\mathbf{M}\mathbf{L}^{-\mathsf{T}}\mathbf{u} = \lambda\,\mathbf{L}\mathbf{L}^\mathsf{T}\mathbf{L}^{-\mathsf{T}}\mathbf{u} = \lambda\,\mathbf{L}\mathbf{u}. \]

Left-multiply by \( \mathbf{L}^{-1} \):

\[ \underbrace{\mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-\mathsf{T}}}_{\mathbf{A}} \mathbf{u} = \lambda\,\mathbf{u}. \]

When \( \mathbf{K}=\mathbf{K}^\mathsf{T} \), the matrix \( \mathbf{A}=\mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-\mathsf{T}} \) is symmetric:

\[ \mathbf{A}^\mathsf{T} = \left(\mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-\mathsf{T}}\right)^\mathsf{T} = \mathbf{L}^{-1}\mathbf{K}^\mathsf{T}\mathbf{L}^{-\mathsf{T}} = \mathbf{A}. \]

Therefore, all \( \lambda \) are real and can be computed using standard symmetric solvers, then recover the mode shapes by \( \boldsymbol{\phi}=\mathbf{L}^{-\mathsf{T}}\mathbf{u} \).

flowchart TD
  A["Given symmetric M (SPD) and K"] --> B["Cholesky: M = L L^T"]
  B --> C["Form A = inv(L) * K * inv(L^T)"]
  C --> D["Solve A u = (w^2) u (symmetric)"]
  D --> E["Recover modes: phi = inv(L^T) u"]
  E --> F["Mass-normalize: phi^T M phi = 1"]
  F --> G["Verify: Phi^T M Phi = I and Phi^T K Phi = diag(w^2)"]
        

• Python: numpy, scipy.linalg.eigh (dense symmetric generalized), scipy.sparse.linalg.eigsh (large sparse), plus matplotlib for plotting.
• C++: Eigen (generalized self-adjoint eigensolver), and optionally Spectra for large sparse eigenproblems.
• Java: EJML for Cholesky + symmetric eigendecomposition (via the reduction in Section 5).
• MATLAB: built-in eig(K,M) for generalized eigenpairs; state-space analysis via eig(A).
• Simulink: implement \( \mathbf{M}\ddot{\mathbf{q}}+\mathbf{K}\mathbf{q}=0 \) using a State-Space block (with \( \mathbf{x}=[\mathbf{q};\dot{\mathbf{q}}] \)) or integrator blocks with matrix gains; compute frequencies from \( \mathbf{A} \) eigenvalues for the undamped case.
• Wolfram Mathematica: Eigensystem[{K, M}] solves the generalized problem directly.

This script uses a Cholesky reduction (Section 5) to compute modes, verifies orthogonality, and simulates free response using modal superposition.

Chapter13_Lesson1.py

"""
Chapter 13 - Vibrations and Multi-Degree-of-Freedom (MDOF) Systems
Lesson 1: Natural Frequencies and Normal Modes (Eigenvalue Problems)

This script:
  1) Builds (M, K) for a 3-DOF undamped mass-spring chain.
  2) Solves the symmetric generalized eigenproblem: K phi = (w^2) M phi.
  3) Mass-normalizes modes so that Phi^T M Phi = I and Phi^T K Phi = diag(w^2).
  4) Simulates free vibration via modal superposition for given initial conditions.
"""

import numpy as np

def mdof_chain_mk(masses, springs):
    """
    3-DOF chain (generalizable):
      wall --k1-- m1 --k2-- m2 --k3-- m3 -- (free end)
    K is assembled with end conditions matching the picture above.
    """
    n = len(masses)
    M = np.diag(masses).astype(float)

    # Springs: k1 between wall and m1; k2 between m1 and m2; ...; kn between m(n-1) and mn
    # Here springs length should be n (k1..kn)
    if len(springs) != n:
        raise ValueError("springs must have length n (k1..kn).")

    K = np.zeros((n, n), dtype=float)
    for i in range(n):
        # left spring contribution
        K[i, i] += springs[i]
        if i > 0:
            K[i, i] += springs[i-1]
            K[i, i-1] -= springs[i-1]
            K[i-1, i] -= springs[i-1]
    return M, K

def generalized_modes(M, K):
    """
    Solve K phi = (w^2) M phi for symmetric (M,K) with M SPD and K SPSD/SPD.

    Uses dense approach:
      - Convert to standard symmetric eigenproblem via Cholesky: M = L L^T
      - Solve (L^-1 K L^-T) u = (w^2) u
      - Recover phi = L^-T u
    """
    # Cholesky factor
    L = np.linalg.cholesky(M)

    # A = L^-1 K L^-T (symmetric)
    Linv = np.linalg.inv(L)
    A = Linv @ K @ Linv.T

    # Standard symmetric eigenproblem
    w2, U = np.linalg.eigh(A)  # ascending
    # Numerical cleanup for tiny negative values
    w2 = np.maximum(w2, 0.0)
    w = np.sqrt(w2)

    # Recover modes in physical coordinates
    Phi = np.linalg.solve(L.T, U)  # L^T Phi = U

    # Mass-normalize: phi_i^T M phi_i = 1
    for i in range(Phi.shape[1]):
        mi = Phi[:, i].T @ M @ Phi[:, i]
        Phi[:, i] = Phi[:, i] / np.sqrt(mi)

    return w, Phi

def check_orthogonality(M, K, w, Phi):
    Mt = Phi.T @ M @ Phi
    Kt = Phi.T @ K @ Phi
    I = np.eye(Mt.shape[0])
    D = np.diag(w**2)
    err_M = np.linalg.norm(Mt - I, ord=np.inf)
    err_K = np.linalg.norm(Kt - D, ord=np.inf)
    return err_M, err_K, Mt, Kt

def modal_free_response(M, K, w, Phi, q0, v0, t):
    """
    Undamped free vibration:
      q(t) = Phi * eta(t), with Phi^T M Phi = I
      eta_i(t) = eta0_i cos(w_i t) + (etaDot0_i / w_i) sin(w_i t)
    where eta0 = Phi^T M q0 and etaDot0 = Phi^T M v0.
    """
    q0 = np.asarray(q0, dtype=float).reshape(-1)
    v0 = np.asarray(v0, dtype=float).reshape(-1)
    eta0 = Phi.T @ M @ q0
    etaDot0 = Phi.T @ M @ v0

    eta = np.zeros((len(w), len(t)))
    for i, wi in enumerate(w):
        if wi < 1e-12:
            eta[i, :] = eta0[i] + etaDot0[i] * t
        else:
            eta[i, :] = eta0[i]*np.cos(wi*t) + (etaDot0[i]/wi)*np.sin(wi*t)

    q = Phi @ eta
    return q, eta

def main():
    # Example parameters
    masses  = [2.0, 1.5, 1.0]          # kg
    springs = [200.0, 300.0, 250.0]    # N/m

    M, K = mdof_chain_mk(masses, springs)
    w, Phi = generalized_modes(M, K)

    print("Natural frequencies (rad/s):")
    for i, wi in enumerate(w, start=1):
        print(f"  w{i} = {wi:.6f}")

    err_M, err_K, Mt, Kt = check_orthogonality(M, K, w, Phi)
    print("\nPhi^T M Phi:\n", Mt)
    print("\nPhi^T K Phi:\n", Kt)
    print(f"\nErrors (inf-norm): err_M={err_M:.3e}, err_K={err_K:.3e}")

    # Free response example
    t = np.linspace(0.0, 5.0, 1501)
    q0 = [0.02, 0.0, -0.01]   # meters
    v0 = [0.0, 0.0, 0.0]      # m/s
    q, eta = modal_free_response(M, K, w, Phi, q0, v0, t)

    # Plot (optional)
    try:
        import matplotlib.pyplot as plt
        plt.figure()
        for i in range(q.shape[0]):
            plt.plot(t, q[i, :], label=f"q{i+1}(t)")
        plt.xlabel("t (s)")
        plt.ylabel("displacement (m)")
        plt.title("3-DOF undamped free vibration (modal superposition)")
        plt.legend()
        plt.grid(True)
        plt.show()
    except Exception as e:
        print("Plotting skipped:", e)

if __name__ == "__main__":
    main()
      

For symmetric \( \mathbf{K} \) and SPD \( \mathbf{M} \), Eigen provides a direct solver GeneralizedSelfAdjointEigenSolver that computes \( \omega^2 \) and mode shapes.

Chapter13_Lesson1.cpp

// Chapter 13 - Vibrations and Multi-Degree-of-Freedom (MDOF) Systems
// Lesson 1: Natural Frequencies and Normal Modes (Eigenvalue Problems)
//
// Build:
//   - Requires Eigen (header-only): https://eigen.tuxfamily.org/
//   - Example (g++):
//       g++ -O2 -std=c++17 Chapter13_Lesson1.cpp -I path/to/eigen -o mdof_modes
//
// This program:
//   1) Builds (M, K) for a 3-DOF undamped mass-spring chain.
//   2) Solves K phi = (w^2) M phi using Eigen's GeneralizedSelfAdjointEigenSolver.
//   3) Mass-normalizes modes and verifies Phi^T M Phi = I and Phi^T K Phi = diag(w^2).

#include <iostream>
#include <vector>
#include <cmath>
#include <Eigen/Dense>

static void mdofChainMK(const std::vector<double>& masses,
                        const std::vector<double>& springs,
                        Eigen::MatrixXd& M,
                        Eigen::MatrixXd& K)
{
    const int n = static_cast<int>(masses.size());
    if (static_cast<int>(springs.size()) != n) {
        throw std::runtime_error("springs must have length n (k1..kn).");
    }

    M = Eigen::MatrixXd::Zero(n, n);
    K = Eigen::MatrixXd::Zero(n, n);

    for (int i = 0; i < n; ++i) M(i, i) = masses[i];

    for (int i = 0; i < n; ++i) {
        K(i, i) += springs[i];
        if (i > 0) {
            K(i, i)     += springs[i-1];
            K(i, i-1)   -= springs[i-1];
            K(i-1, i)   -= springs[i-1];
        }
    }
}

int main()
{
    try {
        std::vector<double> masses  = {2.0, 1.5, 1.0};
        std::vector<double> springs = {200.0, 300.0, 250.0};

        Eigen::MatrixXd M, K;
        mdofChainMK(masses, springs, M, K);

        Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd> solver(K, M);
        if (solver.info() != Eigen::Success) {
            std::cerr << "Eigen solve failed.\n";
            return 1;
        }

        Eigen::VectorXd w2 = solver.eigenvalues();    // ascending (w^2)
        Eigen::MatrixXd Phi = solver.eigenvectors();  // columns are modes

        // Clean small negative numerical noise
        for (int i = 0; i < w2.size(); ++i) if (w2(i) < 0.0 && std::abs(w2(i)) < 1e-12) w2(i) = 0.0;

        Eigen::VectorXd w(w2.size());
        for (int i = 0; i < w2.size(); ++i) w(i) = std::sqrt(std::max(0.0, w2(i)));

        // Mass-normalize columns: phi_i^T M phi_i = 1
        for (int i = 0; i < Phi.cols(); ++i) {
            double mi = Phi.col(i).transpose() * M * Phi.col(i);
            Phi.col(i) /= std::sqrt(mi);
        }

        std::cout << "Natural frequencies (rad/s):\n";
        for (int i = 0; i < w.size(); ++i) {
            std::cout << "  w" << (i+1) << " = " << w(i) << "\n";
        }

        Eigen::MatrixXd Mt = Phi.transpose() * M * Phi;
        Eigen::MatrixXd Kt = Phi.transpose() * K * Phi;

        std::cout << "\nPhi^T M Phi (should be I):\n" << Mt << "\n";
        std::cout << "\nPhi^T K Phi (should be diag(w^2)):\n" << Kt << "\n";

        // Infinity norm errors
        Eigen::MatrixXd I = Eigen::MatrixXd::Identity(Mt.rows(), Mt.cols());
        Eigen::MatrixXd D = w2.asDiagonal();
        double errM = (Mt - I).cwiseAbs().maxCoeff();
        double errK = (Kt - D).cwiseAbs().maxCoeff();

        std::cout << "\nErrors (max abs entry): err_M=" << errM << ", err_K=" << errK << "\n";
    }
    catch (const std::exception& e) {
        std::cerr << "ERROR: " << e.what() << "\n";
        return 1;
    }
    return 0;
}
      

EJML does not always expose a one-call symmetric generalized solver for \( (\mathbf{K},\mathbf{M}) \), but the reduction in Section 5 is numerically standard and keeps symmetry.

Chapter13_Lesson1.java

/*
Chapter 13 - Vibrations and Multi-Degree-of-Freedom (MDOF) Systems
Lesson 1: Natural Frequencies and Normal Modes (Eigenvalue Problems)

This example uses EJML (Efficient Java Matrix Library) to solve the generalized symmetric eigenproblem:
    K phi = (w^2) M phi
by converting it to a standard symmetric eigenproblem via Cholesky:
    M = L L^T
    A = L^{-1} K L^{-T}  (symmetric)
    A u = (w^2) u
    phi = L^{-T} u

Dependency (Gradle):
    implementation 'org.ejml:ejml-all:0.43'
*/

import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;
import org.ejml.dense.row.decomposition.chol.CholeskyDecompositionInner_DDRM;
import org.ejml.interfaces.decomposition.EigenDecomposition_F64;
import org.ejml.dense.row.factory.DecompositionFactory_DDRM;

import java.util.Arrays;

public class Chapter13_Lesson1 {

    static DMatrixRMaj diag(double[] d) {
        int n = d.length;
        DMatrixRMaj M = new DMatrixRMaj(n, n);
        for (int i = 0; i < n; i++) M.set(i, i, d[i]);
        return M;
    }

    // Build 3-DOF chain: wall-k1-m1-k2-m2-k3-m3-(free end)
    static DMatrixRMaj buildK(double[] springs) {
        int n = springs.length;
        DMatrixRMaj K = new DMatrixRMaj(n, n);

        for (int i = 0; i < n; i++) {
            K.add(i, i, springs[i]);
            if (i > 0) {
                K.add(i, i, springs[i-1]);
                K.add(i, i-1, -springs[i-1]);
                K.add(i-1, i, -springs[i-1]);
            }
        }
        return K;
    }

    static class Modes {
        double[] w;       // natural frequencies (rad/s)
        DMatrixRMaj Phi;  // columns are mass-normalized modes
        double[] w2;
    }

    static Modes generalizedModes(DMatrixRMaj M, DMatrixRMaj K) {
        int n = M.numRows;

        // Cholesky: M = L L^T (lower)
        CholeskyDecompositionInner_DDRM chol = new CholeskyDecompositionInner_DDRM(true);
        if (!chol.decompose(M.copy())) {
            throw new RuntimeException("Cholesky failed: M not SPD?");
        }
        DMatrixRMaj L = chol.getT(null); // lower-triangular

        // Compute Linv = inv(L)
        DMatrixRMaj Linv = new DMatrixRMaj(n, n);
        CommonOps_DDRM.invert(L, Linv);

        // A = Linv * K * Linv^T
        DMatrixRMaj temp = new DMatrixRMaj(n, n);
        DMatrixRMaj A = new DMatrixRMaj(n, n);
        CommonOps_DDRM.mult(Linv, K, temp);
        CommonOps_DDRM.multTransB(temp, Linv, A); // temp * Linv^T

        // Symmetric eigen decomposition
        EigenDecomposition_F64<DMatrixRMaj> eig = DecompositionFactory_DDRM.eig(n, true);
        if (!eig.decompose(A)) {
            throw new RuntimeException("Eigen decomposition failed.");
        }

        // Extract eigenpairs (not guaranteed sorted)
        double[] w2 = new double[n];
        DMatrixRMaj U = new DMatrixRMaj(n, n); // eigenvectors u
        for (int i = 0; i < n; i++) {
            w2[i] = eig.getEigenvalue(i).getReal();
            DMatrixRMaj ui = eig.getEigenVector(i);
            for (int r = 0; r < n; r++) U.set(r, i, ui.get(r, 0));
        }

        // Sort by w2 ascending
        Integer[] idx = new Integer[n];
        for (int i = 0; i < n; i++) idx[i] = i;
        Arrays.sort(idx, (a, b) -> Double.compare(w2[a], w2[b]));

        double[] w2s = new double[n];
        DMatrixRMaj Us = new DMatrixRMaj(n, n);
        for (int col = 0; col < n; col++) {
            int j = idx[col];
            w2s[col] = Math.max(0.0, w2[j]);
            for (int r = 0; r < n; r++) Us.set(r, col, U.get(r, j));
        }

        // phi = L^{-T} u => solve L^T phi = u
        DMatrixRMaj Phi = new DMatrixRMaj(n, n);
        DMatrixRMaj Lt = new DMatrixRMaj(n, n);
        CommonOps_DDRM.transpose(L, Lt);
        DMatrixRMaj Ltinv = new DMatrixRMaj(n, n);
        CommonOps_DDRM.invert(Lt, Ltinv);
        CommonOps_DDRM.mult(Ltinv, Us, Phi);

        // Mass-normalize columns: phi_i^T M phi_i = 1
        DMatrixRMaj tempv = new DMatrixRMaj(n, 1);
        for (int i = 0; i < n; i++) {
            DMatrixRMaj phi = CommonOps_DDRM.extract(Phi, 0, n, i, i+1);
            CommonOps_DDRM.mult(M, phi, tempv);
            double mi = dot(phi, tempv);
            scaleColumn(Phi, i, 1.0 / Math.sqrt(mi));
        }

        double[] w = new double[n];
        for (int i = 0; i < n; i++) w[i] = Math.sqrt(w2s[i]);

        Modes out = new Modes();
        out.w = w;
        out.w2 = w2s;
        out.Phi = Phi;
        return out;
    }

    static double dot(DMatrixRMaj a, DMatrixRMaj b) {
        double s = 0.0;
        for (int i = 0; i < a.numRows; i++) s += a.get(i, 0) * b.get(i, 0);
        return s;
    }

    static void scaleColumn(DMatrixRMaj A, int col, double s) {
        for (int r = 0; r < A.numRows; r++) A.set(r, col, A.get(r, col) * s);
    }

    static void printMatrix(String name, DMatrixRMaj A) {
        System.out.println(name + ":");
        for (int i = 0; i < A.numRows; i++) {
            for (int j = 0; j < A.numCols; j++) {
                System.out.printf("%12.6f ", A.get(i, j));
            }
            System.out.println();
        }
    }

    public static void main(String[] args) {
        double[] masses  = {2.0, 1.5, 1.0};
        double[] springs = {200.0, 300.0, 250.0};

        DMatrixRMaj M = diag(masses);
        DMatrixRMaj K = buildK(springs);

        Modes modes = generalizedModes(M, K);

        System.out.println("Natural frequencies (rad/s):");
        for (int i = 0; i < modes.w.length; i++) {
            System.out.printf("  w%d = %.6f%n", i+1, modes.w[i]);
        }

        // Verify orthogonality: Phi^T M Phi = I and Phi^T K Phi = diag(w^2)
        DMatrixRMaj Mt = new DMatrixRMaj(3, 3);
        DMatrixRMaj Kt = new DMatrixRMaj(3, 3);

        DMatrixRMaj temp = new DMatrixRMaj(3, 3);
        CommonOps_DDRM.multTransA(modes.Phi, M, temp);     // Phi^T M
        CommonOps_DDRM.mult(temp, modes.Phi, Mt);          // (Phi^T M) Phi

        CommonOps_DDRM.multTransA(modes.Phi, K, temp);     // Phi^T K
        CommonOps_DDRM.mult(temp, modes.Phi, Kt);          // (Phi^T K) Phi

        printMatrix("\nPhi^T M Phi (should be I)", Mt);
        printMatrix("\nPhi^T K Phi (should be diag(w^2))", Kt);
    }
}
      

MATLAB directly supports \( \mathbf{K}\boldsymbol{\phi}=\omega^2\mathbf{M}\boldsymbol{\phi} \) using eig(K,M). As a consistency check, you can also build the undamped state-space matrix

\[ \mathbf{x}=\begin{bmatrix}\mathbf{q}\\ \dot{\mathbf{q}}\end{bmatrix}, \quad \dot{\mathbf{x}}= \underbrace{\begin{bmatrix} \mathbf{0} & \mathbf{I}\\ -\mathbf{M}^{-1}\mathbf{K} & \mathbf{0} \end{bmatrix}}_{\mathbf{A}} \mathbf{x}. \]

For the undamped case, \( \operatorname{eig}(\mathbf{A})=\{\pm j\omega_i\} \). In MATLAB: abs(imag(eig(A))) returns the natural frequencies (with duplicates).

Chapter13_Lesson1.m

% Chapter 13 - Vibrations and Multi-Degree-of-Freedom (MDOF) Systems
% Lesson 1: Natural Frequencies and Normal Modes (Eigenvalue Problems)

clear; clc;

m = [2.0; 1.5; 1.0];
k = [200.0; 300.0; 250.0];

M = diag(m);
K = zeros(3,3);

% Assemble chain: wall-k1-m1-k2-m2-k3-m3-(free end)
for i=1:3
    K(i,i) = K(i,i) + k(i);
    if i > 1
        K(i,i)   = K(i,i)   + k(i-1);
        K(i,i-1) = K(i,i-1) - k(i-1);
        K(i-1,i) = K(i-1,i) - k(i-1);
    end
end

% Generalized eigenproblem
[Phi, D] = eig(K, M);     % columns of Phi are modes; D diag of w^2 (not necessarily sorted)
w2 = diag(D);

% Sort ascending
[w2s, idx] = sort(max(w2,0));
Phi = Phi(:, idx);
w = sqrt(w2s);

% Mass-normalize
for i=1:3
    mi = Phi(:,i)'*M*Phi(:,i);
    Phi(:,i) = Phi(:,i)/sqrt(mi);
end

disp('Natural frequencies (rad/s):');
disp(w.');

Mt = Phi'*M*Phi;
Kt = Phi'*K*Phi;

disp('Phi^T M Phi (should be I):');
disp(Mt);

disp('Phi^T K Phi (should be diag(w^2)):');
disp(Kt);

% State-space cross-check
A = [zeros(3) eye(3); -M\K zeros(3)];
lam = eig(A);
w_check = sort(unique(round(abs(imag(lam)), 10)));
disp('Frequencies from state-space eig(A):');
disp(w_check.');

% Free response via modal superposition
t = linspace(0,5,1501);
q0 = [0.02; 0.0; -0.01];
v0 = [0.0; 0.0; 0.0];

eta0 = Phi'*M*q0;
etaDot0 = Phi'*M*v0;

eta = zeros(3, numel(t));
for i=1:3
    if w(i) < 1e-12
        eta(i,:) = eta0(i) + etaDot0(i)*t;
    else
        eta(i,:) = eta0(i)*cos(w(i)*t) + (etaDot0(i)/w(i))*sin(w(i)*t);
    end
end

q = Phi*eta;

figure;
plot(t, q(1,:), t, q(2,:), t, q(3,:));
grid on;
xlabel('t (s)');
ylabel('displacement (m)');
title('3-DOF undamped free vibration (modal superposition)');
legend('q1','q2','q3');
      

Simulink implementation (no figures). Create a State-Space block with matrices \( \mathbf{A} \) above, \( \mathbf{B}=\mathbf{0} \), \( \mathbf{C}=[\mathbf{I}\;\mathbf{0}] \), \( \mathbf{D}=\mathbf{0} \), and set the initial condition \( \mathbf{x}(0)=[\mathbf{q}(0);\dot{\mathbf{q}}(0)] \). The outputs are displacements \( \mathbf{q}(t) \). Compare oscillation frequencies with \( \omega_i \).

Mathematica solves generalized eigenproblems directly via Eigensystem[{K, M}]. We then mass-normalize mode vectors and verify \( \mathbf{\Phi}^\mathsf{T}\mathbf{M}\mathbf{\Phi}=\mathbf{I} \).

Chapter13_Lesson1.nb

(* Define M and K for a 3-DOF chain: wall-k1-m1-k2-m2-k3-m3 *)
m = {2.0, 1.5, 1.0};
k = {200.0, 300.0, 250.0};

M = DiagonalMatrix[m];
K = ConstantArray[0.0, {3, 3}];
Do[
  K[[i, i]] += k[[i]];
  If[i > 1,
    K[[i, i]] += k[[i - 1]];
    K[[i, i - 1]] -= k[[i - 1]];
    K[[i - 1, i]] -= k[[i - 1]];
  ];
, {i, 1, 3}];

{vals, vecs} = Eigensystem[{K, M}];
vals = Map[Max[#, 0.0] &, vals];
omega = Sqrt[vals];

(* Sort by eigenvalue *)
perm = Ordering[vals];
vals = vals[[perm]];
omega = omega[[perm]];
vecs = vecs[[perm]];

(* Mass-normalize modes *)
Phi = Map[#/Sqrt[#.M.#] &, vecs];
PhiMat = Transpose[Phi];

Morth = Chop[PhiMat.M.Transpose[PhiMat]];  (* Phi^T M Phi *)
Korth = Chop[PhiMat.K.Transpose[PhiMat]];  (* Phi^T K Phi *)

{omega, Morth, Korth}
      

Problem 1 (2-DOF analytic eigenproblem): Consider two masses \(m\) connected in a line by three springs: wall–\(k\)–mass1–\(k\)–mass2–\(k\)–wall. (a) Write \( \mathbf{M} \) and \( \mathbf{K} \). (b) Compute the natural frequencies \( \omega_1,\omega_2 \) and one eigenvector for each mode.

Solution:

The generalized coordinates are \( \mathbf{q}=[q_1\;q_2]^\mathsf{T} \). The matrices are

\[ \mathbf{M}=\begin{bmatrix}m & 0\\ 0 & m\end{bmatrix},\quad \mathbf{K}=\begin{bmatrix}2k & -k\\ -k & 2k\end{bmatrix}. \]

Solve \( \det(\mathbf{K}-\omega^2\mathbf{M})=0 \):

\[ \det\!\begin{bmatrix}2k-m\omega^2 & -k\\ -k & 2k-m\omega^2\end{bmatrix} =(2k-m\omega^2)^2-k^2=0. \]

Hence \( 2k-m\omega^2=\pm k \), giving

\[ \omega_1^2=\frac{k}{m},\quad \omega_2^2=\frac{3k}{m} \;\;\Longrightarrow\;\; \omega_1=\sqrt{\frac{k}{m}},\quad \omega_2=\sqrt{\frac{3k}{m}}. \]

For \( \omega_1^2=k/m \), solve \( (\mathbf{K}-\omega_1^2\mathbf{M})\boldsymbol{\phi}_1=\mathbf{0} \):

\[ \begin{bmatrix}k & -k\\ -k & k\end{bmatrix}\boldsymbol{\phi}_1=\mathbf{0} \;\;\Longrightarrow\;\; \boldsymbol{\phi}_1\propto \begin{bmatrix}1\\ 1\end{bmatrix}. \]

For \( \omega_2^2=3k/m \), \( \boldsymbol{\phi}_2\propto [1\;-1]^\mathsf{T} \).

Problem 2 (orthogonality proof): Assume \( \mathbf{M}\succ 0 \), \( \mathbf{K}=\mathbf{K}^\mathsf{T} \), and two eigenpairs \( \mathbf{K}\boldsymbol{\phi}_i=\lambda_i\mathbf{M}\boldsymbol{\phi}_i \), \( \mathbf{K}\boldsymbol{\phi}_j=\lambda_j\mathbf{M}\boldsymbol{\phi}_j \) with \( \lambda_i\neq\lambda_j \). Prove \( \boldsymbol{\phi}_i^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_j=0 \).

Solution:

As shown in Section 4, symmetry of \( \mathbf{K} \) implies \( \boldsymbol{\phi}_j^\mathsf{T}\mathbf{K}\boldsymbol{\phi}_i=\boldsymbol{\phi}_i^\mathsf{T}\mathbf{K}\boldsymbol{\phi}_j \). Using the eigenrelations yields \( (\lambda_i-\lambda_j)\boldsymbol{\phi}_i^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_j=0 \), hence \( \boldsymbol{\phi}_i^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_j=0 \).

Problem 3 (Rayleigh-quotient bound): Let \( \lambda_1=\min \lambda_i \) be the smallest eigenvalue of \( \mathbf{K}\boldsymbol{\phi}=\lambda\mathbf{M}\boldsymbol{\phi} \) with \( \mathbf{M}\succ 0 \). Show that for any nonzero \( \mathbf{x} \), \( \lambda_1 \le \mathcal{R}(\mathbf{x}) \), where \( \mathcal{R}(\mathbf{x})=\frac{\mathbf{x}^\mathsf{T}\mathbf{K}\mathbf{x}}{\mathbf{x}^\mathsf{T}\mathbf{M}\mathbf{x}} \).

Solution:

Expand \( \mathbf{x} \) in the mass-orthonormal eigenbasis: \( \mathbf{x}=\sum_{i=1}^n c_i\boldsymbol{\phi}_i \) with \( \boldsymbol{\phi}_i^\mathsf{T}\mathbf{M}\boldsymbol{\phi}_j=\delta_{ij} \). Then

\[ \mathbf{x}^\mathsf{T}\mathbf{M}\mathbf{x}=\sum_{i=1}^n c_i^2,\quad \mathbf{x}^\mathsf{T}\mathbf{K}\mathbf{x}=\sum_{i=1}^n c_i^2\lambda_i. \]

Therefore

\[ \mathcal{R}(\mathbf{x})=\frac{\sum_{i=1}^n c_i^2\lambda_i}{\sum_{i=1}^n c_i^2} \ge \frac{\lambda_1\sum_{i=1}^n c_i^2}{\sum_{i=1}^n c_i^2}=\lambda_1. \]

Thus \( \mathcal{R}(\mathbf{x}) \) is an upper bound on the smallest eigenvalue, and the bound is tight when \( \mathbf{x} \) aligns with \( \boldsymbol{\phi}_1 \).

Problem 4 (Cholesky reduction equivalence): Let \( \mathbf{M}=\mathbf{L}\mathbf{L}^\mathsf{T} \) and define \( \mathbf{A}=\mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-\mathsf{T}} \). Prove that \( \lambda \) is an eigenvalue of \( (\mathbf{K},\mathbf{M}) \) iff it is an eigenvalue of \( \mathbf{A} \).

Solution:

(Forward) If \( \mathbf{K}\boldsymbol{\phi}=\lambda\mathbf{M}\boldsymbol{\phi} \), set \( \mathbf{u}=\mathbf{L}^\mathsf{T}\boldsymbol{\phi} \). Then

\[ \mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-\mathsf{T}}\mathbf{u} = \mathbf{L}^{-1}\mathbf{K}\boldsymbol{\phi} = \lambda\,\mathbf{L}^{-1}\mathbf{M}\boldsymbol{\phi} = \lambda\,\mathbf{L}^{-1}\mathbf{L}\mathbf{L}^\mathsf{T}\boldsymbol{\phi} = \lambda\,\mathbf{u}. \]

Hence \( \mathbf{A}\mathbf{u}=\lambda\mathbf{u} \). (Reverse) If \( \mathbf{A}\mathbf{u}=\lambda\mathbf{u} \), set \( \boldsymbol{\phi}=\mathbf{L}^{-\mathsf{T}}\mathbf{u} \) and reverse the algebra to get \( \mathbf{K}\boldsymbol{\phi}=\lambda\mathbf{M}\boldsymbol{\phi} \).

Problem 5 (rigid-body mode as zero eigenvalue): Assume \( \mathbf{M}\succ 0 \) and \( \mathbf{K}\succeq 0 \) with \( \mathbf{K}\mathbf{r}=\mathbf{0} \) for some \( \mathbf{r}\neq\mathbf{0} \). Show that \( \lambda=0 \) is an eigenvalue of the generalized problem and interpret \( \omega \).

Solution:

Taking \( \boldsymbol{\phi}=\mathbf{r} \), the relation \( \mathbf{K}\boldsymbol{\phi}=\lambda\mathbf{M}\boldsymbol{\phi} \) becomes \( \mathbf{0}=\lambda\mathbf{M}\mathbf{r} \). Since \( \mathbf{M}\mathbf{r}\neq\mathbf{0} \) (SPD), we must have \( \lambda=0 \), hence \( \omega=0 \). Physically, this indicates a rigid-body motion that stores no elastic energy (no restoring force).

We derived the undamped MDOF vibration equation \( \mathbf{M}\ddot{\mathbf{q}}+\mathbf{K}\mathbf{q}=0 \), formulated the generalized eigenproblem \( \mathbf{K}\boldsymbol{\phi}=\omega^2\mathbf{M}\boldsymbol{\phi} \), and proved that (for symmetric \( \mathbf{M}\succ 0 \), \( \mathbf{K}\succeq 0 \)) eigenvalues are real and nonnegative, with distinct modes orthogonal under the mass inner product. We also presented the Cholesky reduction to a standard symmetric eigenproblem and implemented the full workflow in Python, C++, Java, MATLAB/Simulink, and Mathematica.