Modal eigenvalue problem and shift-invert Lanczos#

Free-vibration modal analysis solves the generalised symmetric-definite eigenvalue problem

(1)#\[\mathbf{K}\, \boldsymbol{\phi} \;=\; \omega^{2}\, \mathbf{M}\, \boldsymbol{\phi},\]

with \(\mathbf{K}\) the assembled stiffness matrix (Global assembly) and \(\mathbf{M}\) the mass matrix. Both are sparse, symmetric, and positive (semi-)definite — \(\mathbf{K}\) becomes positive-definite once the rigid- body modes are removed by Dirichlet BCs (Boundary-condition elimination).

This chapter develops the algorithm scipy.sparse.linalg.eigsh() runs underneath femorph-solver’s Model.modal_solve(): an Arnoldi / Lanczos iteration on the shift-inverted operator that converts the smallest-eigenvalue search into the largest- magnitude search ARPACK was designed for.

Why not just solve \(\mathbf{K}^{-1}\, \mathbf{M}\)?#

The straightforward power iteration on \(\mathbf{K}^{-1}\, \mathbf{M}\) extracts the largest eigenvalue of that operator — equivalently the smallest \(\omega^{2}\) of (1). But power iteration converges at a rate set by the spectral gap \(|\lambda_2 / \lambda_1|\), which for the lowest few modes is typically near 1. Lanczos accelerates this by projecting onto a Krylov subspace and computing all extreme eigenvalues simultaneously; ARPACK is the canonical implementation (Lehoucq, Sorensen & Yang 1998).

Shift-invert Lanczos#

Define the shifted-inverted operator

(2)#\[\mathbf{S}_{\sigma} \;=\; (\mathbf{K} - \sigma\, \mathbf{M})^{-1}\, \mathbf{M}.\]

The eigenvalues of \(\mathbf{S}_{\sigma}\) are \(\mu = 1 / (\omega^{2} - \sigma)\). Eigenvalues \(\omega^{2}\) near \(\sigma\) map to the largest-magnitude eigenvalues of \(\mathbf{S}_{\sigma}\) — exactly what Lanczos converges to fastest. The default shift in femorph-solver is \(\sigma = 0\), so \(\mathbf{S}_{0} = \mathbf{K}^{-1}\, \mathbf{M}\) and the lowest \(n_\mathrm{modes}\) are returned by ascending \(\omega^{2}\).

Implementation:

Factor \(\mathbf{K} - \sigma\, \mathbf{M}\) once through the registered linear backend (Pardiso / CHOLMOD / MUMPS / SuperLU auto-chain — see femorph_solver.solvers.linear.list_linear_solvers()).
Each Lanczos iteration applies \(\mathbf{S}_{\sigma}\, \mathbf{v}\) by solving \((\mathbf{K} - \sigma \mathbf{M})\, \mathbf{w} = \mathbf{M}\, \mathbf{v}\) against the cached factor.
ARPACK’s reverse-communication interface returns \((\mu, \mathbf{q})\) pairs; femorph-solver back- transforms via \(\omega^{2} = \sigma + 1/\mu\) and returns the original eigenvectors \(\boldsymbol{\phi} = \mathbf{q}\) (Lanczos eigenvectors of \(\mathbf{S}_{\sigma}\) and of the original GEVP (1) are the same up to back-transformation of the eigenvalue).

The whole flow factors \(\mathbf{K} - \sigma \mathbf{M}\) once and re-uses the factor for every Lanczos matrix-vector product — a critical optimisation because the factor cost dominates a sparse direct solve. See Ericsson & Ruhe (1980) for the original treatment.

Cluster handling near \(\sigma\)#

Eigenvalues clustered near the shift produce ill-conditioned \(\mathbf{S}_{\sigma}\) factors and slow Lanczos convergence. The standard remedy is shift-and-invert with a slightly off-cluster shift: pick \(\sigma\) near but not at the cluster. ARPACK exposes this via the sigma parameter to eigsh; femorph-solver currently uses \(\sigma = 0\) (suitable for static-stress-free models where rigid-body modes have already been removed by BCs) but the API allows an explicit override. When a frequency-shifted sweep is needed (e.g., to extract modes near 1 kHz on a model whose first 50 modes are below 100 Hz), set Model.modal_solve(n_modes=N, sigma=2 * pi * 1000).

Real-symmetric vs complex-Hermitian#

For an isotropic linear-elastic structure with conforming mesh, both \(\mathbf{K}\) and \(\mathbf{M}\) are real symmetric. scipy.sparse.linalg.eigsh (the real-symmetric Lanczos) is the natural choice and femorph-solver’s modal solver dispatches there.

The cyclic-symmetry path (Global assembly) reduces the global GEVP to a per-harmonic complex-Hermitian sub-problem \(\mathbf{K}_{k}\, \boldsymbol{\phi}_{k} = \omega^{2}\, \mathbf{M}_{k}\, \boldsymbol{\phi}_{k}\). Rather than calling the complex eigs (which is slower and lacks positive-definite support), femorph-solver applies the Grimes–Lewis–Simon (1994) real 2n-augmentation:

\[\begin{split}\tilde{\mathbf{K}}_{k} = \begin{bmatrix} \mathrm{Re}(\mathbf{K}_{k}) & -\mathrm{Im}(\mathbf{K}_{k}) \\ \mathrm{Im}(\mathbf{K}_{k}) & \mathrm{Re}(\mathbf{K}_{k}) \end{bmatrix}, \qquad \tilde{\mathbf{M}}_{k} = \begin{bmatrix} \mathrm{Re}(\mathbf{M}_{k}) & -\mathrm{Im}(\mathbf{M}_{k}) \\ \mathrm{Im}(\mathbf{M}_{k}) & \mathrm{Re}(\mathbf{M}_{k}) \end{bmatrix},\end{split}\]

which are real-symmetric and share the same eigenvalue spectrum (each appearing with multiplicity 2). The same eigsh shift-invert path runs unchanged on the doubled matrices; the eigenvectors recombine into a complex pair on the way out.

Mode normalisation#

Two conventions for normalising \(\boldsymbol{\phi}\) see common use:

Mass-orthonormal — \(\boldsymbol{\phi}_{i}^{\!\top}\, \mathbf{M}\, \boldsymbol{\phi}_{j} = \delta_{ij}\). Used by Bathe §8 and by femorph-solver’s modal API.
Unit-amplitude — \(\max_n |\phi_{i,n}| = 1\). Used by some legacy MAPDL pipelines.

mass-orthonormal is the natural choice for modal-superposition post-processing because the modal mass matrix becomes the identity and the harmonic / transient response decouples per mode.

References#

Parlett, B. N. (1998) The Symmetric Eigenvalue Problem, SIAM, §13 (Lanczos), §14 (shift-invert).
Lehoucq, R. B., Sorensen, D. C. and Yang, C. (1998) ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, SIAM Software Environments and Tools 6.
Ericsson, T. and Ruhe, A. (1980) “The spectral transformation Lanczos method for the numerical solution of large sparse generalized symmetric eigenvalue problems,” Mathematics of Computation 35, 1251–1268 (the original shift-invert formulation).
Saad, Y. (2011) Numerical Methods for Large Eigenvalue Problems, 2nd ed., SIAM, §6 (Lanczos for symmetric problems), §13 (shift-and-invert strategies).
Bathe, K.-J. (2014) Finite Element Procedures, 2nd ed., §11 (eigensolution methods), §12 (mode-superposition dynamics).
Wilkinson, J. H. (1965) The Algebraic Eigenvalue Problem, Oxford UP — foundational treatment of symmetric-definite eigenproblem stability.
Grimes, R. G., Lewis, J. G. and Simon, H. D. (1994) “A shifted block Lanczos algorithm for solving sparse symmetric generalized eigenproblems,” SIAM J. Matrix Analysis and Applications 15, 228–272 (real 2n augmentation for the complex-Hermitian cyclic-symmetry sub-problem).