All Generalized Eigenvalue Forms — GVLRG, GVCRG, GVLSF, EVCSF

This example solves the structural dynamics eigenvalue problem \(K \mathbf{v} = \lambda M \mathbf{v}\) for a 4-DOF spring-mass chain using four IMSL-compatible routines:

The stiffness matrix is the tridiagonal \(K = \text{tridiag}(-1,\,2,\,-1)\) and the mass matrix is the identity.

Example Code

"""IMSL generalized eigenvalue examples: gvlrg, gvcrg, gvlsf, evcsf.

Structural dynamics problem: K*v = λ*M*v
- K = 4×4 tridiagonal stiffness matrix (diag=2, off-diag=-1)
- M = 4×4 identity mass matrix

Demonstrates:
- gvlrg  — generalized eigenvalues only (general solver)
- gvcrg  — generalized eigenvalues + eigenvectors
- gvlsf  — generalized eigenvalues for symmetric-definite pair (ascending)
- evcsf  — eigenvalues+vectors of K alone; orthonormality check

Outputs:
- Printed tables of eigenvalues and orthonormality error
- SVG plot: mode shapes and frequency spectrum
- Save SVG to test_output/example_imsl_all_generalized.svg
- Console output captured to docs/_generated_text/all_generalized_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 gvlrg, gvcrg, gvlsf, evcsf


def _tridiag(n: int, diag: float, offdiag: float) -> np.ndarray:
    """Build an n×n symmetric tridiagonal matrix.

    Args:
        n (int): Size of the matrix.
        diag (float): Value on the main diagonal.
        offdiag (float): Value on the super- and sub-diagonals.

    Returns:
        np.ndarray: n×n symmetric tridiagonal matrix.
    """
    a = np.zeros((n, n))
    np.fill_diagonal(a, diag)
    np.fill_diagonal(a[1:], offdiag)
    np.fill_diagonal(a[:, 1:], offdiag)
    return a


def run_demo_all_generalized() -> dict:
    """Run all generalized eigenvalue demonstrations.

    Solves the structural vibration problem K·v = λ·M·v for a
    4-DOF spring-mass chain.

    Args:
        None

    Returns:
        dict: Keys ``gvlrg_eigs``, ``gvcrg_eigs``, ``gvlsf_eigs``,
            ``evcsf_eigs``, ``mode_shapes`` (np.ndarray), ``plot_path`` (str).
    """
    n = 4
    K = _tridiag(n, 2.0, -1.0)  # stiffness
    M = np.eye(n)                # mass (identity)

    # ------------------------------------------------------------------
    # gvlrg — general generalized eigenvalues (no vectors)
    # ------------------------------------------------------------------
    res_gvlrg = gvlrg(K, M)
    eigs_gvlrg = np.sort(res_gvlrg.eigenvalues.real)

    print("\nIMSL GVLRG Example: Generalized Eigenvalues K·v = λ·M·v")
    print("K = 4×4 tridiagonal(2, -1),  M = I")
    print("-" * 50)
    print(f"{'Index':<8} {'λ (real part)':>16} {'λ (imag part)':>16}")
    print("-" * 50)
    for i, ev in enumerate(sorted(res_gvlrg.eigenvalues, key=lambda z: z.real)):
        print(f"  {i+1:<6} {ev.real:>16.8f} {ev.imag:>+16.2e}")

    # ------------------------------------------------------------------
    # gvcrg — general generalized eigenvalues + eigenvectors
    # ------------------------------------------------------------------
    res_gvcrg = gvcrg(K, M)
    # Sort by real part of eigenvalue
    idx_sort = np.argsort(res_gvcrg.eigenvalues.real)
    eigs_gvcrg = res_gvcrg.eigenvalues[idx_sort]
    mode_shapes = res_gvcrg.eigenvectors[:, idx_sort].real

    print("\nIMSL GVCRG Example: Generalized Eigenvalues+Vectors")
    print("-" * 50)
    for i, ev in enumerate(eigs_gvcrg):
        print(f"  Mode {i+1}: λ={ev.real:.6f}, shape={mode_shapes[:, i]}")

    # ------------------------------------------------------------------
    # gvlsf — symmetric-definite generalized eigenvalues (ascending)
    # ------------------------------------------------------------------
    res_gvlsf = gvlsf(K, M)
    freqs_gvlsf = np.sqrt(np.maximum(res_gvlsf.eigenvalues.real, 0))

    print("\nIMSL GVLSF Example: Symmetric-Definite Generalized Eigenvalues")
    print("-" * 50)
    print(f"Eigenvalues (ascending): {res_gvlsf.eigenvalues}")
    print(f"Natural frequencies ω=√λ: {freqs_gvlsf}")

    # Verify gvlsf agrees with gvlrg
    diff = np.max(np.abs(np.sort(res_gvlsf.eigenvalues)
                         - np.sort(eigs_gvlrg)))
    print(f"Max difference |gvlsf - gvlrg| = {diff:.2e}")

    # ------------------------------------------------------------------
    # evcsf — symmetric eigenvalues+vectors of K alone
    # ------------------------------------------------------------------
    res_evcsf = evcsf(K)
    V = res_evcsf.eigenvectors

    print("\nIMSL EVCSF Example: Symmetric Eigenvalues+Vectors of K")
    print("-" * 50)
    print(f"Eigenvalues of K (ascending): {res_evcsf.eigenvalues}")
    orth_err = np.max(np.abs(V.T @ V - np.eye(n)))
    print(f"Orthonormality error ‖V^T V - I‖_max = {orth_err:.2e}")

    # ------------------------------------------------------------------
    # Plot: mode shapes (left) and frequency spectrum (right)
    # ------------------------------------------------------------------
    output_dir = Path("test_output")
    output_dir.mkdir(parents=True, exist_ok=True)
    plot_path = output_dir / "example_imsl_all_generalized.svg"

    node_pos = np.arange(1, n + 1)
    colors = ["#be123c", "#0369a1", "#059669", "#d97706"]

    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 5))

    # Mode shapes
    for k in range(n):
        shape = mode_shapes[:, k]
        shape = shape / np.max(np.abs(shape))  # normalise to ±1
        ax1.plot(node_pos, shape, marker="o", color=colors[k],
                 label=f"Mode {k+1}  λ={eigs_gvcrg[k].real:.3f}")
    ax1.axhline(0, color="gray", linewidth=0.6, linestyle="--")
    ax1.set_xlabel("Node index")
    ax1.set_ylabel("Normalised displacement")
    ax1.set_title("GVCRG: Mode Shapes  K·v = λ·M·v")
    ax1.legend(fontsize=8)
    ax1.set_xticks(node_pos)

    # Frequency spectrum from gvlsf
    ax2.bar(np.arange(1, n + 1), freqs_gvlsf, color="#0369a1",
            alpha=0.8, edgecolor="#1e3a5f")
    for i, f in enumerate(freqs_gvlsf):
        ax2.text(i + 1, f + 0.02 * max(freqs_gvlsf), f"{f:.3f}",
                 ha="center", va="bottom", fontsize=9)
    ax2.set_xticks(np.arange(1, n + 1))
    ax2.set_xticklabels([f"Mode {i}" for i in range(1, n + 1)])
    ax2.set_ylabel("Natural frequency  ω = √λ")
    ax2.set_title("GVLSF: Natural Frequency Spectrum")

    fig.suptitle("Generalized Eigenvalue Problems — IMSL Routines", fontsize=13,
                 fontweight="bold")
    fig.tight_layout()
    fig.savefig(plot_path, format="svg")
    plt.close(fig)
    print(f"\nPlot saved: {plot_path}")

    return {
        "gvlrg_eigs": res_gvlrg.eigenvalues,
        "gvcrg_eigs": res_gvcrg.eigenvalues,
        "gvlsf_eigs": res_gvlsf.eigenvalues,
        "evcsf_eigs": res_evcsf.eigenvalues,
        "mode_shapes": mode_shapes,
        "plot_path": str(plot_path),
    }


if __name__ == "__main__":
    run_demo_all_generalized()

Plot Output

Mode shapes from GVCRG and natural frequency spectrum from GVLSF

Left: normalised mode shapes for each of the four vibration modes (GVCRG). Right: natural frequencies \(\omega = \sqrt{\lambda}\) (GVLSF).

Console Output

IMSL GVLRG Example: Generalized Eigenvalues K┬╖v = ╬╗┬╖M┬╖v
K = 4×4 tridiagonal(2, -1),  M = I
--------------------------------------------------
Index       ╬╗ (real part)    ╬╗ (imag part)
--------------------------------------------------
  1            0.38196601        +0.00e+00
  2            1.38196601        +0.00e+00
  3            2.61803399        +0.00e+00
  4            3.61803399        +0.00e+00

IMSL GVCRG Example: Generalized Eigenvalues+Vectors
--------------------------------------------------
  Mode 1: ╬╗=0.381966, shape=[0.37174803 0.60150096 0.60150096 0.37174803]
  Mode 2: ╬╗=1.381966, shape=[ 0.60150096  0.37174803 -0.37174803 -0.60150096]
  Mode 3: ╬╗=2.618034, shape=[ 0.60150096 -0.37174803 -0.37174803  0.60150096]
  Mode 4: ╬╗=3.618034, shape=[-0.37174803  0.60150096 -0.60150096  0.37174803]

IMSL GVLSF Example: Symmetric-Definite Generalized Eigenvalues
--------------------------------------------------
Eigenvalues (ascending): [0.38196601 1.38196601 2.61803399 3.61803399]
Natural frequencies ω=√λ: [0.61803399 1.1755705  1.61803399 1.90211303]
Max difference |gvlsf - gvlrg| = 3.11e-15

IMSL EVCSF Example: Symmetric Eigenvalues+Vectors of K
--------------------------------------------------
Eigenvalues of K (ascending): [0.38196601 1.38196601 2.61803399 3.61803399]
Orthonormality error ΓÇûV^T V - IΓÇû_max = 6.66e-16

Plot saved: test_output\example_imsl_all_generalized.svg