ABSTRACT BOUNDS on THE DISCRETIZATION ERROR FOR THE NUMERICAL SOLUTION OF THE DIRICHLET PROBLEM FOR LAPLACE'S EQUATION BY Robert L. Meyer Let u(x,y) denote the continuous solution to a Dirichlet problem for Laplace's equation defined on a region R in the (x,y)-plane which is bounded by a piecewise analytic curve C. Suppose the data prescribed on C is continuous. Let Gh be an approximating mesh of width h which covers R, and let uh(P) denote the discrete solution of a finite difference analogue of the continuous problem. The function e(P) = u(P) - uh(P), defined for points P E Gh’ is called the discretization error for the Dirichlet problem for Laplace's equation. Bounds on the absolute value of the discretization error, ‘e(P)‘, are Obtained in this paper. The absolute value of the discrete Laplacian of the contin- uous function, denoted by ‘A:u(P)‘, forms a part of the bound for \e(P)|. In Chapter 2, bounds for ‘A:u(P)‘ are obtained from integral and Fourier series representations of the function u(x,y). These bounds are an improvement of any previous bounds in that they are in terms of known data and do not involve derivatives of the unknown function u(x,y). In Chapter 3, the bounds for ‘A;u(P)‘ are used to obtain bounds on \e(P)‘ for a finite difference analogue of a Dirichlet problem defined on a polygon. In Chapter 4, the Schwarz Reflection Principle is used Robert L. Meyer to extend the usefulness of the bounds for \A:u(P)‘ and to obtain some bounds on ‘€(P)| for a Dirichlet problem defined on a region R bounded by a more general piecewise analytic curve. The Appendix contains an example of some computed bounds for ‘e(P)‘ for a Dirichlet problem defined on a square. BOUNDS ON THE DISCRETIZATION ERROR FOR THE NUMERICAL SOLUTION OF THE DIRICHLET PROBLEM FOR LAFLACE'S EQUATION BY \ Robert L: Meyer A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Mathematics 1969 ACKNOWLEDGMENTS I wish to express my sincere thanks to my major professor, Dr. Norman L. Hills, for the role he has played in the develop- ment of this thesis. I became interested in the problem discussed in this thesis largely because of his excellent presentation of the basic related material in courses at Michigan State University. I am most grateful to him for the encouragement and assistance he has given me. My sincere thanks also to Professors J.S. Frame, W.T. Sledd, and R.H. Wasserman who gave me some helpful suggestions during the later stages of the preparation of this thesis and served as members of my guidance committee. I am grateful also to the other members of the Mathematics Department who expressed an interest in this problem and encouraged me during the time in which the research was conducted. ii Chapter I. Introduction .................... . ...... . .......... II. Bounds for the Discrete Laplacian of a Harmonic Function 00....-OOOOIOOOOOOCOOOOIOOOOOO 00000000 III. The Discretization Error for a Dirichlet Problem Defined on a Polygon ........... ....... IV. The Discretization Error for a Dirichlet Problem Defined on a Region Bounded by a Piecewise Analytic Curve ... ............ . ...... APPENDIX ....... . .. ........... . .. . ............... . BIBLIOGRAPHY ............ . ...... . .. . ................ TABLE OF CONTENTS iii 38 55 65 89 CHAPTER I INTRODUCTION Let R be a bounded simply connected region in the (x,y)- plane whose boundary C is a contour. Let f(x,y) be a function which is defined and continuous on C. The Dirichlet problem for Laplace's equation is the problem of determining a function u(x,y) defined and continuous on G = R+C, harmonic in R, and identical with f(x,y) on C. Denote this problem by 2 2 AU(x,y) = L; + 5421' = 0 (x,y) E R, BX BY (1.1) U(X.y) = f(x,y) (my) 6 C. The existence of a unique solution u(x,y) to (1.1) has been established by a variety of means [1, p. 7]. However, the exact solution can be determined only in special cases. If R is a rectangle, the solution can be given explicitly by Fourier series. If R is the disk ‘2‘ = |reie| < a, the solution is given by the Poisson Integral Formula 2 2 1 an a - rO i¢ (1.2) u(rO,Oo = 33.} L32 2 “(39 )d¢ o +r -2ar cos(e -¢) 0 o o 16 o for all ‘z \ = ‘r e 1 < a. o 0 Since exact solutions to (1.1) are generally not obtain- able it is necessary to seek approximate solutions to the problem. In this paper approximations to (1.1) by the method of finite differences are considered. The basic idea is to "discretize" (1.1) by replacing the set C = R + C by a discrete point set, and the differential equation by a difference equation. The earliest investigations in this area were concerned with proving the existence and uniqueness of a discrete solution which converges to the solution of the continuous problem as the mesh constant decreases to zero, (Phillips and Wiener, [l3]; Courant, Friedrichs, and Lewy, [14]; Gerschgorin, [15]). More recently the development of high speed computers has greatly improved the actual computing of the discrete solution and has produced renewed interest in the problem of obtaining known bounds on the discret- ization error, i.e., known bounds on the difference between the continuous solution and a discrete solution. In this paper bounds are obtained on the discretization error. The bounds given here differ from those obtained in previous investigation in that they are in terms of known data and do not require bounds on the derivatives of the unknown function u(x,y). Following Bramble and Hubbard [2, p. 314] let G be covered by a grid consisting of two families of lines where each family consists of equally spaced lines parallel to one of the coordinate axes. The intersections of grid lines will be called grid or mesh points. Let Rh be the set of all grid points whose nearest neighbors in the x and y directions lie in R. Let G: denote the grid points in R which are not in Rh' Points of inter- section of C with the grid will be denoted by Ch. Let the continuous Laplacian Au = uxx+uyy be replaced by * the discrete Laplacian Ahu(x,y) appearing in [2]. For the grid point P with coordinate (x,y) shown in Figure 1.1, / O (X aymzh) + P(XFY) A (X-d3h,)’) (X'erh,y) ”(x,y-a h) 1 4 Figure 1.1 * 2 "(X+Ulh:Y) u(x,y+u2h) U(x-a3h,y) u(x,y-a4h) (1'3) A “(XOY) = """ + + + h 2 a1(al+03) a2(az+04) a3(al+u3) 04(az+a4) h 1 l a) ‘ + U(X,y)] 2 (U103 0’20 where O < oi s l, i = 1,2,3,4. For P E Rh all neighbors are at distance h from P, and * AhU(X:Y) reduces to the usual 5-point approximation to the Laplacian, 1 (1-4) AhU(X.y) = ‘3iu(x+h,y) + U(x,y+h) + U(x-h,y) + u(x,y-h)-4u(x,y)]- h A function which assumes values only at the mesh points of the grid is called a mesh function. An example of a discrete analogue to the continuous problem (1.1) is the problem of finding a mesh * function uh(P) defined for P 6 Ch = Rh+Ch+Ch which satisfies (1.5) Ahuh(P) = 0 P 6 Rh, * _ * Ahuh(P) " 0 P6 Ch: uh(P) = f(P) P 6 Ch- Problem (1.5) is a system of simultaneous linear equations for the determination of the mesh function uh(P). It is well Iquwn1that (1.5) has a unique solution for h sufficiently small, [13]. The mesh function e(P) = u(P) - uh(P) is the discret- ization error for the Dirichlet problem for Laplace's equation. Bounds for the discretization error for elliptic boundary value problems usually involve the expressions ‘A:u(P)‘ and ‘Ahu(P)‘. Upper bounds for these terms have been obtained by expanding the function u(x,y) in a Taylor series about the grid point P. Typical of the bounds obtained by this procedure are those given by Bramble and Hubbard [2, p. 314], * * (1.6) \Ahu(P)‘ s 2/3 M3h P 6 ch, (1.7) ‘Ahu(P)‘ g 1/6 M4h2 P E Rh’ where i (1.8) M, =sup $192 1 =O,1,...,j]. PER 5x152: J Bounds on |A:u(P)| obtained in this manner require that the boundary data f(x,y) and the boundary curve C be smooth enough that the partial derivatives of the unknown function u(x,y) are bounded in R. However, a large class of problems involve a function u(x,y) defined on a region R bounded by a piecewise analytic curve C with analytic boundary data on each of the segments of C. In this case the derivatives Of u(x,y) are generally unbounded in any neighborhood of a corner point which is common to two analytic arcs on which the analytic boundary data is given. Bounds for the derivatives of u(x,y) are not known even on a closed subregion of G which excludes corner points. Laasonen has Obtained some information on the degree of convergence of ‘e(P)‘ for several finite difference analogues of Laplace's equation defined on a region bounded by a piecewise analytic arc and having analytic data on each of the arcs [7, 8]. Laasonen's approach requires information on the derivatives of u(x,y). In [6], Laasonen showed that in the vicinity of a corner point formed by the intersection of two analytic arcs the deriv- atives of u(x,y) have order O(r-O) where r is the distance from the corner point and o is a positive real number which depends on the prescribed boundary data on the arcs, the angle of intersection, and the degree of the derivative. The papers of Laasonen do not contain bounds on |e(P)] because of the presence of these asymptotic estimates for the derivatives of u(x,y). Other investigations in this area have yielded in some special cases bounds on |e(P)‘ which depend only on the known boundary data. Bounds have been obtained for |e(P)‘ in terms of known boundary data in the case of a harmonic function defined in a rectangular region [3, p. 309]. Also, Bramble and Hubbard have obtained bounds for \e(P)| for a function u(x,y) which satisfies Poisson's equation on a region which is bounded by a smooth curve [4]. The technique of Bramble and Hubbard requires bounds for expressions which involve derivatives of the unknown function u(x,y). They obtain bounds for these expressions when the boundary and data are sufficiently smooth. In this paper the boundary Of the region is assumed to be a piecewise analytic curve and the data is assumed to have a polynomial approximation on each segment of the boundary. In Chapter 2 new bounds are obtained for ‘A:u(P)‘. These bounds are in terms of known data; they do not involve derivatives Of the unknown function u(x,y). In Chapter 3 these bounds are used to obtain some computable and a priori bounds for ‘e(P)‘ for a discrete analogue of Laplace's equation with continuous poly- nomial data defined on a polygon. (An a priori bound for ‘e(P)‘ is a bound for |e(P)\ which is in terms of known data and does not contain explicitly the inverse of the matrix of coefficients associated with equation (1.5). The a priori bound for \e(P)| is generally a very rough estimate of the bounds for ‘3(P)‘ which can be computed with the known bounds for \A:u(P)| obtained in Chapter 2. However, the a priori bound has the advantage that it gives information on the behaviour of |e(P)‘ as the mesh constant h is decreased.) In Chapter 4 computable and a priori bounds are obtained for a discrete analogue of Laplace's equation with continuous polynomial data defined on a region bounded by a more general analytic curve. The Appendix contains some lengthy calculations and a sample problem. CHAPTER II BOUNDS FOR THE DISCRETE LAPLACIAN OF A HARMONIC FUNCTION In this section new bounds are obtained for the discrete Laplacian A:u(P) of a harmonic function u(x,y). Unlike the bounds which appear in the literature, the bounds given here are in terms of known data and do not involve derivatives of the unknown function u(x,y). Discrete Laplacian for a Function Harmonic l2.§,DiSk Using complex notation, u(x,y) = u(z), where z = x+iy, suppose u(z) is harmonic for ‘2‘ < a, and is continuous on ‘2‘ s a. The Poisson Integral Formula for a disk centered at the origin gives the value of u(zo) in terms of the values of u(z) on ‘2‘ = a. The Poisson Formula given in (1.2) may be written 1 2n z+zo i¢ (2.1) u(zo) = EE'I Re(z_z )u(ae )d¢ o 0 i9 i where z = r e 0, z = ae ¢, r < a. o o o * In the evaluation of Ahu(P) from (2.1) it may be assumed without loss of generality that the grid point P is at the origin and the grid is oriented as shown in Figure 2.1. 8 AV 1’ (0,0211) P(0,0) —; ¢ 7X (-a3h,0) (alh,0) 11(O,-a4h) Figure 2.1 19. J From (2.1), for any grid point zj = r,e where J ‘zj‘ < a, z = aei¢ 2n 2+2 2n _L i =L - i¢ (2.2) u(zj) - 2n g Re(;:;§)u(ae ¢)d¢ 2n £ P(z,zj)u(ae )d¢. With P at the origin, equation (1.3) becomes * u(alh) u(iazh) U(-a3h) u(-iaAh) = —- —-———— — + -———"— _ l +’ l u(O) (“ 1°’3 “2% . Since the discrete Laplacian is a linear Operator, 2n __ . gum = g—ny <§P>u(ae1¢)dq>. O (2.4) h *_ 2 P(z,alh) ‘P(z,-03h) P(z,iazh) P ,O‘ = — + + w ere Ah (2 ) 2 a1(al+a3) 013 ((214013) 012(012-i-Or4) h P(Z"i°’4h) 1 1 F(z,0)] + - + —- . oz4 (“2“4) (011013 0201 By the Schwarz inequality, 2n '_ % 2n . % |§1u(9)| s Egg [AEF(2,0)]2dQ] - [if]; o2(ae1¢)d¢] . Suppose ‘u(ae1¢)‘ s M for O 5 ¢ 5 2n. Then 2n _ % (2.5) ‘A:u(P)‘ sB-fi-j' [AEP(Z:O)]2qu . M. O From (2.5) N. Hills obtained the following result [5]. If u(x,y) is harmonic for \z-zo‘ < a, continuous on lz-zo‘ S a, and [u(z)] S M for points on ‘z-zo\ = a, then at a grid point P with coordinate z , 2 (2.6) [Ahu(P)‘ 5553—13—11 1/88 -h where h is the mesh size and a > h is the radius of the largest disk with center at 20 and contained in G. In the * following, a bound is obtained for [Ahu(P)\ which reduces to (2.6) for the case 01 = l, i = 1,2,3,4. As above, suppose the grid point is at the origin and the grid is oriented as indicated in Figure 2.1. ‘_ ‘P(z,a h) E(z,ia h) ‘P(z,-a h) (2.7) AZF(z,O) = £- 1 2 3 hz a1(al+a3) a2(a2+04) +'a3(a1+u3) 'P(z,-ia h) _. + 4 - ( 1 + 1 ) P(z,O) , 04(az+04) 0103 0204 Z _ 2+2. 1+(-1) 2 P(z,z)=Re[—l:]=Re[ z ]=ReE1+—————]. j ”1 fl :1 1-(z ) 1 - (z > z If |-—1| < 1. Z 10 (2.8) P(z,zj)=1m[1—%L—+ 2 k2 4(11) fink],— *— (2.9) ohms) = -2- D‘IN f_——] ...: + N -((21 + 1 )1 °’l°’3 “2% L. Simplification of (2.9) gives r 3 4 *—- _ 4 1 2 2 h 3 3 h 7 (2-10) AhP(z,0) - ‘5 a m (011-113) 3 cos 3(1) + (almig—4 cos 4(1) + h 1 3 a a 5 6 4 4 h 5 5 h (011-013); COS 5O + (011-1013)} cos 6g) +...] 3 4 +1 2 2 h +03 h azhl‘ [(az-a a4)—3-(- sin 3¢) + (0/: +014)? cos 4d) + a a 5 h6 ' (cg-a:)h§ sin 5¢ +-(a H+u )2 6(- cos 6¢) + .] L a j where ¢ is the variable of integration in (2.5). From (2.5) and (2.10), (2.11) ‘A:u(P)‘2 s 8M h4 2 ll 2 (01‘1013) a 2 (a2+u4) [118 fl 3 3 3 3 12 h 5 5 5 5 -a12 (a1+°‘3) (“24%) 16 h 7 7 7 7 Tami-m3) (“21%) or 2 (2.12) ‘A:u(P)‘2 s d2(h,a,ozi) - M , where d2(h,a,ai) is defined by (2.1 ). For the special case ai = l, i = 1,2,3,4, (2.11) has the closed form (2.6). If a1 = 03, a2 = ah, i.e., the grid is rectangular, then from (2.11), 2 4 4 2 2 2 k * 2/2 h M °‘1 0’2 °’1°’2 ”fans 2 444+ 444* 4224 a (a “(1111 ) (a 'azh ) (3 'kl’lazh ) After further algebraic simplification (Appendix I), the formula above reduces to * 2/2 h2M(ozi-I-Ot:) (ai-a§)2h4a4 55 (2-13) IAhuml ‘- ____.___ 1 + 4 a 4 4 4 4 [a8_aaa4h8 (a -oz1h )(a ~azh ) l 2 12 * From (2.11) it is clear that |Ahu(F)| is 0(h) when 0’1 1 0/3, or a2 # a4. * The grid point P 6 Ch will be said to be of Type 1 if it has at least one neighbor in each of the x and y directions, e.g., a1 = a2 = l, O < a3 5 1, O < a S 1. In this case a 4 closed form upper bound for (2.11), independent of ai’ is . k (2.14) [Ahu(P)‘< 4 h M [l +8h::] , (Appendix I). a/a4_h4 a To simplify the discussion, it is assumed in the remainder of this paper that all P E c: are of Type 1. The above results are summarized in the following Theorem I. Let P be a grid point with coordinate 20. Suppose u(z) is harmonic in the disk ‘z-zo‘ < a, continuous on ‘z—zo‘ s a, and satisfied ‘u(z)| S M on Iz-zo‘ = a. Then if h < a, ‘A:u(P)\ S d(h,a,ai) ° M where d(h,a,ai) is defined in (2.11). If all neighbors of P are at distance h from P, 2 (2.15) \A u(P)] s fl—h—M. h f 8 8 a -h If the mesh is rectangular; i.e., a1 = a3, a2 - a4, ‘ * ‘ Z/2 M 112612-1012) (02 _a:)2haa: % (2.16) u(P) s 1 + z. a. z. a. 4 Ah (as: a4_a4h8 (a -0111. )(a ‘02 h ) 0’1 2 If P has at least one neighbor at distance h in the x and y directions, 2 $5 (2.17) \A:u(1>)| < iii-[:1 +§h—2] . a/aA-ha 13 The results of Theorem 1 may be readily applied to a function which is harmonic in a general region. Suppose the curve C forms the boundary of a region R in the (x,y)-plane, and suppose Au(x,y) = O in R and ‘u(x,y)‘ s M on C. By the maximum principle for harmonic functions [u(x,y)] s M for all points (x,y) in R. Let P be a grid point in R, and let a denote the radius of the largest disk which is centered at P and con- tained in R. If a > h, the discrete Laplacian A:u(P) satisfies the bounds given in Theorem 1. The importance of the new bounds given in (2.15) to (2.17) is that they are in terms of known data and do not involve deriv- atives of the unknown function u(x,y). Grid Points Near the Boundary The a in the denominator of the expressions of Theorem 1 may always be taken to be the radius of the largest disk with center at P which is contained in G. Hence, the bounds for lA;u(P)l are best for those grid points which lie deep inside the region R. On the other hand, if the grid point P is near the boundary, the bounds are poor; if the grid point is within h units of the boundary, the bounds are not applicable. A technique involving the Schwarz Reflection Principle may be used to extend and improve the results of Theorem 1 for mesh points near the boundary when the boundary is an analytic curve and the prescribed data is a polynomial. The procedure for the important special case when the boundary is a straight line is discussed below. The more complicated problem of a general analytic curve is discussed in Chapter 4. 14 The reflection principle as it applies to a disk with center at the origin may be stated as follows: If u(x,y) is harmonic in D+(p) = {(x,y)‘x2 +-y2 < p2, y > O}, and u(x,0) = O for -p s x s p, then u(x,y) is harmonic in the disk D(P) = [(X,y)‘x2 +y2 < p2]. Moreover, for (x,y) E D-(p) = 2 {(x,y)1x2 + y2 < p . y < 0}. U(X,_y) = -u(x.-y). Zero Data If the mesh point P with coordinate 20 is near a straight line V1V2 oriented to coincide with the x-axis, and if u(x,O) = O on V1V2, then the a in the bound for \A:u(P)‘ given by Theorem 1 may be taken to be the radius of the largest disk centered at 20 and contained in G +-G where G is the set in which u(x,y) R R is known to be harmonic by reflection in V The M in the 1V2' formula is not changes since [u(z)] s M for all z in the disk lz-zo‘ S a, Figure 2.2. Au=O 111111111 1111111111’ \ v U1(x,0)=0 v ’x G R Figure 2.2 The reflection principle may also be used to improve bounds * for \Ahu(P)‘ obtained from Theorem 1 when the grid point P is near a corner formed by the intersection of two straight lines. 15 Suppose u(r,9) is harmonic in the sector of angle a = n/m, m an integer, with vertex at the origin, and sides (r,O), (r,n/m), O s r s 9. Suppose also that |u(r,9)‘ s M in the sector and u(r,O) = u(r,n/m) = O for O S r s p. By repeated application of the reflection principle u(r,e) is harmonic in the disk D(p) = {(r,O)‘O s r < p, 0 S 6 s 2n}, and ‘u(r,9)| S M in D(p). If the angle of the sector is n/m where m is not an integer, u(r,9) may be reflected across the sides of the sector, but in this case the origin is a branch point of the function. Non-Zero Data The improvements on Theorem 1 obtained above for grid points in a region whose boundary contains straight lines and corners require that the harmonic function u(x,y) have zero boundary data on the lines. By subtracting from u(x,y) an appropriate simple harmonic function, it is possible to obtain improved results for ‘A:u(P)‘ when the grid point P is near a line with non-zero data. In the technique described below it is assumed that the boundary is a straight line and the prescribed data is a polynomial. The prOblemsof non-zero data at corners and more general boundary data are discussed later. Suppose the region R is oriented so that the grid point P has coordinates (O,yo) and the boundary segment V coincides 1V2 with the x-axis, -x1 S x 5 x2, x1,x2 > 0, (Figure 2.3). Suppose the prescribed data for u(x,y) on VIV2 is u(x,O) = PN(x) = 2 N + . ' , . ao 81x +-azx +3..+-aNx Let p > 0, p < mln(x1 x2) Let GR denote the region in which u(x,y) is known to be harmonic by l6 * reflection across lines with zero boundary data. Let G = G +-GR. Suppose that u(x,y) is harmonic in D+(p), and |u(x,y)| S M in D+(p), where D+(p) = [(x,y)‘x2 + y2 < p2, y > O}. f-“"‘-_ _'"“_"'—"l I y G l ' A R I ‘ I ' l /\ ”=0 ."'7” ’7r14717’ 4'77rif7/f7777 I j / u=0 / j 11: r/ [RO’YU) G 5 [[111 [L111/ 41/! l/Jl’//1 1 Il/A «....LX v 0x 0) V ‘V(x’0) 1 1’ 2 2’ u(x,0)=PN(x) Figure 2.3 a = p-yo Define AO(X:Y) = 80, n An(x,y) = anRe(z ), n = 1,2, N, DN(X,y) A0(x,y) + A1(x,y) +...+ AN(x,y). The function DN(x,y) is harmonic in D+(p) and assumes the values DN(x,0) = PN(x). The function W1(x,y) = u(x,y) - DN(x,y) is harmonic in D+(p) and is zero for -p s x S p, y = 0. Hence by the reflection principle, w1(x,y) is harmonic in D(p) = [(x,y)]x2 +-y2 é p2]. Since DN(x,-y) = DN(x,y), the function W1(x,y) assumes the values w1 . (x,y) e D+

W1(x,O) = 0 , -p s x s p , y = 0 w1(x.y) = -u . (x,y) e D’ h. The disk ‘z-yo‘ < a is con- N * * tained in D(p). ‘AhDN(P)| s nEo‘AhAnfl)‘. Bounds for each of * the expressions ‘AhAn(P)‘ may be computed directly. For a given * mesh constant ho, and for a given mesh point P, constants K and KP may be determined such that for all h 5 ho’ |A:DN(P)\ < KZh, \AhDN(P)\ < Kphz, (Appendix II). Finally, applying Theorem 1 to W1(x,y), 2 35 (2.18) 1A:u(P)I < 4h(M+M(p)) (1 + gig) + K*h P e 0:, 4 4 a P a -h 2 (2.19) \Ahu(F)\ < 4/2 h (“+M(°)) + K h2 P e Rh. Ja8-h8 P Bounds for ‘A:u(P)‘ and \Aho(F)| given by (2.18) and (2.19) may be computed for several values of p, and an "optimum" may be selected. An alternate form of (2.18) and (2.19)vnjl.be obtained next. Let W2(x,y) =‘W1(x,y) +-DN(x,y). The function W2(x,y) is harmonic * in G + D(p) with values given by 18 w2(x.y) = u(x,y) (x,y) e 6*. W2(XaY) = Pfi(x) -p S x S p, y = O, w2(XFY) = ‘U(X,'Y) + 2DN(an) (X,Y) 6 D-(p)- Hence, if a is the radius of the largest disk centered at P * and contained in G +-D(p), (Figure 2.4), then 2112.5“ 4h (M+2M(o)) 2) P 6 CE, * (2.20) [Ahu(P)‘ < (l + a -h a 2 (2.21) [Ahu(1>)\ s 4/2 h (“4'2“”) P E Rh. a8-h The advantage of the expressions in (2.20) and (2.21) is that it is possible to have a > p. For example, for the grid point P in Figure 2.4, a = p2 2 . ‘V \\\\§_\ \‘ C l O 11.41/11] 1111/11” OV>X Figure 2.4 a =./ 2 72 Remark: The function DN(x,y) defined above to remove the non- zero boundary data is not unique. DN(x,y) could be replaced by any function of the form DN(x,y) + F(x,y) where F(x,y) is 19 + harmonic in D (p) and is zero for -p s x s p, y = 0. Uniform Bound for ‘A:u(P)] and \Ahu(P)|. Let G be a polygon and let 5’ be a subset of G which excludes some neighborhood of each of the vertices, and let ho be a mesh constant chosen small enough so that all points in C are at distance greater than ho [from all vertices of the polygon. From the discussion above it follows that if the boundary data on G is a polynomial, then real numbers B: and BI may be determined so that for all h 5 ho, and for all P E 5.0 Gh’ * * 9: (2.22) ‘Ahu(P)‘ < BIh P 6 Ch’ (2.23) ‘Ahu(P)| < 31112 p e Rn' The bounds (2.22) and (2.23) will be used in Chapter 3 in an expression for an‘g priori bound on ‘e(P)‘ for a discrete analogue of a Dirichlet problem defined on a polygon. Discrete Laplacian for 3 Function Harmonic i 3 Sector It is possible to start with other representations of a harmonic function and obtain bounds for ‘A:u(P)‘ in a manner similar to the one described above. The Fourier Series represent- ation of a function which is harmonic in a sector and zero on both sides of the sector yields an expression for ‘A:u(P)‘ which is useful for polygonal regions. The exact expression is not in closed form, but closed form upper bounds may be computed. Suppose u(r,6) is harmonic in a sector with vertex at the origin and sides (r,O) and (r,n/m), O s r s p, m > 1/2. Suppose 20 also that u(r,e) is continuous for O S r S p, O S 9 S n/m; u(r,O) = u(r,n/m) = 0 for O S r S p, and u(peie) is pre- scribed for 0 S e S n/m. Then for all points (r0,eo) in the sector the function u(ro,60) has the Fourier Series representa- tion m r km (2.24) u(r ,9 ) = 2 a <-2) sin km 9 where o o k=1 k p o 2m‘n/m i¢ ak=r£ sin km¢ u(pe )d¢, ro< p. In complex form, a, km ie (2.25) u(zo) = 2 ak(l) Im(zo)km where 20 = roe 0, ‘20] < P. on km (2.26) A:u(zo) = Zlak(%) A:[Im(zo) |. k: For the grid point 20 indicated in Figure 2.5, 20 = roe 1W . i? _ 1W 21 z alhe , 22 z lazhe , 23 z a3he , 24 = 20 - iaahelw, where W is the angle between the grid line z z and the x-axis. Suppose that [2}] < 1. o Z M A71u111111\\\11\\\\\\\1\L1 ‘X I r u(r,0)=O (P50) Figure 2.5 21 * km 2 (zo-lulhelv)km (z +iazhe1¢)km (2°27) Ah[Im(zo )] = —2- Im a (0’ +0) + C20! +0 ) h l l 3 °’2 2 4 (zo-theiw)km (zo-iaaheiw)km 1 1 2km + + - ( + ) 0 a3(o'l+n/3) 014012-1114) 011013 01204 = Em. = km km-l = kn(km-l)...(km-(n-1)),... Let Ck1 1 , Ck2 2! ,..., Ckn n! Then (20+alhei*)km zzm m alheiw n (2.28) =——1+zc( ) . “Hal“? O’1("’1"‘°’3) n=1 k“ zo (7.0+iozzhei‘i)km 21:“ so 10/ he 1 n =—— - .. < ) , C”2("2""’4) “2(“2+"4) n=1 1‘“ 20 (z -(:1I3hei¢)km 2km 0° -a3hei¢ n clam) =o'(;+cr)1+zckn( z ) . 3 l 3 3 l 3 n=1 0 (z 'iUAheiw)km 2km w -ia heiv n Note that if m is an integer each series in (2.28) is finite and it is not necessary to require ‘2—] < 1. This Special behaviour for m an integer is discussed in detail later. Substituting (2.28) into (2.27), F km 1 2 iv 3 iU 4 0 he 2 2 he 3 3 61le Ck3( z0 ) (“1"”3) + Ck4( 20) (“1“3) i¢5 i¢65 ] + Ck5(h£——) (oi-03) + Ck6(E§-) (a +a5)+..] 20 o l 3 J * km 2 - (2.29) [Im(z )=-—Im km . . Ah 0 1 h2 +zo c be” 3 , 2_ 2 + c hew 4 3 3 0'2'1'014 k3( 20) (‘1)(“2 0’4) k4( z ) (“24%) ii: 5 iv '7 . 4 4 h 5 5 ‘ . + ck5(hezo) 1(02-014) + Ck6< : ) (-l)(02+u/4)+J 4 O 22 If An = (km-n)9o + nw, the expression (2.29) may be written F "1 km 2 2 T 2r 0! 1! ) ( 02 2) * km 0 h 3 l 3 . 2 +014 (2.30) Ah[Im(Zo )] = hz Ck3(r—-o ) W Sln A3 + 7‘:— COS A3 ..J 3 3 O h 465W oil-+63) ((124114) —— ________ ‘ + —_ ' + C 4(r )4 a +0 Sin A4 a +0 Sin A4 1 3 2 4 J 4 4 4 4 n (-'—) ——-——- sin A + ———- cos A r0 a1+u3 5 012+Ot4 5 a a 5 5 4 4 (h 6 [(011%) (‘az'afl 1 + C “) -—-—-—- sin A +--——————- sin A k6 ro al+u3 6 a2+u4 6 L J + L .J . Let S(k) denote the expression in the outer brackets in (2.30). From (2.26) and (2.30), * 2 m to km (2.31) Ahu(zo) = ;§.kilak(3—) S(k). Applying the Schwarz inequality to (2.31) (2 32) \Atho )I2 54— (282) ' (<320 :9 2km(S(k))2) h4 (k=1 k k=l(p ) By Parseval's equality, Hm'n/m 2 (2.33) = Ef-I u 2(pe i¢)d¢ s 2M where M is a bound for o ’ITMS ...: me W E \u(pe )l. 0 s o s m- 19 Let P denote a mesh point with coordinate 20 = roe 0. Then 2 co 2km * 2 8M 2 (2.34) \Ahu(P)| 5—2; 2(4") (S(k)) h k=l\p As a special case of (2.34), if ai = l, i = 1,2,3,4, then 23 h 4 h 8 , S(k) - ck4(ro) sin A4 + ck8(ro) 8111 A8 +... . From (2.30) and (2.33), 2 m r 2km 2 32M 4 . (2.35) IAhu(P)| s h“ 21(‘3-9) [ck43;‘;) 8111 A4 2 h8 . +... + Ck8(ro) Sin A8 1] . Two types of closed form Upper bounds for ‘A:u(P)\ and |Ahu(P)[ are derived below. The first is easy to compute, but it is not very accurate. The second is more accurate, but it requires some preliminary calculations. If P E c: is of Type 1, the term S(k) satisfies h 3 h h 2 (2.36) \S(k)\ < 4 (r0) - [10161 + ‘Ck4‘(ro) + ‘Ck5‘ (:3) Therefore from (2.34) 2 2 m r 2km (2.37) MEMPHZ < lagyi z (.9) [1016‘ + [ck4\(£-) 0 r k=l\p 0 2 h 2 +\ck5\(ro) +...] h 1 As shown in Appendix III, if ;—'S 23 then °° h n 3 2 2km ° [ 2 ‘C —') ] < e for all integers k and all m > k. —3 (2.38) \A:u(P)\2 < T E 0 <=128M____6__2__h2(:g_p e)2m[pzm 92m m-(ro e)2m ’5 -——(:”.>[.: e. 1 . p r 2km 128M2h2 m (roe) l 24 Similarly, an upper bound for ‘Ahu(P)‘ is 2 r e m 2m % r 4/2 h M.( 0 Q, o l, h l (2.40) ‘Ahu(P)‘ < 4 p 2m 2m , for p < e’ r S 2 r p -(r e) o o o r h—-s %‘ are not the best possible bounds. In this case, and in similar cases which will appear later (The bounds 32-< %, and in this paper, the question of best possible bounds is not significant and the indicated bounds are simple or convenient rather than the best possible.) The general behaviour of \A:u(P)\ and \Ahu(P)‘ for various values of h, re, and m is apparent from (2.39) and (2.40). More * accurate bounds for ‘Ahu(P)| and \Ahu(P)‘ may be computed using 19 the techniques described below. It is assumed that P(roe 0) r satisfies 32-5 %3 E— S‘%. These bounds will be used in the re- mainder of this paper. * Suppose P 6 Ch is of Type 1, then 01+03>1, q2+04>13 from (2-30), 2m Pk -1 2(k-l)m 2 2 r o r * 2 32M h o o (2.41) u(P) < -——- (—) z (—) (lck3| + \CRS‘TZ + |ck7|T4 +...) 2 + 2(‘Ck4‘T + \ck6|T3 +... ) 2 ’ k-l 32M2h2(ro)m ° (22)“ )m + -_—__—. — k 4 2 r6 p k: p o o 2 2 . [lcml + |ck4|T + \ck5|T +...] L . . . . h 1 where R0 18 a large p051t1ve integer, T = ;-'s'§. 0 Similarly, 25 2 4 2m ko'l 2(k l)m 32M 8h 0 (2.42) \ u(P)‘2< ___(—) 2 (.2) Ah r8 p k=l p 0 2 .[\ck4| + [ckgl'r4 +|ck12\'r8 +...] :l 2 m 2 k-l +__,_...22(.)2 2(2)‘ ”“ r p k=kop O '[lck4| + ‘ckal'r4 + \cklleB +...] 2] Define 4 —9 (‘ck3\ + |ck5|T2 + ‘Ck7‘T +...) (2.43) w*(m) - l M kO-1(r )2(k-1)m 1 5 + 2<|Ck4lT + ‘Ck6‘T3 + ‘CkB‘T +...) * an r 2(k-l)m 2 2 (2.44) R (m) = 4 kzk(p—°) [lck3| + le4|T + Icksh +...] 0 ko'1 r0 2(k-1)m 2 0-45) won) =k§1(;-) [\de ”C228” ”0121le + J ’ a. r 2(k-l)m 4 8 2 (2-46) R(m) = kill?) UCM‘ + |Ck8|T + |ck12|T +...] Then (2.41) and (2.42) may be written 2 2 2m 2r (2.47) \AEU(P)|2 < 32M6h (3) [W* (m) + R* (m )]= --331—h-(p- “321“) [6*(m)] r. r O '1 p O p O O (321%?4 o _ 32M2h4 ro M (2-48) ‘AhU(P)|2< ——("m [W (m) + R(m) J - T(3—) [G(m)] . r8 r r For given values of k0,39, and m, bounds may be computed * * for W (m), W(m), R (m), and R(m). Some computed bounds for * C (m) and G(m) are given in Table l for various values of m 26 r and ‘-9. The values given in Table 1 are valid for all grid points satisfying E—-s %; k0 was chosen so that kom 2 150. For the special case where m is an integer, the expressions (2.47) and (2.48) may be written without ro in the denominator. For example, if m = 1 in equation (2.41), then 2 r 6 6 k '1 2k-6 * 2 3 2M _9) L) O (2.49) ‘Ahu(P)\ < -;Z—(p (r k 0 =3 2 4 2. (\Ck3| + [ck5|T + \ck7|T +...) + 2(‘ckalT + |ck6|T3 + ‘Ck8‘T5 +.. ) _J m 2k-6 + 22M2(:Q)6 L)6 4 (:2) 4 (r z p h p o k=ko 2 ' 0le +1Ck4‘T + leS‘TZ +...] 2 2 (2.50) |A:u

|2< ””2 “2 The" (1) + 12* (1)] = 32” 6h [6* (1)] 96 Similarly, (2.51) lahu(P)\2< 32M——2— 24 [6(1)] 98 In (2.50) and (2.51) the series involving T terminate and it is * not necessary to require T = %—'< 1. Some values of C (m) and o G(m) for m = 1,2,3,4 are given in Table 2. In these computa- tions it was assumed T s 1; k0 was chosen so that kom 2 150. 27 * Table l. Bounds For C (m), and G(m). 22-22 2922 c*h2 +‘AhDN (P)| P e Rh 0 3. r0 33 Let h = ho be an initial value for the mesh constant. Con- * stants B (m) and B(m) may be determined so that for all mesh 1 points P E S with 2"3'3, and for all h s ho, o F" * 1 B BThh §-< m2< 3 r O * * 5 * * (2.75) |Ahu(P)| < B (3)h(1 + b ‘log rol) m = 3 . P 6 ch 0 s, * LB (m)h m > 3 * where b is a known positive constant, and B m h2 1 4% '2'< '“< 4 r O (2.76) ‘Ahu(P)l < B(4)h2(l + b‘log ro|) m = 4 P e Rh n s, 2 B(m)h m>4 where b is a known positive constant. The inequality (2.75) is proved below. The proof of (2.76) is similar. The bounds in (2.75) and (2.76) will be used in Chapter 3 in an expression for an‘g priori bound for the discretization error for a discrete analogue of a Dirichlet problem defined on a polygon. Also, the bounds for ‘A:u(P)‘ written in the form of (2.75) and (2.76) may be readily compared with some known results concerning the discrete Laplacian of a harmonic function, (Remark 3). Proof of (2.75). If sin 11% r 0 for n = 1,2,...,N, the be- * haviour of ‘AhDN(x,y)‘ is determined by the behaviour of * n * n * ‘AhRe(z )| and ‘Ah1m(z )|. In this case a constant KD may 34 * be determined such that for all P 6 Ch 0 S, and for all h 5 ho, * * |AhDN(P)| < KDh. 1:. — 1 4 * n *h If sin n m - 0 on y for n 2 , then ‘Ahlm(z log z)‘ < Ln h l * for all ‘;- S‘2 (Ln a known constant), and again a constant * ° * , h 1 KB may be determined so that for all P 6 Ch 0 S Wlth ;—’$'§s * 'k and for all h s ho , [AhDN (P)|2< Kbh. * Case 1 ) Suppose -'< 3, i.e., m.> 3. In this case Im(znlog z) is a part of DN(x,y) only for n 2 4. Hence in the sector S where 2h s r0 5 R * (2.77) ‘A:u(P)‘ < 110?? + th < hERm-3T*(m) + [(3] = B*(m)h. r o 11' * Case 2 ) Suppose mi: The worst possible behaviour of gold * ‘AhDN(P)| occurs when Im(z3log z) is a part of DN(x,y). In this case (2.78) lA:u(P)‘ < T*(3)h + th + Clh‘log ro| where Cl is a known positive constant. (2.79) ‘A:u(P)‘ < hEr*(3) + K: + Cl‘1og r01] C * * * = B (3)h[1 + b Slog ro‘] where b = ___1__ T*(3)+K; * n n . n Case 3 ) 3'- 2’m "#'n. Again, the worst possible behaviour of DN(P) occurs when Im(zzlog z) is a part of DN(x,y). This can n=_31 [E TI ‘1 occur only for m 2 . For all angles m> 2 , m f n, 9: (2.82) ‘A:u(P)‘ < Lg? + x; g- = 3-m [T*(m) + Kgri'm]. r r Since m < 2, * h * * 2- B (2.83) |Ahu(P)| < 7:; [T (m) + KDR "j = 0 r0 Case 6*) %= 11. This is the only case where Im(z log 2) may be a a t f D ( ) S' e S *Im(z 10 z)‘ < * E— p r o N x,y . inc Ah g KD r2 , o * T* (l)hr: (2.84) \A:U(P)‘ '% -- if u(x,y) is continuous at the corner, then u(x,y) can be written in the form u(x,y) = u1(x,y) + u2(x,y) where all partial derivatives of u1(x,y) remain bounded as the point (x,y) approaches the corner. If r denotes the distance from the corner, then u2(x,y) has the order of magnitude m C (8(r ) a = H-, m not an integer, m (2.85) u2(x,y) =4 0(rmlog r) a = 3', m an integer. L The relation in (2.85) may be indefinitely formally differ- entiated. From (2.85) and (1.6) Wasow's result implies that in the vicinity of a corner of the type described above * h (2.86) |Ahu(P)‘ = 0 ( 3-m) m not an integer, r * h log r . (2.87) ‘Ahu(P)\ = O ( 3-m ) m an integer. r It follows from Remark 2 above that for the special case of the intersection of two straight lines, where the boundary data is continuous at the vertex and is a polynomial on each of 37 * the sides, the asymptotic estimates for \Ahu(P)‘ in (2.86) and (2.87) may be replaced by known bounds. CHAPTER III THE DISCRETIZATION ERROR FOR A DIRICHLET PROBLEM DEFINED ON A POLYGON Let G be a polygon in the (x,y)-plane with interior R and boundary C. In sections 3.0 to 3.3 computable and a priori bounds on the discretization error |e(P)| are obtained for a Dirichlet problem defined on G. 3.0 Error From Polynomial Approximation of Original Data By the superposition principle for harmonic functions it may be assumed without loss of generality that any Dirichlet prob- lem with continuous boundary data defined on a polygon is zero at all the vertices. Let fk(x,y), k = 1,2,...,ko, denote the prescribed con- tinuous data on the ko sides of a polygon, where fk(x,y) is zero at all the vertices of the polygon. (If the original data fk(x,y) is continuous on a side k, but has a discontinuous derivative at some point (xo,yo) in the interior of a side k it is possible to introduce a new "vertex" at (xo,yo) with angle Ef-==n, and hence a new "side" of the polygon.) Tte problem of determining a polynomial approximation to fk(x,y) is not considered in this paper. It is assumed that on each side k a polynomial Pk(x,y) and a positive real number ek have been determined so that Pk(x,y) is zero at the vertices, has continuous derivatives of all orders on the side k, and 38 39 satisfies. ‘fk(x,y) - Pk(x,y)‘ = IRk(x,y)\ S ek for all (X.Y) on the kth side of the polygon. Let uP(x,y) denote the solution to Laplace's equation on the polygon with Pk(x,y) boundary data on each of the k sides, and let uR(x,y) denote the solution to Laplace's equa- tion with Rk(x,y) data on each of the k sides. From the maximum principle for harmonic functions, ‘uR(x,y)| s eo = 3:? ekk for all (x,y) in the polygon. If uPh(Q), defined on some ,..., o mesh G , is a discrete approximation to uP(x,y), and bounds for h ‘e(Q)‘ = ‘uP(Q) - uPh(Q)‘ have been determined on G then b) (3.0.1) |u(Q) - uPh(Q)| ‘U(Q) - upm) + uP(Q) - uPh(Q)l V\ M» - uP\ + lqu) - uPh(Q)| V\ eO + |e(Q)| for all Q e Gh' 3.1 Computable Bounds For The Discretization Error Denote the RC vertices of the polygon by V1,V2,V3,...,Vk , o and denote the interior angle of the polygon at Vk by {g; , >1 “‘k 2' Let u(x,y) denote the continuous solution to the follow— ing problem defined on G: (3.1.1) Au(x,y) = 0 (x,y) E R, U(x,y) = f(x,y) (x,y) 6 C- In accordance with the discussion above, it is assumed that the data f(x,y) is continuous and is zero at the vertices Vk’ k = 1’2’°'°’ko° On each of the RD sides of the polygon, f(x,y) 40 is assumed to be a polynomial with continuous derivatives of all orders. Let h be given and let G be a square mesh of width h h it that covers 9. By the methods of Chapter 2 a bound for \Ahu(p)| may be computed at all grid points of Gb n R whose distance from a vertex is greater than h. In this section the boundary data is interpolated in the vicinity of each vertex V and the point sets k, * * . Rh and Ch are defined so that for all P E Rh + Ch the distance from P to the nearest vertex is at least 3h. (The procedure used in section 3.3 to determine anla priori bound for ‘e(P)‘ . h_ ‘1 * . . requires r S 3 for all P E Rh + Ch’ where rO is the distance 0 h 1 * . h . —- — + . from P to Vk T e requirement r s 3 for P E Rh Ch is imposed here in order to simplify the discussion in section 3.3.) Let Chk denote the set of all grid points in Gh n R whoie o . - - I = l distance ro from Vk satisfies r0 < 3h, and let Ch kilchk. Denote by Rh those grid points in (Ch 0 R) - CA with all * neighbors in R - Cg. Let Ch = (Ch O R) - CA - Rh’ and let Ch be the set of all points in Gh n C which have a neighbor in * Ch’ (Figure 3.1). Rh _ o * __ Ch"A I Chk D ch—x V k Figure 3.1 41 The following discrete analogue of (3.1.1) will be con- sidered; * * Ahuh(P) = 0 P E Rh + Ch (3.1.2) uh(P) = f(P) P 6 ch uh(P) = f(Vk) =0 P E Chk' The assignment of boundary values for P E C' is known hk as boundary interpolation of degree zero. (In some problems it is possible to compute \A;u(P)| for all points P E G, and hence it is not necessary to interpolate boundary data. For example, the set CA is empty in the sample problem for the square given in the Appendix.) The system of linear equations (3.1.2) has a unique solution uh(P) for h sufficiently small. A procedure will be given for computing a bound for ‘e(P)‘ = |u(P) - uh(P)‘ for all mesh points in Rh + CE. Suppose the linear system (3.1.2) involves m points in * R + C and n points in C + C'. The system (3.1.2) may be h h h h written in the form (3.1.3) A uh = A Bh where A is an m X m matrix, A is an m X n matrix, and Bh is the n X 1 matrix of boundary values of uh for points of Ch +'Cfi. The continuous solution u(x,y) satisfies the equation 2 * . (3.1.4) A u = h (Ahu) + A B * where (Ahu) is an m x 1 matrix whose elements are the discrete Laplacian of u(x,y) evaluated at each of the m points of 42 Rh +-C:. The matrix B in (3.1.4) is the n X 1 ‘matrix of boundary values for u(x,y). Let 6(Q) = B(Q) - Bh(Q) = u(Q) - Uh(Q) for Q 6 Ch + Cfi. (3(Q) = 0 except at points Q 6 Ch.) From (3.1.3) and (3.1.4), 2 -l * -l- * (3.1.5) e(P) = h A (Ahu) + A A e(Q) P E Rh + Ch , or 2 -l * -1. * (3.1.6) 9(1)) = h (-A )(-Ahu) +A A e(Q) P E Rh + ch It is known that all elements of (-A-1) and (A-lA) are non- negative [1, p. 25]. Hence ‘e(P)‘ satisfies (3.1.7) ‘e(P)‘ s (-A'1)[h2\A:u\ + (-A)|e(Q)\] P e Rh + Ch’ * From Chapter 2, bounds for ‘Ahu(Q)\ may be computed for all * Q E Rh + Ch. Bounds for ‘e(Q)|, Q 6 Ca will be given in the next section. Then a bound for ‘e(P)‘ may be computed from * (3.1.7) for any point P E Rh + Ch. 3.2 Bounds For Data Interpolation Error It may be assumed without loss of generality that the vertex Vk is at the origin and the sides of the polygon lie denote the radius along the rays (r,O), (r,E-), r 2 0. Let p mk k of the largest sector F with vertex VR and sides (r,O), k (r,E-), 0 s r 5 pk and such that F is contained in the polygon. k Suppose ‘u(r,e)| s M on the boundary of the polygon. Case 1.) u(r,9) = 0 for (r,O), (r,§;), 0 s r 3 pk. In this case the Fourier Series representation of u(r,e) in Fk is 43 r jmk (3.2.1) u(r,e)= 2 aj(— k) sin jmke where j= -l E. 2mk mk aj = TI u(pk,¢)sin jmk¢ d¢. o By the Schwarz inequality, and the Parseval equality, a: % co" j 2 15 2 r mk (3.2.2) u(r,e) S 2 2 ‘— sin ' 9 ‘ ‘ [Haj] [i=1((°k) Jmk 0) J m ij 35 f2 2 5— k] . s M[j=1(pk) 2mk m p % (3.2.3) ‘u(r,e)‘ s/Z M ('3‘) k[_z_fiL—2m—k:] , £- < l. k pk k-t pk From (3.2.3) it is evident that a positive constant Mk’ inde- pendent of h, may be determined such that “‘k (3.2.4) \e| =1u - uhml = mm s Mkh Q e cgk. Case 2.) N u(r,0) = alr + azr2 +...+ aNkr 1‘ (3.2.5) n 2 Nk u(r,'n'l‘;) = blr + b2r +...+ka1' 0 s r 5 pk. Let DN (r,e) denote the data removal function defined in k Chapter 2 (p.2fl ), and let M(pk) = max‘D (r,e)|. In this case Fk Nk the function W(r,e) = u(r,e) - DN (r,9) has the Fourier Series k representation a: r jmk (3-2-6) V(ra9)= 2 b. (——) sinjmke, where j=1 J pk 1T— 2m "‘k b] =-'k‘f W(Pk,¢)sinjmk¢ d¢. o 44 Hence w MY‘ (3.2.7) ‘W(r,e)l Sf? (M “NPR” (i’) k 2m 2m From (3.2.7) 3 positive constant A independent of h may be k determined such that m. (3.2.8) |W(Q)| s Akh 1‘ , Q g Chk' . I For all pOints Q E Chk’ (3.2.9) |e(Q>l = mm s \w\ +1sz (on. ‘k Let h = ho be chosen initially. The worst possible behaviour of DN (Q) occurs only when E-'= n and Im(z log 2) is a part k of DV (r,e). In this case, there are positive constants M L k and b0 such that for all h S ho, k (3.2.10) |DNk(Q)\ s Mk(h‘log h‘ + boh)’ Q 6 Chk‘ In all other cases a positive constant M may be determined k such that for all h 3 ho (3.2.11) \DNk(Q)‘ s Mkh , Q 6 Chk‘ Summarizing, if Ef'f n, a positive constant fig may be k determined such that for all h s ho, a — k ' (3.2.12) \e(Q)| s Mkh Q 6 chk, where ak 2 min(1,mk). If Ef'= n, and Im(z log 2) is a part k of D (r,O), positive constants fik and b may be determined Nk l 45 such that for all h s ho, (3.2.13) |e(Q)| s iik(h\1og h| + blh) Q 6 Chk' 3.3 A Priori Bound For Legal For any value of h, a bound for the discretization error ‘e(P)\ = ‘u(P) - uh(P)‘ for the discrete problem (3.1.2) may be * computed from (3.1.7) for all points P E Rh + Ch' In this section an a priori bound is obtained for the error ‘e(P)‘ computed by (3.1.7). This bound will also give information on the behaviour of ‘3(P)‘ for decreasing h. Suppose that an initial mesh constant h = ho < l is chosen sufficiently small so that all P E c: are of Type 1, * (p. 12). The sets Rh and Ch defined above will be Sub- divided. For k = 1,2,...,k , let 0 pk be the radius of the largest sector Fk which has its vertex at V is contained in k, the polygon, and whose sides coincide with the boundary of the P polygon. Let Rk =‘Sh , and let Sk be the sector with radius Rk’ vertex Vk, and whose sides coincide with the sides of the * polygon. For P E Rh + C , define th = {P‘P E Rh 0 Sk}, * * Chk - {P1P 6 Ch H 8k}' LGC RhI and C111 denote those pOints * * 1n Rb and Ch which are not in some th, Chk respectively. * Suppose that h is small enough so that the sets th, Chk’ * RhI’ and ChI are not empty, and ‘log Rkl < ‘log 3h|. * From Chapter 2 (p. 19), positive real numbers BI and BI may be determined so that for all h 3 ho, 46 'k 'k * (3.3.1) ‘Ahu(P)‘ < BIh P 6 ch (3.3.2) ‘Ahu(P)\ < BIhZ P e R * Also from Chapter 2 (p. 33), positive real numbers Bk(mk)’ * Bk(mk)’ b , and b, may be determined so that for all mesh points P E Sk’ F‘ * h skimk) _ < 3 3‘mk 2 < mk to * P E Chk. (3.3.3) \ *u(P)\ < 8*(3)h(1 + h*|1og r |) m = 3 Ah k 0 k * h > 3 L Bk(mk) mk " Bk(m )h2 k —~ 4 "ZTE;“' 2 < mk‘< 1.0 P e th. (3.3.4) thhu(p)\ < Bk(4)h2(l + b‘log ro‘) m = 4 2 _Bk(mk)h mk > 4 * Recall that in (3.3.3) and (3.3.4) Bk(mk) and Bk(mk) involve * * the expressions G (m ), G(m ), l D (P)I, and ‘ D (P)‘. If k k Ah Nk Ah Nk r these expressions are computed using h = ho, %- S-%, and 32 s-% * o the bounds Bk(mk) and B(mk) are valid for all mesh constants h 5 ho, and for all mesh points in Sk' Let Gh(PJQ) denote the "Green's Matrix" for the mesh Gh' Gh(PJQ) is uniquely defined by the following properties: 47 (3.3.5) AZGh(P,Q) = - 915493- P e Rh + (3:. Gh(P,Q) = 6(P,Q) P 6 ch + Ch’ * Q E Rh + Ch + C' +'C where 6(PJQ) is the Kronecker delta h h’ [2, p. 315]. The proofs of the next five lemmas follow trivially from [2, pp. 315-318]. Lemma 3.1. (Green's Third Identity) Let W(P) be an arbitrary * * mesh function defined on Rh + Ch +'Cfi + Ch. For any P E Rh + Ch W(P) =h2 2 *Gh(P,Q)[-A:W(Q)] + 2 Gh(P,Q)W(Q). ' QERh+Ch QECh'LCh Lemma 3;2. P z 0 P e + 0* + c' + C P 1 P , * Lemma §;§. z * Gh( ,Q) S c Rh + Ch . QECh * Lemma 3:3. 2 ' Gh(PJQ) = 1 P E Rh + Ch. Qech+ch Lemma 3.5. Let d denote the diameter of the circumscribing circle of G. Then h2 C (P ) s QE- P E + 0* zm*h ’Q 16 Ph h' Q h h From Lemma 3.1 the mesh function e(P) = u(P) - uh(P) has the representation (3.3.6) em =h2 2: *Gh(P,Q)[-A:6(Q)]+ z 'Gh(P.Q)e(Q) QERh+Ch QECh+Ch P E Rh + C:. 48 Since (item) = aim) - efiuhm = gum). and so) = 0. Q 6 Ch, and Gh(PgQ) 2 0, equation (3.3.6) implies (3.3.7) \e(P)| s h2 2 * Gh(P.Q)\A:U(Q)| QeRh+Ch + 2 Gh(P,Q)\e(P)| P e Rh + c2. Qect" Applying Lemma 3.4 to the sums over elements of C6, (3.3.8) \e(P>\ sh?‘ z: *ch(P,Q)\A:u(Q>| +max \e(Q>| . Q6 +ch QECQ * PERh+Ch. An a priori bound for max [3(Q)\ may be obtained from Qécg the results of section 3.2. Let m = min (mk) k=l,...,ko gon is convex (m < l) a positive constant MC may be deter- If the poly- mined such that for all h 5 ho (3.3.9) max |e(Q)I 3. In this case from (3.3.3), for all k hsh, O * 2 * B (m )h h B ( )h (3.3.16) h2 2* Gh(P,Q)\A:u(Q)‘ < h2 ma); ifi— s _Ln.1k__ 3- QEChk Qechk r0 k (3h) "‘k * mk Bk(mk)h * = T ’ P E Rh + Ch' 3 k Bound For Sum Over th Let k be a fixed integer, 1 s k 3 k0. The sum in (3.3.11) over points in th also depends on the angle E— at the vertex V . k Case 1.) Suppose g%-< Eu From equation (3.3.4) and Lemma 3.5, k for all h 5 ho, 2 d2 815%)th2 (3.3.17) h z c (P,Q)\Ahu(Q)| s - max [A u(Q)| < . Q6 h 16 Q6 h 16 th th P c Rh Ch. Case 2.) Suppose TET; = 121-. If Im(z4log z) is a part of DN (r,9), then there are constants 8(4) and b such that for k * h 5 ho, P 6 Rh +Ch 51 2 2 d (3.3.18) h QEER Gh(P,Q)‘Ahu(Q)A s 16 qtgax |Ahu| hk th Bk(4)h2d2 < 16 [1 + b 1og|3h|] Case 3.) Suppose a; >‘E'. The technique used in (1) and (2) above gives a bound for h2 2 G (P,Q)|Ahu(Q)‘ which does not h Qeth tend to zero with decreasing h if mk s 2. A better bound for this sum will be obtained with the aid of the following Lemma which is proved in Appendix IV. Lemma 3.6. Let r denote the distance from a mesh point Q to the vertex Vk. A constant 1 for all h s h , o k may be determined so that 3319 h2 C p 1—<:13 P’ +c* (") 2 h(’Q) 3 h 3 CR1} h' QEth * From Lemma 3.6, for h 5 ho, P E Rh + Ch’ 2 2 2 Bk(mk)h (3.3.20) h z ch(P,Q)|Ahu(Q)| 0) if 2 n. The function gk(x,y) is uniformly bounded on the set 61 From (3.3.23) Laasonen concluded 54 that the degree of convergence of \Z(P)| on any fixed subregion G of G is at least 001) if T-'- 0. CHAPTER IV THE DISCRETIZATION ERROR FOR A DIRICHLET PROBLEM DEFINED ON A REGION BOUNDED BY A PIECEWISE ANALYTIC CURVE In this Chapter the Schwarz Reflection Principle is used to obtain some computable and a priori bounds on the discretiza- tion error for a Dirichlet problem defined on a region bounded by a piecewise analytic arc and having continuous prescirbed poly- nomial boundary data. The results are limited by the lack of accurate a priori information regarding the radius of univalency and the distortion properties of certain analytic functions which arise in applications of the reflection principle. Let R be a bounded, simply connected region in the z = x+iy plane with boundary curve C. Suppose the curve C consists of k0 2 l analytic arcs, T k = l,...,ko. The are k’ F defined by parametric equations 2 = x(t) + iy(t) for t real, ta S t S tb is an analytic are means; (i) x(t) and y(t) may be expanded in Taylor series with real coefficients in a neighborhood of any point t, ta S t S tb’ (ii) for all t, ta S t S tb, [x'(t)]2 + [y'(t)]2 # 0. In this Chapter it is further assumed that the Taylor series expansions for x(t) and y(t) have only a finite number of non-zero terms. If RC > 1, denote the points of intersection of the arcs by V1,V ..,V Denote the angle of intersection at Vk by 2’. k . o , k = l,...,ko. Let f(x,y) be a function which is continuous '51 “‘k 55 56 on C and is a polynomial fk(x,y) with continuous derivatives of all orders on each of the RO arcs. Let u(x,y) denote the continuous solution to the following Dirichlet problem defined on G = R + C; AU(X:Y) 0 (X330 6 R: (4.1) . U(X.y) = f(x,y) (x,y) E C- As in Chapter 3 it is assumed without loss of generality that f(x,y) is zero at all the vertices Vk' In general the function u(x,y) defined by (4.1) cannot be determined exactly and it is necessary to seek an approximate solution to the problem. Let k = 1,2,...,ko, be a set of 6k, positive real numbers, and let G, denote the subset of G which of V k = l,...,k . Let 6k k’ 0 — * Gh = Rh + Ch + Ca + Ch be a square mesh defined on G where the excludes all points of G within * point sets Rh’ C , Cg, and C are defined in a manner analogous h to the correSponding sets for the polygon of Chapter 3. The dis- * tance from any point P E Rh +Ch to a vertex Vk is greater than c = min 8k . Suppose also that the mesh constant h is k=l,...,ko * chosen small enough so that all P 6 Ch are of Type 1. Let uh(P) denote the mesh function on Eh defined by * P 0 P R + 0* Ahuh( ) - E h h, (4.2) uh(P) = f(P) P E Ch, uh(P) = f(vk) = 0 p e Chk’ k = 1,2,...,ko. where Vk is the vertex nearest the point P E Chk' 57 From Theorem 1 (p. 12), if \u(x,y)‘ S M on C, (4.3) [Ah(P)|S Mtg—__— “:4 P E Pm, fa -h 2 3 4 h M 8h * (4.4) |A:U(P)|< 1 + -—§ P 6 Ch’ a/T 4 4 a -h where a is the radius of the largest disk with center at P x and contained in G. The bounds for [Ahu(P)‘ given in (4.3) * and (4.4) cannot be used for any point P E Rh +Ch within h units of the boundary C. However, the Schwarz Reflection Prin- ciple may be used to extend the usefulness of (4.3) and (4.4) * to all points P E Rh + Ch provided h is chosen sufficiently * small. With these known bounds for [Ahu(P)\, a bound for * ‘e(P)‘ - |u(P) - uh(P)\ may be computed for all P E Rh + Ch' An'a priori bound for ‘3(P)‘ may be given in terms of * max [Ahu(P)‘, max [Ahu(P)‘, and a term which arises from the data * h PERh interpolation at points P E Cfi' Schwarz Reflection Principle For An Analytic Arc Let F be one of the analytic arcs which form the boundary C. Suppose F is defined by z = g(t) = x(t) + iy(t), ta S t S tb. Let T = t +’iT be a complex variable in the (t,T)-plane. Let b' By hypothesis g'(ti) f 0, hence there is a positive real number ri such that ti be a real number satisfying ta < ti < t the function gCT) is univalent in the disk D(ri) = {T‘T-ti[ S ri}. The function gCT) maps the disk conformally onto a certain neighborhood of 21 = g(ti), and the diameter of the disk is mapped 58 onto a segment p1 of the arc P. Let p: denote the image + under gCT) of the upper half-disk D (rj) = {TlT-ti‘ S ri, T > O}, + and suppose pi is contained in the set C, (Figure 4.1). T=t+iT CH) ,. Figure 4.1 The following definition is due to H.A. Schwarz; Two points z = gCT) and z* = gCT) are said to be inverse points with respect to the analytic curve P if they are the images of two points T = t + i? and T = t - iT of the disk D(ri). Clearly z and 2* are on opposite sides of the curve F, and it may be shown that the operation of inversion is independent of the choice of the parameter [10, p. 87]. The familiar inversions of a point in a straight line or a circle are Special cases of the above definition. Using the definition of inverse points given above, the Schwarz Reflection Principle for a general analytic arc may be stated as follows: Suppose u(x,y) = u(z) is harmonic in a region R which has the analytic arc F as a part of its boundary. Suppose u(z) = 0 on P. Then the function u(z) is harmonic 59 in a region which contains F as a part of its interior. The values of u(z) in the region are given by u(z) = -u(z*) where z and 2* are inverse points with respect to the curve F. * Bound For [Ahu(P)‘. With ti’ 21, and r1 defined as above, let 6i denote a positive real number such that the image of the disk D(ri) = {T‘T-ti‘ S ri} contains the disk D(éi) = [z||z-zi| S 61}, (Figure 4.1). (The determination of acceptable values for ri and 6i is not discussed in this paper. Positive lower bounds for the numbers ri and 61 may be obtained from some general theorems on conformal mappings (for example, [12, p. 171, 215]), but these values which apply to a general analytic function are much smaller than the values which can be computed for ri and 6i from the known function g(T). A useful criteria for com- puting a radius of univalency ri is given in [11].) The function u(x,y) is harmonic in p:, and u(x,y) = f(x,y) for (x,y) E pi<: F. In the T = t + iT plane, let G(t,¢) = u(x(t),y(t)), and f(t) = f(x(t),y(t)) denote the functions which correspond to u(x,y) and f(x,y) in the z = x + iy plane. Since f(x,y) is a polynomial in x and y, and x(t) and y(t) are polynomials in t, there exists an integer N such that for (x,y) 6 pi, _ _ ~ - _ _ N (4.5) u(x,y) - f(x,y) - f(t) - a0 + 31(t ti) +---+'3N(t ti) where a j = 0,1,...,N are real constants. Define j, N . = + - +. . .+ - . (4 6) I),(t,T) a alRe(T t.) a Re(T t.) 60 The function §i(t,T) = fi(t,T) - Di(t,T) is harmonic in the half disk D+(ri) = {T‘T-ti‘ S ri, T > O}, and is zero on the diameter of the disk. Hence, by the reflection principle, fii(t,T) is harmonic in the whole disk D(ri) = {T“T-ti‘ S ri} with values given by §i(t.'r) = G(m) - 13102.10 «r > 0. (4.7) = O T = O, A c> = -u(t1-T) +'Di(t1-T) T If the function Di(t,T) = Di(t,-T) is now added to the function §i(t,T) in the disk D(ri), the resulting function Wi(t,w) is harmonic in D(ri) with values given by fii(t,¢) = fi(t,¢) T > 0, (4.8) = “f(t) 'r = 0, = -U(t,-T) + 251(C,T) T < 0- Also, in the disk D(ri) the function Wi(t,T) satisfies \wi(t,T)\ s |u(t,-r)| + 2\fii(t,-r)|, N (4.9) s M + 2{|ao| +|a1\ri +...+ \aN\ri}, = M +~M(ri) where M is a bound for |u(x,y)‘ in G. Let Wi(x,y) denote the function which corresponds to Wi(t,m) under the conformal mapping 2 ‘ g(T). The function Wi(X,y) is harmonic in the disk D(Gi) {Zi‘Z-Zi‘ S 81}, and \wi(x,y)| an +M(ri) in D(bi). 61 If P is a grid point near zi, then 4/2 hzmmri» (4.10) A U(P) S P E a 1 h 1 ___8 8 R. a -h 4hm+M(r.)) 2 t (4.11) [a:u(P)[ s 1 [1 + $815] P e 0:, a/a4_h4 a where a is the radius of the largest disk with center at P and contained in G U D(6i). If u(x,y) = 0 for (x,y) E pi, then M(ri) = 0. With formulas of the type (4.3),(4.4), (4.10), and (4.11), a bound for the discrete Laplacian ‘A:u(P)[ may be computed for all points P E Rh + C:° As described in section (3.1), [g(P)‘ satisfies * (4.12) \e

\ s (-A’1>[h2\ofiu«2)\ + (-A>\e|1 P e Rh + cIn * A bound for ‘e(P)‘, P E Rh + Ch may be computed from (4.12) with the bounds for \e(Q)|, Q 6 Ch obtained below. Bounds For Data Interpolation Error Suppose the vertex Vk is at the origin, and let pk > ek be a real number chosen small enough so that the points on Pk or Fk+l are the only points of the boundary C within pk of Vk' Let F be a sector with center at Vk, with radius and k pk with angle 0 S 0 S Bk’ where Bk z'fl- is chosen so that Fk contains all points on the curves Pk, Fk+l Within pk of Vk. Let S = G n Fk’ (Figure 4.2). Let grad f(x,y) denote the k 62 Figure 4.2 gradient of the function f(x,y), and define Mi = max\grad fk(P)" P . v = P e Fk (1 Pk, Mk+1 mgx‘grad f (P)|, P 6 Pk 0 Pk“. Let k+l Mk = max (M£,Mi+1). Suppose |u(x,y)‘ S M for all points in G. Case 1.) Suppose Bk.< n. Define B Mkr cos(9 - ;E M r cos(e - SE) (4.13) Bk(r,0) = B B cos‘ak) pk cos(§k) The function Bk(r,0) is harmonic, and ‘u(r,9)‘ S Bk(r,e) for all points on the boundary of Sk' By the maximum principle, ‘u(r,6)‘ S Bk(r,6) for all points in S In particular, for k. 8 Q 6 Chk’ 414 |6Q)|= ( - ) _ ‘U«Q)‘ S ‘BkGQ)‘ 5 Mk8 +' M e = M e k k cos(§—) pk cos(§—) where ‘Mk is independent of ek' Case 2.) Suppose n S Bk.< 2n. Let a < 1 be a positive real k < flu Define number such that akBk 2 63 B B k k a p cos a 0 - “' 0 cos a 0 "" k (4.15) Bk(r,6) = Mk 15- , k k; 2 ) +-M‘%- .k 5: B 2 1 pk (“k k) k (‘ k k) cos 2 cos 2 The function Bk(r,e) is harmonic, and ‘u(r,9)‘ S Bk(r,e) on the boundary of 8 As above, for Q E C' k' hk’ ".9 a _ a k M k k (4.16) [g(Q)|s + e = e , O’k (“haw “k “kek k Mk k p cos 2 p cos( 2 ) k k : where Mk is independent of ek’ 0 < ak_< 1- Note that if k0 = l, i.e., the boundary C is a single analytic curve, there is no need to interpolate boundary values and the set Cg is empty. Bounds of the type (4.14) or (4.16) may be determined at each of the vertices V . Then a bound for [3(P)‘ may be com- k * puted from (4.12) for all points P E Rh + Ch. A Priori Bound For IESPZI k = l,...,k, It is clear that for each analytic arc Pk’ o a pOSitive integer j, and point sets Tk = {tkl’tk2’°"’tkj}’ RR = {rkl’rk2’°"’rkj}’ Zk = {zkl’zk2’°"’zkj}’ and AR = {6k1’6k2’°'°’6kj} may be determined, along with a real, positive number 5k S so = min 6k , such that for any k=1,000,k -' 0 20 E Pk.“ G the disk ‘z-zo| S 6k is contained in the set x * J G U G, where G = U D(5 ,). Let 5 = min (6 ) , Pk Pk i=1 1“ k=1, . .1.<,ko 771k = max (M(r i))’ 771 = max Mk . It follows from the discussion i=1....,§ k=1,...,k 0 above that if h is chosen to satisfy h < 6, then 64 2 (4.17) max lAhU(P)\ s “[2 h (Mfim) , PERh /58_h8 8h2 3 (4.18) max *|Ahu(P)|< <9—h—(Mim [1 +7] . PCCh MGZ-hh 5 From Lemmas (3.1) to (3.5), for P E Rh + CE, (4.19) \e(P>| shz z Gh,(PQ)\AhU(Q)| +1122 0 h,\e(Q>\. 06cm; (4.20) \e(P)| s g max‘Ahu(Q)‘ + h2 ma3£|A:u(Q)‘ + max if? QGR1n QECh k=1, . . . ,ko (4.21) [G(P)‘< _/2 h2 (M47712 +‘flfiflfilzzl[1 +91%]: max 171k SZR [88—- h8 V544,“ 6. k=1, .,k0 Clearly, the a priori bound (4.21) is a very rough approximation to the bounds for ‘3(P)\ which can be computed from (4.12). The results of this chapter are summarized in the following Theorem 4, Let u(x,y) denote the continuous solution to (4.1), and let uh(P) denote the discrete solution defined by (4.2). A real number 6 may be determined which depends only on the in- version properties of the curve C such that if h < 5, bounds on ‘AEu(P)‘ and ‘Ahu(P)| of the type (4.10) and (4.11) may be determined at all points P E Rh + CZ. A bound for the dis-- cretization error |e(P)| = [u(P) - uh(P)‘ may be computed at all points P E Rh +-C: by (4.12). The computed bound for |e(P)| will satisfy (4.21). APPENDIX APPENDIX * I. Bounds For ‘Ahu(P)|. These are some details of the calculations associated with expressions (2.13) and (2.14), page 11. Expression (2.13) To derive (2.13) from 2 or4 a4 (12012 25 * Z/Z h M 1 2 2 1 2 ““1“” ”hum” ‘ 2 4 44+ 4 44+ 4 22h4 ’ a a -alh -a2h a +d102 2_ 2_ _ h“ let a1 - x, a2 - y, t (a) . Then 04 04 2 azaz l 2 1 2 =‘l_ x y “'1'” 444+444+4224 4 + a-ah a-ah a+aah a l-xt l-yt 1 2 1 2 +111 1+xyt x2 2 2 x (x+y)2(l-xyt)2 (A.1.3) 2 +-“1—2 +—11+xyt = l-x t l-y t (l-xzt)(1-y2t)(l-x2y2t) = 9:322 1+ 3302:; (l-xzyzt) <1-x2c)<1-y2t> From (A.1.l) and (A.1.3), 2 2 * 2/2 h M(ai + a3) ((1% - 0:) hl‘a4 35 (2.13) \A u(P)| s 1 + h [W (a4-o,“h‘*) (.4. 4h“) a -a102h 1 0[2 65 66 Expression (2.14) From (2.11), 6 h8 h10 T * 2 1 h 2 2 2+ 3 3 2+ 4 4 2 | NM —(Cr-a) — 71 * Suppose P 6 Ch is of type 1, i.e., a + a > 1, 3 02 'l' 34 > 1. From (A.2.10) * 2 r: F 2 h (A.2.11) |AO\< 2 [3(7) + h o * 8h 1 (4.2.12) |Ao\ < 3m 1.1.] Suppose P E Rh’ i.e., a = a = a = a = 1. In this case write (A.2.l3) Ah(Im(znlog 2)) = A0 +A +...+A +...+A. From (4.2.10), _ Im zo ’helw 4 l he1w 8 1 heiw 12 A0 ’ 2 ' -1'( z ) - 2 z -.3 Z -... h o o o , r“ 4 4 8 h h lel < 42" (—) [1 + (_) + (L) +...] r r r h o o o , h2 1 (A.2.14) \Ao\ < 4m h 4 r l - '- o ( .) . * The term Ak for k = 1 is 72 r iv iw 2 * _ 2h. n-l e _ l_ 2 2 he (A.2.15) A1 2 n Im(zo ). (r—‘l'd 2(01 0(3) ( z ) h l 3 o ' 3 1W 4 1 3 3 hell 1 he + 3(“1+”3) ( z ) ' 4 z ) + o o + eil 1( 2 02) (heil)2 —— -—a - _— az+u4 2 2 4 z +.l 3 3 heil 3 1, hew 4 3(02+u4) z ) - 41 ( z ) +'°' L o 0 Note that the terms involving log zo have sum equal to zero. * From (A.2.15), if P 6 Ch’ r- ‘1 2n r“ 1h 1 h 2 2 h 3 2 h 4 m <——,—° 3(r) +5(:‘) + <2) (*) + l h o o o 2 3 4 l h 2 h 2 h +[2(r)+3(r)+4(r)+ ] o o o , L. J) * 4 n h 1 (A.2.l6) \A1\ < r3'“ - E_ o r O c For P E Rh’ 1 h o 2 3 z ) 7 z 11 ( z ) -1 4 2 2n rH 3 -'n h o fl_h__ l, l = 3 l (A'2'17) ‘All < h B(r ) 2 h 4 4-n h o r l - (‘-) o l - ('—' r r o o * Bounds for the remaining AR and Ak terms may be obtained * in a similar manner. The behaviour of ‘Ah1m(q?10g 23‘ described in Chapter 2 is evident from these expressions. 73 m n-3 Appendix III. Bound For 2 C T -— — “=3 kn Let k be a positive integer, and m a real number, a V NIr—t It is required to show that (A.3.1) [ck3\ + leéJT + \ckssz +...+\ckn|1‘"'3 < ekm where ‘C \ = l(km)(km-l);..(km-(n-1)L , T S 1 . kn n. 2 For given k,m, there is an integer Nk Such that -lsk . Nk m<< Nk Case 1) km > 3. In this case Nk 2 4. For all n, n = 3,4,...,N \km - (n-1)| = km - (n-l) < km. Hence k ‘C ‘ 1 um>1< mi k3 3! 3! N k km-l ... km- -1 ‘Cm‘=\(m)< )'< mk)‘<_§fll_':. k (N18- (Nk). Since lK—Efili< 1 for n > Nk’ \km-Nk| < C ' -——————-< C \CkCNk'l'l)‘ ‘ kal Nk+1 l kal ‘CkCNk+P)‘ < \cka\ for all P = 1,2,... From the above inequalities, 74 2 n-3 ‘CkB‘ + |ck4\1‘ +|ck511 +...+|ckn\1‘ +.. N (km)3 (km) k 2 Nk‘3 < +...+ [1+T+T +...]T 3! (NR)! Nk-B 1 <-2>...<- 2 >- . h 2r cosA +a h m) -1 -1 “’1 o o 1 ) (A.4.3) ‘zohulhe ‘ r0 [:1 + C1 ( ro)( r o C a h 2 2r cosA +0 h 2 2 o o l +.. _ +... 2! ( r') ( r ) ‘J O 0 ) ‘2 +0 hew‘"1 f-1 2r cosA +a h o l o 1 fl_ 0 o l (A.4.4) ( +u ) = a +u ——'+ C1 r 0’1“1 3 1 3 0’1 to 0 +13; E: 2r0cosAO+a1h)2 + 2! 0[1 2 r r o 0 Similarly, ‘2 -a hewr1 r.1 -2r cosA +0 h o 3 _ o l h o o 3 M» (... ——..— —+c. — ) Q’3 °’1 3 C“1 3 °’3 r6 to C 2 -2r cosA +a h 2 2 o o _- L 3 + 2! a3 2 ( r ) +... r O O . W -1 . + (A 4 6) \z lazhe | _ O .1__+ c L (2r SinAo-l-ozzh) 012(0' ) 0124114 012 1 r0 r0 + C_2_ Q Hi 2rOSinAo+oz2h 2 + 2! 2 r2 r0 0 ) -1 _ ‘z -i(14he ¢| 1 1 h -2r sinA +0! h (A47) (aw) " 3U —+C1—( ° 0“) 4 2 4 0’2 4 014 ro r6 2 . +‘E; a h (-2roSlnAo+a4h) 2 + , A— 2. 2 r r0 0 From (A 4 1) and (A 4 4) - (A 4 7), * 21f)“- h 2 C2 h 2 .1 (A.4.8) Ahgk(Qo) = hz 2C1(-r-) + 27 (T) 4 + o o a) ((0111-03)cosAo+(a2-a4)sinAO) + 2 a3+o/3 a3+u3 L l 3 + 2 4 . (r . a +0! cr +0, + (higher order terms) 0 1 3 2 4 L. J) P * l (AA-9) Ahgkmo) = —3- O h 1+ — - & . 3(1'0) ((01 (1(3)cosAo-|-(or2 0’4)SlnAo) 3 3 2r0’1___:"i_3 “2+“: 3 h z + +4( J L“ 7W3 o’2+0’4 J h -2 L+ 2 - (higher order terms)-(r ) J O O The discussion below would be simplified if a positive * constant 3 could be determined such that Ah(%-\) >s—3 for O ~k . h 1 0 all QO E Ch Wlth :- s 5. However, no general statement, i.e., 9 valid for all Q0 6 Ch with E— s 2%, can be made regarding the sign of the term in (A.4.9) involving (oz1 - (13) cos A0 and (a2 - a4) sin A0. However, it is possible to obtain the desired 78 bound for h2 2 Gh(PgQ) lg by knowing only that QEth r * 88 * (4.4.10) lAhgka < :5 Q 6 ch and l (A.4.1l) Ahgk(Q) > 9—3- Q 6 Pm r . The inequality (A.4.10) may be obtained by bounding the higher order terms in (A.4.9) by a geometric series. From (A.4.8), using ‘§¥l < % for all N, E— s-% , O ‘CB‘h3731C4‘h474 (A.4.12) ‘higher order terms! < 4 —3-!- (r) (-3— + 777(7) (3) o 0 From (A.4.9) and (A.4.12), * (A.4.13) \Ahgk(qo)\ < 33 [1 + 2 + % + 2(42)] < g3- . r O 0 If P E Rh’ the ambiguous term described above is zero since a1 = 03 and 02 = a4. As shown below, if the terms in- volving C1 to C27 in (A.4.4) to (A.4.7) are added exactly, and the remaining terms are bounded by a geometric series, then (A.4.ll) follows. (Proof of A.4.ll) For n = 1,2,..., define n _ n _ n . n (A.4.l4) S - (ZrOCOSAO+h) + ( 2rocosAo+h) + (2rOSInAO+h) + (-2r sinA +h)“, O O +..] 79 Then from (A.4.4) - (A.4.7), with ai = l, i = 1,2,3,4, -l e. k k 1 _r0 3 h s _ 3 2k-1 (A.4.15) Ahq) - 2 2: ,(r) ‘1: where ck- <-e)<-2)...(-1—12 ) h k=1 0 r . 0 5c k k 2 _E. E_. .§_ = .E_ (A.4.16) E ,(r) k (r) [4c1+4c2] k—l o r o O + L4E2CMC#C( 4A + '4A) (r0). 2 3 3 4 cos 0 Sln o ] +L6zc+zc+3c 4A.+1“A (ro)°(3 3 4 3 5‘COS o S“ 0)) 8 10 h 1 2 h 1 +(ro)'(6c4+3cs) +(ro)'(30c). Using the bound — S cos4Ao + sinaAo S l, 5 C k k 2 4 6 S (4.4.17) 2 45 1‘— —2 L + 11.. . -19. + L(§_5_) k! r k r r 8 r 2 k=1 0 r0 0 O O + 1.1.8 _105 11 + 1.)“) _162 (r ) 16 6 ) (r 5 ) 0 o For E—'s‘% o 5 Ck h kSk h 2 h 2 (A.4.18) 1551;;(T) 72(r) [1 - .58] =(;—). (.42). - r0 0 O The terms from k = 6 to k = 27 will be considered in pairs; the terms for k greater than 27 will be bounded by a geometric series. k k k+ k+ (A419) 2735 LkS—= 2 g(Lr. §—+Ck+1 (L)1——S 1 ' ' k=6 k! r k k! ro rk (k+1)! r rk+l ° ro k=6,8,10,.. ,26 0 0 0 II M O‘r—fi CD .75le CA ”I h 1‘ sk (2mg 31“” h.) k 2(k+l) k+1 r r 0 a 2 o o 80 .k_ k! “1th (1) ' (k-j)!j! ’ k __ k k k-l (A.4.20) S - ((ZrOcosAo) + (1) (2rocosAo) h +. . .+ k k-l k (k- 1) (2rocosAO)h + h ) k k k-l + ((-2rocc sAo) + (1‘) (-2rocosAo) h +. . .+ k k-l k (k_1)(-2rocosAo)h + h ) . k k . k-l + ((2rOSinAO) + (1) (2ros inAo) h +. . .+ k . k-l k‘) (k_1)(2roSinAo)h + h . k k . k-l + ((-2roslnAo) + (1)( ZrOSinAo) h +. . .+ k . k-l k) (k_1)(-2rOSinAo)h +h . Denote coskAo + sinkAO by TR. For k even, 0 < TR 3 l, and k 2 k (A.4.21) S—k = 2. 2ka + 2 (‘32ka 2 k 2(L) +...+ 4F.) r r r o o o k+1 3 k+l S __ . k+l k _ k+12k-2Tk-2 h h (A.4.22) WET 2 (1)2k'1‘ (2) + 2(3 )2 (17) +...+ 46) r O O O o For all j, j = 1,2,3,..., (kT1)/( k) = k—,+1- . Hence, for h 1 J j- -1 J all :- S '5' , o C k k k+l 1 (A413) 3 “15' (L) g: ’ 2311:?) :R+1 (3m 2 9 0 O .. 7(3) 25114—4515) ””2053 <‘;(1>(-‘———§‘f;'3) + 2(3glk (1 ' 226:1)9) 81 For k = 6 and 8, all terms on the right in (A.4.23) are positive. For k = 10,12,14,...,26 the only negative terms is _ (2k+1) = 17-2k the coefficient (1 m 18 . th __52k+11,sk+1 h. k 2(k+1) k+1 (r _) 6170 0 r0 Hence, c k _k .h. (A.4.24) z k, (r ) k=6,8,10,?..,2 C k k 17- 2k 2 -2 E ET'(:O)R2 ZKT ‘ 18 l k=lO,12,...,26 Ck 1 For k210, ‘k-fl <3. C k 2 k k h 17- 82k h (A.4.25) 2 k, 2 (to) |1| < ( 04) ( ) k=10,12,...,26 ° Hence (A 4 26) 2 35' h— 2 §E - $3k+llsk+l 5— > -( 08) h— 2 ° ' k! to k 2(k+1) k+1 r0) ' (r ) k=6,8,...,26 r0 r0 0 C h k Sk The terms k! ;-' —k is positive for k even. Hence 0 r O (A 4 27) Ex, E”k g5 > _- 1929‘ h. 29 S29 + |c31| h 31 E3: ° ' 3 k! r k 29: (r ) "26 31: r l 31 k=27 o r o r o r 0 O 0 SR 1 k 7 k h 1 -E1‘< (2 +-§) = (g) . For -— 3'3, the term on the to 2 I.0 right in (A.4.27) is less than 4° L) -1- -2(Z-) (7)27 ——1—-J< .05 b— r 3 3 9 r o 1 Z 2 o (9).J Combining (A.4.18), (A.4.26), and (A.4.27) l. 1 1 4 (A.4.28) Ah(ro) > 3 [.42 - .08 - .05] > 3 . r r O O 2 .4 82 From Lemma 3.1, the mesh function gk(P) defined on has the representation 9: + + ' + P E Rh Ch Ch Ch (12.4.29) gk(P) r122 Gh E- -Ahgk(Q)] + 2&2 I: 6,, (P Q>E- -Ahgk(0>fl QEthh J'i‘k QEthh +h2 2 ch (P.Q)E-A g. (0)] + 2[h2 2 MG (P.Q)E- -Ahg (0)1] Ah k QERhIh J'J‘k QéCzhj + r122 *hG (P,Q)E- -Ahgk (0)] + h2 2 ”G (P Q)E- 'Ahgk ((2)1 qeca QEChkh * + 2 G (P.Q)g (Q) . P E + C Qecfiflhh k Rh h Consider the signs of the terms in (A.4.29). gk(P) > 0 * for all P 6 Rh + Ch +’C' + C . The first three terms on the h h right in (A.4.29) are negative since Ch(PJQ) 2 O for all P41 and -AhngQ)‘< - ll%-. The last term on the right is positive, and the sign of the remaining terms is not known. Hence, in the * worst case, for P E Rh + Ch’ (A.4.30) 4lh3 2 ch (P,Q)— 2 (1122 *Gh(P,Q) §§)+h2 2 Ch (P Q)— QEthh r3 .x 1(’3590) V2050) Figure A.5.l < 1 The problem solved by discrete methods is the problem for u10(x,y) indicated in Figure A.S.2. In this case the non-zero data was approximated by the polynomial f10(x) = (l-CO) - mfi+mfi- -mfl 2 l n 4 l 10! Where (30:1 ’ (2) 2. + (2) 4! n 101 (2 ) 0! so that f10(%) = f10(-%) = 0. =0 u =0 R V4.._u__._lQ_ ____WV3 V4 2. fl.--V3 ulO—O Au10=0 lu10=0 UR: . AuR=O uR=O i j ; u10(x,0)=f10(x)i } uR(x.0>=R ! ‘ L 4‘ 2" P v 2 —\7 ‘7 W1 T 2 1 2 Figure A.S.2 85 2 4 10 Note that f(x) = (1 - Lg-i-(L-l-‘(EQ— -...- LY%%—_ - Co) + (Co + R11(x,0)) = f10(x) + R(x) where R11(x,0) is the Taylor Series Remainder; ‘R11(x,0)‘ s (£912 i%T-< .0000005. The function uR(x,y) with boundary data indicated in Figure A.5.2 satisfies |uR(x,y)\ s |R11(x,0) + col < .000001 for all points on the boundary of the square; from the maximum principle for harmonic functions, ‘u(P) - u10(P)‘ < .000001 for all points P in the square. This error, introduced by replacing the function u(x,y) by u10(x,y), is negligible compared to the computed bounds for the discretization error for ulo(x,y) obtained below. Let G denote a square mesh of width h on the square h with grid lines parallel to the coordinate axes and such that V ,V 1 2, and the origin are grid points. Since all angles are %, it is not necessary to interpolate boundary data at the vertices. For all P E Gh’ all neighbors of P are a distance * h from P. With Rh’ Ch’ and C defined as in Chapter 1, let h u10 (P) denote the solution to the linear system h S l ‘ O * (A. . ) Ahu10h(P) — P E Rh + Ch’ u10h(P) = f10(P) P 6 Ch. Let 310(P) = u10 (P) - ulo(P). Using the techniques of Chapter 2, h * bounds for |Ahu10(P)\ were computed for all points P E Rh + Ch For mesh points P whose distance rO from V (or V l l 2 by symmetry) satisfies ro s 3, the best bounds for ‘Ahu10(P)] were obtained usin the results of Theorem 2 with E-= E, = l, 8 m 2 p i.e., expressions of the form 86 *5 (11.5.2) \Ahu10(1>)| < M2(G(2)) h2(M+M(1)) + ‘AhD10(P)\. In (A.5.2) D10(x,y) is the function which removes the non-zero data from the side Vin; M(l) is a bound for 1D10(X,y)‘ in the sector of radius 1, center at V1; and M is a bound for \u10(x,y)‘ in the square. The value of, G(2) depends on the value of re, (cf. Table 2, p. 27). At the grid point (-.4, .l) the bound ‘Ahu10(P)‘ < .116 was computed from (A.5.2). Let P with coordinate (xo,yo) be a mesh point in the square. It follows from the reflection principle that the function u10(x,y) is harmonic in the disk (x-xo)2 +-(y-yo)2 < y:. Hence, from Theorem 1, for all mesh points P with yO > h, 2 (A.5.3) ‘Ahu10(P)l s 51.2.1114 , (cf. p. 12), xy:_h8 The expression (A.5.3) gave the best bounds for \Ahu10(P)‘ for points in the upper-half of the square. For example, at (0,.7), the bound ‘Ahu10(P)‘ < .015 was computed from (A.5.3). For mesh points P "near" the x-axis the best bounds for ‘Ahu10(P)| were obtained from Theorem 1 and the data removal technique for a straight line, i.e., from expressions of the form M2 >h2 / 8_h8 a (11.5.4) \Ahu10(P)‘ s , (cf. p. 18), At the mesh point (0,.3) the bound 1Ahu10(P)‘ < .185 was computed from (A.5.4). * For most mesh points P E Rh + Ch several of the expressions for ‘Ahu10(P)‘ were applicable - in this case the minimum of the 87 computed bounds was used. An upper bound for the expression ‘u(P) - 1110 (P)‘ is h given by (A.5.5) |u(P) - u10h(P)| s |u(P) - u10(P)| + \u10(P) - u10h(P)l < .000001 + 1310(P)| Some typical computed values for u10 (P): 1610(P)‘, and u(P) h are given in Table A.5.1 for h = .025. Round-off error in this problem is negligible compared to the discretization error [3, p. 319]. 88 Table A.5.l. Computed Bounds for Problem A.5.l, h = .025. - i (x,y) u10h(P) u(P) \u(P) - u10h(P)| 1610(P)| - - __Ji. _.2 (-.4,.1) .225373 .225337 .000036 .00217 (—.2,.1) .590034 .539941 .000093 .00552 ( 0,.1) .729322 .729207 .000115 .00718 (-.4,.2) .164132 .164081 .000051 .00359 (-.2,.2) .429705 .429571 .000134 .00908 ( 0,.2) .531145 .530979 .000166 .01078 (-.4,.3) .119208 .119153 .000055 .00431 (-.2,.3) .312092 .311947 .000145 .01021 ( 0,.3) : .385767 .385588 .000179 .01151 (-.4,.4) é .086134 .086082 .000052 .00395 (-.2,.4) E .225503 .225366 .000137 .00892 ( 0,.4) E .278738 .278568 .000169 .01032 ' (-.3,.5) 1 .117213 .117127 .000086 .00538 (—.1,.5) 1 .189655 .189515 .000140 .00800 I (-.3,.6) % .082241 .082171 .000070 .00396 ‘ (-.l,.6) .133070 .132955 .000115 .00601 (-.3,.7) .055445 .055291 .000054 .00276 (-.1,.7) .089712 .089626 .000086 .00426 (-.3,.8) .034160 .034124 .000036 .00174 (-.1,.3) .055272 .055215 .000057 .00271 (-.3,.9) .016271 .016253 .000018 .00085 -.....—_. - BIBLIOGRAPHY 10. ll. 12. BIBLIOGRAPHY Greenspan, D., Introductory Numerical Analysis of Elliptic Boundary Value Problems, Harper and Row, N.Y., 1965. Bramble, J.H., and Hubbard, B.E., "0n the Formulation of Finite Difference Analogues of the Dirichlet Problem for Poisson's Equation", Numerische Mathematik, 4(1962), 313-327. Forsythe, G., and Wasow, W., Finite Difference Methods for Partial Differential Equations, Wiley, N.Y., 1960. Bramble, J.H., and Hubbard, B.E., "A Priori Bounds on the Discretization Error in the Numerical Solution of the Dirichlet Problem", Contributions £2 Differential Eguations, 2(1963), 229-253. Hills, N. Notices gf_the American Math. Soc. vol. 13, no. 15, (1966), 612. Laasonen, P., "On the Behaviour of the Solution of the Dirichlet Problem at Analytic Corners", Ann. Acad. Sci. Fenn..A;I., 241(1957), 1-13. __ , "0n the Degree of Convergence of Discrete Approximations for the Solutions for the Dirichlet Problem", Ann. Acad. Sci. Fenn. A.I., 246(1957), l-19. , "0n the Truncation Error of Discrete Approxima- tions to the Solutions of Dirichlet Problems in a Domain With Corners", g, Assoc. Comput. Mach., vol. 5, no. 1, (1958), 32-38. Wasow, W., "Asymptotic Development of the Solution of Dirichlet's Problem at Analytic Corners", Duke Math. 1., 24(1957), 47-56. Caratheodary, C., Conformal Representation, Cambridge University Press, No. 28, 1963. Kaplan, W., "Close to Convex Schlicht Functions", Michigan Math. Journal, 1(1952), 169-184. Nehari, Z., Conformal Mapping, McGraw-Hill, New York, 1952. 89 90 13. Phillips, H.B., and Wiener, N., "Nets and the Dirichlet Problem", J, Math. Phys., 2(1923), 105-124. 14. Courant, R., Friedrichs, K.0., and Lewy, H., "fiber die Partiellen Differenzengleichungen des Mathematischen Physik", Math. Ann., 100(1928), 32-74. 15. Gerschgorin, S., "Fehlerabschfitzung ffir das Differenzen- verfahren zur Lfisung Partieller Differentialgleichungen, g, Angew. Math. Mech., 10(1930), 373-382.