Jeremy Knight Math 639 2/17/11 Algorithm 9.3 Inverse Power Method To approximate an eigenvalue and an associated eigenvector of the LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2JVEibkYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi1RIn5GJy9GM1Enbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRj0vJSlzdHJldGNoeUdGPS8lKnN5bW1ldHJpY0dGPS8lKGxhcmdlb3BHRj0vJS5tb3ZhYmxlbGltaXRzR0Y9LyUnYWNjZW50R0Y9LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1GNjYtUSgmdGltZXM7RidGOUY7Rj5GQEZCRkRGRkZIL0ZLUSwwLjIyMjIyMjJlbUYnL0ZORlNGK0Y5 matrix A given a a nonzero vector x. Digits:=14:
<Text-field style="Heading 1" layout="Heading 1">InversePowerMethod(n,A,q,x0,TOL,N,display)</Text-field> InversePowerMethod:=proc(n,A,q,x0,TOL,N,display) #n=dimension,A=Matrix,x=vector,TOL=tolerance,N=max iterations local i,j,k,p,x,y,yp,ERR,A2,mu,mu0,mu1, muhat, temp,MVProduct,EuclidNorm,IdentityMatrix,euclid_norm,gauss: #option trace; ##Load Subroutines MVProduct:=proc(n,A,x) #subroutine for calculating y=Ax for Matrix A and vector x local i,j,p; p:=Array(1..n); for i from 1 to n do for j from 1 to n do p(i):=p(i)+A(i,j)*x(j): end do: end do: return(p): end proc: EuclidNorm:=proc(n,x) #subroutine to compute euclidean norm and identify the index of the norm #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: IdentityMatrix:=proc(n) local i,j,M; M:=Array(1..n,1..n): for i from 1 to n do M(i,i):=1: end do: return(M); end proc: gauss:=proc(aa, bb, n) local i, j, k, p, m, tmp, maxval, tmpval, s, x, nrow, a, b; a:=Array(1..n,1..n); b:=Array(1..n); x:=Array(1..n); nrow:=Array(1..n); # copy matrix into a local variable so original # matrix is not changed by the procedure for i from 1 to n do for j from 1 to n do a[i,j]:= aa[i,j]; end do; end do; # copy vector for same reason for i from 1 to n do b[i]:= bb[i]; end do; # initialize nrow for i from 1 to n do nrow[i]:= i; end do; s:=Array(1..n); # initialize scale factors for i from 1 to n do maxval:=abs(a[i,1]); for j from 2 to n do tmpval:=abs(a[i,j]); if ( tmpval > maxval ) then maxval:=tmpval; end if; end do; s[i]:=maxval; end do; for i from 1 to n-1 do # find pivot row maxval:= abs(a[nrow[i],i])/s[nrow[i]]; p:=i; for j from i+1 to n do tmpval:= abs(a[nrow[j],i])/s[nrow[j]]; if (tmpval > maxval) then p:=j; maxval:=tmpval; end if; end do; if ( maxval = 0 ) then return(0); end if; # switch rows if ( p <> i ) then tmp:=nrow[p]; nrow[p]:=nrow[i]; nrow[i]:=tmp; end if; # zero out column i for k from i+1 to n do m:=a[nrow[k],i]/a[nrow[i],i]; a[nrow[k],i]:=0; for j from i+1 to n do a[nrow[k],j]:=a[nrow[k],j]-m*a[nrow[i],j]; end do; b[nrow[k]]:=b[nrow[k]]-m*b[nrow[i]]; end do; if ( a[nrow[n],n] = 0 ) then return(0); end if; end do; # backward substitution x[n]:= b[nrow[n]]/a[nrow[n],n]; for j from n-1 to 1 by -1 do x[j]:=b[nrow[j]]; for k from n to j+1 by -1 do x[j]:=x[j]- a[nrow[j],k]*x[k]; end do; x[j]:=x[j]/a[nrow[j],j]; end do; return(eval(x)); end: ##End subroutine section k:=1: mu0:=0: mu1:=0: if display=1 then print ("k, x, mu, muhat"): end if: #find the smallest integer p with 1<=p<=n and abs(x[p])=norm(x) euclid_norm,p:=EuclidNorm(n,x0): x:=x0/x0[p]: A2:=A-q*IdentityMatrix(n): while k<=N do y:=gauss(A2,x,n): if y=0 then print("q is an eigenvalue", q): break: end if: 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); temp,p:=EuclidNorm(n,y): ERR,temp:=EuclidNorm(n,x -(y/y[p])): x:=y/y[p]: if ERR<TOL and k>=4 then mu:=1/mu+q: print("tolerance reached"); print("eigenvalue",mu); print("eigenvector",x); break: end if; if display=1 then print(k-1, x, (1.0/mu+q), (1.0/muhat+q)): end if; k:=k+1; mu0:=mu1: mu1:=mu: end do: if k=N then print("Maximum number of iterations exceeded"): end if: end proc:
A:=Array([[-4.0,14,0],[-5.0,13,0],[-1,0,2]]); x0:=Array([1.0,1.0,1.0]); q:=19.0/3; LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKi9VT14j LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKidmbThE JCIrTExMTGohIio= InversePowerMethod(3,A,q,x0,.00001,10,1); UTNrLH5+eCx+fm11LH5+bXVoYXQ2Ig== NiYiIiEtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqK1IhKVsjJCIrdyI9PT0nISIqJCIiIkkpaW5maW5pdHlHRiY= NiYiIiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqL01UXiMkIituOEM8ZyEiKiQiKz4sJEg1J0Yu NiYiIiMtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqZ2NWXiMkIitfX3IsZyEiKiQiKykqKXo0KydGLg== NiYiIiQtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqKyJbKVsjJCIrRDk8K2chIiokIisjeTQrKydGLg== NiYiIiUtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqY0RaXiMkIitJciwrZyEiKiQiKyYzKysrJ0Yu UTJ0b2xlcmFuY2V+cmVhY2hlZDYi NiRRK2VpZ2VudmFsdWU2IiQiKyw8KytnISIq NiRRLGVpZ2VudmVjdG9yNiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEdGJCIqRyxcXiM= LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= TTdSMApJNlJUQUJMRV9TQVZFLzI1MTM2NDIwNFgsJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIyoiJCIkJCEjUyEiIiQhI11GKUYpIiM5IiM4IiIhRi5GLiIiI0YmTTdSMApJNlJUQUJMRV9TQVZFLzI1MTM2NjU5NlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiIzUhIiJGJ0YnRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI0ODgwMzkwMFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKysrKys1ISIqJCIrbXNzc3MhIzUkIStZPi9lPkYsRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI1MTQxMzQwNFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKysrKys1ISIqJCIrKlJzXjooISM1JCErMC5fXUNGLEYmTTdSMApJNlJUQUJMRV9TQVZFLzI1MTQzNTY2MFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKysrKys1ISIqJCIrQEIzV3IhIzUkISsnKVFBJlwjRixGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI0ODg0ODEwMFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKysrKys1ISIqJCIrKWV6SDkoISM1JCErKCpRYCpcI0YsRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI1MTQ3MjU1NlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKyoqKioqKioqKiohIzUkIit2JHBHOShGKSQhK2xVJioqXCNGKUYmTTdSMApJNlJUQUJMRV9TQVZFLzI1MTQ5MDEyOFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKysrKys1ISIqJCIrXyRlRzkoISM1JCErKlwmKioqXCNGLEYm