Frequency analysis (structsolve.freq)

structsolve.freq(K, M, tol=0, sparse_solver=True, silent=False, sort=True, num_eigvalues=25, num_eigvalues_print=5, skip_null_cols=False, symmetric=False, check_rtol=1e-06, max_dense_size=2000, linear_solver='splu')

Frequency analysis

Calculates the eigenvalues lambda2 and eigenvectors u of the free-vibration problem of K and M, i.e. the eigenpairs with (K + lambda2*M) @ u = 0, where lambda2 = -omega_n**2 and omega_n is the natural frequency in rad/s. The null rows and columns of K are removed from both matrices before solving the eigenvalue problem, unless skip_null_cols=True.

Solver policy. As in lb(), with sparse_solver=True the sparse solver runs first for every problem size, and its eigenpairs are verified:

• the relative residual of each eigenpair must not be larger than check_rtol, unless its backward error is at the round-off, see check_rtol;

• with symmetric=True, no eigenvalue omega_n**2 lower than the highest one returned may be missing, see _check_freq().

When the verification fails or the sparse solver raises an error, e.g. ARPACK did not converge, the dense solver is used as a last resort if the problem has at most max_dense_size dofs, with a linear_buckling.DenseFallbackWarning giving the reasons. It returns the num_eigvalues lowest eigenvalues with symmetric=True, otherwise the num_eigvalues eigenvalues closest to the shift, as the sparse solver, and its eigenpairs are verified as well. Otherwise a RuntimeError is raised. Until version 0.5.3, the sparse solver only issued a warning when the residual check failed.

With Intel MKL 2024.2.0 to 2025.0.0, ncv of scipy.sparse.linalg.eigsh() is capped, see lb() and structsolve.arpackutils.

Parameters:

Ksparse_matrix

Stiffness matrix. Should include initial stress stiffness matrix, aerodynamic matrix and so forth when applicable.

Msparse_matrix

Mass matrix.

tolfloat, optional

A tolerance value passed to scipy.sparse.linalg.eigs() and scipy.sparse.linalg.eigsh().

sparse_solverbool, optional

With True, the default, the sparse solver always runs first, whatever the size of the problem: scipy.sparse.linalg.eigs(), or scipy.sparse.linalg.eigsh() with symmetric=True, in shift-invert mode with a negative shift estimated from the matrices, calculating the num_eigvalues eigenvalues closest to the shift. Its eigenpairs are verified, and a dense solver is only used as a last resort, see above. With False, the full problem is explicitly solved with a dense solver, scipy.linalg.eig(), or scipy.linalg.eigh() with symmetric=True, without trying the sparse solver, which calculates all eigenvalues, e.g. for the kinetic criterion of stability, and is a reference for validations.

silentbool, optional

A boolean to tell whether the log messages should be printed.

sortbool, optional

Sort the eigenvalues and the eigenvectors together by increasing natural frequency, keeping only those with a natural frequency larger than 1e-6 rad/s, e.g. without the rigid-body modes. The sort key is the real part of \(\omega_n = \sqrt{-\lambda^2}\) rounded to one decimal, then its imaginary part rounded to one decimal, then, between equal rounded keys, the exact real and imaginary parts of \(\omega_n\), then those of \(\omega_n^2 = -\lambda^2\), which separate two \(\omega_n^2\) differing at the round-off that have the same \(\omega_n\) after the square root. Only identical eigenvalues keep the order of the solver, such that a multiple eigenvalue returned with round-off differences, in an order that may depend on the LAPACK build, is sorted by increasing exact value. With False the eigenpairs are returned in the order of the solver. Until version 0.8.0, only the eigenvectors were sorted and filtered, and the eigenpairs with equal rounded keys kept the order of the solver.

num_eigvaluesint, optional

Number of calculated eigenvalues with the sparse solver, limited to the size of M minus 2.

num_eigvalues_printint, optional

Number of eigenvalues to print.

skip_null_colsbool, optional

If True, skip the removal of null columns from the matrices. Use only when K and M are known to be non-singular.

symmetricbool or None, optional

False, the default, uses the general solvers described above, valid for unsymmetric matrices, e.g. K with an aerodynamic matrix or with the load stiffness of follower loads. True uses the solvers of symmetric matrices, scipy.sparse.linalg.eigsh() in shift-invert mode or scipy.linalg.eigh(), which require symmetric K and M and a positive definite M, and return real eigenvalues and eigenvectors normalized with respect to M. None selects them when linear_buckling.is_symmetric() is true for both matrices.

check_rtolfloat or None, optional

Maximum relative residual of an eigenpair, \(||K u + \lambda^2 M u||/(||K u|| + |\lambda^2| ||M u||)\). For a rigid-body mode, whose \(||K u|| + |\lambda^2| ||M u||\) is not larger than linear_buckling.RIGID_BODY_RTOL times the matrix-norm scale \((||K||_F + |\lambda^2| ||M||_F) ||u||\), the residual is relative to this scale, see linear_buckling._relative_residuals(). An eigenpair with a larger relative residual is also accepted when its normwise backward error \(||K u + \lambda^2 M u||/((||K||_F + |\lambda^2| ||M||_F) ||u||)\) is not larger than linear_buckling.BACKWARD_ERROR_RTOL, the machine precision, e.g. an eigenvalue close to zero of a model pre-stressed close to buckling, see linear_buckling._backward_errors(). Any other larger residual of the sparse solver is a failed verification, see above. For sparse_solver=False, a warning is issued. None skips the check.

max_dense_sizeint, optional

Maximum size of the dense problem solved as a last resort when the sparse solver failed, see above. Larger problems raise a RuntimeError instead, and 0 disables the dense fallback.

linear_solverstr, optional

'splu' (default) for the sparse factorizations of SuperLU, or 'pardiso' for those of PARDISO, which requires the shared library mkl_rt of Intel oneMKL, see structsolve.pardiso, with symmetric=True (or None detecting symmetric matrices) only; see the Notes.

Returns:

lambda2ndarray

Complex array with the eigenvalues \(\lambda^2 = -\omega_n^2\), real with symmetric=True, sorted and filtered as described in sort, otherwise in the order returned by the eigenvalue solver.

eigvecsndarray

The \(i^{th}\) eigenvector is eigvecs[:, i], the eigenvector of lambda2[i], with the size of the original matrices. The sparse solver returns eigenvectors normalized with respect to M, the dense solver returns eigenvectors with unit Euclidean norm.

Raises:

RuntimeError

When the sparse solver failed and the dense fallback was not possible or failed too, listing the attempted solvers and their failures.

Warns:

.linear_buckling.DenseFallbackWarning

When the sparse solver failed and the dense solver was used.

Notes

The free-vibration eigenvalue problem is

\[(\boldsymbol{K} + \lambda^2 \boldsymbol{M}) \boldsymbol{u} = \boldsymbol{0}\]

where \(\lambda^2 = -\omega_n^2\) and \(\omega_n\) is the natural frequency in rad/s. The sparse solvers solve \(\boldsymbol{K} \boldsymbol{u} = \omega_n^2 \boldsymbol{M} \boldsymbol{u}\) in shift-invert mode, which returns the \(\omega_n^2\) nearest the shift \(\sigma\) ([1], Chapter 3, Section 3.2, Eq. (3.2.2)). The negative shift is \(\sigma = -|\boldsymbol{y}^\top \boldsymbol{K} \boldsymbol{y} / \boldsymbol{y}^\top \boldsymbol{M} \boldsymbol{y}|\), the Rayleigh quotient of \(\boldsymbol{y} = \boldsymbol{K}^{-1} \boldsymbol{M} \boldsymbol{x}\) for a fixed random vector \(\boldsymbol{x}\), i.e. one step of inverse iteration towards the lowest \(\omega_n^2\), or \(\sigma = -1\) when K is singular or the solution is inaccurate. With all \(\omega_n^2 \ge 0\), the eigenvalues nearest a negative shift are the lowest ones.

Verification. With symmetric=True, the number of negative pivots of \(\boldsymbol{K} - s \boldsymbol{M}\), for \(s\) slightly lower than the highest \(\omega_n^2\) returned, is the number of eigenvalues below \(s\) (Sylvester’s law of inertia), see _check_freq().

PARDISO. With linear_solver='pardiso' and symmetric matrices, K of the estimate of the shift, \(\boldsymbol{K} - \sigma \boldsymbol{M}\) of the shift-invert mode, passed to scipy.sparse.linalg.eigsh() as OPinv (also in the passes of arpackutils.capped_eigsh()), and \(\boldsymbol{K} - s \boldsymbol{M}\) of the inertia count are factorized by PARDISO, the symmetric indefinite factorization of static._factorize(), see lb(); ARPACK remains the eigenvalue solver. The count of negative eigenvalues of PARDISO is exact also when SuperLU would need a row interchange. As in lb(), SuperLU is used instead for a matrix whose factorization by PARDISO is inaccurate or, for the count, keeps perturbed pivots. The general solver of unsymmetric matrices, scipy.sparse.linalg.eigs(), and the dense solvers are unchanged.

scipy.sparse.linalg.eigsh() uses ARPACK’s symmetric drivers, which return wrong eigenpairs when linked against Intel MKL 2024.2.0 to 2025.0.0 and ncv > 32, see lb() and structsolve.arpackutils, where ncv is capped accordingly.

Non-conservative systems, e.g. under follower loads, may lose stability by flutter, which the static criterion of lb() does not detect: with K the tangent stiffness matrix at a load level, including the load stiffness, the eigenvalues \(\lambda^2\) of the general solvers are real and negative below the flutter load and become complex conjugate pairs above it, when two natural frequencies coalesce, a double root of the frequency equation. This is the dynamic method of [2] (Introduction, Sections 6 and 7, pp. 12-17), where \(\lambda\) is the characteristic exponent \(s = i\Omega\). A real positive \(\lambda^2\), i.e. a real positive \(\lambda\), indicates divergence ([2], p. 17).

References

[1]

Lehoucq, R.B., Sorensen, D.C., Yang, C., ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, SIAM, Philadelphia, 1998. https://doi.org/10.1137/1.9780898719628

[2] (1,2)

Bolotin, V.V., Nonconservative Problems of the Theory of Elastic Stability, Pergamon Press, Oxford, 1963.

See also Eigenvalue solvers for the solver policy.