Asymptotic expansions of the error of
چکیده
We present and analyse a family of fully discrete spline Galerkin methods for the solution of boundary integral equations with logarithmic kernels. Following the analysis of Galerkin methods of a previous paper, we show the existence of asymptotic expansions of the error. In this work we deal with the numerical solution of the integral equation Z 1 0 log jx( ) x(s)j2g(s)ds = u0(s); 8s; where u0 is a given 1 periodic function and x is a smooth 1 periodic parametrization of the boundary of a simply connected bounded open set in the plane. We will suppose that the logarithmic capacity of that curve is not equal to one, in order to avoid problems of existence and uniqueness of solution. In Crouzeix and Sayas (1995) the existence of asymptotic expansions of the error of Galerkin methods with smoothest splines was shown. However, the practical impossibility of evaluating exactly in most cases the coe cients and the right-hand of the linear system obtained by the Galerkin method demands the use of approximate numerical integration in the process. We will prove that the formulae transforming the Galerkin method into the so-called Galerkin collocation method (see Hsiao, Kopp and Wendland (1980 and 1985)) lead to a formulation of the problem where there exists still an asymptotic expansion of the error in powers of the discretization parameter h. In fact, as in Crouzeix and Sayas (1995), we will be able to prove expansions like kg h Qhg M X k=f(m)hkQhfkk ChM+1; (1) where Qh denotes an interpolation operator onto the discrete space, f(m) is an integer depending on the degree of the splines used as test-trial functions for the numerical method, and g h is the numerical approximation of g. Throughout what follows, C will denote a certain constant, independent of the parameter h, possibly di erent in each occurrence. The symbol O(hk) denotes a function of the parameter h such that divided by hk remains bounded as h goes to zero. 1 The Galerkin method Let us consider a simply connected bounded open set inR2, whose boundary is given by an in nitely often di erentiable 1 periodic non-singular parametrization x. We suppose that the logarithmic capacity of (see Sloan (1992) for example) is not equal to 1. We consider the periodic Sobolev spaces H t equipped with the norm kukt := 0@jû(0)j2 +X k 6=0 jkj2tjû(k)j21A 1 2 ; where û(k) := Z 1 0 u(s)e 2k {sds: We denote by ( : ; : )0 both the inner product in H0, that is, in L2(0; 1), and its extension as the duality bracket between H t and H t for any t. Let V : H t ! H t+1 be the single layer operator, de ned as the continuous extension of V g(s) := Z 1 0 log jx(s) x( )j2g( )d 1 and let b denote the bounded symmetric bilinear form on H 12 b(g; r) := (V g; r)0: To solve approximately the variational problem (P )8><>: nd g 2 H 12 such that b(g; r) = (u0; r)0; 8r 2 H 1 2 ; we consider the space Vh of the periodic smoothest splines of degree m with nodes on the points si+ 12 := i+ 12 h; i 2 Z; where h := 1=N and N is a positive integer. Denoting Ii := (si 12 ; si+ 1 2 ), if m = 0, Vh := fqh 2 H0 : qhjIi 2 P0;8ig; and if m 1 Vh := fqh 2 Cm 1 : qhjIi 2 Pm;8ig; being Ck the space of real 1 periodic k times di erentiable functions and Pk the set of the polynomials of degree lesser than or equal to k. The Galerkin method is the sequence of nite-dimensional variational problems (Ph)8><>: nd gh 2 Vh such that b(gh; rh) = (u0; rh)0; 8rh 2 Vh; uniquely solvable for h small enough. With the application of numerical integration in coe cients and right-hand side of a linear system equivalent to (Ph), this will be transformed into a new problem (P h )8><>: nd g h 2 Vh such that bh(g h; rh) = (u0; rh)h; 8rh 2 Vh; where bh( : ; : ) and ( : ; : )h are respective approximations of the bilinear form b and of the inner product in L2(0; 1), which will be speci ed in the next sections. Let si := hi for all i 2 Z. The interpolation operator Qh maps u 2 C0 to the unique Qhu 2 Vh such that Qhu(zi) = u(zi); 8i; where zi := ( si; if m is even, si 1 2 ; if m is odd. Consequently Qh interpolates on the midpoints of the grid for m even and on the nodes for m odd. 2 Once we have an expansion of bh(Qhg g h; rh) similar to that of Proposition 1 in Crouzeix and Sayas (1995), the uniform inf-sup condition will be used to prove expansions like (1). Since for all rh 2 Vh bh(Qhg g h; rh) = = b(Qhg gh; rh) + (u0; rh)0 (u0; rh)h + bh(Qhg; rh) b(Qhg; rh) (2) the study of such an expansion can be divided into three parts, the rst of which is that of the Galerkin method. The second part includes the in uence of the approximations of the right-hand side of the linear system, and the third one that of the approximations of the coe cients of the matrix (which require the evaluation of double integrals). Some additional notations about the spline space are called for. Let 0 be the characteristic function of the interval ( 1 2; 1 2) and m := 0 m 1 for all m 1. Then for m 1, m is a non-negative function with bounded support [ m+1 2 ; m+1 2 ], unit integral, of class Cm 1, piecewise polynomial of degree m, and so on. These and other properties can be reviewed in Aubin (1972) or Hsiao, Kopp and Wendland (1980). For i 2 Z, and h small enough, we de ne i to be the 1 periodic extension of the function m s zi h ; de ned in a unit length interval around zi. Obviously f igNi=1 is a basis of Vh. Let us nally consider a quadrature formula ZR m(t)f(t)dt ' Lm(f) := l X i= l cif(xi); where for all i xi = x i; ci = c i; and l is a certain positive or null integer. We will demand that the formula Lm be of degree m+ 1 at least, that means, that it be exact for polynomials of degree lesser than or equal to m + 1. However, the above demanded symmetry of nodes and coe cients implies that Lm will always have odd degree. For instance, for m = 0, we can take a single point formula (l = 0) L0(f) := f(0); which is exact in P1. For m = 1, the following three-point formula L1(f) := 1 12f( 1) + 5 6f(0) + 1 12f(1) is exact not only in P2, but also in P3. A three-point formula will also be su cient for m = 2, L2(f) := 1 8f( 1) + 3 4f(0) + 18f(1): 3 2 Approximation of the right-hand side The problem arising in the computation of the right-hand side is the evaluation of the integrals Z 1 0 u0(s) j(s)ds = h ZR u0(zj + th) m(t)dt; which, applying the quadrature formula, can be substituted by hLm(u0(zj + h : )) = h l X i= l ciu0(zj + hxi): This is equivalent to de ning an approximate linear form in Vh, (u0; rh)h := h N X j=1 rh;jLm(u0(zj + h : )); (3) where rh(s) = N X j=1 rh;j j(s): Because of their periodicity, the functions i are equal modulo N , so we can associate to rh 2 Vh the sequence frh;jgj2Z, with rh;j = rh;k for j k, modulo N . Consequently the summation in (3) can be done on any set of N consecutive indices. We denote [x] the largest integer which is lesser than or equal to a real number x. Proposition 1 For all u 2 C1, there exists a sequence of functions fungn 1 C1 such that for all rh 2 Vh, (u; rh)0 (u; rh)h = M X n=[m+1 2 ]+1 h2n(un; rh)0 +O(h2M+2)krhk 1 2 : Proof: For simplicity of notations we only prove the case m even. Since ZR m(t)dt = 1; we can write (u; rh)0 (u; rh)h = h N Xj=1 rh;j (zj; h) (4) where (s; h) := ZR m(t)0@u(s+ ht) l X i= l ciu(s+ hxi)1A dt: Let us remark that is in nitely often di erentiable in both variables, is 1 periodic in s, and satis es that for all s (s; h) = (s; h) (5) 4 because of the symmetry of nodes and coe cients. Moreover, since Lm is exact in Pm+1, we have that for k m+ 1 ZR m(t)0@tk l X j= l cixki1A dt = 0; from where @k @hk (x; 0) = 0; (6) for all k m+ 1. Denoting fk(s) := 1 (2k)! @2k @h2k (s; 0); by Taylor expansions in (4) and applying (5) and (6), it follows that (u; rh)0 (u; rh)h = M X k=m2 +1 h2k[fk; rh]h +O(h2M+2)krhk 1 2 ; (7) where [ ; rh]h := h N X j=1 rh;j (zj): (8) In Crouzeix and Sayas (1995) (see formula (29)) the existence of a sequence of real numbers fDjg, depending only on m (even or odd), such that for all 2 C1 and for all M [ ; rh]h = ( ; rh)0 + M Xl=1 h2lDl( (2l); rh)0 +O(h2M+2)krhk 1 2 (9) is proven. With uk(s) := k m2 1 X n=0 Dnf (2n) k n (s); the statement of the proposition is simply a consequence of (7) and (9). 3 Approximation of the bilinear form. Let us denote Jj := zj m+ 1 2 h; zj + m+ 1 2 h : Then the restriction to Vh of the bilinear form b can be written b(gh; rh) = N X i;j=1 bi;jgh;irh;j where bi;j := ZJi ZJj log jx( ) x(s)j2 i( ) j(s)d ds 5 and fgh;ig; frh;ig are the respective coe cients of gh and rh in the basis f ig of Vh. Our aim is to calculate approximately the coe cients bi;j. For simplicity of notations we will consider the usual sequence of coe cients frh;jgj2Z for each rh 2 Vh. From the periodicity of the functions i, bi;j can be de ned for all i; j 2 Z, modulo N . Let us de ne the index set I(N; j) := 8><>: fj p; j p + 1; : : : ; j + pg; if N = 2p + 1, fj p; j p + 1; : : : ; j + p 1g; if N = 2p: Thus, the bilinear form can also be written b(gh; rh) = N Xj=10@ X i2I(N;j)bi;jgh;i1A rh;j : With this choice of indices, ji jj N=2. Denoting F ( ; s) := log jx( ) x(s)j2 ( s)2 ; we can write, for all j and for all i 2 I(N; j), bi;j = ZJi ZJj F ( ; s) i( ) j(s)d ds+ h2(log h2 +Bi j;m); where Bk;m := ZR ZR log (k + (t u))2 m(t) m(u)dtdu; (10) or also, bi;j = h2 ZR ZR F (zi + ht; zj + hu) m(t) m(u)dtdu+ h2(log h2 +Bi j;m): (11) Let us notice that for all k, B k;m = Bk;m. Since F ( ; s) 2 C1(D), where D := f( ; s) : j sj 2 3g; we can easily apply the bidimensional quadrature formula proceeding from Lm and de ne an approximation of bi;j, for all j and for all i 2 I(N; j) i;j = h2 l X k;k0= l ckck0F (zi + hxk; zj + hxk0) + h2(log h2 +Bi j;m) (12) where the Bk;m are calculated exactly as shown in Hsiao, Kopp and Wendland (1980). That leads us to de ne an approximate bilinear form in Vh by bh(gh; rh) = N Xj=10@ X i2I(N;j) i;jgh;i1A rh;j: (13) 6 It is important to see that bh( : ; : ) is symmetric. If N is odd, then j 2 I(N; i) if and only if i 2 I(N; j), so the symmetry is obvious. If N = 2p, a more careful look is needed. If we de ne j+p;j by (12), then it is easy to see that j p;j = j+p;j, because of the properties of the coe cients and nodes of the quadrature formula Lm. Hence, we can de ne k := 8><>: 1 2 ; if k = p or k = p, 1; if jkj 0 such that for all h h1, jb(gh; rh) bh(gh; rh)j Chm+1kghk 1 2 krhk 1 2 ; for all gh; rh 2 Vh. Proof: Subtracting (11) from (12), we get that bh(gh; rh) b(gh; rh) = h2 N X j=10@ X i2I(N;j) gh;irh;j (zi; zj;h)1A (15) 7 where ( ; s;h) := ZR ZR(F ( +ht; s+hu) l X k;k0= lckck0F ( +hxk; s+hxk0) m(t) m(u)dtdu: The function belongs to C1(D [ h0; h0]) and satis es ( + 1; s+ 1;h) = ( ; s;h); ( ; s; h) = ( ; s;h); s+ 1 2 ; s;h = s 12 ; s;h and for all k m+ 1 @k @hk ( ; s; 0) 0: This last property is a consequence of the degree of the quadrature formula in the sense given in its de nition. Hence, jb(gh; rh) bh(gh; rh)j Chm+4 N X j=1 jrh;jj N X j=1 jgh;jj from where the result follows readily since h N X j=1 jrh;jj Ch 1 2krhk 1 2 for all rh 2 Vh. From the previous lemma, it follows that (at least for h small enough) the family of approximate bilinear forms bh satis es uniformly the Brezzi-Babu ska condition. Therefore (P h ) is uniquely solvable and its solution is uniformly bounded in norm H 12 by the data. The convergence of this method was proven in Hsiao, Kopp and Wendland (1980). We will deal with this aspect in the following section, before giving the asymptotic expansions of the error. Proposition 3 There exists ffngn m2 +1 C1 such that for h = 1=(2p + 1), and for all rh 2 Vh, b(Qhg; rh) bh(Qhg; rh) = M X n=m2 +1h2n(fn; rh)0 +O(h2M+2)krhk 12 : Proof: Let be the function de ned in the proof of the previous lemma and let k( ; s) := 1 (2k)! @2k @h2k ( ; s; 0); 8 for k m2 + 1, which satisfy k( + 1; s + 1) = k( ; s): Then, if faigi2ZN are the coe cients of Qhg in the usual basis of Vh, it follows bh(Qhg; rh) b(Qhg; rh) = h2 N Xj=10@ j+p X i=j p airh;j (zi; zj;h)1A = = M X k=m2 +1 h2k+20@ N Xj=1 j+p X i=j p airh;j k(zi; zj)1A +O(h2M+ 3 2 )krhk 1 2 ; because of the properties of . From properties of the operator Qh (see Sayas (1994)), we know also that there exists a sequence of numbers Cl independent of h and g such that ai = M Xl=0Clh2lg(2l)(zi) +O(h2M+2); uniformly in i 2 ZN . Thus, denoting for k m2 + 1 k( ; s) := k m2 1 Xl=0 Clg(2l)( ) k l( ; s); we have that bh(Qhg; rh) b(Qhg; rh) = = M X k=m2 +1h2k+20@ N Xj=1 j+p X i=j p k(zi; zj)rh;j1A +O(h2M+ 3 2 )krhk 1 2 ; (16) for all rh 2 Vh. Set now k;0(s) := Z s+ 12 s 12 k( ; s)d = Z 12 12 k( + s; s)d : Since for all ( ; s) 2 D, k( + 1; s + 1) = k( ; s); it follows that the functions k;0 are 1 periodic, and given that the functions k are smooth in D, so are k, from where k;0 2 C1. Remark now that h j+p X i=j p k(zi; zj) is the composite midpoint rule applied to compute k;0(zj). By the asymptotic expansion of the error in this formula h j+p X i=j p k(zi; zj) = k;0(zj) + M 0 X n=1 h2n k;n(zj) +O(h2M 0+2); (17) 9 where k;n(s) := 2n Z 1 2 12 @2n @ 2n k( + s; s)d and 2n are some real coe cients. Note that k;n are also smooth and 1-periodic. Applying (17) with M 0 = M k to (16), and denoting 'n(s) := n m2 1 X k=0 n;n k(s); it follows that for all rh 2 Vh, bh(Qhg; rh) b(Qhg; rh) = M X n=m2 +1h2n['n; rh]h +O(h2M+2)krhk 1 2 ; where the discrete product [ : ; : ]h is de ned by (8). Finally, (9) gives the result. Proposition 4 There exists f ngn m2 +1 C1 such that for h = 1=2p, and for all rh 2 Vh, b(Qhg; rh) bh(Qhg; rh) = M X n=m2 +1 h2n( n; rh)0 +O(h2M+2)krhk 1 2 : Proof: Let k be the same functions as those de ned in the proof of Proposition 3, and let us remark that since k s+ 1 2 ; s = k s 12 ; s ; then k s+ 12 ; s = k s 12 ; s for all s and consequently k(zj p; zj) = k(zj+p; zj); for all j. Hence h j+p 1 X i=j p k(zi; zj) is the composite trapezoidal rule applied to compute k;0(zj) (being this function the same as in the previous proposition). The remainder of the proof of Proposition 3 is valid, changing the coe cients appearing in the asymptotic expansion of the error for the midpoint rule by those for the trapezoidal rule. Remark. For m odd, we have that for all rh; gh 2 Vh jb(gh; rh) bh(gh; rh)j Chm+2kghk 1 2 krhk 1 2 ; and consequently existence and uniqueness of solution follows. Propositions 3 and 4 also hold here, changing only the lower limit of the sum, which now is n = (m + 3)=2. The expansions for N even and odd are also di erent. 10 4 Asymptotic expansion of the error In Hsiao, Kopp and Wendland (1980), the following convergence result is proven for m = 0; 1 and 2, although it can be extended to arbitrary m. Theorem (convergence). Let u0 2 Hm+2, g be the unique solution of (P ) and g h the unique solution of (P h ). Thenkg h gks Chm+1 skgkm+1 (18) for all s 2 [ 1 2;m], if m 1, and for s = 0 if m = 0. Let us simply remark that the original result in that paper is stated for a slightly modi ed problem. Propositions 1, 3, 4 and the decomposition of bh(Qhg g h) given in (2) prove the following result. Proposition 5 For all positive integer M , there exist functions ek, as smooth as desired, such that bh(Qhg g h; rh) = M X k=2[m+2 2 ]hk(ek; rh)0 +O(hM+1)krhk 1 2 (19) for all rh 2 Vh, if h = 1=(2p + 1). A similar expansion also holds when h = 1=2p. Let us consider the Sobolev space Wm;1 := fu 2 L1 : u(j) 2 L1; 1 j mg; with the essential supremum norm k : km;1. Theorem 6 For all positive integer M , there exist functions fk, as smooth as desired, such that kg h Qhg M X k=2[m+2 2 ]hkQhfkkm;1 ChM+1 (20) if h = 1=(2p + 1). A similar expansion also holds when h = 1=2p. Proof: The proof of Theorem 2 of Crouzeix and Sayas (1995) is valid here, applying now Proposition 5 and the properties of the family of approximate bilinear forms bh. Since the convergence theorems give the same orders of convergence in the same `strong' norms for the Galerkin and the Galerkin collocation methods, the other results of this section do not di er essentially from those given for the Galerkin method. Corollary 7 Let h be the set of the midpoints of the grid. Then for all 0 l m, there exists a sequence of constants fFk;l(z) : z 2 hg, independent of h = 1=(2p + 1), such that sup z2 h j(g h)(l)(z) g(l)(z) M X k=2[m l+2 2 ]hkFk;l(z)j = O(hM+1) 11 for all positive integer M . A similar expansion also holds for all h = 1=(2p). Let h be the set of the nodes of the grid. Then for all 0 l m 1, there exists a sequence of constants fF k;l(z) : z 2 hg, independent of h = 1=(2p + 1), such that sup z2 h j(g h)(l)(z) g(l)(z) M X k=2[m l+2 2 ]hkF k;l(z)j = O(hM+1) for all positive integer M . A similar expansion also holds for all h = 1=(2p). Finally, what is said about superconvergence in L1 norm for the Galerkin method is true here too. Likewise it is possible to build an interpolate solution of higher order when m is even as was pointed out there. 5 Numerical results We have tested the fully discrete Galerkin collocation method with piecewise constant functions (that is, m = 0), with the midpoint rule as numerical quadrature (see examples given at the beginning), on some problems whose exact solution is known. We show here simply an example where it can be seen that kQhg g hk1 Ch2; (21) and that Richardson extrapolation can be applied as a means of accelerating the con- vergence. Remark that from the last corollary it follows that we have `true' asymptotic expansions (with all coe cients independent of the numerical parameters) when restrict- ing ourselves to a nite set of points, all belonging as midpoints to all meshes under consideration. The number of points where extrapolation is to be applied restricts the size of the di erent grids where the solutions are calculated. As the simplest test we show here results on a single point. Let be the ellipse (x=R)2 + y2 = 1. Then with data u0(x; y) := x3 3xy2; (x; y) 2 ; the solution of the variational problem is the product of q(x; y) = 3(R + 1)2 2R2 x q x2 R4 + y2 (R 2(R + 1)y2); where x = (x; y) is substituted by the usual parametrization of , by the norm of the parametrization jx0(s)j2. We choose R to be equal to 2: Figure 1 shows the midpoint convergence of the solution computed for N = 10; 20; 30; : : : ; 400. We have plotted there the decimal logarithm of the number of intervals against log10(EN ) where EN stands for the maximum error in the midpoints, that is, (21). Applying linear regression to that set of points we obtain an estimate of the order of the method. The slope of the regression line is in this case 1.962. 12 Figure 1: convergence We next show a table of the errors obtained by extrapolation when calculating the solution in (R; 0). The solution has been computed for N = 20; 40; 80; 160 and 320; that is in a Romberg set, making midpoints coincide. Richardson extrapolation has been applied to solutions in consecutive grids. The table shows the errors: columns show di erent stages of the extrapolation process, rows show errors in correspondence with the nest grid used for each computation. N = 20 1:35 N = 40 3:53 10 1 2:0 10 2 N = 80 9:04 10 2 2:65 10 3 1:73 10 4 N = 160 2:28 10 2 3:66 10 4 3:93 10 5 3:04 10 5 N = 320 5:75 10 3 4:86 10 5 3:31 10 6 9:09 10 7 4:44 10 8 Numerical results show how Richardson extrapolation gives a remarkable increase of precision, in addition to faster convergence. Other numerical tests for this method can be seen in Sayas (1994). We wish nally to indicate that when solving some di erential problems by a bound- ary element method, the approximate solution calculated by some postprocessing of the numerical solution of the BEM can inherit the existence of asymptotic expansions and thus the possibility of applying directly extrapolation. References Aubin, J.P. (1972): Approximation of elliptic boundary-value problems, Wiley-Inter- science. 13 Crouzeix, M. and F.J.Sayas (1995): Asymptotic expansions of the error of Boundary Element Methods I: spline Galerkin methods, submitted for publication. Hsiao, G.C, P.Kopp and W.L.Wendland (1980): A Galerkin collocation Method for some Integral Equations of the First Kind, Computing 25, 89-130. Hsiao, G.C, P.Kopp and W.L.Wendland (1985): Some applications of a Galerkin-Col- location Method for Boundary Integral Equations of the First Kind, Math. Meth. in the Appl. Sci. 6, 280-325. Le Roux, M.-N. (1974): R esolution num erique du probl eme du potentiel dans le plan par une m ethode variationelle d' el ements nis, Th ese de troisi eme cycle, Universit e de Rennes, France. Sayas, F.J. (1994): Asymptotic expansion of the error of some boundary element meth- ods, Tesis doctoral, Universidad de Zaragoza, Spain. Sloan, I.H. (1992): Error analysis of boundary integral equations, in Acta Numerica, Cambridge University Press. 14