Consider an undamped linear MDOF vibration model with \( n \) generalized coordinates collected in \( \mathbf{x}(t)\in\mathbb{R}^n \). The equations of motion are

\[ \mathbf{M}\,\ddot{\mathbf{x} }(t) + \mathbf{K}\,\mathbf{x}(t) = \mathbf{f}(t), \quad \mathbf{M}=\mathbf{M}^\top \succ 0,\;\; \mathbf{K}=\mathbf{K}^\top \succeq 0. \]

The modal decoupling objective is to find an invertible coordinate transform \( \mathbf{x}=\mathbf{\Phi}\,\mathbf{q} \) such that the dynamics in \( \mathbf{q}(t) \) become (ideally) diagonal, i.e., a set of independent scalar oscillators:

\[ \ddot{q}_i(t) + \omega_i^2\,q_i(t) = p_i(t), \quad i=1,\dots,n, \]

where \( \omega_i \) are natural frequencies and \( p_i(t) \) are modal forces. This transformation is a change of basis in the configuration space: instead of describing motion using physical coordinates \( \mathbf{x} \), we describe motion as a combination of mode shapes (columns of \( \mathbf{\Phi} \)) weighted by modal coordinates \( \mathbf{q} \).

From Lesson 2 (natural frequencies and normal modes), free vibration \( \mathbf{f}(t)=\mathbf{0} \) is sought in the harmonic form \( \mathbf{x}(t)=\boldsymbol{\phi}\,e^{j\omega t} \), yielding the generalized eigenproblem:

\[ \left(\mathbf{K}-\omega^2\mathbf{M}\right)\boldsymbol{\phi}=\mathbf{0} \quad \Longleftrightarrow \quad \mathbf{K}\boldsymbol{\phi} = \lambda\,\mathbf{M}\boldsymbol{\phi},\;\; \lambda=\omega^2. \]

For symmetric \( \mathbf{M}\succ 0 \) and symmetric \( \mathbf{K} \), the eigenvalues \( \lambda_i \) are real and (for typical structural systems with sufficient constraints) nonnegative. Let \( \{(\lambda_i,\boldsymbol{\phi}_i)\}_{i=1}^n \) denote eigenpairs with \( \lambda_i=\omega_i^2 \). We form the modal matrix

\[ \mathbf{\Phi} = \begin{bmatrix} \boldsymbol{\phi}_1 & \boldsymbol{\phi}_2 & \cdots & \boldsymbol{\phi}_n \end{bmatrix}. \]

The key property enabling decoupling is orthogonality with respect to \( \mathbf{M} \) (and also \( \mathbf{K} \)), proven next.

Theorem (Mass orthogonality for distinct eigenvalues): If \( \mathbf{K}\boldsymbol{\phi}_i=\lambda_i\mathbf{M}\boldsymbol{\phi}_i \) and \( \mathbf{K}\boldsymbol{\phi}_j=\lambda_j\mathbf{M}\boldsymbol{\phi}_j \) with \( \lambda_i \neq \lambda_j \), then \( \boldsymbol{\phi}_i^\top \mathbf{M}\boldsymbol{\phi}_j = 0 \). Moreover, \( \boldsymbol{\phi}_i^\top \mathbf{K}\boldsymbol{\phi}_j = 0 \).

Proof: Start with the two eigen-relations:

\[ \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 by \( \boldsymbol{\phi}_j^\top \) and the second by \( \boldsymbol{\phi}_i^\top \):

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

Because \( \mathbf{K}=\mathbf{K}^\top \), we have \( \boldsymbol{\phi}_j^\top\mathbf{K}\boldsymbol{\phi}_i = (\boldsymbol{\phi}_i^\top\mathbf{K}\boldsymbol{\phi}_j)^\top = \boldsymbol{\phi}_i^\top\mathbf{K}\boldsymbol{\phi}_j \), and because \( \mathbf{M}=\mathbf{M}^\top \), \( \boldsymbol{\phi}_j^\top\mathbf{M}\boldsymbol{\phi}_i = \boldsymbol{\phi}_i^\top\mathbf{M}\boldsymbol{\phi}_j \). Subtract the two equations:

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

Since \( \lambda_i\neq\lambda_j \), it follows that \( \boldsymbol{\phi}_i^\top\mathbf{M}\boldsymbol{\phi}_j = 0 \). Finally, \( \boldsymbol{\phi}_i^\top\mathbf{K}\boldsymbol{\phi}_j = \lambda_j \boldsymbol{\phi}_i^\top\mathbf{M}\boldsymbol{\phi}_j = 0 \), proving stiffness orthogonality as well. ∎

Repeated eigenvalues: if \( \lambda_i=\lambda_j \), the above argument does not force orthogonality, but one can always choose a basis inside that eigenspace to be \( \mathbf{M} \)-orthonormal via a Gram–Schmidt process using the inner product \( \langle \mathbf{u},\mathbf{v}\rangle_M = \mathbf{u}^\top\mathbf{M}\mathbf{v} \).

The eigenvectors can be scaled arbitrarily. A particularly useful scaling is mass normalization:

\[ \boldsymbol{\phi}_i^\top\mathbf{M}\boldsymbol{\phi}_i = 1 \quad \Longrightarrow \quad \mathbf{\Phi}^\top\mathbf{M}\mathbf{\Phi} = \mathbf{I}. \]

Under this choice, stiffness orthogonality becomes:

\[ \mathbf{\Phi}^\top\mathbf{K}\mathbf{\Phi} = \begin{bmatrix} \omega_1^2 & 0 & \cdots & 0\\ 0 & \omega_2^2 & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \omega_n^2 \end{bmatrix} \;=\;\mathbf{\Omega}^2. \]

So the modal matrix simultaneously diagonalizes \( \mathbf{M} \) and \( \mathbf{K} \) under congruence transformations (transpose times matrix times itself). This is the mathematical core of decoupling.

A useful “inverse transform” identity follows immediately. If \( \mathbf{x}=\mathbf{\Phi}\mathbf{q} \) and \( \mathbf{\Phi}^\top\mathbf{M}\mathbf{\Phi}=\mathbf{I} \), then left-multiplying by \( \mathbf{\Phi}^\top\mathbf{M} \) yields \( \mathbf{q}=\mathbf{\Phi}^\top\mathbf{M}\mathbf{x} \). Thus, modal coordinates are simply physical coordinates projected using the \( \mathbf{M} \)-inner product.

Apply the transformation \( \mathbf{x}(t)=\mathbf{\Phi}\mathbf{q}(t) \). Since \( \mathbf{\Phi} \) is constant (time-invariant), we have \( \dot{\mathbf{x} }=\mathbf{\Phi}\dot{\mathbf{q} } \) and \( \ddot{\mathbf{x} }=\mathbf{\Phi}\ddot{\mathbf{q} } \). Substitute into \( \mathbf{M}\ddot{\mathbf{x} }+\mathbf{K}\mathbf{x}=\mathbf{f}(t) \):

\[ \mathbf{M}\mathbf{\Phi}\ddot{\mathbf{q} }(t) + \mathbf{K}\mathbf{\Phi}\mathbf{q}(t) = \mathbf{f}(t). \]

Now left-multiply by \( \mathbf{\Phi}^\top \):

\[ \mathbf{\Phi}^\top\mathbf{M}\mathbf{\Phi}\ddot{\mathbf{q} }(t) + \mathbf{\Phi}^\top\mathbf{K}\mathbf{\Phi}\mathbf{q}(t) = \mathbf{\Phi}^\top\mathbf{f}(t). \]

Using the mass-normalization identities from Section 4, we obtain the decoupled modal equations:

\[ \ddot{\mathbf{q} }(t) + \mathbf{\Omega}^2\,\mathbf{q}(t) = \mathbf{p}(t), \qquad \mathbf{p}(t)=\mathbf{\Phi}^\top\mathbf{f}(t). \]

Component-wise, this is: \( \ddot{q}_i(t) + \omega_i^2 q_i(t) = p_i(t) \), where \( p_i(t) = \boldsymbol{\phi}_i^\top\mathbf{f}(t) \). The system is now a set of independent scalar forced oscillators (undamped).

flowchart TD
  A["Physical model: M xdd + K x = f(t)"] --> B["Solve EVP: K phi = lambda M phi"]
  B --> C["Assemble Phi = [phi1 ... phin]"]
  C --> D["Mass-normalize: Phi^T M Phi = I"]
  D --> E["Transform: x = Phi q"]
  E --> F["Decoupled: qdd + Omega^2 q = Phi^T f(t)"]
        

With mass-normalized modes, the kinetic and potential energies separate into sums over modes. Recall: \( T=\tfrac{1}{2}\dot{\mathbf{x} }^\top\mathbf{M}\dot{\mathbf{x} } \), \( V=\tfrac{1}{2}\mathbf{x}^\top\mathbf{K}\mathbf{x} \). Substitute \( \mathbf{x}=\mathbf{\Phi}\mathbf{q} \) and \( \dot{\mathbf{x} }=\mathbf{\Phi}\dot{\mathbf{q} } \):

\[ T=\tfrac{1}{2}\dot{\mathbf{q} }^\top\left(\mathbf{\Phi}^\top\mathbf{M}\mathbf{\Phi}\right)\dot{\mathbf{q} } =\tfrac{1}{2}\dot{\mathbf{q} }^\top\dot{\mathbf{q} } =\tfrac{1}{2}\sum_{i=1}^n \dot{q}_i^2, \]

\[ V=\tfrac{1}{2}\mathbf{q}^\top\left(\mathbf{\Phi}^\top\mathbf{K}\mathbf{\Phi}\right)\mathbf{q} =\tfrac{1}{2}\mathbf{q}^\top\mathbf{\Omega}^2\mathbf{q} =\tfrac{1}{2}\sum_{i=1}^n \omega_i^2 q_i^2. \]

Hence each modal coordinate behaves like an independent oscillator with “unit modal mass” and stiffness \( \omega_i^2 \). The physical motion is a modal superposition: \( \mathbf{x}(t)=\sum_{i=1}^n \boldsymbol{\phi}_i q_i(t) \).

For forcing, \( p_i(t)=\boldsymbol{\phi}_i^\top\mathbf{f}(t) \) is the generalized force under the \( \mathbf{M} \)-orthonormal basis. If a particular mode shape has small projection on the forcing direction, it receives little energy (small modal participation).

Each modal equation is a scalar undamped forced oscillator: \( \ddot{q}_i + \omega_i^2 q_i = p_i(t) \). For general \( p_i(t) \), the solution can be written via the impulse response \( h_i(t)=\frac{1}{\omega_i}\sin(\omega_i t) \):

\[ q_i(t)=q_i(0)\cos(\omega_i t)+\frac{\dot{q}_i(0)}{\omega_i}\sin(\omega_i t) +\int_{0}^{t}\frac{1}{\omega_i}\sin\!\big(\omega_i(t-\theta)\big)\,p_i(\theta)\,d\theta. \]

For harmonic forcing \( p_i(t)=\hat{p}_i\sin(\Omega t) \) with \( \Omega \neq \omega_i \), the steady-state amplitude is governed by the undamped frequency response:

\[ q_i^{ss}(t)=\frac{\hat{p}_i}{\omega_i^2-\Omega^2}\sin(\Omega t), \qquad \Omega \neq \omega_i. \]

The singularity as \( \Omega \to \omega_i \) is the classical undamped resonance artifact. In practice, damping (next lesson) regularizes this behavior.

flowchart TD
  S["Given M, K, f(t), x0, xd0"] --> E["Generalized EVP (K, M) -> lambdas, phis"]
  E --> N["Normalize modes: Phi^T M Phi = I"]
  N --> T["Transform ICs: q0 = Phi^T M x0, qd0 = Phi^T M xd0"]
  T --> D["Decoupled ODEs: qdd_i + wi^2 qi = pi(t)"]
  D --> R["Reconstruct: x(t) = Phi q(t)"]
  R --> C["Check: compare vs direct integration (optional)"]
        

Below we implement a complete example: a 3-DOF mass–spring chain with boundary springs, compute \( \omega_i \) and \( \mathbf{\Phi} \), mass-normalize the modes, simulate the decoupled modal ODEs, and reconstruct \( \mathbf{x}(t) \).

8.1 Python (NumPy/SciPy): generalized symmetric eigenproblem + ODE integration

Library notes: Use scipy.linalg.eigh for symmetric generalized EVP and scipy.integrate.solve_ivp for ODEs.

File: Chapter13_Lesson3.py

"""
Chapter13_Lesson3.py
System Dynamics — Chapter 13, Lesson 3
Modal Coordinates and Decoupling of MDOF Systems (Undamped Case)

Dependencies:
  numpy, scipy, matplotlib
"""
import numpy as np
from numpy.linalg import norm
from scipy.linalg import eigh
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


def mass_normalize(Phi: np.ndarray, M: np.ndarray) -> np.ndarray:
    """Scale columns of Phi so that Phi^T M Phi = I."""
    Phi_n = Phi.copy()
    for i in range(Phi.shape[1]):
        mi = Phi[:, i].T @ M @ Phi[:, i]
        Phi_n[:, i] = Phi[:, i] / np.sqrt(mi)
    return Phi_n


def build_3dof_chain(m=(1.0, 1.2, 0.9), k=(1200.0, 900.0, 700.0, 1100.0)):
    """
    wall --k1-- m1 --k2-- m2 --k3-- m3 --k4-- wall
    """
    m1, m2, m3 = m
    k1, k2, k3, k4 = k
    M = np.diag([m1, m2, m3])
    K = np.array([
        [k1 + k2,   -k2,       0.0],
        [-k2,     k2 + k3,   -k3],
        [0.0,       -k3,     k3 + k4]
    ], dtype=float)
    return M, K


def modal_matrices(M: np.ndarray, K: np.ndarray):
    """Solve K phi = lambda M phi; return omega, Phi_mass_norm, Phi^T K Phi."""
    lam, Phi = eigh(K, M)  # ascending eigenvalues
    lam = np.maximum(lam, 0.0)
    omega = np.sqrt(lam)
    Phi = mass_normalize(Phi, M)
    Omega2 = Phi.T @ K @ Phi
    return omega, Phi, Omega2


def simulate_modal(omega, Phi, force_fun, t_span=(0.0, 8.0), n_eval=2000):
    """
    Modal ODE (mass-normalized):
      qdd + omega^2 q = p(t),  p(t) = Phi^T f(t)
    """
    n = Phi.shape[0]

    def ode(t, z):
        q = z[:n]
        qd = z[n:]
        p = Phi.T @ force_fun(t)
        qdd = -(omega**2) * q + p
        return np.hstack([qd, qdd])

    t_eval = np.linspace(t_span[0], t_span[1], n_eval)
    z0 = np.zeros(2*n)
    sol = solve_ivp(ode, t_span, z0, t_eval=t_eval, rtol=1e-8, atol=1e-10)

    q = sol.y[:n, :]
    x = Phi @ q
    return sol.t, x, q


def main():
    M, K = build_3dof_chain()
    omega, Phi, Omega2 = modal_matrices(M, K)

    print("Natural frequencies (rad/s):", omega)
    print("Check Phi^T M Phi ~ I:\\n", Phi.T @ M @ Phi)
    print("Check Phi^T K Phi ~ diag(omega^2):\\n", Omega2)

    F0 = 10.0
    Omega = 0.9 * omega[0]  # excite near mode 1
    f = lambda t: np.array([F0*np.sin(Omega*t), 0.0, 0.0])

    t, x, q = simulate_modal(omega, Phi, f)

    plt.figure()
    plt.plot(t, x[0, :], label="x1(t)")
    plt.xlabel("t [s]"); plt.ylabel("displacement")
    plt.grid(True); plt.legend(); plt.title("DOF 1 displacement via modal coordinates")

    plt.figure()
    for i in range(3):
        plt.plot(t, q[i, :], label=f"q{i+1}(t)")
    plt.xlabel("t [s]"); plt.ylabel("modal coordinate")
    plt.grid(True); plt.legend(); plt.title("Modal coordinates")
    plt.show()


if __name__ == "__main__":
    main()
      

8.2 C++ (Eigen): generalized self-adjoint eigenproblem + RK4

Library notes: Eigen provides GeneralizedSelfAdjointEigenSolver when \( \mathbf{M} \) and \( \mathbf{K} \) are symmetric.

File: Chapter13_Lesson3.cpp

/*
Chapter13_Lesson3.cpp
Modal Coordinates and Decoupling of MDOF Systems (Undamped Case)

Dependency: Eigen (header-only)
*/
#include <Eigen/Dense>
#include <iostream>
#include <cmath>

using Eigen::MatrixXd;
using Eigen::VectorXd;

static void massNormalize(MatrixXd& Phi, const MatrixXd& M) {
    for (int i = 0; i < Phi.cols(); ++i) {
        double mi = Phi.col(i).transpose() * M * Phi.col(i);
        Phi.col(i) /= std::sqrt(mi);
    }
}

static void build3DOF(MatrixXd& M, MatrixXd& K) {
    double m1 = 1.0, m2 = 1.2, m3 = 0.9;
    double k1 = 1200.0, k2 = 900.0, k3 = 700.0, k4 = 1100.0;

    M = MatrixXd::Zero(3,3);
    M(0,0)=m1; M(1,1)=m2; M(2,2)=m3;

    K = MatrixXd::Zero(3,3);
    K(0,0)=k1+k2; K(0,1)=-k2;
    K(1,0)=-k2;   K(1,1)=k2+k3; K(1,2)=-k3;
    K(2,1)=-k3;   K(2,2)=k3+k4;
}

static VectorXd force(double t, double F0, double Omega) {
    VectorXd f(3);
    f << F0*std::sin(Omega*t), 0.0, 0.0;
    return f;
}

// RK4 for z=[q;qd], z'=[qd; -w^2 q + p(t)]
static void rk4_step(VectorXd& q, VectorXd& qd,
                     const VectorXd& omega,
                     const VectorXd& p, double dt)
{
    const int n = q.size();
    auto acc = [&](const VectorXd& q_in, const VectorXd& p_in) {
        VectorXd qdd(n);
        for (int i = 0; i < n; ++i) qdd(i) = -omega(i)*omega(i)*q_in(i) + p_in(i);
        return qdd;
    };

    VectorXd k1_q = qd;
    VectorXd k1_qd = acc(q, p);

    VectorXd q2 = q + 0.5*dt*k1_q;
    VectorXd qd2 = qd + 0.5*dt*k1_qd;
    VectorXd k2_q = qd2;
    VectorXd k2_qd = acc(q2, p);

    VectorXd q3 = q + 0.5*dt*k2_q;
    VectorXd qd3 = qd + 0.5*dt*k2_qd;
    VectorXd k3_q = qd3;
    VectorXd k3_qd = acc(q3, p);

    VectorXd q4 = q + dt*k3_q;
    VectorXd qd4 = qd + dt*k3_qd;
    VectorXd k4_q = qd4;
    VectorXd k4_qd = acc(q4, p);

    q  += (dt/6.0)*(k1_q  + 2.0*k2_q  + 2.0*k3_q  + k4_q);
    qd += (dt/6.0)*(k1_qd + 2.0*k2_qd + 2.0*k3_qd + k4_qd);
}

int main() {
    MatrixXd M, K;
    build3DOF(M, K);

    Eigen::GeneralizedSelfAdjointEigenSolver<MatrixXd> ges(K, M);
    VectorXd lambda = ges.eigenvalues();
    MatrixXd Phi = ges.eigenvectors();

    VectorXd omega = lambda.cwiseMax(0.0).cwiseSqrt();
    massNormalize(Phi, M);

    std::cout << "Natural frequencies (rad/s): " << omega.transpose() << "\\n";
    std::cout << "Phi^T M Phi:\\n" << (Phi.transpose()*M*Phi) << "\\n";

    double F0 = 10.0;
    double Omega = 0.9 * omega(0);

    VectorXd q = VectorXd::Zero(3);
    VectorXd qd = VectorXd::Zero(3);

    double t0 = 0.0, tf = 8.0, dt = 1e-3;
    int N = static_cast<int>((tf-t0)/dt);

    for (int step = 0; step <= N; ++step) {
        double t = t0 + step*dt;
        VectorXd p = Phi.transpose() * force(t, F0, Omega);

        if (step % 1000 == 0) {
            VectorXd x = Phi * q;
            std::cout << "t=" << t << "  x1=" << x(0) << "  q1=" << q(0) << "\\n";
        }
        rk4_step(q, qd, omega, p, dt);
    }
    return 0;
}
      

8.3 Java (EJML): Cholesky reduction to standard eigenproblem + RK4

Library notes: For symmetric generalized EVP, a stable approach is: \( \mathbf{M}=\mathbf{L}\mathbf{L}^\top \), then solve \( \mathbf{A}=\mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-\top} \).

File: Chapter13_Lesson3.java

/*
Chapter13_Lesson3.java
Modal Coordinates and Decoupling of MDOF Systems (Undamped Case)

Dependency: EJML (org.ejml:ejml-all)
*/
import org.ejml.data.DMatrixRMaj;
import org.ejml.dense.row.CommonOps_DDRM;
import org.ejml.dense.row.decomposition.DecompositionFactory_DDRM;
import org.ejml.interfaces.decomposition.CholeskyDecomposition_F64;
import org.ejml.interfaces.decomposition.EigenDecomposition_F64;

import java.util.Arrays;

public class Chapter13_Lesson3 {

    static void build3DOF(DMatrixRMaj M, DMatrixRMaj K) {
        double m1 = 1.0, m2 = 1.2, m3 = 0.9;
        double k1 = 1200.0, k2 = 900.0, k3 = 700.0, k4 = 1100.0;

        M.reshape(3,3); CommonOps_DDRM.fill(M,0);
        M.set(0,0,m1); M.set(1,1,m2); M.set(2,2,m3);

        K.reshape(3,3); CommonOps_DDRM.fill(K,0);
        K.set(0,0,k1+k2); K.set(0,1,-k2);
        K.set(1,0,-k2);   K.set(1,1,k2+k3); K.set(1,2,-k3);
        K.set(2,1,-k3);   K.set(2,2,k3+k4);
    }

    static DMatrixRMaj invLower(DMatrixRMaj L) {
        int n = L.numRows;
        DMatrixRMaj invL = CommonOps_DDRM.identity(n);
        CommonOps_DDRM.solve(L, invL, invL); // L * invL = I
        return invL;
    }

    static class ModalResult { double[] omega; DMatrixRMaj Phi; }

    static ModalResult modal(DMatrixRMaj M, DMatrixRMaj K) {
        int n = M.numRows;

        CholeskyDecomposition_F64<DMatrixRMaj> chol = DecompositionFactory_DDRM.chol(n, true);
        if (!chol.decompose(M.copy())) throw new RuntimeException("M not SPD.");
        DMatrixRMaj L = chol.getT(null);

        DMatrixRMaj invL = invLower(L);

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

        // Eigen-decomposition of A (symmetric)
        EigenDecomposition_F64<DMatrixRMaj> eig = DecompositionFactory_DDRM.eig(n, true);
        if (!eig.decompose(A)) throw new RuntimeException("eig failed");

        double[] lambda = new double[n];
        DMatrixRMaj Y = new DMatrixRMaj(n,n);
        for (int i=0;i<n;i++){
            lambda[i] = eig.getEigenvalue(i).getReal();
            DMatrixRMaj yi = eig.getEigenVector(i);
            for(int r=0;r<n;r++) Y.set(r,i, yi.get(r,0));
        }

        // Recover generalized eigenvectors: Phi = inv(L^T) * Y = (invL^T)*Y
        DMatrixRMaj invLT = new DMatrixRMaj(n,n);
        CommonOps_DDRM.transpose(invL, invLT);
        DMatrixRMaj Phi = new DMatrixRMaj(n,n);
        CommonOps_DDRM.mult(invLT, Y, Phi);

        // Mass-normalize: phi_i^T M phi_i = 1
        DMatrixRMaj MPhi = new DMatrixRMaj(n,n);
        CommonOps_DDRM.mult(M, Phi, MPhi);
        for(int i=0;i<n;i++){
            double mi = 0;
            for(int r=0;r<n;r++) mi += Phi.get(r,i)*MPhi.get(r,i);
            double s = 1.0/Math.sqrt(mi);
            for(int r=0;r<n;r++) Phi.set(r,i, Phi.get(r,i)*s);
        }

        double[] omega = new double[n];
        for(int i=0;i<n;i++) omega[i] = Math.sqrt(Math.max(lambda[i],0));

        ModalResult res = new ModalResult();
        res.omega = omega; res.Phi = Phi;
        return res;
    }

    static double[] force(double t, double F0, double Omega){
        return new double[]{F0*Math.sin(Omega*t), 0.0, 0.0};
    }

    static double[] pModal(DMatrixRMaj Phi, double[] f){
        int n = Phi.numRows;
        double[] p = new double[n];
        for(int i=0;i<n;i++){
            double s=0;
            for(int r=0;r<n;r++) s += Phi.get(r,i)*f[r];
            p[i]=s;
        }
        return p;
    }

    static void rk4(double[] q, double[] qd, double[] omega, double[] p, double dt){
        int n=q.length;
        java.util.function.BiFunction<double[], double[], double[]> acc = (qq,pp) -> {
            double[] qdd = new double[n];
            for(int i=0;i<n;i++) qdd[i] = -omega[i]*omega[i]*qq[i] + pp[i];
            return qdd;
        };

        double[] k1q = Arrays.copyOf(qd,n);
        double[] k1qd = acc.apply(q,p);

        double[] q2=new double[n], qd2=new double[n];
        for(int i=0;i<n;i++){ q2[i]=q[i]+0.5*dt*k1q[i]; qd2[i]=qd[i]+0.5*dt*k1qd[i]; }
        double[] k2q = Arrays.copyOf(qd2,n);
        double[] k2qd = acc.apply(q2,p);

        double[] q3=new double[n], qd3=new double[n];
        for(int i=0;i<n;i++){ q3[i]=q[i]+0.5*dt*k2q[i]; qd3[i]=qd[i]+0.5*dt*k2qd[i]; }
        double[] k3q = Arrays.copyOf(qd3,n);
        double[] k3qd = acc.apply(q3,p);

        double[] q4=new double[n], qd4=new double[n];
        for(int i=0;i<n;i++){ q4[i]=q[i]+dt*k3q[i]; qd4[i]=qd[i]+dt*k3qd[i]; }
        double[] k4q = Arrays.copyOf(qd4,n);
        double[] k4qd = acc.apply(q4,p);

        for(int i=0;i<n;i++){
            q[i]  += (dt/6.0)*(k1q[i]  + 2*k2q[i]  + 2*k3q[i]  + k4q[i]);
            qd[i] += (dt/6.0)*(k1qd[i] + 2*k2qd[i] + 2*k3qd[i] + k4qd[i]);
        }
    }

    static double[] xReconstruct(DMatrixRMaj Phi, double[] q){
        int n=Phi.numRows;
        double[] x=new double[n];
        for(int r=0;r<n;r++){
            double s=0;
            for(int i=0;i<n;i++) s += Phi.get(r,i)*q[i];
            x[r]=s;
        }
        return x;
    }

    public static void main(String[] args){
        DMatrixRMaj M=new DMatrixRMaj(3,3);
        DMatrixRMaj K=new DMatrixRMaj(3,3);
        build3DOF(M,K);

        ModalResult mr = modal(M,K);
        System.out.println("Natural frequencies (rad/s): "+Arrays.toString(mr.omega));

        double F0=10.0;
        double Omega=0.9*mr.omega[0];

        double[] q={0,0,0}, qd={0,0,0};

        double t0=0, tf=8, dt=1e-3;
        int N=(int)Math.round((tf-t0)/dt);

        for(int step=0; step<=N; step++){
            double t=t0+step*dt;
            double[] p = pModal(mr.Phi, force(t,F0,Omega));

            if(step%1000==0){
                double[] x = xReconstruct(mr.Phi, q);
                System.out.printf("t=%.3f  x1=%.6f  q1=%.6f%n", t, x[0], q[0]);
            }
            rk4(q, qd, mr.omega, p, dt);
        }
    }
}
      

8.4 MATLAB + Simulink: eig(K,M) + ode45 + programmatic model build

Library notes: MATLAB natively supports generalized EVP with eig(K,M). Simulink can represent the decoupled system with a single State-Space block.

File: Chapter13_Lesson3.m

% Chapter13_Lesson3.m
% Modal Coordinates and Decoupling of MDOF Systems (Undamped Case)

clear; clc;

m1 = 1.0; m2 = 1.2; m3 = 0.9;
k1 = 1200; k2 = 900; k3 = 700; k4 = 1100;

M = diag([m1 m2 m3]);
K = [k1+k2, -k2,   0;
     -k2,  k2+k3, -k3;
      0,    -k3,  k3+k4];

[Phi,Lam] = eig(K,M);
lam = diag(Lam);
[lam, idx] = sort(lam,'ascend');
Phi = Phi(:,idx);
omega = sqrt(max(lam,0));

for i=1:size(Phi,2)
    mi = Phi(:,i)'*M*Phi(:,i);
    Phi(:,i) = Phi(:,i)/sqrt(mi);
end

Omega2 = Phi'*K*Phi;

F0 = 10; Omega = 0.9*omega(1);
f = @(t)[F0*sin(Omega*t); 0; 0];
p = @(t)Phi'*f(t);

n = 3;
ode = @(t,z)[ z(n+1:end);
             - (omega.^2).*z(1:n) + p(t) ];

tspan = [0 8];
z0 = zeros(2*n,1);

opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,z] = ode45(ode,tspan,z0,opts);
q = z(:,1:n).';
x = Phi*q;

figure; plot(t,x(1,:), 'LineWidth',1.2); grid on;
xlabel('t [s]'); ylabel('x_1'); title('x_1(t) via modal simulation');

% Optional Simulink auto-build:
mdl = 'Chapter13_Lesson3_Simulink';
if bdIsLoaded(mdl); close_system(mdl,0); end
new_system(mdl); open_system(mdl);

A = [zeros(n), eye(n);
     -diag(omega.^2), zeros(n)];
B = [zeros(n,1); Phi.'*[F0;0;0]];
C = [ [1 0 0]*Phi, zeros(1,n) ];
D = 0;

add_block('simulink/Sources/Sine Wave',[mdl '/Sine'], 'Amplitude','1', 'Frequency',num2str(Omega));
add_block('simulink/Continuous/State-Space',[mdl '/ModalSS']);
set_param([mdl '/ModalSS'],'A','A','B','B','C','C','D','D');
add_block('simulink/Sinks/Scope',[mdl '/Scope']);

add_line(mdl,'Sine/1','ModalSS/1');
add_line(mdl,'ModalSS/1','Scope/1');

assignin('base','A',A); assignin('base','B',B);
assignin('base','C',C); assignin('base','D',D);

set_param(mdl,'StopTime','8');
save_system(mdl);
      

8.5 Wolfram Mathematica: Eigensystem[{K,M}] + NDSolve

Library notes: Mathematica supports generalized EVP directly via Eigensystem[{K,M}].

File: Chapter13_Lesson3.nb

(* Chapter13_Lesson3.nb (plain-text Notebook expression) *)

Notebook[{
  Cell["Chapter 13, Lesson 3 — Modal Coordinates and Decoupling of MDOF Systems", "Title"],

  Cell["Build a 3-DOF mass-spring system", "Section"],
  Cell[BoxData@ToBoxes@
    HoldForm[
      m = {1.0, 1.2, 0.9};
      k = {1200., 900., 700., 1100.};
      M = DiagonalMatrix[m];
      {k1,k2,k3,k4} = k;
      K = { {k1+k2, -k2, 0.},
           {-k2, k2+k3, -k3},
           {0., -k3, k3+k4} };
    ], "Input"],

  Cell["Generalized eigenproblem: K.phi == lambda M.phi", "Section"],
  Cell[BoxData@ToBoxes@
    HoldForm[
      {vals, vecs} = Eigensystem[{K, M}];
      ord = Ordering[vals];
      vals = vals[[ord]]; vecs = vecs[[ord]];
      omega = Sqrt[Max[vals, 0]];
    ], "Input"],

  Cell["Mass-normalize modes: Phi^T M Phi = I", "Section"],
  Cell[BoxData@ToBoxes@
    HoldForm[
      Phi = Transpose[vecs];
      Phi = MapIndexed[
        Function[{v, idx}, v/Sqrt[v.M.v]],
        Transpose[Phi]
      ] // Transpose;
      Simplify[Transpose[Phi].M.Phi]
    ], "Input"],

  Cell["Decoupled modal ODE and reconstruction x = Phi q", "Section"],
  Cell[BoxData@ToBoxes@
    HoldForm[
      F0 = 10.; \[CapitalOmega] = 0.9*omega[[1]];
      f[t_] := {F0*Sin[\[CapitalOmega] t], 0., 0.};
      p[t_] := Transpose[Phi].f[t];

      q = Array[q, 3];
      eqs = Table[
        q[i]''[t] + omega[[i]]^2 q[i][t] == p[t][[i]],
        {i, 1, 3}
      ];
      ics = Join[
        Table[q[i][0] == 0, {i,1,3}],
        Table[q[i]'[0] == 0, {i,1,3}]
      ];

      sol = NDSolve[Join[eqs, ics], Table[q[i], {i,1,3}], {t,0,8}][[1]];
      qvec[t_] := Table[q[i][t], {i,1,3}] /. sol;
      xvec[t_] := Phi.qvec[t];

      Plot[Evaluate[xvec[t][[1]]], {t,0,8}, PlotRange->All,
           AxesLabel->{"t","x1"}, PlotLabel->"x1(t) via modal coordinates"]
    ], "Input"]
}]
      

Problem 1 (Orthogonality derivation): Let \( \mathbf{K}\boldsymbol{\phi}_i=\lambda_i\mathbf{M}\boldsymbol{\phi}_i \) and \( \mathbf{K}\boldsymbol{\phi}_j=\lambda_j\mathbf{M}\boldsymbol{\phi}_j \), with \( \mathbf{M}=\mathbf{M}^\top \), \( \mathbf{K}=\mathbf{K}^\top \). Show that if \( \lambda_i\neq\lambda_j \), then \( \boldsymbol{\phi}_i^\top\mathbf{M}\boldsymbol{\phi}_j=0 \).

Solution: Multiply \( \mathbf{K}\boldsymbol{\phi}_i=\lambda_i\mathbf{M}\boldsymbol{\phi}_i \) on the left by \( \boldsymbol{\phi}_j^\top \) and \( \mathbf{K}\boldsymbol{\phi}_j=\lambda_j\mathbf{M}\boldsymbol{\phi}_j \) on the left by \( \boldsymbol{\phi}_i^\top \). By symmetry, \( \boldsymbol{\phi}_j^\top\mathbf{K}\boldsymbol{\phi}_i=\boldsymbol{\phi}_i^\top\mathbf{K}\boldsymbol{\phi}_j \) and \( \boldsymbol{\phi}_j^\top\mathbf{M}\boldsymbol{\phi}_i=\boldsymbol{\phi}_i^\top\mathbf{M}\boldsymbol{\phi}_j \), leading to \( (\lambda_i-\lambda_j)\boldsymbol{\phi}_i^\top\mathbf{M}\boldsymbol{\phi}_j=0 \). Hence \( \boldsymbol{\phi}_i^\top\mathbf{M}\boldsymbol{\phi}_j=0 \). ∎

Problem 2 (Explicit 2-DOF modal decoupling): Consider \( \mathbf{M}=m\mathbf{I}_2 \) and

\[ \mathbf{K}= \begin{bmatrix} 2k & -k\\ -k & 2k \end{bmatrix}. \]

(a) Compute \( \omega_1,\omega_2 \) and mode shapes. (b) Choose a mass-normalized modal matrix \( \mathbf{\Phi} \). (c) Show that \( \mathbf{\Phi}^\top\mathbf{K}\mathbf{\Phi}=\mathbf{\Omega}^2 \).

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

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

Thus \( \omega_1^2=\frac{k}{m} \) (plus sign) and \( \omega_2^2=\frac{3k}{m} \) (minus sign). Associated eigenvectors are \( \boldsymbol{\phi}_1 \propto [1\; 1]^\top \), \( \boldsymbol{\phi}_2 \propto [1\; -1]^\top \). Mass-normalize: for \( \mathbf{M}=m\mathbf{I} \), enforce \( \boldsymbol{\phi}_i^\top m\mathbf{I}\boldsymbol{\phi}_i=1 \). Taking

\[ \boldsymbol{\phi}_1=\frac{1}{\sqrt{2m} }\begin{bmatrix}1\\1\end{bmatrix},\quad \boldsymbol{\phi}_2=\frac{1}{\sqrt{2m} }\begin{bmatrix}1\\-1\end{bmatrix}, \quad \mathbf{\Phi}=[\boldsymbol{\phi}_1\;\boldsymbol{\phi}_2], \]

one verifies directly that \( \mathbf{\Phi}^\top\mathbf{M}\mathbf{\Phi}=\mathbf{I} \) and \( \mathbf{\Phi}^\top\mathbf{K}\mathbf{\Phi}=\mathrm{diag}(\omega_1^2,\omega_2^2) \), hence decoupling follows from Section 5. ∎

Problem 3 (Modal harmonic steady-state): Suppose for some mode \( i \) you have \( \ddot{q}_i+\omega_i^2 q_i=\hat{p}_i\sin(\Omega t) \) with \( \Omega \neq \omega_i \). Find a sinusoidal steady-state solution.

Solution: Try \( q_i^{ss}(t)=A\sin(\Omega t) \). Then \( \ddot{q}_i^{ss}=-A\Omega^2\sin(\Omega t) \), and substitution gives \( (-A\Omega^2+\omega_i^2 A)\sin(\Omega t)=\hat{p}_i\sin(\Omega t) \). Therefore \( A=\frac{\hat{p}_i}{\omega_i^2-\Omega^2} \), so \( q_i^{ss}(t)=\frac{\hat{p}_i}{\omega_i^2-\Omega^2}\sin(\Omega t) \). ∎

Problem 4 (Energy separation): Assume mass normalization \( \mathbf{\Phi}^\top\mathbf{M}\mathbf{\Phi}=\mathbf{I} \) and \( \mathbf{\Phi}^\top\mathbf{K}\mathbf{\Phi}=\mathbf{\Omega}^2 \). Show that the total energy \( T+V \) equals \( \tfrac{1}{2}\sum_i \dot{q}_i^2 + \tfrac{1}{2}\sum_i \omega_i^2 q_i^2 \).

Solution: Substitute \( \mathbf{x}=\mathbf{\Phi}\mathbf{q} \) and \( \dot{\mathbf{x} }=\mathbf{\Phi}\dot{\mathbf{q} } \). Then

\[ T=\tfrac{1}{2}\dot{\mathbf{q} }^\top(\mathbf{\Phi}^\top\mathbf{M}\mathbf{\Phi})\dot{\mathbf{q} } =\tfrac{1}{2}\dot{\mathbf{q} }^\top\dot{\mathbf{q} }=\tfrac{1}{2}\sum_i \dot{q}_i^2, \]

\[ V=\tfrac{1}{2}\mathbf{q}^\top(\mathbf{\Phi}^\top\mathbf{K}\mathbf{\Phi})\mathbf{q} =\tfrac{1}{2}\mathbf{q}^\top\mathbf{\Omega}^2\mathbf{q}=\tfrac{1}{2}\sum_i \omega_i^2 q_i^2. \]

Hence \( T+V \) is the sum of independent modal energies. ∎

Problem 5 (Reduced modal model via projection): Let \( \mathbf{\Phi}_r=[\boldsymbol{\phi}_1\cdots\boldsymbol{\phi}_r] \) contain the first \( r \) mass-normalized modes (\( r<n \)) and approximate \( \mathbf{x}\approx \mathbf{\Phi}_r\mathbf{q}_r \). Derive the reduced equations governing \( \mathbf{q}_r(t) \).

Solution: Substitute \( \mathbf{x}\approx\mathbf{\Phi}_r\mathbf{q}_r \) into \( \mathbf{M}\ddot{\mathbf{x} }+\mathbf{K}\mathbf{x}=\mathbf{f}(t) \) and left-multiply by \( \mathbf{\Phi}_r^\top \):

\[ \mathbf{\Phi}_r^\top\mathbf{M}\mathbf{\Phi}_r\,\ddot{\mathbf{q} }_r + \mathbf{\Phi}_r^\top\mathbf{K}\mathbf{\Phi}_r\,\mathbf{q}_r = \mathbf{\Phi}_r^\top\mathbf{f}(t). \]

Because the retained modes are mass-normalized and mutually \( \mathbf{M} \)-orthonormal, \( \mathbf{\Phi}_r^\top\mathbf{M}\mathbf{\Phi}_r=\mathbf{I}_r \), and \( \mathbf{\Phi}_r^\top\mathbf{K}\mathbf{\Phi}_r=\mathrm{diag}(\omega_1^2,\dots,\omega_r^2) \). Thus the reduced model is: \( \ddot{\mathbf{q} }_r + \mathbf{\Omega}_r^2 \mathbf{q}_r = \mathbf{\Phi}_r^\top\mathbf{f}(t) \). Reconstruction is \( \mathbf{x}(t)\approx \mathbf{\Phi}_r\mathbf{q}_r(t) \). ∎

We proved the \( \mathbf{M} \)-orthogonality of normal modes for symmetric undamped MDOF systems, introduced mass normalization \( \mathbf{\Phi}^\top\mathbf{M}\mathbf{\Phi}=\mathbf{I} \), and derived the exact decoupling \( \ddot{\mathbf{q} }+\mathbf{\Omega}^2\mathbf{q}=\mathbf{\Phi}^\top\mathbf{f}(t) \). The resulting scalar modal ODEs provide a mathematically clean lens for analysis and simulation. In the next lesson, we will study how damping alters (or preserves) this decoupling.