Mathematical Principles for Scientific Computing and Visualization  
http://www.farinhansford.com/books/scv
by Gerald Farin and Dianne Hansford


This notebook finds the eigenvalues of a symmetric matrix using Jacobi iteration as described in Section 6.4.

Given: a symmetric matrix a.
First we find the eigenvalues using Mathematica's Eigenvalue function.
We run  Jacobi iteration and compare the results to the Mathematica results.

See Chapter 6: Eigen-Problems for details.

Updated July 2008

a={{3.0,-1,1},{-1,4,2},{1,2,-3}};
Print["Matrix A = ",MatrixForm[a]]

Print["Eigenvalues (from Mathematica's function) = ",Eigenvalues[a]]

(* Compute eigenvalues with Jacobi method *)

(* Largest off-diagonal element in matrix a is a_{2,3} *)
i=2;
j=3;
(* Compute angle *)
tan2al=2a[[i,j]]/(a[[i,i]]-a[[j,j]])//N;
alpha=ArcTan[tan2al]/2;
c=Cos[alpha];
s=Sin[alpha];
u={{1.0,0,0},{0,1,0},{0,0,1}};
u[[i,j]]=-s;
u[[j,i]]=s;
u[[i,i]]=c;
u[[j,j]]=c;

(* Form new matrix a^(1) -- first iteration *)
a=Transpose[u].a.u;

Print["Jacobi eigenvalue approximations are on the diagonal of the resulting matrix"]
Print["Stop when off-diagonal elements are smaller than a given tolerance"]
Print["1st iteration = ",MatrixForm[a]]

(* ITERATION 2 *)
u={{1.0,0,0},{0,1,0},{0,0,1}};
(* largest off-diagonal element from iteratioin 1 is (1,3) *)
i=1;
j=3;
tan2al=2a[[i,j]]/(a[[i,i]]-a[[j,j]])//N;
alpha=ArcTan[tan2al]/2;
c=Cos[alpha];
s=Sin[alpha];
u[[i,j]]=-s;
u[[j,i]]=s;
u[[i,i]]=c;
u[[j,j]]=c;
a=Transpose[u].a.u;
Print["2nd iteration = ",MatrixForm[a]]

(* ITERATION 3 *)
u={{1.0,0,0},{0,1,0},{0,0,1}};
(* largest off-diagonal element from iteration 2 is (1,2) *)
i=1;
j=2;
tan2al=2a[[i,j]]/(a[[i,i]]-a[[j,j]])//N;
alpha=ArcTan[tan2al]/2;
c=Cos[alpha];
s=Sin[alpha];
u[[i,j]]=-s;
u[[j,i]]=s;
u[[i,i]]=c;
u[[j,j]]=c;
a=Transpose[u].a.u;
Print["3rd iteration = ",MatrixForm[a]]

(* ITERATION 4 *)

u={{1.0,0,0},{0,1,0},{0,0,1}};
i=2;
j=3;
tan2al=2a[[i,j]]/(a[[i,i]]-a[[j,j]])//N;
alpha=ArcTan[tan2al]/2;
c=Cos[alpha];
s=Sin[alpha];
u[[i,j]]=-s;
u[[j,i]]=s;
u[[i,i]]=c;
u[[j,j]]=c;
a=Transpose[u].a.u;
Print["4th iteration = ",MatrixForm[a]]

(* ITERATION 5 *)
u={{1.0,0,0},{0,1,0},{0,0,1}};
i=1;
j=3;
tan2al=2a[[i,j]]/(a[[i,i]]-a[[j,j]])//N;
alpha=ArcTan[tan2al]/2;
c=Cos[alpha];
s=Sin[alpha];
u[[i,j]]=-s;
u[[j,i]]=s;
u[[i,i]]=c;
u[[j,j]]=c;
a=Transpose[u].a.u;
Print["5th iteration = ",MatrixForm[a]]

Print["Let's say we met stopping criterion of all off-diagonal elements < tol"]
Print["Eigenvalue estimates: ",a[[1,1]],"  ",a[[2,2]],"  ",a[[3,3]]]
Print["These match the Mathematica eigenvalues exactly to 6 digits"]