Elliptic Equations and Iterative Methods

3 . 6 . 1 I N T R O D U C T I O N

T h e treatment of initial-value problems involved the formation of a difference operator that permitted the initial conditions to be extended into the domain of interest of the problem. Such a procedure is generally not useful for elliptic equations where the boundary conditions are given over the entire region of interest. T h e numerical solution of elliptic equations is usually accomplished by solution of simultaneous equations with a variety of methods.

One possible means of solving a set of simultaneous equations is by the Gauss reduction scheme. Unfortunately the reduction process for Ν equations in Ν unknowns requires approximately N³ operations.

Furthermore, a certain amount of round-off in each operation may cause the solution to degenerate for large N. On the other hand, a direct reduction procedure is determinate in that a fixed number of steps are needed (in theory) to find the solution.

A n alternative approach to the solution of elliptic equations is an iterative procedure. In general, iterative methods require an infinite number of steps to solve a problem exactly. However, for practical purposes it is usually possible to terminate an iteration after a finite number of steps which are fewer in number than those required for reduction methods. Furthermore, iterative procedures have certain advantages with respect to round-off over direct reduction. T o make these matters clearer, we consider a simple introductory example.

W e desire the solution of the problem

(3.6.1)

with boundary conditions

M O ) =

^{j o};

M

^{f l}

) =

ya ·

W e divide the interval 0 to α into Κ subintervals of equal width h as in Fig. 3.6.1. For simplicity we replace the second derivative with a second

*| h k k+l x=a

F I G . 3.6.1. Simple mesh for the numerical solution of Laplace's equation.

central difference. W e have then

yic-i - 2yk + yk 1+ = 0, 1 < A < * - 1 (3.6.2) T h e set of equations (3.6.2) could be solved by reduction (in this case by matrix factorization, see Section 1.14) in a straightforward manner.

T h e iterative solution of the equations is achieved by first assuming a trial solution at all the interior points of the mesh. W e then write Eq.

(3.6.2) in the form

. yk +1 + yk-j

yk ='• (3.6.3)

Equation (3.6.3) provides an algorithm for computing values of yk in terms of its nearest neighbors. L e t us denote the values of yk for the initial trial with a superscript 0. T o compute the " n e w " values of yk

from the " o l d , " we might consider the rule

and in general

yk — ό

ν yl+l + yl-\

yk =

(3.6.4)

(3.6.5) T h e following important questions now occur:

1. W i l l the iteration rule (3.6.5) ever yield a solution?

2. H o w will we know when we have achieved a solution ? 3. H o w long will it take ?

W e assume that the difference approximation itself has been adequate for the problem.

T o help answer these questions, let us denote the actual solution to the difference equation as y%. W e define the initial error, say e£ , as

<l=yl-yt

(3-6.6)

and after p steps

<ΐ=νΙ-νΐ·

(3-6.7)

T h e answer to the first question can be given in terms of the errors. Thus, if the sequence β£ , , ... , e£ approaches zero, for all k, then the iteration does yield a solution. T o answer the second question is somewhat more difficult in that we never know the errors (else we would also know the solution). T h e usual test for determining when the solution has been achieved is based upon the following criterion. W e define a function, say rA., as

n=y'^

^l

-y^

(3.6.8)

In terms of the errors we have

rl = (3·⁶·9)

If r£ is small, for all k, then one assumes the approximate solution has been found. Obviously this criterion is not always valid since a small difference of errors does not necessarily imply small errors. A practical way to be reasonably sure of convergence is to choose a sufficiently small criterion and to iterate a few times perhaps beyond the criterion.

T h i s second feature guards against the rate of convergence as a function of the iteration index having a minimum. For example, the convergence rate may initially increase only to later decrease. A minimum might result while still far from converged.

T h e question of how long an iteration will take is very important, and many different iteration methods have been devised to hasten the rate at which a solution is obtained. W e shall introduce a variety of methods in later sections, and one of the criteria for judging the merit of a method will be the rate of convergence.

3 . 6 . 2 S T A B I L I T Y OF I T E R A T I O N S

In this section we shall attempt to provide a unified basis for the study of iterative methods. W e assume the set of simultaneous equations are written in the form

A x = y, ( 3 . 6 . 1 0 )

where we further assume a solution does exist but otherwise do not restrict the matrix A . T h e exact solution to the problem, say x^{0 0}, is

x⁰⁰ = A !y. ( 3 . 6 . 1 1 )

W e begin the iteration by considering a trial solution x ° and then operate on x ° to produce a new trial x¹. W e assume the iteration can be written

x¹ = B x ° + z, ( 3 . 6 . 1 2 )

where the matrix Β and the vector ζ are taken independent of the iteration index. Such an iteration is called stationary. T h e iterative procedure can be extended to successive trials in the form

XM I = B x " + z. ( 3 . 6 . 1 3 )

T h e matrix Β is called the iteration matrix, and the ability of achieving a solution and the corresponding rate of convergence are intimately related to the properties of the iteration matrix.

T h e properties we desire of our iteration are that the sequence of vectors x^p approach x^x, and further, that iteration with x⁰⁰ reproduces itself. In mathematical terms the last requirement is

x⁰⁰ = B x * + z . ( 3 . 6 . 1 4 )

Assuming Eq. ( 3 . 6 . 1 4 ) is valid, and defining the error vector e^p as

c" = x " - x^{0 0}, ( 3 . 6 . 1 5 )

we have from Eq. ( 3 . 6 . 1 3 )

c" = B e * -¹. ( 3 . 6 . 1 6 ) Equation ( 3 . 6 . 1 6 ) states that the error vector obeys the homogeneous

form of the iteration equation. N o w , if the sequence of vectors x ^ is to approach x^{0 0}, then the sequence of vectors must approach zero.

From Eq. ( 3 . 6 . 1 6 ) we have

€" = B e " -¹ = B"€°. ( 3 . 6 . 1 7 )

W e require

(3.6.18) lim €*> = lim B^p€° -> 0.

L e t us now assume the matrix Β has a complete set of eigenvectors, say , and corresponding eigenvalues . Since the eigenvectors are complete, we may expand e° in the form

(3.6.19)

(3.6.20)

(3.6.21)

(3.6.22) Operation on ε° by Β yields c¹ as

By induction we have

In order for the error vector to vanish, we must require²

T h e above result is very important for our future considerations and is also quite general. W e shall refer to Eq. (3.6.22) as the stability condition for iterative methods. I f an eigenvalue, say Xr, was of magnitude greater than unity, the iteration would diverge, and in our terminology is unstable. Even if the initial error vector were orthogonal to the eigen­

vector, say er, of the eigenvalue Xr , round-off in computations would introduce components along er and ultimately ruin the iteration.

T h e stability condition also provides a means of computing the convergence rate of an iteration. I f the eigenvalues are ordered such that

(3.6.23) then for sufficiently large p, the error is approximately

T h e term e^lnXi is the decay factor for the errors. In particular we define (3.6.24)

2 The requirement of completeness of the eigenvectors of Β to derive condi­

tion (3.6.22) is over restrictive. By consideration of the Jordon canonical form, the same result can be achieved by any real matrix. See problem 7.

as the convergence rate of the iteration. For small λχ the convergence rate is large, meaning a rapid reduction of the error with increasing number of iterations. T h e worth of an iteration scheme is partly measured in terms of the factor v. W e shall consider detailed examples later. W e note, however, that in general ν decreases for increasing number of unknowns for most iterations. T h a t is, the larger the problem, the longer it takes to solve. In fact, for many methods convergence is reciprocally related to TV², where Ν is the number of unknowns. T h a t the relation goes as N² can be heuristically seen by the following argu­

ment. Since the eigenvalues must be distributed between —1 and + 1 , an increase in Ν should move the ones o£ largest magnitude toward the end points of the interval. Further, since there are more points to be solved, the number of operations grows with N. T h u s we expect the convergence to be related somehow to TV².

T h e above argument is strictly heuristic and not necessarily true.

Later examples will be considered to illustrate the general behavior of the convergence factor.

T h u s far our results have been derived by strictly algebraic considera­

tions. It is interesting and instructive to consider a geometric inter­

pretation of iterations. T o this end, we introduce the residual vector r^p defined as

r" = x^{p +1} - x". (3.6.25)

In terms of the error vector, Eq. (3.6.25) becomes

r* = ( B - iy. (3.6.26) Since the residual vector represents the error vector in a transformed

space, the convergence criteria for the residuals are the same as for the errors. Furthermore, the asymptotic behavior is the same.

T h e residuals are calculable at any stage of the iteration. F r o m the defining equation for the iteration (3.6.13), we have

r" = ( B - I ) x » + z. (3.6.27) W e interpret Eq. (3.6.25) in the form

xp+i = vx + VP (3.6.28) as stating that the vector x^p is corrected by addition of a vector r^p

which in turn is defined by an algorithm (3.6.27). Different iteration methods consist of different algorithms for computing the correction vector. Notice that if r^p is chosen to change only one component of x^p, the vector r^p+1 may still have all of its components changed since

x^p+l is transformed by an operator that is not diagonal. In any event, convergence of the solution requires that x^p approach a limit vector and r^p approach the null vector.

T h e change of the trial vector x^,J by the correction vector r^p is called an iterate or one iteration. T h e change of each component separately is called a relaxation or displacement. Thus, after relaxing each component of xpy once and only once, we complete one iteration.

3 . 6 . 3 S T A T I O N A R Y I T E R A T I O N S

In this section we shall discuss several very common iteration methods and consider their stability properties. For each of the iteration methods introduced, we shall apply the method to the simple Laplace equation in the square.³ That is, we consider the equation

- § - + - p - =' 0 ° < * < β · 0 O < « ( 3 - 6 . 2 9 ) with boundary conditions

φ(0,γ) =φ(α,ν) = 0 ,

φ(χ, 0 ) = 0 , φ{χ, a) = f(x). ( 3 . 6 . 3 0 ) For all iteration methods we use the difference equation

& + - # = 0 , hx = hy = h. ( 3 . 6 . 3 1 ) h^{2 1} h²

x y

A. Method of Simultaneous Displacements

L e t us assume we are solving the matrix equation

A x = y, ( 3 . 6 . 3 2 ) where we assume a solution does exist. T h e matrix A is written in the

form

A = L + D + U , ( 3 . 6 . 3 3 ) where L is strictly lower triangular, U strictly upper triangular, and D

diagonal. W e assume D Φ 0. W e now define the iteration

( L + U ) x " + D x "⁺¹ - y, ( 3 . 6 . 3 4 )

3 The numerical solution of the two-dimensional Laplace equation is perhaps the most thoroughly studied problem in iterative analysis. T h e problem is some­

times referred to as the model problem.

or

XP + I = _ D" i( L + U ) x " + D V · (3.6.35) T h e method of simultaneous displacements then consists in solving the difference equation by means of the iteration in Eq. (3.6.35).

W e next observe that the method of iteration in consistent with a vector being a solution of our original Eq. (3.6.32) for if at some stage the correct answer were found, then, apart from round-off errors, iterating the solution does nothing more than give it back to us again.

In particular, if at some stage x^IJ were the correct solution

A - ' y , (3.6.36)

then by substituting this solution into the iteration scheme (3.6.35), we find that

χρ+ 1 = - l y > A

and we recover the original solution.

T h e important question of whether or not an arbitrary starting vector, say x°, will approach A^{- 1} y depends upon the iteration matrix

— D^{_ 1}( L + U ) . T h e eigenvalue spectrum of the iteration matrix must be such that | λί | < 1 for all /, where the are eigenvalues of the matrix.

Thus, the roots of the equation

or, equivalently

- D - H L + U ) - λί I = 0

I L + AD + U I = 0,

(3.6.37)

(3.6.38) must have magnitude less than unity.

In general the solution of (3.6.38) is quite difficult. However, for many cases of interest, it is possible to determine the nature of the eigenvalue spectrum without solving the secular equation. W e assume that the problem is well posed in the sense that A is irreducible, otherwise we factor A into two separate problems. T h e iteration matrix can be written

D~HL + U ) =

By Gersehgorin's theorem, the spectral radius of D^{_ 1}( L + U ) is bounded by

! Am a x| <

max2) IS ·

^(3-6-40)

Gerschgorin's theorem may be strengthened if the matrix (3.6.39) is irreducible. I f the maximum value of the sum on the right-hand side of Eq. (3.6.40) is denoted p, then Gerschgorins' theorem states

I

^{Am a x}

I

< P- (3.6.41)

If, for any i, and for an irreducible matrix

X ^af < P , (3.6.42)

then

I Am ax I < P, (3.6.43)

i.e., there is strict inequality.⁴ A s a consequence, we see that if the matrix A is such that

I *u ati I (j Φ i) (3.6.44)

j

with inequality for some i,then the method of simultaneous displacements converges.

I f the diagonal elements of A are such that (3.6.42) is true, then we say A has diagonal dominance. Consequently, for an irreducible matrix with diagonal dominance, the method of simultaneous displacements converges. Fortunately most elliptic difference equations of reactor interest have this property. N o t e that the condition is sufficient for convergence but not necessary. For instance, the matrix

A = [-J j (3.6.45)

does not have diagonal dominance, yet the iteration matrix

- D - H L + U ) = β/5 ^ (3.6.46)

has eigenvalues with magnitude less than unity.

4 For proof see Reference 8, Chapter 1.

W e now consider application of the method of simultaneous dis­

placements to the model problem. T h e equation for the error function is the homogeneous form of the equation with homogeneous boundary conditions.⁵ T h e iteration algorithm for the error e?k is, f r o m E q . (3.6.34) e?+ 1.fc + €?_l i fc - fcft¹ + €,%+1 + ej^ - 2efX^l = 0, · (3.6.47) or

pfl 1 rv J r V ι V ι Ρ ι V η (3.6.48) T h e stability of the iteration can be studied by examining the behavior of the various eigenfunctions of the iteration operator. F r o m past results we know the error (with homogeneous boundary conditions) must be of the form

€( m , n ) = w) m sin s

Κ 1 < k < Κ - 1 (3.6.49) where (m, ri) are the indices of the eigenfunction. If we assume e j ^^{, n)} is of the form (3.6.49) and insert in (3.6.48), we have

p +1 _ A(tn, n) After an elementary reduction we have

p+i Αί \ · n7TJ · k m n T h e scale factor ζ determines the rate of growth or decay for the errors.

For this case the errors all decay with differing rates for the different harmonics. T h e largest value of ζ occurs for m = η = 1 or m — η = Κ — 1. For such a case we have

ζ = cos

Κ 1

-2K* ' (3.6.53)

T h e iteration converges with the convergence rate _2

• The argument below follows that of Frankel, Reference 12.

(3.6.54)

Notice that the convergence rate is proportional to K², the number of unknowns. For a mesh of 20 by 20 points, the convergence rate is approximately

υ = 0.01234

and hence about 80 iterations are required to reduce the error by a factor of e.

Notice that the method of simultaneous displacements requires the vectors x^p+1 and x'^J in the computation simultaneously. T h i s means that for a vector of Ν dimensions, 2N numbers must be retained. For large problems this condition may be important in working with digital computers.

T h e method of simultaneous displacements is easily described in terms of the residuals. By Eqs. (3.6.35) and (3.6.28) we have

Γ" = - D - ^ A x " - y ] . (3.6.55) T h e quantity Ax^p — y is a measure of the error in a given iteration.

Some authors call Αχ^μ — y the residual.

A s a final note we remark that the method of simultaneous displace­

ments is also called the Jacobi method or the Richardson method.

B. Method of Successive Displacements

T h e method of successive displacements (also called the Gauss-Seidel method or the Liebmann method) is quite similar to the previous scheme. W e again consider the problem of solving the matrix equation (3.6.32). In this case we assume an iteration of the form

(L + D ) x "⁺¹ + U x " = y (3.6.56) or

x "⁺¹ = - ( L + D ^ U x " + ( L + D ) - ^ , (3.6.57) which is the iteration algorithm that defines the method of successive

displacements.

A s before, it is easily seen that if at some state x^ = A^{- 1} y , then

xp + i = A "¹ y .

W e now consider the convergence properties of the iteration and derive expressions for the eigenvalues of the iteration operator. Assume I D I ^ 0, and since L is strictly lower triangular, the matrix ( L + D )^{- 1} exists and is of the form

( L + D )¹ = L ' + D -1 > (3.6.58)

where L ' is strictly lower triangular. Convergence of the iteration depends upon the eigenvalues of the iteration matrix, i.e., the roots λ of

| - ( L + D )¹! ! — AII = 0 (3.6.59) or

I AL + AD + U j = 0. (3.6.60) T h e roots of Eq. (3.6.60) must all have magnitude less than unity for

convergence. It is usually difficult to solve for the roots explicitly;

however, the rule concerning diagonal dominance of irreducible matrices applies for the method of successive displacement also.

In the method of successive displacements, we are always using the latest computed values for the unknown and hence the name successive displacement. It is interesting to compare the convergence rates of the method of simultaneous and successive displacements for the same problem. W h e n the iteration matrices are non-negative, we shall prove shortly that the two methods converge or diverge together. Further, if they converge, then the successive displacements technique converges more rapidly.

T o prove this result we consider first the method of simultaneous displacement. W e assume that the iteration matrix is nonnegative which is possible if the matrices L and U have nonnegative elements and D has all negative nonzero components, or conversely, for instance. For the method of simultaneous displacements, the iteration matrix is then

- D - X L + U ) = R + T , (3.6.61) where R is strictly lower triangular and Τ strictly upper triangular.

Similarly, the iteration matrix for the method of successive displacements is

- ( L + D ) "^XU == (I - R ) *T. (3.6.62) T h e proof (from Reference 13) of divergence is straightforward. L e t A be

the positive eigenvalue of R - f Τ of greatest magnitude, and similarly σ be the largest eigenvalue of ( I — R )^{_ 1}T . L e t ζ be the eigenvector corresponding to A; that is,

( R + T ) z = λζ. (3.6.63) Therefore

( i - A )^{- 1} ( R + T ) z = λ ( I - ζ. (3.6.64)

N o w the matrix ( I — R / λ )^{- 1} can be written

( ' - - I T

where R is assumed (m + 1) by (m + 1). T h u s terms beyond the mth power of R vanish, since R is strictly lower triangular. N o t e that all the elements of the sum are nonnegative in view of the hypothesis. Using the expansion (3.6.65) in Eq. (3.6.64), we have

( i - - ^ - )^{_ 1}T z = λ ζ . (3.6.66) T h u s λ is an eigenvalue of ( I — R / A )^{_ 1}T . F r o m the properties of

non-negative matrices, if λ > 1, then

( i - Τ < (I - R ) " ^ . (3.6.67) Consequently σ > λ > 1, and the two iterations diverge together. I f

λ = 1, then σ = 1, while for λ < 1 , σ < λ < 1 , which proves the result.

Conversely we could reverse the arguments. L e t σ be the largest positive eigenvalue of ( I — R )^{_ 1}T and ζ the corresponding eigenvector.

T h u s

(I - R ^ T z = σζ, (3.6.68) or

( a R + T ) z = σζ. (3.6.69) T h e r e f o r e σ is also an eigenvalue of a R + T . If σ > 1, then a R + Τ >

R + Τ and σ > λ > 1. For σ = 1, then λ = 1. Finally if σ < 1

aR + Τ < R + Τ . (3.6.70) Therefore

σ < λ < 1. (3.6.71) T h i s last result is the important result as it shows that for nonnegative

iteration matrices the method of successive displacements asymptotically converges faster than the method of simultaneous displacements.

For matrices which are not nonnegative, it is possible for one method to work and not the other, and vice versa (see problem 9).

W e now illustrate the method by considering the model

^ = W l r (

³

·

⁶

·

^{7 2}

)

Stability requires that | λ/ η >η | < 1 for all m, n. Using Eq. (3.6.72) in the difference equation (3.6.56) and factoring, we have

= i l & i . * + W- i , * + ⁶^+i + ^A

m.„«,VJ- (

³

·

⁶

·

^{7 3}

)

Following Frankel (Ref. 12) we assume eigenfunctions of the form

€( m , n ) = Ai ns H^Li Bk _{s i n}

ig*

^{9 (3.6.74)}

where A and Β are assumed constants to be determined. Expanding the trigonometric functions in Eq. (3.6.73) yields the equation

+ + XminBk-i)A> cos _ em sin —

* ^J * ^K (3.6.75)

+ 5 - Am,n^ ' i ) S * sin — sin cos g ?

-+ ( B ^ - Xm,nB«-i)A> sin ^ sin ψ cos .

Since the error must be zero on the boundary, viz. at k or j = 0 or k or

= K, the terms in the cosine must vanish. Hence we require

A M - X m, n A i - i = 0

(3.6.76) - λ„,,„£*-ι = 0.

Therefore,

A* = B * = X m .t (3.6.77) n

Equation (3.6.75) then becomes

ν 1 / ηπ τηπ \²

λτη,η = 4 (cos -γ- + cos - ^ - j . (3.6.78)

T h e maximum value is again found for m = η = 1. Expanding we have

Am,„ ^ 1 (3.6.79)

T h e asymptotic decay rate is then

« = - £ Γ · (3-6.80) Notice that this rate is twice as large as for the method of simultaneous

displacements (see Eq. (3.6.54)). T h u s , we expect the method of suc­

cessive displacements to take roughly half as lorrg for the model problem as the method of simultaneous displacements. T h i s result is consistent with the general result for nonnegative irreducible matrices since the Laplace difference equation (3.6.31) leads to a nonnegative irreducible iteration matrix.

In terms of the residuals the method of successive displacements becomes

χν+ι = XP + TPY (3.6.81)

where

Γ" = — (L + D ^ X A x " - y ) . (3.6.82) It is interesting to write out the component form of the residual for the

model problem: we have

= iKi.*

+ * Γ ί * + + ^ , t - i - 4 * g - (3.6.83) A similar result applies for the method of simultaneous displacements where all terms with superscript p + 1 are replaced by like terms with superscript/). Notice that the residual can be interpreted as the inbalance between the function xpk and the value of the difference relation operating on the function at the pth iterate.

C. Successive Over-Relaxation

For both of the previous methods, the iteration algorithm could be written

xp+i = xv + rp (3.6.84)

with different r^p for different methods. From the discussion above, we interpret the residual as correcting the function at each point, say j, k> so as to satisfy the difference equation. Obviously, if any neighboring

point to j, k is changed, the residual also changes, as illustrated by Eq.

(3.6.84). W e might anticipate further changes in the function by over-correcting (or perhaps under-over-correcting) in hopes of speeding conver­

gence of the iteration. T h e iteration might then be written

xp + i = x/ > + Ρ 5 (3.6.85) α Γ

where a is a real number. For α > 1 we speak of over-relaxation; for

where a is a real number. For α > 1 we speak of over-relaxation; for

In document NUMERICAL SOLUTIONS OF EQUATIONS Ill (Pldal 25-63)