The standard linear viscously-damped MDOF model is \( \mathbf{x}(t)\in\mathbb{R}^n \): the displacement vector, \( \mathbf{M}\in\mathbb{R}^{n\times n} \): symmetric positive definite mass matrix, \( \mathbf{K}\in\mathbb{R}^{n\times n} \): symmetric stiffness matrix, and \( \mathbf{C}\in\mathbb{R}^{n\times n} \): viscous damping matrix.

\[ \mathbf{M}\ddot{\mathbf{x} }(t) + \mathbf{C}\dot{\mathbf{x} }(t) + \mathbf{K}\mathbf{x}(t) = \mathbf{f}(t). \]

The key structural question is whether the same real modal basis that diagonalizes \( \mathbf{M} \) and \( \mathbf{K} \) also diagonalizes \( \mathbf{C} \). This determines whether the system decouples into independent SDOF modal equations.

flowchart TD
  A["Start: M xdd + C xd + K x = f(t)"] --> B["Compute undamped modes: K phi = (w^2) M phi"]
  B --> C["Form modal transform: x = Phi q"]
  C --> D["Check modal damping matrix: Cm = Phi^T C Phi"]
  D -->|"Cm is diagonal (or approx. diagonal)"| E["Classical / proportional damping: \nreal modes decouple"]
  D -->|Cm has strong off-diagonals| F["Non-classical damping: \nuse state-space, complex modes"]
  

In practice, classical damping is a modeling assumption that enables tractable analysis and controller-oriented reduced models; non-classical damping is common when damping devices are localized (dashpots, joints, interfaces), or when multiple dissipation mechanisms combine.

Recall from Lessons 2–3 the undamped generalized eigenproblem and modal matrix: \( \mathbf{K}\boldsymbol{\phi}_r = \omega_r^2 \mathbf{M}\boldsymbol{\phi}_r \): \( \boldsymbol{\Phi}=[\boldsymbol{\phi}_1,\dots,\boldsymbol{\phi}_n] \):. Under standard assumptions ( \( \mathbf{M}\succ 0 \), \( \mathbf{K}\succeq 0 \)), choose mass-normalized modes:

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

Use the modal coordinate transformation \( \mathbf{x}=\boldsymbol{\Phi}\mathbf{q} \):. Substituting into the damped equation and premultiplying by \( \boldsymbol{\Phi}^T \) gives:

\[ \ddot{\mathbf{q} } + \underbrace{\left(\boldsymbol{\Phi}^T\mathbf{C}\boldsymbol{\Phi}\right)}_{\mathbf{C}_m} \dot{\mathbf{q} } + \boldsymbol{\Omega}^2 \mathbf{q} = \mathbf{p}(t), \qquad \mathbf{p}(t)=\boldsymbol{\Phi}^T\mathbf{f}(t). \]

Therefore, modal decoupling holds iff \( \mathbf{C}_m=\boldsymbol{\Phi}^T\mathbf{C}\boldsymbol{\Phi} \) is diagonal (or sufficiently close to diagonal).

A widely used sufficient condition is Rayleigh (proportional) damping:

\[ \mathbf{C} = \alpha \mathbf{M} + \beta \mathbf{K}, \]

with scalars \( \alpha\ge 0 \) and \( \beta\ge 0 \). Then

\[ \mathbf{C}_m = \boldsymbol{\Phi}^T(\alpha\mathbf{M}+\beta\mathbf{K})\boldsymbol{\Phi} = \alpha(\boldsymbol{\Phi}^T\mathbf{M}\boldsymbol{\Phi}) + \beta(\boldsymbol{\Phi}^T\mathbf{K}\boldsymbol{\Phi}) = \alpha \mathbf{I} + \beta \boldsymbol{\Omega}^2, \]

which is diagonal. Hence each modal coordinate satisfies an independent damped SDOF equation:

\[ \ddot{q}_r + (\alpha + \beta \omega_r^2)\dot{q}_r + \omega_r^2 q_r = p_r(t). \]

Writing the standard second-order form \( \ddot{q}_r + 2\zeta_r\omega_r \dot{q}_r + \omega_r^2 q_r = p_r(t) \), we identify the modal damping ratio:

\[ \zeta_r = \frac{\alpha}{2\omega_r} + \frac{\beta \omega_r}{2}. \]

Mode shapes with proportional damping. The mode shapes remain the undamped real vectors \( \boldsymbol{\phi}_r \); damping only changes the modal poles (eigenvalues), not the real modal basis.

For free vibration in the r-th decoupled mode ( \( p_r(t)=0 \)), assume \( q_r(t)=e^{\lambda t} \). Substitution yields:

\[ \lambda^2 + 2\zeta_r\omega_r \lambda + \omega_r^2 = 0. \]

The roots are the modal poles:

\[ \lambda_{r,\pm} = -\zeta_r\omega_r \pm j\,\omega_r\sqrt{1-\zeta_r^2}, \qquad 0 \,\, < \,\, \zeta_r \,\, < \,\, 1. \]

Define \( \omega_{d,r} \): damped natural frequency:

\[ \omega_{d,r} = \omega_r\sqrt{1-\zeta_r^2}. \]

The modal time response is:

\[ q_r(t) = e^{-\zeta_r\omega_r t}\left(A_r\cos(\omega_{d,r} t) + B_r\sin(\omega_{d,r} t)\right). \]

Energy interpretation. In a viscously damped mode, the instantaneous power dissipated is \( P_r(t)=2\zeta_r\omega_r \dot{q}_r^2 \) (for mass-normalized modal coordinate), so damping ratio quantifies per-cycle energy loss in the modal subsystem.

When \( \mathbf{C} \) is not proportional (or not simultaneously diagonalizable with \( \mathbf{M},\mathbf{K} \)), the modal damping matrix \( \mathbf{C}_m=\boldsymbol{\Phi}^T\mathbf{C}\boldsymbol{\Phi} \) contains off-diagonal terms that couple modal velocities. Then the real undamped modes do not yield independent SDOF equations.

For free vibration, a direct exponential trial \( \mathbf{x}(t)=\boldsymbol{\phi}e^{\lambda t} \) gives the quadratic eigenvalue problem (QEP):

\[ \left(\lambda^2\mathbf{M} + \lambda\mathbf{C} + \mathbf{K}\right)\boldsymbol{\phi}=\mathbf{0}, \qquad (\boldsymbol{\phi}\neq \mathbf{0}). \]

In general, the eigenvalues \( \lambda \) are complex and the corresponding eigenvectors \( \boldsymbol{\phi} \) are complex (complex mode shapes). These complex modes capture both spatial deformation and phase differences between DOFs induced by damping.

A standard linearization converts the second-order system to a first-order state-space system. Let \( \mathbf{z}=\begin{bmatrix}\mathbf{x}\\\dot{\mathbf{x} }\end{bmatrix}\in\mathbb{R}^{2n} \). Then:

\[ \dot{\mathbf{z} } = \mathbf{A}\mathbf{z}, \qquad \mathbf{A} = \begin{bmatrix} \mathbf{0} & \mathbf{I}\\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{C} \end{bmatrix}. \]

The state-space eigenproblem \( \mathbf{A}\mathbf{v}=\lambda \mathbf{v} \) produces the same eigenvalues as the QEP. Writing \( \lambda=\sigma + j\omega_d \), we define:

\[ \omega_n = \sqrt{\sigma^2 + \omega_d^2}, \qquad \zeta = -\frac{\sigma}{\omega_n}. \]

Mode shapes with non-proportional damping. The displacement component of \( \mathbf{v}=\begin{bmatrix}\boldsymbol{\phi}\\\lambda\boldsymbol{\phi}\end{bmatrix} \) is complex, so the physical motion is recovered from the real part of a modal superposition. This is the minimum machinery needed for the next lesson on experimental modal analysis, where measured FRFs naturally exhibit complex modes.

In computational work, damping modeling typically follows a staged workflow: first compute undamped modes, then assess whether proportional damping is a defensible approximation, otherwise move to state-space complex modes.

flowchart TD
  S["Inputs: M, C, K"] --> U["Undamped eig: K Phi = (w^2) M Phi"]
  U --> N["Normalize: Phi^T M Phi = I"]
  N --> T["Compute Cm = Phi^T C Phi"]
  T --> Q["Is Cm nearly diagonal?"]
  Q -->|yes| P["Use modal SDOF set: \nqdd + 2 zeta w qd + w^2 q = p(t)"]
  Q -->|no| X["Build state matrix A \nand compute complex eig"]
  X --> R["Extract zeta, wn, wd \nfrom lambda = sigma + j wd"]
  R --> M1["Construct responses via \nstate-space or complex modal expansion"]
        

“Nearly diagonal” depends on the objective: for low-order control models, small off-diagonal terms may be ignored; for accurate prediction near closely-spaced modes, non-proportional damping can be decisive.

File: Chapter13_Lesson4.py

"""
Chapter13_Lesson4.py
Damping in MDOF Systems and Mode Shapes with Damping

This script demonstrates:
1) Undamped modal analysis: K phi = (w^2) M phi, with mass-normalization.
2) Rayleigh (proportional) damping: C = alpha*M + beta*K, and modal damping ratios.
3) Non-proportional damping: state-space eigenanalysis -> complex poles and complex mode shapes.
4) Free response simulation for the non-proportional case.

Dependencies:
  pip install numpy scipy matplotlib
"""

import numpy as np

try:
    from scipy.linalg import eigh, eig
    from scipy.integrate import solve_ivp
except ImportError as e:
    raise SystemExit("SciPy is required. Install with: pip install scipy") from e


def chain_matrices_3dof(m1=1.0, m2=1.0, m3=1.0, k1=2000.0, k2=1500.0, k3=1000.0):
    """3-DOF shear/chain model with base fixed.
    DOF i is the displacement of mass i.
    Springs: k1 (ground-m1), k2 (m1-m2), k3 (m2-m3).
    """
    M = np.diag([m1, m2, m3])

    K = np.array([
        [k1 + k2, -k2,      0.0],
        [-k2,     k2 + k3, -k3 ],
        [0.0,     -k3,      k3 ]
    ], dtype=float)
    return M, K


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


def rayleigh_from_two_targets(omega_i, zeta_i, omega_j, zeta_j):
    """Solve for Rayleigh coefficients alpha, beta using:
         zeta_r = alpha/(2*omega_r) + beta*omega_r/2
    at two frequencies (omega_i, omega_j).
    """
    # alpha + beta*omega^2 = 2*zeta*omega
    A = np.array([[1.0, omega_i**2],
                  [1.0, omega_j**2]], dtype=float)
    b = np.array([2.0*zeta_i*omega_i, 2.0*zeta_j*omega_j], dtype=float)
    alpha, beta = np.linalg.solve(A, b)
    return alpha, beta


def build_state_matrix(M, C, K):
    """Build first-order state matrix A for:
       M xdd + C xd + K x = 0
       z = [x; xd], z_dot = A z
    """
    n = M.shape[0]
    Minv = np.linalg.inv(M)
    Z = np.zeros((n, n))
    I = np.eye(n)
    A = np.block([
        [Z, I],
        [-Minv @ K, -Minv @ C]
    ])
    return A


def modal_damping_matrix(Phi, C):
    return Phi.T @ C @ Phi


def summarize_poles(lam, max_modes=6):
    """Return a list of (sigma, wd, wn, zeta) for conjugate pairs with wd>0."""
    out = []
    # keep only positive imaginary part poles (each mode once)
    lam_pos = [x for x in lam if np.imag(x) > 1e-8]
    lam_pos.sort(key=lambda x: np.imag(x))
    for x in lam_pos[:max_modes]:
        sigma = np.real(x)
        wd = np.imag(x)
        wn = np.sqrt(sigma**2 + wd**2)
        zeta = -sigma / wn if wn > 0 else np.nan
        out.append((sigma, wd, wn, zeta))
    return out


def main():
    np.set_printoptions(precision=5, suppress=True)

    # --------------------------
    # Example system (3-DOF)
    # --------------------------
    M, K = chain_matrices_3dof(m1=1.2, m2=1.0, m3=0.8, k1=2500.0, k2=1800.0, k3=1200.0)
    n = M.shape[0]

    # --------------------------
    # (A) Undamped modal analysis
    # --------------------------
    # For symmetric M, K with M positive definite:
    # scipy.linalg.eigh solves K v = lam M v
    lam, Phi = eigh(K, M)  # lam = omega^2
    omega = np.sqrt(lam)
    Phi = mass_normalize_modes(M, Phi)

    print("Undamped natural frequencies (rad/s):", omega)
    print("Check Phi^T M Phi (should be I):\n", Phi.T @ M @ Phi)
    print("Check Phi^T K Phi (should be diag(omega^2)):\n", Phi.T @ K @ Phi)

    # ------------------------------------
    # (B) Proportional damping (Rayleigh)
    # ------------------------------------
    # Target damping ratios at two modes (typical structural values)
    zeta_1 = 0.02  # 2% at mode 1
    zeta_3 = 0.05  # 5% at mode 3
    alpha, beta = rayleigh_from_two_targets(omega[0], zeta_1, omega[2], zeta_3)

    C_ray = alpha*M + beta*K
    Cm_ray = modal_damping_matrix(Phi, C_ray)

    print("\nRayleigh damping coefficients:")
    print("alpha =", alpha, " beta =", beta)
    print("Modal damping matrix Phi^T C Phi (Rayleigh) ~ diagonal:\n", Cm_ray)

    # Modal damping ratios predicted by Rayleigh:
    zeta_modal = 0.5*(alpha/omega + beta*omega)
    print("Modal damping ratios (Rayleigh):", zeta_modal)

    # ------------------------------------------
    # (C) Non-proportional damping (example)
    # ------------------------------------------
    # Build a localized damping that is not representable as alpha*M + beta*K:
    # Dashpot between DOF1 and DOF3 (non-adjacent coupling) + dashpot to ground at DOF2.
    c13 = 45.0
    c2g = 35.0
    C_np = np.zeros((n, n), dtype=float)
    # between 1 and 3
    C_np[0, 0] += c13
    C_np[2, 2] += c13
    C_np[0, 2] -= c13
    C_np[2, 0] -= c13
    # DOF2 to ground
    C_np[1, 1] += c2g

    # Compare modal damping matrix (should have off-diagonals)
    Cm_np = modal_damping_matrix(Phi, C_np)
    print("\nNon-proportional C (example):\n", C_np)
    print("Modal damping matrix Phi^T C Phi (non-proportional) has off-diagonals:\n", Cm_np)

    # State-space eigenanalysis => complex poles and complex mode shapes
    A = build_state_matrix(M, C_np, K)
    lamA, vA = eig(A)

    print("\nComplex poles (sigma + j*wd) for first few modes (wd>0):")
    poles = summarize_poles(lamA, max_modes=3)
    for i, (sigma, wd, wn, zeta) in enumerate(poles, start=1):
        print(f"Mode {i}: sigma={sigma:+.6f}, wd={wd:.6f}, wn={wn:.6f}, zeta={zeta:.6f}")

    # Extract one complex mode shape (displacement part) corresponding to smallest wd>0
    # Find index of pole with smallest positive imaginary part
    idx_candidates = np.where(np.imag(lamA) > 1e-8)[0]
    idx = idx_candidates[np.argmin(np.imag(lamA[idx_candidates]))]
    lam_sel = lamA[idx]
    v_sel = vA[:, idx]
    phi_c = v_sel[:n]  # displacement component
    # normalize by maximum magnitude for display
    phi_c = phi_c / np.max(np.abs(phi_c))

    print("\nExample complex mode shape (normalized, displacement part):")
    print("lambda =", lam_sel)
    print("phi_complex =", phi_c)

    # ------------------------------------------
    # (D) Free response simulation (non-prop)
    # ------------------------------------------
    def zdot(t, z):
        x = z[:n]
        xd = z[n:]
        xdd = np.linalg.solve(M, -C_np @ xd - K @ x)
        return np.hstack([xd, xdd])

    # Initial condition: small displacement in DOF1, rest zero
    x0 = np.array([0.01, 0.0, 0.0])
    v0 = np.zeros(n)
    z0 = np.hstack([x0, v0])

    t_end = 5.0
    sol = solve_ivp(zdot, (0.0, t_end), z0, max_step=1e-2, rtol=1e-8, atol=1e-10)

    # Optional: print a few samples
    print("\nFree response simulation completed.")
    print("t grid size:", sol.t.size)
    print("x(t_end) =", sol.y[:n, -1])

    # If matplotlib is present, plot displacements
    try:
        import matplotlib.pyplot as plt
        plt.figure()
        for i in range(n):
            plt.plot(sol.t, sol.y[i, :], label=f"x{i+1}(t)")
        plt.xlabel("t [s]")
        plt.ylabel("displacement")
        plt.title("Non-proportional damping free response (3-DOF)")
        plt.legend()
        plt.grid(True)
        plt.show()
    except Exception:
        pass


if __name__ == "__main__":
    main()
      

The script computes (i) undamped modes, (ii) Rayleigh damping coefficients from two target modal damping ratios, (iii) a non-proportional damping case with complex modes from state-space eigenanalysis, and (iv) a free-response simulation.

File: Chapter13_Lesson4.cpp

/*
Chapter13_Lesson4.cpp
Damping in MDOF Systems and Mode Shapes with Damping

Demonstrates:
1) Undamped generalized eigenproblem: K phi = (w^2) M phi
2) Rayleigh damping: C = alpha*M + beta*K -> diagonal modal damping matrix
3) Non-proportional damping: state-space eigenanalysis -> complex poles

Dependencies:
  - Eigen (https://eigen.tuxfamily.org/)
Build (example):
  g++ -O2 -std=c++17 Chapter13_Lesson4.cpp -I /path/to/eigen -o Chapter13_Lesson4
*/

#include <iostream>
#include <vector>
#include <complex>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>

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

static void chain_matrices_3dof(MatrixXd& M, MatrixXd& K) {
    const double m1 = 1.2, m2 = 1.0, m3 = 0.8;
    const double k1 = 2500.0, k2 = 1800.0, k3 = 1200.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;
}

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

static void rayleigh_from_two_targets(double omega_i, double zeta_i,
                                      double omega_j, double zeta_j,
                                      double& alpha, double& beta) {
    // alpha + beta*omega^2 = 2*zeta*omega at two omegas
    Eigen::Matrix2d A;
    A << 1.0, omega_i*omega_i,
         1.0, omega_j*omega_j;
    Eigen::Vector2d b;
    b << 2.0*zeta_i*omega_i,
         2.0*zeta_j*omega_j;
    Eigen::Vector2d x = A.fullPivLu().solve(b);
    alpha = x(0);
    beta  = x(1);
}

static MatrixXd build_state_matrix(const MatrixXd& M, const MatrixXd& C, const MatrixXd& K) {
    const int n = static_cast<int>(M.rows());
    MatrixXd Minv = M.inverse();

    MatrixXd A = MatrixXd::Zero(2*n, 2*n);
    A.block(0, n, n, n) = MatrixXd::Identity(n, n);
    A.block(n, 0, n, n) = -Minv * K;
    A.block(n, n, n, n) = -Minv * C;
    return A;
}

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

    // Undamped generalized eigenproblem (symmetric):
    // Use K' = L^{-1} K L^{-T} with M = L L^T (Cholesky), then eig(K')
    Eigen::LLT<MatrixXd> llt(M);
    MatrixXd L = llt.matrixL();
    MatrixXd Linv = L.inverse();
    MatrixXd Kt = Linv * K * Linv.transpose(); // symmetric

    Eigen::SelfAdjointEigenSolver<MatrixXd> es(Kt);
    VectorXd lam = es.eigenvalues();           // omega^2
    MatrixXd Y = es.eigenvectors();            // in transformed coordinates
    VectorXd omega = lam.array().sqrt();

    // Back-transform mode shapes: phi = L^{-T} y
    MatrixXd Phi = L.transpose().inverse() * Y;
    Phi = mass_normalize_modes(M, Phi);

    std::cout << "Undamped natural frequencies (rad/s):\n" << omega.transpose() << "\n\n";
    std::cout << "Phi^T M Phi (should be I):\n" << (Phi.transpose()*M*Phi) << "\n\n";
    std::cout << "Phi^T K Phi (should be diag(omega^2)):\n" << (Phi.transpose()*K*Phi) << "\n\n";

    // Rayleigh damping parameters from two target zetas (mode 1 and 3)
    double alpha=0.0, beta=0.0;
    rayleigh_from_two_targets(omega(0), 0.02, omega(2), 0.05, alpha, beta);
    MatrixXd C_ray = alpha*M + beta*K;
    MatrixXd Cm_ray = Phi.transpose() * C_ray * Phi;

    std::cout << "Rayleigh coefficients: alpha=" << alpha << ", beta=" << beta << "\n";
    std::cout << "Modal damping matrix Phi^T C Phi (Rayleigh) ~ diagonal:\n" << Cm_ray << "\n\n";

    // Modal damping ratios: zeta_r = 0.5*(alpha/omega_r + beta*omega_r)
    VectorXd zeta = 0.5*(alpha*omega.cwiseInverse() + beta*omega);
    std::cout << "Modal damping ratios (Rayleigh):\n" << zeta.transpose() << "\n\n";

    // Non-proportional damping example: damper between DOF1 and DOF3, plus damper to ground at DOF2
    MatrixXd C_np = MatrixXd::Zero(3,3);
    double c13 = 45.0;
    double c2g = 35.0;
    C_np(0,0) += c13; C_np(2,2) += c13;
    C_np(0,2) -= c13; C_np(2,0) -= c13;
    C_np(1,1) += c2g;

    MatrixXd Cm_np = Phi.transpose() * C_np * Phi;
    std::cout << "Non-proportional C:\n" << C_np << "\n\n";
    std::cout << "Phi^T C Phi (non-proportional) has off-diagonals:\n" << Cm_np << "\n\n";

    // State-space eigenanalysis for complex poles
    MatrixXd A = build_state_matrix(M, C_np, K);
    Eigen::EigenSolver<MatrixXd> ces(A);
    Eigen::VectorXcd evals = ces.eigenvalues();

    // Print a few poles with positive imaginary part (one per conjugate pair)
    std::vector<std::complex<double>> poles;
    poles.reserve(evals.size());
    for (int i=0; i<evals.size(); ++i) {
        if (std::imag(evals(i)) > 1e-8) poles.push_back(evals(i));
    }
    std::sort(poles.begin(), poles.end(),
              [](auto a, auto b){ return std::imag(a) < std::imag(b); });

    std::cout << "Complex poles (sigma + j*wd), first three modes:\n";
    for (int i=0; i<std::min<int>(3, (int)poles.size()); ++i) {
        double sigma = std::real(poles[i]);
        double wd = std::imag(poles[i]);
        double wn = std::sqrt(sigma*sigma + wd*wd);
        double zeta_np = (wn > 0.0) ? (-sigma/wn) : 0.0;
        std::cout << "Mode " << (i+1) << ": sigma=" << sigma << ", wd=" << wd
                  << ", wn=" << wn << ", zeta=" << zeta_np << "\n";
    }

    return 0;
}
      

This example uses a Cholesky-based symmetric transformation for the generalized eigenproblem and a state-space eigenanalysis for complex poles under non-proportional damping.

File: Chapter13_Lesson4.java

/*
Chapter13_Lesson4.java
Damping in MDOF Systems and Mode Shapes with Damping

Demonstrates:
1) Undamped generalized eigenproblem (symmetric) via Cholesky transform:
      K phi = (w^2) M phi
   with M = L L^T, transform to A = L^{-1} K L^{-T}, solve A y = (w^2) y, then phi = L^{-T} y.
2) Rayleigh damping coefficients from two target modal damping ratios.
3) Non-proportional damping: state-space eigenanalysis to obtain complex poles.

Library:
  EJML (Efficient Java Matrix Library)
  Maven coordinates (example):
    org.ejml:ejml-all:0.43

This is a single-file demo for teaching; in production, organize into packages.
*/

import org.ejml.data.DMatrixRMaj;
import org.ejml.data.Complex_F64;
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.ArrayList;
import java.util.Comparator;

public class Chapter13_Lesson4 {

    static class PoleInfo {
        double sigma, wd, wn, zeta;
        PoleInfo(double sigma, double wd, double wn, double zeta) {
            this.sigma = sigma; this.wd = wd; this.wn = wn; this.zeta = zeta;
        }
    }

    static void chainMatrices3DOF(DMatrixRMaj M, DMatrixRMaj K) {
        double m1 = 1.2, m2 = 1.0, m3 = 0.8;
        double k1 = 2500.0, k2 = 1800.0, k3 = 1200.0;

        M.reshape(3,3);
        K.reshape(3,3);
        CommonOps_DDRM.fill(M, 0.0);
        CommonOps_DDRM.fill(K, 0.0);

        M.set(0,0, m1); M.set(1,1, m2); M.set(2,2, m3);

        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);
    }

    static DMatrixRMaj massNormalizeModes(DMatrixRMaj M, DMatrixRMaj Phi) {
        int n = Phi.getNumRows();
        int r = Phi.getNumCols();
        DMatrixRMaj PhiN = Phi.copy();
        DMatrixRMaj tmp = new DMatrixRMaj(n,1);
        for (int i = 0; i < r; i++) {
            // mi = phi_i^T M phi_i
            CommonOps_DDRM.extract(PhiN, 0, n, i, i+1, tmp, 0, 0);
            DMatrixRMaj Mphi = new DMatrixRMaj(n,1);
            CommonOps_DDRM.mult(M, tmp, Mphi);
            double mi = CommonOps_DDRM.dot(tmp, Mphi);
            double scale = 1.0/Math.sqrt(mi);
            for (int k=0; k<n; k++) {
                PhiN.set(k, i, PhiN.get(k, i)*scale);
            }
        }
        return PhiN;
    }

    static double[] rayleighFromTwoTargets(double omegaI, double zetaI, double omegaJ, double zetaJ) {
        // alpha + beta*omega^2 = 2*zeta*omega
        double A11 = 1.0, A12 = omegaI*omegaI;
        double A21 = 1.0, A22 = omegaJ*omegaJ;
        double b1 = 2.0*zetaI*omegaI;
        double b2 = 2.0*zetaJ*omegaJ;

        double det = A11*A22 - A12*A21;
        double alpha = ( b1*A22 - A12*b2 )/det;
        double beta  = ( A11*b2 - b1*A21 )/det;
        return new double[]{alpha, beta};
    }

    static DMatrixRMaj buildStateMatrix(DMatrixRMaj M, DMatrixRMaj C, DMatrixRMaj K) {
        int n = M.getNumRows();
        DMatrixRMaj Minv = new DMatrixRMaj(n,n);
        CommonOps_DDRM.invert(M, Minv);

        DMatrixRMaj A = new DMatrixRMaj(2*n, 2*n);
        CommonOps_DDRM.fill(A, 0.0);

        // top-right I
        for (int i=0; i<n; i++) A.set(i, n+i, 1.0);

        // bottom-left -Minv*K
        DMatrixRMaj MinvK = new DMatrixRMaj(n,n);
        CommonOps_DDRM.mult(Minv, K, MinvK);
        CommonOps_DDRM.scale(-1.0, MinvK);
        CommonOps_DDRM.insert(MinvK, A, n, 0);

        // bottom-right -Minv*C
        DMatrixRMaj MinvC = new DMatrixRMaj(n,n);
        CommonOps_DDRM.mult(Minv, C, MinvC);
        CommonOps_DDRM.scale(-1.0, MinvC);
        CommonOps_DDRM.insert(MinvC, A, n, n);

        return A;
    }

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

        // Cholesky of M: M = L L^T
        CholeskyDecomposition_F64<DMatrixRMaj> chol = DecompositionFactory_DDRM.chol(M.numRows, true);
        if (!chol.decompose(M.copy())) {
            throw new RuntimeException("Cholesky failed; M must be SPD.");
        }
        DMatrixRMaj L = chol.getT(null); // lower triangular

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

        // Transform A = Linv * K * Linv^T (symmetric)
        DMatrixRMaj tmp = new DMatrixRMaj(3,3);
        DMatrixRMaj A = new DMatrixRMaj(3,3);
        CommonOps_DDRM.mult(Linv, K, tmp);
        DMatrixRMaj LinvT = new DMatrixRMaj(3,3);
        CommonOps_DDRM.transpose(Linv, LinvT);
        CommonOps_DDRM.mult(tmp, LinvT, A);

        // Eigen decomposition of symmetric A -> omega^2
        EigenDecomposition_F64<DMatrixRMaj> eigSym = DecompositionFactory_DDRM.eig(3, true);
        if (!eigSym.decompose(A)) {
            throw new RuntimeException("Eigen decomposition failed.");
        }

        double[] omega = new double[3];
        DMatrixRMaj Y = new DMatrixRMaj(3,3); // eigenvectors
        for (int i=0; i<3; i++) {
            Complex_F64 ev = eigSym.getEigenvalue(i);
            omega[i] = Math.sqrt(ev.getReal());
            DMatrixRMaj vi = eigSym.getEigenVector(i);
            for (int r=0; r<3; r++) Y.set(r, i, vi.get(r,0));
        }

        // Back-transform modes: Phi = L^{-T} Y
        DMatrixRMaj LT = new DMatrixRMaj(3,3);
        CommonOps_DDRM.transpose(L, LT);
        DMatrixRMaj LTinv = new DMatrixRMaj(3,3);
        CommonOps_DDRM.invert(LT, LTinv);
        DMatrixRMaj Phi = new DMatrixRMaj(3,3);
        CommonOps_DDRM.mult(LTinv, Y, Phi);

        Phi = massNormalizeModes(M, Phi);

        System.out.println("Undamped omegas (rad/s):");
        for (double w : omega) System.out.printf("%.6f  ", w);
        System.out.println("\n");

        // Rayleigh damping from mode 1 and mode 3 targets
        double[] ab = rayleighFromTwoTargets(omega[0], 0.02, omega[2], 0.05);
        double alpha = ab[0], beta = ab[1];

        DMatrixRMaj C_ray = new DMatrixRMaj(3,3);
        DMatrixRMaj alphaM = M.copy();
        DMatrixRMaj betaK  = K.copy();
        CommonOps_DDRM.scale(alpha, alphaM);
        CommonOps_DDRM.scale(beta, betaK);
        CommonOps_DDRM.add(alphaM, betaK, C_ray);

        // Cm = Phi^T C Phi
        DMatrixRMaj PhiT = new DMatrixRMaj(3,3);
        CommonOps_DDRM.transpose(Phi, PhiT);
        DMatrixRMaj Cm = new DMatrixRMaj(3,3);
        DMatrixRMaj tmp2 = new DMatrixRMaj(3,3);
        CommonOps_DDRM.mult(PhiT, C_ray, tmp2);
        CommonOps_DDRM.mult(tmp2, Phi, Cm);

        System.out.printf("Rayleigh alpha=%.6f, beta=%.6f%n", alpha, beta);
        System.out.println("Modal damping matrix (Rayleigh), Phi^T C Phi:");
        Cm.print();

        // Non-proportional damping example
        DMatrixRMaj C_np = new DMatrixRMaj(3,3);
        CommonOps_DDRM.fill(C_np, 0.0);
        double c13 = 45.0, c2g = 35.0;
        C_np.set(0,0, C_np.get(0,0)+c13);
        C_np.set(2,2, C_np.get(2,2)+c13);
        C_np.set(0,2, C_np.get(0,2)-c13);
        C_np.set(2,0, C_np.get(2,0)-c13);
        C_np.set(1,1, C_np.get(1,1)+c2g);

        // Build state matrix and compute eigenvalues (complex poles)
        DMatrixRMaj Astate = buildStateMatrix(M, C_np, K);
        EigenDecomposition_F64<DMatrixRMaj> eig = DecompositionFactory_DDRM.eig(Astate.getNumRows(), false);
        if (!eig.decompose(Astate)) throw new RuntimeException("State eigen decomposition failed.");

        ArrayList<Complex_F64> poles = new ArrayList<>();
        for (int i=0; i<eig.getNumberOfEigenvalues(); i++) {
            Complex_F64 ev = eig.getEigenvalue(i);
            if (ev.getImaginary() > 1e-8) poles.add(ev);
        }
        poles.sort(Comparator.comparingDouble(Complex_F64::getImaginary));

        System.out.println("Complex poles (sigma + j*wd), first three modes:");
        for (int i=0; i<Math.min(3, poles.size()); i++) {
            double sigma = poles.get(i).getReal();
            double wd    = poles.get(i).getImaginary();
            double wn    = Math.sqrt(sigma*sigma + wd*wd);
            double zeta  = (wn > 0) ? (-sigma/wn) : 0.0;
            System.out.printf("Mode %d: sigma=%+.6f, wd=%.6f, wn=%.6f, zeta=%.6f%n",
                    (i+1), sigma, wd, wn, zeta);
        }
    }
}
      

The Java workflow uses EJML. For the undamped generalized eigenproblem, we apply the symmetric Cholesky transform \( \mathbf{A}=\mathbf{L}^{-1}\mathbf{K}\mathbf{L}^{-T} \), then back-transform eigenvectors. For non-proportional damping, we build the state matrix and extract complex poles.

File: Chapter13_Lesson4.m

% Chapter13_Lesson4.m
% Damping in MDOF Systems and Mode Shapes with Damping
%
% Demonstrates:
% 1) Undamped modes: K*phi = (w^2)*M*phi (generalized eigenproblem)
% 2) Rayleigh damping: C = alpha*M + beta*K -> diagonal modal damping matrix
% 3) Non-proportional damping: state-space eigenanalysis -> complex poles
% 4) Free response simulation with ode45
%
% This is designed as a teaching script (single file).

clear; clc;

% --------------------------
% Example system (3-DOF chain)
% --------------------------
m1 = 1.2; m2 = 1.0; m3 = 0.8;
k1 = 2500; k2 = 1800; k3 = 1200;

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

% --------------------------
% (A) Undamped modal analysis
% --------------------------
% eig(K,M) solves K*V = M*V*D
[V,D] = eig(K,M);
lam = diag(D);              % omega^2
omega = sqrt(lam);

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

disp('Undamped natural frequencies (rad/s):');
disp(omega.');

disp('Check Phi''*M*Phi (should be I):');
disp(Phi'*M*Phi);

disp('Check Phi''*K*Phi (should be diag(omega^2)):');
disp(Phi'*K*Phi);

% ------------------------------------
% (B) Rayleigh damping from two targets
% ------------------------------------
zeta1 = 0.02;   % target at mode 1
zeta3 = 0.05;   % target at mode 3

% alpha + beta*w^2 = 2*zeta*w at two modes
A = [1 omega(1)^2; 1 omega(3)^2];
b = [2*zeta1*omega(1); 2*zeta3*omega(3)];
ab = A\b;
alpha = ab(1); beta = ab(2);

C_ray = alpha*M + beta*K;
Cm_ray = Phi'*C_ray*Phi;

fprintf('Rayleigh coefficients: alpha=%.6f, beta=%.6f\n', alpha, beta);
disp('Modal damping matrix (Rayleigh) Phi''*C*Phi ~ diagonal:');
disp(Cm_ray);

zeta_modal = 0.5*(alpha./omega + beta.*omega);
disp('Modal damping ratios (Rayleigh):');
disp(zeta_modal.');

% ------------------------------------------
% (C) Non-proportional damping (example)
% ------------------------------------------
c13 = 45;  % dashpot between DOF1 and DOF3
c2g = 35;  % dashpot DOF2 to ground
C_np = zeros(3);
C_np(1,1) = C_np(1,1) + c13;
C_np(3,3) = C_np(3,3) + c13;
C_np(1,3) = C_np(1,3) - c13;
C_np(3,1) = C_np(3,1) - c13;
C_np(2,2) = C_np(2,2) + c2g;

disp('Non-proportional damping matrix C_np:');
disp(C_np);

Cm_np = Phi'*C_np*Phi;
disp('Modal damping matrix (non-proportional) Phi''*C_np*Phi (has off-diagonals):');
disp(Cm_np);

% State-space eigenanalysis
n = size(M,1);
Z = zeros(n);
I = eye(n);
Astate = [Z I; -M\K -M\C_np];
lambda = eig(Astate);

% Extract one pole per conjugate pair (imag>0)
lambda_pos = lambda(imag(lambda) > 1e-8);
[~,idx] = sort(imag(lambda_pos));
lambda_pos = lambda_pos(idx);

disp('Complex poles (sigma + j*wd), first three modes:');
for i=1:min(3, numel(lambda_pos))
    sigma = real(lambda_pos(i));
    wd = imag(lambda_pos(i));
    wn = sqrt(sigma^2 + wd^2);
    zeta = -sigma/wn;
    fprintf('Mode %d: sigma=%+.6f, wd=%.6f, wn=%.6f, zeta=%.6f\n', i, sigma, wd, wn, zeta);
end

% ------------------------------------------
% (D) Free response simulation (non-prop)
% ------------------------------------------
x0 = [0.01; 0; 0];
v0 = [0; 0; 0];
z0 = [x0; v0];

tspan = [0 5];
odefun = @(t,z) [ z(n+1:end);
                 -M\(C_np*z(n+1:end) + K*z(1:n)) ];

opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t,z] = ode45(odefun, tspan, z0, opts);

x = z(:,1:n);
figure;
plot(t,x,'LineWidth',1.2);
grid on;
xlabel('t [s]');
ylabel('displacement');
title('Non-proportional damping free response (3-DOF)');
legend('x1','x2','x3');
      

In Simulink, represent the non-proportionally damped system using the first-order state equation \( \dot{\mathbf{z} }=\mathbf{A}\mathbf{z} \) with one State-Space block, where \( \mathbf{z}=[\mathbf{x};\dot{\mathbf{x} }] \).

File: Chapter13_Lesson4.nb

(* Chapter13_Lesson4.nb
   Damping in MDOF Systems and Mode Shapes with Damping

   This notebook is provided as a plain-text Notebook expression.
   It demonstrates:
   1) Undamped generalized eigenpairs using Eigensystem[{K,M}]
   2) Rayleigh damping coefficient identification from two modal damping ratios
   3) Non-proportional damping: state-space eigenvalues (complex poles)
*)

Notebook[{
  Cell["Chapter 13 - Lesson 4: Damping in MDOF Systems and Mode Shapes with Damping", "Title"],

  Cell["(1) Define M and K for a 3-DOF chain", "Section"],
  Cell[BoxData@ToBoxes[
    (
      M = DiagonalMatrix[{1.2, 1.0, 0.8}];
      K = { {2500. + 1800., -1800., 0.},
           {-1800., 1800. + 1200., -1200.},
           {0., -1200., 1200.} };
      {M, K}
    )
  ], "Input"],

  Cell["(2) Undamped modes: K.phi = (w^2) M.phi (mass-normalized)", "Section"],
  Cell[BoxData@ToBoxes[
    Module[{vals, vecs, omega, Phi, MM, KK},
      {vals, vecs} = Eigensystem[{K, M}];  (* vals = omega^2 *)
      omega = Sqrt[vals];

      (* Mass-normalize each eigenvector: phi -> phi / Sqrt[phi^T M phi] *)
      Phi = Transpose @ Map[
        Function[phi, phi/Sqrt[phi.M.phi]],
        vecs
      ];

      MM = Transpose[Phi].M.Phi;
      KK = Transpose[Phi].K.Phi;

      <|"omega" -> omega, "Phi" -> Phi, "PhiTMPhI" -> MM, "PhiTKPhi" -> KK|>
    ]
  ], "Input"],

  Cell["(3) Rayleigh damping: C = alpha M + beta K; identify alpha,beta from two targets", "Section"],
  Cell[BoxData@ToBoxes[
    Module[{zeta1=0.02, zeta3=0.05, omega, A, b, alpha, beta, C, Cm},
      omega = Sqrt[ Eigenvalues[{K, M}] ];

      (* alpha + beta*w^2 = 2*zeta*w at modes 1 and 3 *)
      A = { {1., omega[[1]]^2}, {1., omega[[3]]^2} };
      b = {2 zeta1 omega[[1]], 2 zeta3 omega[[3]]};
      {alpha, beta} = LinearSolve[A, b];

      C = alpha M + beta K;
      Cm = Transpose[Phi].C.Phi;

      <|"alpha" -> alpha, "beta" -> beta, "Cm" -> Cm|>
    ]
  ], "Input"],

  Cell["(4) Non-proportional damping and state-space poles", "Section"],
  Cell[BoxData@ToBoxes[
    Module[{Cnp, c13=45., c2g=35., n, Z, I, Astate, lam, lamPos, info},
      Cnp = ConstantArray[0., {3, 3}];

      (* Damper between DOF1 and DOF3 *)
      Cnp[[1,1]] += c13; Cnp[[3,3]] += c13;
      Cnp[[1,3]] -= c13; Cnp[[3,1]] -= c13;

      (* DOF2 to ground *)
      Cnp[[2,2]] += c2g;

      n = Length[M];
      Z = ConstantArray[0., {n, n}];
      I = IdentityMatrix[n];

      Astate = ArrayFlatten[{
        {Z, I},
        {-Inverse[M].K, -Inverse[M].Cnp}
      }];

      lam = Eigenvalues[Astate];
      lamPos = SortBy[Select[lam, Im[#] > 10^-8 &], Im];

      info = Table[
        With[{sigma = Re[lamPos[[i]]], wd = Im[lamPos[[i]]]},
          With[{wn = Sqrt[sigma^2 + wd^2], zeta = -sigma/Sqrt[sigma^2 + wd^2]},
            <|"sigma" -> sigma, "wd" -> wd, "wn" -> wn, "zeta" -> zeta|>
          ]
        ],
        {i, 1, Min[3, Length[lamPos]]}
      ];

      <|"Cnp" -> Cnp, "poles" -> info|>
    ]
  ], "Input"]
}]
      

The notebook is provided as a plain-text Notebook expression so it can be stored in version control and reconstructed by Mathematica.

Problem 1 (Proof of decoupling under Rayleigh damping): Assume \( \mathbf{M}\succ 0 \), \( \mathbf{K} \) symmetric, and \( \mathbf{C}=\alpha\mathbf{M}+\beta\mathbf{K} \). Show that in mass-normalized modal coordinates \( \mathbf{x}=\boldsymbol{\Phi}\mathbf{q} \), the damping matrix becomes diagonal and the modal equations decouple.

Solution: Using mass-normalized modes, we have

\[ \boldsymbol{\Phi}^T\mathbf{M}\boldsymbol{\Phi}=\mathbf{I},\qquad \boldsymbol{\Phi}^T\mathbf{K}\boldsymbol{\Phi}=\boldsymbol{\Omega}^2. \]

Compute the modal damping matrix:

\[ \mathbf{C}_m=\boldsymbol{\Phi}^T\mathbf{C}\boldsymbol{\Phi} =\boldsymbol{\Phi}^T(\alpha\mathbf{M}+\beta\mathbf{K})\boldsymbol{\Phi} =\alpha\mathbf{I}+\beta\boldsymbol{\Omega}^2, \]

which is diagonal. Therefore the transformed equation \( \ddot{\mathbf{q} }+\mathbf{C}_m\dot{\mathbf{q} }+\boldsymbol{\Omega}^2\mathbf{q}=\mathbf{p}(t) \) splits into independent scalar ODEs, one per mode.

Problem 2 (Identify Rayleigh coefficients from two modal damping ratios): Suppose the damping ratios at two modes are known: \( \zeta_i \) at \( \omega_i \) and \( \zeta_j \) at \( \omega_j \). Derive closed-form expressions for \( \alpha \) and \( \beta \).

Solution: From Rayleigh damping,

\[ \zeta_r=\frac{\alpha}{2\omega_r}+\frac{\beta\omega_r}{2} \quad\Longleftrightarrow\quad 2\zeta_r\omega_r=\alpha+\beta\omega_r^2. \]

Apply to modes \( i \) and \( j \):

\[ \alpha+\beta\omega_i^2=2\zeta_i\omega_i,\qquad \alpha+\beta\omega_j^2=2\zeta_j\omega_j. \]

Subtract the equations to solve for \( \beta \), then back-substitute for \( \alpha \):

\[ \beta = \frac{2(\zeta_j\omega_j-\zeta_i\omega_i)}{\omega_j^2-\omega_i^2},\qquad \alpha = 2\zeta_i\omega_i - \beta\omega_i^2. \]

Problem 3 (Damped modal frequency and logarithmic decrement): For one decoupled modal equation \( \ddot{q}+2\zeta\omega \dot{q}+\omega^2 q=0 \) with \( 0<\zeta<1 \), derive \( \omega_d \) and the logarithmic decrement \( \delta \).

Solution: The characteristic equation is

\[ \lambda^2+2\zeta\omega\lambda+\omega^2=0 \quad\Longrightarrow\quad \lambda=-\zeta\omega\pm j\,\omega\sqrt{1-\zeta^2}. \]

Hence \( \omega_d=\omega\sqrt{1-\zeta^2} \). For successive peaks separated by one period \( T_d=2\pi/\omega_d \), the amplitude ratio is \( e^{\zeta\omega T_d} \), giving:

\[ \delta \equiv \ln\!\left(\frac{x_k}{x_{k+1} }\right) = \zeta\omega T_d = \frac{2\pi\zeta}{\sqrt{1-\zeta^2} }. \]

Problem 4 (Equivalence of QEP and state-space eigenproblem): Starting from \( \mathbf{M}\ddot{\mathbf{x} }+\mathbf{C}\dot{\mathbf{x} }+\mathbf{K}\mathbf{x}=\mathbf{0} \), show that the exponential trial \( \mathbf{x}=\boldsymbol{\phi}e^{\lambda t} \) yields the QEP, and then show that the same eigenvalues arise from the state matrix \( \mathbf{A}=\begin{bmatrix}\mathbf{0}&\mathbf{I}\\-\mathbf{M}^{-1}\mathbf{K}&-\mathbf{M}^{-1}\mathbf{C}\end{bmatrix} \).

Solution: Substituting \( \mathbf{x}=\boldsymbol{\phi}e^{\lambda t} \) gives \( \dot{\mathbf{x} }=\lambda\boldsymbol{\phi}e^{\lambda t} \) and \( \ddot{\mathbf{x} }=\lambda^2\boldsymbol{\phi}e^{\lambda t} \), hence:

\[ (\lambda^2\mathbf{M}+\lambda\mathbf{C}+\mathbf{K})\boldsymbol{\phi}=\mathbf{0}. \]

Now define \( \mathbf{z}=\begin{bmatrix}\mathbf{x}\\\dot{\mathbf{x} }\end{bmatrix} \). For a modal solution, let \( \mathbf{x}=\boldsymbol{\phi}e^{\lambda t} \), so \( \dot{\mathbf{x} }=\lambda\boldsymbol{\phi}e^{\lambda t} \), and therefore \( \mathbf{z}=\begin{bmatrix}\boldsymbol{\phi}\\\lambda\boldsymbol{\phi}\end{bmatrix}e^{\lambda t} \). Substituting into \( \dot{\mathbf{z} }=\mathbf{A}\mathbf{z} \) yields:

\[ \lambda\begin{bmatrix}\boldsymbol{\phi}\\\lambda\boldsymbol{\phi}\end{bmatrix} = \begin{bmatrix} \mathbf{0} & \mathbf{I}\\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{C} \end{bmatrix} \begin{bmatrix}\boldsymbol{\phi}\\\lambda\boldsymbol{\phi}\end{bmatrix} = \begin{bmatrix} \lambda\boldsymbol{\phi}\\ -\mathbf{M}^{-1}(\mathbf{K}+\lambda\mathbf{C})\boldsymbol{\phi} \end{bmatrix}. \]

The first block row is an identity. The second block row implies \( \lambda^2\mathbf{M}\boldsymbol{\phi}+\lambda\mathbf{C}\boldsymbol{\phi}+\mathbf{K}\boldsymbol{\phi}=\mathbf{0} \), i.e., exactly the QEP. Hence the eigenvalues coincide.

Problem 5 (Extract damping ratio from a complex pole): A non-proportionally damped mode yields the pole \( \lambda = -0.4 + j\,12 \). Compute \( \omega_n \), \( \omega_d \), and \( \zeta \).

Solution: Here \( \sigma=-0.4 \) and \( \omega_d=12 \). Then

\[ \omega_n=\sqrt{\sigma^2+\omega_d^2}=\sqrt{0.4^2+12^2}=\sqrt{144.16}\approx 12.0067, \qquad \zeta=-\frac{\sigma}{\omega_n}=\frac{0.4}{12.0067}\approx 0.0333. \]

The damped frequency is \( \omega_d=12 \) rad/s (already given as the imaginary part).

We formulated damped MDOF dynamics in matrix form and showed that modal decoupling hinges on whether \( \boldsymbol{\Phi}^T\mathbf{C}\boldsymbol{\Phi} \) is diagonal. Rayleigh damping \( \mathbf{C}=\alpha\mathbf{M}+\beta\mathbf{K} \) guarantees classical damping, preserves real mode shapes, and yields independent modal damping ratios \( \zeta_r=\alpha/(2\omega_r)+(\beta\omega_r)/2 \). For non-proportional damping, we introduced the quadratic eigenvalue problem and the state-space matrix formulation that naturally produces complex poles and complex mode shapes.