Просмотреть файл в отдельном окне: 1d679630442188a66314ee40b7a8c508.pdf

7

Eigenvalue and

Eigenvector Problems

7.1 INTRODUCTION

There is a class of physical problems which lead to a governing ordinary differential

equation containing an unknown parameter. As an example, consider the buckling

of a slender rod subjected to an axial load P shown in Figure 1. The deflected shape

y(x) is governed by the equation:1

d2y M

Py

=

=?

dx 2 EI

EI

(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 equal

to -Py) of the rod at the section x. If the rod is supported at both ends such that the

boundary conditions are:

y(x = 0) = y(x = L ) = 0

(2)

The unknown parameter appearing in Equation 1 is P which is the load axially

applied causing the rod to buckle. The problem is then to find P and the corresponding

buckled shape y(x). If the value of EI is a constant for all x, this problem can be

solved analytically. The buckling load can be shown to be P = 2EI/L2. For the

general case when EI is the function of x, numerical method has to be applied to

obtain 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 involving

the buckled shape evaluated at N selected stations between the end supports of the

rod will be of the standard form:

([A] ? ?[I]){Y} = {0},

© 2001 by CRC Press LLC

or,

[A]{Y} = ?{Y}

(3)

FIGURE 2. Another example of eigenvalue and eigenvector problem — the vibration of

three masses connected by three springs.

where the matrix [A] will depend on the distances between the stations, {Y} contains

the buckled amount of the rod at the stations, and is related to the unknown

buckling load P. Equation 3 can be interpreted as knowing a matrix [A] and trying

to 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 eigen

means proper. and {Y} in Equation 3 are called the eigenvalue and eigenvector

of [A], respectively. If N is the order of the matrix [A], there are N sets of eigenvalues

and eigenvectors. In Section 7.3, how a polynomial from which all eigenvalues of

a 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 of

these three masses is subjected to some disturbance such as the case when the mass

m3 is pulled down by a certain distance and then released, the whole system will

then be vibrating! One will be interested in knowing at what frequency will they be

oscillating up and down. To formulate the analysis, let us denote the displacements

of the masses as xi(t) for i = 1 to 3 which are functions of time t. If the elastic

constants of the three springs are denoted as ki for i = 1,2,3, it can be shown2 by

application of the Newton’s laws of motion that the governing differential equations

for the displacements are:

m1

© 2001 by CRC Press LLC

d 2 x1

+ ( k1 + k 2 )x1 ? k 2 x 2 = 0

dt 2

(4)

m2

d 2x2

? k 2 x1 + ( k 2 + k 3 )x 2 ? k 3x3 = 0

dt 2

(5)

d 2 x3

? k 3x 2 + k 3x 3 = 0

dt 2

(6)

m3

If 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 ?t

for

i = 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 = 0

2

1

2

1

1

2 2

2

2

1

2

3

2

2

3

3

and

[? k c + ( k

2 2

3

) ]

? m 3? 2 C3 sin ?t = 0

In matrix form, the above equations can be written as:

? k1 + k 2 ? m1? 2

?

?k2

?

?

0

?

?k2

k 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 have

led to no motion at all, this leaves the only choice of requiring that the coefficient

matrix be singular. In other words, its determinant must be equal to zero. The

resulting equation is a cubic polynomial and enables us to solve for three roots which

are the squared values of the frequencies (2) of the vibrating system. For each

frequency, we next need to know the associated amplitudes of the vibration. Equation

8 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 6

may not involve with an unknown parameter as in the buckling problem, but once

© 2001 by CRC Press LLC

the common frequency ? is introduced for the free vibration it becomes a standard

eigenvalue and eigenvector problem described by the matrix [A].

Hence, the question becomes how to find the eigenvalues and their associated

eigenvector of a prescribed matrix. The methods of solution are to be discussed in

Sections 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 eigenvalue

of a matrix is provided and program EigenvIt also will be presented. Prior to these

discussions, in the next section we will first concern with how the matrices connected

with the buckling and vibration problems are to be derived and demonstrate in

advance how the programs CharacEq, Bairstow, EigenVec, and EigenvIt are to

be employed for obtaining the eigenvalues and eigenvectors of these matrices.

7.2 PROGRAMS EigenODE.Stb AND EigenODE.Vib —

FOR SOLVING STABILITY AND VIBRATION PROBLEMS

In order to obtain numerical solution of the buckling load and shape of the rod

shown in Figure 1 in Section 7.1, the central-difference method introduced in

Chapter 4 can be applied to approximate the second derivative term appearing in

Equation 1 there. At a typical location along the rod, say x = xj, Equation 1 can be

approximated as:

Py j

d 2 y y j?1 ? 2 y j + y j+1

=?

=?

2

2

dx

h

(EI) j

(1)

where (EI)j is the rigidity of the rod at xj and yj ? y(x = xj) etc. This approach requires

that the rod be investigated at N stations between the two supports which are labeled

as x0 and xN + 1. These stations are equally spaced so that the increment of x (stepsize

h) is simply h = x = L/(N + 1). As a result of such arrangement, the boundary

conditions 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 = 1

and j = N, based on Equation 1 and the boundary conditions (2), the two simplified

equations 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 LLC

Also 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 jth

equation by –(EI)j/h2 for j = 1,2,…,N, Equations 3 to 5 can be further simplified into

the 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 elements

are to be calculated with the formulas:

?2( EI) i h 2 ,

?

a ij = ??( EI)i h 2 ,

?

?

??0,

for i = j

for i = j ? 1 or i = j + 1

(6)

elsewhere

As 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-between

points since the stepsize h = L/(N + 1) = 1/3. [A] is of order 2 by 2 and having

elements a11 = a22 = 18 and a12 = a21 = –9 according to Equation 7. The eigenvalues

of [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 reduced

to 100x(9.87–9.4)/9.87 = 4.76%.

PROGRAM EIGENODE.STB

Buckling problem belongs to a general class of stability problems, for which a

program called EigenODE.Stb is developed to demonstrate how different increments

or different number of stations can be adopted to continue improving the solution

of eigenvalues and eigenvectors with the aid of programs EigenVec, EigenvIt and

Bairstow. The following shows the interactive application of this program.

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Applications

When program EigenODE.Stb is run for the buckling problem, the screen will

show 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 stations

between the supported ends of the rod are stored in the vector {Y} and is equal

to the buckling load P. The resulting display is:

© 2001 by CRC Press LLC

To find the eigenvalues of the above listed matrices, program CharactEq can

be applied by interactively specifying the elements in these matrices to obtain the

respective 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 = 0

5

4

3

2

7

The eigenvalues for these equations can be found by application of the program

Bairstow. 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 eigenvalue

in magnitude found for the fourth equation is 9.64569. It indicates that if the rod is

partitioned into finer and finer increments, the numerical solution continue to

improve in predicting the first buckling load from 9, 9.3726, 9.5492, to 9.64569 and

converging to the exact value of P = = 2 = 9.8696. For further improvement, the

derivation 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.VIB

For a better understanding of the vibration problem also introduced in

Section 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

?10

20

?10

0 ?

?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 LLC

program 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 VERSION

Sample Application

To 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 to

generate the matrix [A] and then the program CharacEq is used to obtain a characteristic equation 3 – 502 – 600 — 1000 = 0. The program Bairstow enables its

roots 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 and

EigenvIt (to be introduced in Section 7.5) for inverting the matrix [C] and then

© 2001 by CRC Press LLC

iteratively finding the smallest eigenvalue in magnitude (which is related to the

lowest natural frequency of the vibration). The resulting display from using these

two programs is:

The smallest eigenvalue in magnitude of the matrix [A] is therefore equal to

1/0.50489, or, 1.9806 same as obtained by application of the programs CharacEq

and Bairstow.

MATLAB APPLICATIONS

A file EigenvIt.m for MATLAB has been developed and is listed and discussed

in 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 guessed

eigenvector V0, and tries to find the eigenvector Eigenvec and eigenvalue Lambda

iteratively until the sum of the absolute values of the differences of the components

© 2001 by CRC Press LLC

of two consecutive guessed eigenvectors is less than the specified tolerance Tol. The

number of iterations is limited by the user to be no more than NT times. The reader

should refer to the program EigenvIt for more details, here provide a simple example

of using EigenvIt.m:

The display indicates that for the specified matrix [A] of order equal to 3, the

largest 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 sum

of the absolute values of the differences of the corresponding components of the

guessed eigenvectors obtained during the seventh and eighth iterations is less than

the 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 Sample

Applications:

Notice that inv.m of MATLAB has been applied to find the inverse of [A] and

using it for EigenvIt.m to find the eigenvalue and eigenvector by iteration.

© 2001 by CRC Press LLC

MATHEMATICA APPLICATIONS

For 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.166667

In[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 above

derived 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 CHARACTERISTIC

EQUATION OF A SPECIFIED SQUARE MATRIX

The program CharacEq is designed to generate the coefficients of the characteristic

equation of an interactively specified square matrix by use of the Feddeev-Leverrier

method. Such a characteristic equation is needed in the stability, vibration, and other

so-called eigenvalue problems.3 Readers interested in these problems should also

refer to the discussions on the programs EigenODE and EigenVec. The former

discusses how the square matrix is to be generated by finite-difference approximation

of ordinary differential equation. The latter program delineates how the eigenvectors

are to be found by a modified Gaussian elimination method for each eigenvalue and

how the eigenvalues are to be solved from the characteristic equation by the program

Bairstow. Here for derivation of the characteristic equation, let us denote the specified

square 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 recursive

formulas:

pk =

1

Trace of [B]k

k

for 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, of

a square matrix is the sum of the diagonal elements. A specific, numerical example

will help further explain the details involved in applying the formulas presented

above. Consider a square matrix:

? 0

[A] = ???10

?? ?2

2

?1

4

3?

2 ??

7??

(5)

Then, [B]1 = [A] and p1 = Trace([B]1) = 0–1+7 = 6. The other p’s and [B]’s are

to 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

?4

0

6

0

3?

2 ??

7??

? ?6

??10

?

?? ?2

2

?7

4

3?

2 ??

1 ??

7 ?

?30 ?? p2 = Trace ([B]2 ) 2 = ( ?26 ? 5 + 9) 2 = ?11

9 ??

? 0

B

=

A

B

+

11

I

=

[ ]3 [ ]([ ]2 [ ]) ???10

?? ?2

?6

= ??0

??0

2

?1

4

2

?1

4

3?

2 ??

7??

? ?15

? 66

?

???42

?2

6

?4

7 ?

?30 ??

20 ??

0?

0 ?? p3 = Trace ([B]3 ) 3 = (6 + 6 + 6) 3 = 6

6 ??

and finally, the characteristic equation is:

(?1)3 (?3 ? p1?2 ? p2 ? ? p3 ) = ? ?3 + 6?2 ? 11? + 6 = 0

© 2001 by CRC Press LLC

Both QuickBASIC and FORTRAN version of the program CharacEq have

been made available for derivation of the characteristic equation based on the

Feddeev-Leverrier method. The program listings are presented below along with

some sample applications.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

The 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 LLC

Sample Application

© 2001 by CRC Press LLC

MATLAB APPLICATION

MATLAB 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 how

to specify a square matrix of order 3, how poly.m is to be called, and the resulting

display:

For the FORTRAN sample problem, we can have:

Here, we can apply plot.m and polyval of MATLAB to graphically explore the

roots 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 line

intercepts the polynomial curve, it helps indicate where the real roots are. To actually

calculate the values of all roots, real or complex, the roots.m of MATLAB can be

applied as follows:

© 2001 by CRC Press LLC

FIGURE 3.

These results complement well with those presented in Figure 3.

MATHEMATICA APPLICATIONS

For finding the characteristic equation of a given matrix, Mathematica’s function Det which derives the determinant of a specified matrix can be employed. To

do 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 LLC

Output[1] =

m = {{0,2,3), {–10,–1,2}, {–2,4,7}}

Notice that the elements in each row are separated by comma and enclosed by

a pair of braces, and rows are separated also by comma. Next, we derive the

characteristic equation of the matrix m.

Input[2]: =

Det[m — x IdentityMatrix[3]]

Output[2] =

6 – 11 X + 6 X2 — X3

Input[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 – X3

We 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 LLC

Output[6] =

Notice that the graph intercepts the x axis at x = 2, x = 3, and x = 4.

7.4 PROGRAM EIGENVEC — SOLVING EIGENVECTOR

BY GAUSSIAN ELIMINATION METHOD

The program EigenVec is designed to solve for the associated eigenvector {V} when

an eigenvalue of a given square matrix [A] is specified. Eigenvalue and eigenvector

problems are discussed in the programs CharacEq and EigenODE. Here, we

describe how the Gaussian Elimination method can be modified for finding the

eigenvector {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. This equation has nontrivial solution

only if the determinant of the coefficient matrix [A]–[I] is equal to zero. In other

words, the linear algebraic equations represented by Equation 1 are not all independent. The number of equations which are dependent on the other equations, is equal

to the multiplicity of the specified . For example, if the matrix [A] is of order N

and if the multiplicity of is M which means M characteristic roots are equal to ,

then there are M equations in Equation 1 are dependent on the other N-M equations.

© 2001 by CRC Press LLC

When Gaussian Elimination method is applied for solving {V} from Equation

1, the normalization of the last equation cannot be carried out if has a multiplicity

equal to 1 even with the pivoting provision in the program. This is because one of

the N equation is dependent on the other N–1 equations. But, it suggests that we

may assign the last component of {V} to be equal to an arbitrary constant c and

express the other components of {V} in terms of c. This concept can be extended

to the case when ? has a multiplicity of M. Since only N-M equations of (1) are

independent, there are M independent solutions of {V}. To obtain the first solution,

we assign the last component of {V} a value c1 and the other last M–1 components

of {V} equal to zero and then proceed to express the first N-M components of {V}

in terms of c1. To obtain the second solution, we assign the next to the last component

of {V} a value c2 and the other last M–1 components of {V} equal to zero and

express the first N-M components of {V} in terms of c2, and so on. The M solution

of {V} can thus be expressed in terms of ci for i = 1,2,…,M.

The program EigenVec is developed from modifying the program Gauss by

following the above-explained procedure. This program requires the user to interactively specify the order N of [A], the elements of [A], a specified value of and

its multiplicity M. The results produced by the program EigenVec are the M set of

eigenvectors {V}. Both FORTRAN and QuickBASIC versions of this programs are

listed below along with sample applications.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

© 2001 by CRC Press LLC

MATLAB APPLICATIONS

MATLAB has a file called eig.m which can be applied for finding the eigenvalues and normalized vectors of a specified square matrix. To do so, we first

interactively specify the elements of a matrix [A} and then ask for the eigenvalues

Lambda and normalized eigenvectors EigenVec by entering (such as for the sample

problem in the FORTRAN and QuickBASIC versions)

>> A = [3,0,2;0,5,0;2,03]; [EigenVec,Lambda] = eig(A)

It results in a display on the screen:

Notice that the eigenvalues are listed in the diagonal of the matrix Lambda and

the corresponding normalized eigenvectors are listed in the matrix EigenVec as

columns. To list the eigenvalues in a vector Lambda, we could enter:

>> [Lambda] = eig(A)

The resulting display is:

Lambda =

5

1

5

MATHEMATICA APPLICATIONS

Mathematica has functions Eigenvalues and Eigenvectors which can be

applied to find the eigenvalues and eigenvectors, respectively, for a specified matrix

as illustrated by the following example:

© 2001 by CRC Press LLC

In[1]: = a = {{3,0,2},{0,5,0},{2,0,3}}; MatrixForm[a]

Out[1]//MatrixForm: =

3

0

2

0

5

0

2

0

3

In[2]: = Eigenvalues[a]

Out[2]: =

{1, 5, 5}

In[3]: = Eigenvectors[a]

Out[3]: =

{–1, 0, 1}, {1, 0, 1}, {0, 1, 0}}

Notice that the computed eigenvectors are not normalized.

As another example, consider the matrix M generated in the program EigenODE

for the buckling problem when the number of stations is equal to 5. To obtain the

eigenvalues, the interactive application of Mathematica goes as:

In[4]: = Eigenvalues[M]

Out[4]: =

{134.354, 108., 72., 36., 9.64617}

Notice that the smallest eigenvalue is equal to 9.64617 which predicts the lowest

buckling load. Since the exact solution is 9.8696, this further indicates that by

continuously increasing the number of stations the smallest eigenvalue in magnitude

will eventually converge to the expected value.

PRINCIPAL STRESSES

AND

PLANES

As another example of solving the eigenvalues and eigenvectors, consider the

problem of determining the principal stresses at a point within a two-dimensional

body which is subjected to in-plane loadings. If the normal stresses (x and ?y) and

shear stresses (xy = yx), Figure 5, at that point are known, it is a common practice

to graphically determine the principal stresses and principal planes, on which the

principal stresses act by use of Mohr’s circle.4 But, here we demonstrate how the

principal stresses and principal planes can be solved as the eigenvalues and eigenvectors, respectively, of a matrix [A] constructed using the values of x, y, and xy

as follows:

© 2001 by CRC Press LLC

FIGURE 5. If the normal stresses (x and ?y) and shear stresses (xy = yx), are known, it

is a common practice to graphically determine the principal stresses and principal planes, on

which the principal stresses act by use of Mohr’s circle.

? ?x

[A] = ??

?

? xy ?

? y ??

yx

(2)

For the three-dimensional cases, the normal stresses x, y, and z, and shear

stresses xy, yz, and zx (yx = xy, zy = yz, and xz = zx) are involved, Figure 6.

Again, the Mohr’s circle method can be applied to graphically solve for the principal

stresses and the principal planes, on which they act.6 But, as an extension of Equation

2, these principal stresses and principal planes can be determined as the eigenvalues

and eigenvectors, respectively, of a matrix constructed using the values of the normal

and shear stresses as follows:

? ?x

?

[A] = ?? yx

??

? zx

? xy

?y

? zy

? xz ?

?

? yz ?

? z ??

(3)

Presented below are MATLAB solutions of two problems: (a) a two-dimensional

case of x = 50, y = –30, and xy = yx = –20, and (b) a three-dimensional case of

x = 25, y = 36, z = 49, xy = yx = –12, yz = zy = 8, and zx = xz = –9, all in N/cm2.

© 2001 by CRC Press LLC

FIGURE 6. For the three-dimensional cases, the normal stresses x, y, and z, and shear

stresses xy, yz, and zx (yx = xy, zy = yz, and xz = zx) are involved.

Notice that for Problem (a), the result indicates that maximum principal stress

equal to 54.7214 N/cm2 is on a plane having an outward normal vector whose

directional cosines are equal to –0.9732 and 0.2298. That is to say this principal

plane has an outward normal vector making an angle of max = 166.7° (cosmax =

–0.9732 and cos[90°–max] = 0.2298) measured counterclockwise from the x-axis.

The minimum principal stress is found to be equal to –34.7217 N/cm2 which is on

a plane having an outward normal vector whose directional cosines are equal to

–0.2298 and –0.9732, or at an angle equal to min = –103.3° (cosmin = 0.2298 and

© 2001 by CRC Press LLC

cos[90°– min] = –0.9732). The two principal planes are perpendicular to each other.

This can also be proven by taking the dot product of the two normalized eigenvectors:

(–0.9732i + 0.2298j)•(–0.2298i- 0.9732j) = 0.

Similar observation can be made from the results for Problem (b). The principal

stresses are equal to 16.9085, 34.6906, and 58.4009 N/cm2 and they on the planes

having outward normal vectors n1 = 0.8632i + 0.4921j + 0.1192k, n2 = 0.3277i –

0.7218j + 0.6096k, and n3 = –0.3860i + 0.4867j + 0.7837k, respectively. It is easy

to prove that these principal planes are indeed orthogonal by showing that n1•n2 =

n2•n3 = n1•n3 = 0.

QUADRATIC FORMS

AND

CANONICAL TRANSFORMATION

Another interesting application of the procedure involved in solving eigenvalues

and eigenvectors of a square matrix is the canonical transformation of quadratic

forms6 for Consider a surface described by the equation:

25x 2 + 36y 2 + 49z 2 ? 24 xy ? 18xz + 16yz = 100

(4)

The left-hand side is called a quadratic form in x, y, and z. The surface is an

ellipsoid but it is difficult to make out what are the values of its major and minor

axes, for which a transformation of the coordinate system is necessary for changing

the quadratic form into a canonical one. By canonical form, it means that only the

squared terms should remain. To find what transformation is needed, we first write

Equation 4 in matrix form as:

{V}T [A]{V} = 100

(5)

where:

?x ?

{V} = ??y??

?? z ??

and

? 25

[A] = ???12

?? ?9

?12

36

8

?9?

8 ??

49 ??

(6,7)

It can be shown6 that if we find the eigenvalues and eigenvectors of [A] which

in fact are already available in the answer to the previously discussed Problem (b),

the coordinate system x-y-z can be transformed into another set of x -y -z system

by using the so-called normalized modal matrix, [Q], formed with the normalized

eigenvectors as its columns. That is:

?x? ? ?.8623

V

=

{ ?} ??y??? = ??.4921

?? z? ?? ??.1192

© 2001 by CRC Press LLC

.3277

?.7218

.6096

?.3860 ?

.4867 ??

.7837 ??

?x ?

?y ? = Q {V}

? ? [ ]

?? z ??

(8)

Since [Q] has the property of [Q]–1 = [Q]T, the quadratic form can then be written

as {V} T [A]{V} = ([Q] –1 {V }) T [A]([Q] –1 {V }) = {V } T [Q][A][Q] T {V } =

{V }[D]{V } where [D] is a diagonal matrix having the eigenvalues of [A] (16.9085,

34.6906, and 58.4009) along its diagonal. Hence, the resulting canonical form for

the surface defined by Equation 4 can now be expressed as:

16.9085x? 2 + 34.6906y? 2 + 58.4009z? 2 = 10 2

(9)

z -axis and minor axes equal to 10/(16.9085)_ and 10/(34.6906)_ along x- and yaxes, respectively. According to Equation 8, the orientations of the x -, y -, and z axes can be determined by using the rows of [Q]. For example, a unit vector I' along

x -axis is equal to 0.9623i + 0.3277j–0.3860k where i, j, and k are unit vectors along

x-, y-, and z-axes, respectively. That means, the angles between x -axis and x-, y-,

and z-axes by be obtained easily as:

?x?x = cos?1 (.8623) = 30.42 o , ?x?y = cos?1 (.3277) = 70.87o

and

?x?z = cos?1 ( ?.3860) = 112.7o

Similarly, we can find y x = cos–1(.4921) = 60.52°, y y = cos–1(-.7218) = 136.2°,

and y z = cos –1 (.4867) = 60.88°, and z x = cos –1 (.1192) = 83.15°, z y =

cos–1(.6096) = 52.44°, and z z = cos–1(.7837) = 38.40°.

To verify the above assertions, MATLAB is again applied to obtain:

© 2001 by CRC Press LLC

The arc-cosine function acos of MATLAB is employed above and pi ( =

3.14159) also has been used for converting the results in radians into degrees.

7.5 PROGRAM EIGENVIT — ITERATIVE SOLUTION

OF THE EIGENVALUE AND EIGENVECTOR

The program EigenvIt is designed to iteratively solve for the largest eigenvalue in

magnitude max and its associated eigenvector {V} of a given square matrix [A].

Eigenvalue and eigenvector problems are discussed in the programs CharacEq,

EigenODE, and EigenVec. Since the eigenvector {V} satisfies the matrix equation:

[A]{V} = ?{V}

(1)

which indicates that if we make a successful guess of {V} then when it is multiplied

by the matrix [A] the product should be a scaled version of {V} and the scaling

vector is the eigenvalue. Of course, it is not easy to guess correctly what this vector

{V} is. But, we may devise a successive guessing scheme and hope for eventual

convergence toward the needed solution. In order to make the procedure better

organized, let us use normalized vectors, that is, to require all guessed eigenvectors

to have a length equal to unity.

The iterative scheme may be written as, for k = 0,1,2,…

[A]{V( k )} = ?( k ) {V( k +1)}

(2)

where k is the iteration counter and {V(0)} is the initial guess of the eigenvector for

[A]. The iteration is to be terminated when the differences in every components

(denoted in lower case of V) of {V(k + 1)} and {V(k)} are sufficiently small. Or,

mathematically

v(i k +1 ? v(i k ) < ?

(3)

for i = 1,2,…,N and N being the order of [A].

in Equation 3 is a predetermined

tolerance of accuracy. As can be mathematically proven.7 this iterative process leads

to the largest eigenvalue in magnitude, max, and its associated eigenvectors. If it is

the smallest eigenvalue in magnitude, min, and its associated eigenvector {V} that

are necessary to be found, the iterative procedure can also be applied but instead of

© 2001 by CRC Press LLC

[A] its inverse [A]–1 should be utilized. The iteration should use the equation, for

k = 0,1,2,…

[A]?1{V( k )} = ?( k ) {V( k +1)}

(4)

When max is found, the required smallest eigenvalue in magnitude is to be

computed as:

? min = 1 ? max

(5)

The program EigenvIt has been developed following the concept explained

above. Both QuickBASIC and FORTRAN versions are made available and listed

below along with sample applications.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

MATLAB APPLICATION

A EigenvIt.m file can be created and added to MATLAB m files for iterative

solution of the eigenvalue largest in magnitude and its associated normalized eigenvector of a given square matrix. It may be written as follows:

© 2001 by CRC Press LLC

The arguments of EigenvIt are explained in the comment statements which in

MATLAB start with a character %. As illustrations of how this function can be

applied, two examples are given below.

Notice that the first attempt allows 10 iterations but the answer has been obtained

after 8 trials whereas the second attempt allowing only two trials fails to converge

and V is printed as a blank vector.

In fact, the iteration can be carried out without a M file. To resolve the above

problem, MATLAB commands can be repeatedly entered by interactive operations

as follows:

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Notice that the command format compact makes the printout lines to be closely

spaced. The iteration converges after 7 trials. This interactive method of continuous

iteration for finding the largest eigenvalue and its associated eigenvector is easy to

follow, but the repeated entering of the statement “>> Ntry = Ntry + 1…V =

V/Lambda” in the interactive execution is a cumbersome task. To circumvent this

situation, one may enter:

>> A = [2,0,3;0,5,0;3,0,2]; V = [1;0;0]; format compact

>> for Ntry = 1:100, Ntry, V = A*V; Lambda = sqrt(V’*V), V = V/Lambda, pause, end

The pause command enables each iterated result to be viewed. To continue the

trials, simply press any key. The total number of trials is arbitrarily limited at 100;

the actual need is 7 trials as indicated by the above printed results. Viewer can

terminate the iteration by pressing the and keys simultaneously after

satisfactorily seeing the 7th, converged results of Lambda = 5 and V = [0.7071; 0;

0.7071] being displayed on screen.

To iteratively find the smallest eigenvector and its associated, normalized eigenvector of [A] according to Equations 4 and 5, we apply the iterative method to [A]–1

by entering

>> Ainv = inv(A); V = [1;0;0];

>> for Ntry = 1:100, Ntry, V = Ainv*V; Lambda = 1/sqrt(V’*V), V = V*Lambda, pause, end

© 2001 by CRC Press LLC

The resulting display is:

For saving space, the last four iterations are listed in the column on the right.

The components of the last two iterated normalized eigenvectors are numerically

equal but differ in sign, it suggests that the eigenvalue is actually equal to –1.0000

instead of 1.0000.

MATHEMATICA APPLICATIONS

The While command of Mathematica can be effectively applied here for iteration of eigenvalue and eigenvector. It has two arguments, the first is a testing

expression when the condition is true the statement(s) specified in the second argument should then be implemented. When the condition is false, the While statement

© 2001 by CRC Press LLC

should be terminated. To obtain the largest eigenvalue in magnitude, we proceed

directly with a given matrix as follows:

Input[1]: = A = {{2,0,3},{0,5,0},{3,0,2}}; V = {1,0,0}; I = 0;

Input[2]: = While[i

Eigenvalue and

Eigenvector Problems

7.1 INTRODUCTION

There is a class of physical problems which lead to a governing ordinary differential

equation containing an unknown parameter. As an example, consider the buckling

of a slender rod subjected to an axial load P shown in Figure 1. The deflected shape

y(x) is governed by the equation:1

d2y M

Py

=

=?

dx 2 EI

EI

(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 equal

to -Py) of the rod at the section x. If the rod is supported at both ends such that the

boundary conditions are:

y(x = 0) = y(x = L ) = 0

(2)

The unknown parameter appearing in Equation 1 is P which is the load axially

applied causing the rod to buckle. The problem is then to find P and the corresponding

buckled shape y(x). If the value of EI is a constant for all x, this problem can be

solved analytically. The buckling load can be shown to be P = 2EI/L2. For the

general case when EI is the function of x, numerical method has to be applied to

obtain 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 involving

the buckled shape evaluated at N selected stations between the end supports of the

rod will be of the standard form:

([A] ? ?[I]){Y} = {0},

© 2001 by CRC Press LLC

or,

[A]{Y} = ?{Y}

(3)

FIGURE 2. Another example of eigenvalue and eigenvector problem — the vibration of

three masses connected by three springs.

where the matrix [A] will depend on the distances between the stations, {Y} contains

the buckled amount of the rod at the stations, and is related to the unknown

buckling load P. Equation 3 can be interpreted as knowing a matrix [A] and trying

to 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 eigen

means proper. and {Y} in Equation 3 are called the eigenvalue and eigenvector

of [A], respectively. If N is the order of the matrix [A], there are N sets of eigenvalues

and eigenvectors. In Section 7.3, how a polynomial from which all eigenvalues of

a 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 of

these three masses is subjected to some disturbance such as the case when the mass

m3 is pulled down by a certain distance and then released, the whole system will

then be vibrating! One will be interested in knowing at what frequency will they be

oscillating up and down. To formulate the analysis, let us denote the displacements

of the masses as xi(t) for i = 1 to 3 which are functions of time t. If the elastic

constants of the three springs are denoted as ki for i = 1,2,3, it can be shown2 by

application of the Newton’s laws of motion that the governing differential equations

for the displacements are:

m1

© 2001 by CRC Press LLC

d 2 x1

+ ( k1 + k 2 )x1 ? k 2 x 2 = 0

dt 2

(4)

m2

d 2x2

? k 2 x1 + ( k 2 + k 3 )x 2 ? k 3x3 = 0

dt 2

(5)

d 2 x3

? k 3x 2 + k 3x 3 = 0

dt 2

(6)

m3

If 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 ?t

for

i = 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 = 0

2

1

2

1

1

2 2

2

2

1

2

3

2

2

3

3

and

[? k c + ( k

2 2

3

) ]

? m 3? 2 C3 sin ?t = 0

In matrix form, the above equations can be written as:

? k1 + k 2 ? m1? 2

?

?k2

?

?

0

?

?k2

k 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 have

led to no motion at all, this leaves the only choice of requiring that the coefficient

matrix be singular. In other words, its determinant must be equal to zero. The

resulting equation is a cubic polynomial and enables us to solve for three roots which

are the squared values of the frequencies (2) of the vibrating system. For each

frequency, we next need to know the associated amplitudes of the vibration. Equation

8 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 6

may not involve with an unknown parameter as in the buckling problem, but once

© 2001 by CRC Press LLC

the common frequency ? is introduced for the free vibration it becomes a standard

eigenvalue and eigenvector problem described by the matrix [A].

Hence, the question becomes how to find the eigenvalues and their associated

eigenvector of a prescribed matrix. The methods of solution are to be discussed in

Sections 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 eigenvalue

of a matrix is provided and program EigenvIt also will be presented. Prior to these

discussions, in the next section we will first concern with how the matrices connected

with the buckling and vibration problems are to be derived and demonstrate in

advance how the programs CharacEq, Bairstow, EigenVec, and EigenvIt are to

be employed for obtaining the eigenvalues and eigenvectors of these matrices.

7.2 PROGRAMS EigenODE.Stb AND EigenODE.Vib —

FOR SOLVING STABILITY AND VIBRATION PROBLEMS

In order to obtain numerical solution of the buckling load and shape of the rod

shown in Figure 1 in Section 7.1, the central-difference method introduced in

Chapter 4 can be applied to approximate the second derivative term appearing in

Equation 1 there. At a typical location along the rod, say x = xj, Equation 1 can be

approximated as:

Py j

d 2 y y j?1 ? 2 y j + y j+1

=?

=?

2

2

dx

h

(EI) j

(1)

where (EI)j is the rigidity of the rod at xj and yj ? y(x = xj) etc. This approach requires

that the rod be investigated at N stations between the two supports which are labeled

as x0 and xN + 1. These stations are equally spaced so that the increment of x (stepsize

h) is simply h = x = L/(N + 1). As a result of such arrangement, the boundary

conditions 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 = 1

and j = N, based on Equation 1 and the boundary conditions (2), the two simplified

equations 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 LLC

Also 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 jth

equation by –(EI)j/h2 for j = 1,2,…,N, Equations 3 to 5 can be further simplified into

the 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 elements

are to be calculated with the formulas:

?2( EI) i h 2 ,

?

a ij = ??( EI)i h 2 ,

?

?

??0,

for i = j

for i = j ? 1 or i = j + 1

(6)

elsewhere

As 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-between

points since the stepsize h = L/(N + 1) = 1/3. [A] is of order 2 by 2 and having

elements a11 = a22 = 18 and a12 = a21 = –9 according to Equation 7. The eigenvalues

of [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 reduced

to 100x(9.87–9.4)/9.87 = 4.76%.

PROGRAM EIGENODE.STB

Buckling problem belongs to a general class of stability problems, for which a

program called EigenODE.Stb is developed to demonstrate how different increments

or different number of stations can be adopted to continue improving the solution

of eigenvalues and eigenvectors with the aid of programs EigenVec, EigenvIt and

Bairstow. The following shows the interactive application of this program.

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Applications

When program EigenODE.Stb is run for the buckling problem, the screen will

show 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 stations

between the supported ends of the rod are stored in the vector {Y} and is equal

to the buckling load P. The resulting display is:

© 2001 by CRC Press LLC

To find the eigenvalues of the above listed matrices, program CharactEq can

be applied by interactively specifying the elements in these matrices to obtain the

respective 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 = 0

5

4

3

2

7

The eigenvalues for these equations can be found by application of the program

Bairstow. 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 eigenvalue

in magnitude found for the fourth equation is 9.64569. It indicates that if the rod is

partitioned into finer and finer increments, the numerical solution continue to

improve in predicting the first buckling load from 9, 9.3726, 9.5492, to 9.64569 and

converging to the exact value of P = = 2 = 9.8696. For further improvement, the

derivation 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.VIB

For a better understanding of the vibration problem also introduced in

Section 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

?10

20

?10

0 ?

?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 LLC

program 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 VERSION

Sample Application

To 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 to

generate the matrix [A] and then the program CharacEq is used to obtain a characteristic equation 3 – 502 – 600 — 1000 = 0. The program Bairstow enables its

roots 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 and

EigenvIt (to be introduced in Section 7.5) for inverting the matrix [C] and then

© 2001 by CRC Press LLC

iteratively finding the smallest eigenvalue in magnitude (which is related to the

lowest natural frequency of the vibration). The resulting display from using these

two programs is:

The smallest eigenvalue in magnitude of the matrix [A] is therefore equal to

1/0.50489, or, 1.9806 same as obtained by application of the programs CharacEq

and Bairstow.

MATLAB APPLICATIONS

A file EigenvIt.m for MATLAB has been developed and is listed and discussed

in 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 guessed

eigenvector V0, and tries to find the eigenvector Eigenvec and eigenvalue Lambda

iteratively until the sum of the absolute values of the differences of the components

© 2001 by CRC Press LLC

of two consecutive guessed eigenvectors is less than the specified tolerance Tol. The

number of iterations is limited by the user to be no more than NT times. The reader

should refer to the program EigenvIt for more details, here provide a simple example

of using EigenvIt.m:

The display indicates that for the specified matrix [A] of order equal to 3, the

largest 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 sum

of the absolute values of the differences of the corresponding components of the

guessed eigenvectors obtained during the seventh and eighth iterations is less than

the 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 Sample

Applications:

Notice that inv.m of MATLAB has been applied to find the inverse of [A] and

using it for EigenvIt.m to find the eigenvalue and eigenvector by iteration.

© 2001 by CRC Press LLC

MATHEMATICA APPLICATIONS

For 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.166667

In[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 above

derived 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 CHARACTERISTIC

EQUATION OF A SPECIFIED SQUARE MATRIX

The program CharacEq is designed to generate the coefficients of the characteristic

equation of an interactively specified square matrix by use of the Feddeev-Leverrier

method. Such a characteristic equation is needed in the stability, vibration, and other

so-called eigenvalue problems.3 Readers interested in these problems should also

refer to the discussions on the programs EigenODE and EigenVec. The former

discusses how the square matrix is to be generated by finite-difference approximation

of ordinary differential equation. The latter program delineates how the eigenvectors

are to be found by a modified Gaussian elimination method for each eigenvalue and

how the eigenvalues are to be solved from the characteristic equation by the program

Bairstow. Here for derivation of the characteristic equation, let us denote the specified

square 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 recursive

formulas:

pk =

1

Trace of [B]k

k

for 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, of

a square matrix is the sum of the diagonal elements. A specific, numerical example

will help further explain the details involved in applying the formulas presented

above. Consider a square matrix:

? 0

[A] = ???10

?? ?2

2

?1

4

3?

2 ??

7??

(5)

Then, [B]1 = [A] and p1 = Trace([B]1) = 0–1+7 = 6. The other p’s and [B]’s are

to 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

?4

0

6

0

3?

2 ??

7??

? ?6

??10

?

?? ?2

2

?7

4

3?

2 ??

1 ??

7 ?

?30 ?? p2 = Trace ([B]2 ) 2 = ( ?26 ? 5 + 9) 2 = ?11

9 ??

? 0

B

=

A

B

+

11

I

=

[ ]3 [ ]([ ]2 [ ]) ???10

?? ?2

?6

= ??0

??0

2

?1

4

2

?1

4

3?

2 ??

7??

? ?15

? 66

?

???42

?2

6

?4

7 ?

?30 ??

20 ??

0?

0 ?? p3 = Trace ([B]3 ) 3 = (6 + 6 + 6) 3 = 6

6 ??

and finally, the characteristic equation is:

(?1)3 (?3 ? p1?2 ? p2 ? ? p3 ) = ? ?3 + 6?2 ? 11? + 6 = 0

© 2001 by CRC Press LLC

Both QuickBASIC and FORTRAN version of the program CharacEq have

been made available for derivation of the characteristic equation based on the

Feddeev-Leverrier method. The program listings are presented below along with

some sample applications.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

The 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 LLC

Sample Application

© 2001 by CRC Press LLC

MATLAB APPLICATION

MATLAB 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 how

to specify a square matrix of order 3, how poly.m is to be called, and the resulting

display:

For the FORTRAN sample problem, we can have:

Here, we can apply plot.m and polyval of MATLAB to graphically explore the

roots 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 line

intercepts the polynomial curve, it helps indicate where the real roots are. To actually

calculate the values of all roots, real or complex, the roots.m of MATLAB can be

applied as follows:

© 2001 by CRC Press LLC

FIGURE 3.

These results complement well with those presented in Figure 3.

MATHEMATICA APPLICATIONS

For finding the characteristic equation of a given matrix, Mathematica’s function Det which derives the determinant of a specified matrix can be employed. To

do 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 LLC

Output[1] =

m = {{0,2,3), {–10,–1,2}, {–2,4,7}}

Notice that the elements in each row are separated by comma and enclosed by

a pair of braces, and rows are separated also by comma. Next, we derive the

characteristic equation of the matrix m.

Input[2]: =

Det[m — x IdentityMatrix[3]]

Output[2] =

6 – 11 X + 6 X2 — X3

Input[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 – X3

We 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 LLC

Output[6] =

Notice that the graph intercepts the x axis at x = 2, x = 3, and x = 4.

7.4 PROGRAM EIGENVEC — SOLVING EIGENVECTOR

BY GAUSSIAN ELIMINATION METHOD

The program EigenVec is designed to solve for the associated eigenvector {V} when

an eigenvalue of a given square matrix [A] is specified. Eigenvalue and eigenvector

problems are discussed in the programs CharacEq and EigenODE. Here, we

describe how the Gaussian Elimination method can be modified for finding the

eigenvector {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. This equation has nontrivial solution

only if the determinant of the coefficient matrix [A]–[I] is equal to zero. In other

words, the linear algebraic equations represented by Equation 1 are not all independent. The number of equations which are dependent on the other equations, is equal

to the multiplicity of the specified . For example, if the matrix [A] is of order N

and if the multiplicity of is M which means M characteristic roots are equal to ,

then there are M equations in Equation 1 are dependent on the other N-M equations.

© 2001 by CRC Press LLC

When Gaussian Elimination method is applied for solving {V} from Equation

1, the normalization of the last equation cannot be carried out if has a multiplicity

equal to 1 even with the pivoting provision in the program. This is because one of

the N equation is dependent on the other N–1 equations. But, it suggests that we

may assign the last component of {V} to be equal to an arbitrary constant c and

express the other components of {V} in terms of c. This concept can be extended

to the case when ? has a multiplicity of M. Since only N-M equations of (1) are

independent, there are M independent solutions of {V}. To obtain the first solution,

we assign the last component of {V} a value c1 and the other last M–1 components

of {V} equal to zero and then proceed to express the first N-M components of {V}

in terms of c1. To obtain the second solution, we assign the next to the last component

of {V} a value c2 and the other last M–1 components of {V} equal to zero and

express the first N-M components of {V} in terms of c2, and so on. The M solution

of {V} can thus be expressed in terms of ci for i = 1,2,…,M.

The program EigenVec is developed from modifying the program Gauss by

following the above-explained procedure. This program requires the user to interactively specify the order N of [A], the elements of [A], a specified value of and

its multiplicity M. The results produced by the program EigenVec are the M set of

eigenvectors {V}. Both FORTRAN and QuickBASIC versions of this programs are

listed below along with sample applications.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

© 2001 by CRC Press LLC

MATLAB APPLICATIONS

MATLAB has a file called eig.m which can be applied for finding the eigenvalues and normalized vectors of a specified square matrix. To do so, we first

interactively specify the elements of a matrix [A} and then ask for the eigenvalues

Lambda and normalized eigenvectors EigenVec by entering (such as for the sample

problem in the FORTRAN and QuickBASIC versions)

>> A = [3,0,2;0,5,0;2,03]; [EigenVec,Lambda] = eig(A)

It results in a display on the screen:

Notice that the eigenvalues are listed in the diagonal of the matrix Lambda and

the corresponding normalized eigenvectors are listed in the matrix EigenVec as

columns. To list the eigenvalues in a vector Lambda, we could enter:

>> [Lambda] = eig(A)

The resulting display is:

Lambda =

5

1

5

MATHEMATICA APPLICATIONS

Mathematica has functions Eigenvalues and Eigenvectors which can be

applied to find the eigenvalues and eigenvectors, respectively, for a specified matrix

as illustrated by the following example:

© 2001 by CRC Press LLC

In[1]: = a = {{3,0,2},{0,5,0},{2,0,3}}; MatrixForm[a]

Out[1]//MatrixForm: =

3

0

2

0

5

0

2

0

3

In[2]: = Eigenvalues[a]

Out[2]: =

{1, 5, 5}

In[3]: = Eigenvectors[a]

Out[3]: =

{–1, 0, 1}, {1, 0, 1}, {0, 1, 0}}

Notice that the computed eigenvectors are not normalized.

As another example, consider the matrix M generated in the program EigenODE

for the buckling problem when the number of stations is equal to 5. To obtain the

eigenvalues, the interactive application of Mathematica goes as:

In[4]: = Eigenvalues[M]

Out[4]: =

{134.354, 108., 72., 36., 9.64617}

Notice that the smallest eigenvalue is equal to 9.64617 which predicts the lowest

buckling load. Since the exact solution is 9.8696, this further indicates that by

continuously increasing the number of stations the smallest eigenvalue in magnitude

will eventually converge to the expected value.

PRINCIPAL STRESSES

AND

PLANES

As another example of solving the eigenvalues and eigenvectors, consider the

problem of determining the principal stresses at a point within a two-dimensional

body which is subjected to in-plane loadings. If the normal stresses (x and ?y) and

shear stresses (xy = yx), Figure 5, at that point are known, it is a common practice

to graphically determine the principal stresses and principal planes, on which the

principal stresses act by use of Mohr’s circle.4 But, here we demonstrate how the

principal stresses and principal planes can be solved as the eigenvalues and eigenvectors, respectively, of a matrix [A] constructed using the values of x, y, and xy

as follows:

© 2001 by CRC Press LLC

FIGURE 5. If the normal stresses (x and ?y) and shear stresses (xy = yx), are known, it

is a common practice to graphically determine the principal stresses and principal planes, on

which the principal stresses act by use of Mohr’s circle.

? ?x

[A] = ??

?

? xy ?

? y ??

yx

(2)

For the three-dimensional cases, the normal stresses x, y, and z, and shear

stresses xy, yz, and zx (yx = xy, zy = yz, and xz = zx) are involved, Figure 6.

Again, the Mohr’s circle method can be applied to graphically solve for the principal

stresses and the principal planes, on which they act.6 But, as an extension of Equation

2, these principal stresses and principal planes can be determined as the eigenvalues

and eigenvectors, respectively, of a matrix constructed using the values of the normal

and shear stresses as follows:

? ?x

?

[A] = ?? yx

??

? zx

? xy

?y

? zy

? xz ?

?

? yz ?

? z ??

(3)

Presented below are MATLAB solutions of two problems: (a) a two-dimensional

case of x = 50, y = –30, and xy = yx = –20, and (b) a three-dimensional case of

x = 25, y = 36, z = 49, xy = yx = –12, yz = zy = 8, and zx = xz = –9, all in N/cm2.

© 2001 by CRC Press LLC

FIGURE 6. For the three-dimensional cases, the normal stresses x, y, and z, and shear

stresses xy, yz, and zx (yx = xy, zy = yz, and xz = zx) are involved.

Notice that for Problem (a), the result indicates that maximum principal stress

equal to 54.7214 N/cm2 is on a plane having an outward normal vector whose

directional cosines are equal to –0.9732 and 0.2298. That is to say this principal

plane has an outward normal vector making an angle of max = 166.7° (cosmax =

–0.9732 and cos[90°–max] = 0.2298) measured counterclockwise from the x-axis.

The minimum principal stress is found to be equal to –34.7217 N/cm2 which is on

a plane having an outward normal vector whose directional cosines are equal to

–0.2298 and –0.9732, or at an angle equal to min = –103.3° (cosmin = 0.2298 and

© 2001 by CRC Press LLC

cos[90°– min] = –0.9732). The two principal planes are perpendicular to each other.

This can also be proven by taking the dot product of the two normalized eigenvectors:

(–0.9732i + 0.2298j)•(–0.2298i- 0.9732j) = 0.

Similar observation can be made from the results for Problem (b). The principal

stresses are equal to 16.9085, 34.6906, and 58.4009 N/cm2 and they on the planes

having outward normal vectors n1 = 0.8632i + 0.4921j + 0.1192k, n2 = 0.3277i –

0.7218j + 0.6096k, and n3 = –0.3860i + 0.4867j + 0.7837k, respectively. It is easy

to prove that these principal planes are indeed orthogonal by showing that n1•n2 =

n2•n3 = n1•n3 = 0.

QUADRATIC FORMS

AND

CANONICAL TRANSFORMATION

Another interesting application of the procedure involved in solving eigenvalues

and eigenvectors of a square matrix is the canonical transformation of quadratic

forms6 for Consider a surface described by the equation:

25x 2 + 36y 2 + 49z 2 ? 24 xy ? 18xz + 16yz = 100

(4)

The left-hand side is called a quadratic form in x, y, and z. The surface is an

ellipsoid but it is difficult to make out what are the values of its major and minor

axes, for which a transformation of the coordinate system is necessary for changing

the quadratic form into a canonical one. By canonical form, it means that only the

squared terms should remain. To find what transformation is needed, we first write

Equation 4 in matrix form as:

{V}T [A]{V} = 100

(5)

where:

?x ?

{V} = ??y??

?? z ??

and

? 25

[A] = ???12

?? ?9

?12

36

8

?9?

8 ??

49 ??

(6,7)

It can be shown6 that if we find the eigenvalues and eigenvectors of [A] which

in fact are already available in the answer to the previously discussed Problem (b),

the coordinate system x-y-z can be transformed into another set of x -y -z system

by using the so-called normalized modal matrix, [Q], formed with the normalized

eigenvectors as its columns. That is:

?x? ? ?.8623

V

=

{ ?} ??y??? = ??.4921

?? z? ?? ??.1192

© 2001 by CRC Press LLC

.3277

?.7218

.6096

?.3860 ?

.4867 ??

.7837 ??

?x ?

?y ? = Q {V}

? ? [ ]

?? z ??

(8)

Since [Q] has the property of [Q]–1 = [Q]T, the quadratic form can then be written

as {V} T [A]{V} = ([Q] –1 {V }) T [A]([Q] –1 {V }) = {V } T [Q][A][Q] T {V } =

{V }[D]{V } where [D] is a diagonal matrix having the eigenvalues of [A] (16.9085,

34.6906, and 58.4009) along its diagonal. Hence, the resulting canonical form for

the surface defined by Equation 4 can now be expressed as:

16.9085x? 2 + 34.6906y? 2 + 58.4009z? 2 = 10 2

(9)

z -axis and minor axes equal to 10/(16.9085)_ and 10/(34.6906)_ along x- and yaxes, respectively. According to Equation 8, the orientations of the x -, y -, and z axes can be determined by using the rows of [Q]. For example, a unit vector I' along

x -axis is equal to 0.9623i + 0.3277j–0.3860k where i, j, and k are unit vectors along

x-, y-, and z-axes, respectively. That means, the angles between x -axis and x-, y-,

and z-axes by be obtained easily as:

?x?x = cos?1 (.8623) = 30.42 o , ?x?y = cos?1 (.3277) = 70.87o

and

?x?z = cos?1 ( ?.3860) = 112.7o

Similarly, we can find y x = cos–1(.4921) = 60.52°, y y = cos–1(-.7218) = 136.2°,

and y z = cos –1 (.4867) = 60.88°, and z x = cos –1 (.1192) = 83.15°, z y =

cos–1(.6096) = 52.44°, and z z = cos–1(.7837) = 38.40°.

To verify the above assertions, MATLAB is again applied to obtain:

© 2001 by CRC Press LLC

The arc-cosine function acos of MATLAB is employed above and pi ( =

3.14159) also has been used for converting the results in radians into degrees.

7.5 PROGRAM EIGENVIT — ITERATIVE SOLUTION

OF THE EIGENVALUE AND EIGENVECTOR

The program EigenvIt is designed to iteratively solve for the largest eigenvalue in

magnitude max and its associated eigenvector {V} of a given square matrix [A].

Eigenvalue and eigenvector problems are discussed in the programs CharacEq,

EigenODE, and EigenVec. Since the eigenvector {V} satisfies the matrix equation:

[A]{V} = ?{V}

(1)

which indicates that if we make a successful guess of {V} then when it is multiplied

by the matrix [A] the product should be a scaled version of {V} and the scaling

vector is the eigenvalue. Of course, it is not easy to guess correctly what this vector

{V} is. But, we may devise a successive guessing scheme and hope for eventual

convergence toward the needed solution. In order to make the procedure better

organized, let us use normalized vectors, that is, to require all guessed eigenvectors

to have a length equal to unity.

The iterative scheme may be written as, for k = 0,1,2,…

[A]{V( k )} = ?( k ) {V( k +1)}

(2)

where k is the iteration counter and {V(0)} is the initial guess of the eigenvector for

[A]. The iteration is to be terminated when the differences in every components

(denoted in lower case of V) of {V(k + 1)} and {V(k)} are sufficiently small. Or,

mathematically

v(i k +1 ? v(i k ) < ?

(3)

for i = 1,2,…,N and N being the order of [A].

in Equation 3 is a predetermined

tolerance of accuracy. As can be mathematically proven.7 this iterative process leads

to the largest eigenvalue in magnitude, max, and its associated eigenvectors. If it is

the smallest eigenvalue in magnitude, min, and its associated eigenvector {V} that

are necessary to be found, the iterative procedure can also be applied but instead of

© 2001 by CRC Press LLC

[A] its inverse [A]–1 should be utilized. The iteration should use the equation, for

k = 0,1,2,…

[A]?1{V( k )} = ?( k ) {V( k +1)}

(4)

When max is found, the required smallest eigenvalue in magnitude is to be

computed as:

? min = 1 ? max

(5)

The program EigenvIt has been developed following the concept explained

above. Both QuickBASIC and FORTRAN versions are made available and listed

below along with sample applications.

QUICKBASIC VERSION

© 2001 by CRC Press LLC

Sample Application

FORTRAN VERSION

© 2001 by CRC Press LLC

Sample Application

MATLAB APPLICATION

A EigenvIt.m file can be created and added to MATLAB m files for iterative

solution of the eigenvalue largest in magnitude and its associated normalized eigenvector of a given square matrix. It may be written as follows:

© 2001 by CRC Press LLC

The arguments of EigenvIt are explained in the comment statements which in

MATLAB start with a character %. As illustrations of how this function can be

applied, two examples are given below.

Notice that the first attempt allows 10 iterations but the answer has been obtained

after 8 trials whereas the second attempt allowing only two trials fails to converge

and V is printed as a blank vector.

In fact, the iteration can be carried out without a M file. To resolve the above

problem, MATLAB commands can be repeatedly entered by interactive operations

as follows:

© 2001 by CRC Press LLC

© 2001 by CRC Press LLC

Notice that the command format compact makes the printout lines to be closely

spaced. The iteration converges after 7 trials. This interactive method of continuous

iteration for finding the largest eigenvalue and its associated eigenvector is easy to

follow, but the repeated entering of the statement “>> Ntry = Ntry + 1…V =

V/Lambda” in the interactive execution is a cumbersome task. To circumvent this

situation, one may enter:

>> A = [2,0,3;0,5,0;3,0,2]; V = [1;0;0]; format compact

>> for Ntry = 1:100, Ntry, V = A*V; Lambda = sqrt(V’*V), V = V/Lambda, pause, end

The pause command enables each iterated result to be viewed. To continue the

trials, simply press any key. The total number of trials is arbitrarily limited at 100;

the actual need is 7 trials as indicated by the above printed results. Viewer can

terminate the iteration by pressing the and keys simultaneously after

satisfactorily seeing the 7th, converged results of Lambda = 5 and V = [0.7071; 0;

0.7071] being displayed on screen.

To iteratively find the smallest eigenvector and its associated, normalized eigenvector of [A] according to Equations 4 and 5, we apply the iterative method to [A]–1

by entering

>> Ainv = inv(A); V = [1;0;0];

>> for Ntry = 1:100, Ntry, V = Ainv*V; Lambda = 1/sqrt(V’*V), V = V*Lambda, pause, end

© 2001 by CRC Press LLC

The resulting display is:

For saving space, the last four iterations are listed in the column on the right.

The components of the last two iterated normalized eigenvectors are numerically

equal but differ in sign, it suggests that the eigenvalue is actually equal to –1.0000

instead of 1.0000.

MATHEMATICA APPLICATIONS

The While command of Mathematica can be effectively applied here for iteration of eigenvalue and eigenvector. It has two arguments, the first is a testing

expression when the condition is true the statement(s) specified in the second argument should then be implemented. When the condition is false, the While statement

© 2001 by CRC Press LLC

should be terminated. To obtain the largest eigenvalue in magnitude, we proceed

directly with a given matrix as follows:

Input[1]: = A = {{2,0,3},{0,5,0},{3,0,2}}; V = {1,0,0}; I = 0;

Input[2]: = While[i