218 7 Eigensystems
We will choose the angleθin order to haveb 23 =b 32 = 0. We get the new symmetric matrix
B=
a 11 a 12 c−a 13 s a 12 s+a 13 c
a 12 c−a 13 s a 22 c^2 +a 33 s^2 − 2 a 23 sc (a 22 −a 33 )sc+a 23 (c^2 −s^2 )
a 12 s+a 13 c(a 22 −a 33 )sc+a 23 (c^2 −s^2 ) a 22 s^2 +a 33 c^2 + 2 a 23 sc
.
Note thata 11 is unchanged! As it should.
We have then
b 11 =a 11
b 12 =a 12 cosθ−a 13 sinθ, 16 = 2 , 16 = 3
b 13 =a 13 cosθ+a 12 sinθ, 16 = 2 , 16 = 3
b 22 =a 22 cos^2 θ− 2 a 23 cosθsinθ+a 33 sin^2 θ
b 33 =a 33 cos^2 θ+ 2 a 23 cosθsinθ+a 22 sin^2 θ
b 23 = (a 22 −a 33 )cosθsinθ+a 23 (cos^2 θ−sin^2 θ)We will fix the angleθso thatb 23 = 0.
We get then a new matrix
B=
b 11 b 12 b 13
b 12 b 22 0
b 13 0 a 33
.
We repeat assuming thatb 12 is the largest non-diagonal matrix element and get a new matrix
C=
c−s 0
s c 0
0 0 1
b 11 b 12 b 13
b 12 b 22 0
b 13 0 b 33
c s 0
−s c 0
0 0 1
.
We continue this process till all non-diagonal matrix elements are zero. It is easy to con-
vince oneself that when performing the above operations, the matrix elementb 23 which was
previously set to zero may become different from zero. This is one of the problems which
slows down the Jacobi procedure. We leave this experience tothe reader in form of a large
numerical project at the end of this chapter.
An implementation of the above algorithm, normally referred to as the classical Jacobi
algorithm, is exposed partially in the code here.
http://folk.uio.no/compphys/programs/chapter07/cpp/jacobi.cpp
/
Jacobi's method for finding eigenvalues
eigenvectors of the symetric matrix A.
The eigenvalues of A will be on the diagonal
of A, with eigenvalue i being A[i][i].
The j-th component of the i-th eigenvector
is stored in R[i][j].
A: input matrix (n x n)
R: empty matrix for eigenvectors (n x n)
n: dimention of matrices
/
#include
#include
#include"jacobi.h"
voidjacobi_method (double A,doubleR,intn )