Iterative Methods for Solving Large Linear Systems (I)

(1)

師大

Iterative Methods for Solving Large Linear Systems (I)

Tsung-Ming Huang

Department of Mathematics National Taiwan Normal University

October 25, 2011

(2)

師大

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

(3)

師大

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⁻¹N x + M⁻¹b ≡ T x + f.

This suggests an iterative process

x_k+1= T x_k+ f = M⁻¹N x_k+ M⁻¹b, (3) where x₀ is given. Then the solution x of (1) is determined by iteration.

(4)

師大

Remark 1

(a) Define ε_k= x_k− x. Then

ε_k+1= x_k+1− x = M⁻¹N x_k+ M⁻¹b − M⁻¹N x − M⁻¹b

= (M⁻¹N )ε_k = (M⁻¹N )^kε₀

which implies that ρ(M⁻¹N ) < 1 if and only if {ε_k} → 0.

(b) Let rk= b − Axk. Then,

xk+1 = M⁻¹N xk+ M⁻¹b

= M⁻¹(M − A)x_k+ M⁻¹b

= x_k+ M⁻¹(b − Ax_k)

= x_k+ z_k where M zk= rk.

(5)

師大

Example 1

We consider the standard splitting of A

A = D − L − R, (4)

where A = [a_ij]ⁿ_i,j=1, D = diag(a₁₁, a₂₂, · · · , a_nn),

−L =







0 0

a21 0

... . .. . ..

an1 · · · an,n−1 0





 ,

−R =







0 a₁₂ · · · a_1n 0 . .. ...

. .. a_n−1,n

0 0





 .

(6)

師大

For a_i,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⁻¹(L + R)xk+ D⁻¹b or in formula

x_k+1,j = 1

a_jj(−X

i6=j

a_jix_k,i+ b_j), j = 1, . . . , n, k = 0, 1, . . . .

(7)

師大

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)⁻¹Rxk+ (D − L)⁻¹b or in formula

x_k+1,j = 1

a_jj(−X

i<j

ajix_k+1,i−X

i>j

ajix_k,i+bj), j = 1, . . . , n, k = 1, 2, . . . .

- Total-Step Method = TSM = Jacobi method.

- Single-Step Method = SSM = Gauss-Seidel method.

(8)

師大

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ε₀k = kT^mε₀k

kε₀k ≤ kT^mk. (5)

Here kT^mk^m¹ is a measure for the average of reduction of error εm per iteration step. We call

Rm(T ) = − ln(kT^mk^m¹) = −1

mln(kT^mk) (6)

the average of convergence rate for m iterations.

(9)

師大

The larger is R_m(T ), so the better is convergence rate. Let σ = (kε_mk/kε₀k)^m¹. From (5) and (6) we get

σ ≤ kT^mk^m¹ ≤ e^−R^m^{(T )}, or

σ^1/R^m^{(T )}≤ 1 e.

That is, after 1/R_m(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 kT^mk^m¹ = ρ(T ), it follows

R∞(T ) = lim

m→∞R_m(T ) = − ln ρ(T ).

R∞ is called the asymptotic convergence rate. It holds always R_m(T ) ≤ R∞(T ).

(10)

師大

Example 3

Consider the Dirichlet boundary-value problem (Model problem):

−∆u ≡ −u_xx− u_yy= f (x, y), 0 < x, y < 1, (7) u(x, y) = 0 (x, y) ∈ ∂Ω,

for the unit square Ω := {x, y|0 < x, y < 1} ⊆ R² with boundary ∂Ω.

To solve (7) by means of a difference methods, one replaces the differential operator by a difference operator. Let

Ω_h := {(x_i, yi)|i, j = 1, . . . , N + 1},

∂Ωh := {(x_i, 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 +1}¹ , N ≥ 1, an integer.

(11)

師大

The differential operator −u_xx− u_yy can be replaced for all (x_i, y_i) ∈ Ω_h by the difference operator:

4ui,j− u_i−1,j − u_i+1,j − u_i,j−1− u_i,j+1

h² (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− z_i−1,j− z_i+1,j − z_i,j−1− z_i,j+1= h²fi,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 u_i,j. Let

z = [z_1,1, z_2,1, · · · , z_N,1, z_1,2, · · · , z_N,2, · · · , z_1,N, · · · , z_N,N]^T and

b = h²[f1,1, · · · , f_N,1, f1,2, · · · , f_N,2, · · · , f_1,N, · · · , f_N,N]^T.

(12)

師大

Then (9) is equivalent to a linear system Az = b with the N²× N²matrix.







4 −1 −1

−1 ...

...

... −1

...

−1 4 −1

...

... −1

... ...

...

.. .

.. . −1

.. .

−1 −1 4

.. . ...

... −1

...

.. .

..

. −1

−1 4 −1

..

. −1

.. .

.. . ..

.

.. . −1

−1 −1 4







(11)

(13)

師大

Let A = D − L − R. The matrix J = D⁻¹(L + R) belongs to the Jacobi method (TSM). The N² eigenvalues and eigenvectors of J can be determined explicitly. We can verify at once, by substitution, that N² vectors z^(k,l), k, l = 1, . . . , N with components

z_i,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.

(14)

師大

J thus has eigenvalues λ^(k,l), 1 ≤ k, l ≤ N . Then we have ρ(J ) = λ_1,1 = cos π

N + 1 = 1 −π²h²

2 + O(h⁴) (12)

and

R∞(J ) = − ln(1 − π²h²

2 + O(h⁴)) = π²h²

2 + O(h⁴).

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.

(15)

師大

Some theorems and definitions

ρ(T ): A measure of quality for convergence.

Definition 4

A real m × n-matrix A = (a_ik) is called nonnegative (positive), denoted by A ≥ 0 (A > 0), if a_ik ≥ 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 ⇒ a_ij = 0.

A is not reducible ⇔ A is irreducible.

(16)

師大

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 → P_j 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 = i₁, i₂, · · · , i_s= j such that (a_i₁_,i₂ · · · a_i_s−1_,i_s) 6= 0.

(17)

師大

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).

(18)

師大

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)

(19)

師大

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 b_ij ≤ 0, i 6= j and B⁻¹ exists with B⁻¹≥ 0.

(20)

師大

Sufficient conditions for convergence of TSM and SSM

Theorem 10

Let B be a real matrix with b_ij ≤ 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⁻¹C) < 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⁻¹C) < 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

(21)

師大

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

|a_ij| < |a_ii|, i = 1, . . . , n. (17)

or (Weak Row Sum Criterion):

X

j6=i

|aij| ≤ |aii|, i = 1, . . . , n,

< |ai₀i₀|, at least one i0, (18)

for A irreducible. Then TSM(Jacobi) and SSM(GS) both are convergent.

Proof

(22)

師大

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

x_k+1 = (D − ωL)⁻¹((1 − ω)D + ωR) x_k+ ω(D − ωL)⁻¹b. (19) From Remark 1 (b) the iteration (19) is equivalent to

x_k+1 = x_k+ ωz_k where

(D − ωL)z_k= r_k≡ b − Ax_k.

(23)

師大

Define

Lω := (D − ωL)⁻¹((1 − ω)D + ωR) . We may assume D = I, i.e.,

Lω:= (I − ωL)⁻¹((1 − ω)I + ωR) .

Otherwise, we can let ˜A = D⁻¹A, ˜L = D⁻¹L, ˜R = D⁻¹R. Then it holds that

A = I − ˜˜ L − ˜R.

ω < 1: under relaxation

ω = 1: single-step method (GS) ω > 1: over relaxation.

(24)

師大

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.

(25)

師大

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 + α⁻¹R are independent of α.

(26)

師大

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

(λ²− µ²_i), n = 2r + m (m = 0, possible). (20)

Theorem 16

Let A be 2-consistly ordered, a_ii= 1, ω 6= 0. Then hold:

(i) If λ 6= 0 is an eigenvalue of L_ω and µ satisfies the equation

(λ + ω − 1)²= λµ²ω², (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ω.

(27)

師大

Remark 3

If ω = 1, then λ = µ², and ρ((I − L)⁻¹R) = (ρ(L + R))². Proof: We first prove the identity

det(λI − sL − rR) = det(λI −√

sr(L + R)). (22)

Since both sides are polynomials of the form λⁿ+ · · · and sL + rR =√

sr(r s

rL +r r

sR) =√

sr(αL + α⁻¹R), 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.

(28)

師大

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)²− λµ²_iω²), (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).

(29)

師大

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)² = λω²µ²_i, for one i (from (23)),

= λω²µ², (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.

(30)

師大

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 − ρ²(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)

(31)

師大

Theorem 18 One has in general,

ρ(Lω) =

( ω − 1, for ωb≤ ω ≤ 2 1 − ω + ¹₂ω²µ²+ ωµ

q

1 − ω + ¹₄ω²µ², for 0 < ω ≤ ω_b (24)

(32)

師大

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)ⁿ, it follows immediately that

ρ(Lω) = max

i |λ_i(Lω)| ≥ |ω − 1|.

(33)

師大

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 λ⁽¹⁾_i (ω, µ_i) and λ⁽²⁾_i (ω, µ_i) of L_ω, which are obtained by solving the quadratic equation (21) in λ.

Geometrically, λ⁽¹⁾_i (ω) and λ⁽²⁾_i (ω) are obtained as abscissae of the points of intersection of

the straight line g_ω(λ) = λ + ω − 1 ω and

the parabola m_i(λ) := ±√ λµ_i

(see Figure 2). The line g_ω(λ) has the slope 1/ω and passes through the point (1,1). If gω(λ) ∩ mi(λ) = φ, then λ⁽¹⁾_i (ω) and λ⁽²⁾_i (ω) are conjugate complex with modulus |ω − 1| (from (21)).

(34)

師大

Evidently

ρ(L_ω) = max

i (|λ⁽¹⁾_i (ω)|, |λ⁽²⁾_i (ω)|) = max(|λ⁽¹⁾(ω)|, |λ⁽²⁾(ω)|), where λ⁽¹⁾(ω), λ⁽²⁾(ω) 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.

(35)

師大

1 1

( )

gω λ ( )

b gω λ

( )

m λ

( ) i

m λ

λ

(1)

i

λ

( 2 )

i

λ θ

Figure: Geometrical view of λ⁽¹⁾_i (ω) and λ⁽²⁾_i (ω).

(36)

師大

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 N²× N² 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)⁻¹R.

(37)

師大

From the Remark of Theorem 16 and (12) follows that ρ(H) = ρ(J )² = cos² π

N + 1.

According to Theorem 17 the optimal relaxation parameter ω_b and ρ(L_ω_b) are given by

ω_b = 2

1 +q

1 − cos² _{N +1}^π

= 2

1 + sin_{N +1}^π and

ρ(Lω_b) = cos²_{N +1}^π (1 + sin_{N +1}^π )².

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 ).

(38)

師大

Now for small z one has ln(1 + z) = z − z²/2 + O(z³) and for large N we have

cos

π

N + 1

= 1 − π²

2(N + 1)² + O( 1 N⁴).

Thus that

ln ρ(J ) = π²

2(N + 1)² + O( 1 N⁴).

Similarly,

ln ρ(L_ω_b) = 2[ln ρ(J ) − ln(1 + sin π N + 1)]

= 2[− π²

2(N + 1)² − π

N + 1 + π²

2(N + 1)² + O( 1 N³)]

= − 2π

N + 1+ O( 1

N³) (for large N ).

and

k = k(N ) ≈ 4(N + 1)

π .

(39)

師大

The optimal relaxation method is more than N times as fast as the Jacobi method. The quantities

R_J := − ln 10

ln ρ(J ) ≈ 0.467(N + 1)². (26) R_H := 1

2R_J ≈ 0.234(N + 1)² (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.

(40)

師大

SSOR (Symmetric Successive Over Relaxation):

A is symmetric and A = D − L − L^T. Let

M_ω: = D − ωL,

N_ω: = (1 − ω)D + ωL^T, and

M_ω^T = D − ωL^T, N_ω^T = (1 − ω)D + ωL.

Then from the iterations

Mωx_i+1/2 = Nωxi+ ωb, M_ω^Txi+1 = N_ω^Tx_i+1/2+ ωb, follows that

x_i+1 = M_ω^−TN_ω^TM_ω⁻¹N_ω x_i+ ˜b

≡ Gx_i+ ω M_ω^−TN_ω^TM_ω⁻¹+ M_ω^−T b

≡ Gx_i+ M (ω)⁻¹b.

(41)

師大

It holds that

((1 − ω)D + ωL) (D − ωL)⁻¹+ I

= (ωL − D − ωD + 2D)(D − ωL)⁻¹+ I

= −I + (2 − ω)D(D − ωL)⁻¹+ I

= (2 − ω)D(D − ωL)⁻¹, Thus

M (ω)⁻¹ = ω D − ωL^T−1

(2 − ω)D(D − ωL)⁻¹, then

M (ω) = 1

ω(2 − ω)(D − ωL)D⁻¹ D − ωL^T

(29)

≈ (D − L)D⁻¹ D − L^T , (ω = 1).

(42)

師大

Appendix

Proof of Theorem

(1) =⇒ (2): Let e = (1, · · · , 1)^T. Since B⁻¹≥ 0 is nonsingular it follows ν = B⁻¹e > 0 and Bν = B(B⁻¹e) = 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 −¹_sC). For ρ(¹_sC) < 1 and from Theorem 2.6 (I −¹_sC)⁻¹ follows that there exists a series expansion

∞

P

ν=0

(¹_sC)^k. Since the terms in sum are nonnegtive, we get B⁻¹ = ¹_s(I − ¹_sC)⁻¹ ≥ 0.

(2) =⇒ (4): From Bν = Dν − Cν > 0 follows Dν > Cν ≥ 0 and d_i > 0, for i = 1, 2, · · · , n. Hence D⁻¹ ≥ 0 and ν > D⁻¹Cν ≥ 0. From (15) follows that ρ(D⁻¹C) < 1.

(4) =⇒ (5): Trivial.

(43)

師大

(5) =⇒ (1): Since ρ(D⁻¹C) < 1, it follows from Theorem 2.6 that (I − D⁻¹C)⁻¹ exists and equals to

∞

P

k=0

(D⁻¹C)^k. Since the terms in sum are nonnegative, we have (I − D⁻¹C)⁻¹ is nonnegative and

B⁻¹= (I − D⁻¹C)⁻¹D⁻¹ ≥ 0.

(2) =⇒ (6): Trivial.

(6) =⇒ (5): Consider the decomposition B = D − C, with d_i = b_ii. Let {I = i | d_i≤ 0}. From d_iνi−P

k6=icikνk≥ 0 follows c_ik = 0 for i ∈ I, and k 6= i. Since Bν ≥ 0, 6= 0 =⇒ I 6= {1, · · · , n}. But B is irreducible

=⇒ I = ∅ and d_i > 0. Hence for Dν >, 6= Cν also ν >, 6= D⁻¹Cν and (16) show that ρ(D⁻¹C) < 1.

return

(44)

師大

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

|a_ik||x_k| ≤

n

X

k=1

cik|x_k|.

Thus,

ρ(A)|x| ≤ C|x|.

If |x| > 0, then from (15) we have ρ(A) ≤ ρ(C). Otherwise, let

I = {i | x_i 6= 0} and C_I be the matrix, which consists of the ith row and ith column of C with i ∈ I. Then we have ρ(A)|xI| ≤ C_I|x_I|. Here |x_I| consists of ith component of |x| with i ∈ I. Then from |x_I| > 0 and (15) follows ρ(A) ≤ ρ(C_I). We now fill C_I with zero up to an n × n matrix ˜C_I. Then ˜CI ≤ C. Thus, ρ(C_I) = ρ( ˜CI) ≤ ρ(C) (by Theorem ??(3)).

return

(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|)⁻¹|R|) < 1.

Now Lemma 11 shows that

ρ(L + R) ≤ ρ(|L| + |R|) < 1.

So TSM is convergent. Similarly,

ρ((I − L)⁻¹R) = ρ(R + LR + · · · + Lⁿ⁻¹R)

≤ ρ(|R| + |L||R| + · · · + |L|ⁿ⁻¹|R|)

= ρ((I − |L|)⁻¹|R|) < 1.

So SSM is convergent.

return