Ch-07 (523180), страница 2
Текст из файла (страница 2)
The program Bairstow enables itsroots to be found as equal to 1.9806, 15.550, and 37.470.It is also of interest to show an application of the programs MatxInvD andEigenvIt (to be introduced in Section 7.5) for inverting the matrix [C] and then© 2001 by CRC Press LLCiteratively finding the smallest eigenvalue in magnitude (which is related to thelowest natural frequency of the vibration). The resulting display from using thesetwo programs is:The smallest eigenvalue in magnitude of the matrix [A] is therefore equal to1/0.50489, or, 1.9806 same as obtained by application of the programs CharacEqand Bairstow.MATLAB APPLICATIONSA file EigenvIt.m for MATLAB has been developed and is listed and discussedin the program EigenvIt.
This function is in the form of [EigenVec,Lambda] =EigenvIt(A,N,V0,NT,Tol). It accepts a matrix [A] of order N, an initial guessedeigenvector V0, and tries to find the eigenvector Eigenvec and eigenvalue Lambdaiteratively until the sum of the absolute values of the differences of the components© 2001 by CRC Press LLCof two consecutive guessed eigenvectors is less than the specified tolerance Tol. Thenumber of iterations is limited by the user to be no more than NT times. The readershould refer to the program EigenvIt for more details, here provide a simple exampleof using EigenvIt.m:The display indicates that for the specified matrix [A] of order equal to 3, thelargest eigenvalue in magnitude is equal to 5.0000 and its associated eigenvector is[0.7071 0 0.7071]T after 8 iterative steps. The iteration is terminated when the sumof the absolute values of the differences of the corresponding components of theguessed eigenvectors obtained during the seventh and eighth iterations is less thanthe specified tolerance 0.0001.To find the smallest eigenvalue and its associated eigenvector by iteration,EigenvIt.m also can be applied effectively.
Let us use the example in SampleApplications:Notice that inv.m of MATLAB has been applied to find the inverse of [A] andusing it for EigenvIt.m to find the eigenvalue and eigenvector by iteration.© 2001 by CRC Press LLCMATHEMATICA APPLICATIONSFor the buckling problem, Mathematica can be applied as follows:In[1]: = Ns = 5; EI = 1.; L = 1.; H = L/(Ns + 1);In[2]: = (Print[“Number of Station = “, Ns, “ EI = “, EI, “ Length = “, L,“Delta L = “, H)Out[2]: =Number of Station = 5 EI = 1.
Length = 1. Delta L = 0.166667In[3]: = (Do[Do[If[i == j, M[[i,j]] = 2.*EI/H^2,If[i == (j + 1)||i = = (j–1), M[[i,j]] = EI/H^2, 0]],{i,Ns}],{j,Ns}]); MatrixForm[M]Out[3]//MatrixForm: =72.36.0.0.0.36.72.36.0.0.0.36.72.36.0.0.0.36.72.36.0.0.0.36.72.In the next section, we will show how the characteristic equation for the abovederived matrix [M] can be determined by application of Mathematica and subsequently how the eigenvalues and eigenvectors are to be obtained.7.3 PROGRAM CHARACEQ — DERIVATION OF CHARACTERISTICEQUATION OF A SPECIFIED SQUARE MATRIXThe program CharacEq is designed to generate the coefficients of the characteristicequation of an interactively specified square matrix by use of the Feddeev-Leverriermethod.
Such a characteristic equation is needed in the stability, vibration, and otherso-called eigenvalue problems.3 Readers interested in these problems should alsorefer to the discussions on the programs EigenODE and EigenVec. The formerdiscusses how the square matrix is to be generated by finite-difference approximationof ordinary differential equation.
The latter program delineates how the eigenvectorsare to be found by a modified Gaussian elimination method for each eigenvalue andhow the eigenvalues are to be solved from the characteristic equation by the programBairstow. Here for derivation of the characteristic equation, let us denote the specifiedsquare matrix be [A] and its elements be ai,j for i,j = 1,2,…,n with n being the order of[A]. The Feddeev-Leverrier method first express the characteristic equation of [A] as:(−1)n (λn − p1λn −1 − p2 λn −2 − … − p n −1λ − p n ) = 0© 2001 by CRC Press LLC(1)where the coefficients p1 through pn are to be determined by the following recursiveformulas:pk =1Trace of [B]kkfor k = 1, 2,…, n(2)and[B]1 = [A]and[B]k = [A]([B]k −1 − p k −1[I])(3,4)Equation 4 is to be applied for j = 2,3,…,n.
Trace, appearing in Equation 2, ofa square matrix is the sum of the diagonal elements. A specific, numerical examplewill help further explain the details involved in applying the formulas presentedabove. Consider a square matrix: 0[A] = −10 −22−1432 7(5)Then, [B]1 = [A] and p1 = Trace([B]1) = 0–1+7 = 6. The other p’s and [B]’s areto be calculated according to Equations 2 and 4, and finally the characteristic equation is to be expressed according to Equation 1 as: 0[B]2 = [A]([B]1 − 6[I]) = −10 −2−26= 66−42−2−5−406032 7 −6−10 −22−7432 1 7 −30 p2 = Trace ([B]2 ) 2 = ( −26 − 5 + 9) 2 = −119 0B=AB+11I=[ ]3 [ ]([ ]2 [ ]) −10 −26= 002−142−1432 7 −15 66−42−26−47 −30 20 00 p3 = Trace ([B]3 ) 3 = (6 + 6 + 6) 3 = 66 and finally, the characteristic equation is:(−1)3 (λ3 − p1λ2 − p2 λ − p3 ) = − λ3 + 6λ2 − 11λ + 6 = 0© 2001 by CRC Press LLCBoth QuickBASIC and FORTRAN version of the program CharacEq havebeen made available for derivation of the characteristic equation based on theFeddeev-Leverrier method.
The program listings are presented below along withsome sample applications.QUICKBASIC VERSION© 2001 by CRC Press LLCSample ApplicationThe display screen will show the following questions-and-answers and the computed results when the matrix [A] given in (5) is interactively entered as the matrix,for which its characteristic equation is to be obtained:FORTRAN VERSION© 2001 by CRC Press LLCSample Application© 2001 by CRC Press LLCMATLAB APPLICATIONMATLAB has a file called poly.m which can be applied to obtain the characteristic equation of a specified square matrix. The following is an example of howto specify a square matrix of order 3, how poly.m is to be called, and the resultingdisplay:For the FORTRAN sample problem, we can have:Here, we can apply plot.m and polyval of MATLAB to graphically explore theroots of this obtained polynomial P(x) = x3–9x2 + 26x–24 = 0 by interactive entering:>> x = [1:0.05:5]; y = polyval(p,x); plot(x,y), hold>> XL = [1 5]; YL = [0 0]; plot(XL,YL)The resulting curve is shown in Figure 3.
Notice that the added horizontal lineintercepts the polynomial curve, it helps indicate where the real roots are. To actuallycalculate the values of all roots, real or complex, the roots.m of MATLAB can beapplied as follows:© 2001 by CRC Press LLCFIGURE 3.These results complement well with those presented in Figure 3.MATHEMATICA APPLICATIONSFor finding the characteristic equation of a given matrix, Mathematica’s function Det which derives the determinant of a specified matrix can be employed. Todo so, the matrix should be entered first and then Det is to be called next.Input[1]: =m = {{0,2,3), {–10,–1,2}, {–2,4,7}}© 2001 by CRC Press LLCOutput[1] =m = {{0,2,3), {–10,–1,2}, {–2,4,7}}Notice that the elements in each row are separated by comma and enclosed bya pair of braces, and rows are separated also by comma. Next, we derive thecharacteristic equation of the matrix m.Input[2]: =Det[m — x IdentityMatrix[3]]Output[2] =6 – 11 X + 6 X2 — X3Input[3]: =m = {{1,2,3), {–10,0,2}, {–2,4,8}}Output[3] =m = {{1,2,3), {–10,0,2}, {–2,4,8}}Input[4]: =Det[m — x IdentityMatrix[3]]Output[4] =24 – 26 X + 9 X2 – X3We may proceed to solve the characteristic roots as follows:Input[5]: =NSolve[24–26x + 9x^2x^3 = = 0,x]Output[5] ={{x -> 2.}, {x -> –3.}, {x -> 4.}}Again, the polynomial can be plotted with:Input[6]: =Plot[x^3–9x^2 + 26x–24, {x,1,5},Frame->True}, AspectRatio->1]© 2001 by CRC Press LLCOutput[6] =Notice that the graph intercepts the x axis at x = 2, x = 3, and x = 4.7.4 PROGRAM EIGENVEC — SOLVING EIGENVECTORBY GAUSSIAN ELIMINATION METHODThe program EigenVec is designed to solve for the associated eigenvector {V} whenan eigenvalue of a given square matrix [A] is specified.
Eigenvalue and eigenvectorproblems are discussed in the programs CharacEq and EigenODE. Here, wedescribe how the Gaussian Elimination method can be modified for finding theeigenvector {V}. Since the eigenvector {V} satisfies the matrix equation:([A] − λ[I]){V} = {0}(1)where [I] is the identity matrix of same order as [A]. Equation 1 is called homogeneous since the right-hand side is a null vector.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
