Jeremy Knight Math 639 2/17/11 Algorithm 9.2
<Text-field style="Heading 1" layout="Heading 1">SymmetricPowerMethod(n,A,x,TOL,N,display)</Text-field> SymmetricPowerMethod:=proc(n,A,x0,TOL,N,display) #n=dimension, A=matrix, x=initial vector, TOL=tolerance, N=Maximum number of iterations local i,j,L2Norm,MVProduct,k,x,y,u,temp,ERR,mu,mu0,mu1,muhat: #option trace; ##subroutines #subroutine for finding the L2 Norm L2Norm:=proc(n,x) local j,total; total:=0: for j from 1 to n do total:=total+x[j]^2: end do: return (evalf(total^(1/2))); end proc: #subroutine for Matrix Vector product y=Ax MVProduct:=proc(n,A,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: ##end subroutines if display=1 then print("k, x, mu^m, muhat^m, "): end if: k:=1: mu0:=0: mu1:=0; x:=x0/L2Norm(n,x0): while k<=N do y:=MVProduct(n,A,x): mu:=0: for j from 1 to n do mu:=mu+x[j]*y[j]: end do; if L2Norm(n,y)=0 then print("Eigenvector=", x); print("A has eigenvalue 0, select new vector x and restart."); break; end if; if display=1 then print(k, x, mu, muhat); end if; #use Aitken's delta^2 procedure to speed convergence muhat:=mu0-((mu1-mu0)^2/(mu-2*mu1+mu0)): temp:=L2Norm(n,y): ERR:=L2Norm(n,(x-y/temp)): x:=y/temp; if ERR<TOL then print("eigenvalue", mu); print("eigenvector", x); print("tolerance reached"); break: end if; k:=k+1: mu0:=mu1: mu1:=mu: end do; if k=N then print("Maximum number of iterations exceeded"); break; end if: end proc:
A:=Array([[4,-1,1],[-1,3,-2],[1,-2,3]]); x0:=Array([1,0,0]); LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKiNcY0JE LUkmQXJyYXlHJSpwcm90ZWN0ZWRHNiMvSSQlaWRHNiIiKktrYmkj SymmetricPowerMethod(3,A,x0,.00001,20,1); UTprLH5+eCx+fm11Xm0sfn5tdWhhdF5tLH5+NiI= NiYiIiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqM2p2aSMkIiIlIiIhSSZtdWhhdEdGKg== NiYiIiMtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqJSlwLmsjJCIrKioqKioqKipcISIqJCIiIUYw NiYiIiQtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqJUcsP0QkIitxbW1tYyEiKiQiK0tMTExgRi4= NiYiIiUtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqJyo+bmkjJCIrMzQ0NGYhIiokIitVKysrcUYu NiYiIiYtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqTyl6PEQkIisoPVduKGYhIiokIitSIT53LydGLg== NiYiIictSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqK3RoaSMkIitYPzolKmYhIiokIitaRCRIKydGLg== NiYiIigtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqPz5lXiMkIityZWApKmYhIiokIis/Sj0rZ0Yu NiYiIiktSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqc3FYaiMkIitSUWoqKmYhIiokIitZOSwrZ0Yu NiYiIiotSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqJT0xS0UkIitfJTMqKipmISIqJCIrRzIrK2dGLg== NiYiIzUtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqa11vaiMkIis7cigqKipmISIqJCIrUysrK2dGLg== NiYiIzYtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqZ1FjaiMkIit2VSoqKipmISIqJCIrOysrK2dGLg== NiYiIzctSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqU0s0XyMkIitwJikqKioqZiEiKiQiKyYqKioqKioqZkYu NiYiIzgtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqZzQqPUQkIitYJyoqKioqZiEiKiQiKy0rKytnRi4= NiYiIzktSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqK0AmUUUkIis2KioqKioqZiEiKiQiKywrKytnRi4= NiYiIzotSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqK1tuaiMkIit1KioqKioqZiEiKiQiKykqKioqKioqZkYu NiYiIzstSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqI1wjPmojJCIrJioqKioqKipmISIqJCIrKCoqKioqKipmRi4= NiYiIzwtSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqT1hKaiMkIissKysrZyEiKkYs NiYiIz0tSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEc2IiIqN1MmPkQkIissKysrZyEiKiQiKy4rKytnRi4= NiRRK2VpZ2VudmFsdWU2IiQiKywrKytnISIq NiRRLGVpZ2VudmVjdG9yNiItSSZBcnJheUclKnByb3RlY3RlZEc2Iy9JJCVpZEdGJCIqS11kaiM= UTJ0b2xlcmFuY2V+cmVhY2hlZDYi LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic= TTdSMApJNlJUQUJMRV9TQVZFLzI1MjM1NjQ5MlgsJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIyoiJCIkIiIlISIiIiIiRigiIiQhIiNGKUYrRipGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2MjU1NjQzMlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCIiIiIiIUYoRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2Mjc1NjMwOFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiIiIiIiEkRilGKUYqRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2NDAzNjk4NFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKzsvNEclKiEjNSQhKy9FLWRCRikkIisvRS1kQkYpRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI1MjAwMTI4NFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKzVlJ1w7KSEjNSQhKzJIWyMzJUYpJCIrMEhbIzMlRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2MjY3MTk5NlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK1gwcDFyISM1JCErIlEkb3VcRikkIisiUSRvdVxGKUYmTTdSMApJNlJUQUJMRV9TQVZFLzI1MTc3OTgzNlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKyRSbSpwayEjNSQhK2YnUTtSJkYpJCIraCdRO1ImRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2MjYxNzMwMFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKyl6ayRHaCEjNSQhK0ZuaShlJkYpJCIrRm5pKGUmRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI1MTU4MTkyMFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKyplckMmZiEjNSQhKyhwLz5vJkYpJCIrKHAvPm8mRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2MzQ1NzA3MlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKztlTmplISM1JCErLXcvR2RGKSQiKy13L0dkRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2MzIwNjE4NFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKypSPiY9ZSEjNSQhKyg+aTN2JkYpJCIrKD5pM3YmRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2MzY4NTA2NFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK0hMLid6JiEjNSQhK1JWP2lkRikkIitSVj9pZEYpRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2MzU2Mzg2MFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK1pOeCV5JiEjNSQhKy4hZnl3JkYpJCIrLiFmeXcmRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI1MjA5MzI0MFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKypcUiJ6ZCEjNSQhK0JBb3FkRikkIitCQW9xZEYpRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI1MTg5MDk2MFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK1Q5S3dkISM1JCErLkc0c2RGKSQiKy5HNHNkRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2Mzg1MjEwMFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK15AIlx4JiEjNSQhK015enNkRikkIitNeXpzZEYpRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI2MzY3NDgwMFgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK1N1P3VkISM1JCErJUddSngmRikkIislR11KeCZGKUYmTTdSMApJNlJUQUJMRV9TQVZFLzI2MzE5MjQ5MlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK3NdJlF4JiEjNSQhKyZcRUx4JkYpJCIrJlxFTHgmRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2MzMxNDUzNlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKyYpKXlPeCYhIzUkISsmZjlNeCZGKSQiKyZmOU14JkYpRiY=TTdSMApJNlJUQUJMRV9TQVZFLzI1MTk1NDAxMlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiKyJ6IWZ0ZCEjNSQhK1YnZU14JkYpJCIrVydlTXgmRilGJg==TTdSMApJNlJUQUJMRV9TQVZFLzI2MzU3NTAzMlgqJSlhbnl0aGluZ0c2IjYiW2dsISElISEhIiQiJCQiK1FuYXRkISM1JCErbTFbdGRGKSQiK20xW3RkRilGJg==