Ch-07 (Pao - Engineering Analysis)

7Eigenvalue andEigenvector Problems7.1 INTRODUCTIONThere is a class of physical problems which lead to a governing ordinary differentialequation containing an unknown parameter. As an example, consider the bucklingof a slender rod subjected to an axial load P shown in Figure 1. The deflected shapey(x) is governed by the equation:1d2y MPy==−dx 2 EIEI(1)FIGURE 1. The buckling of a slender rod subjected to an axial load P.where EI is the rigidity and M is the internal bending moment (in this case equalto -Py) of the rod at the section x.

If the rod is supported at both ends such that theboundary conditions are:y(x = 0) = y(x = L ) = 0(2)The unknown parameter appearing in Equation 1 is P which is the load axiallyapplied causing the rod to buckle. The problem is then to find P and the correspondingbuckled shape y(x).

If the value of EI is a constant for all x, this problem can besolved analytically. The buckling load can be shown to be P = 2EI/L2. For thegeneral case when EI is the function of x, numerical method has to be applied toobtain approximate solutions.In this chapter, we will apply the finite-difference approximation to solve Equation 1. As will be presented in Section 7.2, the resulting matrix equation involvingthe buckled shape evaluated at N selected stations between the end supports of therod will be of the standard form:([A] − λ[I]){Y} = {0},© 2001 by CRC Press LLCor,[A]{Y} = λ{Y}(3)FIGURE 2. Another example of eigenvalue and eigenvector problem — the vibration ofthree masses connected by three springs.where the matrix [A] will depend on the distances between the stations, {Y} containsthe buckled amount of the rod at the stations, and is related to the unknownbuckling load P.

Equation 3 can be interpreted as knowing a matrix [A] and tryingto find a proper vector {Y} when it is multiplied by [A], a scaled {Y} will result.This becomes the well-known eigenvector and eigenvalue problem because eigenmeans proper. and {Y} in Equation 3 are called the eigenvalue and eigenvectorof [A], respectively. If N is the order of the matrix [A], there are N sets of eigenvaluesand eigenvectors. In Section 7.3, how a polynomial from which all eigenvalues ofa given matrix can be found as roots will be discussed.As another example of eigenvalue and eigenvector problem, consider the vibration of three masses connected by three springs shown in Figure 2.

If any one ofthese three masses is subjected to some disturbance such as the case when the massm3 is pulled down by a certain distance and then released, the whole system willthen be vibrating! One will be interested in knowing at what frequency will they beoscillating up and down. To formulate the analysis, let us denote the displacementsof the masses as xi(t) for i = 1 to 3 which are functions of time t.

If the elasticconstants of the three springs are denoted as ki for i = 1,2,3, it can be shown2 byapplication of the Newton’s laws of motion that the governing differential equationsfor the displacements are:m1© 2001 by CRC Press LLCd 2 x1+ ( k1 + k 2 )x1 − k 2 x 2 = 0dt 2(4)m2d 2x2− k 2 x1 + ( k 2 + k 3 )x 2 − k 3x3 = 0dt 2(5)d 2 x3− k 3x 2 + k 3x 3 = 0dt 2(6)m3If we assume that the masses are vibrating sinusoidally with a common frequency but with different amplitudes Ci, their displacements can then be expressed as:x i (t ) = Ci sin ωtfori = 1, 2, 3(7)Substituting Equation 7 into Equations 4 to 6, we obtain:[(k + k − m ω )C − k c ]sin ωt = 0[− k A + (k + k − m ω )A − k A ]sin ωt = 0212112 2221232233and[− k c + ( k2 23) ]− m 3ω 2 C3 sin ωt = 0In matrix form, the above equations can be written as: k1 + k 2 − m1ω 2−k20−k2k 2 + k 3 − m 2ω 2−k3  C1 00  − k 3  C2  sin ωt = 0 0 k 3 − m 3ω 2  C3 (8)Since the amplitudes C1–3 and sinωt cannot be equal to zero which would haveled to no motion at all, this leaves the only choice of requiring that the coefficientmatrix be singular.

In other words, its determinant must be equal to zero. Theresulting equation is a cubic polynomial and enables us to solve for three roots whichare the squared values of the frequencies (2) of the vibrating system. For eachfrequency, we next need to know the associated amplitudes of the vibration.

Equation8 can be arranged into the standard form, Equation 3 by letting {Y} = [C1 C2 C3]T and:( k1 + k 2 ) m1λ = ω , [A] =  − k 2 m 202− k 2 m1(k 2 + k3 ) m2− k 3 m3− k3 m2 − k 3 m 3 0(9,10)This example shows that the governing ordinary differential Equations 4 to 6may not involve with an unknown parameter as in the buckling problem, but once© 2001 by CRC Press LLCthe common frequency ω is introduced for the free vibration it becomes a standardeigenvalue and eigenvector problem described by the matrix [A].Hence, the question becomes how to find the eigenvalues and their associatedeigenvector of a prescribed matrix. The methods of solution are to be discussed inSections 7.3 and 7.4 where the programs CharacEq and EigenVec are introduced.In Section 7.5, an iterative method for finding the eigenvector when an eigenvalueof a matrix is provided and program EigenvIt also will be presented.

Prior to thesediscussions, in the next section we will first concern with how the matrices connectedwith the buckling and vibration problems are to be derived and demonstrate inadvance how the programs CharacEq, Bairstow, EigenVec, and EigenvIt are tobe employed for obtaining the eigenvalues and eigenvectors of these matrices.7.2 PROGRAMS EigenODE.Stb AND EigenODE.Vib —FOR SOLVING STABILITY AND VIBRATION PROBLEMSIn order to obtain numerical solution of the buckling load and shape of the rodshown in Figure 1 in Section 7.1, the central-difference method introduced inChapter 4 can be applied to approximate the second derivative term appearing inEquation 1 there.

At a typical location along the rod, say x = xj, Equation 1 can beapproximated as:Py jd 2 y y j−1 − 2 y j + y j+1=˙=−22dxh(EI) j(1)where (EI)j is the rigidity of the rod at xj and yj ≡ y(x = xj) etc. This approach requiresthat the rod be investigated at N stations between the two supports which are labeledas x0 and xN + 1. These stations are equally spaced so that the increment of x (stepsizeh) is simply h = x = L/(N + 1). As a result of such arrangement, the boundaryconditions previously defined in Equation 2 now become:y 0 = y N +1 = 0(2)By writing out the equation for the first and last in-between stations, i.e., j = 1and j = N, based on Equation 1 and the boundary conditions (2), the two simplifiedequations are, respectively:h2P y +y =0−2 +(EI)1  1 2(3)h2P y N −1 + −2 +y = 0(EI)N  N(4)and© 2001 by CRC Press LLCAlso Equation 1 can be rearranged into the form, for j = 2,3,…,N–1h2P y j−1 + −2 +y + y = 0(EI) j  j j+1(5)where (EI)j is the rigidity of the beam at the jth station.

By multiplying the jthequation by –(EI)j/h2 for j = 1,2,…,N, Equations 3 to 5 can be further simplified intothe standard matrix form [A–I]{Y} = {0} where {0} is a null vector of order N, = P, {Y} = [y1 y2 … yN]T, and [A] = [aij] which for i,j = 1,2,…,N the elementsare to be calculated with the formulas:2( EI) i h 2 ,a ij = −( EI)i h 2 ,0,for i = jfor i = j − 1 or i = j + 1(6)elsewhereAs a simple numerical example, consider the case of EI = 1, L = 1, and N = 2.We are seeking only the solution of displacements, y1 and y2 at two in-betweenpoints since the stepsize h = L/(N + 1) = 1/3.

[A] is of order 2 by 2 and havingelements a11 = a22 = 18 and a12 = a21 = –9 according to Equation 7. The eigenvaluesof [AλI] can be easily obtained to be λ1 = 9 and λ2 = 27. The exact solution P =π2EI/L2 in this case is P = λ = π2 = 9.87, which indicates that λ1 is off about 10%from the exact value. If N is increased to 3, h = 0.25 and the eigenvalues are λ1 =9.4, λ2 = 32, and λ3 = 54.6. The error in estimating the first buckling load is reducedto 100x(9.87–9.4)/9.87 = 4.76%.PROGRAM EIGENODE.STBBuckling problem belongs to a general class of stability problems, for which aprogram called EigenODE.Stb is developed to demonstrate how different incrementsor different number of stations can be adopted to continue improving the solutionof eigenvalues and eigenvectors with the aid of programs EigenVec, EigenvIt andBairstow.

The following shows the interactive application of this program.FORTRAN VERSION© 2001 by CRC Press LLCSample ApplicationsWhen program EigenODE.Stb is run for the buckling problem, the screen willshow the coefficient [C] in the standard eigenvalue problem of the form [A–I]{Y} ={0} where the values of the buckled shape y(x) computed in 2, 3, 4, and 5 stationsbetween the supported ends of the rod are stored in the vector {Y} and is equalto the buckling load P. The resulting display is:© 2001 by CRC Press LLCTo find the eigenvalues of the above listed matrices, program CharactEq canbe applied by interactively specifying the elements in these matrices to obtain therespective characteristic equations asλ2 – 36λ + 243 = 0,λ3 – 96λ2 – 2560λ – 16384 = 0,λ4 – 200λ3 + 13125λ2 – 312500λ + 1953125 = 0,andλ – 360λ + 46656 – 2612736λ + 5.878656x10 λ – 3.62791x108 = 054327The eigenvalues for these equations can be found by application of the programBairstow.

The sets of eigenvalues for the first three equations are (9 and 27), (9.3726,32, and 54.627), and (9.5492, 34.549, 65.451, and 90.451). The smallest eigenvaluein magnitude found for the fourth equation is 9.64569. It indicates that if the rod ispartitioned into finer and finer increments, the numerical solution continue toimprove in predicting the first buckling load from 9, 9.3726, 9.5492, to 9.64569 andconverging to the exact value of P = = 2 = 9.8696. For further improvement, thederivation of the characteristic equations of order 6 and higher is given as a homework problem for the reader to practice application of the programs EigenODE.Stb,CharacEq, and Bairstow.PROGRAM EIGENODE.VIBFor a better understanding of the vibration problem also introduced inSection 7.1, let us assign values for the spring constants and masses to be k1 = k2 =k3 = 10 lb/ft and m1 = m2 = m3 = 1 lb-sec2/ft.

The matrix [A] becomes: 20[A] = −10 0−1020−100 −10 10 Indeed, [A] is singular. The determinant of [A–λI] gives the characteristic equation 3–502 + 600–1000 = 0. The roots can be obtained by application of the© 2001 by CRC Press LLCprogram Bairstow to be = 2 = 1.98, 15.5, and 32.5. For = 1.98 the frequency is equal to 1.407 radians/second, the amplitude ratios are C2/C1 = 1.80 and C3/C1 =2.25. The program EigenVib has been developed for generation of the matrix [A]when the values of the masses, m’s, and the spring constants, k’s, are provided.FORTRAN VERSIONSample ApplicationTo demonstrate application of the program EigenODE.Vib, the numerical example for the vibration of three masses shown in Figure 1 in Section 7.1 is run togenerate the matrix [A] and then the program CharacEq is used to obtain a characteristic equation 3 – 502 – 600 — 1000 = 0.