White R.E. Computational Mathematics: Models, Methods, and Analysis with MATLAB and MPI

Подождите немного. Документ загружается.

78 CHAPTER 2. STEADY STATE DISCRETE MODELS

usion in 2D. Let x = x({> |) = temperature on a square.

0 =

i + (Nx

{

)

{

+ (Nx

)

and (2.4.3)

x = given on the boundary= (2.4.4)

The discrete models can be either viewed as discrete versions of the Fourier

heat law, or as ﬁnite di

erence approximations of the continuous models.

Discrete Models:

usion in 1D. Let x

approximate x(lk) with k = O@q.

0 =

i + (x

l+1

+ x

l1

)  2x

(2.4.5)

where

l = 1> ===> q  1 and  = N@k

and

> x

= given= (2.4.6)

usion in 2D. Let x

approximate x(lk> mk) with k = O@q = { = |.

0 =

i + (x

l+1>m

+ x

l1>m

)  2x

l>m

(x

l>m+1

+ x

l>m1

)  2x

l>m

(2.4.7)

where

l> m = 1> ===> q  1 and  = N@k

and

0>m

> x

q>m

> x

l>0

> x

l>q

= given= (2.4.8)

The matrix version of the discrete 1D model with

q = 6 is as follows. This

1D model will have 5 unknowns, which we list in classical order from left to

right. The matrix

D will be 5 × 5 and is derived from(2.4.5) by dividing both

sides by

 = N@k

2 1

1 2 1

1 2

= (1@)

The matrix version of the discrete 2D model with q = 6 will have 5

= 25

unknowns. Consequently, the matrix

D will be 25 × 25. The location of its

components will evolve from line (2.4.7) and will depend on the ordering of the

unknowns

. The classical method of ordering is to start with the bottom

grid row (

m = 1) and move from left (l = 1) to right (l = q  1) so that

x =

with X

is a grid row m of unknowns. The ﬁnal grid row corresponds to m = q  1. So,

it is reasonable to think of

D as a 5 × 5 block matrix where each block is 5 × 5

and corresponds to a grid row. With careful writing of the equations in (2.4.7)

one can derive

D as

2.4. SCHUR COMPLEMENT 79

E L

L E L

L E

= (1@)

where

E =

4 1

1 4 1

1 4

and L =

2.4.4 Method

In the above 5×5 block matrix it is tempting to try a block version of Gaussian

elimination. The ﬁrst block row could b e used to eliminate the

L in the block

(2,1) position (block row 2 and block column 1). Just multiply block row 1 by

1

and add the new block row 1 to block row 2 to get

0 (

E  E

1

) L 0 0

where the 0 represents a 5 × 5

}hur pdwul{. If all the inverse matrices of any

subsequent block matrices on the diagonal exist, then one can continue this

until all blocks in the lower block part of

D have been modiﬁed to 5 × 5 zero

matrices.

In order to make this more precise, we will consider just a 2×2 block matrix

where the diagonal blocks are square but may not have the same dimension

D =



E H

I F

= (2.4.9)

In general

D will be q × q with q = n + p, E is n × n, F is p × p, H is n × p

and I is p×n= For example, in the above 5×5 block matrix we may let q = 25,

n = 5 and p = 20 and

F =

E L

L E L



L E L

L E

and H = I

L 0 0 0

If E has an inverse, then we can multiply block row 1 by I E

1

and sub-

tract it from block row 2. This is equivalent to multiplication of

D by a block

elementary matrix of the form



I E

1

80 CHAPTER 2. STEADY STATE DISCRETE MODELS

D{ = g is viewed in block form, then



E H

I F

¸

[



= (2.4.10)

The above blo ck elementary matrix multiplication gives



E H

0 F  I E

1

¸

[



 I E

1

= (2.4.11)

So, if the block upper triangular matrix is nonsingular, then this last block

equation can be solved.

The following basic properties of square matrices play an important role in

the solution of (2.4.10). These prop erties follow directly from the deﬁnition of

an inverse matrix.

Theorem 2.4.1 (Basic Matrix Properties) Let B and C be square matrices

that have inverses. Then the following equalities hold:



E 0

1



1

0 F

1



I L

1



I L



E 0

I F



¸

1

I L

and



I F

1



1

I L

1



1



1

F

1

I E

1

Deﬁnition. Let D have the form in (2.4.9) and E be nonsingular. The Schur

complement of

E in D is F  I E

1

Theorem 2.4.2 (Schur Complement Existence) Consider A as in (2.4.10).

If both B and the Schur complement of B in A are nonsingular, then A is

nonsingular. Moreover, the solution of

D{ = g is given by using a block upper

triangular solve of (2.4.11).

The choice of the blocks

E and F can play a very important role. Often

the choice of the physical object, which is being modeled, suggests the choice

of E and F= For example, if the heat diusion in a thin wire is being modeled,

the unknowns associated with E might be the unknowns on the left side of the

thin wire and the unknowns associated with

F would then be the right side.

Another alternative is to partition the wire into three parts: a small center and

a left and right side; this might be useful if the wire was made of two types of

materials. A somewhat more elaborate example is the model of airﬂow over an

aircraft. Here we might partition the aircraft into wing, rudder, fuselage and

"connecting" components. Such partitions of the physical object or the matrix

are called domain decompositions.

2.4. SCHUR COMPLEMENT 81

2.4.5 Implementation

MATLAB will be used to illustrate the Schur complement, domain decompo-

sition and di

erent ordering of the unknowns. The classical ordering of the

unknowns can be changed so that the "solve" or "inverting" of B or its Schur

complement is a minimal amount of work.

1D Heat Di

usion with n = 6 (5 unknowns).

Classical order of unknowns

> x

gives the coe!cient matrix

D =

2 1

1 2 1

1 2

Domain decomposition order of unknowns is x

; x

> x

; x

> x

= In order to

form the new co e!cient matrix D

, list the equations in the new order. For

example, the equation for the third unknown is

x

+ 2x

 x

= (1@)i

, and

so, the ﬁrst row of the new coe

!cient matrix should be

2 0

1 1 0

The other rows in the new coe!cient matrix are found in a similar fashion so

that

2 1 1

1

1 1 2

1 2 1

1 2

Here E = [2] and F is block diagonal. In the following M

ATLA B calculations

note that

E is easy to invert and that the Schur complement is more complicated

than the F matrix.

Ab = [2];

Ae = [0 - 1 -1 0];

Af = e’;

Ac = [2 -1 0 0;-1 2 0 0;0 0 2 -1;0 0 -1 2];

Aa = [b e;f c]

a =

2 0 -1 -1 0

0 2 -1 0 0

-1 -1 2 0 0

-1 0 0 2 -1

0 0 0 -1 2

Aschurcomp = c - f*inv(b)*e % 4x4 tridiagonal matrix

schurcomp =

2.0 -1.0 0 0

-1.0 1.5 -0.5 0

82 CHAPTER 2. STEADY STATE DISCRETE MODELS

0 -0.5 1.5 -1.0

0 0 -1.0 2.

Ad1 = [1];

Ad2 = [1 1 1 1]’;

Add2 = d2 - f*inv(b)*d1

dd2 =

1.0000

1.5000

1.0000

Ax2 = schurcomp\dd2 % block upper triangular solve

x2 =

2.5000

4.0000

2.5000

Ax1 = i nv(b)*(d1 - e*x2)

x1 =

4.5000

Ax = a\[d1 d2’]’

x =

4.5000

2.5000

4.0000

2.5000

Domain decomposition order of unknowns is

> x

; x

> x

; x

so that the

new coe!cient matrix is

2 1

1 2 1

2 1 1

1 2

1 1 2

Here F = [2] and E is block diagonal. The Schur complement of E will be 1 × 1

and is easy to invert. Also,

E is easy to invert b ecause it is block diagonal. The

following M

AT LAB calculations illustrate this.

Af = [ 0 -1 -1 0];

Ae = f’;

Ab = [2 -1 0 0;-1 2 0 0;0 0 2 -1;0 0 -1 2];

Ac = [2];

Aa = [ b e;f c]

2.4. SCHUR COMPLEMENT 83

a =

2 -1 0 0 0

-1 2 0 0 -1

0 0 2 -1 -1

0 0 -1 2 0

0 -1 -1 0 2

Aschurcomp = c -f*inv(b)*e % 1x1 matrix

schurcomp =

0.6667

Ad1 = [1 1 1 1]’;

Ad2 = [1];

Add2 = d2 -f*inv(b)*d1

dd2 =

Ax2 = schurcomp\dd2 % block upper triangular solve

x2 =

4.5000

Ax1 = inv(b)*(d1 - e*x2)

x1 =

2.5000

4.0000

2.5000

Ax = inv(a)*[d1’ d2]’

x =

2.5000

4.0000

2.5000

4.5000

2D Heat Di

usion with n = 6 (25 unknowns).

Here we will use domain decomposition where the third grid row is listed

last, and the ﬁrst, second, fourth and ﬁfth grid rows are listed ﬁrst in this order.

Each block is 5 × 5 for the 5 unknowns in each grid row, and

l is a 5 ×5 identity

84 CHAPTER 2. STEADY STATE DISCRETE MODELS

matrix

e l

l e l

e l l

l e

l l e

where

e =

4 1

1 4 1

1 4

The E will be the block 4×4 matrix and F = e. The B matrix is block diagonal

and is relatively easy to invert. The C matrix and the Schur complement of

E are 5 × 5 matrices and will be easy to invert or "solve". With this type of

domain decomposition the Schur complement matrix will be small, but it will

have mostly nonzero components. This is illustrated by the following M

ATLAB

calculations.

Aclear

Ab = [4 -1 0 0 0;-1 4 -1 0 0; 0 -1 4 -1 0; 0 0 -1 4 -1;0 0 0 -1 4];

Aii = -eye(5);

Az = zeros(5);

AB = [b ii z z;ii b z z; z z b ii; z z ii b];

Af = [z ii ii z];

Ae = f’;

AC = b;

Aschurcomp = C - f*inv(B)*e % full 5x5 matrix

schurcomp =

3.4093 -1.1894 -0.0646 -0.0227 -0.0073

-1.1894 3.3447 -1.2121 -0.0720 -0.0227

-0.0646 -1.2121 3.3374 -1.2121 -0.0646

-0.0227 -0.0720 -1.2121 3.3447 -1.1894

-0.0073 -0.0227 -0.0646 -1.1894 3.4093

Awhos

Name Size Bytes Class

B 20x20 3200 double array

C 5x5 200 double array

b 5x5 200 double array

e 20x5 800 double array

f 5x20 800 double array

ii 5x5 200 double array

schurcomp 5x5 200 double array

z 5x5 200 double array

2.4. SCHUR COMPLEMENT 85

2.4.6 Assessment

Heat and mass transfer models usually involve transfer in more than one direc-

tion. The resulting discrete models will have structure similar to the 2D heat

diusion model. There are a number of zero components that are arranged in

very nice patterns, which are often block tridiagonal. Here domain decomposi-

tion and the Schur complement will continue to help reduce the computational

burden.

The proof of the Schur complement theorem is a direct consequence of using

a block elementary row op eration to get a zero matrix in the block row 2 and

column 1 position



I E

1

¸

E H

I F



E H

0 F  I E

1

Thus



E H

I F



I E

1

¸

E H

0 F  I E

1

Since b oth matrices on the right side have inverses, the left side, D, has an

inverse.

2.4.7 Exercises

1. Use the various orderings of the unknowns and the Schur complement to

solve

D{ = g where

D =

2 1

1 2 1

1 2

and g =

2. Consider the above 2D heat diusion model for 25 unknowns. Suppose

g is a 25×1 column vector whose components are all equal to 10. Use the Schur

complement with the third grid row of unknowns listed last to solve

D{ = g.

3. Repeat problem 2 but now list the third grid row of unknowns ﬁrst.

4. Give the proofs of the four basic properties in Theorem 2.4.1.

5. Find the inverse of the block upper triangular matrix



E H

0 F

6. Use the result in problem 5 to ﬁnd the inverse of



E H

0 F  I E

1

7. Use the result in problem 6 and the p roof of the Schur complement the-

orem to ﬁ nd the inverse of

86 CHAPTER 2. STEADY STATE DISCRETE MODELS



E H

I F

2.5 Convergence to Steady State

2.5.1 Introduction

In the applications to heat and mass transfer the discrete time-space dependent

models have the form

n+1

= Dx

+ e=

Here x

n+1

is a sequence of column vectors, which could represent approximate

temperature or concentration at time step

n + 1= Under stability conditions on

the time step the time dependent solution may "converge" to the solution of

the discrete steady state problem

x = Dx + e=

In Chapter 1.2 one condition that ensured this was when the matrix products

"converged" to the zero matrix, then x

n+1

"converges" to x. We would like

to be more precise about the term “converge” and to show how the stability

conditions are related to this "convergence."

2.5.2 Vector and Matrix Norms

There are many dierent norms, which are a "measure" of the length of a vector.

A common norm is the Euclidean norm

 ({

{)

Here we will only use the inﬁnity norm. Any real valued function of { 5 R

that satisﬁes the properties 1-3 of subsequent Theorem 2.5.1 is called a norm.

Deﬁnition. The inﬁnity norm of the

q × 1 column vector { = [{

] is a real

number

k{k  max

| =

The inﬁnity norm of an q × q matrix D{ = [d

] is

Dk  max

| =

Example.

Ohw { =



9

and D =

1 3 4

1 3 1

3 0 5

k{k = max{1> 6> 9} = 9 and kDk = max{8> 5> 8} = 8=

2.5. CONVERGENCE TO STEADY STATE 87

Theorem 2.5.1 (Basic Properties of the Inﬁnity Norm) Let

D and E be q×q

matrices and {> | 5 R

. Then

1. k

{k  0, and k{k = 0 if and only if { = 0,

2. k

{ + |k  k{k + k|k >

3. k{k  || k{kwhere  is a real number,

4. k

D{k  kDk k{k and

5. k

DEk  kDk kEk =

Proof. The proofs of 1-3 are left as exercises. The pro of of 4 uses the

deﬁnitions of the inﬁnity norm and the matrix-vector product.

D{k = max

{

 max

| · |{

 (max

|) · (max

|) = kDk k{k =

The proof of 5 uses the deﬁnition of a matrix-matrix product.

DEk  max

 max

| |e

= max

 (max

|)(max

= k

Dk kEk

Prop erty 5 can be generalized to any number of matrix p roducts.

Deﬁnition. Let x

and { be vectors. {

converges to { if and only if each com-

ponent of

{

converges to {

. This is equivalent to

{

 {

= max

{

 {

converges to zero.

Like the geometric series of single numbers the iterative scheme

{

n+1

+ e can be expressed as a summation via recursion

{

n+1

= D{

+ e

= D(D{

n1

+ e) + e

= D

{

n1

+ De + e

= D

(D{

n2

+ e) + De + e

= D

{

n2

+ (D

+ D

+ L)e

n+1

{

+ (D

+ · · · + L)e= (2.5.1)