Pao Y.C. Engineering Analysis: Interactive Methods and Programs with FORTRAN, QuickBASIC, MATLAB, and Mathematica

Подождите немного. Документ загружается.

Since [Q] has the property of [Q]

–1

= [Q]

, the quadratic form can then be written

as {V}

[A]{V} = ([Q]

–1

{V})

[A]([Q]

–1

{V}) = {V}

[Q][A][Q]

{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 deﬁned by Equation 4 can now be expressed as:

(9)

z-axis and minor axes equal to 10/(16.9085)

and 10/(34.6906)

along x- and y-

axes, 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:

and

Similarly, we can ﬁnd 

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:

16 9085 34 6906 58 4009 10

2222

...

′

=xyz

θθ

′

−

′

−

()

cos . . , cos . .

8623 30 42 3277 70 87

′

−

=−

()

cos . .

3860 112 7

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} satisﬁes the matrix equation:

(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,…

(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 sufﬁciently small. Or,

mathematically

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

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

AV V

[]

{}

AV V

kkk

[]

{}

() ()

()

( +

()

−<

[A] its inverse [A]

–1

should be utilized. The iteration should use the equation, for

k = 0,1,2,…

(4)

When 

max

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

computed as:

(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

AV V

kkk

[]

{}

−

() ()

()

λα

min

max

Sample Application

FORTRAN VERSION

Sample Application

MATLAB APPLICATION

A EigenvIt.m ﬁle can be created and added to MATLAB m ﬁles for iterative

solution of the eigenvalue largest in magnitude and its associated normalized eigen-

vector of a given square matrix. It may be written as follows:

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 ﬁrst 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 ﬁle. To resolve the above

problem, MATLAB commands can be repeatedly entered by interactive operations

as follows:

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 ﬁnding 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 <Ctrl> and <Break> 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 ﬁnd the smallest eigenvector and its associated, normalized eigen-

vector 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

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 iter-

ation of eigenvalue and eigenvector. It has two arguments, the ﬁrst is a testing

expression when the condition is true the statement(s) speciﬁed in the second argu-

ment should then be implemented. When the condition is false, the While statement

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<8, I = I + 1; VN = A.V; Lambda = Sqrt[VN.VN];

Print["Iteration # ",I," ","Lambda = ",N[Lambda],

" {V} =",V = N[VN/Lambda]]]

Output[2]: =

Iteration # 1 Lambda = 3.60555 {V} = {0.5547, 0, 0.83205}

Iteration # 2 Lambda = 4.90682 {V} = {0.734803, 0, 0.67828}

Iteration # 3 Lambda = 4.99616 {V} = {0.701427, 0, 0.712741}

Iteration # 4 Lambda = 4.99985 {V} = {0.708237, 0, 0.709575}

Iteration # 5 Lambda = 4.99999 {V} = {0.70688, 0, 0.707333}

Iteration # 6 Lambda = 5. {V} = {0.707152, 0, 0.707062}

Iteration # 7 Lambda = 5. {V} = {0.707098, 0, 0.707016}

Iteration # 8 Lambda = 5. {V} = {0.707109, 0, 0.707105}

Notice that the testing condition is whether the running iteration counter I is

less than 8 or not. This setup enables up to eight iterations to be conducted. The

function N in Input[2] requests the value of the variable inside the brackets to be

given in numerical form. For example, when the value of Lambda is displayed as

sqrt[4], it will be displayed as 2.00000 if the input is N[Lambda]. VA·VB is the dot

product of VA and VB. In Input[2] some sample printouts of the character strings

speciﬁed inside a pair of parentheses are also demonstrated.

For iterating the smallest eigenvalue in magnitude, we work on the inverse of

this given matrix as follows:

Input[3]: = Ainv = Inverse[A]; V = {1,0,0}; I = 0;

Input[4]: = While[i<8, I = I + 1; VN = Ainv.V; Lambda = Sqrt[VN.VN];

Print["Iteration # ",I," ","Lambda = ",N[Lambda],

" {V} = ",V = N[VN/Lambda]]]

Output[4]: =

Iteration # 1 Lambda = 0.72111 {V} = {–0.5547, 0, 0.83205}

Iteration # 2 Lambda = 0.981365 {V} = {0.734803, 0, –0.67828}

Iteration # 3 Lambda = 0.999233 {V} = {–0.701427, 0, 0.712741}

Iteration # 4 Lambda = 0.999969 {V} = {0.708237, –0, 0.709575}

Iteration # 5 Lambda = 0.999999 {V} = {–0.70688, 0, 0.707333}

Iteration # 6 Lambda = 1. {V} = {0.707152, 0, –0.707062}

Iteration # 7 Lambda = 1. {V} = {–0.707098, 0, 0.707016}

Iteration # 8 Lambda = 1. {V} = {0.707109, 0, –0.707105}