logo
Tsung-Ming Huang
Department of Mathematics National Taiwan Normal University, Taiwan
September 12, 2015
logo
Outline
1 Norms of vectors and matrices
2 Eigenvalues and eigenvectors
3 The Jacobi and Gauss-Siedel Iterative Techniques
4 Relaxation Techniques for Solving Linear Systems
5 Error bounds and iterative refinement
6 The conjugate gradient method
2 / 99
logo
Definition 1
k · k : Rn→ R is a vector norm if (i) kxk ≥ 0, ∀ x ∈ Rn,
(ii) kxk = 0 if and only if x = 0,
(iii) kαx| = |α|kxk ∀ α ∈ R and x ∈ Rn, (iv) kx + yk ≤ kxk + kyk ∀ x, y ∈ Rn. Definition 2
The `2 and `∞norms for x = [x1, x2, · · · , xn]T are defined by
kxk2 = (xTx)1/2= ( n
X
i=1
x2i )1/2
and kxk∞= max
1≤i≤n|xi|.
The `2 norm is also called the Euclidean norm.
logo
Theorem 3 (Cauchy-Bunyakovsky-Schwarz inequality) For each x = [x1, x2, · · · , xn]T and y = [y1, y2, · · · , yn]T in Rn,
xTy =
n
X
i=1
xiyi ≤ ( n
X
i=1
x2i
)1/2( n X
i=1
y2i )1/2
= kxk2· kyk2.
Proof: If x = 0 or y = 0, the result is immediate.
Suppose x 6= 0 and y 6= 0. For each α ∈ R, 0 ≤ kx − αyk22 =
n
X
i=1
(xi− αyi)2=
n
X
i=1
x2i − 2α
n
X
i=1
xiyi+ α2
n
X
i=1
yi2,
and 2α
n
X
i=1
xiyi≤
n
X
i=1
x2i + α2
n
X
i=1
yi2= kxk22+ α2kyk22.
4 / 99
logo
Since kxk2> 0and kyk2 > 0, we can let α = kxk2
kyk2 to give
2kxk2
kyk2
n
X
i=1
xiyi
!
≤ kxk22+kxk22
kyk22kyk22= 2kxk22. Thus
xTy =
n
X
i=1
xiyi ≤ kxk2kyk2.
logo
For each x, y ∈ Rn, kx + yk∞ = max
1≤i≤n|xi+ yi| ≤ max
1≤i≤n(|xi| + |yi|)
≤ max
1≤i≤n|xi| + max
1≤i≤n|yi| = kxk∞+ kyk∞
and
kx + yk22 =
n
X
i=1
(xi+ yi)2=
2
X
i=1
x2i + 2
n
X
i=1
xiyi+
n
X
i=1
yi2
≤ kxk22+ 2kxk2kyk2+ kyk22 = (kxk2+ kyk2)2, which gives
kx + yk2 ≤ kxk2+ kyk2.
6 / 99
logo
Definition 4
A sequence {x(k)∈ Rn}∞k=1is convergent to x with respect to the norm k · k if ∀ ε > 0, ∃ an integer N (ε) such that
kx(k)− xk < ε, ∀ k ≥ N (ε).
Theorem 5
{x(k)∈ Rn}∞k=1 converges to x with respect to k · k∞if and only if
k→∞lim x(k)i = xi, ∀ i = 1, 2, . . . , n.
Proof: “⇒” Given any ε > 0, ∃ an integer N (ε) such that
1≤i≤nmax |x(k)i − xi| = kx(k)− xk∞< ε, ∀ k ≥ N (ε).
logo
This result implies that
|x(k)i − xi| < ε, ∀ i = 1, 2, . . . , n.
Hence
k→∞lim x(k)i = xi, ∀ i.
“⇐” For a given ε > 0, let Ni(ε)represent an integer with
|x(k)i − xi| < ε, whenever k ≥ Ni(ε).
Define
N (ε) = max
1≤i≤nNi(ε).
If k ≥ N (ε), then
1≤i≤nmax |x(k)i − xi| = kx(k)− xk∞< ε.
This implies that {x(k)} converges to x with respect to k · k∞.
8 / 99
logo
Theorem 6 For each x ∈ Rn,
kxk∞≤ kxk2 ≤√
nkxk∞. Proof: Let xj be a coordinate of x such that
kxk2∞= |xj|2≤
n
X
i=1
x2i = kxk22,
so kxk∞≤ kxk2 and kxk22 =
n
X
i=1
x2i ≤
n
X
i=1
x2j = nx2j = nkxk2∞,
so kxk2 ≤√
nkxk∞.
logo
Definition 7
A matrix norm k · k on the set of all n × n matrices is a real-valued function satisfying for all n × n matrices A and B and all real number α:
(i) kAk ≥ 0;
(ii) kAk = 0 if and only if A = 0;
(iii) kαAk = |α|kAk;
(iv) kA + Bk ≤ kAk + kBk;
(v) kABk ≤ kAkkBk;
Theorem 8
If k · k is a vector norm on Rn, then kAk = max
kxk=1kAxk
is a matrix norm. 10 / 99
logo
For any z 6= 0, we have x = z/kzk as a unit vector. Hence kAk = max
kxk=1kAxk = max
z6=0
A
z kzk
= max
z6=0
kAzk kzk .
Corollary 9
kAzk ≤ kAk · kzk.
Theorem 10
If A = [aij]is an n × n matrix, then kAk∞= max
1≤i≤n n
X
j=1
|aij|.
logo
Proof: Let x be an n-dimension vector with 1 = kxk∞= max
1≤i≤n|xi|.
Then
kAxk∞ = max
1≤i≤n
n
X
j=1
aijxj
≤ max
1≤i≤n n
X
j=1
|aij| max
1≤j≤n|xj| = max
1≤i≤n n
X
j=1
|aij|.
Consequently,
kAk∞= max
kxk∞=1kAxk∞≤ max
1≤i≤n n
X
j=1
|aij|.
On the other hand, let p be an integer with
n
X
j=1
|apj| = max
1≤i≤n n
X
j=1
|aij|,
12 / 99
logo
and x be the vector with xj =
1, if apj≥ 0,
−1, if apj< 0.
Then
kxk∞= 1 and apjxj = |apj|, ∀ j = 1, 2, . . . , n, so
kAxk∞= max
1≤i≤n
n
X
j=1
aijxj
≥
n
X
j=1
apjxj
=
n
X
j=1
|apj|
= max
1≤i≤n n
X
j=1
|aij|.
This result implies that kAk∞= max
kxk∞=1
kAxk∞≥ max
1≤i≤n n
X
j=1
|aij|.
which gives
kAk∞= max
1≤i≤n n
X|aij|. 13 / 99
logo
Exercise
Page 441: 5, 9, 10, 11
14 / 99
logo
Eigenvalues and eigenvectors
Definition 11 (Characteristic polynomial)
If A is a square matrix, the characteristic polynomial of A is defined by p(λ) =det(A − λI).
Definition 12 (Eigenvalue and eigenvector)
If p is the characteristic polynomial of the matrix A, the zeros of p are eigenvalues of the matrix A. If λ is an eigenvalue of A and x 6= 0 satisfies (A − λI)x = 0, then x is an eigenvector of A corresponding to the eigenvalue λ.
Definition 13 (Spectrum and Spectral Radius)
The set of all eigenvalues of a matrix A is called the spectrum of A.
The spectral radius of A is
ρ(A) = max{|λ|; λis an eigenvalue ofA}.
logo
Theorem 14
If A is an n × n matrix, then (i) kAk2 =p
ρ(ATA);
(ii) ρ(A) ≤ kAkfor any matrix norm.
Proof: Proof for the second part. Suppose λ is an eigenvalue of Aand x 6= 0 is a corresponding eigenvector such that Ax = λx and kxk = 1. Then
|λ| = |λ|kxk = kλxk = kAxk ≤ kAkkxk = kAk, that is, |λ| ≤ kAk. Since λ is arbitrary, this implies that ρ(A) = max |λ| ≤ kAk.
Theorem 15
For any A and anyε > 0, there exists a matrix norm k · k such that
ρ(A) < kAk < ρ(A) + ε.
16 / 99
logo
Definition 16
We call an n × n matrix A convergent if
k→∞lim(Ak)ij = 0 ∀ i = 1, 2, . . . , n and j = 1, 2, . . . , n.
Theorem 17
The following statements are equivalent.
1 Ais aconvergent matrix;
2 lim
k→∞kAkk = 0forsomematrix norm;
3 lim
k→∞kAkk = 0forallmatrix norm;
4 ρ(A) < 1;
5 lim
k→∞Akx = 0foranyx.
logo
Exercise
Page 449: 11, 12, 18, 19
18 / 99
logo
Jacobi and Gauss-Siedel Iterative Techniques
For small dimension of linear systems, it requires for direct techniques.
For large systems, iterative techniques are efficient in terms of both computer storage and computation.
The basic idea of iterative techniques is to split the coefficient matrix A into
A = M − (M − A),
for some matrix M , which is called thesplitting matrix. Here we assume that A and M are bothnonsingular. Then the original problem is rewritten in the equivalent form
M x = (M − A)x + b.
logo
This suggests an iterative process
x(k)= (I − M−1A)x(k−1)+ M−1b ≡ T x(k−1)+ c,
whereT is usually called theiteration matrix. The initial vector x(0)can be arbitrary or be chosen according to certain
conditions.
Two criteria for choosing the splitting matrix M are x(k)is easily computed. More precisely, the system M x(k)= y is easy to solve;
the sequence {x(k)} converges rapidly to the exact solution.
Note that one way to achieve the second goal is to choose M so that M−1approximate A−1.
In the following subsections, we will introduce some of the mostly commonly used classic iterative methods.
20 / 99
logo
Jacobi Method
If we decompose the coefficient matrix A as A = L + D + U,
whereDis thediagonal part,Lis thestrictly lower triangular part, andU is thestrictly upper triangular part, ofA, and chooseM = D, then we derive the iterative formulation for Jacobi method:
x(k)= −D−1(L + U )x(k−1)+ D−1b.
With this method, the iteration matrixTJ = −D−1(L + U )and c = D−1b. Each component x(k)i can be computed by
x(k)i =
bi−
i−1
X
j=1
aijx(k−1)j −
n
X
j=i+1
aijx(k−1)j
,
aii.
logo
a11x(k)1 + a12x(k−1)2 + a13x(k−1)3 + · · · + a1nx(k−1)n = b1
a21x(k−1)1 +a22x(k)2 + a23x(k−1)3 + · · · + a2nx(k−1)n = b2 ... an1x(k−1)1 + an2x(k−1)2 + an3x(k−1)3 + · · · +annx(k)n = bn. Algorithm 1 (Jacobi Method)
Given x(0), tolerance T OL, maximum number of iteration M . Set k = 1.
While k ≤ M and kx − x(0)k2 ≥ T OL Set k = k + 1, x(0) = x.
For i = 1, 2, . . . , n xi=
bi−Pi−1
j=1aijx(0)j −Pn
j=i+1aijx(0)j
,
aii
End For
End While 22 / 99
logo
Example 18
Consider the linear system Ax = b given by
E1 : 10x1 − x2 + 2x3 = 6,
E2 : −x1 + 11x2 − x3 + 3x4 = 25, E3 : 2x1 − x2 + 10x3 − x4 = −11,
E4 : 3x2 − x3 + 8x4 = 15
which has the unique solution x = [1, 2, −1, 1]T. Solving equation Ei for xi, for i = 1, 2, 3, 4, we obtain
x1 = 1/10x2 − 1/5x3 + 3/5,
x2 = 1/11x1 + 1/11x3 − 3/11x4 + 25/11,
x3 = −1/5x1 + 1/10x2 + 1/10x4 − 11/10,
x4 = − 3/8x2 + 1/8x3 + 15/8.
logo
Then Ax = b can be rewritten in the form x = T x + c with
T =
0 1/10 −1/5 0
1/11 0 1/11 −3/11
−1/5 1/10 0 1/10
0 −3/8 1/8 0
and c =
3/5 25/11
−11/10 15/8
and the iterative formulation for Jacobi method is x(k)= T x(k−1)+ c for k = 1, 2, . . . .
The numerical results of such iteration is list as follows:
24 / 99
logo
k x1 x2 x3 x4
0 0.0000 0.0000 0.0000 0.0000 1 0.6000 2.2727 -1.1000 1.8750 2 1.0473 1.7159 -0.8052 0.8852 3 0.9326 2.0533 -1.0493 1.1309 4 1.0152 1.9537 -0.9681 0.9738 5 0.9890 2.0114 -1.0103 1.0214 6 1.0032 1.9922 -0.9945 0.9944 7 0.9981 2.0023 -1.0020 1.0036 8 1.0006 1.9987 -0.9990 0.9989 9 0.9997 2.0004 -1.0004 1.0006 10 1.0001 1.9998 -0.9998 0.9998
logo
Matlab code of Example
clear all; delete rslt.dat; diary rslt.dat; diary on;
n = 4; xold = zeros(n,1); xnew = zeros(n,1); T = zeros(n,n);
T(1,2) = 1/10; T(1,3) = -1/5; T(2,1) = 1/11;
T(2,3) = 1/11; T(2,4) = -3/11; T(3,1) = -1/5;
T(3,2) = 1/10; T(3,4) = 1/10; T(4,2) = -3/8; T(4,3) = 1/8;
c(1,1) = 3/5; c(2,1) = 25/11; c(3,1) = -11/10; c(4,1) = 15/8;
xnew = T * xold + c; k = 0;
fprintf(’ k x1 x2 x3 x4 \n’);
while ( k <= 100 & norm(xnew-xold) > 1.0d-14 ) xold = xnew; xnew = T * xold + c; k = k + 1;
fprintf(’%3.0f ’,k);
for jj = 1:n
fprintf(’%5.4f ’,xold(jj));
end
fprintf(’\n’);
end 26 / 99
logo
Gauss-Seidel Method
When computing x(k)i for i > 1, x(k)1 , . . . , x(k)i−1have already been computed and are likely to be better approximations to the exact x1, . . . , xi−1than x(k−1)1 , . . . , x(k−1)i−1 . It seems reasonable to compute x(k)i using these most recently computed values. That is
a11x(k)1 + a12x(k−1)2 + a13x(k−1)3 + · · · + a1nx(k−1)n = b1
a21x(k)1 +a22x(k)2 + a23x(k−1)3 + · · · + a2nx(k−1)n = b2
a31x(k)1 +a32x(k)2 +a33x(k)3 + · · · + a3nx(k−1)n = b3
... an1x(k)1 +an2x(k)2 +an3x(k)3 + · · · +annx(k)n = bn. This improvement induce the Gauss-Seidel method.
The Gauss-Seidel method setsM = D + Land defines the iteration as
x(k)= −(D + L)−1U x(k−1)+ (D + L)−1b.
logo
That is, Gauss-Seidel method usesTG= −(D + L)−1U as the iteration matrix. The formulation above can be rewritten as
x(k)= −D−1
Lx(k)+ U x(k−1)− b .
Hence each component x(k)i can be computed by
x(k)i =
bi−
i−1
X
j=1
aijx(k)j −
n
X
j=i+1
aijx(k−1)j
,
aii.
For Jacobi method, only the components of x(k−1)are used to compute x(k). Hence x(k)i , i = 1, . . . , n,can be computed in parallel at each iteration k.
At each iteration of Gauss-Seidel method, since x(k)i can not be computed until x(k)1 , . . . , x(k)i−1are available, the method is not a parallel algorithm in nature.
28 / 99
logo
Algorithm 2 (Gauss-Seidel Method)
Given x(0), tolerance T OL, maximum number of iteration M . Set k = 1.
For i = 1, 2, . . . , n xi=
bi−Pi−1
j=1aijxj−Pn
j=i+1aijx(0)j
,
aii
End For
While k ≤ M and kx − x(0)k2 ≥ T OL Set k = k + 1, x(0) = x.
For i = 1, 2, . . . , n xi=
bi−Pi−1
j=1aijxj−Pn
j=i+1aijx(0)j ,
aii End For
End While
logo
Example 19
Consider the linear system Ax = b given by
E1 : 10x1 − x2 + 2x3 = 6,
E2 : −x1 + 11x2 − x3 + 3x4 = 25, E3 : 2x1 − x2 + 10x3 − x4 = −11,
E4 : 3x2 − x3 + 8x4 = 15
which has the unique solution x = [1, 2, −1, 1]T. Gauss-Seidel method gives the equation
x(k)1 = 101x(k−1)2 − 15x(k−1)3 + 35, x(k)2 = 111x(k)1 + 111x(k−1)3 − 113x(k−1)4 + 2511, x(k)3 = −15x(k)1 + 101x(k)2 + 101x(k−1)4 − 1110, x(k)4 = − 38x(k)2 + 18x(k)3 + 158.
30 / 99
logo
The numerical results of such iteration is list as follows:
k x1 x2 x3 x4
0 0.0000 0.0000 0.0000 0.0000 1 0.6000 2.3273 -0.9873 0.8789 2 1.0302 2.0369 -1.0145 0.9843 3 1.0066 2.0036 -1.0025 0.9984 4 1.0009 2.0003 -1.0003 0.9998 5 1.0001 2.0000 -1.0000 1.0000 The results of Example appear to imply that the
Gauss-Seidel method is superior to the Jacobi method.
This is almost always true, but there are linear systems for which the Jacobi method converges and the Gauss-Seidel method does not.
See Exercises 17 and 18 (8th edition).
logo Matlab code of Example
clear all; delete rslt.dat; diary rslt.dat; diary on;
n = 4; xold = zeros(n,1); xnew = zeros(n,1); A = zeros(n,n);
A(1,1)=10; A(1,2)=-1; A(1,3)=2; A(2,1)=-1; A(2,2)=11; A(2,3)=-1; A(2,4)=3; A(3,1)=2; A(3,2)=-1;
A(3,3)=10; A(3,4)=-1; A(4,2)=3; A(4,3)=-1; A(4,4)=8; b(1)=6; b(2)=25; b(3)=-11; b(4)=15;
for ii = 1:n xnew(ii) = b(ii);
for jj = 1:ii-1
xnew(ii) = xnew(ii) - A(ii,jj) * xnew(jj);
end for jj = ii+1:n
xnew(ii) = xnew(ii) - A(ii,jj) * xold(jj);
end
xnew(ii) = xnew(ii) / A(ii,ii);
end
k = 0; fprintf(’ k x1 x2 x3 x4 \n’);
while ( k <= 100 & norm(xnew-xold) > 1.0d-14 ) xold = xnew; k = k + 1;
for ii = 1:n xnew(ii) = b(ii);
for jj = 1:ii-1
xnew(ii) = xnew(ii) - A(ii,jj) * xnew(jj);
end for jj = ii+1:n
xnew(ii) = xnew(ii) - A(ii,jj) * xold(jj);
end
xnew(ii) = xnew(ii) / A(ii,ii);
end
fprintf(’%3.0f ’,k);
for jj = 1:n
fprintf(’%5.4f ’,xold(jj));
end fprintf(’\n’);
end
diary off 32 / 99
logo
Lemma 20
Ifρ(T ) < 1, then(I − T )−1exists and (I − T )−1=
∞
X
i=0
Ti= I + T + T2+ · · · .
Proof: Let λ be an eigenvalue of T , then 1 − λ is an eigenvalue of I − T . But |λ| ≤ ρ(A) < 1, so 1 − λ 6= 0 and 0 is not an eigenvalue of I − T , which means (I − T ) is nonsingular.
Next we show that (I − T )−1= I + T + T2+ · · ·. Since (I − T )
m
X
i=0
Ti
!
= I − Tm+1, and ρ(T ) < 1 implies kTmk → 0 as m → ∞, we have
(I − T ) lim
m→∞
m
XTi
!
= (I − T )
∞
XTi
!
= I.
logo
Theorem 21
Foranyx(0)∈ Rn , the sequence produced by x(k)= T x(k−1)+ c, k = 1, 2, . . . ,
convergesto theuniquesolution of x = T x + c if and only if ρ(T ) < 1.
Proof: Suppose ρ(T ) < 1. The sequence of vectors x(k)produced by the iterative formulation are
x(1)= T x(0)+ c
x(2)= T x(1)+ c = T2x(0)+ (T + I)c x(3)= T x(2)+ c = T3x(0)+ (T2+ T + I)c
... In general
x(k)= Tkx(0)+ (Tk−1+ Tk−2+ · · · + T + I)c.
34 / 99
logo
Since ρ(T ) < 1, limk→∞Tkx(0)= 0for any x(0) ∈ Rn. By Lemma 20,
(Tk−1+ Tk−2+ · · · + T + I)c → (I − T )−1c, as k → ∞.
Therefore
k→∞lim x(k)= lim
k→∞Tkx(0)+
∞
X
j=0
Tj
c = (I − T )−1c.
Conversely, suppose {x(k)} → x = (I − T )−1c. Since
x − x(k)= T x + c − T x(k−1)− c = T (x − x(k−1)) = T2(x − x(k−2))
= · · · = Tk(x − x(0)).
Let z = x − x(0). Then
k→∞lim Tkz = lim
k→∞(x − x(k)) = 0.
It follows from theorem ρ(T ) < 1.
logo
Theorem 22
IfkT k < 1, then the sequencex(k)convergestoxforanyinitial x(0)and
1 kx − x(k)k ≤ kT kkkx − x(0)k
2 kx − x(k)k ≤ 1−kT kkT kk kx(1)− x(0)k.
Proof: Since x = T x + c and x(k) = T x(k−1)+ c, x − x(k)= T x + c − T x(k−1)− c
= T (x − x(k−1))
= T2(x − x(k−2)) = · · · = Tk(x − x(0)).
The first statement can then be derived
kx − x(k)k = kTk(x − x(0))k ≤ kT kkkx − x(0)k.
For the second result, we first show that
kx(n)− x(n−1)k ≤ kT kn−1kx(1)− x(0)k for any n ≥ 1. 36 / 99
logo
Since
x(n)− x(n−1)= T x(n−1)+ c − T x(n−2)− c
= T (x(n−1)− x(n−2))
= T2(x(n−2)− x(n−3)) = · · · = Tn−1(x(1)− x(0)), we have
kx(n)− x(n−1)k ≤ kT kn−1kx(1)− x(0)k.
Let m ≥ k, x(m)− x(k)
=
x(m)− x(m−1) +
x(m−1)− x(m−2)
+ · · · +
x(k+1)− x(k)
=Tm−1
x(1)− x(0)
+ Tm−2
x(1)− x(0)
+ · · · + Tk
x(1)− x(0)
=
Tm−1+ Tm−2+ · · · + Tk
x(1)− x(0) ,
logo
hence
kx(m)− x(k)k
≤
kT km−1+ kT km−2+ · · · + kT kk
kx(1)− x(0)k
=kT kk
kT km−k−1+ kT km−k−2+ · · · + 1
kx(1)− x(0)k.
Since limm→∞x(m) = x, kx − x(k)k
= lim
m→∞kx(m)− x(k)k
≤ lim
m→∞kT kk
kT km−k−1+ kT km−k−2+ · · · + 1
kx(1)− x(0)k
=kT kkkx(1)− x(0)k lim
m→∞
kT km−k−1+ kT km−k−2+ · · · + 1
=kT kk 1
1 − kT kkx(1)− x(0)k.
This proves the second result.
38 / 99
logo
Theorem 23
If A isstrictly diagonal dominant, then both theJacobiand Gauss-Seidelmethodsconvergesforanyinitial vector x(0).
Proof: By assumption, A is strictly diagonal dominant, hence aii6= 0 (otherwise A is singular) and
|aii| >
n
X
j=1,j6=i
|aij|, i = 1, 2, . . . , n.
For Jacobi method, the iteration matrix TJ= −D−1(L + U )has entries
[TJ]ij = (−aaij
ii, i 6= j,
0, i = j.
Hence
kTJk∞= max
1≤i≤n n
X
j=1,j6=i
aij aii
= max
1≤i≤n
1
|aii|
n
X
j=1,j6=i
|aij| < 1,
and this implies that the Jacobi method converges.
logo
For Gauss-Seidel method, the iteration matrix
TG= −(D + L)−1U. Let λ be any eigenvalue of TGand y, kyk∞= 1, is a corresponding eigenvector. Thus
TGy = λy =⇒ −U y = λ(D + L)y.
Hence for i = 1, . . . , n,
−
n
X
j=i+1
aijyj = λaiiyi+ λ
i−1
X
j=1
aijyj. This gives
λaiiyi = −λ
i−1
X
j=1
aijyj−
n
X
j=i+1
aijyj
and
|λ||aii||yi| ≤ |λ|
i−1
X
j=1
|aij||yj| +
n
X
j=i+1
|aij||yj|.
40 / 99
logo
Choose the index k such that |yk| = 1 ≥ |yj| (this index can always be found since kyk∞= 1). Then
|λ||akk| ≤ |λ|
k−1
X
j=1
|akj| +
n
X
j=k+1
|akj|
which gives
|λ| ≤ Pn
j=k+1|akj|
|akk| −Pk−1
j=1|akj| <
Pn
j=k+1|akj| Pn
j=k+1|akj| = 1
Since λ is arbitrary, ρ(TG) < 1. This means the Gauss-Seidel method converges.
The rate of convergence depends on the spectral radius of the matrix associated with the method.
One way to select a procedure to accelerate convergence is to choose a method whose associated matrix has minimal spectral radius.
logo
Exercise
Page 459: 9, 10, 11
42 / 99
logo
Relaxation Techniques for Solving Linear Systems
Definition 24
Suppose ˜x ∈ Rn is an approximated solution of Ax = b. Theresidual vectorrfor ˜xisr = b − A˜x.
Let the approximate solution x(k,i)produced by Gauss-Seidel method be defined by
x(k,i)=h
x(k)1 , . . . , x(k)i−1, x(k−1)i , . . . , x(k−1)n iT
and
r(k)i =h
r(k)1i , r2i(k), . . . , rni(k)iT
= b − Ax(k,i)
be the corresponding residual vector. Then the mth component of ri(k)is
r(k)mi = bm−
i−1
X
j=1
amjx(k)j −
n
X
j=i
amjx(k−1)j ,
logo
or, equivalently,
r(k)mi = bm−
i−1
X
j=1
amjx(k)j −
n
X
j=i+1
amjx(k−1)j − amix(k−1)i ,
for each m = 1, 2, . . . , n.
In particular, the ith component of ri(k)is
r(k)ii = bi−
i−1
X
j=1
aijx(k)j −
n
X
j=i+1
aijx(k−1)j − aiix(k−1)i , so
aiix(k−1)i + r(k)ii = bi−
i−1
X
j=1
aijx(k)j −
n
X
j=i+1
aijx(k−1)j
= aiix(k)i .
44 / 99
logo
Consequently, the Gauss-Seidel method can be characterized as choosing x(k)i to satisfy
x(k)i = x(k−1)i +r(k)ii aii
.
Relaxation method is modified the Gauss-Seidel procedure to x(k)i = x(k−1)i + ωrii(k)
aii
= x(k−1)i + ω aii
bi−
i−1
X
j=1
aijx(k)j −
n
X
j=i+1
aijx(k−1)j − aiix(k−1)i
= (1 − ω)x(k−1)i + ω aii
bi−
i−1
X
j=1
aijx(k)j −
n
X
j=i+1
aijx(k−1)j
(1) for certain choices of positive ω such that the norm of the
residual vector is reduced and the convergence is significantly
logo
These methods are called for
ω < 1: under relaxation, ω = 1: Gauss-Seidel method, ω > 1: over relaxation.
Over-relaxation methods are calledSOR (Successive
over-relaxation). To determine the matrix of the SOR method, we rewrite (1) as
aiix(k)i + ω
i−1
X
j=1
aijx(k)j = (1 − ω)aiix(k−1)i − ω
n
X
j=i+1
aijx(k−1)j + ωbi, so that if A = L + D + U , then we have
(D + ωL)x(k) = [(1 − ω)D − ωU ] x(k−1)+ ωb or
x(k) = (D + ωL)−1[(1 − ω)D − ωU ] x(k−1)+ ω(D + ωL)−1b
≡ Tωx(k−1)+ cω.
46 / 99
logo
Example 25
The linear system Ax = b given by
4x1 + 3x2 = 24,
3x1 + 4x2 − x3 = 30,
− x2 + 4x3 = −24, has the solution [3, 4, −5]T.
Numerical results of Gauss-Seidel method with x(0)= [1, 1, 1]T:
k x1 x2 x3
0 1.0000000 1.0000000 1.0000000 1 5.2500000 3.8125000 -5.0468750 2 3.1406250 3.8828125 -5.0292969 3 3.0878906 3.9267578 -5.0183105 4 3.0549316 3.9542236 -5.0114441 5 3.0343323 3.9713898 -5.0071526 6 3.0214577 3.9821186 -5.0044703
logo
Numerical results of SOR method withω = 1.25and x(0)= [1, 1, 1]T:
k x1 x2 x3
0 1.0000000 1.0000000 1.0000000 1 6.3125000 3.5195313 -6.6501465 2 2.6223145 3.9585266 -4.6004238 3 3.1333027 4.0102646 -5.0966863 4 2.9570512 4.0074838 -4.9734897 5 3.0037211 4.0029250 -5.0057135 6 2.9963276 4.0009262 -4.9982822 7 3.0000498 4.0002586 -5.0003486
48 / 99
logo
Numerical results of SOR method withω = 1.6and x(0)= [1, 1, 1]T:
k x1 x2 x3
0 1.0000000 1.0000000 1.0000000 1 7.8000000 2.4400000 -9.2240000 2 1.9920000 4.4560000 -2.2832000 3 3.0576000 4.7440000 -6.3324800 4 2.0726400 4.1334400 -4.1471360 5 3.3962880 3.7855360 -5.5975040 6 3.0195840 3.8661760 -4.6950272 7 3.1488384 4.0236774 -5.1735127
logo Matlab code of SOR
clear all; delete rslt.dat; diary rslt.dat; diary on;
n = 3; xold = zeros(n,1); xnew = zeros(n,1); A = zeros(n,n); DL = zeros(n,n); DU = zeros(n,n);
A(1,1)=4; A(1,2)=3; A(2,1)=3; A(2,2)=4; A(2,3)=-1; A(3,2)=-1; A(3,3)=4;
b(1,1)=24; b(2,1)=30; b(3,1)=-24; omega=1.25;
for ii = 1:n DL(ii,ii) = A(ii,ii);
for jj = 1:ii-1
DL(ii,jj) = omega * A(ii,jj);
end
DU(ii,ii) = (1-omega)*A(ii,ii);
for jj = ii+1:n
DU(ii,jj) = - omega * A(ii,jj);
end end
c = omega * (DL \ b); xnew = DL \ ( DU * xold ) + c;
k = 0; fprintf(’ k x1 x2 x3 \n’);
while ( k <= 100 & norm(xnew-xold) > 1.0d-14 ) xold = xnew; k = k + 1; xnew = DL \ ( DU * xold ) + c;
fprintf(’%3.0f ’,k);
for jj = 1:n
fprintf(’%5.4f ’,xold(jj));
end fprintf(’\n’);
end diary off
50 / 99
logo
Theorem 26 (Kahan)
If aii6= 0, for each i = 1, 2, . . . , n, then ρ(Tω) ≥ |ω − 1|. This implies that the SOR method can converge only if 0 < ω < 2.
Theorem 27 (Ostrowski-Reich)
If A ispositive definiteand the relaxation parameter ω satisfying 0 < ω < 2, then theSORiterationconvergesforanyinitial vector x(0).
Theorem 28
If A ispositive definiteandtridiagonal, thenρ(TG) = [ρ(TJ)]2< 1and theoptimalchoice of ω for the SOR iteration is
ω = 2
1 + q
1 − [ρ(TJ)]2 .
With this choice of ω,ρ(Tω) = ω − 1.
logo
Example 29 The matrix
A =
4 3 0
3 4 −1
0 −1 4
,
given in previous example, is positive definite and tridiagonal.
Since
TJ = −D−1(L + U ) =
1
4 0 0
0 14 0 0 0 14
0 −3 0
−3 0 1
0 1 0
=
0 −0.75 0
−0.75 0 0.25
0 0.25 0
,
52 / 99
logo
we have
TJ − λI =
−λ −0.75 0
−0.75 −λ 0.25
0 0.25 −λ
, so
det(TJ − λI) = −λ(λ2− 0.625).
Thus,
ρ(TJ) =√ 0.625 and
ω = 2
1 +p1 − [ρ(TJ)]2 = 2 1 +√
1 − 0.625 ≈ 1.24.
This explains the rapid convergence obtained in previous example when using ω = 0.125
logo
Symmetric Successive Over Relaxation (SSOR) Method
Let A be symmetric and A = D + L + LT. The idea is in fact to implement the SOR formulationtwice,one forwardandone backward, at each iteration. That is, SSOR method defines
(D + ωL)x(k−12) = (1 − ω)D − ωLT x(k−1)+ ωb, (2) (D + ωLT)x(k) = [(1 − ω)D − ωL] x(k−12)+ ωb. (3) Define
Mω: = D + ωL,
Nω: = (1 − ω)D − ωLT. Then from the iterations (2) and (3), it follows that
x(k) = Mω−TNωTMω−1Nω x(k−1)+ ω Mω−TNωTMω−1+ Mω−T b
≡ T (ω)x(k−1)+ M (ω)−1b.
54 / 99
logo
But
((1 − ω)D − ωL) (D + ωL)−1+ I
= (−ωL − D − ωD + 2D)(D + ωL)−1+ I
= −I + (2 − ω)D(D + ωL)−1+ I
= (2 − ω)D(D + ωL)−1. Thus
M (ω)−1 = ω D + ωLT−1
(2 − ω)D(D + ωL)−1, then the splitting matrix is
M (ω) = 1
ω(2 − ω)(D + ωL)D−1 D + ωLT . The iteration matrix is
T (ω) = (D + ωLT)−1[(1 − ω)D − ωL] (D + ωL)−1(1 − ω)D − ωLT .
logo
Exercise Page 467: 2, 8
56 / 99
logo
Error bounds and iterative refinement
Example 30
The linear system Ax = b given by
1 2
1.0001 2
x1
x2
=
3
3.0001
has the unique solution x = [1, 1]T.
The poor approximation ˜x = [3, 0]T has the residual vector r = b − A˜x =
3
3.0001
−
1 2
1.0001 2
3 0
=
0
−0.0002
, so krk∞= 0.0002. Although the norm of the residual vector is small, the approximation ˜x = [3, 0]T is obviously quite poor; in fact, kx − ˜xk∞= 2.
logo
The solution of above example represents the intersection of the lines
`1: x1+ 2x2 = 3 and `2 : 1.0001x1+ 2x2 = 3.0001.
`1 and `2 are nearly parallel. The point (3, 0) lies on `1 which implies that (3, 0) also lies close to `2, even though it differs significantly from the intersection point (1, 1).
58 / 99
logo
Theorem 31
Suppose that ˜xis an approximate solution of Ax = b, A is nonsingular matrix and r = b − A˜x. Then
kx − ˜xk ≤ krk · kA−1k and if x 6= 0 and b 6= 0,
kx − ˜xk
kxk ≤ kAk · kA−1kkrk kbk. Proof: Since
r = b − A˜x = Ax − A˜x = A(x − ˜x) and A is nonsingular, we have
kx − ˜xk = kA−1rk ≤ kA−1k · krk. (4) Moreover, since b = Ax, we have
kbk ≤ kAk · kxk.
logo
It implies that
1
kxk ≤ kAk
kbk. (5)
Combining Equations (4) and (5), we have kx − ˜xk
kxk ≤ kAk · kA−1k kbk krk.
Definition 32 (Condition number)
The condition number of nonsingular matrix A is κ(A) = kAk · kA−1k.
For any nonsingular matrix A,
1 = kIk = kA · A−1k ≤ kAk · kA−1k = κ(A).
60 / 99
logo
Definition 33
A matrix A iswell-conditionedifκ(A)is close to1, and is ill-conditionedwhen κ(A) is significantly greater than 1.
In previous example, A =
1 2
1.0001 2
. Since
A−1 =
−10000 10000 5000.5 −5000
, we have
κ(A) = kAk∞· kA−1k∞= 3.0001 × 20000 = 60002 1.
logo
How to estimate the effective condition number in t-digit arithmetic without having to invert the matrix A?
If the approximate solution ˜xof Ax = b is being determined using t-digit arithmetic and Gaussian elimination, then
krk = kb − A˜xk ≈ 10−tkAk · k˜xk.
All the arithmetic operations in Gaussian elimination technique are performed using t-digit arithmetic, but the residual vector r are done in double-precision (i.e., 2t-digit) arithmetic.
Use the Gaussian elimination method which has already been calculated to solve
Ay = r.
Let ˜ybe the approximate solution.
62 / 99
logo
Then
˜
y ≈ A−1r = A−1(b − A˜x) = x − ˜x and
x ≈ ˜x + ˜y.
Moreover,
k˜yk ≈ kx − ˜xk = kA−1rk
≤ kA−1k · krk ≈ kA−1k(10−tkAk · k˜xk) = 10−tk˜xkκ(A).
It implies that
κ(A) ≈ k˜yk k˜xk10t. Iterative refinement
In general, ˜x + ˜yis a more accurate approximation to the solution of Ax = b than ˜x.
logo
Algorithm 3 (Iterative refinement)
Given tolerance T OL, maximum number of iteration M , number of digits of precision t.
Solve Ax = b by using Gaussian elimination in t-digit arithmetic.
Set k = 1 while ( k ≤ M )
Compute r = b − Ax in 2t-digit arithmetic.
Solve Ay = r by using Gaussian elimination in t-digit arithmetic.
If kyk∞< T OL, then stop.
Set k = k + 1 and x = x + y.
End while
64 / 99
