師大
Iterative Methods for Solving Large Linear Systems (I)
Tsung-Ming Huang
Department of Mathematics National Taiwan Normal University
October 25, 2011
師大
Outline
1 Jacobi and Gauss-Seidel methods
2 Successive Over-Relaxation (SOR) Method
3 Symmetric Successive Over Relaxation
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 2 / 45
師大
General procedures for the construction of iterative methods
Given a linear system of nonsingular A
Ax = b, (1)
we consider the splitting of A
A = M − N (2)
with M nonsingular. Then (1) is equivalent to M x = N x + b, or x = M−1N x + M−1b ≡ T x + f.
This suggests an iterative process
xk+1= T xk+ f = M−1N xk+ M−1b, (3) where x0 is given. Then the solution x of (1) is determined by iteration.
師大
Remark 1
(a) Define εk= xk− x. Then
εk+1= xk+1− x = M−1N xk+ M−1b − M−1N x − M−1b
= (M−1N )εk = (M−1N )kε0
which implies that ρ(M−1N ) < 1 if and only if {εk} → 0.
(b) Let rk= b − Axk. Then,
xk+1 = M−1N xk+ M−1b
= M−1(M − A)xk+ M−1b
= xk+ M−1(b − Axk)
= xk+ zk where M zk= rk.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 4 / 45
師大
Example 1
We consider the standard splitting of A
A = D − L − R, (4)
where A = [aij]ni,j=1, D = diag(a11, a22, · · · , ann),
−L =
0 0
a21 0
... . .. . ..
an1 · · · an,n−1 0
,
−R =
0 a12 · · · a1n 0 . .. ...
. .. an−1,n
0 0
.
師大
For ai,i6= 0, i = 1, . . . , n, D is nonsingular. If we choose M = D and N = L + R
in (2), we then obtain the Jacobi Method (Total-step Method):
xk+1= D−1(L + R)xk+ D−1b or in formula
xk+1,j = 1
ajj(−X
i6=j
ajixk,i+ bj), j = 1, . . . , n, k = 0, 1, . . . .
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 6 / 45
師大
Example 2
If D − L is nonsingular in (4), then we choose M = D − L, N = R
as in (2) are possible and yields the so-called Gauss-Seidel Method (Single-Step Method):
xk+1= (D − L)−1Rxk+ (D − L)−1b or in formula
xk+1,j = 1
ajj(−X
i<j
ajixk+1,i−X
i>j
ajixk,i+bj), j = 1, . . . , n, k = 1, 2, . . . .
- Total-Step Method = TSM = Jacobi method.
- Single-Step Method = SSM = Gauss-Seidel method.
師大
We consider the following points on Examples 1 and 2:
(i) flops counts per iteration step.
(ii) Convergence speed.
Let k · k be a vector norm, and kT k be the corresponding operator norm.
Then kεmk
kε0k = kTmε0k
kε0k ≤ kTmk. (5)
Here kTmkm1 is a measure for the average of reduction of error εm per iteration step. We call
Rm(T ) = − ln(kTmkm1) = −1
mln(kTmk) (6)
the average of convergence rate for m iterations.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 8 / 45
師大
The larger is Rm(T ), so the better is convergence rate. Let σ = (kεmk/kε0k)m1. From (5) and (6) we get
σ ≤ kTmkm1 ≤ e−Rm(T ), or
σ1/Rm(T )≤ 1 e.
That is, after 1/Rm(T ) steps in average the error is reduced by a factor of 1/e. Since Rm(T ) is not easy to determine, we consider m → ∞. Since
m→∞lim kTmkm1 = ρ(T ), it follows
R∞(T ) = lim
m→∞Rm(T ) = − ln ρ(T ).
R∞ is called the asymptotic convergence rate. It holds always Rm(T ) ≤ R∞(T ).
師大
Example 3
Consider the Dirichlet boundary-value problem (Model problem):
−∆u ≡ −uxx− uyy= f (x, y), 0 < x, y < 1, (7) u(x, y) = 0 (x, y) ∈ ∂Ω,
for the unit square Ω := {x, y|0 < x, y < 1} ⊆ R2 with boundary ∂Ω.
To solve (7) by means of a difference methods, one replaces the differential operator by a difference operator. Let
Ωh := {(xi, yi)|i, j = 1, . . . , N + 1},
∂Ωh := {(xi, 0), (xi, 1), (0, yj), (1, yj)|i, j = 0, 1, . . . , N + 1}, where
xi= ih, yj = jh, i, j = 0, 1, . . . , N + 1, h := N +11 , N ≥ 1, an integer.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 10 / 45
師大
The differential operator −uxx− uyy can be replaced for all (xi, yi) ∈ Ωh by the difference operator:
4ui,j− ui−1,j − ui+1,j − ui,j−1− ui,j+1
h2 (8)
up to an error τi,j. For small h one can expect that the solution zi,j, for i, j = 1, . . . , N , of the linear system
4zi,j− zi−1,j− zi+1,j − zi,j−1− zi,j+1= h2fi,j, i, j = 1, . . . , N, (9) z0,j = zN +1,j = zi,0= zi,N +1 = 0, i, j = 0, 1, . . . , N + 1,
obtained from (8) by omitting the error τi,j, agrees approximately with the ui,j. Let
z = [z1,1, z2,1, · · · , zN,1, z1,2, · · · , zN,2, · · · , z1,N, · · · , zN,N]T and
b = h2[f1,1, · · · , fN,1, f1,2, · · · , fN,2, · · · , f1,N, · · · , fN,N]T.
師大
Then (9) is equivalent to a linear system Az = b with the N2× N2matrix.
4 −1 −1
−1 ...
...
...
... −1
...
−1 4 −1
−1 4 −1
...
... −1
... ...
...
.. .
.. . −1
.. .
−1 −1 4
.. . ...
... −1
...
...
...
...
...
...
.. .
..
. −1
−1 4 −1
..
. −1
.. .
.. . ..
.
.. . −1
−1 −1 4
(11)
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 12 / 45
師大
Let A = D − L − R. The matrix J = D−1(L + R) belongs to the Jacobi method (TSM). The N2 eigenvalues and eigenvectors of J can be determined explicitly. We can verify at once, by substitution, that N2 vectors z(k,l), k, l = 1, . . . , N with components
zi,j(k,l) := sin kπi
N + 1sin lπj
N + 1, 1 ≤ i, j ≤ N, satisfy
J z(k,l) = λ(k,l)z(k,l) with
λ(k,l) := 1
2(cos kπ
N + 1 + cos lπ
N + 1), 1 ≤ k, l ≤ N.
師大
J thus has eigenvalues λ(k,l), 1 ≤ k, l ≤ N . Then we have ρ(J ) = λ1,1 = cos π
N + 1 = 1 −π2h2
2 + O(h4) (12)
and
R∞(J ) = − ln(1 − π2h2
2 + O(h4)) = π2h2
2 + O(h4).
These show that (i) TSM converges;
(ii) Diminution of h will not only enlarge the flop counts per step, but also the convergence speed will drastically make smaller.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 14 / 45
師大
Some theorems and definitions
ρ(T ): A measure of quality for convergence.
Definition 4
A real m × n-matrix A = (aik) is called nonnegative (positive), denoted by A ≥ 0 (A > 0), if aik ≥ 0 (> 0), i = 1, . . . , m, k = 1, . . . , n.
Definition 5
An m × n-matrix A is called reducible, if there is a subset
I ⊂ ˜N ≡ {1, 2, . . . , n}, I 6= φ, I 6= ˜N such that i ∈ I, j 6∈ I ⇒ aij = 0.
A is not reducible ⇔ A is irreducible.
師大
Remark 2
G(A) is the directed graph associated with the matrix A. If A is an n × n-matrix, then G(A) consists of n vertices P1, · · · , Pn and there is an (oriented) arc Pi → Pj in G(A) precisely if aij 6= 0.
It is easily shown that A is irreducible if and only if the graph G(A) is connected in the sense that for each pair of vertices (Pi, Pj) in G(A) there is an oriented path from Pi to Pj. i.e., if i 6= j, there is a sequence of indices i = i1, i2, · · · , is= j such that (ai1,i2 · · · ais−1,is) 6= 0.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 16 / 45
師大
Theorem 6 (Perron-Frobenius) Let A ≥ 0 irreducible. Then
(i) ρ = ρ(A) is a simple eigenvalue;
(ii) There is a positive eigenvector z associated to ρ, i.e., Az = ρz, z > 0;
(iii) If Ax = λx, x ≥ 0, then λ = ρ, x = αz, α > 0;
(iv) A ≤ B, A 6= B =⇒ ρ(A) < ρ(B).
師大
Theorem 7
Let A ≥ 0, x > 0. Define the quotients:
qi(x) ≡ (Ax)i
xi
= 1 xi
n
X
k=1
aikxk, for i = 1, . . . , n.
Then
1≤i≤nmin qi(x) ≤ ρ(A) ≤ max
1≤i≤nqi(x). (13)
If A is irreducible, then it holds additionally, either
q1 = q2= · · · = qn (then x = µz, qi= ρ(A)) or
1≤i≤nmin qi(x) < ρ(A) < max
1≤i≤nqi(x). (14)
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 18 / 45
師大
Theorem 8
The statements in Theorem 7 can be formulated as: Let A ≥ 0, x > 0.
(13) corresponds:
Ax ≤ µx ⇒ ρ ≤ µ, Ax ≥ νx ⇒ ν ≤ ρ.
(15)
Let A ≥ 0, irreducible, x > 0. (14) corresponds :
Ax ≤ µx, Ax 6= µx ⇒ ρ < µ,
Ax ≥ νx, Ax 6= νx ⇒ ν < ρ. (16)
Definition 9
A real matrix B is called an M -matrix if bij ≤ 0, i 6= j and B−1 exists with B−1≥ 0.
師大
Sufficient conditions for convergence of TSM and SSM
Theorem 10
Let B be a real matrix with bij ≤ 0 for i 6= j. Then the following statements are equivalent.
(i) B is an M −matrix.
(ii) There exists a vector v > 0 so that Bv > 0.
(iii) B has a decomposition B = sI − C with C ≥ 0 and ρ(C) < s.
(iv) For each decomposition B = D − C with D = diag (di) and C ≥ 0, it holds: di> 0, i = 1, 2, . . . , n, and ρ(D−1C) < 1.
(v) There is a decomposition B = D − C, with D = diag(di) and C ≥ 0 it holds: di> 0, i = 1, 2, . . . , n and ρ(D−1C) < 1.
Further, if B is irreducible, then (vi) is equivalent to (i)-(v).
(vi) There exists a vector v > 0 so that Bv ≥ 0, 6= 0.
Proof
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 20 / 45
師大
Lemma 11
Let A be an arbitrary complex matrix and define |A| = [|aij|]. If |A| ≤ C, then ρ(A) ≤ ρ(C). Especially ρ(A) ≤ ρ(|A|).
Proof
Theorem 12
Let A be an arbitrary complex matrix. It satisfies either (Strong Row Sum Criterion):
X
j6=i
|aij| < |aii|, i = 1, . . . , n. (17)
or (Weak Row Sum Criterion):
X
j6=i
|aij| ≤ |aii|, i = 1, . . . , n,
< |ai0i0|, at least one i0, (18)
for A irreducible. Then TSM(Jacobi) and SSM(GS) both are convergent.
Proof
師大
Relaxation Methods (Successive Over-Relaxation (SOR) Method) Consider the parametrized linear system ωAx = ωb and consider the splitting
ωA = ωD − ωL − ωR + D − D
= (D − ωL) − ((1 − ω)D + ωR) ≡ M − N.
From (3) we have the iteration
xk+1 = (D − ωL)−1((1 − ω)D + ωR) xk+ ω(D − ωL)−1b. (19) From Remark 1 (b) the iteration (19) is equivalent to
xk+1 = xk+ ωzk where
(D − ωL)zk= rk≡ b − Axk.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 22 / 45
師大
Define
Lω := (D − ωL)−1((1 − ω)D + ωR) . We may assume D = I, i.e.,
Lω:= (I − ωL)−1((1 − ω)I + ωR) .
Otherwise, we can let ˜A = D−1A, ˜L = D−1L, ˜R = D−1R. Then it holds that
A = I − ˜˜ L − ˜R.
ω < 1: under relaxation
ω = 1: single-step method (GS) ω > 1: over relaxation.
師大
We now try to choose an ω such that ρ(Lω) is small as possible. But this is only under some special assumptions possible. we first list a few
qualitative results about ρ(Lω).
Theorem 13
Let A = D − L − L∗ be Hermitian and positive definite. Then the relaxation method is convergent for 0 < ω < 2.
Theorem 14
Let A be Hermitian and nonsingular with positive diagonal. If SSM converges, then A is positive definite.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 24 / 45
師大
Determination of the Optimal Parameter ω for 2-consistly Ordered Matrices
For an important class of matrices the more qualitative assertions of Theorems 13 and 14 can be considerably sharpened. This is the class of consistly ordered matrices. The optimal parameter ωb with
ρ(Lωb) = min
ω ρ(Lω) can be determined. We consider A = I − L − R.
Definition 15
A is called 2-consistly ordered, if the eigenvalues of αL + α−1R are independent of α.
師大
If A is 2-consistly ordered, then L + R and −(L + R) (α = −1) has the same eigenvalues. The nonzero eigenvalues of L + R appear in pairs.
Hence
det(λI − L − R) = λm
r
Y
i=1
(λ2− µ2i), n = 2r + m (m = 0, possible). (20)
Theorem 16
Let A be 2-consistly ordered, aii= 1, ω 6= 0. Then hold:
(i) If λ 6= 0 is an eigenvalue of Lω and µ satisfies the equation
(λ + ω − 1)2= λµ2ω2, (21) then µ is an eigenvalue of L + R (so is −µ).
(ii) If µ is an eigenvalue of L + R and λ satisfies the equation (21), then λ is an eigenvalue of Lω.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 26 / 45
師大
Remark 3
If ω = 1, then λ = µ2, and ρ((I − L)−1R) = (ρ(L + R))2. Proof: We first prove the identity
det(λI − sL − rR) = det(λI −√
sr(L + R)). (22)
Since both sides are polynomials of the form λn+ · · · and sL + rR =√
sr(r s
rL +r r
sR) =√
sr(αL + α−1R), if sr 6= 0, then sL + rR and √
sr(L + R) have the same eigenvalues. It is obviously also for the case sr = 0. The both polynomials in (22) have the same roots, so they are identical.
師大
For
det(I − ωL) det(λI − Lω) = det(λ(I − ωL) − (1 − ω)I − ωR)
= det((λ + ω − 1)I − ωλL − ωR) = Φ(λ) and det(I − ωL) 6= 0, λ is an eigenvalue of Lω if and only if Φ(λ) = 0.
From (22) follows
Φ(λ) = det((λ + ω − 1)I − ω
√
λ(L + R)) and that is (from (20))
Φ(λ) = (λ + ω − 1)m
r
Y
i=1
((λ + ω − 1)2− λµ2iω2), (23)
where µi is an eigenvalue of L + R. Therefore, if µ is an eigenvalue of (L + R) and λ satisfies (21), so is Φ(λ) = 0, then λ is eigenvalue of Lω. This shows (ii).
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 28 / 45
師大
Now if λ 6= 0 an eigenvalue of Lω, then one factor in (23) must be zero.
Let µ satisfy (21). Then
(i) µ 6= 0: From (21) follows λ + ω − 1 6= 0, so
(λ + ω − 1)2 = λω2µ2i, for one i (from (23)),
= λω2µ2, (from (21)).
This shows that µ = ±µi, so µ is an eigenvalue of L + R.
(ii) µ = 0: We have λ + ω − 1 = 0 and 0 = Φ(λ) = det((λ + ω − 1)I − ω√
λ(L + R)) = det(−ω√
λ(L + R)), i.e., L + R is singular, so µ = 0 is eigenvalue of L + R.
師大
Theorem 17
Let A = I − L − R be 2-consistly ordered. If L + R has only real eigenvalues and satisfies ρ(L + R) < 1, then it holds
ρ(Lωb) = ωb− 1 < ρ(Lω), for ω 6= ωb, where
ωb = 2
1 +p1 − ρ2(L + R) (solve ωb in (21)).
0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2
0.9975 0.998 0.9985 0.999 0.9995 1
ω
spectral radius
Figure: figure of ρ(Lωb)
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 30 / 45
師大
Theorem 18 One has in general,
ρ(Lω) =
( ω − 1, for ωb≤ ω ≤ 2 1 − ω + 12ω2µ2+ ωµ
q
1 − ω + 14ω2µ2, for 0 < ω ≤ ωb (24)
師大
Remark: We first prove the following Theorem proposed by Kahan: For arbitrary matrices A it holds
ρ(Lω) ≥ |ω − 1|, for all ω. (25) Since det(I − ωL) = 1 for all ω, the characteristic polynomial Φ(λ) of Lω is
Φ(λ) = det(λI − Lω) = det((I − ωL)(λI − Lω))
= det((λ + ω − 1)I − ωλL − ωR).
For
n
Q
i=1
λi(Lω) = Φ(0) = det((ω − 1)I − ωR) = (ω − 1)n, it follows immediately that
ρ(Lω) = max
i |λi(Lω)| ≥ |ω − 1|.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 32 / 45
師大
Proof of Theorem: By assumption the eigenvalues µi of L + R are real and −ρ(L + R) ≤ µi≤ ρ(L + R) < 1. For a fixed ω ∈ (0, 2) (by (25) in the Remark it suffices to consider the interval (0,2)) and for each µi there are two eigenvalues λ(1)i (ω, µi) and λ(2)i (ω, µi) of Lω, which are obtained by solving the quadratic equation (21) in λ.
Geometrically, λ(1)i (ω) and λ(2)i (ω) are obtained as abscissae of the points of intersection of
the straight line gω(λ) = λ + ω − 1 ω and
the parabola mi(λ) := ±√ λµi
(see Figure 2). The line gω(λ) has the slope 1/ω and passes through the point (1,1). If gω(λ) ∩ mi(λ) = φ, then λ(1)i (ω) and λ(2)i (ω) are conjugate complex with modulus |ω − 1| (from (21)).
師大
Evidently
ρ(Lω) = max
i (|λ(1)i (ω)|, |λ(2)i (ω)|) = max(|λ(1)(ω)|, |λ(2)(ω)|), where λ(1)(ω), λ(2)(ω) being obtained by intersecting gω(λ) with m(λ) := ±√
λµ, with µ = ρ(L + R) = maxi|µi|. By solving (21) with µ = ρ(L + R) for λ, one verifies (24) immediately, and thus also the remaining assertions of the theorem.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 34 / 45
師大
1 1
( )
gω λ ( )
b gω λ
( )
m λ
( ) i
m λ
λ
(1)
i
λ
( 2 )
i
λ θ
Figure: Geometrical view of λ(1)i (ω) and λ(2)i (ω).
師大
Application to Finite Difference Methods: Model Problem
We consider the Dirichlet boundary-value problem (Model problem) as in Example 3. We shall solve a linear system Az = b of the N2× N2 matrix A as in (11).
To Jacobi method: The iterative matrix is J = L + R = 1
4(4I − A).
It is easily seen that A is 2-consistly ordered (Exercise!).
To Gauss-Seidel method: The iterative matrix is H = (I − L)−1R.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 36 / 45
師大
From the Remark of Theorem 16 and (12) follows that ρ(H) = ρ(J )2 = cos2 π
N + 1.
According to Theorem 17 the optimal relaxation parameter ωb and ρ(Lωb) are given by
ωb = 2
1 +q
1 − cos2 N +1π
= 2
1 + sinN +1π and
ρ(Lωb) = cos2N +1π (1 + sinN +1π )2.
The number k = k(N ) with ρ(J )k= ρ(Lωb) indicates that the k steps of Jacobi method produce the same reduction as one step of the optimal relaxation method. Clearly,
k = ln ρ(Lωb)/ ln ρ(J ).
師大
Now for small z one has ln(1 + z) = z − z2/2 + O(z3) and for large N we have
cos
π
N + 1
= 1 − π2
2(N + 1)2 + O( 1 N4).
Thus that
ln ρ(J ) = π2
2(N + 1)2 + O( 1 N4).
Similarly,
ln ρ(Lωb) = 2[ln ρ(J ) − ln(1 + sin π N + 1)]
= 2[− π2
2(N + 1)2 − π
N + 1 + π2
2(N + 1)2 + O( 1 N3)]
= − 2π
N + 1+ O( 1
N3) (for large N ).
and
k = k(N ) ≈ 4(N + 1)
π .
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 38 / 45
師大
The optimal relaxation method is more than N times as fast as the Jacobi method. The quantities
RJ := − ln 10
ln ρ(J ) ≈ 0.467(N + 1)2. (26) RH := 1
2RJ ≈ 0.234(N + 1)2 (27)
RLωb := − ln 10
ln ρ(Lωb) ≈ 0.367(N + 1) (28) indicate the number of iterations required in the Jacobi, the Gauss-Seidel method, and the optimal relaxation method, respectively, in order to reduce the error by a factor of 1/10.
師大
SSOR (Symmetric Successive Over Relaxation):
A is symmetric and A = D − L − LT. Let
Mω: = D − ωL,
Nω: = (1 − ω)D + ωLT, and
MωT = D − ωLT, NωT = (1 − ω)D + ωL.
Then from the iterations
Mωxi+1/2 = Nωxi+ ωb, MωTxi+1 = NωTxi+1/2+ ωb, follows that
xi+1 = Mω−TNωTMω−1Nω xi+ ˜b
≡ Gxi+ ω Mω−TNωTMω−1+ Mω−T b
≡ Gxi+ M (ω)−1b.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 40 / 45
師大
It holds that
((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
M (ω) = 1
ω(2 − ω)(D − ωL)D−1 D − ωLT
(29)
≈ (D − L)D−1 D − LT , (ω = 1).
師大
Appendix
Proof of Theorem
(1) =⇒ (2): Let e = (1, · · · , 1)T. Since B−1≥ 0 is nonsingular it follows ν = B−1e > 0 and Bν = B(B−1e) = e > 0.
(2) =⇒ (3): Let s > max(bii). It follows B = sI − C with C ≥ 0. There exists a ν > 0 with Bν = sν − Cν (via (2)), also sν > Cν. From the statement (15) in Theorem 8 follows ρ(C) < s.
(3) =⇒ (1): B = sI − C = s(I −1sC). For ρ(1sC) < 1 and from Theorem 2.6 (I −1sC)−1 follows that there exists a series expansion
∞
P
ν=0
(1sC)k. Since the terms in sum are nonnegtive, we get B−1 = 1s(I − 1sC)−1 ≥ 0.
(2) =⇒ (4): From Bν = Dν − Cν > 0 follows Dν > Cν ≥ 0 and di > 0, for i = 1, 2, · · · , n. Hence D−1 ≥ 0 and ν > D−1Cν ≥ 0. From (15) follows that ρ(D−1C) < 1.
(4) =⇒ (5): Trivial.
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 42 / 45
師大
(5) =⇒ (1): Since ρ(D−1C) < 1, it follows from Theorem 2.6 that (I − D−1C)−1 exists and equals to
∞
P
k=0
(D−1C)k. Since the terms in sum are nonnegative, we have (I − D−1C)−1 is nonnegative and
B−1= (I − D−1C)−1D−1 ≥ 0.
(2) =⇒ (6): Trivial.
(6) =⇒ (5): Consider the decomposition B = D − C, with di = bii. Let {I = i | di≤ 0}. From diνi−P
k6=icikνk≥ 0 follows cik = 0 for i ∈ I, and k 6= i. Since Bν ≥ 0, 6= 0 =⇒ I 6= {1, · · · , n}. But B is irreducible
=⇒ I = ∅ and di > 0. Hence for Dν >, 6= Cν also ν >, 6= D−1Cν and (16) show that ρ(D−1C) < 1.
return
師大
Proof of Lemma 11 There is a x 6= 0 with Ax = λx and |λ| = ρ(A).
Hence
ρ(A)|xi| = |
n
X
k=1
aikxk| ≤
n
X
k=1
|aik||xk| ≤
n
X
k=1
cik|xk|.
Thus,
ρ(A)|x| ≤ C|x|.
If |x| > 0, then from (15) we have ρ(A) ≤ ρ(C). Otherwise, let
I = {i | xi 6= 0} and CI be the matrix, which consists of the ith row and ith column of C with i ∈ I. Then we have ρ(A)|xI| ≤ CI|xI|. Here |xI| consists of ith component of |x| with i ∈ I. Then from |xI| > 0 and (15) follows ρ(A) ≤ ρ(CI). We now fill CI with zero up to an n × n matrix ˜CI. Then ˜CI ≤ C. Thus, ρ(CI) = ρ( ˜CI) ≤ ρ(C) (by Theorem ??(3)).
return
T.M. Huang (NTNU) Iterative Methods for LS October 25, 2011 44 / 45
師大
Proof of Theorem ?? Let A = D − L − R. From (17) and (18) D must be nonnsingular and then as in Remark 9.3 we can w.l.o.g. assume that D = I. Now let B = I − |L| − |R|. Then (17) can be written as Be > 0.
From Theorem 10(2) and (1) follows that B is an M -matrix.
(18) can be written as Be ≥ 0, Be 6= 0. Since A is irreducible, also B, from Theorem 10 (6) and (1) follows that B is an M -matrix.
Especially, from theorem 10(1), (4) and Theorem ?? follows that ρ(|L| + |R|) < 1 and ρ((I − |L|)−1|R|) < 1.
Now Lemma 11 shows that
ρ(L + R) ≤ ρ(|L| + |R|) < 1.
So TSM is convergent. Similarly,
ρ((I − L)−1R) = ρ(R + LR + · · · + Ln−1R)
≤ ρ(|R| + |L||R| + · · · + |L|n−1|R|)
= ρ((I − |L|)−1|R|) < 1.
So SSM is convergent.
return