Adaptation of the Lanczos Algorithm for the Solution of Buckling Eigenvalue Problems

Abstract

An adaptation of the conventional Lanczos algorithm is proposed to solve the general symmetric eigenvalue problem Kϕ = λK _Gϕ in the case when the geometric stiffness matrix K_G is not necessarily positive-definite. The only requirement for the new algorithm to work is that matrix K must be positive-definite. Firstly, the algorithm is presented for the standard situation where no shifting is assumed. Secondly, the algorithm is extended to include shifting since this procedure may be important for enhanced precision or acceleration of convergence rates. Neither version of the algorithm requires matrix inversion, but more resources in terms of memory allocation are needed by the version with shifting.

Keywords
Lanczos algorithm; general eigenproblem; shear buckling

1 INTRODUCTION

One of the most efficient methods for extracting eigenpairs (eigenvalues and eigenvectors) was proposed by Lanczos ( Lanczos, 1950 ). Lanczos algorithm is considered today to be the most efficient (if not the most efficient) numerical method for computing extreme eigenvalues and associated eigenvectors ( Paige, 1971 ; Paige, 1972 ) for large problems. The finite arithmetic of computers introduce some difficulties in the procedure since it may lead to loss of orthogonality of the eigenvectors being computed. However, the loss of orthogonality can be treated through corrective schemes ( Parlett and Scott, 1979 ) already available.

In structural modal analysis the eigenvalue problem to be solved is usually not in the standard form Aϕ = λϕ, where A is a symmetric matrix. Because of inertial effects, in addition to the stiffness matrix K, there is always a mass matrix M involved, such that the eigenvalue problem is stated as Kϕ = λ Mϕ, known as the generalized form. In modal analyses the stiffness matrix is guaranteed to be positive semi-definite, whereas the mass matrix is positive-definite ( Ericsson and Ruhe, 1980 ; Ramaswamy, 1980 ; Bathe 1996 ). Nevertheless, the traditional linearized buckling analysis requires computation of a geometric stiffness matrix K_G, which is dependent on the prebuckling stress distribution ( Zienkiewicz, 1991 ). In the linearized buckling analysis matrix K_G replaces matrix M. However, the positive-definiteness of K_G is not guaranteed what brings apparently unsurmountable an obstacle to the conventional Lanczos algorithm applied to the generalized eigenproblem. Jones and Patrick ( Jones and Patrick, 1989 ) propose the use of a spectral transformation to solve buckling eigenproblems, but it relies on a shifted operator and does not address the basic difficulty of taking the square root of a possibly negative number, except by taking preventive action of selecting appropriate initial guess vectors that must be preconditioned. However, if K_G is indefinite, their method crashes.

Given the incapacity of the conventional Lanczos algorithm to handle indefinite matrices K_G this work proposes an adaptation of the algorithm to compute eigenvalues (critical loads) in linearized buckling problems. The strategy consists basically in relying on the positive-definiteness of K, despite the fact that K_G is indefinite. Initially, a version of the adapted algorithm is proposed without shifting that is sufficient to obtain the lowest positive eigenvalues usually sought. Subsequently, a version of the algorithm with shifting is proposed in order to either enhance precision or accelerate rates of convergence. It is shown that the version with shifting requires more computer resources in terms of memory allocation.

2 THE ADAPTED LANCZOS ALGORITHM

The generalized linear eigenvalue buckling problem is

$K\varphi =\lambda K_G\varphi$ (1)

where K is the stiffness matrix, ϕ is the eigenvector, λ is the eigenvalue and K_G is the geometric stiffness matrix. K is symmetric and positive-definite. K_G is symmetric but may not be positive-definite. In classical linearized buckling analysis when there is only compressive prebuckling loads matrix K_G can be assumed to be positive-definite. However, this is not always the case. Structural components subject to shear loads invariably lead to indefinite K_G matrices. Two simple examples of these situation are panels subject to shear and closed box beams subject to torsion, both of which consist of practical cases encountered in aircraft wing design. Nonetheless, even compressive loads may result in indefinite K_G matrices. Almeida and Hansen (2000) investigated rectangular composite plates with circular cutouts and observed local buckling modes related to negative eigenvalues. The absolute values of those negative eigenvalues were high in comparison to the positive eigenvalues. Even though these negative eigenvalues are not of practical relevance, they cripple the traditional Lanczos algorithm that will eventually be required to compute the square root of a negative number as will be shown. Figure 1 presents the traditional Lanczos algorithm proposed by Bathe ( Bathe, 1996 ) to solve the generalized eigenvalue posed in Eq. (1) .

Figure 1
Algorithm 1.

The algorithm in Fig. 1 generates a sequence of vectors ϕ_i, at least theoretically, orthonormal to K_G. As q in algorithm 1 increases, the symmetric tridiagonal matrix

$T=\left[\begin{array}{ccc}\begin{array}{ccc}{\alpha }_1& {\beta }_2& \\ {\beta }_2& {\alpha }_2& {\beta }_3\\ & & \end{array}& \begin{array}{cc}& \\ & \\ \ddots & \end{array}& \begin{array}{cc}& \\ & \\ & \end{array}\\ \begin{array}{ccc}& & \\ & & \end{array}& \begin{array}{cc}& {\beta }_{q-1}\\ & \end{array}& \begin{array}{cc}{\alpha }_{q-1}& {\beta }_q\\ {\beta }_q& {\alpha }_q\end{array}\end{array}\right]$ (2)

will possess eigenvalues that converge to the inverse of the eigenvalues of the original problem. In theory, it is shown that, when q is exactly the dimension of the eigenvalue problem, matrix T given in Eq. (2) will have eigenvalues that are exactly the inverse of those of the original eigenproblem.

Observe that, in step 3 of algorithm 1, matrix products like ϕ^TK_Gϕ must be computed and, subsequently, their square roots must be taken. In the case of buckling problems K_G is indefinite, therefore, ϕ^TK_G ϕ may be negative, fatally crippling the algorithm. If one insists in using the conventional algorithm and takes, instead, the square root of the absolute value of ϕ^TK_Gϕ, the algorithm flies off and converges to absolutely meaningless results. If the inverted eigenvalue problem K_G ϕ = υKϕ is solved instead, where υ = 1/λ, matrix products of the form ϕ^TKϕ would be computed and, because K is positive-definite, it would lead necessarily to a positive number. Since the conventional Lanczos algorithm computes approximated values to the lowest eigenvalues, the lowest υ’s are obtained and they correspond to the highest λ’s (recall that υ = 1/λ). Hence, solving the inverted eigenproblem is not an option if the lowest λ’s are to be obtained.

The fundamental aspect to be noted in proposing a new variant of the Lanczos algorithm is to realize that, since K is positive-definite, its inverse, K⁻¹ , is also positive-definite. One can therefore define the transformation

$\varphi =K^{-1}x$ (3)

that is always unique because K is positive-definite, therefore invertible. Substitution of Eq. (3) into Eq. (1) and pre-multiplication by (K_G)⁻¹ yields

$K_G^{-1}x=\lambda K^{-1}x$ (4)

There is apparently one inconsistency in Eq. (2): K_G is indefinite and may be singular, such that (K_G)⁻¹ may not even exist. Let us forget this inconsistency for now and employ the conventional Lanczos algorithm to the eigenproblem posed in Eq. (4) .

The difficulty of having to take the square root of a negative number is now resolved since the product x^TK⁻¹ x results in a positive number. However, the algorithm presented in Fig. 2 is apparently highly awkward because both K⁻¹ and ( K_G) ⁻¹ are required. If one carefully examines the algorithm it will be clear that inversions are unnecessary. The inverse of K_G is required only in the evaluation of $K_G^{-1}{\overline{x}}_{i+1}=K^{-1}{x}_i$ that can also be recast as ${\overline{x}}_{i+1}=K_G{K}^{-1}{x}_i$ . Therefore, (K_G) ⁻¹ is not actually needed. As for K⁻¹, notice that only products of the form K⁻¹x are required. These can be more effectively evaluated if recast as Kϕ = x. Numerical solution of this system is achieved through traditional LU decomposition of K (Crout, Cholesky, etc.). Before beginning algorithm 2, matrix K can be LU decomposed and, since it is symmetric, only L must be stored. Refinement of algorithm 2 results in algorithm 3 presented in Fig. 3 .

Figure 2
Algorithm 2

Figure 3
Algorithm 3

Vector v is only auxiliary. Computation of v in step 5 ( $v=K^{-1}{\tilde{x}}_{i+1}$ ), is obtained through solution of the system $Kv={\tilde{x}}_{i+1}$ , which is not troublesome since matrix K has already been LU decomposed. A particularly effective way of selecting the starting vector ϕ ₁ is to make it equal to the main diagonal of matrix K. Since the transformation proposed in Eq. (3) was adopted, the eigenvectors of the original and transformed problems are related by ϕ_i = K⁻¹ x_i.

The eigenproblem with shifting is posed in Eq. (5)

$\left(K-\sigma K_G\right)\varphi =\mu K_G\varphi$ (5)

where σ is the shift and the eigenvalues µ to be determined satisfy µ = λ − σ. If the transformation proposed in Eq. (3) is used in Eq. (5) , pre-multiplication by (K_G)⁻¹ leads to

$K_G^{-1}\left(K-\sigma K_G\right)K^{-1}x=\mu K^{-1}x$ (6)

The particularly complicated matrix product $K_G^{-1}\left(K-\sigma K_G\right)K^{-1}$ is purposely not computed. Both K⁻¹ and (K _G)⁻¹ in this product are not actually evaluated in the algorithm proposed in Fig. 4 . However, matrix, (K − σK _G), is now present and it must be LU decomposed. The adapted Lanczos algorithm with shifting is:

Figure 4
Algorithm 4

In algorithm 4 it is suggested that the two sequences of eigenvectors x_i and ϕ_i are stored. This procedure precludes the necessity to solve systems of the form Kϕ = x. However, if computer memory is scarce, then either x or ϕ only may be stored. Storing the sequence of ϕ_i makes more sense since these are the eigenvectors of the original eigenvalue problem. Computer memory may also become an issue because, in addition to K and K_G, matrix K − σK_G decomposed must also be stored. Another source of hardship may be the decomposition of K − σK_G. If the shift σ is near one of the eigenvalues of the problem stated in Eq. (1) , K − σK_G may become singular, therefore not invertible.

3 NUMERICAL EXAMPLES

A simple eigenvalue problem was proposed by Ramaswamy (1980): find all the eigenvalues of Eq. (1) for matrices K and K_G given by

$\begin{array}{ccc}K=\left[\begin{array}{cc}\begin{array}{ccc}1& 0& 0\\ 0& 3& 0\\ 0& 0& 5\end{array}& \begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{cc}4& 0\\ 0& 2\end{array}\end{array}\right]& and& K_G=\left[\begin{array}{cc}\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& -1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}\end{array}\right]\end{array}$ (7)

Since K_G is indefinite, this simple problem helps illustrating the obstacles encountered by the conventional Lanczos algorithm presented in Fig. 1 . The starting vector is chosen as

$\varphi ={\left\{\begin{array}{cc}\begin{array}{ccc}1/\sqrt{5}& 1/\sqrt{5}& 1/\sqrt{5}\end{array}& \begin{array}{cc}1/\sqrt{5}& 1/\sqrt{5}\end{array}\end{array}\right\}}^T$ (8)

Blunt application of algorithm 1, requires, in the first iteration, evaluation of β₂² = −0.1181, halting the algorithm. Stubbornly taking the square root of the absolute value of ${({\tilde{\varphi }}_{i+1}^T{K}_G{\tilde{\varphi }}_{i+1})}^{1/2}$ to compute β_i+1 does not help; it does not fix the algorithm. In this simple example this procedure would lead to the wrong set of eigenvalues 1.3409, 0.3819, −0.8515, 0.0392 and −0.0809. On the other hand, algorithm 3 produces the correct eigenvalues and eigenvectors placed columnwise:

$\begin{array}{ccc}\left\{\begin{array}{c}1.0000\\ 2.0000\\ -5.0000\\ 4.0000\\ 3.0000\end{array}\right\}& and& \left[\begin{array}{ccccc}1.0000& 0.0000& 0.0003& 0.0006& 0.0005\\ -0.0003& 0.0000& 0.0000& 0.0000& 0.5774\\ -0.0002& 0.0002& 0.4472& 0.0000& 0.0000\\ -0.0003& 0.0006& 0.0000& 0.5000& 0.0000\\ 0.0000& -0.7071& 0.0003& 0.0008& -0.0001\end{array}\right]\end{array}$ (9)

If one makes the geometric stiffness matrix K_G indefinite by, for instance assuming that

$K_G=\left[\begin{array}{cc}\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& -1\end{array}& \begin{array}{cc}0& 0\\ 0& 0\\ 0& 0\end{array}\\ \begin{array}{ccc}0& 0& 0\\ 0& 0& 0\end{array}& \begin{array}{cc}1& 0\\ 0& 1\end{array}\end{array}\right]$ (10)

then the new variant of the Lanczos algorithm delivers the solution

$\begin{array}{ccc}\left\{\begin{array}{c}1.0000\\ -5.0000\\ 2.0000\\ 4.7661\times {10}^{16}\\ 4.0000\end{array}\right\}& and& \left[\begin{array}{ccccc}1.0000& 0.0000& -0.0023& 0.0000& -0.0007\\ 0.0000& 0.0000& 0.0000& 0.5774& 0.0000\\ 0.0000& -0.4472& 0.0000& 0.0000& -0.0005\\ 0.0004& -0.0005& 0.0000& 0.0000& 0.5000\\ -0.0016& 0.0000& -0.7071& 0.0000& 0.0000\end{array}\right]\end{array}$ (11)

Notice the “infinite” eigenvalue obtained numerically (4.7661×10 ¹⁶) and the correct eigenvector associated to it (column 4).

A more practical example consists of a 40 cm × 40 cm square composite plate whose laminate is just a single layer of T300-5208 graphite/epoxy (properties in Tab. 1 ) of 0.15 mm thickness. The single layer is oriented according to angle θ as shown in Fig. 5 . Simply supported boundary conditions are imposed along the four edges and pure shear loading is applied.

Table 1
Mechanical properties of the T300-5208 graphite/epoxy.

Property	Value
Longitudinal modulus of elasticity E1	154.0 GPa
Transverse modulus of elasticity E2	11.13 GPa
In-plane Poisson ratio ν12	0.304
In-plane shear modulus G12	6.98 GPa
Transverse shear modulus G13 = G 12	6.98 GPa
Transverse shear modulus G23	3.36 GPa

Figure 5
Single ply composite square plate.

The composite plate is modeled with a 4×4 mesh of bicubic elements ( Faria, 2000 ), and the stiffness and geometric stiffness matrices are computed. The Lanczos algorithm 3 for θ = 45^o leads to the critical loads 20.91 N/m (lowest positive eigenvalue) and −3.17 N/m (highest negative eigenvalue). Along the plate diagonal, where x = y, the pure shear load induces compressive normal stresses. Therefore, a high critical buckling load is expected (20.91 N/m). When the layer angle is θ = −45^o the critical loads become 3.17 N/m (lowest positive eigenvalue) and −20.91 N/m (highest negative eigenvalue).

Figure 6 shows how the lowest positive eigenvalue (λ_pos) and highest negative eigenvalue (λ_neg) vary with the layer angle θ. λ_pos and λ_neg have the same absolute values for θ = 0^o and θ = ±90^o. The maximum | λ| is obtained for θ = ±45^o . If the plate is required to equally support either positive or negative shear loads then the optimal strategy would be to choose θ = 0^o or θ = 90^o.

Figure 6
Critical load variation.

4 CONCLUSIONS

A deeper investigation of the Lanczos algorithm applied to the solution of the eigenvalue problem stated in Eq. (3) reveals that the principal contribution to vector x_{i +1} is from ϕ_i, where ϕ_i is the eigenvector related to the ith least eigenvalue in absolute value. Hence, vector x₂ consists mostly of ϕ₁; vector x₃ consists mostly of ϕ₂; so on so forth, where the eigenvalues are ordered such that |λ₁| < |λ ₂| < ... < |λ_n|. Thus, one expects that the procedure converges quickly to the least eigenvalues in absolute value and related eigenvectors. If a structure is subject to loadings in either direction (positive or negative), the relevant critical buckling load is the one associated to the least eigenvalue in absolute value. Therefore, it is clear that eigensolvers of practical relevance must be able to compute eigenvalues of either sign.

• http://dx.doi.org/10.1590/1679-78253932

References

• Almeida, S.F.M., Hansen, J.S., (2000). Natural frequencies of composite plates with thermal residual stresses. International Journal of Solids and Structures 36: 3517-3539.
• Bathe, K-J., (1996). Finite Element Procedures, Prentice Hall (New Jersey).
• Ericsson, T., 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.
• Faria, A.R., (2000). Buckling optimization of composite plates and cylindrical shells: uncertain loading combinations, Ph.D. Thesis, University of Toronto.
• Jones, M.T., Patrick, M.L., (1989) The use of Lanczo’s method to solve the large generalized symmetric definite N eigenvalue problem. NASA Contractor Report 181914, ICASE Report No. 89-69.
• Lanczos, C., (1950). An iteration method for the solution of the eigenvalue problem of linear differential and integral operators. Journal of Research of the Natural Bureau of Standards 45: 255-282.
• Paige, C.C., (1971). The computation of eigenvalues and eigenvectors of very large sparse matrices, Ph.D. Thesis, University of London.
• Paige, C.C., (1972). Computational variants of the Lanczos method for the eigenproblem. Journal of the Institute of Mathematics and its Applications 10: 373-381.
• Parlett, B.N., Scott, D.S., (1979). The Lanczos algorithm with selective orthogonalization. Mathematics of Computation 33: 217-238.
• Ramaswamy, S., (1980). On the effectiveness of the Lanczos method for the solution of large eigenvalue problems. Journal of Sound and Vibration 73:405-418.
• Zienkiewicz, O.C., (1991). The Finite Element Method, Vol. 2: Solid and Fluid Mechanics, Dynamics and Non-linearity. 4th ed., McGraw-Hill: London

Publication Dates

• Publication in this collection
  2018

History

• Received
  14 Apr 2017
• Reviewed
  30 May 2017
• Accepted
  22 Sept 2017