alg093.txt

Numerical Anaysis 8th Edition Burden and Faires (Maple Source)

> restart;
> # INVERSE POWER METHOD ALGORITHM 9.3
> #
> # To approximate an eigenvalue and an associated eigenvector of the
> # n by n matrix A given a nonzero vector x:
> #
> # INPUT:   Dimension n; matrix A; vector x; tolerance TOL;
> #          maximum number of iterations N.
> #
> # OUTPUT:  Approximate eigenvalue MU; approximate eigenvector x
> #          or a message that the maximum number of iterations was
> #          exceeded.
> alg093 := proc() local MULTIP, SOLVE, OK, AA, NAME, INP, N, I, J, A, X, TOL, NN, FLAG, OUP, Q, S, K, LP, AMAX, B, YMU, ERR, T;
> MULTIP := proc(N,OK,NROW,Q,A) local I, M, IMAX, J, IP, L1, L2, JJ, I1, J1, K;
> # Procedure MULTIP determines the row ordering and multipliers for the
> # matrix (A-Q*I)
> for I from 1 to N do
> NROW[I-1] := I;
> od; OK := TRUE;
> I := 1;
> M := N - 1;
> while I <= M and OK = TRUE do
> IMAX := I;
> J := I+1;
> for IP from J to N do
> L1 := NROW[IMAX-1];
> L2 := NROW[IP-1];
> if abs(A[L2-1,I-1]) > abs(A[L1-1,I-1]) then
> IMAX := IP;
> fi;
> od;
> if abs(A[NROW[IMAX-1]-1,I-1]) <= 1.0e-20 then
> OK := FALSE;
> printf(`'A - Q * I is singular, Q = %.8e is an eigenvalue\n`, Q);
> else
> JJ := NROW[I-1];
> NROW[I-1] := NROW[IMAX-1];
> NROW[IMAX-1] := JJ;
> I1 := NROW[I-1];
> for JJ from J to N do
> J1 := NROW[JJ-1];
> A[J1-1,I-1] := A[J1-1,I-1] / A[I1-1,I-1];
> for K from J to N do
> A[J1-1,K-1] := A[J1-1,K-1] - A[J1-1,I-1] * A[I1-1,K-1];
> od;
> od;
> fi;
> I := I+1;
> od;
> if abs(A[NROW[N-1]-1,N-1]) <= 1.0e-20 then
> OK := FALSE;
> printf(`A - Q * I is singular, Q = %.8e is an eigenvalue\n`, Q);
> fi;
> end;
> SOLVE := proc(N,B,A,Y,NROW) local M, I, J, I1, JJ, J1, N1, L, K, N2, KK;
> # Procedure SOLVE solves the linear system (A-Q*I)*Y=X given a new
> # vector X and returns the solution in Y
> M := N - 1;
> for I from 1 to M do
> J := I+1;
> I1 := NROW[I-1];
> for JJ from J to N do;
> J1 := NROW[JJ-1];
> B[J1-1] := B[J1-1] - A[J1-1,I-1] * B[I1-1];
> od;
> od;
> N1 := NROW[N-1];
> Y[N-1] := B[N1-1] / A[N1-1,N-1];
> L := N - 1;
> for K from 1 to L do J := L - K + 1;
> JJ := J + 1;
> N2 := NROW[J-1];
> Y[J-1] := B[N2-1];
> for KK from JJ to N do
> Y[J-1] := Y[J-1] - A[N2-1,KK-1] * Y[KK-1];
> od;
> Y[J-1] := Y[J-1] / A[N2-1,J-1];
> od;
> end;
> printf(`This is the Inverse Power Method.\n`);
> OK := FALSE;
> printf(`The array will be input from a text file in the order:\n`);
> printf(`A(1,1), A(1,2), ..., A(1,n), A(2,1), A(2,2), ..., A(2,n),\n`);
> printf(`..., A(n,1), A(n,2), ..., A(n,n)\n\n`);
> printf(`Place as many entries as desired on each line, but separate `);
> printf(`entries with\n`);
> printf(`at least one blank.\n`);
> printf(`The initial approximation should follow in same format.\n\n\n`);
> printf(`Has the input file been created? - enter Y or N.\n`);
> AA := scanf(`%c`)[1];
> if AA = "Y" or AA = "y" then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example:   A:\\\DATA.DTA\n`);
> NAME := scanf(`%s`)[1];
> INP := fopen(NAME,READ,TEXT);
> OK := FALSE;
> while OK = FALSE do
> printf(`Input the dimension n.\n`);
> N := scanf(`%d`)[1];
> if N > 0 then
> for I from 1 to N do
> for J from 1 to N do
> A[I-1,J-1] := fscanf(INP, `%f`)[1];
> od;
> od;
> for I from 1 to N do
> X[I-1] := fscanf(INP, `%f`)[1];
> od;
> fclose(INP);
> while OK = FALSE do
> printf(`Input the tolerance.\n`);
> TOL := scanf(`%f`)[1];
> if TOL > 0 then
> OK := TRUE;
> else
> printf(`Tolerance must be positive number.\n`);
> fi;
> od;
> OK := FALSE;
> while OK = FALSE do
> printf(`Input maximum number of iterations `);
> printf(`- integer.\n`);
> # Use NN in place of N.
> NN := scanf(`%d`)[1];
> if NN > 0 then
> OK := TRUE;
> else
> printf(`Number must be positive integer.\n`);
> fi;
> od;
> else
> printf(`The dimension must be a positive integer.\n`);
> fi;
> od;
> else 
> printf(`The program will end so the input file can be created.\n`);
> fi;
> if OK = TRUE then
> printf(`Choice of output method:\n`);
> printf(`1. Output to screen\n`);
> printf(`2. Output to text file\n`);
> printf(`Please enter 1 or 2.\n`);
> FLAG := scanf(`%d`)[1];
> if FLAG = 2 then
> printf(`Input the file name in the form - drive:\\name.ext\n`);
> printf(`for example   A:\\OUTPUT.DTA\n`);
> NAME := scanf(`%s`)[1];
> OUP := fopen(NAME,WRITE,TEXT);
> else
> OUP := default;
> fi;
> fprintf(OUP, `INVERSE POWER METHOD\n\n`);
> # Step 1
> # Q could be input instead of computed.
> Q := 0;
> S := 0;
> for I from 1 to N do
> S := S + X[I-1] * X[I-1];
> for J from 1 to N do
> Q := Q + A[I-1,J-1] * X[I-1] * X[J-1];
> od;
> od;
> Q := Q / S;
> printf(`Q is %.8e\n`, Q);
> printf(`Input new Q? Enter Y or N.\n`);
> AA := scanf(`\n%c`)[1];
> if AA = "Y" or  AA = "y" then
> printf(`input new Q\n`);
> Q := scanf(`%f`)[1];
> fi;
> fprintf(OUP, `Iteration  Eigenvalue  Eigenvector\n`);
> # Step 2
> K := 1;
> for I from 1 to N do
> A[I-1,I-1] := A[I-1,I-1] - Q;
> od;
> # Call subroutine to compute multipliers M(I,J) and upper triangular
> # matrix for the matrix (A-Q*I)
> MULTIP(N, OK, NROW, Q, A);
> if OK = TRUE then
> # Step 3
> LP := 1;
> for I from 2 to N do
> if abs(X[I-1]) > abs(X[LP-1]) then
> LP := I;
> fi;
> od;
> # Step 4
> AMAX := X[LP-1];
> for I from 1 to N do
> X[I-1] := X[I-1] / (AMAX);
> od;
> # Step 5
> while K <= NN and OK = TRUE do 
> # Steps 6 and 7
> for I from 1 to N do
> B[I-1] := X[I-1];
> od;
> # Subroutine SOLVE returns the solution of (A-Q*I)*Y=b in Y
> SOLVE(N, B, A, Y, NROW);
> # Step 8
> YMU := Y[LP-1];
> # Steps 9 and 10
> LP := 1;
> for I from 2 to N do
> if abs(Y[I-1]) > abs(Y[LP-1]) then
> LP := I;
> fi;
> od;
> AMAX := Y[LP-1];
> ERR := 0;
> for I from 1 to N do;
> T := Y[I-1] / AMAX;
> if abs(X[I-1] - T) > ERR then
> ERR := abs(X[I-1] - T);
> fi;
> X[I-1] := T;
> od;
> YMU := 1 / YMU + Q;
> # Step 11 
> fprintf(OUP, `%3d %12.8f\n`, K, YMU);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I-1]);
> od;
> fprintf(OUP, `\n`);
> if ERR < TOL then 
> OK := FALSE;
> fprintf(OUP, `Eigenvalue = %12.8f`, YMU);
> fprintf(OUP, ` to tolerance = %.10e\n`, TOL);
> fprintf(OUP, `obtained on iteration number = %d\n\n`, K);
> fprintf(OUP, `Unit eigenvector is :\n`);
> for I from 1 to N do
> fprintf(OUP, ` %11.8f`, X[I-1]);
> od;
> fprintf(OUP, `\n`);
> else
> # Step 12
> K := K+1;
> fi;
> od;
> if K > NN then
> fprintf(OUP,  `No convergence in %d iterations\n`,NN);
> fi;
> fi;
> if OUP <> default then
> fclose(OUP):
> printf(`Output file %s created successfully`,NAME);
> fi;
> fi;
> RETURN(0);
> end;
> alg093();