Generalized Eigenvalue Example — Mechanical Vibration

This example solves the generalized eigenvalue problem

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

for a 3-DOF spring-mass chain where K is the stiffness matrix and M is the diagonal mass matrix. The natural frequencies are \(\omega_i = \sqrt{\lambda_i}\) (rad/s) and the eigenvectors are the vibration mode shapes.

eigensystems.gvcsf() is used because the pair (K, M) is symmetric-definite, which guarantees real, positive eigenvalues and M-orthonormal mode shapes.

Example Code

"""IMSL generalized eigenvalue example: mechanical vibration natural frequencies.

Demonstrates:
- Formulate the generalized eigenvalue problem Ax = λBx for a 3-DOF spring-mass system.
  A = stiffness matrix [[4,-1,0],[-1,4,-1],[0,-1,4]]
  B = mass matrix      [[2,0,0],[0,2,0],[0,0,2]]
- Solve using :func:`eigensystems.gvcsf` (symmetric-definite pair).
- Compute natural frequencies as ω = sqrt(λ).
- Plot mode shapes (eigenvectors) for each natural frequency.
- Save to test_output/example_imsl_generalized_eigen.svg

Outputs:
- Printed table of eigenvalues and natural frequencies.
- SVG figure saved to test_output/.
- Console output captured to docs/_generated_text/generalized_eigen_output.txt
"""
from __future__ import annotations

from pathlib import Path

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np

from eigensystems import gvcsf


def run_demo_generalized_eigen() -> dict:
    """Run IMSL GVCSF example: natural frequencies of a 3-DOF spring-mass system.

    Solves the generalized eigenvalue problem K*phi = lambda*M*phi where K is
    the stiffness matrix and M is the mass matrix.  Natural frequencies are
    omega = sqrt(lambda) and the eigenvectors are the vibration mode shapes.

    Args:
        None

    Returns:
        dict: Result dictionary with keys ``eigenvalues`` (np.ndarray),
            ``natural_frequencies`` (np.ndarray), ``eigenvectors`` (np.ndarray),
            and ``plot_path`` (str).
    """
    # Stiffness matrix (tri-diagonal spring chain)
    K = np.array([
        [ 4.0, -1.0,  0.0],
        [-1.0,  4.0, -1.0],
        [ 0.0, -1.0,  4.0],
    ])

    # Mass matrix (diagonal, equal lumped masses)
    M = np.array([
        [2.0, 0.0, 0.0],
        [0.0, 2.0, 0.0],
        [0.0, 0.0, 2.0],
    ])

    result = gvcsf(K, M)
    eigenvalues = result.eigenvalues.real
    eigenvectors = result.eigenvectors.real
    natural_frequencies = np.sqrt(np.abs(eigenvalues))

    print("\nIMSL GVCSF Example: Generalized Eigenvalue Problem — 3-DOF Spring-Mass System")
    print("  Solve K*phi = lambda*M*phi   (omega = sqrt(lambda))")
    print("-" * 62)
    print(f"{'Mode':<8} {'Eigenvalue λ':>16} {'Nat. freq ω (rad/s)':>22}")
    print("-" * 62)
    for i, (lam, omega) in enumerate(zip(eigenvalues, natural_frequencies)):
        print(f"  Mode {i+1:<3}  {lam:>16.6f}  {omega:>22.6f}")
    print("-" * 62)

    print("\n  Mode shape vectors (columns = eigenvectors):")
    headers = [f"Mode {i+1}" for i in range(len(eigenvalues))]
    header_line = f"  {'DOF':<6}  " + "  ".join(f"{h:>12}" for h in headers)
    print(header_line)
    print("  " + "-" * (len(header_line) - 2))
    for dof in range(eigenvectors.shape[0]):
        vals = "  ".join(f"{eigenvectors[dof, m]:>12.6f}" for m in range(eigenvectors.shape[1]))
        print(f"  DOF {dof+1:<3}  {vals}")

    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_generalized_eigen.svg"

    n_modes = eigenvectors.shape[1]
    dof_labels = [f"DOF {i+1}" for i in range(eigenvectors.shape[0])]
    colors = ["#be123c", "#0369a1", "#059669"]

    fig, axes = plt.subplots(1, n_modes, figsize=(12, 5), sharey=False)
    fig.suptitle(
        "IMSL GVCSF: Mode Shapes of a 3-DOF Spring-Mass System\n"
        "K·φ = λ·M·φ  →  ω = √λ",
        fontsize=11,
    )

    for m, ax in enumerate(axes):
        mode_vec = eigenvectors[:, m]
        x_pos = np.arange(len(mode_vec))
        ax.bar(x_pos, mode_vec, color=colors[m % len(colors)], alpha=0.85,
               edgecolor="black", linewidth=0.8)
        ax.axhline(0, color="gray", linewidth=0.8, linestyle="--")
        ax.set_xticks(x_pos)
        ax.set_xticklabels(dof_labels, fontsize=9)
        ax.set_ylabel("Amplitude")
        ax.set_title(
            f"Mode {m+1}\nλ = {eigenvalues[m]:.4f}\nω = {natural_frequencies[m]:.4f} rad/s",
            fontsize=9,
        )

    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)
    print(f"\nPlot saved: {plot_path}")

    return {
        "eigenvalues": eigenvalues,
        "natural_frequencies": natural_frequencies,
        "eigenvectors": eigenvectors,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_generalized_eigen()

Plot Output

The three vibration modes of the spring-mass chain. Mode 1 is an in-phase motion (lowest frequency); Mode 3 is the highest-frequency alternating pattern.

Console Output

IMSL GVCSF Example: Generalized Eigenvalue Problem ΓÇö 3-DOF Spring-Mass System
  Solve K*phi = lambda*M*phi   (omega = sqrt(lambda))
--------------------------------------------------------------
Mode         Eigenvalue λ    Nat. freq ω (rad/s)
--------------------------------------------------------------
  Mode 1            1.292893                1.137055
  Mode 2            2.000000                1.414214
  Mode 3            2.707107                1.645329
--------------------------------------------------------------

  Mode shape vectors (columns = eigenvectors):
  DOF           Mode 1        Mode 2        Mode 3
  ------------------------------------------------
  DOF 1       -0.353553     -0.500000      0.353553
  DOF 2       -0.500000      0.000000     -0.500000
  DOF 3       -0.353553      0.500000      0.353553

Plot saved: test_output\example_imsl_generalized_eigen.svg