Linear buckling analysis (`structsolve.linear_buckling`)#

structsolve.lb(K, KG, tol=0, sparse_solver=True, silent=False, num_eigvalues=25, num_eigvalues_print=5, skip_null_cols=False, max_dense_size=2000, check_rtol=0.001, symmetric=None, linear_solver='splu')#

Linear buckling analysis

Calculates the load multipliers lambda and the buckling modes u of the linear buckling problem of K and KG, i.e. the eigenpairs with (K + lambda*KG) @ u = 0, where lambda is the load multiplier of the load state that generated KG. It can also be used for more general eigenvalue analyses, e.g. if K is the tangent stiffness matrix of a given load state. The null rows and columns of K are removed from both matrices before solving the eigenvalue problem, unless skip_null_cols=True.

Solver policy. With sparse_solver=True the sparse solvers run first for every problem size, in the order given below, and each result must pass the verification: the relative residual of each eigenpair must not be larger than check_rtol and, for symmetric matrices with a positive definite K, no positive load multiplier below the highest one returned, and no negative one closer to zero than the returned ones, may be missing, see _check_both_signs(). Only when all of them failed (verification, ARPACK error or no convergence, singular factorization, no positive real load multiplier…) a dense solver is used, if the dense problem has at most max_dense_size dofs, with a DenseFallbackWarning that gives the reasons of the failures. Its result is verified as well. A RuntimeError listing every attempted solver is raised when the problem is too large for the dense fallback, or when the dense solution also failed.

Symmetric matrices. The sparse solvers are:

scipy.sparse.linalg.eigsh() in Cayley mode, with the shift estimated from the matrices by _estimate_sigma();
the same with a 10 times larger shift.

Each one computes the lowest positive load multipliers and, when fewer than num_eigvalues exist, the negative ones closest to zero, i.e. the critical load multipliers of the reversed load, with -KG, see _eigsh_both_signs(). A reference load without any positive load multiplier, e.g. a tensile load, is solved by the sparse solvers too. The dense fallback condenses out the dofs where KG has null rows and solves the condensed problem with scipy.linalg.eigh(), see _eigh_condensed().

Unsymmetric matrices, e.g. KG with the load stiffness of a follower pressure, may have complex load multipliers. The sparse solver is a search of shifts along the positive real axis with scipy.sparse.linalg.eigs(), see _eigs_shift_search(), which finds the lowest positive real load multiplier. The returned eigenpairs are the num_eigvalues load multipliers nearest the last shift of the search, which is zero unless the multipliers of smallest modulus did not include a positive real one, not necessarily the load multipliers of smallest modulus. The dense fallback, _eig_condensed(), keeps the dofs where either the row or the column of KG is not null, and returns the num_eigvalues first load multipliers in the order described in Returns; when none of them is positive and real, they are returned with a warning. Eigenvalues whose imaginary part is below 1e-8 of their modulus are returned as real, with real eigenvectors. The arrays are real when the returned eigenpairs are real, even if other eigenvalues of the solver are complex.

Intel MKL. When SciPy is linked against Intel MKL 2024.2.0 to 2025.0.0, ncv of scipy.sparse.linalg.eigsh() is capped at 32 and more than 15 eigenpairs are computed in several passes, see structsolve.arpackutils, where the cap can be overridden, e.g. with the environment variable STRUCTSOLVE_ARPACK_MAX_NCV=0 to lift it with MKL 2025.0.1 or newer. The dense solvers of unsymmetric matrices call scipy.linalg.eig() with complex matrices.

Parameters:

Ksparse_matrix: Stiffness matrix. Should include all constant terms of the initial stress stiffness matrix, aerodynamic matrix and so forth when applicable.
KGsparse_matrix: Initial stress stiffness matrix that multiplies the load multiplier \(\lambda\) of the eigenvalue problem.
tolfloat, optional: A float tolerance passed to scipy.sparse.linalg.eigsh() and scipy.sparse.linalg.eigs().
sparse_solverbool, optional: With True, the default, the sparse solvers always run first, whatever the size of the problem, and a dense solver is only used as a last resort, see above. num_eigvalues load multipliers are returned. With False, the full problem is explicitly solved with a dense solver, without trying the sparse solvers: scipy.linalg.eigh(), which requires a positive definite K, or scipy.linalg.eig() for unsymmetric matrices, returning all the load multipliers. It is a reference for validations and for the callers that need all of them.
silentbool, optional: A boolean to tell whether the log messages should be printed.
num_eigvaluesint, optional: Number of load multipliers calculated with the sparse solver, and number of returned eigenvectors. Fewer are returned by the sparse solver of symmetric matrices when fewer finite load multipliers exist, the infinite ones, e.g. of the null rows of KG, being discarded.
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 is known to be non-singular.
max_dense_sizeint, optional: Maximum size of the dense problem solved as a last resort, when all the sparse solvers failed, i.e. the number of non-null rows of KG, or of non-null rows or columns for unsymmetric matrices, see above. Larger problems raise a RuntimeError instead, and 0 disables the dense fallback. Until version 0.5.3, it was the size below which the dense solver was tried first.
check_rtolfloat, optional: Maximum relative residual ||K u + lambda KG u|| / (||K u|| + |lambda| ||KG u||) of the eigenpairs of the sparse solver, see above and _relative_residuals().
symmetricbool or None, optional: Whether K and KG are symmetric. None checks it with is_symmetric(), True skips the check and uses the solvers of symmetric matrices, False forces the solvers of unsymmetric matrices, required e.g. when KG includes the load stiffness of follower loads, see above.
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, in the sparse solvers of symmetric matrices only; see the Notes.

Returns:

eigvalsndarray: The load multipliers \(\lambda\), calculated as -1/eigval from the eigenvalues eigval of KG u = eigval K u. The positive load multipliers come first in increasing order, followed by the negative ones in order of increasing modulus, i.e. the critical load multipliers of the reversed load first, for both solvers. For unsymmetric matrices the array is complex when some eigenvalues are complex, which come after the real ones by increasing modulus. The load multipliers of the null eigval, |eigval| not larger than NULL_MU_RTOL times the largest |eigval|, infinite up to the round-off, come last; the dense solvers also take as null the eigval not determined by the round-off, see linear_buckling._null_mu_pencil, returned as -inf.
eigvecsndarray: The \(i^{th}\) eigenvector is eigvecs[:, i], with the size of the original matrices. Only num_eigvalues eigenvectors are returned.

Raises:

RuntimeError: When no eigenvalue solver passed the verification, with the list of the attempted solvers and their failures.

Warns:

DenseFallbackWarning: When the sparse solvers failed and the dense solver was used.
RuntimeWarning: For unsymmetric matrices, when the dense fallback found no positive real load multiplier.

Notes

The eigenvalue problem is

\[(\boldsymbol{K} + \lambda \boldsymbol{K}_G) \boldsymbol{u} = \boldsymbol{0}\]

where \(\lambda\) is the load multiplier of the load state that generated \(\boldsymbol{K}_G\). It is solved as \(\boldsymbol{K}_G \boldsymbol{u} = \mu \boldsymbol{K} \boldsymbol{u}\), with \(\lambda = -1/\mu\); the null eigenvalues \(\mu\), e.g. of the dofs where \(\boldsymbol{K}_G\) is null, give infinite load multipliers.

Solver policy. The dense fallback keeps structsolve robust where ARPACK is unreliable, normal runs never pay for a dense solution.

Symmetric matrices. The shift of the Cayley mode is estimated by power iterations on \(\boldsymbol{K}^{-1} \boldsymbol{K}_G\), see _estimate_sigma() and _eigsh_cayley(). The dense fallback condenses out the dofs where KG has null rows, typically the in-plane dofs, with a sparse factorization of K, see _eigh_condensed(). The eigenpairs are verified as follows:

the relative residual \(||K u + \lambda K_G u|| / (||K u|| + |\lambda| ||K_G u||)\) of each eigenpair must not be larger than check_rtol, see _relative_residuals();
when K is positive definite, no load multiplier below the highest positive one returned may be missing: the number of negative pivots of \(\boldsymbol{K} + s \boldsymbol{K}_G\), for \(s\) slightly below that multiplier, is the number of load multipliers in \((0, s)\) (Sylvester’s law of inertia). The factorization is done with SuperLU, without row interchanges, or with PARDISO, see below. When a row interchange was needed and the count is unknown, \(\boldsymbol{K} + s \boldsymbol{K}_G\) must be positive definite for \(s\) slightly below the lowest positive load multiplier. See _check_eigenpairs();
the same for the negative load multipliers with \(-\boldsymbol{K}_G\), i.e. the positive load multipliers of the reversed load, see _check_both_signs().

The Cayley mode selects all the negative \(\mu\), i.e. the positive load multipliers, before the null and the positive \(\mu\). The load multipliers of each sign are first counted with the inertia of \(\boldsymbol{K} + s \boldsymbol{K}_G\), for \(|s|\) equal to INERTIA_COUNT_FACTOR times the lowest modulus that a load multiplier may have, and at most as many are requested, since ARPACK may not converge when it must return eigenvalues of the cluster of the null \(\mu\). When fewer than num_eigvalues positive load multipliers exist, the Cayley mode with \(-\boldsymbol{K}_G\) gives the negative load multipliers closest to zero, see _eigsh_both_signs().

PARDISO. With linear_solver='pardiso', the factorizations of the sparse solvers of symmetric matrices are done by PARDISO, through static._factorize(): the symmetric indefinite factorization \(\boldsymbol{L} \boldsymbol{D} \boldsymbol{L}^\top\) of the matrix scaled to a unit diagonal, with its fallback on perturbed pivots. ARPACK remains the eigenvalue solver: the operator \((\boldsymbol{K}_G - \sigma \boldsymbol{K})^{-1}\) of the Cayley mode, also of the passes of arpackutils.capped_eigsh(), is passed to scipy.sparse.linalg.eigsh() as OPinv, each solution refined by GMRES when its residual is too large, see arpackutils.pardiso_solve(); K of the power iterations of _estimate_mu_max() and the matrices \(\boldsymbol{K} + s \boldsymbol{K}_G\) of the inertia counts are factorized by PARDISO too. The number of negative pivots is then the number of negative eigenvalues given by PARDISO, exact also when SuperLU would need a row interchange, for which its count is unknown. A factorization of PARDISO is used only when it is accurate: the relative residual of its solution for a random right-hand side must not be larger than arpackutils.PARDISO_ACCURACY_RTOL, and no pivot may remain perturbed for the counts; otherwise SuperLU is used for that matrix, with a RuntimeWarning for the operator of ARPACK, see arpackutils.pardiso_solve() and _negative_pivots(). The symmetric indefinite factorization of PARDISO (MKL 2025.0, default parameters, iparm(2) = 2, iparm(10) = 8) of the Ritz models of the tests, e.g. cylinder_5panels (1140 dofs) and plate_ssss_Nxx (300 dofs), perturbed 5 to 80 pivots and gave residuals of 1 to 50, also for a positive definite K, and its fallback iparm(24) = 1 of static._factorize() reported no perturbed pivot with residuals of 1 to 2e7 and a wrong inertia; for these models 'pardiso' thus falls back to SuperLU. The solvers of unsymmetric matrices and the dense fallbacks, also the sparse factorization of K of _eigh_condensed(), keep SuperLU and LAPACK.

Unsymmetric matrices, e.g. KG with the load stiffness of a follower pressure that is not conservative, may have complex eigenvalues. The search of shifts of _eigs_shift_search() finds the lowest positive real load multiplier even when many negative or complex multipliers have a smaller modulus, verifies the residuals, and checks the parity of the number of real multipliers below the critical one. The dofs kept by the dense fallback, where either the row or the column of KG is not null, are for follower loads almost all of them. When no returned load multiplier is positive and real, no divergence load exists, e.g. for a system that only loses stability by flutter, see freq(). The inertia check requires symmetric matrices and is skipped.

Intel MKL. ARPACK, used by scipy.sparse.linalg.eigsh(), returns wrong eigenpairs at random, without any error, or raises ArpackError -8 when linked against Intel MKL 2024.2.0 to 2025.0.0, e.g. in some Anaconda builds of SciPy: the dsteqr routine of these MKL versions returns wrong eigenvectors for matrices larger than 32 x 32, which ARPACK uses with the size ncv. When such a LAPACK is detected, ncv is capped at 32 and more than 15 eigenpairs are computed in several passes, combined with a Rayleigh-Ritz projection, see structsolve.arpackutils. The non-symmetric drivers of scipy.sparse.linalg.eigs() are not affected.

With Intel MKL 2025.0.0, the real LAPACK drivers dgeev and dggev of scipy.linalg.eig() crash the Python process for some matrices when computing the eigenvectors, with any number of threads. The dense solvers of unsymmetric matrices therefore call scipy.linalg.eig() with complex matrices, i.e. zgeev and zggev, which are not affected.

structsolve.linear_buckling.is_symmetric(A, rtol=1e-10)#

Tell whether ||A - A^T|| <= rtol ||A||, in the Frobenius norm

A cheap check, of the order of the number of non-zero entries of A, used by lb() to select the solvers of symmetric or unsymmetric matrices.

Parameters:

Aarray or sparse matrix: Square matrix, converted to CSR format.
rtolfloat, optional: Relative tolerance, dimensionless. Default 1e-10.

Returns:

symmetricbool: True when the Frobenius norm of A - A^T is not larger than rtol times the Frobenius norm of A.

Linear buckling analysis (structsolve.linear_buckling)#

Linear buckling analysis (`structsolve.linear_buckling`)#