Jeremy Knight Math 639 2/17/11 Algorithm 9.1: Power Method
<Text-field style="Heading 1" layout="Heading 1">PowerMethod(n,A,x,TOL,N,display)</Text-field> PowerMethod:=proc(n,A,x0,TOL,N,display) #n=dimension, A=Array, x=vector, TOL=tolerance, N=max iteraions #This version has been modified to use Aitkens delta^2 procedure local i,j,k,p,x,xp,y,mu,yp,EuclidNorm,ERR,temp,mu0,mu1,muhat: #option trace; k:=1: mu0:=0: mu1:=0: x:=x0: if display=1 then print("iterations, eigenvalue, eigenvector, muhat"); end if; #subroutine for finding the Euclidean Norm and returning the norm and the index of the norm. EuclidNorm:=proc(n,x) #n=vector size, x= 1 by n array local euclid_norm,j,p: euclid_norm:=abs(x[1]): p:=1; for j from 2 to n do if abs(x[j])>euclid_norm then euclid_norm:=abs(x[j]): p:=j: end if: end do: return(euclid_norm,p); end proc: #find the smallest integer p with 1<=p<=n and abs(x[p])=norm(X) xp,p:=EuclidNorm(n,x); #normalize x x:=x/x[p]: while k<=N do #Calculate next eigenvector approximate y:=Array(1..n): for i from 1 to n do for j from 1 to n do y[i]:=y[i]+A[i,j]*x[j]: end do: end do: mu:=y[p]: #use Aitken's delta^2 procedure to speed convergence muhat:=mu0-((mu1-mu0)^2/(mu-2*mu1+mu0)): #find the smallest integer p with 1<=p<=n and abs(y[p])=norm(y) yp,p:=EuclidNorm(n,y): if y[p]=0 then print("Eigenvector",evalf(x)); print("A has the eigenvalue 0, select new vector x and restart"); break; end if: #Modify with abs(x)-abs(y/y[p] to find an accurate error measurement when signs alternate in eigenvector. ERR,temp:=EuclidNorm(n,abs(x)-abs(y/y[p])): x:=y/y[p]; #check for tolerance reached if ERR<TOL and k>=4 then print("The procedure was successful"); print("eigenvalue approximate (muhat)=",evalf(muhat)); print("eigenvector approximate=",evalf(x)); break; end if; if display=1 then print(k, evalf(mu), evalf(x), evalf(muhat)); end if; k:=k+1: mu0:=mu1: mu1:=mu: end do; if k=N then print("The maximum number of iterations exceeded"); end if: if ERR>TOL then print("tolerance not reached"); print("eigenvalue approximate=",evalf(muhat)); print("eigenvector approximate=",evalf(x)); end if; end proc: