UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

On local mode analysis in multi-grid methods Reimers, Mark Allan 1983

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata

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

ON LOCAL MODE ANALYSIS IN MULTI-GRID METHODS By MARK ALLAN REIMERS B. Sc., The University 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 thesis as conforming to the required standard THE UNIVERSITY OF August © M a r k Allan BRITISH COLUMBIA 1983 Reimers, 1983 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. I t i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. requirements for an advanced degree at the University Department of The University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 DE-6 (.3/81) i i ABSTRACT Local Mode Analysis was introduced by Achi Brandt as a h e u r i s t i c p r a c t i c a l tool for determining the expected convergence rate of a Multi-Grid algorithm. It i s shown that this tool may be j u s t i f i e d rigorously. Furthermore, computations analogous to actual relaxation processes y i e l d results much l i k e Local Mode Analysis predictions. This analysis is also useful for a h e u r i s t i c understanding of certain relaxations independently of the Multi-Grid context. i i i TABLE OF CONTENTS Abstract i i L i s t of Tables v L i s t of Figures . . . . . v i Acknowledgement v i i Chapter 1 - Introduction 1 . 1 The Subject 1 1.2 Outline of this Thesis 4 Chapter 2 - Background 2.1 F i n i t e Difference Solution 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 Multi-Grid (MG) Method 11 2.3 Local Mode Analysis 17 Chapter 3 - A J u s t i f i c a t i o n of Local Mode Analysis 3.1 The Meaning of Local Mode Analysis 21 3.2 Accuracy of Local Mode Analysis 23 3.3 Examples 35 3.4 Comments and Miscellaneous Consequences 41 Chapter 4 - A Heuristic Look at Relaxations using Local Mode Analysis 4. 1 The SOR Method 46 4.2 Alternating Direction Implicit Iteration 51 4.3 Smoothing Property of Relaxations 55 i v Chapter 5 - Discussion 5.1 A p p l i c a b i l i t y of Local Mode Analysis 57 5.2 The Potential of LMA Outside the MG Context..58 5.3 Possible 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 Fourier Components of Gauss-Seidel Relaxation of Laplace's Equation in the Unit Square 38 Table II - Line Relaxation 39 Table III - D i s t r i b u t i v e Relaxation of the Cauchy-Riemann Equations 41 Table IV - Best Values of co 50 v i FIGURES Figure 1 - Orderings: (a) Lexicographic (b) Red-Black.10 Figure 2 - Relation of Fine (f) and Coarse (C) Grids..13 Figure 3 - Numbering of Grid Points 24 Figure 4 - Entries in R1 Corresponding to Grid Points on the Boundary 32 Figure 5 - Contour Lines of |M(0)| for SOR Relaxation.48 ACKNOWLEDGMENT I would l i k e to thank Robert Miura for many encouraging talks, even though my work was outside his f i e l d of research, and Uri Ascher, who f i r s t exposed me to Multi-Grid, and Marc Lemieux, a very f r i e n d l y office-mate. Most of a l l I would l i k e to thank Mary Regester who put up with me during a l l the consummate aggravation of producing this work. 1 CHAPTER ONE INTRODUCTION 1.1 The Subject Numerical techniques for the solution of p a r t i a l d i f f e r e n t i a l equations (PDEs) are finding increasing application in science and industry. Much current research i s devoted to finding fast and robust algorithms for this purpose. E l l i p t i c PDEs arise from steady state problems, and commonly the f i r s t step in their solution is the translation of the continuous problem into a large sparse system of algebraic equations by the use of a f i n i t e difference or f i n i t e element d i s c r e t i z a t i o n . These systems are frequently solved by f i r s t approximating the solution very roughly, and then reducing the error in th i s aproximation by an i t e r a t i v e method (relaxation scheme). In the Multi-Grid (MG) method, relaxation techniques are used, not to eliminate the error in an approximation to the solution, but only to smooth th i s error. (For a description of relaxation techniques and the MG method, the reader may refer to Chapter Two.) Since the p r a c t i c a l introduction of the MG method by Achi Brandt in 1971, i t has shown great early promise and i s now being investigated systematically. A c r u c i a l problem for applications is the estimation of the rate of convergence of a MG algorithm. This requires an estimation of the e f f i c i e n c y of the smoothing performed by the relaxat ion. Local Mode Analysis (LMA) or Local Fourier Analysis was introduced by Brandt (see Brandt(1977)) to r e a l i s t i c a l l y 2 estimate the convergence rate of the smoothing step in a Multi-Grid algorithm. Like the von Neumann s t a b i l i t y analysis for parabolic systems, LMA treats a model system without boundaries, i . e . , on a gri d which i s i n f i n i t e in extent. LMA can be used to obtain good estimates of the (rapid) convergence rate of the MG algorithm as a whole. Brandt has not given a rigorous j u s t i f i c a t i o n of this procedure (on t h i s topic, see Brandt (1982)). Other mathematicians have t r i e d to develop rigourous convergence estimates for MG methods, but they have not found LMA useful for their theories. Most such rigourous convergence estimates have been too conservative (by several orders of magnitude) to be of any p r a c t i c a l value. Therefore LMA remains the p r a c t i c a l tool for r e a l i s t i c estimation and decision 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 of LMA to real problems has been neither rigourously j u s t i f i e d nor disproven. It i s not even clear how the analysis of an i n f i n i t e model system ought to be brought to bear on the f i n i t e system at hand. Thus far mathematical work related to LMA has taken three di r e c t i o n s . Trottenberg and Stueben (1982) carried out a rigourous Fourier analysis of certain special MG situations ( e s s e n t i a l l y D i r i c h l e t problems in a square using Red-Black or Richardson relaxation for smoothing). They derived convergence rates for the smoothing step which turned out to be i d e n t i c a l with convergence rates derived via an extension of LMA for the same problems. They claimed that the identit y of these results gave c r e d i b i l i t y to the application of LMA in more general 3 settings. A second avenue of the mathematical work has been the elaboration of Brandt's o r i g i n a l notion of an i n f i n i t e g r i d model. Hemker (1980a) described in d e t a i l the Fourier analysis of i n f i n i t e g r i d analogues of prolongation and r e s t r i c t i o n operators. Further Hemker (1980b, 1981), Mol (1981), and de Vries (1982) analysed an i n f i n i t e g r i d model of the Incomplete L-U Decomposition (ILU) relaxation process. However, they did not rigourously j u s t i f y the use of these estimates for f i n i t e g r i d problems. It i s worth noting that de Vries (1982) went further than Brandt's o r i g i n a l intent regarding the range of a p p l i c a b i l t y of LMA. de Vries computed the LMA convergence rate estimate for a certain low frequency Fourier component. He then compared th i s with the actual spectral radius of the relaxation operator on a grid whose mesh spacing was determined by the frequency. These two values were close, e s p e c i a l l y for fine meshes, which suggested that LMA might be used to estimate the o v e r a l l performance of relaxation schemes. This idea extends the o r i g i n a l horizons of the method, which was envisaged to be r e a l i s t i c only for high frequencies (see Brandt 1977, 1982). A t h i r d mathematical dir e c t i o n related to LMA i s mentioned by Trottenberg and Stueben (1982), MO1(1981), and Hackbusch (1983). LMA may be interpreted as a rigourous analysis of a f i n i t e discrete problem which arises from a d i f f e r e n t i a l equation problem in a rectangle with doubly periodic boundary conditions. These authors imply that the s i m i l a r i t y between the periodic 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 for the l a t t e r . 4 Hackbusch (1983) went on to relate the periodic and D i r i c h l e t cases e x p l i c i t l y . He treated the MG solution of singular perturbation problems in one dimension, which are to be smoothed by Gauss-Seidel relaxation. He derived an equation involving the periodic relaxation operator, the D i r i c h l e t relaxation operator, and a rank one perturbation matrix. He used t h i s to l i n k the size of the D i r i c h l e t relaxation operator to the size of the periodic operator in a special norm related to the high frequency components. This estimate is further used in a rigourous estimate of the convergence rate of the MG algorithm for t h i s problem. Hackbusch i s the only author who has made use of LMA type analysis in a rigourous treatment of MG. 1.2 Outline of t h i s Thesis In t h i s work I treat two-dimensional D i r i c h l e t problems on grids which approximate a r b i t r a r y domains. The discrete systems a r i s i n g from these problems are relaxed by any of a large class of i t e r a t i v e schemes to be described. In Chapter Three, I interpret LMA as an approximate description of the relaxation of Fourier components on the f i n i t e grids. This description i s asymptotically correct in the l i m i t of very fine grids. To show t h i s , I perturb the matrix describing the actual relaxation to obtain a matrix which operates on Fourier components exactly as described by LMA. In the case of rectangular domains, this perturbed operator i s i d e n t i c a l with the relaxation operator for a problem with periodic boundary conditions. I show that t h i s perturbation i s asymptotically small in norm, and hence validate LMA on fine grids. This 5 v a l i d a t i o n i s i l l u s t r a t e d w i t h s e v e r a l examples worked out in d e t a i l , to show the a c c u r a c y of LMA. LMA i d e a l i z a t i o n s for s e l e c t e d F o u r i e r components are compared wi th a c t u a l r e l a x a t i o n s over a range of problems and meshes. In Chapter F o u r , I i n v e s t i g a t e , as d i d de V r i e s (1982), the p o s s i b i l i t y of u s i n g LMA to d e s c r i b e the o v e r a l l convergence r a t e of r e l a x a t i o n s o u t s i d e the MG c o n t e x t . The u s e f u l n e s s of LMA in p r e d i c t i n g the optimum r e l a x a t i o n parameter for the SOR scheme i s i l l u s t r a t e d over a range of problems and meshes. LMA i s a l s o used to h e u r i s t i c a l l y e x p l a i n the e f f i c i e n c y of the ADI method a p p l i e d to L a p l a c e ' s e q u a t i o n and the d e g r a d a t i o n i n performance for non-symmetric o p e r a t o r s . Chapter F i v e i s a d i s c u s s i o n of some of the ideas beh ind t h i s work, as w e l l as a proposed e x t e n s i o n . I hope t h i s work w i l l b u i l d c o n f i d e n c e i n LMA and suggest i t s p o s s i b l e uses as a t h e o r e t i c a l t o o l . 6 CHAPTER TWO 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 E q u a t i o n s In p r a c t i c e the 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 E q u a t i o n s (PDEs) i s f r e q u e n t l y o n l y f e a s i b l e by n u m e r i c a l methods. Commonly the (continuum) PDE i s t r a n s l a t e d i n t o a system of d i s c r e t e e q u a t i o n s ( d i s c r e t i z a t i o n ) . These e q u a t i o n s may be s o l v e d e i t h e r by a d i r e c t method or by an i t e r a t i v e a l g o r i t h m . A r e c e n t development has been the M u l t i - G r i d (MG) method i n which t h e r e i s i n t e r a c t i o n between the i t e r a t i v e s o l u t i o n of the e q u a t i o n s and the d i s c r e t i z a t i o n i t s e l f . The MG method has proven q u i t e e f f e c t i v e i n t r e a t i n g many d i f f e r e n t t y p e s of problems. In t h i s work one f a c e t of the MG t e c h n i q u e , L o c a l Mode A n a l y s i s (LMA), i s j u s t i f i e d and a p p l i e d t o the s t u d y of i t e r a t i v e methods. C o n s i d e r a t w o - d i m e n s i o n a l f i n i t e domain fl and an 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 o p e r a t o r L. We w i l l c o n s i d e r o n l y l i n e a r o p e r a t o r s L. The problem i s t o f i n d a f u n c t i o n u on 0 which s a t i s f i e s (2-1 ) Lu = f ( x , y ) , where f i s a p r e s c r i b e d f u n c t i o n on fl and u and/or i t s d e r i v a t i v e s s a t i s f y some p r e s c r i b e d l i n e a r boundary c o n d i t i o n s on 9 f l , t he boundary of fl. Our 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 t h i s problem proceeds as f o l l o w s . P l a c e a d i s c r e t e s e t of N p o i n t s on fl, spaced 7 evenly in x and in y to form a rectangular g r i d . Frequently the same spacing i s used in both co-ordinate directions and thi s spacing i s c a l l e d the mesh size h. There are several ways to approximate the boundary of a non-rectangular domain. We w i l l discuss one of them in Chapter Three. The g r i d points are ordered so that numbers assigned to them form the components of a vector U which we c a l l a g r i d vector. Now we approximate the continuum equation for u by a system of N algebraic equations which determine the values of U at the N gri d points. The components of U are to approximate the values of u on the corresponding gri d points. We construct an equation at each i n t e r i o r grid point (x,y) by approximating the derivatives that appear in Lu(x,y) by f i n i t e differences, e.g., (2-2) 9 2u/9x 2(x,y) = (1/h 2)(u(x+h,y)-2u(x,y)+u(x-h,y)) +0(h 2). These difference approximations to the derivatives are used to build up an approximation L*U(x,y) to Lu(x,y). The discrete equation corresponding to a grid point (x,y) i s then (2-3) L*U(x,y) = f ( x , y ) . Thus for each i n t e r i o r g r i d point we have an equation involving the component of U corresponding to that point and the components of U corresponding to neighbouring points. At grid points (x,y) on or near the boundary of fi, similar approximate equations are set up, but now they must take into account the boundary conditions imposed on u. We obtain a 8 linear N x N system of equations AX = b, whose solution X = U approximates the solution u of the PDE. The approximation generally becomes better as N increases. For further information on the d i s c r e t i z a t i o n process see Young (1971 ) . Now we .must find an e f f i c i e n t algorithm for solving AX = b. This i s the most costly part of the f i n i t e difference solution of the PDE. For the large systems (N > 10 3) that are desirable to ensure accuracy, the usual Gauss elimination (GE) algorithm turns out to be too space- and time-consuming. Other d i r e c t methods, based on the Fourier transform or on reduction (transformation of the large sparse system into a smaller dense system), may be used to obtain exact solutions of some special discrete systems. These methods are faster than GE and the development of large core computers has made them p r a c t i c a l for the solution of larger systems. In general large systems can be handled more quickly i f they are solved only approximately using i t e r a t i v e methods (or relaxation techniques). A common class of relaxation schemes are the fixed point i t e r a t i o n s . It i s these methods with which we s h a l l be concerned. A fixed point i t e r a t i o n for the system AX = b obtains a matrix B which approximates A and which i s e a s i l y inverted. Then the function (2-4) h(X) = X - B"1AX + B" 1b has a fixed point X = U the desired discrete solution. Choose any grid vector X° and construct a sequence of i t e r a t e s X°, X 1 = h(X°), X 2 = h(X 1),... Let e be the 9 difference between an it e r a t e X and U, e = X - U. Then the difference e' between the next approximation X' and U i s given by (2-5) e' = (I - B"1 A)e. If the spectral radius of the amplification matrix I - B"'A i s less than 1 then c l e a r l y the errors in the sequence of iterates w i l l go to 0 and the sequence w i l l converge to X = U. There are simple s u f f i c i e n t conditions which are known to guarantee convergence of such an i t e r a t i o n scheme (see Young, 1971). For example: l e t A be s p l i t into A = B - (B - A); i f B"1 and B - A are non-negative (a "regular" s p l i t t i n g ) , then the i t e r a t i o n w i l l converge. Most of the common fixed point i t e r a t i o n s produce regular s p l i t t i n g s as outlined below. In the Jacobi (or simultaneous displacement) method, B i s the diagonal of A (see Young, 1971). In the Gauss-Seidel (G-S) method, B i s the lower triangular part of A (including the diagonal). For l i n e relaxation, B is the lower tr i a n g l e and one of the non-zero superdiagonals of A. In Successive Over-Relaxation (SOR), B i s D + a>L where D i s the diagonal and L the lower triangle not including the diagonal of A. Incomplete L-U Decomposition (ILU) (see de Vries, 1982) finds matrices L, U, and R such that L i s lower triangular, U is upper triangular and R i s "small" in some sense and such that A = LU + R; B i s LU and i s therefore e a s i l y inverted. For sparse matrices A a r i s i n g from PDEs L , U and R should be sparse. The Alternating-10 Direction Implicit (ADI) method is not based on a simple s p l i t t i n g but rather on a sequence of two s p l i t t i n g s . It w i l l be discussed further l a t e r . The Conjugate-Gradient (CG) it e r a t i o n i s not of the type described above. The p r a c t i c a l implementation of a relaxation scheme en t a i l s more than just a knowledge of the s p l i t t i n g on which i t i s based. We w i l l discuss the implementation of the simplest widely used relaxation, the Gauss-Seidel scheme. This method w i l l be used frequently in examples. The gri d points are ordered in some manner. In that order each approximate solution value i s adjusted to s a t i s f y the one corresponding discrete equation. Depending on the ordering of the points, the determination of the value of the next iterate at a grid point may involve using the newly determined values at one or more neighbouring points as well as values of the 1 1 --5-1 5 — 8 I I I I 2-12--6-16 . I I I I 9 — 3-1 3 ~ 7 I I I I 1-10--4-14 (a) (b) Figure 1. Orderings: (a) Lexicographic and (b) Red-Black current i t e r a t e . The most common orderings are lexicographic and red-black (see Figure 1). Lexicographic relaxation starts at one corner of the grid and proceeds along one grid l i n e . When the end of the l i n e i s reached, relaxation starts again at the f i r s t position of the next l i n e , and so on. u n t i l a l l points on a l l grid l i n e s have been modified. This constitutes 13-14-15-16 I I I I 9-10-11-12 11 one i t e r a t i o n . Red-Black relaxation divides the grid points into two ordered subsets related to each other l i k e the red squares and the black squares of a checkerboard. A l l the points of the red subset are relaxed in their order, and then a l l the points of the black subset. 2.2 The Multi-Grid (MG) Method The best of the relaxation schemes mentioned above are better than Gauss Elimination, but they are s t i l l quite i n e f f i c i e n t , and therefore s t i l l c o s t l y . Typ i c a l l y on the f i r s t few passes of a relaxation scheme, the error in the approximate solution (as measured by the size of the residuals r = b - AX) i s decreased substantially (on the order of 10% to 50%). However, subsequent passes are not so e f f i c i e n t and long before the desired tolerance on the residuals i s reached, they are being reduced by only a small f r a c t i o n on each i t e r a t i o n . Very many relaxation sweeps are thus required to solve the equations. For most of these schemes the theory shows that the spectral radius of the amplification matrix i s of order 1 -k 0(h ), where h is the mesh size. Depending on the scheme used and the s p e c i f i c problem, the power k of h varies from 3/2 to 5/2. In practice the convergence rate asymptotically approaches the spectral radius, since the components of the error along the eigenvectors of the amplification matrix, with the largest eigenvalues, gradually become dominant. That i s , the quickly converging error components (small eigenvalues) are almost eliminated by the f i r s t few it e r a t i o n s , but the slowly converging components are worn down only gradually and 1 2 thus stand out after many passes. If we are to' more e f f i c i e n t l y solve the problem we need to find some other way to eliminate these components. What do the slowly converging error components look lik e ? The residuals after several applications of nearly any relaxation scheme are smoother than the residuals for the i n i t i a l guess X°. That i s , there i s more co r r e l a t i o n between the residuals at neighbouring points, and fewer sign changes as one moves along a row of the gr i d . Thus i t i s natural to conjecture that the slowly converging error components are also smooth. Numerical computation of the eigenvalues of the amplification matrix for simple model problems bears out th i s i n t u i t i o n . The largest eigenvalues are generally associated with eigenvectors, whose values at neighbouring gri d points are strongly correlated. On the other hand, the eigenvectors corresponding to small eigenvalues vary quickly along g r i d l i n e s . These observations suggest that the relaxation schemes might be better understood by decomposing the error e = X - U in an it e r a t e X, into a f i n i t e Fourier-type ser i e s . A he u r i s t i c analysis suggests, and this thesis makes i t rigourous, that while most relaxation schemes are e f f i c i e n t at reducing the high frequency Fourier components of error, they are not e f f i c i e n t at reducing the low frequency components. Thus we cannot in general eliminate the slowly converging error components from one relaxation scheme simply by using a di f f e r e n t scheme. The idea behind the MG algorithm i s to use the relaxation methods to do only what they do well - reducing the high 13 frequency components of the error, and to find some other method to eliminate the low frequency components. The problem then is to find a way to eliminate a "smooth" or slowly varying error in an approximate solution. This may be approached as follows. Suppose several relaxation sweeps have been done on the discrete system AX = b obtained from (2-1), y i e l d i n g an approximate solution U 1, and suppose that the residuals r 1 = b - AU1 are smooth. The determination of the exact discrete solution U is equivalent to the determination of the smooth error vector e 1 = U 1 - U which solves Ae 1 = r 1 . Consider Ae 1 = r 1 as the d i s c r e t i z a t i o n of a PDE problem Le = r in J2, where e i s to s a t i s f y homogeneous boundary conditions. We might expect then that e w i l l be smooth r e l a t i v e to the grid we are using for Lu = f. Then an accurate solution for e may be obtained with a much coarser g r i d . Thus we set up a f i n i t e difference system A*e* = r* by d i s c r e t i z i n g Le = r on a coarser grid (see Figure 2). The solution e* of thi s system should be a good —-C f C f I I I I — f f f f I I I I — c f c f Figure 2. Relation of Fine (f) and Coarse (C) Grids approximation on the points of the coarse g r i d to the smooth function e which i s approximated on the fine grid by e 1. Since e i s smooth, values of e on the fine g r i d points may be obtained accurately simply by interpolation from the coarse 14 gr i d points. Thus a good approximation for e 1 can be obtained from e*, and this smooth error can be almost eliminated. Why would this procedure (called Coarse Grid Correction (CGC) ) be more e f f i c i e n t than simply eliminating e 1 by many more relaxations on the o r i g i n a l fine grid? It i s because each relaxation of A*e* = r* i s much less expensive than a relaxation of Ae 1 = r 1 (which i s equivalent to a further relaxation of AX = b). In addition fewer i t e r a t i o n s w i l l be required to solve the coarse gr i d system to a given tolerance than w i l l be required to solve the fine g r i d system to the same tolerance. Thus a smooth error can be eliminated more quickly t h i s way. Of course the solution of the CGC problem may be accomplished in the same way; a few relaxation sweeps and then transfer to a yet coarser g r i d . Thus a recursive multiple-grid algorithm i s devised. On each g r i d the rough component of the error i s reduced by relaxation and the smooth ( r e l a t i v e to the resolution of that grid) component is reduced via the next coarser g r i d . These ideas may be implemented in several ways. The algorithm sketched below is referred to as an adaptive V-cycle correction scheme (see Brandt, 1977, 1982) Let the sequence of 1 2 n successively finer grids be G , G ,... ,G ; the desired n solution i s a g r i d vector on G . Let the mesh sizes be 1 n m m+1 h > ... > h ; usually h = 2h . There are operators which transform gr i d vectors on one g r i d to g r i d vectors on another k k-1 g r i d . Prolongation I acts on g r i d vectors on G and k-1 smoothly interpolates to y i e l d values at a l l the g r i d points k k-1 in G . The reverse procedure, r e s t r i c t i o n I , interpolates k 1 5 k-1 values at g r i d points in G using the data supplied by the k values of a g r i d vector on G . Both operations find a good representation on the target grid of a function which i s represented by a vector on the o r i g i n a l g r i d . k k k k We solve a problem A X = b on each grid G as follows. The system i s relaxed several times u n t i l the error k in the last approximate solution Z i s judged to be "smooth". This i s determined by measuring the rate of decrease of the residuals in some suitable norm. As long as the error has large non-smooth components, each relaxation should e f f i c i e n t l y reduce the residuals. An e f f i c i e n c y c r i t e r i o n K (0 < K < 1) must be set. As soon as the size of the residuals after a relaxation i s no less than K times the size before, then the error i s judged to be smooth. Ty p i c a l l y less than 5 i t e r a t i o n s are done. After these relaxation sweeps we set up the CGC problem k-1 k k k k k-1 k-1 k on G . Let r = b - A Z , and b =1 r . The system k-1 k to be solved using grid G is then k-1 k-1 k-1 (2-6) A X = b k-1 k-1 The operator A used for G may be taken to be the k-1 f i n i t e difference approximation on G to the underlying d i f f e r e n t i a l operator L , or, (more suitably for theoretical 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 s u f f i c i e n t number of relaxation sweeps. Using either method the amount of work done i s small on the coarsest g r i d . On the other grids (2-6) i s solved by smoothing steps and a CGC. It 1 6 is only necessary to solve the CGC problem to an accuracy commensurate with the size of the remaining non-smooth error k on G k-1 k-1 k-1 When a solution X = U i s obtained on G then l e t k k k-1 k k e = 1 U . An improved solution Y on G is obtained via k k k-1 k Y = Z - e. This eliminates most of the smooth part of the k k k error in Z . Since the error Y - U is now dominated by i t s non-smooth component, t h i s error may be further e f f i c i e n t l y reduced by relaxation steps. Thus the non-smooth component of the error in the n o r i g i n a l approximation on G has been cut down by relaxations, and the smooth part has been almost eliminated by the CGC. The Multi-Grid procedure outlined above is complicated and although the MG idea seems promising, there i s s t i l l controversy among those who have applied i t . Brandt f i r s t implemented the scheme in the early seventies (see Brandt, 1977) and suggested that i t was an e f f i c i e n t (rapidly converging) algorithm for numerical solution of PDEs. He promoted i t vigorously and many d i f f i c u l t problems were solved e f f i c i e n t l y via MG. However, slow convergence of the method for other problems has led to much doubt regarding i t s robustness. There i s l i t t l e theory which can help decide whether the poor performance of some MG applications i s the result of inappropriate implementation, or whether i t r e f l e c t s a fundamental l i m i t a t i o n of the method. Convergence proofs for the MG algorithm (given already in the early s i x t i e s by Bakhvalov (see Brandt, 1977)) show that a reduction of one 1 7 order of magnitude in the o v e r a l l error can be achieved in 0(N) operations (N i s the number of grid p o i n t s ) 1 . However the constant in this 0(N) i s so large as to have no bearing on convergence rates in practice. Brandt's claim that the MG algorithm w i l l converge quickly i s based not on rigorous theory but rather on a h e u r i s t i c l o c a l Fourier analysis of relaxation (and also r e s t r i c t i o n and prolongation). With suitable r e s t r i c t i o n s and prolongations, this analysis suggests that MG w i l l converge rapidly provided that the non-smooth error components on each grid are substantially reduced by relaxation. In t h i s thesis we w i l l be concerned with the application of t h i s Local Mode Analysis (LMA) to relaxations. 2.3 Local Mode Analysis For a given problem and a s p e c i f i e d relaxation scheme LMA i s used to estimate the rate of reduction of the various Fourier modes which comprise the error. Brandt's o r i g i n a l idea was that relaxation at one g r i d point strongly affects the relaxation only at nearby gri d points, i . e . , "relaxation i s a l o c a l process" (see Brandt, 1977). Thus he imagined a relaxation acting on an i n f i n i t e g r i d and asked how relaxation would reduce a p a r t i c u l a r Fourier component. We w i l l i l l u s t r a t e LMA by considering Gauss-Seidel relaxation of Laplace's equation in some domain. The equation l i n k i n g the new values X' to the old values X during a sweep over lexicographically ordered g r i d points i s 1 Straightforward use of a good relaxation w i l l result in convergence in 0(N 2) (or possibly N to the power 3/2) operat ions. 18 (2-7) 4X'(i,j) -X'(i-1,j) -X'(i,j-1) -X(i+1,j) -X(i,j+1) = 0. Note that when the point ( i , j ) i s reached, newly computed values X' at points ( i - 1 , j ) and ( i , j - l ) are already in place. The desired solution U also s a t i s f i e s this equation and upon subtracting we obtain the error transformation equation, (2-8) 4e'(i,j) - e ' ( i - 1 , j ) - e ' ( i , j - 1 ) -e(i+1,j) - e ( i , j + l ) = 0. Brandt (1977) suggested that one can "analyse i t [2-8] in the i n t e r i o r of the grid by ( l o c a l l y ) expanding the error in Fourier s e r i e s . ... Thus to study the d = (d*,d2) Fourier component of the error before and after the relaxation sweep, we put (2-9) e ( i , j ) = A ( 0 ) e x p ( i ( 0 1 i + 0 2 j ) ) ( i = /-I) and e ' ( i , j ) = A*(0)exp(t(0 1i + 0 2 j ) ) . " Here 01 and 82 must be interpreted as wave numbers in the x and y d i r e c t i o n s , respectively. Using the id e n t i t y exp( c0(k±1 )) = exp( t0k)exp(± it9) the error transformation equation (2-8) becomes (2-10) (4 - exp(-i0 1) - exp(-i0 2) )A'(0) + (-expd .0 1 ) - exp(t0 2) )A(0) = 0. The convergence rate of the 9 \ 8 2 component i s then given by the r a t i o 19 (2 -11) M ( 0 ) = A' ( 0 ) / A ( 0 ) = (- exp( <,en - exp( id2) ) / ( 4 - e x p ( - i 0 1 ) - e x p ( - i 0 2 ) ) . Brandt expected this estimate of the convergence rate to be good for the high frequencies (which are the only ones of interest in the MG context where the role of relaxation i s smoothing only), because for these frequencies the interactions of the solution value at an i n t e r i o r point with the boundary conditions and with far away points are small. For the usual case where the mesh size in the finer g r i d i s half the mesh size in the coarser grid, the boundary between "high frequencies" and "low frequencies" i s taken as max ( | 9 1 | , | 92 \ ) = 7r/2 , since the ( 0 1 , 0 2 ) component i s indistinguishable from the (9 1 -TT, 92-IT) component i f only values at every second grid point in the fine g r i d are considered. The exact form of u.(9) i s usually ignored and an estimate of the convergence rate of the high frequency components (called the smoothing factor) i s taken as (2 -12) u = max{|/z(0)| ; max ( | 9 1 | , | 92 | ) > TT /2}. This e f f e c t i v e rate of smoothing can be used to estimate the e f f i c i e n c y of the algorithm as a whole, and to estimate the t o t a l work required for a MG solution. In practice i t estimates accurately the high e f f i c i e n c y of many MG implementat ions. Despite i t s p r a c t i c a l success so far, many people who have worked on MG methods have found the LMA procedure 20 dubious. Brandt himself has never come forth with a j u s t i f i c a t i o n for 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 w i l l consider the v a l i d i t y of Local Mode Analysis (LMA) for problems with D i r i c h l e t boundary conditions. We w i l l take the point of view that LMA should approximately describe the transformation of a Fourier mode error under relaxation. We w i l l show that this approximation becomes better as the grid i s refined. We w i l l give several concrete examples of how good the approximation i s in p r a c t i c a l situations. 3.1 The Meaning of Local Mode Analysis. In explanations of Local Mode Analysis (e.g., Brandt, 1977) the error is supposed to be written as a " l o c a l Fourier series" near a gr i d point ( i , j) . This notion of " l o c a l Fourier series" arises from an attempt to make LMA look plausible in the case of general boundary conditions and variable c o e f f i c i e n t s , by ignoring the boundary conditions and 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 must . f i r s t dispose of the notion of a " l o c a l Fourier series" in order to make a rigourous, global, interpretation of the analysis. We f i r s t discuss this in the context of a one-dimensional continuum, and l a t e r we w i l l discuss the multi-dimensional discrete case by analogy. Consider a function f(x) defined on some i n t e r v a l I containing a point x = a. Consider what would be required of a " l o c a l Fourier s e r i e s " (3-1) 2 L d(0)exp(i0(x-a)) e 22 for f(x) near x = a. We might require i t to converge in a neighbourhood of x = a and to describe the behaviour of f close to a, perhaps converging at a to the value f ( a ) . But these ' l o c a l ' requirements actually come from the idea of a Taylor series. 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 * a+c ma+c (3-2) d(0) = l e x p ( - 0 i x ) £ ( x ) d x /J |exp(-10x)| 2dx. J a-c / ' a-c Note that (3-3) lim d(0) = exp(- i 0 a ) f ( a ) c-K) which is independent of the behaviour of f(x) near a, so that the l i m i t of d(0) as c goes to 0 is not a determination of the l o c a l importance of the 0-component. We must then pick some c > 0, but which one? This brings us to the heart of the problem. A Fourier series can be defined only for a function and an i n t e r v a l , not a function and an in d e f i n i t e neighbourhood of a point. In the case of a discrete point set in real N-space, e.g., a g r i d covering a domain 1, consider a g r i d vector v and a s p e c i f i c point a. The c o e f f i c i e n t of a Fourier mode in a 1 The discrete Fourier modes are then vectors with components equal to the continuum Fourier functions evaluated at points of t h i s ordered set. 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 vectors with summation over a l l grid points. <v,w> = £v(i,j)w(i,j) 23 " l o c a l expansion" of v near x = a w i l l depend on which points are considered to form the neighbourhood of a. This neighbourhood i s not well-defined. It is also not clear which discrete frequencies are to be selected for the expansion. A discrete Fourier series can only be defined for values of a gri d vector corresponding to a s p e c i f i c rectangular subset of the g r i d . There i s then, no unique way to form a " l o c a l Fourier s e r i e s " in the usual sense that ' l o c a l ' objects are defined in mathematics. The concept of " l o c a l Fourier series" i s r e a l l y a conflation of the idea of a Fourier series and that of a Taylor s e r i e s . While the existence of such a " l o c a l Fourier series", and of " l o c a l Fourier modes", i s suggested by the language in some of the papers written on MG methods, th i s is misleading. We suggest that LMA be reinterpreted as a l o c a l analysis of the action of a relaxation scheme on Fourier modes over the whole gr i d . We w i l l show that the a t t r i t i o n rate of errors having the form of a Fourier component (over the whole grid) i s accurately computed by LMA for fine meshes. The interest of these special forms of error i s only that they w i l l be used to give an acceptable interpretation of LMA. 3.2 Accuracy of Local Mode Analysis We develop an interpretation of LMA and j u s t i f y i t in the following s i t u a t i o n . We consider a second order linear e l l i p t i c PDE problem in two dimensions with constant c o e f f i c i e n t s of the form (3-4) Lu = a9 2u/9x 2 + b3 2u/9y 2 + c3u/3x + d9u/9y + eu = f 24 in a domain fl with smooth boundary 3fl. The solution u is assumed to s a t i s f y D i r i c h l e t boundary conditions on 3fl. A grid G covers fl and without loss of generality we assume that the g r i d spacings in the x di r e c t i o n and in the y di r e c t i o n are equal. We consider zero-order approximations to the boundary and to the boundary conditions. That i s , the domain i s approximated by l i n e segments along grid l i n e s , and the D i r i c h l e t data are assigned to gri d points on these l i n e s . Grid points w i l l be doubly indexed ( i , j ) (in 2D) as i f the gri d were part of a larger rectangular grid running p a r a l l e l to the co-ordinate axes (see Figure 3). The f i r s t (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 ) Figure 3 - Numbering of Grid Points index i numbers the values of x and j numbers the values of y. The equations are doubly indexed ( i , j ) l i k e the grid points, so that equation ( i , j) involves values of the solution at points ( i , j ) , ( i - 1 , j ) , (i+1,j), ( i , j - l ) , ( i , j + l ) ( i f a l l these are in the gri d i n t e r i o r ) . This notation makes discussion of LMA easier. If the usual centered differences are used to approximate the derivatives in t h i s PDE, then only the values at a point and i t s four nearest neighbours are involved in the difference equation corresponding to that point. Thus a " f i v e point 25 formula" may be used in the i n t e r i o r of fl. A system of f i n i t e difference equations i s obtained, AX = f, where each equation corresponds to a gr i d point.. In case the equation at a point involves boundary points, (at most three), the fixed values at those points are to be included in the grid vector f. We w i l l suppose that the equations have been multiplied through by h 2, where h i s the mesh si z e . We w i l l j u s t i f y LMA for linear relaxations of AX = f which meet the following conditions. F i r s t the relaxations must be of SPLITTING TYPE, i . e . , the c o e f f i c i e n t matrix A i s rewritten as A = L - R with L non-singular. We fi n d the next iterate X' from the current i t e r a t e x by solving (3-5) LX' = RX + f. This s p l i t t i n g need not be regular (L~ 1 and R non-negative), and i t need not be the same at each i t e r a t i o n . Secondly, we w i l l only discuss s p l i t t i n g s for which each s p l i t equation involves only the same fi v e points as the o r i g i n a l equation. This means (3-6) a(m,n) =0 => l(m,n) = r(m,n) = 0. F i n a l l y , in order that LMA i t s e l f be feasible we must require that every l i n e in LX' = RX + f which corresponds to a point which i s not next to a boundary point i s of the form 26 ( 3 - 7 ) l 1 X ' ( i , j ) + l 2 X ' ( i - 1 , j ) + l 3 X ' ( i , j - 1 ) + 1 " X ' ( i + 1 , j ) + 1 5 X ' ( i , j + 1 ) = r ' X l i , ] ) + r 2 X ( i - 1 , j ) + r 3 X ( i , j - 1 ) + r'X(i+l,j) + r 5 X ( i , j + 1 ) + f ( i , j ) , where 1 1 , . . . , 1 5 , r 1 , . . . , r 5 are constant c o e f f i c i e n t s . This condition w i l l be subsequently referred to as 'homogeneity'. Many common linear relaxation schemes s a t i s f y these conditions, including Jacobi, Gauss-Seidel, SOR, l i n e versions of these, and ADI (see Chapter Two). One i t e r a t i o n of ADI consists of two steps such as are considered here, with L and R interchanged. The ILU method does not s a t i s f y the homogeneity condition and cannot, s t r i c t l y speaking, be analysed by LMA. In practice an i d e a l i z a t i o n of the ILU method is analysed t h i s way, but "we w i l l not j u s t i f y that here. Red-Black and Zebra relaxations do not s a t i s f y these conditions and i t is known (S. McCormick and K. Stueben, private communication) that LMA applied simple-mindedly yi e l d s misleading estimates of smoothing rates in these cases. In j u s t i f y i n g LMA we w i l l consider the reduction of error from one iterate X to the next iterate X' which i s given by (see notation in Chapter Two) ( 3 - 8 ) Le' = Re. It i s too complicated to consider an ar b i t r a r y error e, and write a Fourier series for i t (which amounts to taking a f i n i t e Fourier transform), and then ask what happens to the values of the Fourier transform under a relaxation step. Rather we suppose that the error consists of a multiple of 27 only one Fourier component \p(8) , as defined below, and show that af t e r a relaxation step the new error i s U(8)\}J{8) plus something small. On a grid with the index convention, the value of a discrete exponential Fourier component, \p(8) (with frequencies 8 = ( 0 1 , 0 2 ) , in two dimensions), at a grid point ( i , j ) i s given by ( 3 - 9 ) .//(0)(i,j) = e x p { i ( i 0 1 + j 0 2 ) } . We need only consider -n < 8},82 <n. As mentioned in Chapter Two, LMA treats the Fourier component ^(8) as an eigenvector of the amplification matrix of the relaxation. As noted in Chapter One, t h i s i s only the case for certain relaxations of problems with periodic boundary conditions on a rectangle. However, as w i l l be seen, for many relaxations, the Fourier components are in fact "close to being eigenvectors", because they are eigenvectors of an idealized relaxation that i s "close" to the actual relaxation. The idealized relaxation w i l l be i d e n t i c a l with the actual relaxation except at points next to a boundary point. Since the proportion of such points in the grid becomes smaller as the mesh becomes f i n e r , the idealized relaxation gives a good approximation to the actual relaxation in the l i m i t of fine grids. The estimate of the accuracy of LMA which w i l l be developed makes such arguments rigourous. For a s p e c i f i c Fourier component, \p(8) given by ( 3 - 9 ) , we w i l l construct a system L*\p* = R*\p(8), similar to L\//' = R \ M # ) , for which 4/(8) i s an eigenvector of both L* and R*. F i r s t note that every l i n e in L\p' = R\p(8) which 28 corresponds to a point which i s not next to a boundary point is of the form (3-10) l V ( i , j ) + l V ( i " 1 , j ) + 1 3 V>' ( i , j - 1 ) + I V ( i + 1 , j) + I V ( i , j + 1 ) = r V 0 ) ( i , j ) + r V S ) (i-1 , j) + r V 6) ( i , j-1 ) + r V 0 ) ( i + 1 , j) + r V 0 ) ( i , j + 1 ) , (cf. (3-7)). Using (3-9) the right hand side s i m p l i f i e s to (3-11) ( r 1 + r 2 e x p ( - i 0 M + r 3 e x p ( - i 0 2 ) + r"expU0 1) + r 5exp ( i 0 2 ) )tf ( e)(i,j) = r * i / / ( e ) ( i , j ) . Now i f \p(8) i s substituted for \p' on the l e f t side of (3-10) then i t s i m p l i f i e s to (3-12) ( l 1 + l 2 e x p ( - i 0 1 ) + l 3 e x p ( - i 0 2 ) + l a e x p ( i 0 1 ) + l 5 e x p ( i 0 2 ) )<//(0)(i,j) = l * ^ ( 0 ) ( i , j ) . Now we construct R* so that values of [ R*\//( d) ] ( i , j ) which correspond to points ( i , j ) next to a boundary are equal to r*i//( 8) ( i , j ) , l i k e values corresponding to other i n t e r i o r points. We construct L* likewise so that for every i n t e r i o r point ( i , j ) [L*\p(8) ] ( i , j) = l * ^ ( 0 ) ( i , j ) . We do th i s by including an appropriate multiple of e x p ( ± i # 1 ) or e x p ( ± t # 2 ) , for every equation involving a boundary point, to 29 make i t l i k e the equation for points away from the boundary 2. S p e c i f i c a l l y , suppose we are constructing an equation at a l e f t edge point ( i , j ) where a term 12\^' (i-1 , j ) and a term r2\//( 8) ( i-1 , j ) are missing in (3-10). Then to form the ( i , j ) equation in L*\J/* = R*\p(8) , add l 2 e x p ( - i d 1 ) i//* ( i , j ) and r 2exp( - c 9 1 ) \//( 8) ( i , j ) to the l e f t and right sides, respect i v e l y . 3 If we denote these extra elements along the diagonal by L1 = L* - L and R1 = R* - R, then LI i s 0 except for diagonal entries corresponding to boundary points. LEMMA 1: The Fourier component ip{6) is an eigenvector of the system ^L*e = R*e with eigenvalue M(#). PROOF: If ( i , j ) is an i n t e r i o r grid point the ( i , j ) entries in L*\jj(8) and L\p{8) are both equal to ( 6) ( i , j ) . At a point next to a boundary point the additional entries in L* are chosen to make L*\//( 6) ( i , j ) equal to l * i / > ( 8) ( i , j ) . Sim i l a r l y , for i n t e r i o r points and points near the boundary, the ( i , j ) entry in R*4/{8) i s r * t f ( 0 ) ( i , j ) . If we do a LMA and put 2 If Q i s a rectangle, and the values of 0 1, 82 are r e s t r i c t e d to those values that would occur in a Fourier transform, then the matrices L* and R* describe the relaxation of a problem with periodic boundary conditions. However the construction used here i s not r e s t r i c t e d to rectangular domains. 3 E.g., for the Gauss-Seidel relaxation of Laplace's equation,the l e f t edge equation which has the form 4e' ( i , j ) - e'(i , j - 1 ) = e(i+1,j) + e ( i , j + l ) , becomes, (4 - e x p ( - i f l 1 ) ) e * ( i , j ) - e * ( i , j - l ) = e ( i + 1,j) + e ( i , j + 1). 30 (3-13) v = A ( 0)exp( i ( 0 1 i + 0 2 j ) ) , and v' = A ' ( 0 ) e x p ( i ( 0 1 i + 0 2 j ) ) , where v and v' are the errors before and after relaxation, then (3-14) A ' ( 0 ) ( 1 1 + l 2 e x p ( - t 0 1 ) + l 3 e x p ( - t 0 2 ) + l " e x p ( i 0 1 ) + l 5 e x p ( i 0 2 ) ) e x p ( i ( 0 1 i + 0 2 j ) ) = A ( 0 ) ( r 1 + r 2 e x p ( - i 0 1 ) + r 3 e x p ( - i 0 2 ) + r " e x p ( i 0 1 ) + r 5 e x p ( i 0 2 ) ) e x p ( i ( 0 1 i + 0 2 j ) ) . Thus, (3-15) /x ( 0 ) = A ' ( 0 ) / A ( 0 ) = r * / l * and therefore (3-16) L * M ( 0 ) < / > ( 0 ) = R*\p{9). QED Now l e t us consider the difference between i / / * = M ( 0 ) * M 0 ) which solves (3-17) L*\p* = R*tf(0) and the grid vector which i s the actual solution of (3-18) LiT = R\p(6) . The difference between L and L* occurs only at points next to a boundary point. For thi s difference to be small then, we need a lemma to guarantee a small number of grid points next to a boundary point. LEMMA 2: Let the domain 0 be f i n i t e with piecewise smooth boundary 3S2. Then for h s u f f i c i e n t l y small there i s a number a, such that for any grid of mesh size h, the number N(b) of 31 grid points in the i n t e r i o r of fl, one of whose four nearest neighbours l i e s on or outside 9fl, is bounded by a/N, where N is the t o t a l number of g r i d points in intfi. PROOF: This is i n t u i t i v e l y obvious although there are subtle points in a rigourous j u s t i f i c a t i o n . For technical d e t a i l s see Appendix A. For ease of analysis we w i l l consider only relaxations where, in the equation corresponding to an i n t e r i o r grid point ( i , j ) , the values at each of the neighbouring points ( i - 1 , j ) , ( i , j-1) , ( i + 1 , j ) , ( i , j + 1) occur either on the l e f t side or on the right side of t h i s equation, but not both. Formally th i s says that only one of the members of each of the pairs ( 1 2 , r 2 ),...,(1 5, r 5 ) is not zero. A l l of the common relaxation schemes s a t i s f y this condition. Now we are ready to consider the rel a t i o n s h i p between equations (3-17) and (3-18). THEOREM 1: If 90 i s smooth and the induced 2-norm of L~ 1 is bounded independently of the number of grid points N, then there is a constant A independent of the grid such that (3-19) - u(d),p(6)\ < A|i//(0)|/VN PROOF: From (3-16) we find that (3-20) R\l/(6) = (R* - R1 )i/>(0) = u(6) (L + L1 )\jj(e) - R 1 ^ ( 0 ) . Thus from (3-18) 3 2 ( 3 - 2 1 ) = L " yR\p{9) = L " 1 (ix(d)L^(d) + u(d)L1\p(6) ) - L-1R1i// = fi(d)4/(6) + L " 1 ( M(0)L1 tV(0) - RI<//(0) ). This equation shows that v//(0) i s "almost an eigenvector of the relaxation" i f the term L~ 1 ( M ( 0) L11// ( 0) - R1iM0) ) i s of small norm. We w i l l show thi s by showing that L1\J/(0) and R1<H0) are of small norm. Consider the entries of R1i/>(0). This grid vector w i l l be zero at grid points not next to a boundary. Corresponding to a grid point exactly one of whose four nearest neighbours i s on the boundary, the vector w i l l have one of the following values: r 2 e x p ( - i 0 1 ) , r 3 e x p ( - i 0 2 ) , r * e x p ( i 0 1 ) , or r 5 e x p ( i 0 2 ) . For the case of a point with two nearest neighbours on the boundary, e.g., with points to the right and on the bottom, one obtains the entry r 3 e x p ( - i 0 2 ) + r a e x p ( i 0 1 ) (see Figure 4). At other 'corner' points and at points with three nearest __ 0 0 r*exp( i 0 1 ) __ 0 0 r"exp( t0 1 ) — 0 r 3 e x p ( - i 0 2 ) — r 3 e x p ( - c 0 2 ) + r " e x p ( 1 0 1 ) I Figure 4 - Entries in R1 Corresponding to Grid Points on the Boundary neighbours on the boundary, the entry w i l l s i m i l a r l y be the sum of two (or three) of the preceeding numbers. L1\^(0) w i l l be s i m i l a r . Let N(e), N(c) and N(w) be the number of points with one (edges), two (corners), and three (wedges) neighbouring boundary points, respectively. Let 33 a 0 = m a x ( | r 2 | , . . . , | r 5 | , | l 2 | , . . . , | l 5 | ) . Then in the grid vector M(0)L1I//( 0) - R1\//(0) there are N(e) entries of modulus bounded by a 0 (since |ju(0)| ^ 1), N(c) entries of modulus bounded by 2a°, and N(w) entries of modulus no larger than 3a°. Thus taking the 2-norm by summing the squares of the entries, (3-22) | M(0)L1<//(0) " R1i/>(0) | £ aV(N(e) + 4N(c) + 9N(w)). From Lemma 2 there is a number a such that for any grid 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- 1|a°VN 3v/a. If N i s the number of i n t e r i o r g r i d points then the norm of \p(8) , each of whose entries has modulus one, i s i/N. Since 3|L _ 1| a 0 /a i s bounded independently of N the result follows. QED This i s a bound on the r e l a t i v e norm of the difference between an actual relaxation and the LMA i d e a l i z a t i o n . We may also obtain simple lower bounds on this difference, i f for some b(e) > 0, N(e) ^ b(e)/N for N s u f f i c i e n t l y large. This occurs for c i r c u l a r and rectangular domains. Obviously for rectangular domains which have been rotated through 7r/4 radians there are no edge points so that the following simple argument won't go through. Let b° be the minimum of the non-zero c o e f f i c i e n t s | r21, ...,|r5|,112|,...,115|. Then the g r i d . vector 34 u(8)L]\p(6) - R1\/>(0) has at least N(e) entries of size at least |M(0)|b° > 0 (since each of the four nearest neighbours i s involved in the relaxation at a point). Thus since N(e) ^ b(e)/N (3-25) |M(0)L1<M0) - R1<M0)| > |y(0)|bVN(e) > \u(6) |bVb(e) VN. Now (3-26) min |L~ 1x|/|x| = (max |Lx|/|x|)" 1 =1/|L|. Thus from (3-20) and (3-26), (3-27) |*« - M ( 0 ) ^ ( 0 ) | ^ (1/VN) ( |bVb(e)/|L| ) \*(6) | .QED In interpreting i n e q u a l i t i e s (3-24) and (3-27) i t i s important to r e a l i z e that the estimate above i s NOT the difference between the norm of the result of relaxing the vector ii, and the norm of u(6)\}/. These norms may. be closer than t h i s estimate would suggest, as w i l l be seen in the tables. It is rather the norm of the vector difference between these two which i s considered. This difference may have a s i g n i f i c a n t component normal to \/>, which means that t h i s component "feeds" other Fourier components during relaxation, e s p e c i a l l y on small grids. It seems from some sample calculations that usually the projection of this vector difference - u(8)\p(8), is largest on Fourier components whose f i r s t frequency l i e s in the range (6 1-0(h),d 1+0(h)) and whose second frequency l i e s in the range (0 2-O(h),0 2+O(h)). 35 (0(h) means of bounded by ch where c i s a positive constant.) Except for simple cases, some of which w i l l be discussed l a t e r , the estimates (3-24) and (3-27) above are too cumbersome to be of much use in practice. However (3-24) can be used to j u s t i f y LMA asymptotically." This gives some psychological grounds for confidence in i t s p r a c t i c a l use on fine grids. On the other hand, (3-27) suggests that LMA i s not such a good estimate for very coarse grids. The examples w i l l i l l u s t r a t e how good LMA i s in some t y p i c a l s i t u a t i o n s . 3.3 Examples 3.3.1 Gauss-Seidel Relaxation of Laplace's Equation in the Unit Square with D i r i c h l e t Boundary Conditions. The gri d i s as follows. We may suppose that the boundaries coincide with grid l i n e s . The mesh in x and y w i l l be h = l/(n+l), and the values of x and y which occur at grid points w i l l be indexed from 0 through n+1. Points not on the grid boundary w i l l have both indices i , j between 1 and n. The f i n i t e difference system consists of n 2 equations, where equation ( i , j ) has the form (3-28) 4u(i,j) - u(i-1,j) - u ( i , j - l ) -u(i+1,j) -u(i,j+l) = 0, " One unexpected inference we may draw from (3-24) and (3-27) above i s that LMA works well for giving an asymptotic estimate of the a t t r i t i o n rate of the low frequency Fourier components. This_ i s contrary to Brandt's expectation (1977) that LMA predictions would only be good for high frequency components. This fact suggests the use of LMA for the understanding of the ove r a l l performance (as d i s t i n c t from smoothing rate) of relaxation schemes. This w i l l be discussed l a t e r . 36 where i , j run from 1 to n. If any indices in this equation take on values 0 or n+1 then the values w i l l be understood to be taken from the boundary data. Relaxation w i l l proceed from the bottom l e f t corner to the right (lexicographic ordering). The equation to be solved to obtain the next iterate x' from the current iterate x i s (3-29) 4 x ' ( i , j ) - x ' ( i - 1 , j ) - x ' ( i , j - l ) = x ( i+ 1 , j ) +x(i,j+ 1 ) . Then the generic equation transforming the errors from one it e r a t i o n to the next i s (3-30) 4 e ' ( i , j ) - e ' ( i - 1 , j ) - e ' ( i , j - l ) = e ( i+ 1 , j ) +e(i,j+l). We estimate the difference between the result of relaxing a Fourier component \p(6) and the LMA estimate u(&)\ls(6) by following through the steps of the theorem. F i r s t we estimate |L1i//(0)| and |R1i//(0)|. The nxn grid vector L1i//(0) i s 0 , except for 2(n-l) entries of modulus 1 along the top and right sides, and an entry of modulus between 0 and 2 in the bottom l e f t corner. R]\p(6) i s 0 except for 2(n-l) entries of modulus one along the top and right sides, and one of modulus between 0 and 2 in the top right corner. Then M ( 0 ) L 1 ^ ( 0 ) + R1^(0) w i l l have entries of modulus |M ( 0 ) | along the bottom and l e f t sides and of modulus 1 along the top and right sides. The sum of the moduli of the corner entries i s bounded by 6 ( 1 + | M ( 0 ) | ) . Thus 37 (3-31) (1+|M(0)I)/(2n-4) < \u(6)L^(8) + R1<//(c9) | < ( ]+u(6) | V(2n + 2) . The 2-norm of L i s bounded by 6 and that of L" 1 by 1/3. These estimates are obtained in Appendix B. The estimates of the difference between \//' and u(6)\j/(6) become (3-32) (l/6n)/(2n-4) < - »(d)t(6) \/\4>(6) | (1 + |M(0) I ) < (1/2n)/(2n+2). For an n=l00 by 100 gri d (10,000 i n t e r i o r points), (3-33) .0233 < -u(e)xP(6) \/\>l>(6) \ < .1421. This should be compared with Table I. The values in Table I were obtained as follows. For an (n+2)x(n+2) gri d , values were assigned to the nxn i n t e r i o r points as per the d e f i n i t i o n of \p(9) for various 0 1,0 2. One Gauss-Seidel relaxation sweep was performed, taking the values of boundary points to be 0. The norm of the relaxed vector was computed, as well as i t s projection on the o r i g i n a l Fourier component . These were both normalized, using |^(0)| = n. Table I i l l u s t r a t e s that as n gets larger, the LMA estimate better approximates the relaxed vector. 38 Table I - A t t r i t i o n Rates of Fourier Components of Gauss-Seidel Relaxation of Laplace's Equation in the Unit Square Grid (nxn) (6\d2)/ir Relative norm Projection |M(0)| \r\/\1>ie)\ (<//' ,tf(0) ) / \ $ ( 6 ) | 2 5x5 20x20 100x100 1 , 1 1/5,3/5 1/5,1 3/5,3/5 1 , 1 4,4 1/5,1/2 1/5,1 1/2,1/2 1 . 1 1/5,1/5 1/5,1/2 1/5,1 1/2,1/2 1 , 1 .6174 .3829 .1839 .3625 .3028 .7201 .4864 . 1 589 .4360 .3256 .7464 .4973 .1488 .4449 .3318 5859 3453 1 1 22 3303 2731 71 04 4768 1 369 4279 3187 7444 4953 1 442 4433 3304 .7529 .4232 . 1 460 .4004 .3333 .7529 .5000 . 1 460 .4472 .3333 .7529 .5000 . 1 460 .4472 .3333 3.3.2 Line Relaxation of .01(3 2U/3X 2) + (3 2u/3y 2) = 0 in the Unit Square. Table II was obtained l i k e the preceeding table using y - l i n e Gauss-Seidel relaxation on a 100x100 g r i d . In a sweep of y - l i n e relaxation, the values of the new ite r a t e along an entire g r i d l i n e of points whose second index j i s the same are determined simultaneously, using values of the new iterate at points whose second index i s j-1 and values of the current it e r a t e for points whose j index i s j+1. The LMA estimate i s less accurate here than for the point G-S relaxation on a comparable g r i d . This i s because the norm of the matrix L~ 1 i s much larger (although s t i l l bounded independently of h). 39 Table II - Line Relaxation (0',0 2 ) / 7 r Norm afte r relax- u(6) ation of Fourier mode .01,.01 .7767 .9094 .50,.01 .3907 .4301 .50,.05 .2121 .2187 .50,.50 .0051 .0049 1.0,1.0 .0025 .0025 3.3.3 Accuracy of Local Mode Analysis Estimates for Di s t r i b u t i v e Gauss-Seidel (DGS) Relaxation of the Cauchy-Riemann Equations We seek solutions u , and v to (3-34) 3u/3x = 3v/3y, -3u/3y = 3v/3x in the unit square. The domain on which a solution i s sought i s parcelled up into square c e l l s of side h. Values of the discrete variable u are sought at the centers of v e r t i c a l c e l l sides and values of v are sought at the centers of horizontal c e l l sides. This placing of discrete variables ensures second order accuracy of the discrete 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 f i r s t equation i s located at c e l l centers, and the second is located at c e l l corners. There i s no natural correspondence between variables and equations. To solve these equations a sweep i s made over a l l the f i r s t equations and then over the second equations. To relax 40 each equation each of the four variables that are involved in that equation are adjusted. The d e t a i l s of the d i s t r i b u t i v e relaxation may be found in Brandt(1979). In that paper he calculates a smoothing rate for residuals. The second column in Table III corresponds to the actual a t t r i t i o n rate of these residuals for the 40x40 g r i d . In the LMA i d e a l i z a t i o n , the amplitudes a 1 , a 2 of Fourier components r 1 = a 1 ( 0 ) e x p ( i ( i 0 1 + j 0 2 ) , and r 2 = a 2 ( 0 ) e x p ( i ( i 0 1 + j 0 2 ) before relaxation are related l i n e a r l y to the amplitudes a 1 , ( 0 ) and a 2'(0) after the sweep. Thus a matrix A(0) may be constructed specifying the linear linkage between the amplitudes of r 1 and r 2 before and after the sweep. In r e a l i t y a s i g n i f i c a n t component orthogonal to the o r i g i n a l vector a r i s e s . The norms of the relaxed vectors are compared below with the i d e a l i z a t i o n A(0) . The f i r s t column of each matrix records the r e l a t i v e norms of the relaxed residuals i f the o r i g i n a l residuals were r 1=\//( 0) , r 2 = 0 , the second column i f they were r 1 =0 , r 2=\i/( 0) . One sweep of DGS was performed taking the i n i t i a l and boundary values of u and v to be 0. After both phases of the relaxation sweep, the residuals were recomputed and the norms of the two residuals taken. These norms were normalized by the vector norms of the o r i g i n a l residuals, |r1|=N, | r 2 j = N-1. The t h i r d column records the actual transformation of errors, which i s too d i f f i c u l t to obtain via LMA. To obtain these values, vectors u and v were constructed, one of which was given values as per ^(0), and the other was set to 0. DGS relaxation was done, taking f1= f2=0, and the norms of the resulting vectors were taken and normalized. As may be seen from Table III the norms of the residuals 41 Table III - D i s t r i b u t i v e Relaxation of the Cauchy-Riemann Equations e\e2 7T/2 , TT/2 7T/40 , 7T 7T/40 , 7T/2 TT/40 , TT/40 LMA Estimate .333 0 0 .333 .440 0 0 .440 .020 0 0 .020 460 0 0 .460 Actual Residual Actual Error .997 0 0 .997 .110 0 0 . 1 07 721 0 0 .720 709 0 0 .716 727 0 0 .722 .997 0 0 .968 .548 .038 .038 .548 .444 .022 .022 .444 .119 .035 .025 .079 .465 .021 .020 .454 .973 .012 .012 .973 after relaxation are considerably d i f f e r e n t from what would be predicted on the basis of LMA. In fact the LMA prediction for the residuals seems to conform more to the norm of the errors after relaxation though the difference is s t i l l marked. In th i s situation then, LMA does not y i e l d very accurate estimates of the e f f e c t of relaxation on the various Fourier components. For this weakly e l l i p t i c system, the relaxation i s complicated and cannot be described in the terms used to prove Theorem I. 3.4 Comments and Miscellaneous Consequences 3.4.1 Alternate Fourier Components In the above discussion only exponential Fourier components have been used. There i s no reason a p r i o r i why we could not have embodied the various frequencies as sine Fourier components. With the usual indexing convention such a component corresponding to frequencies 01 and 62 would be 42 defined as (3-36) u ( 0 ) ( i , j ) = s i n ( 0 1 i ) s i n ( 0 2 j ) . In fact, sine Fourier components l i k e CJ above are the eigenvectors of the Jacobi relaxation of the model problem, whereas exponential components \jj are not eigenvectors for any scheme. However, one serious d i f f i c u l t y a r i s e s ; we cannot do LMA for sine Fourier components for lack of appropriate simple i d e n t i t i e s connecting values at adjacent points. As i t turns out the rates of a t t r i t i o n of sine Fourier components are quite d i f f e r e n t from exponential Fourier components for the same relaxation on the same problem. This casts some doubt on the interchangeability of the two types of Fourier components. 3.4.2 Alternate Norms We have used the 2-norm throughout, because th i s arises from an inner product which must be used to compute orthogonal projections. The reader may check that i f the 1-norm were used then the result of 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 useful bound i s obtained, since the contribution of boundary values to the norm of a gri d vector is not in proportion to their numbers. 3.4.3 Residuals Using LMA i t i s possible to j u s t i f y the i n t u i t i v e l y reasonable claim that the residuals 43 (3-38) r = b - AX = -Ae, r e f l e c t l o c a l errors in the current iterate X more than global errors, so that the size of residuals alone i s not a good estimate of the error in an approximate solution. We must also know whether the error 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 the use of LMA type computations of residuals, and then we use these computations to show that the low frequencies in the error contribute much less to the residuals than the high frequencies do. LEMMA 3: Suppose the difference equation (3-39) a 1 U ( i , j ) + a 2U(i-1,j) + a 3 U ( i , j ~ 1 ) + a"U(i-1,j) + a 5U(i,j+l) = h 2 b ( i , j ) i s to be s a t i s f i e d at every i n t e r i o r point ( i , j ) of a grid G covering a convex domain fl, and that zero order approximations to the boundary of fl and to D i r i c h l e t boundary conditions are made. Let \p(6) be as before. Then as h -> 0 the residuals associated with an error of the form \p{6) are asymptotically given by (3-40) r ( i , j ) = ( a 1 + a 2exp(-i0 1) + a 3exp(-i0 2) + a*exp(i0 1) + a 5exp(i0 2) ) \p(d){i,j). PROOF: Consider the system of equations 44 (3-41 ) ( I + A )X = 0. Suppose we relax t h i s system by solving an equation of the form (3-42) IX' = -AX. This i t e r a t i o n s a t i s f i e s a l l the conditions of Theorem I. If the error in the current iterate X i s \p(6) then the error in the next iterate X' i s r, the residual for AU = h 2b. Theorem I then says that (3-43) | r - M(0)*(0) I = 0 ( 1 / V N ) . But y(0) for this relaxation i s given by (3-40). QED In our description of residual computation we assume that the sum a 1+a 2+a 3+a 4+a 5 is 0. F u l l generality i s restored by adding a term a°U(i,j) i f the function u i t s e l f appears in the d i f f e r e n t i a l equation which underlies (3-39). Consider what happens to .a single Fourier component of the form e x p ( i ( 0 1 i + 0 2j) under the mapping e->r: (3-44) r ( i , j ) = e x p ( i ( 0 1 i + 0 2j))(a° + a 1 + a 2exp(- t0 1) + a 3exp(-i0 2) + a"exp(i0 1) + a 5 e x p ( i 0 2 ) ) If 0 1 and 0 2 are both small, then exp(±i0 1)and exp(±i0 2) are close to 1. At an i n t e r i o r point ( i , j ) therefore the modulus of r ( i , j ) i s close to |a°|. However i f 0 1 and 0 2 are both close to TT then exp(± i0 1 ) and exp(±i0 2) are a l l nearly -1, and 45 | D ( e ) ( i , j ) | i s close to |a° + 2a 1|. Lemma 3 insures that t h i s indicates the size of the global residuals. Thus unless - a 0 i s of the same order as 2a 1 ( i t i s t y p i c a l l y much smaller in pra c t i c e ) , the high frequency amplitudes are very much over-represented in the residual v i s - a - v i s the low frequency components. Since the lowest frequencies are on the order of the mesh size h, and exp(0(h)) i s 1 + 0(h), we see from (3-44) that the amplitude of a very low frequency residual i s only 0(h) of the amplitude of the corresponding low frequency error. On the other hand the amplitude of a very high frequency residual can be many times (for the Gauss-Seidel relaxation of the model problem a 0 + 2a 1 = 8) the amplitude of the high frequency error. 46 CHAPTER FOUR A HEURISTIC LOOK AT RELAXATIONS USING LOCAL MODE ANALYSIS In Chapter Three i t was shown that LMA could be j u s t i f i e d in the l i m i t of very fine grids. It was also shown that the v a l i d i t y of the function M(0) was not r e s t r i c t e d to values of 6\82 bounded away from 0. Here we compute the function u(6) for two common schemes: the SOR method and the ADI method. The form of uid) w i l l help us understand h e u r i s t i c a l l y what these schemes do. The algebraic theory of these schemes for special cases w i l l be adduced to confirm the understanding using LMA. LMA w i l l also be used to discuss why i t i s that most relaxations which meet the conditions outlined in Chapter Three do not e f f i c i e n t l y reduce smooth errors. 4 . 1 The SOR Method The SOR (Successive Over-Relaxation) method has been commonly used for constant c o e f f i c i e n t problems with some success. A parameter co sets the amount of over-relaxation performed. A c r u c i a l question i s how to pick co (0 < co < 2) in order to optimize the performance of the scheme over many it e r a t i o n s , i . e . , minimize the spectral radius of the amplification matrix. Here we study the behaviour of the SOR method applied to the Laplace equation in a domain J2 which has been d i s c r e t i z e d on a grid with mesh size h. The discrete Laplace equation at a point ( 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)). 4 7 Applying the SOR method, the equation of relaxation at a point ( i , j ) i s ( 4 - 2 ) X ' ( i , j ) = (l - c o ) X ( i , j ) + ( cu/4 ) ( X ' ( i - 1 , j ) + X ' ( i , j - 1 ) + X ( i+ 1 , j ) + X ( i , j + 1 ) ) . Subtracting ( 4 - 1 ) from ( 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 ) ) . Let e ( i , j ) = A ( 0 ) e x p ( i ( 0 1 i + 0 2 j ) ) and l e t e ' ( i , j ) = A' ( 0 ) e x p ( i ( 0 1 i + 0 2 j ) ) . Then ( 4 - 3 ) becomes ( 4 - 4 ) A * ( 0 ) = (i-u)A ( 0 ) + (w/4)( A' ( 0 ) e x p ( - i 0 1 ) + A'(0)exp(- i d 2 ) + A ( 0 ) e x p ( i 0 1 ) + A ( 0 ) e x p ( i 0 2 ) ) . The LMA estimate i s then given by ( 4 - 5 ) M ( 0 ) = (exp( it?1 )+exp( i 0 2 ) - 4 ( w - 1 )/w) / (4/co - e x p ( - t 0 1 ) - e x p ( - c 0 2 ) ) . The maximum of | M ( 0 ) | occurs at 0 1 = 82 = 0 and has the value 1 . The minimum occurs at ( 0 , T T) or (T T , 0 ) and has the value CJ - 1 . The p r o f i l e of the surface uid) resembles that for Gauss-Seidel relaxation, with a l o c a l maximum at ( 7 r , 7 r ) , see Figure 5 . 48 If CJ i s near 1 then the values of u(6) are well ca = 1 (Gauss-Seidel) CJ = 1.5 Figure 5 - Contour Lines of |M(0)| for SOR Relaxation d i s t r i b u t e d , but as CJ increases to 2 the values of u(6) cluster together. For model problems i t i s known (Young, 1971) that the eigenvalues of the amplification matrix behave similar i l y . Can u(6) be used to estimate the rate of convergence of the SOR relaxation? Unlike de Vries (1982), I found that M(7rh ,7rh) was not a good estimate of K the spectral radius of the amplification matrix. There i s no other clear choice of Fourier component which should be analogous to the smoothest eigenvector of the relaxation (since M(0,0) = 1 always, i t cannot be used). T y p i c a l l y 1 - K was two or three times bigger than 1 - Ai (7 rh ,7 rh ) for problems in the unit square. Nevertheless, in certain cases i t i s possible to use LMA to estimate the best value of the parameter CJ as described below. We consider Laplace's equation in the following domains: 0 1 i s the unit square; fi2 i s the c i r c l e (x-1/2) 2 + (y-1/2) 2 <l/4; fi3 i s the curved region {(x,y)| y>0, 4/5-(2x-l) 2< y <1-(2x-1) 2}. F i n i t e difference systems were set up for these problems for several values of h, as outlined in Table IV. It was not clear 49 which frequency ought to be analogous to the smoothest eigenvector for the domains fi2 and fl3. Since the number of grid points i s readily computed, the value of M was computed for ft = (7r//N, 7r/v/N) . This is nearly equivalent to ( 7 r h , 7 r h ) in the unit square. The value cjr of CJ which minimizes \u\ i s recorded in Table IV. For comparison purposes the value C J * of CJ which minimizes the actual radius of the amplification matrix i s included. For the unit square C J * i s known from a complete analysis of this model problem (see Young(1971)), and i s given by (4-6) CJ* = 2/( l+sinUh) ) . For the other problems, K ' an approximation to K , was determined using the successive iterates X1,X2,... . It i s defined by (4-7) K ' = |X 5 2 - X 5 1| / |X51 - X 5 0|. An C J * was computed which minimized K ' as a function of C J . We see that C J ' i s a good estimate of C J * p a r t i c u l a r l y , for large N. For fl3 the r a t i o of the number of gri d points next to the boundary to the t o t a l number of i n t e r i o r g r i d points i s rather larger than for the other domains with comparable N. The v a l i d i t y of LMA depends on t h i s r a t i o being small, and so i t i s not surprising that estimates for J23 were not as good as those for fl1 and fl2. Can LMA be used to choose parameters for more general 50 Table IV - Best values of CJ h N C J * . 2 1 6 1 . 2 3 6 1 . 2 6 0 . 1 81 1 . 5 2 4 1 . 5 2 8 . 0 5 361 1 . 7 2 9 1 . 7 3 0 . 0 2 2 4 0 1 1 . 8 8 2 1 . 8 8 2 n 2 . 0 2 1 8 4 7 1 . 8 6 1 . 8 9 . 0 2 4 1 9 1 . 7 3 1 . 5 8 . 0 1 1 7 8 7 1 . 8 6 1 . 7 9 . 0 0 5 7 3 7 1 ' 1 . 9 3 1 . 9 3 problems? I considered extension to two classes of equations; asymmetric constant c o e f f i c i e n t equations and equations with one variable c o e f f i c i e n t . I encountered two d i f f i c u l t i e s with asymmetric operators of the form - 9 2 / 9 x 2 - 9 2 / 9 y 2 - a9/9 x . The function UL{6) computed by LMA for such an operator may have a maximum at (T T , 0 ) rather than at ( 0 , 0 ) for some values of C J . There i s then even less p l a u s i b l i t y in computing u, for a low frequency component. Secondly, the eigenfunctions of the asymmetric operator are skewed exponentially (u(x,y) = exp(-ax/2) s i n ( n 7 r x ) sin ( m 7 r y ) in the unit square, for integers m,n). For h of order 1/a or larger, the discrete eigenfunctions are quite d i s s i m i l a r to the Fourier components, and thus we cannot expect LMA to be v a l i d at a l l . Several methods to estimate optimal parameter values using LMA were t r i e d and none approached the experimentally determined optimal value. Extension to problems with one variable c o e f f i c i e n t was more successful. Numerical experiments were done using the following problem: 51 (4-8) (1/10 +x2 + y 2)9 2u/9x 2 + 9 2u/9y 2 = 0, with D i r i c h l e t boundary conditions on [0,1]X[0,1]. A mesh size of .02 was used. The determination of co' was insensitive to the c o e f f i c i e n t of 9 2u/9x 2 over the range .1 to 2.; CJ' = 1.88. The value of co* was determined experimentally, as outlined previously, to be co* = 1.90. Similar comparisons where a term -u/lOh 2 was included in the discrete operator f a i l e d to give such good agreement between co1 and co* . Apparently LMA has some limited usefulness with respect to choice of parameter for solution of Laplace type equations by SOR relaxation. This i s p a r t i c u l a r l y true for domains with a small r a t i o of boundary length to area, approximated by fine grids. However, there is no obvious way to use LMA to estimate the actual rate of convergence with the optimal parameter. 4.2 Alternating-Direction Implicit Iteration For ADI i t e r a t i o n on an equation AX = h 2b, the matrix A is s p l i t into matrices H 1 and V 1. For any grid vector X, and any gr i d point ( i , j) , ( H 1 X ) ( i , j ) involves only the values of X at the points ( i , j ) , ( i - 1 , j ) , and (i+1,j), and V'X(i,j) involves only the values of X at the points ( i , j) , ( i , j - 1 ) , and ( i , j + l ) . Thus H1 and V 1 represent the f i n i t e - d i f f e r e n c i n g in the x and y d i r e c t i o n s , respectively. A parameter p, which may vary from one i t e r a t i o n to the next, determines the form of the relaxation. Optimal choice of values for p is important" in ensuring quick convergence. One i t e r a t i o n of ADI consists of solving 52 (4-9) (H 1 + p h 2 I )X* = ( p h 2 I - V 1 )X + h 2 b , and then s o l v i n g (4-10) (V 1 + p h 2 I )X' = ( p h 2 I - H 1 )X* + h 2 b , to o b t a i n the next 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 , u s i n g X* as an i n t e r m e d i a t e . I f the two m a t r i c e s H 1 and V 1 commute, then the a l g e b r a i c t h e o r y guarantees r a p i d convergence of the i t e r a t e s (see Young 1971). The m a t r i c e s commute o n l y for the s p e c i a l case of s e p a r a b l e e q u a t i o n s of the form (4-11) ( 3 / 3 x ) ( f ( x ) 3 u / 3 x ) + ( 3 / 3 y ) ( g ( y ) 3 u / 3 y ) + cu = h ( x , y ) on r e c t a n g u l a r domains . I f these c o n d i t i o n s are not met, ( for example, i f l a r g e f i r s t - o r d e r d e r i v a t i v e terms are p r e s e n t ) then convergence i s f r e q u e n t l y slow i n p r a c t i c e . An h e u r i s t i c e x p l a n a t i o n of t h i s phenomenon u s i n g LMA i s now g i v e n . Theorem 1 of Chapter 2 p r o v i d e s some j u s t i f i c a t i o n for t h i s n o n - r i g o r o u s a p p r o a c h . The l o c a l mode a n a l y s i s f or the ADI method a p p l i e d to L a p l a c e ' s e q u a t i o n runs as f o l l o w s . The r e l a x a t i o n a t g r i d p o i n t ( i , j ) c o n s i s t s of s o l v i n g (4-12) ((2 + p h 2 ) X * ( i , j ) - X * ( i - 1 , j ) - X * ( i + 1 , j ) ) = ( ( p h 2 - 2 ) X ( i , j ) + X ( i , j - 1 ) + X ( i , j + 1 ) ) + h 2 b ( i , j ) , 53 and then s o l v i n g (4-13) ((2 + p h 2 ) X ' ( i , j ) - X ' ( i , j - 1 ) - X ' ( i , j + D) = ( ( p h 2 - 2 ) X * ( i r j ) + x * ( i - l , j ) + X * ( i + l , j ) ) + h 2 b ( i , j ) . L e t e* be the e r r o r i n X*. Then the e r r o r t r a n s f o r m a t i o n c o r r e s p o n d i n g t o (4-13) i s (4-14) ((2 + p h 2 ) e * ( i , j ) - e * ( i , j - 1 ) - e * ( i , j + l ) ) = ( ( p h 2 - 2 ) e ( i , j ) + e ( i , j - D + e ( i , j + l ) ). L e t e ( i , j ) = A ( 0 ) e x p ( i ( 0 1 i + 0 2 j ) ) , and l e t e * ( i , j ) = A * ( 0 ) e x p ( i ( e ' i + & 2 j ) ) . Then (4-15) becomes (4-15) A * ( 0 ) ( 2 + ph 2 - 2 c o s 0 1 ) = A ( 0 ) ( -2 + p h 2 + 2 c o s 0 2 ) . U s i n g s i m i l a r arguments we a r r i v e a l s o a t (4-16) A ' ( 0 ) ( 2 + p h 2 - 2cos0 2) = A * ( 0 ) ( -2 + p h 2 + 2COS0 1). Thus (4-17) n(6) V - S V - t V + S V + t where v = p h 2 , s = 2 - 2 c o s 0 1 , and t = 2 - 2 c o s 0 2 . N o t i c e t h a t i t i s p o s s i b l e t o make ix(d) e q u a l t o 0 f o r any 8 by s u i t a b l e c h o i c e of p. For s m a l l p the components w i t h s m a l l 0 1 or s m a l l 62 v a l u e s a r e b e i n g e l i m i n a t e d s i n c e 1 - c o s ( 0 ( h ) ) = 0 ( h 2 ) . As p i s i n c r e a s e d , the z e r o of M(0) 54 o c c u r s a t i n c r e a s i n g l y h i g h e r f r e q u e n c i e s , u n t i l f o r p = 4/h 2 t h e z e r o i s ( 7 T , 7 r ) . T h i s i s t h e same ra n g e f o r p t h a t was a r r i v e d a t by more d i f f i c u l t c a l c u l a t i o n s i n Young ( 1 9 7 1 ) . Now c o n s i d e r t h e PDE (4-18) 9 2 u / 9 x 2 + 9 2 u / 9 y 2 - a9u/9x - b9u/9y = 0. As m e n t i o n e d a b o v e , ADI methods a r e n o t e x p e c t e d t o work a s w e l l f o r t h i s e q u a t i o n . The f i n i t e d i f f e r e n c e f o r m of t h e e q u a t i o n a t g r i d p o i n t ( i , j ) i s (4-19) - 4 U ( i , j ) + d + a h ) U ( i - 1 , j ) + ( 1 + b h ) U ( i , j - 1 ) + ( l - a h ) U ( i + 1 , j ) + ( 1 - b h ) U ( i , j + 1 ) = 0. Then LMA y i e l d s (4-20) M(0) = v ~ s - 2 i a h s i n 6 1 ^ v - t - 2 i b h s i n 6 2 v + s + 2 i a h s i n 6 1 v + t + 2 i b h s i n 6 2 N o t i c e t h a t i t i s n o t p o s s i b l e t o c h o o s e p i n o r d e r t o make u(6) e q u a l t o 0 f o r any g i v e n 0. W h i l e f o r s m a l l p t h e minimum of |ju(0)| s t i l l o c c u r s a t low f r e q u e n c i e s , i t s v a l u e i s o f o r d e r 1 i f ah and bh a r e 0 ( 1 ) . F o r l a r g e r t h e minimum o c c u r s a t h i g h f r e q u e n c i e s and a l t h o u g h i t s v a l u e i s n o t 0, i t w i l l be s m a l l ( 0 ( h ) ) p r o v i d e d t h a t ah and bh a r e n o t 0 ( h " 1 ) . Thus t h e r e i s no c h o i c e of p w h i c h w i l l e f f e c t i v e l y r e d u c e t h e smooth (low f r e q u e n c y ) e r r o r components. Thus a l t h o u g h t h e ADI method w i l l s t i l l e f f e c t i v e l y smooth t h e e r r o r i t s o v e r a l l p e r f o r m a n c e w i l l be s i g n i f i c a n t l y d e g r a d e d from i t s 55 performance on the model problem. 4.3 Smoothing Property of Relaxations Suppose we have a difference equation at every i n t e r i o r point of a g r i d covering a domain fl; (4-21) a 1 U ( i , j ) + a 2U(i-1,j) + a 3 U ( i , j - l ) + a 4U(i+1,j) + a5U(i,j+D = h 2 b ( i , j ) , where a 1,...,a 5 are constant c o e f f i c i e n t s . Suppose there i s no term of degree zero in u in the underlying PDE so that (4-22) a 1 +a2 + a 3 + a• + a 5 = 0. Suppose we have a s p l i t t i n g A = L - R which yie l d s a relaxation scheme for (4-22), (4-23) l 1 X ' ( i , j ) + l 2 X ' ( i - 1 , j ) + l 3 X ' ( i , j - D + 1 " X ' ( i + 1 , j ) + 1 5 X ' ( i , j + i ) = r 1 X ( i , j ) + r 2 X ( i - 1 , j ) + r 3 X ( i , j - l ) + r 4X(i+1,j) + r 5 X ( i , j + l ) + h 2 b ( i , j ) . Then since 1* - r* = a*, * = 1,...,5, we must have (4-24) 1 1 + 1 2 + 1 3 + 1 « + 1 5 = r 1+r 2+r 3+r"+r 5. The reader w i l l r e c a l l from Chapter Three that the function u(6) for the relaxation scheme (4-24) i s 56 (4-25) u(6) = ( r 1 + r 2 e x p ( - i 0 1 ) + r 3 e x p ( - i 0 2 ) + r ' e x p d f l 1 ) + r 5 e x p ( t 0 2 ) ) / ( l 1 + l 2 e x p ( - i 0 1 ) + l 3 e x p ( - t 0 2 ) + l"exp( i0 1 ) + l 5 e x p ( i d 2 ) ). Thus, i f (4-26) l 1 + l 2 + l 3 + 1" + l 5 * 0, 1 (4-27) u(e) - > 1 , as e\e2 -> o. The l i m i t (4-28) says that the function u(6) i s near 1 i f 6 \ 62 are near 0. Thus for the re laxat ion of a constant c o e f f i c i e n t equation by a scheme which s a t i s f i e s the condit ions out l ined in Chapter Three, the low frequency Fourier components in the error are not e f f i c i e n t l y reduced. The best we can hope for in general is e f f i c i e n t smoothing, but not e f f i c i e n t reduction of the error o v e r a l l . 1 Only the ADI re laxat ion of the model problem with a small value of the parameter r comes close to v i o l a t i n g th i s cond i t i on . The reader w i l l r e c a l l that low frequency errors could be e f f i c i e n t l y relaxed in th i s case. 57 CHAPTER FIVE DISCUSSION 5.1 A p p l i c a b i l i t y of Local Mode Analysis The results in section 2.2 j u s t i f y LMA under the following conditions. The relaxation for which LMA i s done must be a linear fixed point i t e r a t i o n which i s 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 into A = L - R. Alternately each relaxation step could be a sequence of several procedures each of which 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 regular, they must be homogeneous, i . e . , each point must be in the same position r e l a t i v e to the set of previously relaxed points. Thus red-black relaxations are excluded. A further condition, which we imposed for ease of analysis, was that in the relaxation equation at a given point, no other neighbouring points occur than those which occur in the o r i g i n a l difference equation at that point. We considered only second order difference schemes based on f i v e point formulae, which excluded relaxations for which the relaxation at point ( i , j ) required values at points other than the four nearest neighbours. This condition may not, in fact, be neccessary for the results but no relaxation scheme i s known to us for which i t i s not true. The condition that the relaxation 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 analysis was meaningful. We do not think that Fourier analysis is relevant to the understanding of schemes such as the Conjugate Gradient Method. The homogeneity condition guaranteed that the 58 cross-coupling of Fourier components undergoing relaxation was asymptotically negligible; i . e . , that each Fourier component behaved l i k e an eigenvector. A l o c a l Fourier analysis is possible for inhomogeneous schemes such as red-black relaxation; then however, because cross-coupling of the Fourier components is important, they must be arranged in invariant subspaces rather than being treated l i k e eigenvectors. Under red-black ordered Gauss-Seidel relaxation, Fourier components with frequencies ( 0 1 , 0 2 ) , ( 0 1 , 7 T - 0 2 ) , ( T T - 0 1 , 0 2 ) , and ( 7 T-0 1 , 7 T-0 2) feed into each other s i g n i f i c a n t l y for a l l mesh sizes, and must be treated as a four-dimensional invariant subspace (see Trottenberg and Stueben, 1982). Although the theorem in Chapter Three was established only for relaxations of scalar problems with D i r i c h l e t boundary conditions which met the requirements described above, some extensions may be contemplated. Scalar problems with Neumann or mixed boundary conditions or higher order approximations to D i r i c h l e t boundary conditions may be treated using the method in Chapter Three, however, the algebra may be much more complicated. It seems straightforward to extend the result to strongly e l l i p t i c systems. Extension to weakly e l l i p t i c systems may be possible, but 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 relaxation schemes. Of course a theorem l i k e that in Chapter Three i s not meaningful for variable c o e f f i c i e n t problems or non-linear problems. 5.2 The Potential of LMA Outside the MG Context The theorem in Chapter Three j u s t i f i e d our explorations 59 in Chapter Four. We saw that LMA could be used to explain the behaviour of the SOR and ADI methods, and also the fact that linear relaxations reduce smooth errors poorly. We speculate that LMA has the potential to be more widely used in order to gain further h e u r i s t i c understanding of such methods. Frequently 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 relaxation scheme for the solution of a problem and wishes to know how appropriate i t i s or how to pick relevant parameters. A rigourous algebraic answer would involve finding the eigenvalues and the eigenspaces of the amplification matrix of the relaxation. This i s generally more work than i s needed to solve the o r i g i n a l problem. On the other hand LMA i s a simple a n a l y t i c a l tool and i t can y i e l d much the same information. 5.3 Possible Developments of the Theorem The theorem in Chapter Three establishes that LMA is a meanigful tool and imparts confidence in i t s use for MG problems. However, we have only shown that LMA y i e l d s a good description of the r e s u l t s of relaxation where the error i s a multiple of a single Fourier component. As w i l l be discussed below i t would be very useful to have a result of this kind for a r b i t r a r y errors. Suppose we have an a r b i t r a r y error e= x-u and decompose i t as (5-1 ) e = £A(0H(0) . e We ask how closely does e' the error after relaxation resemble e* the LMA estimate, which i s given by 60 (5-2) e* = £ ^ ( 0 ^ ( 0 ) ^ ( 0 ) . e An argument along the l i n e s of Chapter Three might be made that the difference e' - e* i s asymptotically small. We w i l l sketch t h i s here. We w i l l consider only rectangular regions so that the Fourier decomposition (5-1) i s well defined and the c o e f f i c i e n t s A(0) form the f i n i t e Fourier transform of e , where 01 and 0 2 take the values 27rk/N for k = 0,1,...,N-1. Note that the low frequency region i s (0 , 7r/2 )U( Ztt/2 , 27r) with these choices of frequencies. This i s equivalent to the previous d e f i n i t i o n . Now matrices R* and L* may be constructed independently of 0 so that (5-3) R*t//(0) = r*(e)xp(6) and L*«//(0) = l*(0)<//(0), for 0 1 , 0 2 = 27rk/N. Then since r*(0)/l*(0) = M(0), formula (3-17) reads (5-4) e' = e* + L" 1(L1e* - R1e), where again L1 = L* - L and R1 = R* - R. Now however we have no guarantee that L i e * and R1e are asymptotically small r e l a t i v e to e. In fact i f e is non zero only on points next to the boundary, L i e * and R1e can be large for any N. Thus the estimate obtained in Chapter Three may not be extended to arbi t r a r y errors. Does th i s fact invalidate convergence estimates of MG algorithms that rely on LMA for estimating the smoothing factor for any error? I think not. The smoothing rate i s taken as M=max|M(0)| for one of 0 \ 0 2 in the range (7r/2,3 7r/2) . Most 61 of the high frequency Fourier components are reduced by a factor much smaller than u. I believe that i t is s t i l l true that the non-smooth part er 1 of an a r b i t r a r y error e i s s t i l l reduced by a factor which is bounded by u + 0 ( i / V N ) . If t h i s could be proved then i t would be possible to show convergence of the MG method in much more generality than previously. Trottenberg and Stuben (1981) have already shown that direct solution of the coarse grid problem (the coarse grid correction step) s i g n i f i c a n t l y reduces the smooth component es = e - er of the error in the fine grid problem. Thus a proof that the relaxation on the fine gr i d strongly reduces the non-smooth error i s a l l that i s needed to show convergence of the whole MG algorithm. So far proofs of the convergence of the MG method have been d i f f i c u l t and r e s t r i c t e d to a few simple cases. They have also been highly a r t i f i c i a l . A proof via the LMA estimate which is used in practice would go a long way toward bringing together theory and ap p l i c a t i o n . 5.4 Summary There seems to be potential for development of LMA beyond the p r a c t i c a l role for which Brandt introduced i t . As we have seen, under suitable conditions, l o c a l Fourier analysis can y i e l d quantitative information about the performance of relaxation schemes outside the MG context. It can also be made rigourous and thus is possibly part of a sound mathematical basis for the MG method. 1 This i s the orthogonal projection of e onto the space spanned by {4/(6) \ 61 , d2 = 2nk/U, one of 6\62 e (n/2 , 3TT/2 )} 62 BIBLIOGRAPHY Brandt, A. "Multi-Level Adaptive Solutions to Boundary-Value Problems", Mathematics of Computation, 31 (1977), 333-390. Brandt, A. "Guide to Multgrid Development", in Multiqrid, Methods, W. Hackbusch and U. Trottenberg, eds., Koeln-Porz, Nov 1981. LNM 960, Springer Verlag (1982), 220-312. Hackbusch, W. "Multi-Grid Convergence for a Singular Perturbation Problem", Proceedings, Copper Mountain, Apr 1983. (to appear). Hemker, P. W. "Fourier Analysis of Grid Functions, Prolongations and Restrictions", Report NW 93/80, Mathematical Centre, Amsterdam (1980). Hemker, P. W. "The Incomplete LU-Decomposition as a Relaxation Method in Multi-Grid Algorithms", in Boundary and Interior  Layers - Computational and Asymptotic Methods,. J. M i l l e r ed., Boole Press, Dublin (1980), 402-406. Hemker, P. W. "An Accurate Method Without Direc t i o n a l Bias for the Numerical Solution of a 2-D E l l i p t i c Singular Perturbation Problem", Report NW 117/81, Mathematical Centre, Amsterdam (1981). 63 Mol, W.J. "Smoothing and Coarse Grid Approximation Properties of M u l t i g r i d Methods", Report NW 110/81, Mathematical Centre, Amsterdam (1981). Trottenberg, U. and K. Stueben. "Multigrid Methods", in Mult i g r i d Methods, W. Hackbusch and U. Trottenberg, eds., Koeln-Porz, Nov 1981. LNM 960, Springer Verlag (1982), 1-176. de Vries, H.B. "The Two-Level Algorithm for the Solution of the I n i t i a l Value Problem for P a r t i a l D i f f e r e n t i a l Equations", Report NW 136/82, Mathematical Centre, Amsterdam (1982). Young, D. M. Iterative Solution of Large Linear Systems. Academic Press, New York, 1971. 64 Appendix A Proof of Lemma 3, Chapter Two Let 7 ( t ) : [ 0,L]-> R 2 be a continuous piecewise d i f f e r e n t i a b l e curve parametrized by arc l e n g t h . For 0 < t < L, l e t H(t,h) be (A- 1 ) H(t,h) = {(x,y)| | ( x , y ) - 7 ( t * ) | < h f o r some t * , 0<t*<t} Define A l ( t ) as the area of H ( t , h ) , denoted A(H(t,h)) Lemma A 1 : If 7 i s a s t r a i g h t l i n e then dAl/dt = 2 h . Proof: Obvious. Lemma A 2 : A(H(L,h) ) = A l ( L ) < T r h 2 + 2 h L Proof: Let D(a,h) be a d i s k of r a d i u s h centered at a po i n t a in R 2. Then c l e a r l y A l ( t ) = A( U D ( 7 ( r ) , h ) , 0<r<t). C l e a r l y A l i s a continuous non-decreasing f u n c t i o n of t . I t i s d i f f e r e n t i a b l e whenever 7 i s . Let t be a p o i n t where 7 i s d i f f e r e n t i a b l e ; then (see F i g u r e ) (A- 2 ) A l ' ( t ) = l i m (A(H(t+e,h)) - A ( H ( t , h ) ) ) / e 6-HD = l i m ( A ( D ( 7 ( t ) + e 7 ' ( t ) , h ) U H ( t , h ) ) - A(H(t,h))) / e e-K) = l i m A ( D ( 7 ( t ) + e 7 ' ( t ) , h ) \ H ( t , h ) ) / e < l i m A ( D ( 7 ( t ) + e 7 ' ( t ) , h ) \ D ( 7 ( t ) , h ) ) / e. e-H3 65 HI D( 7(t),h) ///D( T(t)+e T'(t),h) 7(t) 7 ( t ) + e 7 ' ( t ) Figure This l a s t l i m i t i s independent of 7 (since | 7 ' ( t ) | = 1 was spe c i f i e d ) . If 7 is a straight l i n e the la s t inequality i s an equality and for thi s case A l ' ( t ) = 2h by Lemma A1. Thus (A-3) A l ' ( t ) < 2h. Since Al(0) = 7rh2, the lemma follows. QED Lemma A3: (Lemma 2 of Chapter Three): Let fl be a domain with piecewise smooth boundary 7. There i s a constant a such that i f G i s any s u f f i c i e n t l y fine g r i d , and N is the number of grid points of G which l i e i n t e r i o r to fl, then the number M of such points which l i e next to a grid point which i s outside fl is bounded by (A-4) M < a/N Proof: (I am indebted to David Boyd of the Math. Dept. UBC for his suggestion of using gri d squares) Let A be the area of fl, L be the length of 7, and h be the mesh size of G. Let n be the number of grid squares (each of area h 2) which l i e e n t i r e l y in fl. Let m be the number of squares which intersect 7, or whose boundary 66 intersects 7. Then (A-5) A < (m+n)h2. The m squares which meet 7 must l i e e n t i r e l y in the set H(L,2/2"h) defined above. Now by Lemma A2 (A-6) mh2 < A(H(L,2v/2h) < 8vrh2 + 4v/"2hL. Choose h < /2L/(8TT). Then (A-7) mh2 < 5/2Lh. Now choose h < A/(10/2L). Clearly A(H(L,2/2h)) < A/2, thus by (A-7) and (A-9) (A-8) nh 2 > A/2. Thus (A-9) m < 5/2L/h = (10L//A)(/A//2h) < (10L//A)/n. Now each of the n squares i n t e r i o r to S2 contribute at least one point to N (say the lower l e f t corner). Thus n < N. Each of the M points which l i e next to a point outside £2, i s a corner of one of the m squares. Clearly there are at most three contributions 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" 1 for the Gauss-Seidel Relaxation of the Model Problem 1. Form of L and L" 1 For an nxn grid with n 2 i n t e r i o r points the n 2 by n 2 matrix L has the form (B-1 ) A 0 0 -I A 0 0 -I A 0 0 . . . 0 -I A where A is an nxn matrix of the form (B-2) 4 0 0 0 . . . 0 -1 4 0 0 • 0 -1 4' 0 • . 0 0 . . . 0 - 1 4 The reader may v e r i f y that A - 1 has the form 68 1/4 0 0 0 ... 0 1/16 1/4 0 0 1/64 1/16 1/4 0 (B-3) 0 • • . . . 1/16 1/4 where in the lower triangle the ( i , j) entry i s 4 to the power -1 - i + j . L" 1 has the form A" 1 0 0 0 ... 0 A' 2 A"1 0 0 A" 3 A" 2 A"1 0 (B-4) . . 0 . . . A" 2 A" 1 where the block in position ( i , j ) is A to the power -1-i+j. If we compare the ( i , j ) entry in AA~k with the ( i , j ) entry in A we obtain (B-5) 4A k ( i , j ) = A _ k ( i - 1 , j ) + A " k + 1 ( i , j ) . C learly the elements along the main diagonal of A k are (B-6) A _ k ( i , i ) = 4 • 4 k 69 Using (B-6) and (B-5) the reader may v e r i f y that (k+i-j-1 )! k+i-j - -/4 , i f i ^ j , and (k-1)! ( i - j ) l ' 0, i f i<j. Note that the entries in the rows are equal to the entries in the columns under the correspondence (B-8) A " k ( i , j ) = A'k (n-j + 1 ,n-i + 1 ) . Inspection of the form of (B-4) shows that then (B-9) L " 1 ( i , j ) = L" 1(n 2-j+1,n 2-i+1). 2. Norms The square of the 2-norm induced on L by the 2-norm on gr i d vectors x i s (B-10) |L| 2 = max{(Lx,Lx)| | X| 2=1}. Thus |L|2= /K where K i s the largest eigenvalue of L TL. Let z be an eigenvector of L T L with eigenvalue K. Then (B-11) |K| = |L^Lz| /|z| < |LT| |L| = | L | J L ! , 1 1 1 1 1 where |L| denotes the induced 1-norm which i s the maximum of the sums of the entries in the columns of L, and |L|w denotes (B-7) A " k ( i , j ) = 70 the induced supremum norm which i s the maximum of the sums of the entries in the rows of L . From (B-1) and (B-2) both are 6. Thus K < 36 and and we may write (B-14) | L - ' | < I L - 1 ! . 2 1 The reader may v e r i f y by inspection of (B-4) and (B-7) that the sum of the entries in the f i r s t column of L ~ 1 i s larger than the sum of the entries in any other column. Thus (B-12) | L | 2 < 6 . By a similar argument (B-13) | L - ' | 2 2 < I L - ' I J L - ^ . Because of the correspondence (B-9), (B-15) |L-'| 2 " L, ,m+l l> n (m-j):j! j m=l 4 j=0 C O I 4 m=0 4 m+1 0 0 I — m=0 2 m+2 1/2. 

Cite

Citation Scheme:

        

Citations by CSL (citeproc-js)

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>
                        
                    
IIIF logo 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

Comment

Related Items