Page 125 - Numerical Methods for Chemical Engineering

114 3 Matrix eigenvalue analysis

Applying Gershgorin’s theorem to study the convergence of
iterative linear solvers

As a demonstration of the usefulness of Gershgorin’s theorem, we generate a convergence
criterion for the Jacobi iterative method of solving Ax = b. This example is typical of the
use of eigenvalues in numerical analysis, and also shows why the questions of eigenvector
basis set existence raised in the next section are of such importance.
We develop here a simple alternative to Gaussian elimination for solving a linear system
[0]
Ax = b. It is based on forming an initial guess of the solution, x , and iteratively refining
[1]
[2]
it to form a sequence x , x ,... that hopefully converges to a solution; i.e.,
& [k] &
lim & Ax − b = 0 (3.67)
&
k→∞
We add to each side of Ax = b the term Bx, for some non singular B,
Ax + Bx = b + Bx (3.68)

Upon rearrangement, this yields the identity
Bx = b + (B − A)x (3.69)
[k]
and suggests a rule for forming a new guess x [k+1] from x ,
Bx [k+1] = b + (B − A)x [k] (3.70)

If B = A, this rule yields the exact solution after only one iteration, but it is equivalent
to solving Ax = b directly in the first place. We therefore want to choose some B that
approximates A such that: (1) B − A is small and (2) the update linear system is easy to
solve (e.g. B is triangular). In the Jacobi method, we simply choose B to be the diagonal
part of A,

 
a 11
a 22
 
 
. (3.71)
B =  . 
 . 
a NN
Under what conditions must the Jacobi method converge? To answer this question, let us
define the error vector at iteration k as
−1
ε [k] ≡ x [k] − x x = A b (3.72)
[k]
We want lim k→∞ ε = 0. We obtain a rule for the transformation of the error vector at
each iteration by subtracting (3.69) from (3.70),

Bε [k+1] = (B − A)ε [k] ⇒ ε [k+1] = B −1 (B − A)ε [k] (3.73)
[k]
We now use Gershgorin’s theorem to find when lim k→∞ ε = 0. Let us assume that the
[2]
[1]
matrix B −1 (B − A) has a set of eigenvectors {w , w ,..., w [N] }, satisfying B −1 (B −
N
[ j]
A)w [ j] = λ j w , that is linearly independent and forms a basis for . We can then express