Numerical Methods: Mathematica Notebooks
(c) John H. Mathews, 1996

To accompany the text:

Numerical Methods:  for Mathematics, 
Science, and Engineering, 2nd Ed, 1992
Prentice Hall, Inc.
Englewood Cliffs, New Jersey 07632  U.S.A.
Prentice Hall, Inc.; USA, Canada, Mexico ISBN 0-13-624990-6
Prentice Hall, International Editions:      ISBN 0-13-625047-5

This free software is compliments of the author.
E-mail address:  in%"mathews@fullerton.edu"
http://titan.fullerton.edu/~mathews/

CONTENTS

Chapter 11.	Eigenvalues and Eigenvectors
	Algorithm 11.1	Power Method	
	Algorithm 11.2	Shifted Inverse Power Method	
	Algorithm 11.3	Jacobi Iteration for Eigenvalues and Eigenvectors	
	Algorithm 11.4	Reduction to Tridiagonal Form	
	Algorithm 11.5	The QL Method with Shifts

Algorithm 11.1 (Power Method).  To compute the dominant value  l1
and its associated eigenvector V1 for the  n by n  matrix A. It is assumed
that the n eigenvalues have the dominance property |l1| > |l2| >...> |ln|.
Section	11.2	The Power Method   Page  557

Off[General::spell1];
Clear[PowerMeth,A0,X0,epsilon,max1,show,V];
PowerMeth[A0_,X0_,epsilon_,max1_,show_] :=
Module[{A,Y,c1,cnt,dc,dv,err,state,done,iterating},
  Clear[maxsize,norm];
  maxsize[V_] := Module[{},
    vv = Sort[V];
    vmax = vv[[-1]];
    If[Abs[vv[[1]]]>Abs[vv[[-1]]],vmax=vv[[1]]];
    Return[vmax];];
  norm[V_] := Sqrt[V^2/.List->Plus]//N;
  A=A0//N;
  X=X0//N;
  lambda = 0;
  cnt = 0;
  err = 1;
  done = 0;
  iterating = 1;
  state = iterating;
  While[ (cnt<=max1) && (state==iterating),
    Y = A.X;
    c1 = maxsize[Y];
    dc = Abs[lambda-c1];
    Y = (1/c1) Y;
    dv = norm[X-Y];
    err = Max[{dc,dv}];
    X = Y;
    lambda = c1;
    state = done;
    If[err>epsilon, state = iterating];
    cnt = cnt+1;
    If[show==1, Print[N[{lambda,X},14]]; ]; ];
  V = X;
  Return[N[{lambda,X},14]]; ];
On[General::spell1];

Chapter 11, Example 11.6, Page 555.

Find the dominant eigenvalue and eigenvector for the matrix:

Clear[A,X];
A = {{0  , 11  , -5},
     {-2 ,  17 ,  -7},
     {-4 ,  26 , -10}};
A//TableForm
X = {1,1,1}

sol = PowerMeth[A,X,10^-10,25,1]

Now check our work.

lambda

X

A.X

lambda X

Compare with Mathematica's Eigensystem procedure.
Observe that Mathematica returns unit length eigenvectors.

matset = Eigensystem[N[A]];
matset[[1,1]]
matset[[2,1]]
X/norm[X]

Algorithm 11.2  (Shifted Inverse Power Method).  To compute the
dominant eigenvalue lj and its associated eigenvector Vj for the n by n matrix A.  It is assumed that the n eigenvalues are  l1 < l2 <...< ln  and  a  is a real number such that |lj - a| < |li - a|  for each  i=1,2, ...,j-1,j+1 ...,n.
Section	11.2	The Power Method   Page  558

First, execute the following group of cells:

Off[General::spell1];
Clear[LUfactor,a0,det,row];
LUfactor[a0_] :=
Module[{a,c,i,j,k,p},
  a = a0;
  mna = Dimensions[a];
  n  = mna[[1]];
  m = mna[[2]];
  det = 1;
  Clear[row];
  For[j=1,j<=n,j++, row[j]=j];
  For[p=1,p<=n-1,p++,
  Block[{},
    For[k=p+1,k<=n,k++,
      If[Abs[a[[row[k],p]]]>Abs[a[[row[p],p]]],
         t=row[p]; row[p]=row[k]; row[k]=t; det=-det];];
    det = det a[[row[p],p]];
    If[det == 0,
      Print["The matrix is singular."];
      Goto[GetOut]; ];
    Clear[k];
    For[k=p+1,k<=n,k++,
      rowk = row[k];
      rowp = row[p];
      a[[rowk,p]] = a[[rowk,p]]/a[[rowp,p]];
      For[c=p+1,c<=n,c++,
        a[[rowk,c]] = a[[rowk,c]] - 
                      a[[rowk,p]] a[[rowp,c]];];];];];
  Label[GetOut];
  det = det a[[row[n],n]];
  Return[a]; ]

Clear[SolveLU,lu0,b0]
SolveLU[lu0_,b0_] :=
Module[{a,c,i,j,k,p,sum},
  a = lu0;
  b = b0;
  x = Table[0,{i,n}];
  For[k=1,k<=n,k++,
    If[a[[row[k],k]]==0, Goto[GetOut] ]; ];
  x[[1]] = b[[row[1]]];
  Clear[k];
  For[k=2,k<=n,k++,
    sum = 0;
    rowk = row[k];
    For[c=1,c<=k-1,c++,
      sum = sum + a[[rowk,c]] x[[c]] ];
    x[[k]] = b[[rowk]] - sum; ];
  x[[n]] = x[[n]]/a[[row[n],n]];
  Clear[k];
  For[k=n-1,1<=k,k--,
    sum = 0;
    rowk = row[k];
    For[c=k+1,c<=n,c++,
      sum = sum + a[[rowk,c]] x[[c]] ];
    x[[k]] = (x[[k]] - sum)/a[[rowk,k]]; ];
  Label[GetOut];
  Return[x]; ];

Clear[InvPowerMeth,A0,alpha,X0,V,epsilon,max1,show];
InvPowerMeth[A0_,alpha_,X0_,epsilon_,max1_,show_] :=
Module[{A,Y,c1,cnt,dc,dv,err,state,done,iterating},
Clear[maxsize,norm];
maxsize[V_] := Module[{},
  vv = Sort[V];
  vmax = vv[[-1]];
  If[Abs[vv[[1]]]>Abs[vv[[-1]]],vmax=vv[[1]]];
  Return[vmax];];
norm[V_] := Sqrt[V^2/.List->Plus]//N;
n = Length[A0];
A = A0 - alpha*IdentityMatrix[n];
LU = LUfactor[A];
lambda = 0;
cnt = 0;
err = 1;
done = 0;
iterating = 1;
state = iterating;
X = X0;
While[(cnt<=max1)&&(state==iterating),
  Y = SolveLU[LU,X];  
  c1 = maxsize[Y];  
  dc = Abs[lambda - c1];
  Y = (1/c1)*Y;
  dv = norm[X-Y];
  err = Max[{dc,dv}];
  X = Y;
  lambda = c1; 
  state = done;
  If[err>epsilon, state = iterating];
  cnt = cnt+1;
  lamb = alpha + 1/c1;
  If[show==1, Print[N[{lambda,X},14]]; ]; ];
lambda = alpha + 1/c1;
Return[N[{lambda,X},14]] ];
On[General::spell1];

Chapter 11, Example 11.6, Page 555.

Find the dominant eigenvalue and eigenvector for the matrix:

Clear[alpha,A,X];
A = {{0  , 11  , -5},
     {-2 ,  17 ,  -7},
     {-4 ,  26 , -10}};
A//TableForm
X = {1,1,1}
alpha = 4.2

InvPowerMeth[A,alpha,X,0.0000001,50,1]

Now check our work.

lambda

X

A.X

lambda X

Compare with Mathematica's Eigensystem procedure.
Observe that Mathematica returns unit length eigenvectors.

matset = Eigensystem[N[A]];
matset[[1,1]]
matset[[2,1]]
X/norm[X]

Algorithm 11.3 (Jacobi Iteration for Eigenvalues and
Eigenvectors).  To compute the full set of eigen-pairs  { lj,Vj }j=1j=n
of the  n by n  real symmetric matrix A.  Jacobi iteration is used to find
all eigenvalue-eigenvector pairs.
Section	11.3	Jacobi's Method   Page  571

First, execute the following group of cells:

Off[General::spell1];
Clear[ZeroPQ,p,q];
ZeroPQ[p_,q_] :=
Module[{c,s,i,theta,XP,XQ,YP,YQ},
A = A//N;
theta = (A[[q,q]] - A[[p,p]])/(2*A[[p,q]]);
If[theta<10^19,
  T = 1/(Abs[theta] + Sqrt[theta^2 + 1]),
  T = 1/Abs[2*theta] ];      
If[theta<0, T = -T];
c = 1/Sqrt[T*T+1];
s = c*T;
n = Length[A];
XP = Table[0,{n}]; XQ = Table[0,{n}];
YP = Table[0,{n}]; YQ = Table[0,{n}];
For[i=1,i<=n,i++,
  XP[[i]] = A[[i,p]];
  XQ[[i]] = A[[i,q]]; ];
A[[p,q]] = 0;
A[[q,p]] = 0;
A[[p,p]] = c*c*XP[[p]] + s*s*XQ[[q]] - 2*c*s*XQ[[p]];
A[[q,q]] = s*s*XP[[p]] + c*c*XQ[[q]] + 2*c*s*XQ[[p]];
For[i=1,i<=n,i++,
  If[(i!=p)&&(i!=q),
  A[[i,p]] = c*XP[[i]] - s*XQ[[i]];
  A[[i,q]] = c*XQ[[i]] + s*XP[[i]];
  A[[p,i]] = A[[i,p]];
  A[[q,i]] = A[[i,q]]; ]; ];
For[i=1,i<=n,i++,
  YP[[i]] = V[[i,p]];
  YQ[[i]] = V[[i,q]]; ];
For[i=1,i<=n,i++,
  V[[i,p]] = c*YP[[i]] - s*YQ[[i]];
  V[[i,q]] = s*YP[[i]] + c*YQ[[i]]; ]; ];

Clear[MaxPQ];
MaxPQ :=
Module[{i,j,Pm,Qm},
n = Length[A];
MaxA = Abs[A[[1,n]]];
Pm = 1;
Qm = n ;     
For[i=1,i<=n,i++,
  For[j=i+1,j<=n,j++,        
	If[Abs[A[[i,j]]]>MaxA,
	  Pm = i;
	  Qm = j;
	  MaxA = Abs[A[[i,j]]]; ]; ]; ];      
p = Pm; q = Qm; Apq = A[[p,q]]; ];

RMS[A_] :=
Module[{j,sum,rms},
  sum = 0.0;
  n = Length[A];
  For[j=1,j<=n,j++,        
    sum = sum + A[[j,j]]^2 ];
  rms = Sqrt[sum/n];
Return[rms]; ];

Clear[OldJacobi];
OldJacobi[epsilon_]:=
Module[{Atemp,Stat,Done,Working},
  Atemp = A;
  V = IdentityMatrix[Length[A]];
  Stat = Working;       
  While[Stat==Working,
    MaxPQ;
    ZeroPQ[p,q];
    If[MaxA<epsilon*RMS[A], Stat = Done]; ];
  L = A;
  A = Atemp;
  Print[""];  
  Print[Chop[L]//TableForm]; 
  Print[""];  
  Print[Chop[V]//TableForm]; ];

Clear[CyclicJacobi,epsilon,MaxS];
CyclicJacobi[epsilon_,MaxS_]:=
Module[{Atemp,CountS},
Atemp = A;
CountS = 0;
V = IdentityMatrix[Length[A]];
Stat = Working;       
While[(Stat==Working) && (CountS<MaxS),
  CountS = CountS+1;
  t = RMS[A];
  Stat = Done;
  For[p=1,p<=n-1,p++,        
    For[q=p+1,q<=n,q++,
      If[Abs[A[[p,q]]]/t>epsilon,
        ZeroPQ[p,q];
        Stat = Working; ]; ]; ]; ];
  L = A;
  A = Atemp;
  Print[""];  
  Print[Chop[L]//TableForm]; 
  Print[""];  
  Print[Chop[V]//TableForm]; ];
On[General::spell1];

Chapter 11, Example 11.7, Page 569.

Use the original Jacobi method to find all the eigenvalues and 

eigenvectors of the matrix:

A = {{8,-1,3,-1},{-1,6,2,0},{3,2,9,1},{-1,0,1,7}}; 
A//TableForm

OldJacobi[10^-11]

Check our work by multiplying A by v1, etc.

l1 = L[[1,1]]

v1 = Transpose[V][[1]]

A.v1

l1 v1

Chapter 11, Example 11.7, Page 569.

Use the cyclic Jacobi method to find all the eigenvalues and eigenvectors 

of the matrix:

A = {{8,-1,3,-1},{-1,6,2,0},{3,2,9,1},{-1,0,1,7}}; 
A//TableForm

CyclicJacobi[10^-11,50]

Check our work by multiplying A by v1, etc.

l1 = L[[1,1]]

v1 = Transpose[V][[1]]

A.v1

l1 v1

Compare with Mathematica's Eigensystem procedure.
Observe that Mathematica returns unit length eigenvectors.

sol = Eigensystem[N[A]];
sol[[1]]
Transpose[sol[[2]]]//TableForm

Algorithm 11.4   (Reduction to Tridiagonal Form).  To reduce the
n x n  symmetric matrix A to tridiagonal form by using n-2 Householder
transformations.
Section	11.4	Eigenvalues of Symmetric Matrices   Page  581

First, execute the following group of cells:

Off[General::spell1];
Clear[FormW,k,n];
FormW[k_,n_] :=
Module[{i,j,sum},
  sum = 0;
  For[j=k+1,j<=n,j++,
    sum = sum + A1[[j,k]]^2; ];
  s0 = Sqrt[sum];
  If[A1[[k+1,k]]<0,s0 = -s0;];
  r0 = Sqrt[2*(A1[[k+1,k]] + s0)*s0];
  For[j=1,j<=k,j++,W1[[j]]=0;];
  W1[[k+1]] = (A1[[k+1,k]] + s0)/r0;
  For[j=k+2,j<=n,j++,
    W1[[j]] = A1[[j,k]]/r0;]; ];

Clear[FormV,k,n];
FormV[k_,n_] :=
Module[{i,j,sum},
  For[j=1,j<=k,j++,V1[[j]]=0;];
  For[i=k+1,i<=n,i++,
    sum = 0;
    For[j=k+1,j<=n,j++,
      sum = sum + A1[[i,j]]*W1[[j]];];
    V1[[i]] = sum;  ]; ];

Clear[FormQ,k,n];
FormQ[k_,n_] :=
Module[{j},
  c = 0;
  For[j=k+1,j<=n,j++,
    c = c + W1[[j]]*V1[[j]];];
  For[j=1,j<=k,j++,Q1[[j]]=0;];
  For[j=k+1,j<=n,j++,
    Q1[[j]] = V1[[j]]-c*W1[[j]];]; ];

Clear[FormA,k,n];
FormA[k_,n_] :=
Module[{i,j},
  For[j=k+2,j<=n,j++,
    A1[[j,k]] = 0;
    A1[[k,j]] = 0; ];
  A1[[k+1,k]] = -s0;
  A1[[k,k+1]] = -s0;
  For[j=k,j<=n,j++,
    A1[[j,j]] = A1[[j,j]]-4Q1[[j]]*W1[[j]];];
  For[i=k+1,i<=n,i++,
    For[j=i+1,j<=n,j++,
      A1[[i,j]] = A1[[i,j]]-2W1[[i]]*Q1[[j]]-2Q1[[i]]*W1[[j]];
      A1[[j,i]] = A1[[i,j]]; ]; ]; ];

Clear[Householder,A0,k,n];
Householder[A0_] :=
Module[{k},
  A1 = A0//N;
  n = Length[A1];
  Q1 = Table[0,{i,n}];
  V1 = Table[0,{i,n}];
  W1 = Table[0,{i,n}];
  For[k=1,k<=n-2,k++,
    FormW[k,n];
    FormV[k,n];
    FormQ[k,n];
    FormA[k,n];
    Print[" "]; ];
  Print[" "];
  Print[A1//TableForm]; ];
On[General::spell1];

Chapter 11, Example 11.8, Page 580.

Use Householder's method to reduce the following

symmetric matrix to tridiagonal form:

A = {{4,2,2,1},
     {2,-3,1,1},
     {2,1,3,1},
     {1,1,1,2}};
A//TableForm

Householder[A]

Algorithm 11.5  (The QL Method with Shifts).  Given the symmetric 
tridiagonal matrix A whose diagonal elements stored in  D(1),D(2),...,D(n)
and subdiagonal elements are stored in  E(1),E(2),...,E(n-1).  To construct
the eigenvalues  of  A  and store them in  D(1),D(2),...,D(n).
Section	11.4	Eigenvalues of Symmetric Matrices   Page  587

First, execute the following group of cells:

Off[General::spell1];
Clear[AtoTRI,n];
AtoTRI[n_] :=
Module[{j},
  For[j=1,j<=n,j++,
    D2[[j]] = A2[[j,j]]; ];
  For[j=1,j<=n-1,j++,
    E2[[j]] = A2[[j,j+1]]; ]; ];

Clear[TRItoA,n];
TRItoA[n_] :=
Module[{j},
  For[j=1,j<=n,j++,
    A2[[j,j]] = D2[[j]]; ];
  For[j=1,j<=n-1,j++,
    A2[[j,j+1]] = E2[[j]];
    A2[[j+1,j]] = E2[[j]]; ]; ];

Clear[FindShift,m];
FindShift[m_] :=
Module[{b1,c1,d1},
  b1 = -(D2[[m+1]] + D2[[m]]);
  c1 = D2[[m+1]]*D2[[m]] - E2[[m]]*E2[[m]];
  d1 = Sqrt[Abs[b1*b1 - 4*c1]];      
  If[b1>0, r1 = (-b1 - d1)/2;
           r2 = 2*c1/(-b1 - d1),
           r1 = (-b1 + d1)/2;
           r2 = 2*c1/(-b1 + d1);];
  sh = r1;
  If[Abs[D2[[m]]-r2]<Abs[D2[[m]]-r1],
    sh = r2; ]; ];

Clear[FormL,m,n];
FormL[n_,m_] :=
Module[{j,Pj1,Qj},
  For[j=1,j<=n,j++,
    P2[[j]] = D2[[j]];];  
  For[j=1,j<=n-1,j++,
    Q2[[j]] = E2[[j]];]; 
  Pj1 = P2[[n]];
  Qj = Q2[[n-1]];
  For[j=n-1,m<=j,j--,
    P2[[j+1]] = Sqrt[Pj1^2 + Q2[[j]]^2];
    C2[[j]] = Pj1/P2[[j+1]];
    S2[[j]] = Q2[[j]]/P2[[j+1]];
    Q2[[j]] = C2[[j]]*Qj + S2[[j]]*D2[[j]];
    Pj1 = C2[[j]]*D2[[j]] - S2[[j]]*Qj;
    If[j>m,Qj = C2[[j]]*Q2[[j-1]];]; ];
  P2[[m]] = Pj1; ];

Clear[FormA,m,n];
FormA[n_,m_] :=
Module[{j},
  D2[[n]] = S2[[n-1]]*Q2[[n-1]] + C2[[n-1]]*P2[[n]];
  E2[[n-1]] = S2[[n-1]]*P2[[n-1]];     
  For[j=n-2,m<=j,j--,
    D2[[j+1]]=S2[[j]]*Q2[[j]]+C2[[j]]*C2[[j+1]]*P2[[j+1]];
    E2[[j]] = S2[[j]]*P2[[j]]; ]; 
  D2[[m]] = C2[[m]]*P2[[m]]; ];

Clear[QLIteration,A0,delta,k,m,n];
QLIteration[A0_,delta_] :=
Module[{},
  A2 = A0//N;
  n = Length[A2];
  C2 = Table[0,{i,n}];
  D2 = Table[0,{i,n}];
  E2 = Table[0,{i,n}];
  P2 = Table[0,{i,n}];
  S2 = Table[0,{i,n}];
  Q2 = Table[0,{i,n}];
  AtoTRI[n];
  shift = 0;
  m = 1;
  While[m<=n-2,
    k = 1;
    Cond = 0;
    While[k<50&&Cond==0,
      FindShift[m];
      If[Abs[E2[[m]]]>delta,
      shift = shift + sh;
      For[j=m,j<=n,j++,
        D2[[j]] = D2[[j]] - sh;],
        D2[[m]] = D2[[m]] + shift;
        m = m+1;
        Cond = 1; ];
      FormL[n,m];
      FormA[n,m];
      k = k+1; ]; ];
  FindShift[m];
  D2[[n-1]] = r1 + shift;
  D2[[n]] = r2 + shift;
  For[j=1,j<=n-1,j++, E2[[j]]=0;];
  TRItoA[n]; 
  Return[A2//TableForm]; ];
On[General::spell1];

Chapter 11, Example 11.9, Page 584.

Use the QL method to find the eigenvalues of the symmetric tridiagonal 

matrix:

A1 = {{ 4., -3., 0, 0}, 
 
      {-3.,  2., 3.162277660168379332, 0}, 
 
      { 0 ,  3.162277660168379332, -1.4, -0.2}, 
      { 0 ,  0, -0.2, 1.4}};
A1//TableForm

QLIteration[A1,0.00000001]