ON LOCAL MODE ANALYSIS IN MULTI-GRID METHODS By MARK ALLAN REIMERS B. Sc., The U n i v e r s i t y of Toronto, 1978 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES (Department of Mathematics) We accept t h i s t h e s i s as conforming to the r e q u i r e d standard THE UNIVERSITY OF BRITISH COLUMBIA August ©Mark 1983 A l l a n Reimers, 1983 In p r e s e n t i n g t h i s t h e s i s i n p a r t i a l f u l f i l m e n t of requirements f o r an advanced degree a t the the University of B r i t i s h Columbia, I agree t h a t the L i b r a r y s h a l l make it f r e e l y a v a i l a b l e f o r reference and study. I further agree t h a t p e r m i s s i o n f o r e x t e n s i v e copying o f t h i s t h e s i s f o r s c h o l a r l y purposes may be department or by h i s o r her granted by the head of representatives. my It is understood t h a t copying or p u b l i c a t i o n of t h i s t h e s i s f o r f i n a n c i a l gain s h a l l not be allowed without my permission. Department of The U n i v e r s i t y of B r i t i s h 1956 Main Mall Vancouver, Canada V6T 1Y3 DE-6 (.3/81) Columbia written ii ABSTRACT L o c a l Mode A n a l y s i s was introduced by Achi a h e u r i s t i c p r a c t i c a l t o o l f o r determining convergence r a t e of a that t h i s t o o l may be j u s t i f i e d computations yield Multi-Grid analogous to the algorithm. rigorously. actual i s also useful of c e r t a i n r e l a x a t i o n s context. expected I t i s shown Furthermore, relaxation r e s u l t s much l i k e L o c a l Mode A n a l y s i s This analysis Brandt as processes predictions. f o r a h e u r i s t i c understanding independently of the Multi-Grid iii TABLE OF CONTENTS Abstract ii List of T a b l e s v List of F i g u r e s .....vi Acknowledgement Chapter vii 1 - Introduction 1 . 1 The Subject 1 1.2 O u t l i n e of t h i s T h e s i s 4 Chapter 2 - Background 2.1 F i n i t e D i f f e r e n c e S o l u t i o n of E l l i p t i c P a r t i a l D i f f e r e n t i a l Equations 6 2.2 The M u l t i - G r i d (MG) Method 11 2.3 L o c a l Mode A n a l y s i s 17 Chapter 3 - A J u s t i f i c a t i o n of L o c a l Mode A n a l y s i s 3.1 The Meaning of L o c a l Mode A n a l y s i s 21 3.2 Accuracy of L o c a l Mode A n a l y s i s 23 3.3 Examples 35 3.4 Comments and M i s c e l l a n e o u s Consequences Chapter 4 - A H e u r i s t i c Look a t R e l a x a t i o n s 41 using L o c a l Mode A n a l y s i s 4. 1 The SOR Method 46 4.2 A l t e r n a t i n g D i r e c t i o n I m p l i c i t 4.3 Smoothing Property Iteration of R e l a x a t i o n s 51 55 iv Chapter 5 - D i s c u s s i o n 5.1 A p p l i c a b i l i t y of L o c a l Mode A n a l y s i s 57 5.2 The P o t e n t i a l of LMA Outside the MG Context..58 5.3 P o s s i b l e Developments of the Theorem 59 5.4 Summary 61 Bibliography 62 Appendix A 64 Appendix B 67 V TABLES Table I - A t t r i t i o n of F o u r i e r Components of G a u s s - S e i d e l R e l a x a t i o n of L a p l a c e ' s Equation i n the U n i t Square Table II - L i n e R e l a x a t i o n Table I I I - D i s t r i b u t i v e 39 R e l a x a t i o n of the Cauchy-Riemann Equations Table IV - Best Values of 38 co 41 50 vi FIGURES Figure 1 - Orderings: (a) L e x i c o g r a p h i c F i g u r e 2 - R e l a t i o n of F i n e (b) Red-Black.10 ( f ) and Coarse (C) Grids..13 F i g u r e 3 - Numbering of G r i d P o i n t s 24 F i g u r e 4 - E n t r i e s i n R1 Corresponding to G r i d P o i n t s on the Boundary F i g u r e 5 - Contour L i n e s of 32 |M(0)| f o r SOR Relaxation.48 ACKNOWLEDGMENT I would l i k e t o thank encouraging his Robert Miura f o r many t a l k s , even though my work was o u t s i d e f i e l d of r e s e a r c h , and U r i Ascher, who f i r s t exposed me t o M u l t i - G r i d , and Marc Lemieux, a very friendly office-mate. thank Mary Regester Most of a l l I would l i k e t o who put up with me during a l l the consummate aggravation of producing t h i s work. 1 CHAPTER ONE INTRODUCTION 1.1 The Subject Numerical differential application devoted techniques equations i n science to the algebraic are fast and robust problem equations by finding partial increasing research algorithms is for this problems, and step i n t h e i r s o l u t i o n i s the t r a n s l a t i o n by into the a large use of f i n i t e element d i s c r e t i z a t i o n . These solved of PDEs a r i s e from steady s t a t e first of the continuous (PDEs) solution and i n d u s t r y . Much current finding purpose. E l l i p t i c commonly f o r the sparse system of a f i n i t e d i f f e r e n c e or systems are frequently f i r s t approximating the s o l u t i o n very roughly, and then reducing the e r r o r i n t h i s aproximation by an iterative method ( r e l a x a t i o n scheme). In the M u l t i - G r i d used, not to e l i m i n a t e s o l u t i o n , but only relaxation (MG) method, r e l a x a t i o n techniques are the e r r o r i n an to smooth t h i s e r r o r . approximation (For a d e s c r i p t i o n of techniques and the MG method, the reader may r e f e r to Chapter Two.) Since the p r a c t i c a l i n t r o d u c t i o n method promise by t o the Achi and Brandt i s now in being 1971, of convergence of a MG a l g o r i t h m . of the efficiency of the the i t has shown great investigated c r u c i a l problem f o r a p p l i c a t i o n s of This smoothing early systematically. i s the e s t i m a t i o n requires MG A of the rate an performed estimation by the relaxat ion. Local introduced Mode by Analysis Brandt (LMA) or L o c a l F o u r i e r (see Brandt(1977)) to Analysis was realistically 2 estimate the convergence rate of the smoothing step i n a M u l t i - G r i d a l g o r i t h m . L i k e the von Neumann s t a b i l i t y for parabolic systems, LMA treats a analysis model system without boundaries, i . e . , on a g r i d which i s i n f i n i t e in extent. can of the be used to obtain good estimates convergence r a t e of the MG a l g o r i t h m as a not given whole. LMA (rapid) Brandt has a r i g o r o u s j u s t i f i c a t i o n of t h i s procedure (on t h i s t o p i c , see Brandt (1982)). Other mathematicians convergence LMA have tried develop rigourous e s t i m a t e s f o r MG methods, but they have not found u s e f u l f o r t h e i r t h e o r i e s . Most such rigourous convergence estimates have been too c o n s e r v a t i v e magnitude) the to (by several orders to be of any p r a c t i c a l v a l u e . T h e r e f o r e LMA practical tool for realistic estimation of remains and d e c i s i o n making with regard to MG methods. Despite i t s p r a c t i c a l success, the a p p l i c a b i l i t y to real problems disproven. It i s infinite finite LMA model of LMA has been n e i t h e r r i g o u r o u s l y j u s t i f i e d nor not even system clear ought to how be the analysis of an brought to bear on the system at hand. Thus f a r mathematical work related to has taken three d i r e c t i o n s . Trottenberg and Stueben (1982) c a r r i e d out a rigourous F o u r i e r a n a l y s i s of c e r t a i n s p e c i a l MG D i r i c h l e t problems relaxation situations i n a square using Red-Black (essentially or Richardson f o r smoothing). They d e r i v e d convergence r a t e s f o r the smoothing step which turned out to be identical convergence r a t e s d e r i v e d v i a an extension of LMA f o r the same problems. They claimed that the i d e n t i t y of these r e s u l t s credibility to the application of LMA in with more gave general 3 s e t t i n g s . A second avenue of the mathematical the elaboration of Brandt's of restriction infinite operators. (1981), and de V r i e s the grid Incomplete in detail analogues Further of Hemker L-U Decomposition intent for certain this relaxation low frequency with the operator used to This idea method, which was frequencies A on a the These rate extends two mesh values performance the original be realistic may direction be a n a l y s i s of a f i n i t e d i s c r e t e differential equation p e r i o d i c boundary then of the spacing were of was close, might be relaxation horizons only problem related (1982), interpreted problem in c o n d i t i o n s . These between He of the for high 1982). mathematical LMA these estimate radius whose mentioned by T r o t t e n b e r g and Stueben similarity of component. spectral grid envisaged to (1983). Mol (1982) went f u r t h e r than Fourier overall (see Brandt 1977, third Hackbusch use f o r f i n e meshes, which suggested that LMA estimate schemes. the convergence actual determined by the frequency. especially 1981), r e g a r d i n g the range of a p p l i c a b i l t y de V r i e s computed the LMA compared and problems. of LMA. a prolongation (ILU) r e l a x a t i o n p r o c e s s . i s worth n o t i n g that de V r i e s original Fourier (1982) analysed an i n f i n i t e g r i d model of estimates f o r f i n i t e g r i d Brandt's been infinite the (1980b, However, they d i d not r i g o u r o u s l y j u s t i f y It has o r i g i n a l n o t i o n of an g r i d model. Hemker (1980a) d e s c r i b e d analysis work a which LMA MO1(1981), as imply is and a rigourous arises rectangle authors to from a with doubly that the the p e r i o d i c problem and the (more usual) D i r i c h l e t problem j u s t i f i e s the use of LMA f o r the l a t t e r . 4 Hackbusch (1983) went Dirichlet singular be cases on explicitly. perturbation to He relate the treated problems i n one the MG the relaxation operator, and used this periodic to link operator to the norm related is further rate of relaxation size MG only author who Outline In this description D i r i c h l e t relaxation estimate operator estimate made use of of in the convergence LMA type analysis in a MG. from arbitrary domains. these problems are LMA of F o u r i e r components on is The relaxed i t e r a t i v e schemes to be d e s c r i b e d . interpret f i n e g r i d s . To actual He work I t r e a t two-dimensional D i r i c h l e t problems l a r g e c l a s s of relaxation matrix. of t h i s T h e s i s arising I the D i r i c h l e t high frequency components. T h i s on g r i d s which approximate Three, the equation special rigourous treatment of systems of to a l g o r i t h m for t h i s problem. Hackbusch i s the has 1.2 of a used i n a rigourous the an perturbation s i z e of the p e r i o d i c to the solution derived operator, a rank one the and dimension, which are smoothed by Gauss-Seidel r e l a x a t i o n . He involving periodic discrete by any In Chapter as an approximate d e s c r i p t i o n asymptotically the finite correct i n the of grids. to obtain a matrix which operates on domains, t h i s perturbed operator i s i d e n t i c a l with I LMA. operator for a problem with show small in norm, and that t h i s perturbation hence v a l i d a t e LMA on In the Fourier rectangular conditions. by the exactly relaxation described This components the as the l i m i t of very show t h i s , I p e r t u r b the matrix d e s c r i b i n g relaxation of a periodic is fine case of boundary asymptotically grids. This 5 validation detail, is to i l l u s t r a t e d with show the selected a c c u r a c y of Fourier relaxations over a range possibility convergence rate usefulness of parameter for of and efficiency the of LMA using the idealizations compared as LMA to LMA i s is d i d de the the MG in for actual illustrated also used in Vries (1982), the overall to context. optimum over performance a range Laplace's for The relaxation heuristically ADI method a p p l i e d t o degradation with describe outside predicting SOR scheme of LMA worked out p r o b l e m s and m e s h e s . investigate, in examples are relaxations p r o b l e m s and m e s h e s . the I of the LMA. components In C h a p t e r F o u r , the several of explain equation non-symmetric operators. Chapter this work, as I hope its Five possible is well this as a discussion a proposed work w i l l uses as of build some of the ideas behind i n LMA and suggest extension. confidence a theoretical tool. 6 CHAPTER TWO BACKGROUND 2.1 Finite Differential In Difference Solution (PDEs) is of e l l i p t i c frequently only m e t h o d s . Commonly t h e ( c o n t i n u u m ) PDE system of d i s c r e t e equations may be Elliptic solved either Partial Differential feasible by n u m e r i c a l is translated (discretization). by a d i r e c t in solution MG which there different has method o r by an iterative i s interaction proven quite t e c h n i q u e , L o c a l Mode A n a l y s i s to the study of i t e r a t i v e Consider elliptic only on a 0 which MG and a p p l i e d finite L. domain We fl will and an consider The p r o b l e m i s t o f i n d a f u n c t i o n u is a f(x,y), prescribed function on some p r e s c r i b e d l i n e a r fl and u and/or boundary its conditions 9 f l , t h e b o u n d a r y o f fl. Our as many methods. L. Lu = derivatives satisfy on treating The satisfies (2-1 ) f in itself. (LMA), i s j u s t i f i e d two-dimensional operators (MG) work one f a c e t o f t h e p a r t i a l d i f f e r e n t i a l operator linear where Multi-Grid between t h e i t e r a t i v e effective types of problems. In t h i s a equations of t h e e q u a t i o n s and t h e d i s c r e t i z a t i o n method into These a l g o r i t h m . A r e c e n t d e v e l o p m e n t h a s been t h e method Partial Equations practice the solution Equations of finite follows. difference Place a solution discrete of t h i s problem proceeds s e t o f N p o i n t s on fl, s p a c e d 7 evenly i n x and i n y to form a the same spacing rectangular grid. Frequently i s used i n both c o - o r d i n a t e d i r e c t i o n s and t h i s spacing i s c a l l e d the mesh s i z e h. There are s e v e r a l ways to approximate the boundary of a non-rectangular domain. We w i l l d i s c u s s one of them in Chapter Three. The g r i d p o i n t s are ordered so that numbers assigned to them form the components of a v e c t o r U which we c a l l a g r i d Now we approximate the continuum equation f o r system of U at vector. u by N a l g e b r a i c equations which determine the v a l u e s of the N grid points. The components of U are to approximate the v a l u e s of u on the corresponding g r i d We a c o n s t r u c t an equation at each i n t e r i o r g r i d p o i n t approximating the d e r i v a t i v e s that appear differences, points. (x,y) by i n Lu(x,y) by finite e.g., (2-2) 9 u / 9 x ( x , y ) = (1/h )(u(x+h,y)-2u(x,y)+u(x-h,y)) + 0 ( h ) . 2 2 2 2 These d i f f e r e n c e approximations to the d e r i v a t i v e s are used to build up an approximation L*U(x,y) equation c o r r e s p o n d i n g to a g r i d p o i n t (2-3) Thus L*U(x,y) = for each interior the grid components of points (x,y) U U (x,y) i s then point we the have an equation c o r r e s p o n d i n g to that p o i n t and c o r r e s p o n d i n g to neighbouring p o i n t s . At on or near the boundary approximate equations are set up, but now account discrete f(x,y). grid i n v o l v i n g the component of to L u ( x , y ) . The boundary conditions of fi, s i m i l a r they must take imposed on u. We into obtain a 8 linear N x N X = U system of approximates equations the AX solution approximation g e n e r a l l y becomes b e t t e r further information on = b, whose u of as N solution the PDE. The increases. For the d i s c r e t i z a t i o n process see Young (1971 ) . Now AX we .must = b. find an efficient algorithm T h i s i s the most c o s t l y part of the s o l u t i o n of the PDE. solving f i n i t e difference that are d e s i r a b l e to ensure accuracy, the usual Gauss e l i m i n a t i o n (GE) algorithm turns out For for the large systems (N > 10 ) 3 to be too space- and time-consuming. Other d i r e c t methods, based on the F o u r i e r transform (transformation of dense system), may the large sparse be used to obtain development practical of large core exact relaxation l a r g e systems can techniques). f i x e d point made and them = b obtains (2-4) h(X) has point any X°, grid X 1 if i t e r a t i v e methods (or iteration f o r the a matrix B which approximates A and i n v e r t e d . Then the a fixed quickly i t e r a t i o n s . It i s these methods with which AX iterates some A common c l a s s of r e l a x a t i o n schemes s h a l l be concerned. A f i x e d point Choose has be handled more we easily of f a s t e r than GE computers they are solved only approximately using the solutions smaller for the s o l u t i o n of l a r g e r systems. In general are reduction system i n t o a s p e c i a l d i s c r e t e systems. These methods are the or on + 1 X = U = h(X°), which i s function = X - B" AX vector system the X° X B" b 1 desired and 2 discrete construct = h(X ),... 1 Let solution. a sequence of e be the 9 difference between the d i f f e r e n c e an i t e r a t e e' X and U, e = X - U. between the next approximation X' Then and U i s given by (2-5) e' = (I - B" A)e. 1 If the s p e c t r a l radius of the a m p l i f i c a t i o n matrix is less than 1 then i t e r a t e s w i l l go to clearly 0 and the I - B"'A the e r r o r s i n the sequence sequence will converge of to X = U. There are simple s u f f i c i e n t c o n d i t i o n s which are known to guarantee convergence of such an i t e r a t i o n 1971). For example: l e t B" 1 and B - A A be s p l i t are non-negative into In the J a c o b i i s the d i a g o n a l of (G-S) method, B A = B - (B - A); i f (a " r e g u l a r " s p l i t t i n g ) , the i t e r a t i o n w i l l converge. Most of the i t e r a t i o n s produce scheme (see Young, common fixed (or simultaneous displacement) method, A (see Young, 1971). In the is the t r i a n g l e and one of the lower non-zero S u c c e s s i v e Over-Relaxation (SOR), diagonal d i a g o n a l of Vries, and L the A. Incomplete 1982) triangular i s upper some sense and such that part easily inverted. from PDEs L , and R is D + a>L of A i s the lower of where A. In D is lower t r i a n g l e not i n c l u d i n g the L, U, and R t r i a n g u l a r and A = LU + R; therefore U B B Gauss-Seidel superdiagonals L-U Decomposition finds matrices lower t r i a n g u l a r , U point r e g u l a r s p l i t t i n g s as o u t l i n e d below. ( i n c l u d i n g the d i a g o n a l ) . For l i n e r e l a x a t i o n , B the then For B (ILU) (see such that R is i s "small" in is sparse m a t r i c e s should be sparse. The L de LU and A is arising Alternating- 10 Direction splitting be Implicit iteration it further more on a simple later. The Conjugate-Gradient will (CG) implementation of a relaxation scheme than j u s t a knowledge of the s p l i t t i n g on which i s based. We will simplest widely used method w i l l be used The based i s not of the type d e s c r i b e d above. practical entails i s not but rather on a sequence of two s p l i t t i n g s . I t discussed The (ADI) method discuss the implementation r e l a x a t i o n , the Gauss-Seidel of the scheme. T h i s f r e q u e n t l y i n examples. g r i d p o i n t s are ordered i n some manner. In that order each approximate s o l u t i o n value i s a d j u s t e d t o s a t i s f y the one corresponding d i s c r e t e equation. Depending on the o r d e r i n g of the p o i n t s , the d e t e r m i n a t i o n of the value of the next iterate at a g r i d p o i n t may i n v o l v e using the newly determined values at one or more neighbouring p o i n t s as w e l l as values of the 13-14-15-16 I I I I 9-10-11-12 1 1 --5-1 5 — 8 I I I I I I I I I I I I 2-12--6-16 . 9 — 3-1 3 ~ 7 1-10--4-14 (a) Figure current and (b) 1. O r d e r i n g s : (a) L e x i c o g r a p h i c and (b) Red-Black i t e r a t e . The most common o r d e r i n g s red-black are l e x i c o g r a p h i c (see F i g u r e 1). L e x i c o g r a p h i c r e l a x a t i o n at one corner of the g r i d and proceeds along When relaxation starts at the end of the l i n e i s reached, the f i r s t points one p o s i t i o n of the next l i n e , and so grid starts line. again on. u n t i l a l l on a l l g r i d l i n e s have been m o d i f i e d . T h i s c o n s t i t u t e s 11 one i t e r a t i o n . Red-Black r e l a x a t i o n d i v i d e s into two ordered all grid points subsets r e l a t e d to each other l i k e the red squares and the black points the squares of a checkerboard. A l l the of the red subset are r e l a x e d i n t h e i r order, and then the p o i n t s of the black subset. 2.2 The M u l t i - G r i d The better (MG) Method best of the r e l a x a t i o n schemes than Gauss Elimination, i n e f f i c i e n t , and t h e r e f o r e first few passes approximate of solution a mentioned but still they costly. above are s t i l l quite Typically on r e l a x a t i o n scheme, the e r r o r substantially 50%). However, subsequent i n the (on the order of 10% to passes are not so e f f i c i e n t and long before the d e s i r e d t o l e r a n c e on the r e s i d u a l s i s reached, are being reduced by only a small f r a c t i o n on each many the (as measured by the s i z e of the r e s i d u a l s r = b - AX) i s decreased Very are relaxation sweeps they iteration. are thus r e q u i r e d to s o l v e the equations. For most of these spectral k 0(h ), radius of schemes the shows that the the a m p l i f i c a t i o n matrix i s of order 1 - where h i s the mesh s i z e . Depending on the scheme used and the s p e c i f i c problem, the power to theory 5/2. In approaches practice the k of h v a r i e s convergence the s p e c t r a l r a d i u s , s i n c e the rate from 3/2 asymptotically components of the e r r o r along the e i g e n v e c t o r s of the a m p l i f i c a t i o n matrix, with the largest e i g e n v a l u e s , g r a d u a l l y become dominant. That i s , the q u i c k l y converging e r r o r are almost eliminated by components the f i r s t (small eigenvalues) few i t e r a t i o n s , but the slowly converging components are worn down only g r a d u a l l y and 12 thus stand efficiently out after many passes. If we are solve the problem we need to f i n d some t o ' more other way to e l i m i n a t e these components. What do the slowly converging e r r o r components look The residuals after several r e l a x a t i o n scheme are smoother initial guess X ° . That applications than the of as one moves smooth. Numerical The converging computation strongly between changes error components f o r simple model problems bears out largest are of the e i g e n v a l u e s of the this e i g e n v a l u e s are g e n e r a l l y a s s o c i a t e d with e i g e n v e c t o r s , whose values at are sign the along a row of the g r i d . Thus i t i s n a t u r a l to a m p l i f i c a t i o n matrix intuition. any for i s , there i s more c o r r e l a t i o n c o n j e c t u r e that the slowly also nearly residuals the r e s i d u a l s at neighbouring p o i n t s , and fewer like? neighbouring grid points c o r r e l a t e d . On the other hand, the e i g e n v e c t o r s corresponding to small eigenvalues vary quickly along grid lines. These observations suggest that the r e l a x a t i o n might be b e t t e r understood by decomposing the e r r o r in a an iterate heuristic X, analysis into finite suggests, and Fourier-type this thesis schemes e = X - U series. makes A it r i g o u r o u s , that while most r e l a x a t i o n schemes are e f f i c i e n t at reducing the high frequency F o u r i e r components of e r r o r , are not e f f i c i e n t at reducing the Thus we cannot in general low frequency they components. e l i m i n a t e the slowly converging e r r o r components from one r e l a x a t i o n scheme simply by using different a scheme. The idea behind the MG a l g o r i t h m i s to use the r e l a x a t i o n methods to do only what they do w e l l - reducing the high 13 frequency components of the e r r o r , method then varying error This may way AX = b to eliminate i n an approximate be have obtained from as been some other of the - U 1 which s o l v e s Ae or Suppose the slowly several discrete r 1 = b - AU system discrete the are 1 smooth U error 1 Le = r smooth. solution = r . Consider Ae 1 d i s c r e t i z a t i o n of a PDE problem satisfy on exact e q u i v a l e n t to the d e t e r m i n a t i o n of = U "smooth" (2-1), y i e l d i n g an approximate s o l u t i o n 1 determination a follows. done U , and suppose that the r e s i d u a l s 1 find solution. approached r e l a x a t i o n sweeps e to to e l i m i n a t e the low frequency components. The problem i s to f i n d a The and 1 = r i n J2, where e vector 1 as the i s to homogeneous boundary c o n d i t i o n s . We might expect that e w i l l Lu = f . be smooth r e l a t i v e to the g r i d we are is then using f o r Then an a c c u r a t e s o l u t i o n f o r e may be obtained with a much c o a r s e r g r i d . Thus we set up a f i n i t e d i f f e r e n c e system A*e* = r * by d i s c r e t i z i n g Le = r on a c o a r s e r g r i d 2). The solution e* —-C — — c I of f I this C I I f f I I I c should function Since e f (C) G r i d s on the p o i n t s of the coarse g r i d to which i s approximated e i s smooth, v a l u e s of a good f F i g u r e 2. R e l a t i o n of Fine ( f ) and Coarse approximation be f I f f system (see F i g u r e the smooth on the f i n e g r i d by e on the f i n e g r i d points e . 1 may be obtained a c c u r a t e l y simply by i n t e r p o l a t i o n from the coarse 14 grid p o i n t s . Thus a good approximation for e from e*, and Why (CGC) t h i s smooth e r r o r can would t h i s procedure ) be more e f f i c i e n t can 1 be almost be obtained eliminated. ( c a l l e d Coarse Grid than simply e l i m i n a t i n g Correction e by many 1 more r e l a x a t i o n s on the o r i g i n a l f i n e g r i d ? It i s because each relaxation much of A*e* = r * r e l a x a t i o n of relaxation required than Ae = r 1 of AX 1 is (which is = b ) . In a d d i t i o n to a be required quickly this Of accomplished the to solve the in the solution same way; t r a n s f e r to a yet coarser i s devised. i s reduced by resolution further tolerance f i n e g r i d system to can be be eliminated the more way. course algorithm than a fewer i t e r a t i o n s w i l l to solve the coarse g r i d system to a given will coarser expensive equivalent same t o l e r a n c e . Thus a smooth e r r o r error less of On of a few CGC problem may r e l a x a t i o n sweeps and be then g r i d . Thus a r e c u r s i v e m u l t i p l e - g r i d each g r i d the r e l a x a t i o n and that the grid) the rough component of the smooth ( r e l a t i v e to the component i s reduced v i a the next grid. These ideas may algorithm be implemented in several ways. sketched below i s r e f e r r e d to as an a d a p t i v e c o r r e c t i o n scheme (see Brandt, 1977, 1982) Let 1 2 successively finer g r i d s be G , G ,... ,G n s o l u t i o n i s a g r i d v e c t o r on G . Let the 1 n m m+1 h > ... > h ; u s u a l l y h = 2h . There are transform g r i d v e c t o r s the n V-cycle sequence of ; the mesh The desired sizes be o p e r a t o r s which on one g r i d to g r i d v e c t o r s on another k k-1 grid. Prolongation I acts on g r i d v e c t o r s on G and k-1 smoothly i n t e r p o l a t e s to y i e l d v a l u e s at a l l the grid points k k-1 in G . The reverse procedure, r e s t r i c t i o n I , interpolates k 15 values at g r i d points values of a g r i d vector representation on solve a follows. in the The the has last using k the data s u p p l i e d G . Both operations target on grid of large find the o r i g i n a l g r i d . k k k A X = b on each problem a i s determined by measuring the components, reduce the a r e l a x a t i o n i s no grid G The operator k as l e s s than rate of decrease long as the relaxation X error should K s i z e of the times the K residuals size smooth. T y p i c a l l y l e s s before, than 5 done. After these r e l a x a t i o n sweeps we set up the CGC k-1 k k k k k-1 k-1 k on G . Let r = b - A Z , and b =1 r . The k-1 k to be solved using g r i d G i s then A good efficiency criterion soon as the then the e r r o r i s judged to be i t e r a t i o n s are each r e s i d u a l s . An (0 < K < 1) must be s e t . As (2-6) the s e v e r a l times u n t i l the e r r o r k solution Z i s judged to be approximate non-smooth k-1 by a f u n c t i o n which i s r e s i d u a l s in some s u i t a b l e norm. As efficiently after k-1 system i s r e l a x e d "smooth". T h i s of on the represented by a vector We in G k-1 A = b problem system k-1 k-1 used for G f i n i t e d i f f e r e n c e approximation on k-1 G may k-1 be taken to be to the the underlying differential operator L , or, (more s u i t a b l y f o r t h e o r e t i c a l k-1 k-1 k k purposes), as A =1 A I k-1 k k-1 On G (2-6) may be solved either by GE or by a sufficient number of r e l a x a t i o n sweeps. Using e i t h e r method the amount of work done i s small on the c o a r s e s t other grids (2-6) i s solved by g r i d . On smoothing steps and a CGC. the It 16 i s only necessary commensurate k on G to solve the with CGC problem to an accuracy the s i z e of the remaining non-smooth e r r o r k-1 k-1 k-1 When a s o l u t i o n X = U i s obtained on G then l e t k k k-1 k k e = 1 U . An improved s o l u t i o n Y on G i s obtained v i a k k k-1 k Y = Z e. T h i s e l i m i n a t e s most of the smooth part of the k k k e r r o r i n Z . Since the e r r o r Y - U i s now dominated by i t s non-smooth reduced component, this e r r o r may be f u r t h e r efficiently by r e l a x a t i o n s t e p s . Thus the original non-smooth component n approximation on G of the error has been cut r e l a x a t i o n s , and the smooth part has been almost i n the down by e l i m i n a t e d by the CGC. and The M u l t i - G r i d procedure outlined although seems the MG idea c o n t r o v e r s y among those who implemented 1977) and the scheme suggested converging) in that algorithm have above is promising, there i s s t i l l applied i t . Brandt i t was an f o r numerical efficient solution problems has l e d to much (rapidly of PDEs. He problems were s o l v e d e f f i c i e n t l y v i a MG. However, slow convergence other first the e a r l y s e v e n t i e s (see Brandt, promoted i t v i g o r o u s l y and many d i f f i c u l t for complicated of doubt the method regarding i t s robustness. There i s l i t t l e poor performance inappropriate fundamental the MG Bakhvalov theory which can h e l p decide whether of some MG implementation, limitation algorithm a p p l i c a t i o n s i s the r e s u l t of or whether it reflects a of the method. Convergence p r o o f s f o r (given (see Brandt, the already in 1977)) show that the a e a r l y s i x t i e s by reduction of one 17 order 0(N) of magnitude in the o v e r a l l e r r o r can o p e r a t i o n s (N i s the constant in t h i s 0(N) convergence rates Brandt's quickly claim local number of g r i d p o i n t s ) . However i s so l a r g e as to that on the Fourier prolongations, this we are bearing the algorithm of will but to relaxation a given problem and local estimate the r e l a x a t i o n at one relaxation (see r e l a x a t i o n a c t i n g on an rate w i l l converge will Local of r e d u c t i o n grid Brandt, point 1977). infinite illustrate LMA of the strongly the new Mode values X' Straightforward use convergence i n 0(N ) operat i o n s . 2 of (or various the i.e., "relaxation is a Thus g r i d and he imagined asked how a relaxation component. by considering Gauss-Seidel to the o l d values X during over l e x i c o g r a p h i c a l l y ordered g r i d p o i n t s 1 LMA affects r e l a x a t i o n of Laplace's equation i n some domain. The linking each e r r o r . Brandt's o r i g i n a l idea would reduce a p a r t i c u l a r F o u r i e r We and r e l a x a t i o n . In t h i s t h e s i s only at nearby g r i d p o i n t s , process" also a s p e c i f i e d r e l a x a t i o n scheme modes which comprise the that (and a relaxations. For to on non-smooth e r r o r components on L o c a l Mode A n a l y s i s used converge rather suggests that MG s u b s t a n t i a l l y reduced by (LMA) on With s u i t a b l e r e s t r i c t i o n s analysis 2.3 Fourier was no w i l l be concerned with the a p p l i c a t i o n of t h i s Analysis is have theory analysis prolongation). r a p i d l y provided that MG rigorous r e s t r i c t i o n and grid the 1 in p r a c t i c e . i s based not heuristic be achieved in equation a sweep is a good r e l a x a t i o n w i l l r e s u l t i n possibly N to the power 3/2) 18 (2-7) 4X'(i,j) -X'(i-1,j) -X'(i,j-1) Note that when the p o i n t values The X' at p o i n t s desired we newly computed ( i , j - l ) are a l r e a d y obtain the error transformation (1977) suggested that interior Fourier of the g r i d by s e r i e s . ... Thus component we ( i - 1 , j ) and reached, 0. in place. and upon equation, 4 e ' ( i , j ) - e ' ( i - 1 , j ) - e ' ( i , j - 1 ) -e(i+1,j) - e ( i , j + l ) = Brandt the is s o l u t i o n U a l s o s a t i s f i e s t h i s equation subtracting (2-8) (i,j) -X(i+1,j) -X(i,j+1) = of the to one can 0. "analyse i t [2-8] in ( l o c a l l y ) expanding the e r r o r in study e r r o r before and the d = (d*,d ) Fourier 2 a f t e r the relaxation sweep, put e(i,j) = A(0)exp(i(0 i+0 j)) (2-9) 1 2 (i = /-I) and e ' ( i , j ) = A*(0)exp(t(0 i +0 j ) ) . " 1 Here 0 1 and y and 8 2 must be directions, 2 interpreted as wave numbers respectively. exp( c0(k±1 )) = exp( t0k)exp(± it9) the equation (2-8) (2-10) (4 - e x p ( - i 0 ) - e x p ( - i 0 ) 1 1 the convergence ratio the x Using the identity error transformation becomes + (-expd.0 ) The in 2 )A'(0) - e x p ( t 0 ) )A(0) r a t e of the 2 9\8 2 = 0. component i s then given by 19 (2-11) M (0) = A'(0)/A(0) = (- exp( <,en exp( id ) - ) 2 / ( 4 - exp(-i0 ) - exp(-i0 ) 1 Brandt expected good for interest t h i s estimate the high i n the MG smoothing of the convergence frequencies context only), ). 2 role for boundary conditions and the "high mesh s i z e frequencies" 1 2 and indistinguishable from values second at every c o n s i d e r e d . The estimate of components (2-12) "low interior i n the since the 1 work estimates the with grid is 2 point in of the the ; max ( | 9 | , | 9 2 | ) > the accurately for a the MG is only and are an frequency as TT/2}. solution. high if fine grid high a l g o r i t h m as a whole, and required as component 2 i s u s u a l l y ignored rate 1 taken component -TT, 9 -IT) 1 grid convergence (0 ,0 ) the (9 is r a t e of smoothing can be used to of point finer frequencies" form of u.(9) u = max{|/z(0)| efficiency is frequencies ( c a l l e d the smoothing f a c t o r ) i s taken This e f f e c t i v e total exact the relaxation i n the c o a r s e r g r i d , the boundary between \ ) = 7r/2 , max ( | 9 | , | 9 be with f a r away p o i n t s are s m a l l . For the usual case where the mesh s i z e half of these i n t e r a c t i o n s of the s o l u t i o n value at an the to (which are the only ones of where the because rate estimate to estimate In efficiency practice of many the the it MG implementat i o n s . Despite have worked i t s p r a c t i c a l success on MG methods have so far, found many the people LMA who procedure 20 dubious. Brandt justification himself for has never come forth with a i t or even a statement about when i t might be reasonably expected to work, and when i t might not. 21 CHAPTER THREE A JUSTIFICATION OF LOCAL MODE ANALYSIS In t h i s chapter we Mode Analysis (LMA) w i l l consider for the problems with c o n d i t i o n s . We w i l l take the p o i n t of approximately describe the e r r o r under r e l a x a t i o n . We becomes better concrete as examples will how view show that that the Local boundary LMA should of a F o u r i e r mode this i s r e f i n e d . We good of Dirichlet transformation the g r i d of validity approximation will give approximation several is in practical situations. 1977) 3.1 The Meaning of L o c a l Mode A n a l y s i s . In explanations of Local series" arises case of ( i , j ) . T h i s notion Brandt, by boundary ignoring the v a r i a t i o n of c o e f f i c i e n t s . However we notion rigourous, We of "local Fourier and variable c o n d i t i o n s and must . f i r s t dispose s e r i e s " in order first d i s c u s s t h i s in the context l a t e r we will discuss of a one-dimensional the multi-dimensional case by analogy. Consider a f u n c t i o n f ( x ) d e f i n e d r e q u i r e d of a " l o c a l F o u r i e r s e r i e s " (3-1) of to make a some i n t e r v a l I c o n t a i n i n g a p o i n t x = a. Consider what be the g l o b a l , i n t e r p r e t a t i o n of the a n a l y s i s . continuum, and discrete a Fourier look p l a u s i b l e i n conditions boundary Fourier of " l o c a l from an attempt to make LMA general coefficients, the (e.g., the e r r o r i s supposed to be w r i t t e n as a " l o c a l s e r i e s " near a g r i d point the Mode A n a l y s i s 2Ld(0)exp(i0(x-a)) e on would 22 for f ( x ) near x = a. We might r e q u i r e i t t o converge i n a neighbourhood of x = a and to close these to a, 'local' perhaps describe the behaviour of f converging at a to the value f ( a ) . But requirements a c t u a l l y come from the idea of a T a y l o r s e r i e s . Furthermore, how would the c o e f f i c i e n t s d(0) be determined? If we determine them over the neighbourhood (a~c,a+c) of x = a, then the c o e f f i c i e n t of the 0 component i s given by ma+c * a+c (3-2) d(0) = l e x p ( - 0 i x ) £ ( x ) d x J a-c Note /J |exp(-10x)| dx. 2 / ' a-c that (3-3) l i m d(0) = exp(-i0a)f(a) c-K) which the i s independent of the behaviour of f ( x ) near a, so l i m i t of d(0) as c goes to 0 i s not a d e t e r m i n a t i o n of the local importance of the 0-component. We must then pick some c > 0, but which one? T h i s b r i n g s problem. and that an A us to the heart of the F o u r i e r s e r i e s can be d e f i n e d only f o r a f u n c t i o n interval, not a function and an indefinite neighbourhood of a p o i n t . In the case of a d i s c r e t e p o i n t set i n r e a l N-space, e.g., a g r i d c o v e r i n g a domain , c o n s i d e r a g r i d v e c t o r v 1 a specific point a. and The c o e f f i c i e n t of a F o u r i e r mode i n a The d i s c r e t e F o u r i e r modes are then v e c t o r s with components equal to the continuum F o u r i e r f u n c t i o n s e v a l u a t e d at points of t h i s ordered s e t . The inner product < , > on such a g r i d i s pointwise complex conjugate m u l t i p l i c a t i o n of two v e c t o r s with summation over a l l g r i d p o i n t s . <v,w> = £v(i,j)w(i,j) 1 23 " l o c a l expansion" of v near are considered neighbourhood to form the neighbourhood i s not w e l l - d e f i n e d . d i s c r e t e frequencies discrete x = a w i l l depend on which p o i n t s Fourier are to be series s e l e c t e d f o r the can a. This It i s a l s o not c l e a r which only be d e f i n e d g r i d vector corresponding to a s p e c i f i c the of expansion. f o r values rectangular A of a subset of grid. There is then, s e r i e s " i n the usual mathematics. The series. s e r i e s " , and language unique way sense that to form a " l o c a l ' l o c a l ' o b j e c t s are d e f i n e d i n idea of a While Fourier series the e x i s t e n c e and of " l o c a l F o u r i e r modes", We suggest that LMA that of such a " l o c a l is of a Fourier suggested i n some of the papers w r i t t e n on MG misleading. Fourier concept of " l o c a l F o u r i e r s e r i e s " i s r e a l l y a c o n f l a t i o n of the Taylor no by the methods, t h i s i s be r e i n t e r p r e t e d as a local a n a l y s i s of the a c t i o n of a r e l a x a t i o n scheme on F o u r i e r modes over the whole g r i d . We e r r o r s having the grid) is will show that the a t t r i t i o n form of a F o u r i e r component accurately computed by LMA i n t e r e s t of these s p e c i a l forms of e r r o r w i l l be used to give an a c c e p t a b l e 3.2 We the whole f o r f i n e meshes. is only that i n t e r p r e t a t i o n of The they LMA. Accuracy of L o c a l Mode A n a l y s i s develop an following elliptic (over r a t e of i n t e r p r e t a t i o n of LMA situation. PDE problem We consider in two c o e f f i c i e n t s of the form (3-4) + b3 u/9y Lu = a 9 u / 9 x 2 2 2 2 a and justify second dimensions i t in the order with linear constant + c3u/3x + d9u/9y + eu = f 24 in a domain fl with smooth boundary assumed to s a t i s f y D i r i c h l e t G covers fl and without We and the to l o s s of g e n e r a l i t y we assume that the i n the y direction c o n s i d e r zero-order approximations boundary approximated by line solution u is boundary c o n d i t i o n s on 3fl. A g r i d g r i d spacings i n the x d i r e c t i o n and equal. 3fl. The conditions. segments That to the boundary is, along the grid domain lines, D i r i c h l e t data are assigned to g r i d p o i n t s on these G r i d p o i n t s w i l l be doubly the grid were part of a indexed larger (i,j) (in and Figure is the lines. 2D) as rectangular grid p a r a l l e l to the c o - o r d i n a t e axes (see are 3). if running The first (3,5)--(4,5) (2^4>^T3)4)--( 4,4) ( 1 ,3)<f2,3)--(3,3)--(4,3) (2^>K^3,2)--(4,2) (3,1)--(4, 1 ) F i g u r e 3 - Numbering of G r i d P o i n t s index The i numbers the v a l u e s of x and equations are doubly so t h a t equation points are LMA (i,j), in ( i , j) involves (i-1,j), the g r i d indexed (i+1,j), j numbers the values of y. ( i , j ) l i k e the g r i d p o i n t s , values (i,j-l), of the (i,j+l) solution at ( i f a l l these i n t e r i o r ) . T h i s n o t a t i o n makes d i s c u s s i o n of easier. If the usual centered d i f f e r e n c e s are used to approximate the d e r i v a t i v e s i n t h i s PDE, and then only the values at i t s four nearest neighbours equation corresponding to a point are i n v o l v e d i n the d i f f e r e n c e that point. Thus a " f i v e point 25 formula" may be used i n the i n t e r i o r of fl. A system of d i f f e r e n c e equations i s obtained, corresponds to finite AX = f , where each equation a g r i d point.. In case the equation at a p o i n t i n v o l v e s boundary p o i n t s , (at most t h r e e ) , the f i x e d values a t those p o i n t s are to be i n c l u d e d i n the g r i d v e c t o r f . We will suppose that the equations have been m u l t i p l i e d through by h , 2 where h i s the mesh s i z e . We will justify which meet the f o l l o w i n g must be conditions. matrix A i s A = L - R with L n o n - s i n g u l a r . X' from the c u r r e n t This splitting iterate need not be r e g u l a r i t need not be the same at each only First AX = f of SPLITTING TYPE, i . e . , the c o e f f i c i e n t LX' = RX + f . will r e l a x a t i o n s of relaxations (3-5) and for linear the r e w r i t t e n as iterate LMA discuss splittings x We f i n d the next by s o l v i n g (L~ 1 and R non-negative), iteration. Secondly, f o r which each s p l i t i n v o l v e s only the same f i v e p o i n t s as the original we equation equation. T h i s means (3-6) a(m,n) = 0 Finally, i n order that every l i n e i n => l(m,n) = r(m,n) = 0. that LMA i t s e l f LX' = RX + f be f e a s i b l e we must require which corresponds to a p o i n t which i s not next to a boundary p o i n t i s of the form 26 (3-7) l X'(i,j) + l X'(i-1,j) + l X'(i,j-1) 1 2 3 + 1"X'(i+1,j) + 1 X'(i,j+1) 5 = r'Xli,]) + r X(i-1,j) + r X ( i , j - 1 ) 2 3 + r'X(i+l,j) + r X ( i , j + 1 ) + f(i,j), 5 where 1 , .. . , 1 , r , . . . , r 1 5 condition w i l l Many 1 linear relaxation R of two steps interchanged. homogeneity coefficients. schemes including Jacobi, Gauss-Seidel, of these, and ADI (see Chapter consists constant be subsequently r e f e r r e d to as 'homogeneity'. common conditions, are 5 Two). SOR, l i n e One ILU condition method and does cannot, these versions iteration such as are considered The satisfy of ADI here, with L and not satisfy strictly the speaking, be analysed by LMA. In p r a c t i c e an i d e a l i z a t i o n of the ILU method i s analysed t h i s way, but "we Red-Black and conditions Zebra and i t i s will not relaxations known justify do not (S. McCormick that here. satisfy these and K. Stueben, p r i v a t e communication) that LMA a p p l i e d simple-mindedly y i e l d s misleading estimates of smoothing r a t e s In j u s t i f y i n g LMA we w i l l c o n s i d e r from one iterate (see n o t a t i o n (3-8) It X to the next of e r r o r i t e r a t e X' which i s given by Le' = Re. a Fourier f i n i t e Fourier values the r e d u c t i o n i n Chapter Two) i s too complicated to c o n s i d e r write i n these c a s e s . of series transform), the Fourier for an a r b i t r a r y e r r o r it and (which amounts t o t a k i n g a and then ask what transform e, under Rather we suppose that the e r r o r c o n s i s t s happens to the a relaxation of a step. multiple of This 27 only one component \p(8) , as d e f i n e d below, and Fourier that a f t e r a r e l a x a t i o n step the new error is show U(8)\}J{8) plus something s m a l l . On a g r i d with the index convention, value of a discrete (with f r e q u e n c i e s = 8 point ( i , j ) i s given (3-9) .//(0)(i,j) = We 1 ^(8) dimensions), -n + 1 for grid 2 < 8 ,8 <n. Chapter Two, } in 2 as an eigenvector certain a j0 )}. LMA t r e a t s the F o u r i e r of the a m p l i f i c a t i o n matrix of the r e l a x a t i o n . As noted in Chapter One, case at by exp{i(i0 mentioned component 2 \p(8) F o u r i e r component, i n two (0 ,0 ), need only c o n s i d e r As exponential the relaxations of this is problems only with the periodic boundary c o n d i t i o n s on a r e c t a n g l e . However, as w i l l be for many relaxations, the Fourier components " c l o s e to being e i g e n v e c t o r s " , because they of that an idealized r e l a x a t i o n . The the actual relaxation is except at are i n f a c t are eigenvectors " c l o s e " to the a c t u a l idealized relaxation will relaxation seen, be identical with p o i n t s next to a boundary p o i n t . Since the p r o p o r t i o n of such p o i n t s i n the g r i d becomes s m a l l e r as the mesh becomes f i n e r , gives limit will we a good of f i n e g r i d s . The idealized relaxation to the a c t u a l r e l a x a t i o n i n the estimate of the accuracy of LMA which be developed makes such arguments r i g o u r o u s . For a s p e c i f i c F o u r i e r component, \p(8) will = R*\p(8), construct L\//' = R \ M # ) , and approximation the R*. f o r which First note a system 4/(8) L*\p* given by similar i s an e i g e n v e c t o r that every line in L\p' (3-9), of = both R\p(8) to L* which 28 corresponds to a point which i s not next to a boundary point i s of the form (3-10) lV(i,j) + lV(i"1,j) + 1 V>' ( i , j - 1 ) 3 + I V (i + 1 , j) + I V ( i , j + 1 ) = r V 0 ) ( i , j ) + r V S ) (i-1 , j) + r V + rV 0 ) (cf. ( 3 - 7 ) ) . Using (i + 1 , j) + r V 0 ) + r exp(-i0 ) 2 2 3 + r"expU0 ) + r exp(i0 ) 1 2 5 )tf(e)(i,j) r*i//(e)(i,j). = Now (i, j+ 1 ) , (3-9) the r i g h t hand side s i m p l i f i e s to (3-11) ( r + r e x p ( - i 0 M 1 6) ( i , j-1 ) if \p(8) i s s u b s t i t u t e d f o r \p' on the l e f t s i d e of (3-10) then i t s i m p l i f i e s to (3-12) ( l + l e x p ( - i 0 ) 1 + l exp(i0 ) 1 a = Now correspond point to (i, j), points. 2 3 + l exp(i0 ) 5 2 )<//(0)(i,j) R* so that values points like values We c o n s t r u c t (i,j) of [ R*\//( d) ] ( i , j ) which ( i , j ) next to a boundary are equal to L* [L*\p(8) i n c l u d i n g an a p p r o p r i a t e for l exp(-i0 ) l*^(0)(i,j). we c o n s t r u c t r*i//( 8) + 1 2 corresponding to other l i k e w i s e so that f o r every ] ( i , j) = l * ^ ( 0 ) ( i , j ) . m u l t i p l e of We exp(±i# ) 1 every equation i n v o l v i n g a boundary p o i n t , to do or interior interior this by exp(±t# ), 2 29 make i t l i k e the equation f o r p o i n t s away from the Specifically, left suppose edge p o i n t are 2 in L*\J/* = R*\p(8) If we , to 1 respect i v e l y . 2 c o n s t r u c t i n g an equation a t a 1 \^' (i-1 , j ) and a 2 term m i s s i n g i n (3-10). Then to form the ( i , j ) r e x p ( - c 9 ) \//( 8) ( i , j ) 2 are ( i , j ) where a term r \//( 8) ( i-1 , j ) equation we boundary . add l exp(-id the left 2 ) i//* 1 and ( i, j) and right sides, 3 denote these L1 = L* - L and extra elements R1 = R* - R, then along LI the d i a g o n a l by is 0 except for d i a g o n a l e n t r i e s corresponding t o boundary p o i n t s . LEMMA of 1: the system The F o u r i e r component ^L*e = R*e in L*\jj(8) At a p o i n t next L* are chosen Similarly, the interior L\p{8) M(#). grid point the ( i , j ) are both equal to ( 6) (i,j) . to a boundary p o i n t the a d d i t i o n a l e n t r i e s to make L*\//( 6) ( i , j ) equal t o f o r i n t e r i o r p o i n t s and p o i n t s near ( i , j ) entry i n If and i s an e i g e n v e c t o r with eigenvalue PROOF: I f ( i , j ) i s an entries ip{6) R*4/{8) i s the l * i / > ( 8) in (i, j) . boundary, r*tf(0)(i,j). we do a LMA and put I f Q i s a r e c t a n g l e , and the values of 0 , 8 are r e s t r i c t e d to those v a l u e s that would occur i n a F o u r i e r transform, then the m a t r i c e s L* and R* d e s c r i b e the r e l a x a t i o n of a problem with periodic boundary c o n d i t i o n s . However the c o n s t r u c t i o n used here i s not r e s t r i c t e d t o r e c t a n g u l a r domains. E.g., f o r the G a u s s - S e i d e l relaxation of Laplace's equation,the l e f t edge equation which has the form 2 1 2 3 4 e ' ( i , j ) - e ' ( i , j - 1 ) = e(i+1,j) + e ( i , j + l ) , (4 - e x p ( - i f l ) ) e * ( i , j ) - e * ( i , j - l ) = e ( i + 1,j) + e ( i , j + 1). 1 becomes, 30 v = A(0)exp( i ( 0 i (3-13) 2 = A'(0)exp(i(0 i v' where + 0 j ) ) , and 1 v and v' 0 j)), + 1 2 are the e r r o r s before and a f t e r r e l a x a t i o n , then (3-14) A ' ( 0 ) ( 1 + 1 l 2 e x p ( - 0 t + l"exp(i0 ) 1 = A(0)(r + + 1 r 2 1 + ) l 3 e x p ( - 0 t + l exp(i0 ) 5 exp(-i0 r"exp(i0 ) 1 + r 5 ) + r exp(i0 2 3 ) )exp(i(0 i 1 2 1 2 exp(-i0 2 + 2 j ) ) ) )exp(i(0 i + ) 0 1 0 j)). 2 Thus, (3-15) and /x(0) L*M(0)</>(0) Now let = r*/l* L*\p* = us = consider QED the difference between R*tf(0) the g r i d v e c t o r (3-18) R*\p{9). which s o l v e s M ( 0 ) * M 0 ) = (3-17) and A ' ( 0 ) / A ( 0 ) therefore (3-16) i//* = LiT = which i s the a c t u a l s o l u t i o n of R\p(6) . The d i f f e r e n c e between L and L* occurs only at p o i n t s next to a boundary p o i n t . For t h i s d i f f e r e n c e to be small then, we need a lemma t o guarantee a small number of g r i d to points next a boundary p o i n t . LEMMA 2: L e t the domain 0 be f i n i t e with piecewise boundary 3S2. Then f o r h s u f f i c i e n t l y small there smooth i s a number a, such that f o r any g r i d of mesh s i z e h, the number N(b) of 31 g r i d p o i n t s i n the i n t e r i o r of fl, one of neighbours lies whose four nearest on or o u t s i d e 9fl, i s bounded by a/N, where N i s the t o t a l number of g r i d p o i n t s i n i n t f i . PROOF: T h i s i s i n t u i t i v e l y subtle points in a details see Appendix A. obvious rigourous although there j u s t i f i c a t i o n . For For ease of a n a l y s i s we w i l l consider only are technical relaxations where, i n the equation corresponding to an i n t e r i o r g r i d p o i n t (i,j), the values at each of the neighbouring points ( i - 1 , j ) , ( i , j-1) , ( i + 1 , j ) , ( i , j + 1) occur e i t h e r on the l e f t or on the r i g h t s i d e of t h i s equation, but not both. Formally t h i s says that only one ( 1 , r ),...,(1 , r ) 2 2 5 side 5 of the members of each of the pairs i s not zero. A l l of the common r e l a x a t i o n schemes s a t i s f y t h i s c o n d i t i o n . Now equations we are ready to (3-17) and THEOREM consider the relationship (3-18). 1: If 90 i s smooth and the induced 2-norm of i s bounded independently of the number of g r i d p o i n t s there i s a constant A independent (3-19) - u(d),p(6)\ < PROOF: From (3-16) we (3-20) R\l/(6) = Thus from between = of the g r i d such that A|i//(0)|/VN f i n d that (R* - R1 )i/>(0) u(6) (3-18) (L + L1 )\jj(e) - R1^(0) . N, L~ 1 then 3 2 (3-21 ) = L " R\p{9) y = L " (ix(d)L^(d) = fi(d)4/(6) T h i s equation 1 + L " ( M ( 0 ) L 1 V ( 0 ) - RI<//(0) ). 1 t shows that v//(0) i s "almost an eigenvector of the i f the term L~ ( M ( 0) L11// ( 0) - R1iM0) ) i s of relaxation" norm. ) - L- R1i// + u(d)L1\p(6) 1 1 We w i l l small show t h i s by showing that L1\J/(0) and R1<H0) a r e of small norm. Consider the e n t r i e s of R1i/>(0). T h i s g r i d v e c t o r w i l l be zero at g r i d p o i n t s not next t o a boundary. Corresponding t o a grid p o i n t e x a c t l y one of whose four nearest the boundary, the v e c t o r values: For r exp(-i0 ), 2 have one of 3 e.g., with obtains At other 2 1 3 2 neighbours 2 + r exp(i0 ) a 1 p o i n t s and a t p o i n t s on the (see F i g u r e 4 ) . with __ 0 0 r*exp( i 0 ) __ 0 0 r"exp( t 0 ) —0 5 p o i n t s t o the r i g h t and on the bottom, the entry r e x p ( - i 0 ) 'corner' the f o l l o w i n g r e x p ( - i 0 ) , r * e x p ( i 0 ) , or r e x p ( i 0 ) . 1 the case of a point with two nearest boundary, one will neighbours i s on three nearest 1 1 r exp(-i0 )—r exp(-c0 )+r"exp(10 ) 3 2 3 2 1 I Figure neighbours on 4 - E n t r i e s i n R1 Corresponding t o G r i d P o i n t s on the Boundary the boundary, the entry w i l l sum of two (or three) of the preceeding be similar. with one neighbouring L e t N(e), N(c) (edges), s i m i l a r l y be the numbers. L1\^(0) will and N(w) be the number of p o i n t s two (corners), boundary points, and three (wedges) respectively. Let 33 a = m a x ( | r | , . . . , | r | , | l | , . . . , | l | ) . Then i n the g r i d 0 2 5 2 M(0)L1I//( 0) - R1\//(0) there are N(e) e n t r i e s of modulus by a 0 vector 5 ( s i n c e |ju(0)| ^ 1), bounded N(c) e n t r i e s of modulus bounded by 2a°, and N(w) e n t r i e s of modulus no larger than 3a°. Thus t a k i n g the 2-norm by summing the squares of the e n t r i e s , (3-22) | M(0)L1<//(0) " R1i/>(0) | £ a V ( N ( e ) + 4N(c) + 9N(w)). From Lemma 2 there i s a number a such that f o r any g r i d with h small enough, (3-23) N(e) + N(c) +N(w) < a/N. Thus from (3-20) and (3-22) (3-24) - M ( 0 ) ^ ( 0 ) | ^ | L - | a ° V N 3v/a. 1 If N i s the number of i n t e r i o r g r i d p o i n t s then the \p(8) , of 3|L | a _1 0 norm each of whose e n t r i e s has modulus one, i s i/N. Since /a i s bounded independently of N the r e s u l t follows. QED T h i s i s a bound on the r e l a t i v e norm between an a c t u a l r e l a x a t i o n and the LMA We may difference, also if sufficiently for large. domains. Obviously rotated the b(e) > 0, T h i s occurs 5 lower bounds on N(e) ^ b(e)/N domains this for N rectangular which have been the argument won't go through. minimum | r 1, ...,|r |,11 |,...,11 |. 2 difference idealization. f o r c i r c u l a r and for rectangular f o l l o w i n g simple be some simple the 7r/4 r a d i a n s there are no edge p o i n t s so that through Let b° obtain of 2 5 of the Then non-zero the coefficients grid . vector 34 u(8)L]\p(6) R1\/>(0) has - at l e a s t N(e) |M(0)|b° > 0 ( s i n c e each of the involved N(e) in the e n t r i e s of s i z e at l e a s t four relaxation at nearest a neighbours point). Thus is since ^ b(e)/N (3-25) | (0)L1<M0) - R1<M0)| > M > |y(0)|bVN(e) \u(6) |bVb(e) VN. Now (3-26) min |L~ x|/|x| Thus from (3-20) and (3-27) |*« In = (max 1 - ( 0 ) ^ ( 0 ) | ^ (1/VN) ( realize d i f f e r e n c e between the and ii, the that (3-24) the norm of norm of t a b l e s . It i s r a t h e r two which the on small calculations that difference - whose to first grids. usually u(8)\p(8), frequency l i e s whose second frequency l i e s will be seen means from projection in the the on is the the closer in the d i f f e r e n c e between which seems largest in relaxing components during It is i s NOT T h i s d i f f e r e n c e may \/>, the it These norms may. be as considered. "feeds" other F o u r i e r especially of the norm of the vector is (3-27) above result u(6)\}/. s i g n i f i c a n t component normal component and estimate than t h i s estimate would suggest, these =1/|L|. |bVb(e)/|L| ) \*(6) | .QED M to 1 (3-26), interpreting inequalities important vector |Lx|/|x|)" have a that relaxation, some of this sample this Fourier vector components range ( 6 - 0 ( h ) , d + 0 ( h ) ) and 1 range 1 (0 -O(h),0 +O(h)). 2 2 35 (0(h) means of bounded by ch where c is a positive constant.) Except later, f o r simple cases, some of which w i l l the estimates (3-24) and (3-27) be above cumbersome to be of much use i n p r a c t i c e . However be used to justify LMA asymptotically." discussed are too (3-24) This can g i v e s some p s y c h o l o g i c a l grounds f o r confidence i n i t s p r a c t i c a l use on f i n e g r i d s . On the other hand, (3-27) suggests that LMA i s not such a good estimate f o r very coarse g r i d s . The examples w i l l illustrate how good LMA i s i n some t y p i c a l situations. 3.3 Examples 3.3.1 Gauss-Seidel R e l a x a t i o n of the U n i t Square with D i r i c h l e t The boundaries be grid i s as Laplace's Equation in Boundary C o n d i t i o n s . follows. We may suppose that c o i n c i d e with g r i d l i n e s . The mesh i n x and y the will h = l / ( n + l ) , and the v a l u e s of x and y which occur at g r i d p o i n t s w i l l be indexed g r i d boundary w i l l The finite where equation from 0 through n+1. P o i n t s not have both difference i n d i c e s i , j between on the 1 and n. system c o n s i s t s of n 2 equations, ( i , j ) has the form (3-28) 4 u ( i , j ) - u ( i - 1 , j ) - u ( i , j - l ) -u(i+1,j) - u ( i , j + l ) = 0, " One unexpected i n f e r e n c e we may draw from (3-24) and (3-27) above i s that LMA works w e l l f o r g i v i n g an asymptotic estimate of the a t t r i t i o n r a t e of the low frequency F o u r i e r components. This_ i s c o n t r a r y to Brandt's expectation (1977) that LMA p r e d i c t i o n s would only be good f o r high frequency components. T h i s f a c t suggests the use of LMA f o r the understanding of the overall performance (as d i s t i n c t from smoothing r a t e ) of r e l a x a t i o n schemes. T h i s w i l l be d i s c u s s e d l a t e r . 36 1 to n. I f any where i , j run from in on values 0 or n+1 then the values w i l l take be taken from the right proceed from the bottom left corner to = x(i+1,j) +x(i,j+1). Then the g e n e r i c equation t r a n s f o r m i n g (3-30) be understood to i t e r a t e x' from the c u r r e n t i t e r a t e x i s 4x'(i,j) -x'(i-1,j) -x'(i,j-l) iteration equation ( l e x i c o g r a p h i c o r d e r i n g ) . The equation to be s o l v e d o b t a i n the next (3-29) this the boundary d a t a . Relaxation w i l l to indices the errors from one to the next i s 4e'(i,j) -e'(i-1,j) -e'(i,j-l) = e(i+1,j) +e(i,j+l). We estimate the d i f f e r e n c e between the r e s u l t of r e l a x i n g a Fourier component f o l l o w i n g through |L1i//(0)| for the LMA estimate the steps of the theorem. F i r s t we by estimate except 2 ( n - l ) e n t r i e s of modulus 1 along the t o p and r i g h t sides, an and | R 1 i / / ( 0 ) | . The nxn g r i d v e c t o r L1i//(0) u(&)\ls(6) i s 0, and entry c o r n e r . R]\p(6) along and and \p(6) 2 the of modulus between 0 and 2 i n the bottom l e f t i s 0 except top and r i g h t f o r 2 ( n - l ) e n t r i e s of M(0)L1^(0) + R1^(0) will | M ( 0 ) | along the bottom and l e f t sides and of modulus 1 along the top and r i g h t moduli one s i d e s , and one of modulus between 0 i n the top r i g h t c o r n e r . Then have e n t r i e s of modulus modulus of the corner e n t r i e s s i d e s . The sum of the i s bounded by 6 ( 1 + | M ( 0 ) | ) . Thus 37 (3-31) The (1+|M(0)I)/(2n-4) < \u(6)L^(8) < ( ]+u(6) R1<//(c9) | + | V(2n + 2) . 2-norm of L i s bounded by 6 and that of L " estimates are obtained i n Appendix B. d i f f e r e n c e between \//' and (3-32) The u(6)\j/(6) become \/\4>(6) | (1 + |M(0) 1 by 1/3. These estimates of the (l/6n)/(2n-4) < - »(d)t(6) I) < (1/2n)/(2n+2). For an n=l00 by 100 g r i d (3-33) .0233 < (10,000 i n t e r i o r -u(e)xP(6) points), .1421. \/\>l>(6) \ < T h i s should be compared with Table I. The values in Table I were obtained as f o l l o w s . For an (n+2)x(n+2) g r i d , values were a s s i g n e d points as per Gauss-Seidel to the d e f i n i t i o n of \p(9) the nxn interior f o r v a r i o u s 0 , 0 . One 1 2 r e l a x a t i o n sweep was performed, t a k i n g the values of boundary p o i n t s to be 0. The norm of the r e l a x e d v e c t o r was computed, as w e l l as i t s p r o j e c t i o n on component Table the . These were both normalized, using I illustrates that as n gets original Fourier | ^ ( 0 ) | = n. l a r g e r , the LMA estimate b e t t e r approximates the r e l a x e d v e c t o r . 38 Table I - Attrition Rates of F o u r i e r G a u s s - S e i d e l R e l a x a t i o n of Laplace's Equation Square (6\d )/ir Grid (nxn) R e l a t i v e norm \r\/\1>ie)\ 2 Components of i n the Unit Projection |M(0)| (<//' ,tf(0) ) / \ $ ( 6 ) | 2 5x5 20x20 100x100 1, 1 1/5,3/5 1/5,1 3/5,3/5 1, 1 .6174 .3829 .1839 .3625 .3028 5859 3453 1 1 22 3303 2731 .7529 .4232 . 1 460 .4004 .3333 4,4 1/5,1/2 1/5,1 1/2,1/2 1. 1 .7201 .4864 . 1 589 .4360 .3256 71 04 4768 1 369 4279 3187 .7529 .5000 . 1 460 .4472 .3333 1/5,1/5 1/5,1/2 1/5,1 1/2,1/2 1, 1 .7464 .4973 .1488 .4449 .3318 7444 4953 1 442 4433 3304 .7529 .5000 . 1 460 .4472 .3333 3.3.2 L i n e R e l a x a t i o n of . 0 1 ( 3 U / 3 X ) + ( 3 u / 3 y ) = 0 2 2 2 2 in the Unit Square. Table II was obtained like the preceeding t a b l e u s i n g y - l i n e G a u s s - S e i d e l r e l a x a t i o n on a 100x100 g r i d . of y-line a sweep r e l a x a t i o n , the values of the new i t e r a t e along an e n t i r e g r i d l i n e of p o i n t s whose second are determined at In index j is the s i m u l t a n e o u s l y , u s i n g values of the new p o i n t s whose second index same iterate i s j-1 and values of the c u r r e n t i t e r a t e f o r p o i n t s whose j index i s j+1. The LMA estimate i s l e s s a c c u r a t e here than f o r the p o i n t G-S r e l a x a t i o n on a comparable g r i d . T h i s i s because the of the matrix L~ 1 independently of h ) . is much larger (although s t i l l norm bounded 39 Table II - L i n e R e l a x a t i o n (0',0 )/7r Norm a f t e r r e l a x a t i o n of F o u r i e r mode u(6) .01,.01 .50,.01 .50,.05 .50,.50 1.0,1.0 .7767 .3907 .2121 .0051 .0025 .9094 .4301 .2187 .0049 .0025 2 3.3.3 Accuracy Distributive of Local Gauss-Seidel Mode (DGS) Analysis Estimates f o r Relaxation of the Cauchy-Riemann Equations We seek s o l u t i o n s (3-34) u , and v to 3u/3x = 3v/3y, -3u/3y = 3v/3x in the u n i t square. The domain on which a s o l u t i o n into square c e l l s of s i d e i s sought h. Values of the d i s c r e t e u are sought at the c e n t e r s of v e r t i c a l c e l l of i s parcelled variable s i d e s and v are sought at the c e n t e r s of h o r i z o n t a l c e l l up values sides. This p l a c i n g of d i s c r e t e v a r i a b l e s ensures second order accuracy of the d i s c r e t e equations (3-35) u(x+h/2,y) - u(x-h/2,y) = v(x,y+h/2) - v(x,y-h/2) u(x,y-h/2) - u(x,y+h/2) = v(x+h/2,y) - v(x-h/2,y). The first equation i s l o c a t e d a t c e l l i s l o c a t e d a t c e l l c o r n e r s . There c e n t e r s , and the i s no n a t u r a l second correspondence between v a r i a b l e s and e q u a t i o n s . To first solve these equations a sweep i s made over a l l the equations and then over the second e q u a t i o n s . To relax 40 each equation each of the four v a r i a b l e s that are that equation relaxation are a d j u s t e d . The may be found details in of the for the amplitudes a , a 1 g r i d . In the LMA 1 relaxation (0) 1 , and a a r i s e s . The r 2 the and residuals were to and 2 the linkage amplitudes A(0) may between a f t e r the sweep. In to . The first are vector below column of each matrix residuals if the the second column i f 2 sweep of DGS reality compared relaxed 1 be the the o r i g i n a l r=\//( 0) , r = 0 , 2 was performed taking boundary values of u and v to be 0. A f t e r both of the r e l a x a t i o n sweep, the r e s i d u a l s were recomputed the norms of the two normalized 1 linear the r e l a t i v e norms of the initial |r |=N, 1 linearly orthogonal they were r =0 , r=\i/( 0) . One phases 2 before and component 1 the = a (0)exp(i(i0 +j0 ) 2 related i d e a l i z a t i o n A(0) original the r norms of the r e l a x e d v e c t o r s the records are and significant with 2 specifying 1 these idealization, 2 amplitudes of r column r a t e of a ' ( 0 ) a f t e r the sweep. Thus a matrix constructed he of F o u r i e r components 2 1 before a 40x40 = a ( 0 ) e x p ( i ( i 0 + j 0 ) , and 1 second III corresponds to the a c t u a l a t t r i t i o n residuals r distributive Brandt(1979). In that paper c a l c u l a t e s a smoothing rate f o r r e s i d u a l s . The in Table involved in by the vector | r j = N-1. 2 transformation residuals The norms third taken. These norms were of the o r i g i n a l r e s i d u a l s , column records of e r r o r s , which i s too d i f f i c u l t the actual to o b t a i n v i a LMA. To o b t a i n these v a l u e s , v e c t o r s u and v were c o n s t r u c t e d , one of set which was to 0. DGS given values as per ^ ( 0 ) , and r e l a x a t i o n was done, taking f1= f2=0, norms of the r e s u l t i n g v e c t o r s were taken and As may be seen from Table the other and was the normalized. I I I the norms of the r e s i d u a l s Table III Cauchy-Riemann e\e LMA 2 Distributive Equations Estimate .333 0 0 .333 Relaxation of 41 the Actual Residual Actual Error .110 0 0 . 1 07 .548 .038 .038 .548 7T/2 , TT/2 .440 0 0 .440 721 0 0 .720 .444 .022 .022 .444 .020 0 0 .020 709 0 0 .716 .119 .025 .035 .079 460 0 0 .460 727 0 0 .722 .465 .020 .021 .454 TT/40 , TT/40 .997 0 0 .997 .997 0 0 .968 .973 .012 .012 .973 7T/40 , 7T 7T/40 , 7T/2 after r e l a x a t i o n are c o n s i d e r a b l y d i f f e r e n t p r e d i c t e d on the b a s i s of LMA. the after from what would be In f a c t the LMA prediction for r e s i d u a l s seems to conform more to the norm of the e r r o r s r e l a x a t i o n though the d i f f e r e n c e this situation estimates then, LMA does is not still yield of the e f f e c t of r e l a x a t i o n on the components. For t h i s weakly e l l i p t i c complicated and cannot be d e s c r i b e d marked. very In accurate various Fourier system, the r e l a x a t i o n i s i n the terms used to prove Theorem I. 3.4 Comments and M i s c e l l a n e o u s 3.4.1 In Consequences A l t e r n a t e F o u r i e r Components the above discussion only exponential components have been used. There i s no reason could not have embodied the various F o u r i e r components. With the usual component corresponding to a priori frequencies indexing convention frequencies 0 1 and 6 2 Fourier why as we sine such a would be 42 defined as (3-36) u ( 0 ) ( i , j ) = In fact, sine eigenvectors sin(0 i)sin(0 j). 1 Fourier of the whereas exponential scheme. LMA However, components Jacobi CJ above are one s e r i o u s d i f f i c u l t y a r i s e s ; we for the rates of at adjacent attrition same p o i n t s . As any cannot f o r sine F o u r i e r components f o r lack of a p p r o p r i a t e values the r e l a x a t i o n of the model problem, do simple it turns of sine F o u r i e r components are q u i t e d i f f e r e n t from e x p o n e n t i a l the like components \jj are not e i g e n v e c t o r s i d e n t i t i e s connecting out 2 Fourier components for the r e l a x a t i o n on the same problem. T h i s c a s t s some doubt on i n t e r c h a n g e a b i l i t y of the two 3.4.2 We from an components. A l t e r n a t e Norms have used the 2-norm throughout, because this arises inner product which must be used to compute orthogonal p r o j e c t i o n s . The then types of F o u r i e r the reader may result of check that i f the 1-norm were used the main theorem (whose developement i s independent of the norm used) i s (3-37) I * ' - u(e)4>{6)\ = 0 ( I V (67) | /i/N). If the supremum norm i s used no since the boundary values g r i d vector 3.4.3 contribution i s not of in p r o p o r t i o n useful to t h e i r bound is obtained, to the norm of a numbers. Residuals Using LMA it is possible to justify the intuitively 43 reasonable c l a i m that (3-38) r reflect = b - AX = local errors errors, so that the r e s i d u a l s -Ae, i n the current the estimate of the e r r o r size i t e r a t e X more than of r e s i d u a l s alone i s not a good i n an approximate s o l u t i o n . We must a l s o know whether the e r r o r i s smooth or o s c i l l a t o r y . F i r s t we prove a lemma which j u s t i f i e s type computations of residuals, computations to show that the low contribute frequencies much less to and of LMA use these frequencies in the error residuals than the then use the high do. (3-39) a U ( i , j ) equation + a U(i-1,j) + a U(i,j~1) 1 2 3 + a"U(i-1,j) + a U ( i , j + l ) 5 = h b(i,j) 2 i s t o be s a t i s f i e d at every i n t e r i o r point to the we LEMMA 3: Suppose the d i f f e r e n c e covering global ( i , j ) of a grid G a convex domain fl, and that zero order approximations the boundary of fl and t o D i r i c h l e t boundary c o n d i t i o n s a r e made. Let \p(6) be as b e f o r e . associated with Then as h -> 0 the an e r r o r of the form \p{6) are a s y m p t o t i c a l l y given by (3-40) r ( i , j ) residuals = ( a 1 + a exp(-i0 ) + a exp(-i0 ) 2 1 + a*exp(i0 ) + a e x p ( i 0 ) ) 1 5 2 3 2 \p(d){i,j). PROOF: Consider the system of equations 44 (3-41 ) ( I + A )X = 0. Suppose we r e l a x t h i s system by s o l v i n g an equation of the form (3-42) IX' This = -AX. iteration the e r r o r s a t i s f i e s a l l the c o n d i t i o n s of Theorem i t e r a t e X i s \p(6) then the i n the c u r r e n t I. error If in the next i t e r a t e X' i s r , the r e s i d u a l f o r AU = h b . Theorem I 2 then says that (3-43) | r - M(0)*(0) I = 0(1/VN). But y(0) f o r t h i s r e l a x a t i o n the i s given by (3-40). QED In our d e s c r i p t i o n of r e s i d u a l computation we assume that sum i s r e s t o r e d by a +a +a +a +a 1 2 3 4 i s 0. F u l l 5 generality adding a term a°U(i,j) i f the f u n c t i o n u itself appears differential (3-39). Consider happens equation which u n d e r l i e s to .a exp(i(0 i + 0 j) 1 2 (3-44) r ( i , j ) single Fourier component of i n the the what form under the mapping e->r: = e x p ( i ( 0 i + 0 j))(a° + a 1 2 + a exp(- 0 ) 1 2 1 t + a exp(-i0 ) + a"exp(i0 ) + a e x p ( i 0 ) ) 3 If 0 1 close 2 and 0 to of r ( i , j ) 2 1 a r e both s m a l l , 5 2 then exp(±i0 )and 1. At an i n t e r i o r p o i n t 1 exp(±i0 ) 2 ( i , j ) therefore i s c l o s e to |a°|. However i f 0 1 and 0 the modulus 2 are c l o s e to TT then exp(± i 0 ) and exp(±i0 ) are a l l n e a r l y 1 2 are both -1, and 45 | D ( e ) ( i , j ) | i s c l o s e t o |a° + 2 a | . Lemma 3 i n s u r e s that 1 indicates of the the s i z e of the g l o b a l same practice), order the high over-represented in as 2a amplitude 0(h) of the error. On frequency relaxation amplitudes other very f r e q u e n c i e s are on the of much hand corresponding the amplitude times of the model problem a error. order of i s 1 + 0 ( h ) , we see from (3-44) the r e s i d u a l can be many the high frequency are the r e s i d u a l v i s - a - v i s the low frequency of a very low frequency amplitude the 0 ( i t i s t y p i c a l l y much smaller i n 1 the mesh s i z e h, and exp(0(h)) the r e s i d u a l s . Thus u n l e s s - a i s frequency components. Since the lowest that this 0 residual low of ( f o r the + 2a 1 a i s only frequency very high Gauss-Seidel = 8) the amplitude of 46 CHAPTER FOUR A HEURISTIC LOOK AT RELAXATIONS USING LOCAL MODE ANALYSIS In Chapter in the Three i t was limit of very shown that LMA f i n e g r i d s . I t was v a l i d i t y of the f u n c t i o n M(0) was 6\8 2 two could be justified a l s o shown that the values of bounded away from 0. Here we compute the f u n c t i o n u(6) for common schemes: the SOR of uid) not r e s t r i c t e d method and the ADI w i l l h e l p us understand to method. The form h e u r i s t i c a l l y what these schemes do. The a l g e b r a i c theory of these schemes f o r s p e c i a l cases will be adduced to c o n f i r m the understanding using LMA. also be used to d i s c u s s why relaxations which Three do not efficiently the conditions reduce i t i s that most outlined in Chapter LMA will meet smooth e r r o r s . 4 . 1 The The commonly SOR used SOR Method (Successive for constant success. A parameter performed. A co s e t s the i.e., the performance minimize a m p l i f i c a t i o n matrix. Here we method coefficient problems amount c r u c i a l q u e s t i o n i s how order to o p t i m i z e iterations, Over-Relaxation) the some over radius behaviour of method a p p l i e d t o the Laplace equation i n a domain J2 been d i s c r e t i z e d on a g r i d with mesh s i z e with over-relaxation scheme spectral study the been co (0 < co < 2) i n to p i c k of the of has many of the the SOR which has h. The d i s c r e t e Laplace equation at a p o i n t ( i , j ) i s given by (4-1) U ( i , j ) = ( 1 / 4 ) ( U ( i - 1 , j ) +U(i,j-1) +U(i+1,j) +U(i,j+1)). 47 Applying (i, j) (4-2) the SOR method, the equation of r e l a x a t i o n at a p o i n t is X'(i,j) = (l-co)X(i,j) + (cu/4)( + X(i+1,j) + X(i,j+1) ( 4 - 1 ) from Subtracting X'(i-1,j) + X'(i,j-1) ). ( 4 - 2 ) we obtain the transformation of errors (4-3) e'(i,j) = ( l - u ) e ( i , j ) + (co/4)( e ' ( i - 1 , j ) + e ' ( i , j - 1 ) + e(i+1,j) + e(i,j+1) e(i,j) = A(0)exp(i(0 i+0 j)) Let e'(i,j) 1 (4-4) let = A' ( 0 ) e x p ( i ( 0 i + 0 j ) ) . 1 2 A * ( 0 ) = (i-u)A(0) + ( w / 4 ) ( A' ( 0 ) e x p ( - i 0 ) 1 + A'(0)exp(-id ) + A(0)exp(i0 ) 2 LMA estimate (4-5) + A(0)exp(i0 ) 1 2 ). i s then given by M ( 0 ) = (exp( it? )+exp( i 0 ) - 4 ( w - 1 )/w) 1 / 2 (4/co - e x p ( - t 0 ) 1 - exp(-c0 )). 2 The maximum of | M ( 0 ) | occurs at 0 The minimum occurs a t The and 2 ( 4 - 3 ) becomes Then The ). profile (0,TT) or of the s u r f a c e uid) 1 = 8 (TT,0) 2 = 0 and has the value 1 . and has the value CJ - 1. resembles t h a t f o r Gauss-Seidel r e l a x a t i o n , with a l o c a l maximum a t (7r,7r), see F i g u r e 5 . 48 If CJ i s near 1 then the values ca = 1 (Gauss-Seidel) F i g u r e 5 - Contour distributed, cluster that but u(6) of well CJ = 1.5 L i n e s of | M ( 0 ) | f o r SOR R e l a x a t i o n as CJ i n c r e a s e s to 2 eigenvalues of the u(6) the values of t o g e t h e r . For model problems i t i s known the are (Young, amplification 1971) matrix behave similar i l y . the Can u(6) SOR relaxation? M(7rh,7rh) the be used t o estimate the r a t e of Unlike de Vries was not a good estimate of amplification matrix. There of the relaxation (since be used). T y p i c a l l y 1 - K than 1 - Ai(7rh,7rh) problems Nevertheless, in certain cases 2 in to M(0,0) = 1 the s e v e r a l v a l u e s of systems were h, as o u t l i n e d always, i t square. CJ as d e s c r i b e d below. (x-1/2) 1 + ( y - 1 / 2 ) <l/4; 2 2 2 difference smoothest unit fi i s the curved region { ( x , y ) | y>0, 4 / 5 - ( 2 x - l ) < Finite the i n the f o l l o w i n g domains: 0 circle 3 of i t i s p o s s i b l e to use LMA t o We c o n s i d e r L a p l a c e ' s equation square; fi i s the found that was two or t h r e e times b i g g e r estimate the best value of the parameter i s the u n i t I the s p e c t r a l r a d i u s K analogous cannot for (1982), of i s no other c l e a r c h o i c e of F o u r i e r component which should be eigenvector convergence y <1-(2x-1) }. 2 set up f o r these problems f o r i n Table IV. I t was not clear 49 which frequency ought to be analogous to the smoothest eigenvector f o r the domains fi and fl . Since the number of points readily 2 is computed, the value of M was computed f o r ft = (7r//N, 7r/v/N) . T h i s i s n e a r l y e q u i v a l e n t unit square. The value to the actual of the cj of CJ which minimizes \u\ i s recorded r radius of the CJ* is i n c l u d e d . For the u n i t square analysis in (7rh,7rh) C J * of in Table IV. For comparison purposes the value minimizes grid 3 CJ which a m p l i f i c a t i o n matrix i s known from a complete t h i s model problem (see Young(1971)), and i s given by (4-6) CJ* = 2/( l + s i n U h ) ) . For the other determined K' problems, using the an approximation successive iterates K, to X ,X ,... 1 2 was . It i s d e f i n e d by (4-7) K ' = |X An - X | 52 / |X 5 1 51 - X |. 50 C J * was computed which minimized We large CJ' see that N. For the boundary fl 3 K ' as a f u n c t i o n of C J . i s a good e s t i m a t e of C J * particularly, for the r a t i o of the number of g r i d p o i n t s next to to the t o t a l number of interior grid points r a t h e r l a r g e r than f o r the other domains with comparable validity of LMA fl 1 Can and LMA N. The depends on t h i s r a t i o being s m a l l , and so i t i s not s u r p r i s i n g that estimates f o r for is J2 were not as good as 3 those fl . 2 be used to choose parameters f o r more general 50 Table IV - Best values of C J h n 2 N CJ* .2 16 1 .236 1 .260 .1 .05 81 1 .524 1 .528 361 .02 2401 1 .729 1 .882 1 .730 1 .882 .02 1847 1 .86 1.89 .02 .01 419 1787 1 .73 1 .58 1.79 .005 7371 1 .86 1 .93 ' problems? I c o n s i d e r e d extension to two asymmetric constant c o e f f i c i e n t variable classes of equations and equations with one form -9 /9x 2 two d i f f i c u l t i e s with asymmetric o p e r a t o r s of 2 -9 /9y 2 . The f u n c t i o n UL{6) computed - a9/9x 2 by LMA f o r such an operator may have a maximum at than at ( 0 , 0 ) f o r some plausiblity Secondly, in the computing values of C J . There u, f o r a eigenfunctions of low the frequency the For for integers m,n). all. less component. sin(m7ry) h of order l a r g e r , the d i s c r e t e e i g e n f u n c t i o n s are q u i t e d i s s i m i l a r F o u r i e r components, and thus we cannot rather asymmetric operator a r e sin(n7rx) square, (TT,0) i s then even skewed e x p o n e n t i a l l y (u(x,y) = exp(-ax/2) unit equations; coefficient. I encountered the 1 .93 in 1/a or to the expect LMA to be v a l i d at S e v e r a l methods to estimate optimal parameter values u s i n g LMA were t r i e d and none approached the e x p e r i m e n t a l l y determined optimal v a l u e . Extension to problems with more successful. f o l l o w i n g problem: Numerical one variable experiments were coefficient done was using the 51 (4-8) (1/10 +x + y )9 u/9x 2 2 2 + 9 u/9y 2 2 with D i r i c h l e t boundary c o n d i t i o n s on of .02 was of co* of 9 u/9x 2 was 2 was determined included Apparently SOR of LMA of as The outlined in the d i s c r e t e operator f a i l e d to give 1 and c o * . has some l i m i t e d u s e f u l n e s s with respect parameter ratio i n s e n s i t i v e to the experimentally, to f o r s o l u t i o n of Laplace type equations r e l a x a t i o n . T h i s i s p a r t i c u l a r l y true small size S i m i l a r comparisons where a term such good agreement between co choice co' was mesh over the range .1 to 2.; CJ' = 1.88. 2 p r e v i o u s l y , to be co* = 1.90. -u/lOh [0,1]X[0,1]. A used. The determination of coefficient value = 0, 2 boundary for domains with l e n g t h to area, approximated by g r i d s . However, there i s no obvious way to use LMA to by a fine estimate the a c t u a l r a t e of convergence with the optimal parameter. 4.2 Alternating-Direction Implicit For ADI is s p l i t any at i t e r a t i o n on an equation i n t o matrices H and 1 V . Iteration AX = h b, 2 the matrix For any g r i d v e c t o r 1 X, and g r i d p o i n t ( i , j) , ( H X ) ( i , j ) i n v o l v e s only the values of 1 the p o i n t s ( i , j ) , (i-1,j), and Thus H 1 and V 1 represent the f i n i t e - d i f f e r e n c i n g y d i r e c t i o n s , r e s p e c t i v e l y . A parameter iteration to r e l a x a t i o n . Optimal the next, determines c h o i c e of values for i t e r a t i o n of ADI c o n s i s t s of and the p solving is (i,j + l). i n the x and p, which may e n s u r i n g quick convergence. One X ( i + 1 , j ) , and V ' X ( i , j ) i n v o l v e s only the values of X at the p o i n t s ( i , j) , ( i , j - 1 ) , one A vary form of important" from the in 52 (4-9) (H + ph I 1 and t h e n (4-10) to (V + ph I 1 an If 1971). on the next = (ph I - H 2 iterate two m a t r i c e s guarantees The H 2 X' 1 )X* + h b , 2 from t h e matrices equations rectangular if commute of the current convergence explanation Chapter is of 1 commute, iterate only X, using If frequently for algebraic iterates the (see special Young case of + cu = h ( x , y ) these c o n d i t i o n s are not met, (for derivative are present) then slow phenomenon 2 provides the the form first-order this of then + (3/3y)(g(y)3u/3y) domains. large and V 1 r a p i d convergence (3/3x)(f(x)3u/3x) example, of )X + h b , 1 intermediate. separable (4-11) )X' 2 the theory - V 2 solving obtain X* as )X* = ( p h I 2 in using terms practice. LMA i s some j u s t i f i c a t i o n An heuristic now g i v e n . for this Theorem 1 non-rigorous approach. The l o c a l Laplace's (i,j) (4-12) mode a n a l y s i s equation consists ((2 of r u n s as 2 2 the follows. ADI method The r e l a x a t i o n applied at grid to point solving + ph )X*(i,j) = ((ph for - - X*(i-1,j) 2)X(i,j) - + X(i,j-1) X*(i+1,j)) + X(i,j+1) ) + h b(i,j), 2 53 and t h e n s o l v i n g (4-13) ( ( 2 + p h ) X ' ( i , j ) - X ' ( i , j - 1 ) - X ' ( i , j + D) 2 = ( ( p h - 2 ) X * ( i j ) +x*(i-l,j) +X*(i+l,j)) be the e r r o r 2 r Let e* the corresponding (4-14) error in X*. Then + h b(i,j). 2 transformation t o (4-13) i s ((2 + p h ) e * ( i , j ) - e * ( i , j - 1 ) - e * ( i , j + l ) ) 2 = ((ph - 2)e(i,j) + e(i,j-D 2 Let e(i,j) = A(0)exp(i(0 i let e*(i,j) + 0 j)), 1 2 + e(i,j+l) and = A * ( 0 ) e x p ( i ( e ' i + & j ) ) . Then 2 (4-15) A * ( 0 ) ( 2 + p h Using similar 2 ). (4-15) - 2 c o s 0 ) = A ( 0 ) ( -2 + p h 1 becomes + 2 2cos0 ). 2 a r g u m e n t s we a r r i v e a l s o a t (4-16) A ' ( 0 ) ( 2 + p h 2 - 2cos0 ) = A*(0)( 2 -2 + p h 2 + 2COS0 ). 1 Thus n(6) (4-17) where v with 8 - S V - t V + S V + t = ph , 2 Notice any V that by small s = 2 - 2 c o s 0 , and 1 i t i s p o s s i b l e t o make s u i t a b l e choice 0 1 - cos(0(h)) 1 or small = 0(h ). 2 6 2 As of 2 ix(d) e q u a l p. F o r s m a l l values p t = 2 - 2cos0 . are being i s increased, p to 0 for t h e components eliminated the zero of since M(0) 54 occurs the at zero arrived is at Now by 9 u/9x 2 mentioned for this grid (4-19) (4-20) Notice u(6) the same calculations range for p = for i n Young p 4/h that 2 was (1971). PDE 9 u/9y 2 above, the until ADI The - 2 a9u/9x methods finite - b9u/9y are not = 0. expected d i f f e r e n c e form to of work the as well equation at + d+ah)U(i-1,j) + + (1+bh)U(i,j-1) (1-bh)U(i,j+1) = 0. yields M(0) = that of |ju(0)| 1 i f ah v ~ s - 2iahsin6 v + s + 2iahsin6 i t i s not to 0 still and frequencies bh and provided Thus reduce although its is difficult (l-ah)U(i+1,j) equal (0(h)) + 2 -4U(i,j) LMA This frequencies, (i, j) is + Then more equation. point higher (7T,7r). consider (4-18) As increasingly occurs that overall is smooth the given ADI 0. 2ibhsin6 2ibhsin6 choose While for at low frequencies, 0(1). For large although there the p o s s i b l e to f o r any are ^ v - t v + t + 1 1 ah i t s value and no performance not of will still be 2 in order small p i t s value minimum 0, i t to the minimum i s of occurs will make at be order high small 0(h" ). 1 p frequency) will the i s not are choice (low method bh r p 2 which error effectively components. effectively significantly will smooth degraded the from Thus error its 55 performance on the model problem. 4.3 Smoothing Property of R e l a x a t i o n s Suppose we have a d i f f e r e n c e equation at every interior point of a g r i d c o v e r i n g a domain fl; (4-21) a U ( i , j ) + a U ( i - 1 , j ) + a U ( i , j - l ) 1 2 3 + a U ( i + 1 , j ) + a U(i,j+D = h b ( i , j ) , 4 5 where a , . . . , a 1 are constant c o e f f i c i e n t s . Suppose there 5 term of degree zero i n (4-22) a +a 1 + a 2 2 3 u is no i n the u n d e r l y i n g PDE so that + a• + a = 0. 5 Suppose we have a s p l i t t i n g A = L - R which y i e l d s a r e l a x a t i o n scheme f o r (4-22), (4-23) l X ' ( i , j ) + l X ' ( i - 1 , j ) + l X ' ( i , j - D 1 2 3 + 1 " X ' ( i + 1 , j) + 1 X ' ( i , j + i ) 5 = r X(i,j) + r X(i-1,j) + r X(i,j-l) 1 2 3 + r X(i+1,j) + r X(i,j+l) +h b(i,j). 4 5 2 Then s i n c e 1 * - r * = a*, * = 1,...,5, we must (4-24) 1 1 + 1 2 + 1 The reader w i l l 3 + 1 « + 1 recall 5 = have r +r +r +r"+r . 1 2 3 5 from Chapter Three that the f u n c t i o n u(6) for the r e l a x a t i o n scheme (4-24) i s 56 (4-25) u(6) = ( r + r exp(-i0 ) + r exp(-i0 ) 1 2 + r'expdfl ) ( l 5 2 1 Thus, 3 2 + l exp(id ) 5 ). 2 if (4-26) l 1 + l 2 + l 3 + 1" + l (4-27) u(e) - > 1 , as e\e 2 The limit 6 are 2 + l exp(-t0 ) 1 + l"exp(i0 ) 2 ) 2 + l exp(-i0 ) 1 3 + r exp(t0 ) 1 / 1 (4-28) says that near coefficient 0. Thus 5 * 0, 1 -> o. the for function the u(6) is relaxation near of equation by a scheme which s a t i s f i e s the conditions in e r r o r are not e f f i c i e n t l y in general is r e d u c t i o n of the e r r o r efficient 6 \ constant i n Chapter Three, the low frequency F o u r i e r for if a outlined the 1 components reduced. The best we can hope smoothing, but not efficient overall. Only the ADI r e l a x a t i o n of the model problem with a small value of the parameter r comes c l o s e to violating this condition. The reader will recall that low frequency e r r o r s c o u l d be e f f i c i e n t l y r e l a x e d i n t h i s c a s e . 1 57 CHAPTER FIVE DISCUSSION 5.1 A p p l i c a b i l i t y of L o c a l Mode A n a l y s i s The results following in conditions. section The of the A l t e r n a t e l y each several coefficient each LMA of i s based matrix step could which into be the on a A = L - R. a sequence of i s based on a s p l i t t i n g . Although the s p l i t t i n g ( s ) need not be r e g u l a r , they must be homogeneous, i . e . , each relative to the red-black r e l a x a t i o n s are excluded. A f u r t h e r c o n d i t i o n , which we imposed set point under f o r which LMA i s done i t e r a t i o n which relaxation procedures justify relaxation must be a l i n e a r f i x e d p o i n t splitting 2.2 of must be i n the same p o s i t i o n previously f o r ease of a n a l y s i s , was that equation at a given p o i n t , no other than relaxed points. Thus i n the r e l a x a t i o n neighbouring p o i n t s occur those which occur i n the o r i g i n a l d i f f e r e n c e equation at that p o i n t . We considered based on f i v e point which only second order d i f f e r e n c e schemes formulae, which excluded the r e l a x a t i o n at point relaxations ( i , j ) r e q u i r e d values neighbours. This for at p o i n t s other than the four nearest condition may not, i n f a c t , be neccessary f o r the r e s u l t s but no r e l a x a t i o n scheme i s known to us f o r which i t i s not t r u e . The c o n d i t i o n that the r e l a x a t i o n be based on a s p l i t t i n g of the c o e f f i c i e n t matrix was intended to ensure that Fourier a n a l y s i s was meaningful. We do not think that F o u r i e r analysis is relevant to Conjugate Gradient The the understanding of schemes such as the Method. homogeneity condition guaranteed that the 58 cross-coupling of F o u r i e r components undergoing r e l a x a t i o n asymptotically negligible; i . e . , that each F o u r i e r component behaved l i k e an eigenvector. possible inhomogeneous for relaxation; Fourier then invariant Fourier components (7T-0 significantly for four-dimensional all such they than as is red-black of the must be arranged in being treated like ordered Gauss-Seidel r e l a x a t i o n , frequencies 1 analysis cross-coupling important, with Fourier schemes rather and 2 local because Under red-black (TT-0 ,0 ), 1 is subspaces eigenvectors. A however, components was ,7T-0 ) mesh invariant 1 feed 2 2 into s i z e s , and subspace (0 ,7T-0 ), (0 ,0 ), 2 1 each other must be t r e a t e d as a (see Trottenberg and Stueben, 1982). Although only for the theorem relaxations boundary conditions above, some extensions with Neumann or of which may mixed Chapter Three was scalar met be problems the established with Dirichlet requirements contemplated. boundary approximations to D i r i c h l e t using in described Scalar problems c o n d i t i o n s or higher boundary c o n d i t i o n s may be the method i n Chapter Three, however, the a l g e b r a much more complicated. result elliptic to strongly systems I t seems s t r a i g h t f o r w a r d elliptic may be systems. possible, treated may to extend Extension but order be the to weakly we anticipate d i f f i c u l t i e s with d i s t r i b u t i v e r e l a x a t i o n schemes. Of course a theorem like that in Chapter Three i s not meaningful f o r v a r i a b l e c o e f f i c i e n t problems or n o n - l i n e a r 5.2 The P o t e n t i a l of LMA Outside the MG The theorem i n Chapter Three j u s t i f i e d problems. Context our explorations 59 in Chapter Four. We saw that LMA behaviour of the SOR and ADI linear reduce smooth e r r o r s p o o r l y . We relaxations that LMA gain has could be used to e x p l a i n methods, and a l s o the fact speculate to f u r t h e r h e u r i s t i c understanding of such methods. a s c i e n t i s t considers the use of a p a r t i c u l a r r e l a x a t i o n scheme for the s o l u t i o n of a problem and know how appropriate rigourous eigenvalues the that the p o t e n t i a l to be more widely used in order Frequently A the i t i s or how algebraic and would involve a n a l y t i c a l t o o l and finding i t can the other hand LMA y i e l d much the same The theorem in Chapter Three e s t a b l i s h e s that tool and problems. However, we imparts confidence have only in shown that LMA its LMA is a use for MG yields a good of the r e s u l t s of r e l a x a t i o n where the e r r o r i s a m u l t i p l e of a s i n g l e F o u r i e r component. As w i l l it to information. P o s s i b l e Developments of the Theorem below the i s a simple 5.3 description parameters. i s g e n e r a l l y more work than i s needed the o r i g i n a l problem. On meanigful to the eigenspaces of the a m p l i f i c a t i o n matrix of r e l a x a t i o n . This solve answer to pick r e l e v a n t wishes would be very be discussed u s e f u l to have a r e s u l t of t h i s kind for a r b i t r a r y e r r o r s . Suppose we it have an a r b i t r a r y e r r o r e= x-u and decompose as (5-1 ) e = £ A ( 0 H ( 0 ) . e We ask resemble how e* closely the LMA does e' estimate, the error after relaxation which i s given by 60 (5-2) e* = £^(0^(0)^(0). e An argument along the l i n e s of made that the d i f f e r e n c e will sketch regions so defined that and transform = this of e' - e* here. the the We 0,1,...,N-1. Note that (0 , 7r/2 )U( Ztt/2 , 27r) with these equivalent t o the p r e v i o u s may be c o n s t r u c t e d (5-3) for 1 might the s m a l l . We rectangular (5-1) i s form the f i n i t e low be well Fourier 27rk/N f o r k frequency region c h o i c e s of f r e q u e n c i e s . d e f i n i t i o n . Now m a t r i c e s is This is R* and L* of 0 so that and L*«//(0) = l*(0)<//(0), = 27rk/N. Then 2 only take the values 2 independently R*t//(0) = r*(e)xp(6) 0 ,0 A(0) and 0 1 consider decomposition coefficients e , where 0 Three i s asymptotically will Fourier Chapter since r*(0)/l*(0) = M(0), formula (3-17) reads (5-4) e' = where e* + L " ( L 1 e * - R1e), 1 again no guarantee relative the that Lie* obtained arbitrary Does algorithms factor Lie* to e. In f a c t boundary, estimate as L1 = L* - L and R1 = R* - R. Now however we have and R1e are asymptotically small i f e i s non zero only on p o i n t s next t o and R1e can be l a r g e f o r any N. Thus the i n Chapter Three may not be extended to errors. this fact that r e l y invalidate on LMA convergence estimates for estimating the of MG smoothing f o r any e r r o r ? I think not. The smoothing rate i s taken M=max|M(0)| f o r one of 0 \ 0 2 i n the range (7r/2,3 7r/2) . Most 61 of the high factor frequency F o u r i e r much smaller reduced er of an a r b i t r a r y e r r o r 1 previously. method Trottenberg in and much correction component step) es = e - er Thus a proof that the reduces e So the Stuben of the e r r o r r e l a x a t i o n on non-smooth f a r proofs been d i f f i c u l t than (1981) have already shown been problem error reduces in the the fine i s a l l that restricted highly to a few artificial. together theory and 5.4 coarse the smooth grid strongly i s needed to show algorithm. of the convergence of the and (the f i n e g r i d problem. A MG method proof v i a the LMA estimate toward bringing application. Summary There seems to be p o t e n t i a l f o r development of LMA practical r o l e f o r which Brandt seen, under s u i t a b l e c o n d i t i o n s , quantitative relaxation rigourous and is spanned by introduced local information schemes o u t s i d e the MG the {4/(6) beyond i t . As we have Fourier analysis can the performance of about context. thus i s p o s s i b l y part of b a s i s f o r the MG This have simple c a s e s . They have which i s used i n p r a c t i c e would go a long way 1 is generality significantly convergence of the whole MG yield true u + 0(i/VN). more that d i r e c t s o l u t i o n of the coarse g r i d the a t h i s c o u l d be proved then i t would be p o s s i b l e to show convergence of the MG grid by i t is s t i l l reduced by a f a c t o r which i s bounded by If also are than u. I b e l i e v e that that the non-smooth part still components a I t can a l s o be made sound mathematical method. orthogonal \ 6 1 , d = 2nk/U, 2 p r o j e c t i o n of one of 6\6 2 e e onto the (n/2 space , 3TT/2 )} 62 BIBLIOGRAPHY Brandt, A. " M u l t i - L e v e l Adaptive S o l u t i o n s to Problems", Brandt, A. Methods, Nov Mathematics "Guide of Computation, in Multiqrid, W. Hackbusch and U. T r o t t e n b e r g , eds., Koeln-Porz, Hackbusch, Multgrid (1977), 333-390. Development", 1981. LNM to 31 Boundary-Value 960, S p r i n g e r V e r l a g W. "Multi-Grid P e r t u r b a t i o n Problem", (1982), 220-312. Convergence for a Singular Proceedings, Copper Mountain, Apr 1983. (to appear). Hemker, P. W. "Fourier Analysis of P r o l o n g a t i o n s and R e s t r i c t i o n s " , Report NW Centre, Amsterdam Hemker, P. W. "The - Functions, 93/80, Mathematical (1980). Incomplete LU-Decomposition Method i n M u l t i - G r i d A l g o r i t h m s " , Layers Grid Computational and in as a R e l a x a t i o n Boundary Asymptotic and Methods,. Interior J. Miller ed., Boole P r e s s , D u b l i n (1980), 402-406. Hemker, P. W. "An Accurate Method Without D i r e c t i o n a l Bias f o r the Numerical S o l u t i o n of a 2-D Problem", (1981). Report NW 117/81, Elliptic Singular Mathematical Perturbation Centre, Amsterdam 63 Mol, W.J. "Smoothing and Coarse G r i d Approximation of M u l t i g r i d Methods", Report NW Amsterdam Mathematical Centre, (1981). Trottenberg, U. and M u l t i g r i d Methods, K. Stueben. W. Hackbusch Koeln-Porz, Nov 1981. LNM de V r i e s , H.B. "The Two-Level the I n i t i a l Value Problem Report NW 110/81, Properties 960, "Multigrid and U. Methods", Trottenberg, Springer Verlag Iterative Academic P r e s s , New York, eds., (1982), 1-176. A l g o r i t h m f o r the Solution of f o r P a r t i a l D i f f e r e n t i a l Equations", 136/82, Mathematical Centre, Amsterdam Young, D. M. in Solution 1971. of Large (1982). Linear Systems. 64 Appendix A Proof of Lemma 3 , C h a p t e r 7(t) Let differentiable [0,L]-> : curve R 2 Two be a parametrized continuous by arc piecewise length. For 0 < t < L, l e t H ( t , h ) be (A-1) H(t,h) = {(x,y)| Define A l ( t ) as the a r e a 7 of H ( t , h ) , Lemma A 1 : I f 7 i s a s t r a i g h t Proof: denoted line in is Trh 2 + 2hL continuous non-decreasing differentiable whenever differentiable; then Al'(t) at a point a A l ( t ) = A( U D ( 7 ( r ) , h ) , 0 < r < t ) . 2 (A-2) = 2h. dAl/dt L e t D(a,h) be a d i s k o f r a d i u s h c e n t e r e d R . Then c l e a r l y a then A(H(t,h)) Obvious. Lemma A 2 : A ( H ( L , h ) ) = A l ( L ) < Proof: < h f o r some t * , 0<t*<t} |(x,y)- (t*)| (see 7 function of Clearly A l t. i s . L e t t be a p o i n t where Figure) = l i m (A(H(t+e,h)) - A ( H ( t , h ) ) ) / e 6-HD = l i m (A(D(7(t)+e '(t),h)UH(t,h)) 7 - A(H(t,h))) e-K) = lim A(D( (t)+e7'(t),h)\H(t,h)) 7 / e < lim A(D( (t)+e7'(t),h)\D( (t),h)) e-H3 7 It 7 / e. / e is 7is 65 7(t) HI D( (t),h) 7 ///D( (t)+e '(t),h) T 7(t)+e7'(t) T Figure T h i s l a s t l i m i t i s independent specified). If of 7 ( s i n c e 7 i s a straight l i n e the l a s t e q u a l i t y and f o r t h i s case A l ' ( t ) (A-3) Al'(t) |7'(t)| = 1 was inequality i s an = 2h by Lemma A1. Thus < 2h. Since A l ( 0 ) = 7rh , the lemma f o l l o w s . QED 2 Lemma A3: (Lemma 2 of Chapter with piecewise smooth Three): L e t fl be a boundary 7. There i s a constant such that i f G i s any s u f f i c i e n t l y f i n e g r i d , number of g r i d p o i n t s of number domain and N a i s the G which l i e i n t e r i o r t o fl, then the M of such p o i n t s which l i e next t o a g r i d p o i n t which i s o u t s i d e fl i s bounded by (A-4) M < a/N Proof: his (I am indebted t o David Boyd of the Math. Dept. UBC f o r suggestion of u s i n g g r i d Let be squares) A be the area of fl, L be the l e n g t h of the mesh s i z e of G. L e t n be the number of g r i d (each of area h ) which l i e e n t i r e l y 2 number of squares 7, and which intersect in fl. 7, Let or m h squares be the whose boundary 66 intersects 7. Then (A-5) A < (m+n)h . 2 The squares which meet 7 must m H(L,2/2"h) d e f i n e d above. Now (A-6) mh Choose h < /2L/(8TT). Then (A-7) mh Now 2 2 l i e entirely in the by Lemma A2 < A(H(L,2v/2h) < 8vrh 2 + 4v/"2hL. < 5/2Lh. choose h < A/(10/2L). C l e a r l y A(H(L,2/2h)) < A/2, (A-7) and (A-9) (A-8) n h > A/2. 2 set thus by Thus (A-9) m < 5/2L/h = (10L//A)(/A//2h) < (10L//A)/n. Now one each of the n squares i n t e r i o r to S2 c o n t r i b u t e at p o i n t to N (say the lower l e f t c o r n e r ) . Thus of the M p o i n t s which l i e next to a corner of one of the m point outside least n < N. Each £2, is a squares. C l e a r l y there are at most three c o n t r i b u t i o n s to M from each of the m squares. Thus (A-10) M < 3m < (30L/v/A)/n < (30L//A)/N. QED 67 Appendix B Norms of the matrices L and L " f o r the G a u s s - S e i d e l 1 R e l a x a t i o n of the Model Problem 1. Form of L and For an nxn L" 1 g r i d with n matrix L has the 2 interior points form A 0 0 -I A 0 0 -I A (B-1 ) 0 0 ... 0 where A i s an nxn matrix of the (B-2) form 4 0 0 0 -1 4 0 0 • 4' 0 • 0 -1 . . . 0 . 0 The reader may -I . . . verify that A 0 0-14 - 1 has the form A the n 2 by n 2 68 1/4 1/16 1/64 0 0 0 1/4 0 0 1/16 1/4 ... 0 0 (B-3) 0 • • . . . where 1/16 1/4 i n the lower t r i a n g l e the ( i , j) entry i s 4 to the power -1 - i + j . L" 1 has the form A" A' A" . (B-4) 1 2 3 0 A" A" 1 2 0 0 A" . 1 0 0 0 ... where the block ... A" entry in A (B-5) 4A ( i , j ) = A _ k (i-1,j) + A" C l e a r l y the elements along _k 1 A to the power in AA~ k with -1-i+j. the ( i , j ) we o b t a i n k A 0 A" 2 in position ( i , j ) i s If we compare the ( i , j ) e n t r y (B-6) 0 (i,i) = 4• 4 k k + 1 (i,j). the main d i a g o n a l of A k are 69 Using (B-6) and (B-5) the reader may v e r i f y that (B-7) A " ( i , j ) = k (k+i-j-1 )! k+i-j -/4 , i f i ^ j , and (k-1)! ( i - j ) l ' 0, Note i f i<j. that the e n t r i e s i n the rows are equal t o the e n t r i e s i n the columns under the (B-8) A " ( i , j ) correspondence = A' (n-j + 1 , n - i + 1 ) . k k I n s p e c t i o n of the form of (B-4) shows that then (B-9) L " ( i , j ) = 1 L" (n -j+1,n -i+1). 1 2 2 2. Norms The square grid vectors x of the 2-norm induced on L by the 2-norm on is (B-10) |L| = max{(Lx,Lx)| 2 |X| =1}. 2 Thus |L| = / K where K i s the l a r g e s t eigenvalue of L L . Let z T 2 be an e i g e n v e c t o r of L L with e i g e n v a l u e K. Then T (B-11) | K | = |L^Lz| / | z | < |LT| |L| = | L | J L ! , 1 1 where |L| denotes the induced 1 1 1 1-norm which i s the maximum of the sums of the e n t r i e s i n the columns of L, and |L| w denotes 70 the induced supremum norm which i s the maximum of the sums of the entries i n the rows of L . From (B-1) and (B-2) both a r e 6. Thus K < 36 and (B-12) | L | < 2 6. By a s i m i l a r argument (B-13) | L - ' | < 2 2 IL-'IJL-^. Because of the correspondence (B-9), and we may w r i t e (B-14) | L - ' | < I L - ! . 1 2 1 The reader may v e r i f y by i n s p e c t i o n the sum of the e n t r i e s i n the f i r s t than the sum of the e n t r i e s (B-15) |L-'| " l n column of (B-7) L~ 1 i n any other column. Thus 2 , ,m+l > (m-j):j! m=l 4 j=0 L of (B-4) and j that i s larger CO I4 m+1 m=0 4 00 I — m+2 m=0 2 1/2.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- On local mode analysis in multi-grid methods
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
On local mode analysis in multi-grid methods Reimers, Mark Allan 1983
pdf
Page Metadata
Item Metadata
Title | On local mode analysis in multi-grid methods |
Creator |
Reimers, Mark Allan |
Date Issued | 1983 |
Description | Local Mode Analysis was introduced by Achi Brandt as a heuristic practical tool for determining the expected convergence rate of a Multi-Grid algorithm. It is shown that this tool may be justified rigorously. Furthermore, computations analogous to actual relaxation processes yield results much like Local Mode Analysis predictions. This analysis is also useful for a heuristic understanding of certain relaxations independently of the Multi-Grid context. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2010-04-22 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079365 |
URI | http://hdl.handle.net/2429/24019 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-UBC_1983_A6_7 R44.pdf [ 3.45MB ]
- Metadata
- JSON: 831-1.0079365.json
- JSON-LD: 831-1.0079365-ld.json
- RDF/XML (Pretty): 831-1.0079365-rdf.xml
- RDF/JSON: 831-1.0079365-rdf.json
- Turtle: 831-1.0079365-turtle.txt
- N-Triples: 831-1.0079365-rdf-ntriples.txt
- Original Record: 831-1.0079365-source.json
- Full Text
- 831-1.0079365-fulltext.txt
- Citation
- 831-1.0079365.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079365/manifest