Eigenvalue analyses (`pyfe3d.solver.eigen`)#

Linear buckling, with the constitutive and geometric stiffness matrices, and natural frequency, with the constitutive stiffness and mass matrices, the eigenvalue problems that the element matrices of this package are usually assembled for. Neither problem is safe to hand to scipy.sparse.linalg.eigsh() without preparation:

diagonal_preconditioner() equilibrates the diagonal of the matrices, which does not change the eigenvalues;
estimate_cayley_sigma() estimates, from the matrices themselves, a safe shift for the Cayley mode used for linear buckling;
check_eigenpairs() verifies the result: it checks the residual of every eigenpair and, when the constitutive stiffness matrix is positive definite, that no load multiplier was missed below the highest one found, by counting the negative pivots of K + s*KG just below it, see is_positive_definite().

linear_buckling() and natural_frequency() put them together and are what most users should call. They also solve unsymmetric problems, e.g. with the load stiffness of a follower load or with the aerodynamic matrix of piston theory, and they follow the same solver policy: the sparse solvers run first, their eigenpairs are verified, and a dense solver is used only as a last resort, with a DenseFallbackWarning. complex_frequency() follows the same policy for the free vibrations of damped systems, with the stiffness, damping and mass matrices, a quadratic eigenvalue problem.

Note

ARPACK, used by scipy.sparse.linalg.eigsh(), has been reported to return wrong eigenpairs at random, without raising, when linked against Intel MKL 2024.2.0 to 2025.0.0, which includes some Anaconda builds of SciPy, whenever its parameter ncv exceeds 32. MKL 2025.0.1 or newer is not affected. The faulty versions are detected and ncv is then capped at 32, see pyfe3d.solver.arpackutils. The verification of check_eigenpairs() catches such a failure, which is another reason to keep it on.

Notes#

The eigenvalue problems are linear buckling

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

and natural frequency

\[(\boldsymbol{K}_{C0} - \omega_n^2 \boldsymbol{M}) \boldsymbol{u} = \boldsymbol{0}\]

The reason this module exists is that neither problem is safe to hand to scipy.sparse.linalg.eigsh() without preparation, and the two traps are easy to fall into and silent when they bite:

The matrices are badly scaled. A shell element mixes translations and rotations, so the diagonal of \(\boldsymbol{K}_{C0}\) spans many orders of magnitude, and the drilling stiffness makes that worse: with the fictitious penalty of K6ROT the drilling entries are some six orders of magnitude below the membrane ones, and a condition number of \(10^{14}\) is easy to reach. diagonal_preconditioner() equilibrates the diagonal, which does not change the eigenvalues at all and usually decides whether the iterative solver converges.
The shift of the Cayley mode has to be chosen, not guessed. With eigsh(A=KG, M=K, sigma=sigma, mode='cayley', which='SM') the returned eigenvalues are the ones with the smallest \(|(\mu + \sigma)/(\mu - \sigma)|\), which are the most negative \(\mu\), and therefore the lowest positive load multipliers \(\lambda = -1/\mu\), only while \(\sigma\) is larger than \(|\mu|\) of those critical eigenvalues. Otherwise the eigenvalues closest to \(-\sigma\) come back instead, with no warning. estimate_cayley_sigma() estimates a safe shift from the matrices themselves.

The inertia check of check_eigenpairs(), the number of negative pivots of \(\boldsymbol{K}_{C0} + s \boldsymbol{K}_G\) just below the highest load multiplier found, is the only cheap way to know that no mode was skipped, and it is what distinguishes a converged answer from a plausible-looking one.

Solver policy

With sparse_solver=True the sparse solvers run first for every problem size, and each result must pass the verification. Only when all of them failed, e.g. a failed verification, an ARPACK error or a singular factorization, a dense solver is used, if the dense problem has at most max_dense_size degrees-of-freedom, with a DenseFallbackWarning giving the reasons of the failures; its result is verified as well. A RuntimeError listing every attempted solver is raised otherwise. See _sparse_first().

Order of the load multipliers

linear_buckling() returns the positive load multipliers first, in increasing order, followed by the negative ones closest to zero, i.e. the critical load multipliers of the reversed load, then the complex ones by increasing modulus. For symmetric matrices, the negative ones are computed with \(-\boldsymbol{K}_G\) when fewer than num_eigvalues positive load multipliers exist, see _eigsh_both_signs().

The shift estimate, the verification, the dense fallback, the solvers of unsymmetric matrices, the order of the load multipliers and the cap of ncv depend only on NumPy and SciPy. The dsteqr routine of the Intel MKL versions in the note above returns wrong eigenvectors for matrices larger than 32 by 32, see pyfe3d.solver.arpackutils.

exception pyfe3d.solver.eigen.DenseFallbackWarning#

Issued by linear_buckling(), natural_frequency() and complex_frequency() when the sparse eigenvalue solvers failed and a dense solver was used as a last resort

A subclass of RuntimeWarning, issued with warnings.warn(), whose message gives the reasons of the failures of the sparse solvers. It can be turned into an error with warnings.simplefilter('error', DenseFallbackWarning).

pyfe3d.solver.eigen.check_eigenpairs(K, KG, eigvals, eigvecs, rtol=0.001, min_rel_gap=0.0001, check_inertia=True, K_positive_definite=None)#

Verify the eigenpairs of a linear buckling problem

The eigenpairs, load multipliers and buckling modes, of the linear buckling problem of K and KG, see linear_buckling(), are verified with two checks:

the relative residual of every eigenpair with a finite load multiplier lam and eigenvector u, norm(K @ u + lam*(KG @ u))/(norm(K @ u) + abs(lam)*norm(KG @ u)), see _relative_residuals(), must not exceed rtol;
with check_inertia=True and when K is positive definite, no positive load multiplier is missing below the highest positive one returned, lam_n: the number of negative pivots of K + s*KG, with s = lam_n*(1 - min_rel_gap), must be the number of returned load multipliers in (0, s). When this count is unknown, see _negative_pivots(), the weaker check that K + s*KG with s = lam_1*(1 - min_rel_gap) is positive definite, lam_1 being the lowest positive load multiplier, is done.

The eigenpairs may be real or complex; the inertia check ignores the complex load multipliers and requires symmetric matrices. The eigenpairs with a non-finite load multiplier are skipped. When no load multiplier is finite, the check fails at once. The inertia check is skipped when no load multiplier is positive or when K is not positive definite, see is_positive_definite(); missing negative load multipliers are not detected, see _check_both_signs(). Nothing is raised or printed; the result tells which check failed.

Parameters:

K, KGsparse_matrix: Constitutive and geometric stiffness matrices. Shape (n, n).
eigvalsarray_like: Load multipliers, dimensionless, real or complex. 1-D array of size k.
eigvecsarray_like: Eigenvectors, one per column. Shape (n, k).
rtolfloat, optional: Largest acceptable relative residual, dimensionless. Default 1.e-3.
min_rel_gapfloat, optional: Relative margin below the load multipliers at which the inertia is tested, dimensionless. Default 1.e-4.
check_inertiabool, optional: Whether the inertia check is done, only valid for symmetric matrices. Default True.
K_positive_definitebool, callable or None, optional: Whether K is positive definite, or a function of K returning it, e.g. to verify several solutions with one factorization of K. None uses is_positive_definite(). Default None.

Returns:

errorstr or None: Description of the check that failed, or None when all passed.

See also

linear_buckling

Notes

The eigenpairs are those of \((\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_G) \boldsymbol{u} = \boldsymbol{0}\). The inertia check is the one that catches a solver that converged to an interior part of the spectrum, which is silent otherwise. For symmetric \(\boldsymbol{K}_{C0}\) and \(\boldsymbol{K}_G\) with a positive definite \(\boldsymbol{K}_{C0}\), \(\boldsymbol{K}_{C0} + s \boldsymbol{K}_G\) is congruent to \(\boldsymbol{I} + s \boldsymbol{K}_{C0}^{-1/2} \boldsymbol{K}_G \boldsymbol{K}_{C0}^{-1/2}\), whose eigenvalues are \(1 + s \mu_i = 1 - s/\lambda_i\), with \(\mu_i = -1/\lambda_i\) the eigenvalues of \(\boldsymbol{K}_G \boldsymbol{u} = \mu \boldsymbol{K}_{C0} \boldsymbol{u}\). For \(s > 0\), \(1 - s/\lambda_i < 0\) exactly when \(0 < \lambda_i < s\), hence the number of negative eigenvalues of \(\boldsymbol{K}_{C0} + s \boldsymbol{K}_G\), equal to its number of negative pivots by Sylvester’s law of inertia, is the number of load multipliers in \((0, s)\). This is the Sturm sequence check of the eigenvalue solvers of structural analysis ([1], Section 10.2.2, p. 849, Section 11.4.3, p. 943, and Section 11.6, p. 954), where the number of negative elements of \(\boldsymbol{D}\) in the factorization \(\boldsymbol{L} \boldsymbol{D} \boldsymbol{L}^\top\) of \(\boldsymbol{K} - \mu \boldsymbol{M}\) is the number of eigenvalues of the pencil \((\boldsymbol{K}, \boldsymbol{M})\) smaller than \(\mu\). When only one positive load multiplier was returned, none is expected in \((0, s)\) and the weaker check is done. The gap min_rel_gap keeps \(s\) away from \(\lambda_n\), where \(\boldsymbol{K}_{C0} + s \boldsymbol{K}_G\) is singular.

References

[1]

Bathe, K. J., Finite Element Procedures, 2nd ed., K. J. Bathe, Watertown, MA, 2014, ISBN 978-0-9790049-5-7.

pyfe3d.solver.eigen.complex_frequency(K, C, M, num_eigvalues=6, tol=0, sigma=None, precondition=True, check=True, check_rtol=1e-06, sparse_solver=True, max_dense_size=2000, silent=True)#

Complex eigenvalues of a damped system from the stiffness, damping and mass matrices

Solves the quadratic eigenvalue problem of the free vibrations of a damped system, (lam**2*M + lam*C + K) @ u = 0, whose eigenvalues lam give the damped circular frequencies omegad = abs(lam.imag) and the damping ratios zeta = -lam.real/abs(lam). K, C and M must be real and may be unsymmetric, e.g. K with the aerodynamic matrices of the piston theory and C with its damping matrix, or C with gyroscopic terms. K must have no null row or column, see check_null_rows_cols(), whereas C and M may have null rows and columns. The matrices are optionally preconditioned with diagonal_preconditioner().

The problem is linearized into a generalized eigenvalue problem of size 2*n, the first companion form, solved with scipy.sparse.linalg.eigs() in shift-invert mode, which finds the eigenvalues nearest the shift sigma and factorizes the n by n matrix sigma**2*M + sigma*C + K, see _eigs_qep(). The eigenvalues of real matrices are real or come in complex conjugate pairs, and one eigenvalue per pair, the one with a positive imaginary part, is returned, with its eigenvector. The num_eigvalues eigenvalues nearest the shift are returned, sorted by increasing imaginary part, then by increasing modulus, see _qep_modes(); the real ones, e.g. of overdamped modes, come first.

The relative residual of each eigenpair must not be larger than check_rtol, unless its backward error is at the round-off, see _check_qep(). When the verification fails or the sparse solver raises an error, the dense solver of _dense_qep(), scipy.linalg.eig() on the linearization, is used as a last resort, with a DenseFallbackWarning, see the solver policy of pyfe3d.solver.eigen. Whether an eigenvalue nearer the shift than the returned ones was missed is not verified.

With C null, the eigenvalues are 1j*omegan, with the natural frequencies omegan of natural_frequency().

Parameters:

K, C, Msparse_matrix: Stiffness, damping and mass matrices of the free degrees-of-freedom, i.e. with the constrained ones partitioned out, e.g. KC0uu = KC0[bu, :][:, bu]. K should include the constant terms that apply, e.g. the geometric stiffness matrix of a preload. Real, shape (n, n).
num_eigvaluesint, optional: Number of eigenvalues to compute, one per complex conjugate pair. Default 6.
tolfloat, optional: Tolerance passed to scipy.sparse.linalg.eigs(), 0 meaning machine precision. Default 0.
sigmafloat, complex or None, optional: Shift, in rad/s, near which the eigenvalues are sought; it must not be an eigenvalue. None uses 1., which gives the eigenvalues of lowest modulus, see the Notes. Default None.
preconditionbool, optional: Whether to equilibrate the diagonal of K before solving, the same transformation being applied to C and M. Default True.
checkbool, optional: Whether to verify the eigenpairs of the sparse solver. With False they are returned without verification and without dense fallback. Default True.
check_rtolfloat, optional: Largest acceptable relative residual of the eigenpairs, dimensionless. An eigenpair with a larger one is also accepted when its normwise backward error is not larger than BACKWARD_ERROR_RTOL, the machine precision, see _quadratic_backward_errors(). Default 1.e-6.
sparse_solverbool, optional: With False, the full problem is solved with the dense solver, returning all the finite eigenvalues, one per complex conjugate pair, without verification. Default True.
max_dense_sizeint, optional: Maximum size of the dense problem solved as a last resort, i.e. of the linearization, 2*n; 0 disables the dense fallback. Default 2000.
silentbool, optional: With False, the solver used and the eigenvalues, damped circular frequencies and damping ratios are printed, see pyfe3d.solver.logger. Default True.

Returns:

eigvalsnp.ndarray: Complex eigenvalues lam, in rad/s, with lam.imag >= 0. 1-D array of size k, at most num_eigvalues.
eigvecsnp.ndarray: Complex eigenvectors, one per column, in the same order, of unit Euclidean norm with the component of largest modulus real and positive. Shape (n, k).

Raises:

numpy.linalg.LinAlgError: When K has a null row or column, with the indices of these degrees-of-freedom, see check_null_rows_cols().
ValueError: When K, C or M is complex.
RuntimeError: When check=True and no eigenvalue solver passed the verification.

Warns:

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

Notes

The problem solved is

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

the free vibrations \(\boldsymbol{q} = \boldsymbol{u} e^{\lambda t}\) of \(\boldsymbol{M} \ddot{\boldsymbol{q}} + \boldsymbol{C} \dot{\boldsymbol{q}} + \boldsymbol{K} \boldsymbol{q} = \boldsymbol{0}\). With \(\lambda = -\zeta \omega_n \pm i \omega_d\), \(\omega_d = |Im(\lambda)|\) is the damped circular frequency and \(\zeta = -Re(\lambda)/|\lambda|\) the damping ratio: a single degree-of-freedom, \(m \lambda^2 + c \lambda + k = 0\), has the roots \(\lambda = -\zeta \omega_n \pm i \omega_n \sqrt{1 - \zeta^2}\) for \(\zeta < 1\), with \(\omega_n = \sqrt{k/m} = |\lambda|\) and \(\zeta = c/(2 \sqrt{k m})\). A real \(\lambda\), e.g. of an overdamped mode, gives \(\omega_d = 0\), and a positive \(Re(\lambda)\) a growing motion, i.e. an unstable system, by flutter when \(\lambda\) is complex. \(\zeta\) is not defined for \(\lambda = 0\). With \(\boldsymbol{C} = \boldsymbol{0}\), \(\lambda^2 = -\omega_n^2\) of natural_frequency().

Linearization

With \(\boldsymbol{z} = \{\boldsymbol{u}^\top, \lambda \boldsymbol{u}^\top\}^\top\), the problem is the generalized eigenvalue problem of size \(2n\)

\[\begin{split}\begin{bmatrix} \boldsymbol{0} & \boldsymbol{I} \\ -\boldsymbol{K} & -\boldsymbol{C} \end{bmatrix} \boldsymbol{z} = \lambda \begin{bmatrix} \boldsymbol{I} & \boldsymbol{0} \\ \boldsymbol{0} & \boldsymbol{M} \end{bmatrix} \boldsymbol{z}\end{split}\]

the first companion form \(\boldsymbol{A} - \lambda \boldsymbol{B}\) of [1] (Section 3.4, Eq. (3.6) with \(\boldsymbol{N} = \boldsymbol{I}\), p. 253): the first row is \(\lambda \boldsymbol{u} = \lambda \boldsymbol{u}\), and the second gives the quadratic problem. It is valid for any \(\boldsymbol{C}\), symmetric or not, singular or not, and it requires neither a non-singular \(\boldsymbol{K}\) nor a non-singular \(\boldsymbol{M}\): a singular \(\boldsymbol{M}\), e.g. a lumped mass matrix without rotary inertia, only gives infinite eigenvalues, \(\boldsymbol{B}\) being singular, which are discarded. Its shift-invert operator needs the factorization of the \(n \times n\) matrix \(\boldsymbol{K}_{dyn}(\sigma) = \sigma^2 \boldsymbol{M} + \sigma \boldsymbol{C} + \boldsymbol{K}\) only, see _eigs_qep(); MSC/Nastran uses this linearization with a shift-invert Lanczos algorithm ([1], Section 7.1, p. 278). The symmetric linearizations of symmetric \(\boldsymbol{K}\), \(\boldsymbol{C}\) and \(\boldsymbol{M}\) give pencils that are indefinite, unless the system is overdamped ([1], Section 3.8.2, p. 258), for which the symmetric solver of ARPACK cannot be used, and they exclude the unsymmetric matrices of the piston theory. The eigenvectors \(\boldsymbol{u}\) are taken from the first or the last \(n\) rows of \(\boldsymbol{z}\), see _companion_vectors().

Shift

The default shift \(\sigma = 1\) rad/s is the counterpart of the shift \(\omega_n^2 = -1\) of natural_frequency(), \(\lambda^2 = -\omega_n^2 = 1\): with \(\boldsymbol{C} = \boldsymbol{0}\) the factorized matrix is the same, \(\boldsymbol{K} + \boldsymbol{M}\), and \(|\lambda - 1|^2 = 1 + \omega_n^2\) for \(\lambda = \pm i \omega_n\), such that the lowest natural frequencies are found first. With \(\lambda = \omega_n(-\zeta \pm i \sqrt{1 - \zeta^2})\), \(|\lambda - 1|^2 = 1 + 2 \zeta \omega_n + \omega_n^2\), which increases with \(\omega_n\) for \(\omega_n > -\zeta\). When \(\boldsymbol{M}\) is non-singular, \(\det \boldsymbol{K}_{dyn}(\sigma) = \det \boldsymbol{M} \prod_i (\sigma - \lambda_i)\) over the \(2n\) eigenvalues, such that \(\boldsymbol{K}_{dyn}(1)\) is singular only when \(\lambda = 1\) is an eigenvalue, i.e. never for the rigid-body modes, \(\lambda = 0\), nor for a stable system, \(Re(\lambda) \le 0\). A complex shift gives a complex operator, and the conjugate of each eigenvalue found with \(Im(\lambda) < 0\) is returned ([1], Section 6.1, p. 273).

Verification

The relative residual \(\|(\lambda^2 \boldsymbol{M} + \lambda \boldsymbol{C} + \boldsymbol{K}) \boldsymbol{u}\| / (|\lambda|^2 \|\boldsymbol{M} \boldsymbol{u}\| + |\lambda| \, \|\boldsymbol{C} \boldsymbol{u}\| + \|\boldsymbol{K} \boldsymbol{u}\|)\) of each eigenpair is checked, the rigid-body modes being relative to the matrix-norm scale, see _quadratic_residuals(). Unlike natural_frequency(), no inertia check tells that no eigenvalue nearer the shift was missed: by Sylvester’s law of inertia, the negative pivots of \(\boldsymbol{K} - s \boldsymbol{M}\) count the real eigenvalues below \(s\) of a symmetric pencil with a positive definite \(\boldsymbol{M}\), whereas the eigenvalues of a damped system are complex, and the symmetric linearizations are definite only for overdamped systems, whose eigenvalues are real ([1], Sections 3.8.2 and 3.9, p. 258). No count of the eigenvalues of a quadratic eigenvalue problem inside a region of the complex plane is implemented, and that ARPACK found the eigenvalues of largest \(|1/(\lambda - \sigma)|\) is therefore not verified; the dense solution, sparse_solver=False, is the reference. With \(\boldsymbol{C} = \boldsymbol{0}\) and symmetric matrices, natural_frequency() verifies it with the inertia check.

A rigid-body mode, \(\boldsymbol{K} \boldsymbol{u} = \boldsymbol{0}\), gives \(\lambda = 0\), a double eigenvalue when \(\boldsymbol{C} \boldsymbol{u} = \boldsymbol{0}\), which the solvers return as two eigenvalues null up to the round-off.

References

[1] (1,2,3,4,5)

Tisseur, F. and Meerbergen, K., The quadratic eigenvalue problem, SIAM Review, 43(2), 235-286, 2001. https://doi.org/10.1137/S0036144500381988

pyfe3d.solver.eigen.diagonal_preconditioner(K, floor=1e-30)#

Jacobi preconditioner of a stiffness matrix

Returns the diagonal matrix P with the entries 1/sqrt(max(abs(K[i, i]), floor)), such that P @ K @ P has a diagonal of ones, or minus ones for negative entries, when the moduli of the diagonal entries of K are larger than floor. For a generalized eigenvalue problem both matrices are transformed, P @ A @ P and P @ B @ P, which leaves the eigenvalues unchanged, and the eigenvectors of the original problem are recovered with u = P @ u_scaled.

Parameters:

Ksparse_matrix: Matrix whose diagonal sets the scaling, normally the constitutive stiffness matrix with the boundary conditions already applied. Shape (n, n).
floorfloat, optional: Lower bound for the diagonal entries, which guards against the null diagonal of an unconstrained degree-of-freedom. Same dimension as the diagonal entries of K. Default 1.e-30.

Returns:

Psparse_matrix: Diagonal preconditioner, created with scipy.sparse.diags(), shape (n, n).

Notes

The diagonal matrix \(\boldsymbol{P}\) has the entries

\[P_{ii} = \frac{1}{\sqrt{max(|K_{{C0}_{ii}}|, floor)}}\]

such that \(\boldsymbol{P} \boldsymbol{K}_{C0} \boldsymbol{P}\) has a unit diagonal, up to the sign of the negative entries, which appear e.g. in a stiffness matrix with the load stiffness of a follower load beyond the critical load. For a generalized eigenvalue problem both matrices are transformed, \(\boldsymbol{P} \boldsymbol{A} \boldsymbol{P}\) and \(\boldsymbol{P} \boldsymbol{B} \boldsymbol{P}\), which leaves the eigenvalues unchanged and maps the eigenvectors to \(\boldsymbol{P}^{-1} \boldsymbol{u}\): \(\boldsymbol{P} \boldsymbol{A} \boldsymbol{P} \boldsymbol{y} = w \boldsymbol{P} \boldsymbol{B} \boldsymbol{P} \boldsymbol{y}\) is \(\boldsymbol{A} \boldsymbol{u} = w \boldsymbol{B} \boldsymbol{u}\) with \(\boldsymbol{u} = \boldsymbol{P} \boldsymbol{y}\), \(\boldsymbol{P}\) being non-singular. The equilibration removes the spread of the diagonal entries, which usually reduces the condition number, see the notes of pyfe3d.solver.eigen.

Examples

>>> P = diagonal_preconditioner(Kuu)
>>> eigvals, eigvecs = eigsh(A=P@KGuu@P, M=P@Kuu@P)
>>> eigvecs = P @ eigvecs

pyfe3d.solver.eigen.estimate_cayley_sigma(K, KG, safety=10.0, max_iter=50, rel_tol=0.001)#

Safe shift for the Cayley mode of scipy.sparse.linalg.eigsh()

The largest modulus of the eigenvalues mu of the generalized eigenvalue problem KG @ u = mu * K @ u is estimated with power iterations, and the result is multiplied by safety.

The fallback value 1. is returned when K is singular or not positive definite, or when the linear solution is inaccurate.

K is factorized with scipy.sparse.linalg.splu() and the power iteration starts from a vector drawn from numpy.random.RandomState with the seed 42, such that the result is reproducible. The fallback value 1. is also returned when the ratio of norms is not finite or zero, and when any exception is raised during the estimate.

Parameters:

K, KGsparse_matrix: Constitutive and geometric stiffness matrices, with the boundary conditions already applied. Shape (n, n).
safetyfloat, optional: Multiplier applied to the estimated largest abs(mu). Dimensionless. Default 10..
max_iterint, optional: Maximum number of power iterations. Default 50.
rel_tolfloat, optional: Relative increment below which the power iteration is stopped. Dimensionless. Default 1.e-3.

Returns:

sigmafloat: Shift to be passed to eigsh(..., sigma=sigma, mode='cayley'), of the dimension of the eigenvalues mu.

Notes

The shift has to exceed \(|\mu|\) of the critical eigenvalues of \(\boldsymbol{K}_G \boldsymbol{u} = \mu \boldsymbol{K}_{C0} \boldsymbol{u}\), see the module documentation. The largest \(|\mu|\) is estimated with power iterations on \(\boldsymbol{K}_{C0}^{-1} \boldsymbol{K}_G\), whose ratio of norms in the \(\boldsymbol{K}_{C0}\)-norm grows towards it, and the result is multiplied by safety.

At each iteration \(\boldsymbol{x}\) is normalized such that \(\boldsymbol{x}^\top \boldsymbol{K}_{C0} \boldsymbol{x} = 1\), \(\boldsymbol{y} = \boldsymbol{K}_{C0}^{-1} \boldsymbol{K}_G \boldsymbol{x}\) is computed and the ratio of norms is \(\sqrt{|\boldsymbol{y}^\top \boldsymbol{K}_{C0} \boldsymbol{y}|}\); \(\boldsymbol{y}\) becomes the next \(\boldsymbol{x}\). The estimate is the largest ratio found, and the iterations stop when a ratio exceeds the previous estimate by at most rel_tol times the ratio, or after max_iter iterations. For a symmetric positive definite \(\boldsymbol{K}_{C0}\) and a symmetric \(\boldsymbol{K}_G\), \(\boldsymbol{K}_{C0}^{-1} \boldsymbol{K}_G\) is self-adjoint in the inner product of \(\boldsymbol{K}_{C0}\), such that the ratio tends to its spectral radius \(\max |\mu|\). At the first iteration the relative residual of the linear solution must not exceed 1e-6.

The fallback value 1. is returned when \(\boldsymbol{K}_{C0}\) is singular or not positive definite, or when the linear solution is inaccurate, since in those cases the power iteration says nothing.

pyfe3d.solver.eigen.is_positive_definite(A)#

Tells whether the symmetric sparse matrix A is positive definite

A is factorized with a symmetric permutation and without row interchanges, and it is considered positive definite when all the pivots of the factorization are positive.

The factorization uses scipy.sparse.linalg.splu() with the ordering 'MMD_AT_PLUS_A', diag_pivot_thresh=0. and the symmetric mode. The result is False when the factorization fails because A is exactly singular, when a row interchange was nevertheless made, or when a pivot is not positive. The symmetry of A is not checked.

Parameters:

Asparse_matrix: Symmetric matrix, shape (n, n).

Returns:

positive_definitebool: True when all the pivots are positive.

See also

check_eigenpairs

Notes

The Gaussian elimination of a positive definite matrix never meets a zero or a negative pivot, and by Sylvester’s law of inertia a negative pivot means a negative eigenvalue: without row interchanges, the factorization of the symmetric matrix is \(\boldsymbol{L} \boldsymbol{D}_U \boldsymbol{L}^\top\), with the pivots on the diagonal \(\boldsymbol{D}_U\) of \(\boldsymbol{U}\), and \(\boldsymbol{A}\) has as many negative eigenvalues as \(\boldsymbol{D}_U\) has negative entries, see _negative_pivots().

pyfe3d.solver.eigen.linear_buckling(K, KG, num_eigvalues=6, tol=0, sigma=None, precondition=True, check=True, check_rtol=0.001, symmetric=None, sparse_solver=True, max_dense_size=2000)#

Linear buckling analysis with a verified spectrum

Computes the load multipliers and the buckling modes of the linear buckling problem of K and KG, the lowest positive load multipliers being sought. K must have no null row or column, see check_null_rows_cols(), whereas KG may have null rows and columns. K and KG are symmetric when symmetric=None finds them so, see is_symmetric().

Symmetric matrices. The sparse solvers are scipy.sparse.linalg.eigsh() in Cayley mode, using the shift of estimate_cayley_sigma(), or sigma when given, and then 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(). The eigenpairs are verified with the residual and inertia checks of check_eigenpairs(), for both signs, see _check_both_signs(). The dense fallback condenses out the degrees-of-freedom 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 load, 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. Only the residuals are verified. The dense fallback, _eig_condensed(), keeps the degrees-of-freedom where either the row or the column of KG is not null; when none of its load multipliers is positive and real, they are returned with a RuntimeWarning. Eigenvalues whose imaginary part is below 1e-8 of their modulus are returned as real, with real eigenvectors.

The sparse solvers run first and the dense solver is a last resort, see the solver policy of pyfe3d.solver.eigen.

Parameters:

K, KGsparse_matrix: Constitutive and geometric stiffness matrices of the free degrees-of-freedom, i.e. with the constrained ones partitioned out, e.g. KC0uu = KC0[bu, :][:, bu]. K should include the constant terms that apply, e.g. the geometric stiffness matrix of a preload, and KG is the matrix multiplied by the load multiplier. Shape (n, n).
num_eigvaluesint, optional: Number of load multipliers to compute, at most n - 2 with the sparse solvers. Default 6.
tolfloat, optional: Tolerance passed to scipy.sparse.linalg.eigsh() and scipy.sparse.linalg.eigs(), 0 meaning machine precision. Default 0.
sigmafloat or None, optional: Shift of the Cayley mode of the symmetric sparse solver. With None it is estimated from the matrices, which is the recommended use. Default None.
preconditionbool, optional: Whether to equilibrate the diagonal before solving, see diagonal_preconditioner(). This does not change the eigenvalues. Default True.
checkbool, optional: Whether to verify the eigenpairs of the sparse solvers. With False, the first sparse solver is returned without verification and without dense fallback. Default True.
check_rtolfloat, optional: Largest acceptable relative residual of the eigenpairs, dimensionless. Default 1.e-3.
symmetricbool or None, optional: Whether K and KG are symmetric. None checks it with is_symmetric(), True uses the solvers of symmetric matrices, False the solvers of unsymmetric matrices. Default None.
sparse_solverbool, optional: With True the sparse solvers run first, the dense one being a last resort. With False, the full problem is explicitly solved with a dense solver, scipy.linalg.eigh(), which requires a positive definite K, or scipy.linalg.eig() for unsymmetric matrices, returning all the load multipliers, without verification; it is a reference for validations. Default True.
max_dense_sizeint, optional: Maximum size of the dense problem solved as a last resort, the number of non-null rows of KG, or of non-null rows or columns for unsymmetric matrices; 0 disables the dense fallback. Default 2000.

Returns:

eigvalsnp.ndarray: Load multipliers, dimensionless: the positive ones first in increasing order, followed by the negative ones in order of increasing modulus, then the complex ones by increasing modulus, see _lb_order(). Fewer than num_eigvalues are returned by the symmetric sparse solvers when fewer finite load multipliers exist, the infinite ones, e.g. of the null rows of KG, being discarded. The dense solvers return as -inf, last, the load multipliers of the eigenvalues null up to the round-off, see _null_mu_pencil(). Complex only when complex load multipliers are returned.
eigvecsnp.ndarray: Eigenvectors, one per column, in the same order. Shape (n, k).

Raises:

numpy.linalg.LinAlgError: When K has a null row or column, with the indices of these degrees-of-freedom, see check_null_rows_cols().
RuntimeError: When check=True and 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 problem solved is \((\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_G) \boldsymbol{u} = \boldsymbol{0}\), written as \(\boldsymbol{K}_G \boldsymbol{u} = \mu \boldsymbol{K}_{C0} \boldsymbol{u}\) with \(\lambda = -1/\mu\), whose lowest positive load multipliers are the most negative \(\mu\), returned by the Cayley mode while the shift is large enough, see the notes of pyfe3d.solver.eigen. With precondition=True the problem solved is the one of \(\boldsymbol{P} \boldsymbol{K}_{C0} \boldsymbol{P}\) and \(\boldsymbol{P} \boldsymbol{K}_G \boldsymbol{P}\), with the same eigenvalues, see diagonal_preconditioner(); the verification always uses the original matrices.

A null row or column of \(\boldsymbol{K}_{C0}\) is a degree-of-freedom without stiffness, a modelling error: with the same null row or column of \(\boldsymbol{K}_G\), \(\det(\boldsymbol{K}_{C0} + \lambda \boldsymbol{K}_G) = 0\) for every \(\lambda\), a singular pencil whose computed eigenvalues are arbitrary, and otherwise \(\lambda = 0\) is a spurious load multiplier. It is therefore reported instead of being removed, unlike in structsolve, see check_null_rows_cols(). A null row of \(\boldsymbol{K}_G\), at a degree-of-freedom that the load does not destabilize, only gives an infinite load multiplier.

A failure of the verification usually means that the mesh, and not the solver, needs attention: a spectrum with many nearly coincident eigenvalues, as in the buckling of a cylindrical shell, is the typical case.

Non-conservative systems, e.g. under follower loads, may lose stability by flutter, which this static criterion does not detect, see natural_frequency(). When no returned load multiplier is positive and real, no divergence load was found.

pyfe3d.solver.eigen.natural_frequency(K, M, num_eigvalues=6, tol=0, precondition=True, symmetric=None, check=True, check_rtol=1e-06, sparse_solver=True, max_dense_size=2000)#

Natural frequencies from the stiffness and mass matrices

K must have no null row or column, see check_null_rows_cols(), whereas M may have null rows and columns, e.g. a lumped mass matrix without rotary inertia. The matrices are optionally preconditioned with diagonal_preconditioner(), and the eigenvalues, the squares of the circular frequencies, are computed in shift-invert mode with sigma=-1. and which='LM', with scipy.sparse.linalg.eigsh() for symmetric matrices and scipy.sparse.linalg.eigs() otherwise, e.g. with the aerodynamic matrix of piston theory or the load stiffness of a follower load in K. K and M are symmetric when symmetric=None finds them so, see is_symmetric(). The number of computed frequencies is k = min(num_eigvalues, n - 2).

The eigenpairs are verified, see _check_freq(): the relative residual of each eigenpair must not be larger than check_rtol, unless its backward error is at the round-off, and, for symmetric matrices, no eigenvalue lower than the highest one returned may be missing. When the verification fails or the sparse solver raises an error, the dense solver of _dense_freq() is used as a last resort, with a DenseFallbackWarning, see the solver policy of pyfe3d.solver.eigen.

The eigenpairs are sorted together by increasing frequency, the eigenvectors are transformed back to the original degrees-of-freedom, and negative eigenvalues of symmetric matrices give a zero frequency.

Parameters:

K, Msparse_matrix: Constitutive stiffness and mass matrices of the free degrees-of-freedom, i.e. with the constrained ones partitioned out, e.g. KC0uu = KC0[bu, :][:, bu]. K should include the constant terms that apply, e.g. the geometric stiffness matrix of a preload. Shape (n, n).
num_eigvaluesint, optional: Number of frequencies to compute. Default 6.
tolfloat, optional: Tolerance passed to scipy.sparse.linalg.eigsh() or scipy.sparse.linalg.eigs(), 0 meaning machine precision. Default 0.
preconditionbool, optional: Whether to equilibrate the diagonal of K before solving. Default True.
symmetricbool or None, optional: Whether K and M are symmetric. None checks it with is_symmetric(). Default None.
checkbool, optional: Whether to verify the eigenpairs of the sparse solver. With False they are returned without verification and without dense fallback. Default True.
check_rtolfloat, optional: Largest acceptable relative residual of the eigenpairs, dimensionless. An eigenpair with a larger one is also accepted when its normwise backward error is not larger than BACKWARD_ERROR_RTOL, the machine precision, e.g. an eigenvalue close to zero of a model pre-stressed close to buckling, see _backward_errors(). Default 1.e-6.
sparse_solverbool, optional: With False, the full problem is solved with a dense solver, returning all the frequencies, without verification. Default True.
max_dense_sizeint, optional: Maximum size of the dense problem solved as a last resort; 0 disables the dense fallback. Default 2000.

Returns:

omegannp.ndarray: Circular frequencies in rad/s, in increasing order. 1-D array of size k. For unsymmetric matrices whose eigenvalues \(\omega_n^2\) are complex, the complex \(\omega_n = \sqrt{\omega_n^2}\) (principal root), sorted by increasing real part of \(\omega_n^2\), then by increasing imaginary part.
eigvecsnp.ndarray: Eigenvectors, one per column, in the same order. Shape (n, k).

Raises:

numpy.linalg.LinAlgError: When K has a null row or column, with the indices of these degrees-of-freedom, see check_null_rows_cols().
RuntimeError: When check=True and no eigenvalue solver passed the verification.

Warns:

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

Notes

The problem solved is \((\boldsymbol{K}_{C0} - \omega_n^2 \boldsymbol{M}) \boldsymbol{u} = \boldsymbol{0}\), and the frequencies of symmetric matrices are \(\omega_n = \sqrt{\max( \omega_n^2, 0)}\). In the shift-invert mode, the sparse solvers find the eigenvalues closest to the shift sigma=-1., i.e. the lowest ones when \(\boldsymbol{K}_{C0}\) is positive semi-definite and \(\boldsymbol{M}\) positive definite, factorizing \(\boldsymbol{K}_{C0} + \boldsymbol{M}\), which is not singular when \(\boldsymbol{K}_{C0}\) is singular only because of unconstrained rigid-body modes that have mass.

A null row or column of \(\boldsymbol{K}_{C0}\) is a degree-of-freedom without stiffness, a modelling error: with the same null row or column of \(\boldsymbol{M}\), \(\det(\boldsymbol{K}_{C0} - \omega_n^2 \boldsymbol{M}) = 0\) for every \(\omega_n^2\), a singular pencil whose computed eigenvalues are arbitrary, and otherwise \(\omega_n^2 = 0\) is a spurious frequency. It is therefore reported instead of being removed, unlike in structsolve, see check_null_rows_cols().

A mass matrix without rotary inertia, and one without inertia associated with the drilling rotation, is singular, so the eigenvalues associated with those degrees-of-freedom are infinite. The shift-invert mode used here looks for the eigenvalues closest to zero and is not troubled by them, but a lumped mass matrix will place the massless degrees-of-freedom wherever the factorization puts them, so prefer a consistent mass matrix when the lowest frequencies matter.

Non-conservative systems, e.g. under follower loads, may lose stability by flutter, which the static criterion of linear_buckling() does not detect: with K the tangent stiffness matrix at a load level, including the load stiffness, the eigenvalues \(\omega_n^2\) are real and positive 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 [1] (Introduction, Section 6, pp. 12-15, and Section 7, pp. 15-17), where \(\omega_n\) is the frequency \(\Omega\) of the characteristic exponent \(s = i\Omega\). A real negative \(\omega_n^2\), i.e. a real positive \(s\), indicates divergence ([1], p. 17).

References

[1] (1,2)

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

Eigenvalue analyses (pyfe3d.solver.eigen)#

Notes#

Eigenvalue analyses (`pyfe3d.solver.eigen`)#