A CONVERGENT BOUNDARY INTEGRAL METHOD FOR THREE ...

MATHEMATICS OF COMPUTATION Volume 00, Number 0, Xxxx XXXX, Pages 000{000 S 0025-5718(XX)0000-0

A CONVERGENT BOUNDARY INTEGRAL METHOD FOR THREE-DIMENSIONAL WATER WAVES J. THOMAS BEALE Abstract. We design a boundary integral method for time-dependent, threedimensional, doubly periodic water waves and prove that it converges with O(h3 ) accuracy, without restriction on amplitude. The moving surface is represented by grid points which are transportedaccordingto a computedvelocity. An integral equation arising from potential theory is solved for the normal velocity. A new method is developed for the integration of singular integrals, in which the Green's function is regularized and an ecient local correction to the trapezoidal rule is computed. The sums replacing the singular integrals are treated as discrete versions of pseudodierential operators and are shown to have mapping properties like the exact operators. The scheme is designed so that the error is governed by evolution equations which mimic the structure of the original problem, and in this way stability can be assured. The wavelike character of the exact equations of motion depends on the positivity of the operator which assigns to a function on the surface the normal derivative of its harmonic extension; similarly, the stability of the scheme depends on maintaining this property for the discrete operator. With n grid points, the scheme can be implemented with essentially O(n) operations per time step.

1. Introduction In this paper we design a semidiscrete numerical method of boundary integral type for computing time-dependent, doubly periodic, three-dimensional water waves. We prove that this method converges as long as the actual motion remains smooth, without restriction on its amplitude. While the full equations are dicult to work with analytically, boundary integral methods seem naturally suited to water wave motion. The surface is represented in Lagrangian coordinates, and grid points are moved according to a computed velocity. Because the ow is assumed irrotational, the motion is determined by quantities on the surface alone. However, the normal velocity must be found in terms of singular boundary integrals, using potential theory. The evolution has the character of nonlinear, nonlocal wave motion without dissipation. For this reason it seems that the stability of the numerical method depends critically on its design. In the present work we develop a new, ecient approach for calculating the singular integrals, using the limited information computed on the moving surface. We analyze discrete versions of integral operators such as single and double layer potentials by adapting the viewpoint Received by the editor June 28, 1999. 1991 Mathematics Subject Classi cation. Primary 65M12, 76B15; Secondary 65D30. Key words and phrases. water waves, boundary integral methods, integral operators, quadrature of singular integrals. The author was supported in part by NSF Grant #DMS-9870091. 1

c 1997 American Mathematical Society

2

J. THOMAS BEALE

of pseudodierential operators, as well as potential theory, to the discrete setting. We nd the most important terms in the evolution equations for the errors, and ensure by design of the scheme that the errors grow at a bounded rate. Boundary integral methods have been widely used in ocean engineering and applied mathematics, starting with the work of Longuet-Higgins and Cokelet [21] and Vinje and Brevig [33]. Two-dimensional applications have been extensive; see the recent surveys [27],[32] or the references in [5, 6]. Computations in 3-D include [4, 27, 28, 9, 32, 17]. In [6], T. Hou, J. Lowengrub, and the author proved the convergence of a version of the method in 2-D. Calculations illustrate the advantage of this method in capturing small scale features. It seems that convergence in 3-D is more dicult, largely because quadrature for singular integrals on surfaces is more involved than on curves, and because the choice of quadrature aects the stability of the time-dependent solution. One approach is presented in [19]. We assume as usual that the ow of the water is incompressible, inviscid, and irrotational. For simplicity we assume the uid is in nitely deep and the motion is periodic in both horizontal directions, with period 2. (The formulation can easily be modi ed to account for a horizontal bottom; for more general lower boundaries, see [3].) The motion of the uid is governed by Euler's equations and boundary conditions at the surface. Since the ow is irrotational and incompressible, the velocity is the gradient of a potential  which is harmonic below the surface: (1.1)
 = 0 : At the surface the uid pressure matches the atmospheric pressure, which we assume constant. (We neglect surface tension.) Bernoulli's equation, together with the boundary condition for the pressure, provides us with an evolution equation for  at the surface. With the position x and velocity potential  at the surface treated as functions of Lagrangian coordinates  = (1; 2) and time t, the evolution equations are xt = v ; t = 12 jvj2 ? gx3 (1.2) where v is the uid velocity, g is the acceleration of gravity, x3 is the vertical component of x, and  is adjusted by a spatial constant at each t. The rst equation says that a particle on the surface with xed  moves with the uid velocity. The velocity v = r is determined by  on the surface and the condition (1.1), with the requirement that r is square integrable as x3 ! ?1. In the boundary integral approach, the main eort in solving these equations is in nding the velocity v at the surface using potential theory. There is some choice in how this is done. The tangential part of the velocity can be found from dierentiating  along the surface. For the normal component we use a Fredholm integral equation of the second kind which determines the normal velocity n = @=@n directly from  on the surface, Z  (x ? y) 1 n(y) dS(y) (1.3) 2 n(x) + @G@n(x) Z ?
= n(x)
rG (x ? y)  (n(y)  rT(y) dS(y) Here G is the Green's function, periodic in the horizontal variables, obtained from a sum of images of the standard Green's function G(x) = ?1=4jxj with
G = ; n is the unit normal pointing out of the uid. We assume  is periodic and has square

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

3

integrable gradient below the surface. The rst integral is absolutely convergent, but the second is of principal type. The right side is the normal derivative of the double layer potential due to . This integral equation has been used in electromagnetics. It can be derived by writing  below the surface in terms of boundary integrals in  and n, and then taking the normal derivative at the surface (e.g., see [10], x3.9 and Thm. 2.23, or [20]). Its use for three-dimensional water waves has been suggested in [28]. An alternative approach is to write the velocity potential in terms of a dipole source on the boundary; this was done in [3, 4, 17, 5, 6]. To design a speci c scheme, we have choices to make in writing the integrals, nding derivatives, and especially in the numerical integration of the singular integrals. Analytical considerations in the present work and in [6] indicate that these choices can determine whether or not the scheme is stable. There are well developed methods for quadrature of singular integrals in several dimensions which depend on special coordinates near the singularity; see Lyness [24] for a survey. In the present context the surface and integrand are known only from computed values at the Lagrangian grid points. Here we develop a general approach in which we replace the Green's function G by a smooth version Gh regularized on the scale of the grid size h. We discretize the new integral by the trapezoidal rule. Then, using an asymptotic expansion of the quadrature error like that found by Lyness [23] and Goodman et al. [12] for unregularized singular integrals, we identify the largest error term. Because of the regularization, we can compute this rst error and correct for it; it is expressed through the Poisson Summation Formula in terms of the Fourier transform G^ h of Gh . In our case the correction improves an O(h) error to O(h3 ). (We have avoided using extrapolation as in [23] for the sake of a positivity requirement explained below.) As in [6] we are guided in designing the scheme by the linearization of the exact equations about an arbitrary solution. By preserving the same structure at the discrete level we can ensure that the numerical scheme is stable. Without such considerations we could have a mismatch of terms which destroys the well-posedness of the linearization and therefore the numerical stability. Brie y, this linearization reduces to an equation for a variable u related to the disturbance in the velocity potential on the surface, (1.4)

utt + cu  0 ;

where we have dropped less important terms, and  is the principal part of the Dirichlet-to-Neumann operator on the current surface, i.e., the operator that determines n from . It is important that  and c are positive, as discussed further below. The equation describes nonlocal wave motion; in the special case of small disturbances at equilibrium, it reduces to the usual equation for linearized water waves. The linearization about an arbitrary motion was found in the boundary integral setting in 2-D in [5], and in 3-D with orthogonal coordinates in [18]. It was derived earlier for bounded domains in a dierent way in [2]. A primary concern is that the discrete version of the operator  should be positive, or almost so. Violation of this requirement in high wavenumbers would lead to rapid growth of disturbances, as can be seen from (1.4). For the exact operator , the positivity is related to the fact that the Green's function has transform G^ < 0. This property may not be preserved by the discrete analogue. In the present approach, we choose Gh so that G^ h < 0. The discrete transform is a lattice sum of

4

J. THOMAS BEALE

values of G^ h , which converges because Gh is regularized. In this way we can ensure the necessary positivity. We now describe the method in more detail. We rst rewrite the integrals in (1.3) to reduce the order of singularity, taking advantage of the two identities Z Z  (x ? y) @G (1.5) rG (x ? y)  n(y) dS(y) = 0 @n(y) dS(y) = 0 ; where the integral is over one period of the surface. (Cf. [4], x3 or [10], Thm. 2.1. The use of these identities to reduce the singularity was suggested in [4]). On the right side of (1.3) we use vector identities to convert the integrand, ?
?
(1.6) rG  n  rT = rG
rT n ? (rG
n) rT = (rG  n)  rT ? (rG
n) rT : Then using (1.5) we can write the integral on the right in (1.3) as (1.7)

Z

?




(rG  n(y))  rT(y) ? rT(x) dS(y)

?

Z

?




rG
n(y) rT(y) ? rT(x) dS(y) :

The rst integral is absolutely convergent because of the subtraction. We can replace rG by rT G , the part of the gradient tangential at y and write (cf. [10], p. 35), with y = y(0 ) and X1 ; X2 the tangent vectors Xj = @y=@0j , rG (x ? y)  N(y) = (D01 G )X2 ? (D02 G )X1 (1.8) Here N = X1  X2 ; we assume N points outward, so that n = N=jN j. For the integral on the left in (1.3), we resolve rG into parts normal and tangential at y and use (1.5) again, so that the integral becomes (1.9) Z Z rG
n(y) [n(y)
n(x)n(y) ? n(x)] dS(y) + rT G
n(x)n(y) dS(y) : The second integrand is absolutely integrable since rT G ? n(x) at y = x. We work with grid functions on a grid of size h representing equally spaced values of the Lagrangian variable, j = jh, where j = (j1 ; j2) 2 Z2. We assume h = 2=N with N even. For the discrete values in a fundamental period we take (1.10) I = fj 2 Z2 : ?N=2  j  N=2 ? 1;  = 1; 2g : A grid function f has a discrete Fourier transform which we denote by f, a function of k 2 Z2 with period 2=h in k1; k2,  = (2)?2 X f(jh)e?ikjh h2 ; f(jh) = X f(k)e  ikjh : f(k) (1.11) j 2I

k2I

For properties of this transform see e.g. [15, 31]. We measure grid functions on Ih in the discrete L2h norm, given by X X  j2 : (1.12) jf j2L2h = jf(jh)j2 h2 = (2)2 jf(k) j 2I

k2I

Of the exact solution, we assume that the surface function x(; t) and the velocity potential (; t) are smooth up to some time T, and x(1; 2; t) ? (1; 2; 0) and (1; 2; t) are 2-periodic in each of 1; 2. (If the velocity is periodic but  is

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

5

not, we can choose a frame of reference so that  becomes periodic. The scheme can easily be adapted to the more general case.) We also assume that the Jacobian @x=@ is nondegenerate and x(; t) 6= x(0; t) when  6= 0. This implies that (1.13) jx(; t) ? x(0 ; t)j  c j ? 0j ; 0  t  T : Next we discuss the discrete derivative operator. It can be a dierence operator, but to be general we discuss its Fourier representation. For a function f of j 2 Ih, we assume the discrete partial derivative Dh; in  ,  = 1; 2, has the form  ; ( + 2) = () (1.14) (Dh; f)(k) = h?1(k h)f(k) where  : R ! C. We also assume (1.15)  is C 1 on [?; ] ; (?) = () ; C1 j j  j()j  C2j j on [?; ] ; (1.16) () = i + O( 4) as  ! 0 ; (?) =  () where  is the complex conjugate. The latter says that Dh; is accurate to O(h3) and preserves real-valued functions. Because of (1.15), Dh; has a product rule: for f 2 L2h and u a smooth periodic function, (1.17) Dh; (uf) = uDh; f + B(f) where B is an operator depending on u, bounded on L2h uniformly in h (see [6], Lemma 1). This will be useful later. The requirement () 6= 0 for  6= 0 rules out, e.g., the symmetric fourth order dierence operator. A simple example meeting these conditions is the third order dierence operator coming from Lagrange interpolation
? (1.18) Dh; = (6h)?1 ?2S?1 ? 3I + 6S ? S2 where S is the shift operator, S1 f(j1 h; j2h) = f((j1 +1)h; j2h)). Alternatively, we could de ne  directly, e.g. by setting  = i + (1 ? ) for j j  , where () is some even, real-valued function with   1 near  = 0 and   0 near  = . We  , given by shall refer to the adjoint operator Dh;  f)(k) = h?1  (k h)f(k)  : (1.19) (Dh;  is found by replacing S with S?1 . For f smooth, Of course, for (1.18), Dh; 3  = ?@f=@ + O(h3 ). Dh; f = @f=@ + O(h ) while Dh; Various geometric quantities must be computed from the grid points on the surface. From x(j ); j 2 I, we can use Dh; to compute approximate tangent vectors Xj = (Dh; x)j ,  = 1; 2. (We apply Dh; to x() ? (; 0).) From these we compute the metric tensor g = X
X and the inverse g , as functions of j . We then nd the normal vectors N = X1  X2 and n = N=jN j, and also the dual tangent vectors at j , X =  g X ; X
X =  : (1.20) Given the surface x and the velocity potential  at the surface, at one time for the grid points j , j 2 I, we can advance to the next time using (1.2) once we have computed the velocity v at each j . To do this we write the computed velocity as v = rTh  + wn (1.21)

6

J. THOMAS BEALE

where w is an approximation to n. We can compute the tangential gradient as (e.g., cf. [10], p. 33)

rTh  =

(1.22)

X

(Dh; )X :

 =1;2

We will obtain w from a discrete version of (1.3), rewriting the two integrals as in (1.7)-(1.9). We use the regularized Green's function Gh ; for the normal derivatives we use the analytic gradient of Gh , while for rT Gh (x ? x(0 )), the part of the gradient tangential at x(0), we use the discrete 0-derivative,

rTh Gh (x ? x(0 )) = ?

(1.23)

X

 =1;2

(Dh; Gh (x ? x(0 ))) X :

This expression is used to match (1.22); both appear in the discrete operator  mentioned above, and this symmetry gives it the proper structure. The regularized Green's function Gh is constructed in x2 from a free space version in the form Gh (x) = ?(4jxj)?1s(jxj=h) for x 2 R3 , where s is a function chosen with certain conditions (2.2),(2.3),(2.9) which ensure that Gh is smooth and diers from G itself as an integral operator by O(h3). As above, we also assume that G^ h (k) < 0. A speci c choice (2.19) for s is described in x2, based on the Gaussian function, for which the requirements are satis ed and necessary quantities can be computed explicitly. Since the regularization has length scale h, Gh diers appreciably from G only within distance O(h) of the origin. If, for example, x() = (1; 2; 0) and we approximate the single layer potential due to a function f by applying the trapezoidal rule with Gh in place of G, we have (with a limiting value at ` = j) X

(1.24)

`2Z2

G((jh ? `h; 0))s(j ? `)f(`h) h2

and the s factor amounts to a set of quadrature weights. The quadrature error is O(h); however, it can be improved to O(h3) by a correction proportional to f(jh), as explained in x3. Such corrections appear in the discrete integrals below. The discrete integral equation is obtained by substituting (1.7),(1.9) into (1.3) for the two integrals, replacing G with Gh , applying the trapezoidal rule, and adding the needed quadrature correction. It has the form 1 (1.25) 2 wj + (Kh w)j = fj = Fj
nj ; (1.26) X X (Kh w)j = rGh
N` [(n`
nj )w` ? wj ]h2 + rTh Gh
nj w` jN` jh2 + hj wj (1.27) Fj =

`

X

`

[Dh Gh ; X]`  (rTh` ? rThj )h2 + hj ?

X

rGh
N` (rTh` ? rThj )h2

` ` Here rTh  is found from (1.22) and rTh Gh from (1.23); rGh etc. are evaluated at xj ? x`; and in view of (1.8) we have introduced (1.28) [Dh Gh ; X]` = Dh;1 Gh (xj ? x`)X2` ? Dh;2 Gh (xj ? x`)X1`

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

7

where Dh is with respect to `. The correction terms are given by (1.29) j = j [Ph (Dh; 1 X1 + Dh; 2 X2 )]j
Nj ?
(1.30) j = j X1j  [PhDh; 2 (rTh )]j ? X2j  [Ph Dh; 1 (rTh )]j where   q X (1.31) ? 2 g (j )n n

j = 2(det g (j ))1=2 06=n2Z2

and ? is a function related to Gh ; see (2.14),(2.24). The sum in (1.31) is in nite, but we can choose Gh so that ? decreases rapidly and only a few terms are needed (cf. (2.24)). The derivation of these corrections is explained in x3. A linear operator Ph is inserted in (1.29),(1.30) to smooth out the second derivatives. We assume it has the form p  (Ph f)(k) = (k h)f(k) (1.32) with  chosen so that (1.33) () = 1 + O(j j4) as  ! 0 ; () = 0 for j j > c0 ; c0 > 0 : As part of the convergence argument, we prove that the integral equation can be solved by simple iteration, i.e., 1 w(n+1) + Kh w(n) = f ; (1.34) w(0) = f : 2 The following theorem gives the convergence result for this scheme.

Main Theorem. Suppose an initial state is prescribed, such that the exact equations of motion (1.1),(1.2) have a smooth solution for 0  t  T meeting the conditions above (1.13). Let the numerical scheme be as presented, with the discrete derivative Dh; chosen according to (1.14)-(1.16) and the operator Ph as in (1.32),(1.33). Let Gh be a regularized Green's function satisfying (2.1)-(2.3),(2.9) and G^ h < 0. Then for h suciently small, the scheme has a solution with the prescribed initial state for 0  t  T , and in particular the discrete integral equation (1.25) can be solved. For each time t  T , the surface computed at the grid points j diers from the exact surface by O(h3 ) in L2h . The computed velocity potential at the surface has error O(h3 ) in L2h , and the computed velocity at the surface has error O(h2 ) in L2h .

This scheme is semidiscrete, i.e., discrete in space but not time. We expect the time integration to be done using an ODE solver such as the Runge-Kutta or Adams-Bashforth methods. Higher order versions of this scheme are possible but would require more complicated correction terms. We do not deal in detail here with the ecient computation of the discrete integrals; it is important to note, however, that with n = N 2 points to be tracked, the computation can be done in essentially O(n) operations per time step. The four integrals can be put in matrixvector form, where the matrix comes from either Gh or its gradient. There are two ways to proceed, based on the fast multipole method of Greengard and Rokhlin [14] or Ewald summation. Both have been used for doubly periodic water wave computations [4, 17]. We can choose the regularization of the Green's function, e.g. (2.19), so that it is signi cant only at distance O(h) from the singularity. Then fast summation can be used except for O(1) terms per integral, or O(n) terms in total. An alternative is to adapt the fast version of Ewald summation of Strain [29].

8

J. THOMAS BEALE

As usual the theorem asserts the convergence of the scheme as long as a smooth solution exists. Recently Sijue Wu has proved that a smooth solution of the initial value problem for water waves exists for a time depending on the initial state, in a three-dimensional uid of in nite extent, tending to equilibrium at in nity [35]. (Earlier results were two-dimensional or assumed analytic data.) It is reasonable to expect that a result similar to Wu's holds for the present case of doubly periodic motion. Wu shows in her case that the pressure gradient at the surface is bounded away from zero, (1.35)

c = ?rp
n  c0 > 0

where n is the normal outward from the uid region. Such a result was proved for a bounded domain in [2]. The same condition holds in the present case; this can be seen as in x4 of [35]. The essential reason is that ?p is subharmonic and therefore obeys the strict maximum principle. The coecient (1.35) appears in (1.4). The positivity is important for the well-posedness of the exact equations and consequently also for the numerical stability of the scheme under consideration, as seen in x6. It means that the motion is well-posed even after a wave overturns. (Of course the model ceases to apply once the wave crashes.) Further discussion of the signi cance of (1.35) can be found in [5]. In analyzing the stability of the scheme we need to understand the mapping properties of the discrete integral operators applied to error terms. It is natural to view them as discrete versions of pseudodierential operators; the important parts can be regarded as convolutions and estimated in the Fourier transform. This point of view is especially helpful in establishing the positivity of the operator  using the assumption G^ h < 0. In x4 we develop some basic properties of discrete pseudodierential operators under mild conditions, including a version of Garding's inequality. In [25] general results for such operators were derived assuming a high wavenumber cut-o. Here we avoid this assumption and derive more limited results. In x5 we show that various discrete operators related to Gh are bounded, or gain derivatives, as linear operators on L2h . These are discrete versions of standard properties of single and double layer potentials. In the stability estimates of x6 we use these mapping properties to identify the principal error terms and simplify them, without treating less important terms in detail. The singularity subtraction helps to prevent spurious terms from appearing in the stability analysis. In fact, a scheme in 2-D like that of [6] can be shown to converge if singularity subtraction is used, without smoothing the points xj inside the integrals as was done in [6]. We brie y outline the contents of the remaining sections. In x2 we describe the regularized Green's function Gh and derive needed formulas, construct the periodic version Gh , and prove that the error from replacing G with Gh in a single or double layer potential is O(h3 ). In x3 we analyze the quadrature error for the trapezoidal rule applied to singular integrals with regularization and derive a formula to remove the largest error by local correction. We apply these results to the discrete integrals of (1.26),(1.27). In x4 we prove basic properties of discrete pseudodierential operators, and in x5 the boundedness properties of discrete integral operators. The convergence of the scheme is proved in x6: stability estimates are obtained for the computed velocity, using the results of x4 and x5; a simpli ed equation is found for the growth of the error; and estimates are found for the rate of growth. Some arguments needed for the proof are deferred to x7.

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

9

The operator  is essentially a convolution with rT G, following a derivative; see Lemma 6.8 and (6.31) for the discrete version. The positivity of the discrete  depends on the choice of quadrature for the singular integral, as well as the choice of the two derivatives. It appears dicult to maintain this condition in discretizing the integral without regularizing G; the sign condition on G^ h and the matching of the two derivatives ensure the positivity here. In the stability argument of x6, the positivity of  is needed because it enters the energy estimates; see (6.55). We also need the discrete  to dominate the discrete rst derivative, i.e. the derivative can be written as a bounded operator times  (see (6.49-50)). It is for this latter reason that we require the derivative in (1.14) to have nonzero symbol  for k 6= 0; if  decays in the high wavenumbers, then the symbol of  decays faster, since  has two derivatives. In the 2-D convergence argument of [6], the analogue of rT G was the Hilbert transform H; the discrete H given by alternate quadrature had symbol with modulus one, and the two needed properties of  = HD followed directly. 2. The regularized Green's function In this section we describe the class of regularized Green's functions and nd the Fourier transform of the restriction to a plane in terms of the 3-D transform. We need to compute this transform explicitly for the correction terms in the integral equation (1.25). We nd a speci c Gh , based on the Gaussian function, which has all the needed properties. We construct the periodic Green's function as a sum of images, adjusted by constants. Finally, we prove (Theorems 2.2, 2.3) that the error in a single or double layer potential due to the regularization is O(h3 ). The regularized Green's function Gh will have the form (2.1) Gh (x) = ?(4r)?1 s(r=h) ; r = jxj with s chosen so that s(r) ! 1 rapidly as r ! 1 and Gh is smooth. A convenient assumption for the decay of s ? 1 is (2.2) jDrk (s(r) ? 1) j  Cr?4?k ; r  1 : For Gh to be smooth, we need s(jxj)=jxj to be smooth as a function of x, and thus as a function of jxj2; i.e., we assume (2.3) s(r)=r is a smooth function of r2; 0  r < 1 : Then for h > 0, Gh (x) = h?1G1(x=h). For later use we note (cf. [1]) ( C jxj?jkj?1 ; jxj  h ; jDxk Gh (x)j  Ch (2.4) ?jkj?1 ; jxj  h : Letting =
G1, we nd (2.5) (x) =
G1(x) = ?(4r)?1 s00 (r) ; r = jxj and correspondingly (2.6)
Gh (x) = h (x)  h?3 (jxj=h) ; which approximates the delta function as h ! 0. For accuracy we need a further condition on s. It is well-known (e.g. [16]) that the error from replacing G by Gh in a integral is determined by moments of s ? 1

10

J. THOMAS BEALE

or . Here we are interested in integrals on a surface, and so we consider, for a homogeneous polynomial pn (x) of degree n, the integral Z 1 Z (s(jxj=h) ? 1) jxj?1p (x) dx (G (x) ? G(x))p (x) dx = ? (2.7) h n n 4 R2 R2 which becomes in polar coordinates (2.8)

Cn

Z

0

1

(s(r=h) ? 1) rn dr = Cnhn+1

Z

1

0

(s() ? 1)n d :

The largest error is O(h) with n = 0. For this reason we require of s, (2.9)

Z

0

1

(s(r) ? 1) dr = 0 :

For odd n, the error terms above are zero since the integrand is odd. Thus with condition (2.9) we expect the error from the regularization to be O(h3 ). This is veri ed for single layer potentials in Thm. 2.2. Given an arbitrary s(1) satisfying (2.2),(2.3), we can easily produce a related function satisfying (2.9), as well as (2.2),(2.3), by setting (2.10) s(3) (r) = s(1) (r) + r @r s(1) (r) : Condition (2.9) can be veri ed for s(3) by an integration by parts. The corresponding relationship for the 's is (3) = 4 (1) + r@r (1) . (A similar strategy for satisfying moment conditions was used in [8].) When Gh has the form (2.1), we have a similar expression for rGh , x s~(r=h) ; s~(r) = s(r) ? r @ s(r) : (2.11) rGh (x) = 4r r 3 Here s~ has the same properties (2.2),(2.3) as s. If s satis es the moment condition (2.9), then the same holds for s~; this can be seen by an integration by parts. We need explicit formulas for G^ h , the Fourier transform of Gh in R3, and for the transform of Gh on planes, since we work with integrals on surfaces. We will write the Fourier transform of a function f on Rd as Z Z ^ ikx dk : ^ = (2)?d=2 f(x)e?ikx dx ; f(k) = (2)?d=2 f(k)e (2.12) f(k) The hypotheses imply that G^ h (k) is smooth for k 6= 0 and decays rapidly for large k. More particularly, it will be important to require G^ h (k) < 0. Since ^ = ^(k), it is equivalent to require ^(k) > 0. To nd the transform of ?jkj2G(k) the restriction to a plane, we rst identify (x1; x2) 2 R2 with (x1 ; x2; 0) 2 R3 . Then the Fourier transform of the function on R2 given by x ! Gh (x; 0) is Z 1 ? 1 = 2 G1 (
; 0)^(k1 ; k2) = (2) (2.13) G^ 1(k1 ; k2; k3) dk3 : ?1

Since G^ 1 is radial, this depends only on j(k1; k2)j, and we write G1(
; 0)^(k1; k2) = ?(j(k1; k2)j), with ? de ned by Z 1 ? 1 = 2 ?() = (2) (2.14) G^ 1(; 0; `) d` : ?1

More generally we will use the transform of G1  J : R2 ! R, given by  ! G1 (J), where J : R2 ! R3 is a 3  2 matrix. In our applications J will be the Jacobian

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

11

matrix of a coordinate mapping. Since G1 is radial, G1(J) is a function of jJj2 = jJ J
j2 = jBj2, where B = (J  J)1=2. Thus G1  J() = G1(B; 0), and (2.15) (G1  J)^(k) = (det B)?1 G1 (
; 0)^(B ?1 k) ; k 2 R2 ? f0g : We can now combine this with (2.13), using jB ?1 (k)j = j(J  )?1kj, to obtain ?
(2.16) (G1  J)^(k) = (det J  J)?1=2 ? j(J  )?1kj ; k 2 R2 ? f0g : If the two columns of J are vectors X1 ; X2 2 R3, then (J  J) = X
X = g , the usual metric induced in the coordinate plane. Also j(J  )?1 kj2 = g k k , summed over ; , where (g ) = (g )?1 . Thus the formula above becomes p  (G1  J)^(k) = (det g )?1=2 ? g k k ; k 2 R2 ? f0g : (2.17) Of course (G1  J)^(k) < 0 since G^ 1 < 0; this will be important later. We now describe a speci c choice of Gh which meets all our requirements. We start with the error function Z r 2 (2.18) erf(r) = p2 e?s ds : 0

It ful lls conditions (2.2),(2.3). Using (2.10), we de ne s(r) = erf(r) + r erf 0(r) or (2.19) s(r) = erf(r) + 2?1=2re?r2 so that s satis es the moment condition (2.9). For this s we nd (2.20) s~(r) = erf(r) + 2?1=2(2r3 ? r)e?r2 (2.21) (r) = ?3=2(4 ? 2r2)e?r2 : As for the transform, we have   G^ 1 (k) = ?(2)?3=2 12 + jkj?2 e?jkj2=4 ; k 2 R3 : (2.22) Then from above,  Z 1  1 1 ? 2 + 2 + `2 e?(2 +`2 )=4 d` : (2.23) ?() = ?(2) 2 ?1 Using [13], formula 3.466, we have   p ?() = ?(2)?2 e?2 =4 +  erfc(=2) (2.24) where erfc = 1 ? erf. We could decide to use s(jxj=ah) instead of s((jxj=h) in the de nition of Gh , where a > 0 is some constant. The eect would be to replace Gh by Gah . Then (G1  J) would be replaced by (Ga  J), and (2.17) by p (2.25) (Ga  J)^(k) = a(det g )?1=2 ?(a) ;  = g k k : The condition G^ h < 0 is evident in (2.22) for the speci c choice of Gh given by (2.19). The condition would hold for a variety of other choices as well. Suppose we rst choose a smooth delta function as in (2.6) by taking (1) to be a smooth, rapidly decreasing, radial function, with integral one, such that ^(1) (k) > 0 and ^(1) decreases with jkj. If we de ne s(1) by (2.5), s(3) by (2.10), and Gh by (2.1) using s(3) , it follows that the needed conditions hold, including G^ h < 0.

12

J. THOMAS BEALE

We next discuss the periodic versions of the Green's function and its regularization. Because of the slow decay at in nity, we write the sum with a re ection and with a constant subtracted from each term. We de ne  X0  1 (2.26) G(x ? 2(n; 0)) + G(x + 2(n; 0)) + 21jnj G (x) = 2 n2Z2

where the prime on the sum means that 1=jnj is omitted at n = 0. The gradient is X rG (x) = 21 (2.27) (rG(x ? 2(n; 0)) + rG(x + 2(n; 0))) n2Z2

and we de ne Gh , rGh similarly. Lemma 2.1. The sums for Gh (x) and rGh (x) converge uniformly on bounded sets for xed h, and those for G (x) and rG (x) converge in L1 on bounded sets. G (x) and Gh (x) are periodic and satisfy X X
Gh (x) = (2.28)
G (x) = (x + 2(n; 0)) ; h (x + 2(n; 0)) : n

Proof. For xed x and large y, we have 

n



1 1 ? y
x + O(jyj?3 (2.29) as y ! 1 : G(x + y) = ? 4 jyj jyj3 If we add G(x  y), the second term above is canceled, and so with y = 2(n; 0), we see that the nth term in (2.26) is O(jnj?3). Thus the sum for G (x) converges pointwise away from the singularity, and since the singularity is integrable, the sum also converges in L1 on a bounded set. For the regularized version Gh , we note that (2.30) jGh (x  2(n; 0)) ? G(x  2(n; 0))j = O(jnj?5) as n ! 1 because of (2.2). Thus the sum for Gh converges similarly to G , but Gh is bounded for xed h, so that the convergence is uniform on bounded sets. It is correct to apply
termwise in the distributional sense, so that both equations (2.28) hold as distributions. However, since
Gh is smooth, the second case is true in the classical sense as well. Similar considerations apply to the gradients. It remains to verify that G and Gh are actually periodic. Speci cally, let x be a point which is not a periodic image of 0; we verify that G (x + (2; 0; 0)) = G (x). We will use the notation r(n) = jnj?1 if n 6= 0, r(0) = 0, and also n+ = (n1 + 1; n2; 0), n? = (n1 ? 1; n2; 0). Let BM = f(n1; n2) 2 R2 : jnj j  M; j = 1; 2g. Then (2.31) X? + j?1 + jx ? 2n?j?1 ? 2r(n)
: ?8G (x + (2; 0; 0)) = Mlim j x + 2n !1 BM

We will show below that we can replace 2r(n) above by r(n+ ) + r(n? ). If we do this, the sum over BM can be rewritten as X? X?

jx + 2n+ j?1 ? r(n+ ) + (2.32) jx ? 2n? j?1 ? r(n? ) : BM

BM + In the rst sum we can shift the n back to n, ranging over BM , except for O(M) terms, each of size O(M ?2 ). Since BM is symmetric, n ! ?n, this sum converges to ?4G (x) as M ! 1. The same argument applies to the n? term, and the desired equality is veri ed. The same argument works for Gh .

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

13

We still must show that the replacement above P did not change (2.31); i.e., we must show that the absolutely convergent sum (r(n+ ) + r(n? ) ? 2r(n)) is zero. It is the limit as M ! 1 of X? X?

r(n+ ) ? r(n) + (2.33) r(n? ) ? r(n) : BM

BM

We can shift the index in the second sum and write it as the sum of r(n) ? r(n+ ) over a set which is BM , except for O(M) terms each of size O(M ?2). Within BM , the terms cancel exactly, and therefore the expression above approaches 0 as M ! 1. Finally we estimate the error introduced by the regularization of the Green's function for single or double layer potentials on a doubly periodic surface. Theorem 2.2. Let  ! x be a smooth mapping from R2 to R3 so that x()?(; 0) is doubly periodic with period 2 and (1.13) is satis ed. Let f be a smooth periodic function of , G the periodic Green's function, and Gh the regularized version. Then for any 0 (2.34)

Z

Z



Gh (x(0) ? x())f() d ? G (x(0) ? x())f() d  Ch3

where the integral is over a period square in . Proof. We can take the domain of integration to be

(2.35) B0 = f = (1 ; 2) 2 R2 : j ? (0) j <  ;  = 1; 2g : To separate out the singularity, we introduce a cut-o function  on R2 so that () = 1 near  = 0 and () = 0 for jj > =2. Now de ne for all  Gh0 () =

(2.36) Gh0

X

n2Z

2

Gh (x(0) ? x() ? 2(n; 0))(0 ?  ? 2n) :

Then is periodic in , since adding a period to  only shifts the sum; only one term can be nonzero. Similarly we de ne G0 from G using the same . Let Gh1 = Gh ? Gh0 and G1 = G ? G0 , as functions of . Now Gh1 and G1 are smooth and periodic. We rst argue that their dierence contributes little to the integral. From (2.2), we have (2.37) jGh (x(0) ? x()) ? G(x(0) ? x())j  Ch4jx(0) ? x()j?5 provided jx(0) ? x()j is bounded away from zero. By (1.13), the latter is true if j0 ? j is bounded away from zero, as it is on the support of 1 ? . The same is true with  replaced by  + 2n. Consequently Gh1 ? G1 is O(h4), uniformly in ; the nth term decays like jnj?5, so that the sum of estimates converges. Thus the integral corresponding to the dierence Gh1 ? G1 is also O(h4). We are now reduced to the case where Gh and G are replaced by Gh0 and G0 . We note that for  2 B0 , j0 ?   2nj >  for n 6= 0, so that  = 0 there. Consequently Gh0() = Gh (x(0) ? x())(0 ? ) for  2 B0 , and similarly for G . It remains to obtain a local estimate near the singularity at 0. For convenience we assume from now on that 0 = 0 and x(0) = 0. We also may as well assume

14

J. THOMAS BEALE

that f() has support near  = 0, because of the  factor. The proof now amounts to showing that Z

(2.38) The integral is

R2

Z

(Gh (x()) ? 

G(x())) f() d

 Ch3 :



1 f() d ; r = jx()j : ? 2 (s(r=h) ? 1)) 4r (2.39) R Let (; ) be polar coordinates for ; the integral becomes Z 2 Z 1 (2.40) ?(4)?1 (s(r=h) ? 1)) r f() d d : 0 0 We will change variables in the inner integral from  to r. First we write the Taylor formula for x() as x() = J + q() + O(3 ), where J = (@x=@)(0) and q is bilinear. As a consequence we can write (2.41) r= = jx()j=jj = a() + b() + (; )2 ; where is smooth in (; ), regarded as independent variables. Here (2.42) a() = jJj=jj ; b() = J
q()=jj2jJj : For later use we note that a is even and b is odd in , i.e., a( + ) = a(), b( + ) = ?b(). Assuming the support of the cut-o function  was chosen small enough, we can invert (2.41) to obtain (2.43) =r = a~() + ~b()r + ~ (r; )r2 ; where ~a = 1=a, ~b = ?b=a3 , so that a~; ~b have the same parity as a; b, and ~(r; ) is smooth. We now convert the -integral above to Z 1 (s(r=h) ? 1)) r @ (2.44) @r f() dr : 0 Only the low powers of r will matter, and we use Taylor expansions. We have =r = a~ + ~br + O(r2 ) and @=@r = a~ + 2~br + O(r2). We can write the expression for f() as f = f0 + f1 () + O(2 ), with f1 () odd. Converting from  to r, we have f = f0 + a~f1 r + O(r2). Multiplication of these expressions gives ~ )r2 (2.45) (=r)(@=@r)f() = ~a2f0 + (3~a~bf0 + a~3 f1 )r + f(r; with f~ smooth in r; . When we substitute this into the integral, the constant term contributes zero because of the moment condition (2.9). The term linear in r has a coecient which is odd, so that it also contributes zero after integration in . The remaining part of the r-integral is now Z 1 Z 1 2 3 ~ ~ )(r=h)2 d(r=h) : (s(r=h) ? 1)) f(r; )r dr = h (2.46) (s(r=h) ? 1)) f(r; 0

0

The last integral is bounded uniformly in h, since f~ is bounded and s ? 1 decreases rapidly. The total error has been reduced to the last expression and is therefore O(h3) as claimed. Next we have the analogue of Theorem 2.2 for double layer potentials.

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

15

Theorem 2.3. Let  ! x be the mapping of Theorem 2.2, and let F() be a smooth vector-valued function of . Assume that for some 0 the vector F(0) is parallel to n0 , the normal to the surface at x(0). Then (2.47) Z

rGh (x(0) ? x())
F() d ?

Z



rG (x(0 ) ? x())
F() d

 Ch3

where the integral is over a period square in . Proof. As in Thm. 2.2, we can reduce the estimate to the case where 0 = 0,

x(0) = 0, and F is supported in a small neighborhood of  = 0. We then need to show that Z 3 (2.48) 2 (rGh (x()) ? rG(x()))
F() d  Ch : R The integral isZproportional to (2.49) (~s(r=h) ? 1))r?3 x()
F() d2 ; r = jx()j : R2 Proceeding as before we obtain the r-integral Z 1 (2.50) (~s(r=h) ? 1)) x()
F()r?3 @ @r dr : 0

Expanding x as x() = J + q() + O(3 ) and F() = F0 + F1 + O(2 ), and noting that J
F0 = 0, we have (2.51) x
F() = c2 ()2 + c3()3 + O(4 ) with c2 even and c3 odd. Now using  = a~r + ~br2 + O(r3 ), we get (2.52) x
F = c2a~2 r2 + (c3 ~a3 + 2~a~b)r3 + O(r4) : Multiplying out, we nd (2.53) (x
F)r?3(@=@r) = a~4 c2 + (~a5 c3 + 2~a3~bc2 + 3~a3~bc2 )r + O(r2) : The coecient of r is odd, and the argument proceeds just as before. 3. Quadrature of singular integrals In this section we develop a general approach to quadrature for singular integrals with regularization. We nd an expansion in powers of h for the error in the trapezoidal rule, similar to that of Lyness [23] and Goodman et al. [12] without regularization, adapting the approach of [12]. We show that the largest error can be identi ed in terms of the Fourier transform of the regularized kernel using the Poisson Summation Formula. A correction can then be found to improve the order of accuracy. The application to single and double layer potentials is given in Theorem 3.7. We derive the correction terms in (1.26),(1.27) and prove that the quadrature errors in these equations are O(h3). For unregularized singular integrals, corrections are more dicult to nd, but methods have been developed in [22], [17]. For later reference, we state the Poisson Summation Formula for a smooth, rapidly decreasing function f on Rd : X X ^ + 2n=h) f(jh)e?ikjh hd = (2)?d=2 (3.1) f(k j 2Z d

n2Z d

16

J. THOMAS BEALE

We will need a general estimate for quadrature of mildly regular functions. The following statement was proved in [11],[1] for F independent of h. The proof of Lemma 2.2 in [1] is direct from the Poisson Formula. Lemma 3.1. Let F(x) be a function of x in Rd such that Dx F is in L1(Rd ) for multi-indices with j j = `, where `  d + 1. Then, with a universal constant C`, (3.2)

X F(nh) hd n2Z d

?

F(x) dx Rd

Z

 C ` h`

X

j j=`

jDx F jL1 :

We begin with a treatment of quadrature errors for singular integrals, regularized on the scale of the grid size. We start with a singular kernel K(x), de ned for nonzero x in d-space Rd , which is homogeneous of degree m for integer m, i.e. (3.3) K(ax) = am K(x) ; a > 0 ; x 6= 0 : We also work with modi ed kernels of the form Kh (x) = K(x)s(x=h) where s approaches 1 for large argument, and s is chosen so that Kh is smooth up to x = 0. Note that Kh (x) = hm K1 (x=h). The following two lemmas will be fundamental. They are adapted from the proof of Lemma 1 in [12]. Lemma 3.2. Let K and s be smooth functions on Rd ? f0g such that K is homogeneous of degree m, in the sense of (3.3), and (3.4) jDk s(x)j  C jxj?jkj ; jxj  1 : For h > 0, let Kh (x) = K(x)s(x=h), and assume Kh continues as a smooth function to x = 0, so that Kh (0) = hm K1 (0). Finally, let (x) be a smooth cut-o function with (x)  1 for x in a neighborhood of 0 and  = 0 outside a bounded set. Now approximate the integral

I=

(3.5)

Z

Kh (x)(x) dx = d

Z

K(x)s(x=h)(x) dx

R Rd by the sum of values at nh, where n is a multi-integer

(3.6)

S=

Then as h ! 0,

X

n2Zd

Kh (nh)(nh)hd =

X

n2Zd

K(nh)s(n)(nh)hd :

S ? I = c0 hd+m + O(h` ) ; where ` is large depending on the smoothness of K(x) and s(x) for jxj  1. For example, if m = 1 ? d, we have a rst order error with a much smaller remainder. In some cases we can identify the constant c0 , as explained below. In the second expression for the sum S, the proper limiting value should be used for the n = 0 term, namely K1 (0)hd+m . The proof of the lemma below also shows that the quadrature error for the unmodi ed K, omitting the n = 0 term, has the same behavior in case m  1 ? d, as shown in [12]: (3.7)

(3.8)

X

n6=0

K(nh)(nh)hd ?

Z

Rd

K(x)(x) dx = c00 hd+m + O(h` ) :

In the regularized case, the n = 0 term is O(hd+m ), so that the conclusion is valid whether or not the n = 0 term is included in the quadrature.

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

17

Proof. We use a simpli ed version of the argument in [12], modi ed for the smoothing of the kernel. We can rewrite the integral I using the substitution y = x=h and the homogeneity of K as Z (3.9) I = hd+m I ; I  d K(y)s(y)(yh) dy : R Similarly we can write the sum as X (3.10) S = hd+m S ; S  n K(n)s(n)(nh) : Now h appears inside I and S only in . As in [12] we dierentiate in h and compare the resulting sum and integral. We nd Z Z I 0(h) = d K(y)s(y) y
r(yh) dy = h?(d+m+1) d K(x)s(x=h) x
r(x) dx R R and X X S 0 (h) = K(n)s(n) n
r(nh) = h?(d+m+1) K(nh)s(n) (nh)
r(nh)
hd n

n

The integrand and summand are nonzero only on the support of r. We take this support to be fx : r1  jxj  r2g for some r2 > r1 > 0. We can regard the last sum as a trapezoidal rule approximation to the integral above. To estimate the error using Lemma 3.1 we need to estimate the derivatives of the integrand. Since jxj  r1, derivatives of K are uniformly bounded. Moreover from (3.4) we have (3.11) jDxk (s(x=h))j  C jxj?jkj ; jxj  h : Thus the L1 norm of each derivative of the integrand is bounded uniformly in h, up to some order limited by the smoothness of K and s. Then by Lemma 3.1, (3.12) jS 0 (h) ? I 0(h)j  C` h`?d?m?1 ; h > 0 ; for large `. Now S (h) ? I (h) is smooth for h > 0 and has a limit at h = 0 provided ` is large enough. Then (3.13)

(S ? I )(h) = (S ? I )(0) +

Z

h

0

(S 0 ? I 0) = c0 + O(h`?d?m )

for some constant c0, so that S(h) ? I(h) = c0 hd+m + O(h` ), as claimed. Lemma 3.3. Let K , s, Kh , m,  be as in the previous lemma, and let f(x) be a smooth function on Rd such that f and its derivatives are rapidly decreasing. Assume that K , s and  are even in x. Then the integral Z Z I = Kh (x)f(x) dx = (3.14) K(x)s(x=h)f(x) dx d d R is approximated by the sum

(3.15)

with error

S=

X

n2Zd

Kh (nh)f(nh)hd = ?

R

X

n2Zd

K(nh)s(n)f(nh)hd :


S ? I = hd+m c0 f(0) + C2h2 +


+ C2`h2` + O(h2`+2) where c0 is as in (3.7) and C2 etc. depend on g. If K is odd and s is even, we have a similar expansion but with hd+m replaced by hd+m+1 , and the leading error proportional to rf(0). (3.16)

18

J. THOMAS BEALE

Proof. The error from (1 ? )f is easily seen to be high order. For f we use a

second order Taylor expansion for f, multiplied by , ~ ; (3.17) (x)f(x) = f(0)(x) + (rf(0)
x)(x) + f(x) where f~ is a remainder of compact support. The constant term in f leads to an integral whose quadrature error is the rst term above plus a higher order part, according to the preceding Lemma. The linear term in f leads to an integrand which is odd in x, provided K and s are even; thus the integral for this term is zero, and the same is true for the sum, since the grid respects the symmetry. ~ Next we estimate the quadrature error coming from f(x). The remainder f~ consists of a sum of terms, each of the form (x)b(x), where b is bilinear and is smooth. The integrand arising from such a term is (3.18) F(x; h) = K(x)s(x=h)b(x) (x) = hm+2 g(x=h) (x) where g is the smooth function Ksb = K1 b. In order to apply Lemma 3.1, we estimate x-derivatives of F. Now from the assumption (3.4) on s and the homogeneity of K we have jDxk g(x)j  jxjm+2?k, so that (3.19) jDxk F(x; h))j  C jxjm+2?k ; jxj  h ; k  0 : On the other hand, for jxj  h, we can use the smoothness of K1 to estimate (3.20) jDxk F(x; h)j  Chm+2?k ; jxj  h ; k  0 : Combining these two pointwise estimates and integrating, we get jDxk F jL1  Chm+2+d?k : (3.21) Lemma 3.1 now implies that the quadrature error is O(hm+2+d ), and we have proved (3.16) with ` = 0. The expansion can be carried further in the same way, as allowed by the smoothness, with zero contribution from the odd order terms in the Taylor series for f by symmetry. If K is odd rather than even, the sum and integral from the even terms in the Taylor series drop out, beginning with the constant f(0). Lemma 3.4. With K and Kh as in Lemmas 3.2, 3.3, the constant c0 in (3.7), (3.16) is given by

(3.22)

c0 = (2)d=2

X

n6=0

K^ 1 (2n)

Proof. The hypothesis on Kh implies that Dj K1 is in L1 (Rd ) for j large; conse-

quently K^ 1(k) and its derivatives decrease rapidly for k 6= 0. We use Lemma 3.3 with a convenient choice of f. Suppose f(0) = 1. Then from (3,16) and (3.3), h?(m+d) (S ? I) = hlim (3.23) c0 = hlim !0 !0

 X

Z



K1 (n)f(nh) ? K1 (y)f(yh) dy :

^ is smooth everywhere and zero outside jkj  1. Now let We choose f so that f(k) Fh (x) = K1 (x)f(xh). Then Fh is smooth and rapidly decreasing, and we can apply the Poisson Summation Formula (with k = 0, h = 1 in (3.1)) to obtain d=2 X F^h (2n) : c0 = hlim (2) (3.24) !0 n6=0

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

19

We need to justify passage to the limit to obtain (3.22). Since Fh is a product, F^h ^ an approximate delta function. is the convolution of K^ 1 with (2)?d=2 h?d f(k=h), ^ ^ We can write Fh (k) ? K1 (k) as an integral and show by a standard argument, using the compact support of f^ and the decay of DK^ 1 , that jF^h (k) ? K^ 1 (k)j  Chjkj?d?1 for k  1. Thus the sum in (3.24) tends to the sum with K^ 1 in place of F^h . We now apply the above lemmas to determine quadrature errors for single or double layer potentials on a doubly periodic surface, using the regularized Green's function of (2.1). We rst deal with the neighborhood of the singularity. For convenience we place the singularity at zero in both  and x coordinates. Recall that Gh (x) is a smooth function of x the form r?1s(r=h), r = jxj, and Gh (x) = h?1G1 (x=h). Lemma 3.5. Let  ! x be a smooth mapping from a neighborhood U of 0 in R2 into R3, with 0 ! 0. Let J be the Jacobian matrix (@x=@)(0), so that x = J + O(2 ). Assume that J has full rank and U is small enough so that (3.25) jJj=2  jxj  2jJj ;  2 U : Let Gh be a regularized Green's function obeying (2.1)-(2.3), and let f() be a smooth function with support in U . Then as h ! 0 (3.26) where

(3.27)

X

n2Z

2

Gh (x(nh))f(nh) h2 ? c0 = 2

Z

U

X

n6=0

Gh (x())f() d = c0f(0)h + O(h3) ;

(G1  J)^(2n)

or with ? as in (2.14), g = (J  J) , and g the inverse,

(3.28)

c0 = 2(det g )?1=2

X

n6=0

?

p




? 2 g n n :

Proof. Since G1 is a smooth function of x, it is a smooth function of jxj2, i.e., we can

write G1(x) = g(jxj2) for a smooth function g of one variable, and therefore Gh (x) = h?1g(2 ), where  = jxj=h. For later use we note that the decay assumption (2.2) for s implies that (3.29) jDk g()j  C jj?(k+1=2) ; jj  1 : We want to relate the composition Gh (x()) to the simpler case where x() is replaced by J, a case to which the earlier lemmas apply. To this end we use a Taylor expansion for g as a function of 2 about 21 , and then set  = jx()j=h and 1 = jJj=h: (3.30) g(2 ) = g(21 ) + g0 (21 )(2 ? 21 ) + g~(1 ; )(2 ? 21 )2 ; with some smooth g~. Further, we can write x() in a Taylor expansion about  = 0, beginning with J, to obtain (3.31) 2 ? 21 = h?2 (p3 () + y()) where p3 is trilinear in  and the remainder y has degree 4 or higher. Combining these, we have an expansion of Gh in terms with successively higher homogeneity

20

J. THOMAS BEALE

(1) in  and h together; we write Gh = Gh(?1) + G(0) h + Gh , as a function of , where Gh(?1) = h?1g(21 ) = Gh (J) = ?(4jJj)?1s(jJj=h) ; (3.32) ?3 0 2 G(0) (3.33) h = h g (1 )p3 () ; with nonhomogeneous remainder ?3 0 2 ?1 2 22 (3.34) G(1) h = h g (1 )y + h g~(1 ; )( ? 1 ) : We multiply the expanded integrand by a cut-o function () which is even and equals 1 on the support of f and assess the quadrature errors from the various terms. For the leading term Gh(?1), we have the integrand Gh (J)f()(). Lemma 3.3 applies with m = ?1 and x, Kh replaced by , Gh  J. Recalling (2.17), we obtain the rst order error term as stated. For the next term we note that p3 is odd in  and 21 is even, so that G(0) h () is odd, and the contribution from f(0) 3 drops out. Using Lemma 3.3, we see that the error from G(0) h ()f()() is O(h ). (1) Finally we verify that the quadrature error resulting from Gh is O(h3), considering the two terms separately. The remainder y in (3.31) is a sum of terms of the form y~()p4(), where p4 is homogeneous of degree 4 and y~ is smooth. Thus by Lemma 3.3 the integrand h?3 g0 (21 )~y ()p4() has quadrature error O(h?3+4+2 ) = O(h3 ). For the last term we need to estimate directly. The remainder g~ in g is

g~(1 ; ) =

(3.35)

Z

0

1

tg00(t21 + (1 ? t)2 ) dt :

Using (3.29) and (3.25) we can conclude that (3.36) jDk g~(1 ; )j  Ch5 jj?(5+k) ; jj  c1 h with c1 a constant large enough so that jj  c1 h implies jxj  h. We can write 2 ? 21 as a sum of terms of the form h?2q3 ()r(), where q3 is trilinear. Thus the last term in (3.34) leads to the integrand F() = h?5g~(1 ; )q3()2r()2 (). We estimate the quadrature error using Lemma 3.1. Using (3.36) we nd (3.37) jDk F()j  C jj1?k ; jj  c1h : On the other hand, for jj  c1 h, we have jDk g~j = O(h?k ), and consequently (3.38) jDk F()j  Ch1?k ; jj  c1 h : Combining these two pointwise estimates we have jDk F jL1  Ch3?k , and by Lemma 3.1 the quadrature error is O(h3 ). The next lemma concerns the double layer potential near the singularity. We recall that the gradient of the regularized Green's function, rGh has the form r?3x~s(jxj=h). We assume it is multiplied by a vector-valued function F() which is normal at the singular point, so that the product is integrable. Lemma 3.6. Assume the mapping  ! x, the set U , and the function Gh are as in Lemma 3.5. Let F() be a smooth vector-valued function with support in U such that F(0) is parallel to n0 , the unit normal to the surface x() at x(0) = 0. Then (3.39)

X

n2Z2

rGh (x(nh))
F(nh) h2 ?

Z

U

rGh (x())
F() d = C0h + O(h3) ;

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

21

for some C0 depending on F . In particular, C0 = 0 and the error is O(h3 ) under the stronger assumption (3.40) F()
n0 = O() and F T () = O(2) as  ! 0 : where F T () is the projection of F() orthogonal to n0. Proof. We use an argument like the proof of Lemma 3.5. The function s~ in rGh

(cf. (2.11)) has the form s~() = 3 g(2 ), where g is a smooth function, so that (3.41) rGh (x) = h?3g(2 )x ;  = jxj=h : Let x() = x1 () + x2 () + : : : be the Taylor expansion for x() about  = 0, with x1 = J linear, x2 bilinear, etc. As in (3.30) we expand g(2 ) about 21 , with 1 = jx1j=h, and substitute for 2 ?21 in terms of the x's. We collect terms according to degree in ; h together. In this way we write rGh = K?2 + K?1 + K0 + K1 + ::: as functions of . The rst two are (3.42) K?2 = h?3g(21 )x1 ; (3.43) K?1 = h?3g(21 )x2 + 2h?5g0 (21 )(x1
x2)x1 : We consider rst the quadrature of K?2 ()
F(). (As usual we multiply by a cuto covering the support of F.) The term K?2 ()
F(0) is zero since x1() is tangent and F(0) is normal by hypothesis. Then Lemma 3.3 implies that the quadrature error has the form C0h+O(h3 ), with C0 determined by the scalar product of x1() with the linear term in the expansion of F. Under the hypothesis (3.40), this linear term is parallel to n0 , so that the O(h) error vanishes. For K?1 , we again have an error of the form Ch + O(h3 ), according to Lemma 3.3. The O(h) term comes from the constant F(0), and it vanishes under hypothesis (3.40). It can be seen that K0 is odd in , so that by Lemma 3.3 it contributes an error O(h3). Similar considerations show that higher order terms beyond K0 lead to errors O(h3 ) or smaller. To obtain terms up to K0 , we need to expand g to the second order, x to the third order, and 2 ? 21 to the fourth degree terms, all with remainders. Finally, we can estimate the quadrature errors due to the remainders, in a manner similar to Lemma 3.5, and show that they are O(h3 ). The remainder in g(2 ) has the form g~(; 1 )(2 ? 21 )3 , where g~ is an average of g000 . In this case g(2 ) decays like ?3 , and the assumptions on the decay of s~ imply that jDk g~(; 1 )j  C jj?(9+k) for jj  c1 h. Otherwise the treatment of this term proceeds as before. Now that we have treated the quadrature of the single or double layer potential near the singularity, we can obtain similar results for the integral over the doubly periodic surface. Theorem 3.7. Let  ! x be a smooth mapping from R2 to R3 such that x() ? (; 0) is doubly 2-periodic and (1.13) holds. Let f be a smooth periodic function of , and Gh a regularized periodic Green's function as in (2.1)-(2.3), (2.26). Let j = jh for j 2 Z2 , and xj = x(j ). Then for each j 2 Z2 , (3.44) Z X  2 Gh (xj ? x` )f(` ) h ? Gh (xj ? x(0 ))f(0 ) d0 = j f(j )h + O(h3 ) `2I

22

J. THOMAS BEALE

where I is index set (1.10) and the integral is over one period square. Here

(3.45)

j = 2(det g (j ))1=2

X

06=n2Z2



? 2

q

g (j )n n



with ? as in (2.14) and g the inverse metric tensor. ~ 0) be smooth and periodic in ; 0, and let n(0) be the unit Similarly, let f(; normal vector at x(0). Then the quadrature error for the integral Z

(3.46)

~ j ; 0) d0 : rGh (xj ? x(0 ))
n(0)f(

~ 0) = 0 for  = 0, the error improves to O(h3 ). is O(h). If f(; Proof. We need to reduce the general case to the local estimates of Lemmas 3.5,3.6.

As in the proof of Theorem 2.2 we can introduce a cut-o function () with support near zero and write Gh (xj ? x(0 )) as a sum Gh0 (0)+Gh1 (0), depending on j, where the singularity is localized in Gh0 (0), and Gh1 (0) is smooth, with derivatives bounded uniformly in h. We take the domain of integration to be the period square Sj centered at j . We can now split the single layer integral into two parts. The integral with Gh1 (0) has smooth periodic integrand, and its quadrature is high order accurate. The remaining part of the quadrature error comes from the integral with Gh0 . As in the proof of Theorem 2.2, G0 (0 ) = Gh (xj ? x(0 ))(j ? 0) for 0 2 Sj . Thus the single-layer integral is reduced to the case of Lemma 3.5; we need to choose the support of  so that (3.25) holds, depending on the mapping x(). For the double layer integral we use a similar argument to reduce to the case of Lemma 3.6; the hypothesis (3.40) holds when f~ = 0 for  = 0. We can apply this theorem to the discretization error for the four sums representing integrals in the discrete integral equation (1.25)-(1.27). This error results not only from the quadrature but also from the replacement of the derivative by the discrete operator Dh , and we begin with the latter. If f is a smooth, doubly periodic function of , we can regard the values of Dh f on the grid as restrictions of a smooth function Dh f, de ned for all  using the Fourier series for f. Then Dh f = Df + h3r, where r is a smooth function of , depending on h, with a high norm bounded uniformly in h. Thus Xk(h) , the tangent vector computed using Dh , diers from the exact tangent vector Xk by a smooth function of order h3, and similarly for N, n, Xk . In the same way, ?Dh diers from D to O(h3 ). We rst compare the integral (3.47)

I

Z

rT Gh (xj ? x(0))
n(j )n (0)jN(0)j d0

as in (1.26) with the sum X S  rTh Gh (xj ? x` )
n(jh) (n )` jN`(h) j h2 (3.48) `

and derive a correction so that the remaining error is O(h3 ). (Upper (h) indicates the use of Dh . We have Gh in the integral since we have already treated the error in regularizing G in x2.) We start by writing out rT in the integral and integrating

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

23

by parts to obtain (with sum over k = 1; 2 and Dk the derivative in 0k ) (3.49)

I=

Z

Gh (xj ? x(0))Dk fXk (0)
n(j )n(0 )jN(0)jg d0 :

Using the remarks above, we can replace Xk (0) by Xk(h) (0 ) etc., with pointwise  . Since G is integrable, the error O(h3 ). We also replace Dk outside with ?Dh;k h 3 error in the integral is O(h ). We now have (3.50) Z n o  X (h) (0 )
n(h) n (0)jN (h) (0 )j d0 + O(h3) : I = ? Gh (xj ? x(0))Dh;k k j Similarly we write out rTh in the sum and use summation by parts, obtaining o n X  X (h)
n(h) (n)` jN (h) j h2 (3.51) S = ? Gh (xj ? x`)Dh;k j k` ` `

 acts on the ` variable. With the integral and sum in this form, the above where Dh;k theorem applies directly; up to O(h3 ), the error is the O(h) term with coecient

f(j ) given by (3.52)

n

o

 X (h)
n(h) (n )` jN (h) j : rj = ?Dh;k j k` ` `=j  To simplify this, we can replace ?Dh by D, use the fact that Xkj(h)
n(jh) = 0, and then replace D with ?Dh , committing O(h3) errors, to rewrite (3.52) as  X (h) )j
N (h) (n)j : rj = ?(Dh;k (3.53) j k In summary, the integral (3.47) and the sum (3.48) dier by (3.54) S ? I = j rj h + O(h3) : It is this correction term p that appears in (1.26),(1.29) except for the operator Ph . The smoothing by (k h) contributes an error of order h
h2 = h3 provided () ? 1 = O( 4 ) as in (1.33). The rst sum in (1.27) can be treated analogously. The correction simpli es since the last factor vanishes at ` = j and we obtain (1.30). We can use a similar argument to show that the two sums in (1.26),(1.27) with normal derivative of Gh dier from the corresponding integrals by O(h3). After replacing D with Dh again, we have a sum minus an integral with kernel rGh
N`(h) . We cannot apply the second part of Theorem 3.7 directly, since N`(h) is not exactly normal at x(` ). However, we can write it as the exact N(` ) plus a smooth remainder of order h3. The theorem applies to the main part, since the integrand has a factor which vanishes on the diagonal. For the remainder, an argument like Lemmas 3.5,3.6 shows that the quadrature for rGh has O(1) error, and this is multiplied by O(h3) in our remainder term. The total quadrature error is thus O(h3). 4. Discrete pseudodifferential operators In this section we de ne a class of discrete pseudodierential operators and derive basic properties such as boundedness on L2h . We use a symbol class which allows us to compose operators by multiplying symbols, in some cases, with bounded remainder. We prove a discrete version of Garding's inequality and related facts which will be applied to the operator  in x6. We interpret discrete integral operators

24

J. THOMAS BEALE

in this context and express the symbol in terms of the Fourier transform. This connection will be used in x5. For the usual theory of pseudodierential operators, see e.g. [30]. For an extensive theory of discrete operators with high wavenumber cut-o, see [25]. We will work with operators acting on functions de ned on the grid Ih, with I as in (1.10). The general form of a discrete pseudodierential operator A is XX (4.1) Af(jh) = (2)?2 a(jh; `h; k; h)eik(j ?`)hf(`h) h2 or (4.10)

k2I `2I

Af(x) = (2)?2

X X

k2I y2Ih

a(x; y; k; h)eik(x?y)f(y) h2

x 2 Ih :

We generally assume a(jh; `h; k; h) is extended periodically to all integer j; `; k. We will often omit the h-dependence in writing a, although it will be important that our estimates have constants independent of h. In the special case where a is independent of the second variable `h or y, we can rewrite the operator (4.1) as X  (4.2) Af(jh) = a(jh; k; h)eikjhf(k) k

where f is the discrete Fourier transform (1.11). This is the discrete analogue of the standard form of a pseudodierential operator. On the other hand, if a is independent of the rst variable jh or x, we have from (4.1) X (4.3) (Af)(k) = (2)?2 a(`h; k; h)e?ik`hf(`h) h2 : `

As usual we refer to a as the symbol of the operator A. The following lemma is a direct consequence of the de nitions. Lemma 4.1. Suppose A is an operator of type (4.2) with symbol a(jh; k) and B is an operator of type (4.3) with symbol b(`h; k). Then the composition AB is of type (4.1) with symbol a(jh; k)b(`h; k). The adjoint of A is an operator A of type (4.3) with symbol a(`h; k). Often we will want to assume that the symbol a(jh; `h; k; h) is smooth in jh; `h with order m in k, at least in the sense that (4.4) jD+s a(jh; `h; k; h)j  C(jkj + 1)m with C independent of jh; `h 2 Ih, k 2 I, as well as h. Here s is a multi-index, arbitrary up to some large order, and D+ is a forward divided dierence in jh or `h. The next lemma gives the L2h boundedness of an operator with such a symbol with order m = 0. Lemma 4.2. Assume (4.4) holds with m = 0 and with s large. Then (i) The discrete transform a(; ; k; h) of a with respect to (x; y) satis es (4.5) ja(; ; k; h)j  C(1 + j j2 + jj2)?=2 with  large depending on s. (ii) The operator A de ned by (4.1) is bounded in norm on L2h uniformly in h.

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

25

Proof. For (i), we note that the ve-point Laplacian of a in jh is bounded by

hypothesis, and therefore its transform ?4h?2 (sin2(1 h=2) + sin2(2 h=2))a(; ; k) is also. This is at least (2=)2 j j2jaj in magnitude. Similarly, using powers of the discrete Laplacian in jh or `h, we obtain the estimate (4.5). Part (ii) can be proved by adapting a standard argument, e.g. [30], Prop. II.6.1. With x; y 2 Ih we can write X X (4.6) a(x; y; k) = a(; ; k)eixeiy ; A = A ; (4.7)

Now (4.8)

;2I ;2I XX ? 2 ix i`h ik ( x A f(x) = (2) a(; ; k)e e e ?`h) f(`h) h2 : k `

A f(x) =

X

k2I

 ? ) : ei(k+)xa(; ; k)f(k

Each factor is periodic in k, and we can shift k to rewrite this as X  : A f(x) = ei(k++)x a(; ; k + )f(k) (4.9) k2I

 Thus A f(x) is times a function whose transform is a(; ; k + )f(k). 2 2 ? = 2 It follows from (i) that jA f jL2  C((1 + j j + jj ) jf jL2 . Then jAf jL2 is bounded by the sum of such estimates; we obtain the result provided  is large enough so that the sum over ;  has an upper bound independent of h. For further results we need an assumption that the undivided dierence in k, written as
+(k), reduces the order of the symbol by 1. For this reason we de ne a symbol class somewhat like that of the continuous case. We will say that a(x; y; k; h) is a symbol of class S m , for real m, if a(x; y; k; h) is de ned for all x; y 2 R2 and integer k, is periodic in each variable, and satis es the estimate (4.10) j
n+(k)Ds a(x; k; h)j  C(1 + jkj)m?jnj ; k 2 I ; jnj = 0 or 1 where n and s are multi-indices and Ds is an (x; y)-derivative, with s up to large degree. Note that (4.10) is required to hold even at the border points of I. We are primarily interested in the special case of a symbol for an operator of type (4.2) or (4.3). The following lemma shows that the composition of two operators of type (4.2) and (4.3) with total order  1 can be written as an operator of type (4.2) plus a bounded remainder. Lemma 4.3. Let a(x; k; h), b(y; k; h) be symbols of class S m1 , S m2 , m1 + m2  1. Let A; B be the operators of type (4.2), (4.3) corresponding to a(jh; k; h); b(`h; k; h). Also let M be the operator of type (4.2) with symbol a(jh; k; h)b(jh; k; h). Then M ? AB is bounded as an operator on L2h , uniformly in h. Proof. The operator M ? AB has type (4.1) with symbol a(x; k) (b(x; k) ? b(y; k)). We want to perform a summation by parts using   (4.11)
+ eik(x?y) = ei(x ?y ) ? 1 eik(x?y) ; where
+ is the forward dierence in k . We can write   X (4.12) q (x; y; k) ei(x ?y ) ? 1 b(x; k) ? b(y; k) = ei(+)x

 =1;2

26

J. THOMAS BEALE

with

q1 = b(x1 ; x2i(;xk)1 ??y1 )b(y1 ; x2; k) ; q2 = b(y1 ; x2i(;xk)2 ??y2 )b(y1 ; y2; k) : e ?1 e ?1 The denominator in q vanishes only where x = y modulo 2, and q can be extended there so that q is smooth and periodic in x; y. It can be checked that q satis es the same estimates as b. We now have X (M ? AB)f(x) = (2)?2 a(x; k)q (x; y; k)
+ eik(x?y)f(y) h2 (4.14) (4.13)

;k;y

or with a summation by parts X (4.15) (M ? AB)f(x) = ?(2)?2
? (a(x; k)q (x; y; k)) eik(x?y)f(y) h2 ; ;k;y

s (aq ) is using the periodicity in k. The estimates for a; b; q imply that
? Dx;y uniformly bounded for large s, independent of k; h, and the same holds if the (x; y)derivatives are replaced by dierences. Finally we can conclude from Lemma 4.2 that the operator M ? AB is L2h -bounded, independent of h.

The next three results use this lemma. It is convenient to use an operator E on grid functions representing an absolute rst derivative, de ned as  (4.16) k 2I: (Ef)(k) = (1 + jkj2)1=2f(k) More generally, we de ne E r as multiplication by (1 + jkj2)r=2 in the transform for real r. It is important to us that this factor is itself in class S r , i.e., rst dierences in k of the periodic extension reduce the order in k. Lemma 4.4. Let A be an operator of type (4.2) with symbol a(x; k; h) in the class S 1 . Then there are operators B1 ; B2; B3 , bounded on L2h uniformly in h, so that (4.17) A = B1 E ; A = EB2 ; A = E 1=2B3 E 1=2 : Proof. Let b(x; k; h) = a(x; k; h)(k)?1 for k 2 I, where (k) = (1 + jkj2)1=2, and extend periodically in k. Then b is a symbol of class S 0 . Let B1 be the operator of type (4.2) with symbol b(x; k; h). Then B1 is bounded uniformly in h, according to Lemma 4.2, and A = B1 E. This gives the rst equation. For the second, we start with B0  B1 , a bounded operator of type (4.3) with symbol b(y; k; h). We apply Lemma 4.3 with E; B0 in place of A; B, and m1 = 1, m2 = 0. We conclude that EB0 diers by a bounded operator from the type (4.2) operator with symbol (k)b(x; k; h), namely A. Thus we have A = EB0 + B2 for some B2 , or A = E(B0 + E ?1 B2 ), which is the second equation except for notation. The third equation is proved similarly. Let R be the operator of type (4.2) with symbol r(x; k; h) = a(x; k; h)(k)?1=2 for k 2 I, so that r is of class S 1=2 , A = RE 1=2, and r = (k)1=2 b. We apply Lemma 4.3 to the operators E 1=2, B0 , with B0 as above, and conclude that R = E 1=2 B0 +B4 for some bounded B4 . Then A = (E 1=2 B0 + B4 )E 1=2 = E 1=2 (B0 + E ?1=2 B4 )E 1=2. Next we derive a result like Garding's inequality, relating the positivity of a symbol of class S 1 to the positivity of the operator. The resulting estimate is crucial for the numerical stability estimates of x6.

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

27

Lemma 4.5. Let a(x; k; h) be a symbol of class S 1. Assume that a(x; k; h)  c0jkj for k 2 I and all x, with constant c0 > 0. Let A be the operator of type (4.2) with symbol a(x; k; h). Then A ? A is L2h -bounded, uniformly in h, and for all grid functions f and some C  0, X  j2 ? C jf j2L2 : Re (Af; f)L2h  c20 jkjjf(k) (4.18) h k2I

Proof. We mimic the standard proof (e.g., [30], Lemma II.6.2), using Lemma 4.3.

We may assume that a(x; k; h)  c0(k) for k 2 I, where (k) = (1 + jkj2)1=2 , by adding a constant to a without aecting the conclusion. We set b(x; k; h) = a(x; k; h)=(k) for k 2 I and extend periodically in k. Then b is a symbol of class S 0 . Let a1(x; k; h) = a(x; k; h) ? (c0=2)(k) for k 2 I, and let A = A0 + A1 , where A0 has symbol (c0 =2)(k) and A1 has symbol a1(x; k; h), k 2 I. Then a1(x; k; h) = (k)(b(x; k; h) ? c0 =2) for k 2 I; note the last factor is at least c0=2. We want to approximate A1 by an operator of the form RR , which is necessarily positive. Let r(x; k) = (k)1=2 [b(x; k; h) ? c0 =2]1=2. We check that r is a symbol of class S 1=2 : The bound for Dxs r follows from the one for Dxs b. As for the k-dierence, the estimate for
+(k)Dxs b implies that j
+(k)Dxs [b(x; k; h) ? c0 =2]1=2j  C(k)?1 and therefore j
+(k)Dxs r(x; k; h)j  C(k)?1=2. Thus Lemma 4.3 applies to the operator R with symbol r(x; k; h) and the adjoint R with symbol r(y; k; h). We conclude that A1 = RR + B with B bounded. Since RR and A0 are self-adjoint, A ? A = B  ? B, and the rst statement is proved. Further, we have (4.19) (Af; f)L2h = (A0 f; f)L2h + jRf j2L2h + (Bf; f)L2h and the second statement follows. The next lemma expresses the absolute derivative E in terms of the operator A of Lemma 4.5; this is needed in x6. Lemma 4.6. Let A be as in Lemma 4.5. Then for each h there are operators B4 ; B5, bounded on L2h uniformly in h, so that E = B4 A + B5 . Proof. Let Q be the operator of type (4.3) with symbol b?1 , where b = a= for k 2 I. Lemma 4.3 applies with A; Q in place of A; B, and m1 = 1, m2 = 0. We conclude that AQ diers by a bounded operator from the operator of type (4.2) with symbol ab?1 = ; That is, AQ = E + B1 for some bounded B1 . Now Q is bounded, according to Lemma 4.2. Taking adjoints we get Q A = E + B1 . Then, since A ? A  B3 is a bounded operator, we have B4 A = E + B1 ? B4 B3 , with B4 = Q , or E = B4 A + B5 . In dealing with discrete potentials, we will apply the above to discrete integral operators of the form X Af(jh) = K  (jh; jh ? `h)f(`h) h2 (4.20) `2I

where K  (x; z) is periodic in both x and z, and K  has the form X (4.21) K(x; z + 2n) K  (x; z) = n2Z 2

for a function K(x; z) which is smooth in x and z, periodic in x, and rapidly decreasing in z. For each j, we can regard the sum (4.20) as a discrete convolution

28

J. THOMAS BEALE

of K  (jh;
) with f, evaluated at jh. Then Af is the inverse discrete transform of a product, X  ; Af(jh) = (2)2 K  (jh; k)eikjhf(k) (4.22) k2I

an operator of type (4.2). Furthermore, since K  has the form (4.21), its transform as a discrete periodic function can be written as a discrete transform of K itself, X K  (jh; k) = (2)?2 (4.23) K(jh; mh)e?ikmh h2 : m2Z 2

The latter can be related to the usual Fourier transform of K, as a function of z, using the Poisson Summation Formula (3.1): X ^ k + 2n=h) : K  (jh; k) = (2)?1 (4.24) K(jh; n2Z 2

This relation will allow us to study operators related to Gh by estimating their transforms, and to derive an important positivity property using Lemma 4.5. 5. Estimates for integral operators In this section we derive boundedness properties on L2h of discrete integral operators with kernels related to the regularized Green's function. We expand Gh near the singularity into terms where derivatives of Gh are composed with the linearized coordinate mapping. The Fourier transforms of such terms are estimated. It is then shown using the results of x4 that various operators are bounded, with gain of derivatives as for the original integrals. We consider operators of the form X (5.1) Af(j ) = K(j ; `)f(` ) h2 `2I

acting on grid functions f. The kernel K will be related to Gh (xj ?x` ) = Gh (x(j )? x(`)). In order to say that certain operators gain derivatives in L2h , we use the operator E de ned in (4.16), representing an absolute derivative in a discrete sense. Thus if A is an operator of type (4.2) such that ja(x; k; h)j  C(1 + jkj2)?n=2 for some n  0, with the same estimate for x-derivatives of a, then AE n is bounded on L2h , uniformly in h, according to Lemma 4.2, .i.e., A gains n derivatives. When AE n is bounded in this sense for n  0 we will say that A is of order ?n as an operator on L2h . The following theorem summarizes properties of discrete operators associated with Gh which will be needed for the stability estimates for the scheme. Theorem 5.1. Suppose that the regularized Green's function Gh satis es the assumptions (2.1)-(2.3) and G^ h (k) < 0. Let A be the operator X Af(j ) = Gh (xj ? x` )f(` ) h2 (5.2) `2I

Then A is of order ?1, i.e., AE is bounded on L2h uniformly in h. Speci cally, A has the form A = ?A(1) ? A(2) where A(1) has order ?1, A(2) has order ?2, and A(1) is an operator of type (4.2) whose symbol a1(; k; h) is of class S ?1 , in the sense of (4.10), and also satis es, for some c1 > 0,

(5.3)

a1(; k; h)  c1(1 + jkj2)?1=2 :

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

29

More generally, suppose x~j approximates xj = x(j ), with jx~j ? xj j  c0h for each j , and x~j ? xj is periodic. Let (; 0) be an arbitrary smooth, periodic function, and let Am be the operator of the form (5.1) which has the kernel (j ; `)Dxm Gh (~xj ? x~` ), where Dxm is a derivative of order jmj for multi-index m. Then Am is of order ?1 for m = 0 and order 0 for jmj = 1. For jmj = 2, the operator hAm is of order 0. In the special case that  = 0 on the diagonal  = 0, the operator Am has order jmj? 2 for 0  jmj  2. Also for arbitrary  the operator with kernel @Gh (xj ? x` )=@n is of order ?1, where n denotes either normal vector n(xj ) or n(x`).

We also need a more special result concerning the discrete derivative of Gh , showing an extra gain of a derivative due to a smooth factor  which vanishes on the diagonal. Theorem 5.2. With Gh as above, let Dh be a rst order derivative operator as in (1.14)-(1.16), and let (; 0) be smooth and periodic, with  = 0 when  = 0. Then the following two operators are bounded on L2h uniformly in h, where both Dh and E act with respect to `: X B1 f(j ) = Dh Gh (xj ? x` )(j ; `)Ef(` ) h2 (5.4) (5.5)

`2I

B2 f(j ) =

X

`2I

E 2Gh (xj ? x`)(j ; `)f(` ) h2

To begin deriving these facts, we rst split the Green's function into a far- eld part and a local part, in order to focus attention on the latter. We wish to restrict the local analysis to a small neighborhood of the singularity in which the coordinate mapping is well approximated by its linearization. Let J() be the Jacobian matrix @x=@ at . Because of (1.13) and the smoothness of x(), there is some 0 small enough so that 0 < =2 and (5.6) jJ()( ? 0)j=2  jx() ? x(0)j  2jJ()( ? 0)j ; j ? 0j  0 : Also from (1.13) there exists r0 > 0 with r0 < =2 so that jx()?x(0 )j  r0 implies j ? 0j  0 . We can also assume that jJ()( ? 0)j  r0 implies j ? 0 j  0 . We now choose a cut-o function  : R3 ! R so that (x) = 1 for x near 0 and (x) = 0 for jxj > r0 . Thus (5.6) holds whenever x() ? x(0) is within the support of . We write the regularized Green's function as Gh = Gh + (1 ? )Gh . The periodic version is then split as Gh = Gh0 + Gh1 , as in Thm. 2.2. The far- eld part Gh1 is smooth and periodic. Consequently a discrete integral operator with kernel Dm Gh1 (xj ? x` ) is L2h -bounded, uniformly in h, for Dm an x-derivative of any order m. The local term Gh0 when evaluated on the surface is (5.7) Gh0(x() ? x(0)) = Gh (x() ? x(0)) (x() ? x(0)) when j0 ?  j  ; the remaining terms in the sum are zero because the small support of . Thus the boundedness properties of a discrete operator with kernel Dm Gh (xj ? x` ) reduce to consideration of Dm (Gh0 )(xj ? x` ), where Gh0 (x) = Gh (x)(x). We now consider the boundedness properties of the operator (5.2). Because of the above remarks, we replace Gh by Gh0. We will make a Taylor expansion of the kernel which will enable us to use the representation (4.20) for the resulting terms. We treat xj ? x` as a perturbation of its linearization yj` = Jj (j ? ` ),

30

J. THOMAS BEALE

where Jj = J(j ). Thus we can write Gh0 (xj ? x`) as an expansion in terms Dxm Gh0 (yj`)(zj` )m , summed over m with remainder, where zj` = x(j )?x(`)?yj` , and m is a multi-index. We can further expand zj` in powers of j ? ` , quadratic or higher. We then have a linear combination of terms, a typical one being, for j` ? j j < , (5.8) Dxm Gh0 (yj` )(j ? ` )p = Dxm Gh0 (Jj (j ? `))(j ? ` )p : Here p is another multi-index, with jpj  2jmj. The discrete integral operator with (5.8) as kernel, extended periodically, has the form (4.20), and therefore (4.22), with symbol given by (4.24). Thus we can derive boundedness properties of this operator by estimating the Fourier transform of the kernel p Dxm Gh0 (Jj ) with respect to . The relevant estimates are stated in the following lemma. The rapid decay for kh large is needed for the convergence of the sum (4.24). These estimates and the previous remarks show that the operator with kernel (5.8) gains jpj +1 ?jmj derivatives, i.e., it is of order jmj ? jpj ? 1 on L2h in the sense de ned above. Lemma 5.3. Let Kmp () = pDxm Gh0 (J), where J is the Jacobian matrix at some xed 0 . Assume that jpj + 1  jmj. Then the Fourier transform K^ mp satis es the following estimates uniformly in h: 8 > jkj  1 ; C jkjjmj?jpj?1 ; (5.9) 1  jkj  C0=h ; : Ch?jmj+jpj+1 jkhj?n ; jkj  1=h with n large. Similar estimates hold for the derivatives of K^ mp with respect to 0. The same estimates hold for jkj  1 if Gh0 is replaced by Gh . For p = 0 and jmj = 2, we have the same estimate for jkj > 1=h, and jK^ mp (k)j  Ch?1 for jkj  C0=h. In the rest of this section we present the proof of the Lemma 5.3, then the proofs of Theorems 5.1 and 5.2, and one nal lemma. Proof of Lemma 5.3. In estimating this transform it is helpful to note that the xderivatives on Gh can be replaced by -derivatives, as we now explain. First we extend the Jacobian J : R2 ! R3 at 0 to a nonsingular operator J~ : R3 ! R3 by de ning ~ 1; 2; 3) = J(1; 2) + 3N ; (5.10) J( where N = X1  X2 and X1 ; X2 are the tangent vectors at j , as before. Now since Gh is radial, (rGh )(J) has the direction of J for  2 R2 , and thus is perpendicular to N, so that (rGh )(J) = (Gh  J)1 X1 + (Gh  J)2 X2 . Thus Dx Gh can be replaced with D (Gh  J). Next we note that the matrices of second ~ = J~ (Dx2 Gh )J.~ This shows that Dx2 Gh can derivatives are related by D2 (Gh  J) be expressed in -derivatives of Gh  J~ for  2 R3 . However, for 3 = 0, @(Gh  ~ 3 = 0 and the same is true for the second partials in 3; 1 or 3; 2. Also J)=@ the second partial in 3; 3 can be rewritten in 1; 2-derivatives plus a multiple of ~ since
Gh = h . Thus we have reduced the second x-derivatives to 1; 2h  J, ~ derivatives plus the h term. For higher derivatives, we can convert Dxm Gh (J) m 3 ~ to D (Gh  J) for  2 R , and then repeat the above process, converting to 1; 2-derivatives, until a term with h appears. In summary, a term of the form

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

31

Dxm Gh (J) with  2 R2 can be replaced by a sum of terms of the form Dm (Gh  J), with 1 ; 2-derivatives only, and Dxm?2 h (J). We rst estimate the transform of Fmp () = p Dm (Gh (J)). (We ignore the cut-o  for now and account for its eect later.) First, since
Gh = h , we have G^ h (k) = ? ^h (k)=jkj2 = ? ^(kh)=jkj2 for the 3-D transforms. Now as in x2, we obtain the two-dimensional transform of Gh  J from Z 1 Z 1 G^ h (M(k; `)) d` = C (5.11) (Gh  J)^(k) = C jM(k; `)j?2 ^(M(k; `)h) d`

?1 ?1 2  ? 1 where k 2 R and M = (J~ ) , with J~ as above. To obtain the transform of p Dm (Gh (J)), we need to multiply this by (ik)m and then apply (i@k )p . We

obtain various terms where e.g. q derivatives apply to ^ and the remaining jpj ? q to km jM(k; `)j?2. It is important that jM(k; `)j2 depends on ` only through an `2 term; this follows from the de nition of J.~ Thus @k jM(k; `)j2 in linear in k, without ` dependence. The resulting term has the form C

(5.12)

Z

1

?1

kshq jM(k; `)j?n(@ q ^)(M(k; `)h) d`

with n  2 and the total degree in k has been lowered by jpj ? q, i.e., s ? n = jmj ? 2 ? jpj + q. If we let  = jkhj and  = `h, and convert the integral to , we see that the integral is dominated by (5.13)

C jkjshq+n?1

Z

1

?1

(2 +  2 )?n=2(1 + 2 +  2 )?N=2 d

for some large N; we have used the smoothness of . We rst estimate this assuming k 6= 0 and jkhj  C0. In this case we drop the last factor in the integral (the resulting integral exists since n  2) and convert the integral to =. We nd that the integral is of order ?n+1 , and consequently the entire term above is bounded by C jkjshq+n?1 ?n+1 = jkjjmj?jpj?1(jkhjq ), and in summary (5.14) jF^mp (k)j  C jkjjmj?jpj?1 ; jkhj  C0 ; k 6= 0 : On the other hand, for jkhj  1, we return to (5.13), drop the 1 in the second factor, and nd that the integral is of order jj?(n+N ?1). The expression (5.13) is then bounded by C jjs?(n+N ?1)hq+n?1?s = C jjs?(n+N ?1)h?jmj+jpj+1 . The power of  = jkhj can be made large negative by taking N large, and we obtain (5.15) jF^mp (k)j  Ch?jmj+jpj+1 jkhj?n ; jkhj  1 : We need a similar estimate for the transform of Pmp () = p (Dxm?2 h )(J). Omitting the p factor for the moment, we have (5.16)

?




(Dxm?2 h )  J ^(k) = C

Z

1

?1

jM(k; `)jm?2 ^(M(k; `)h) d` :

To obtain the transform of Pmp , we apply @kp to the above. If q derivatives apply to the ^ factor and jpj ? q to the rst factor, we obtain an expression (5.17)

Z

1

?1

(k; `)hq (@ q ^)(M(k; `)h) d`

32

J. THOMAS BEALE

where  is homogeneous of degree s = jmj ? 2 ? jpj + q. We now rewrite this with  = jkhj,  = `h and estimate by (5.18)

hq?s?1

Z

1

?1

(2 +  2 )s=2(1 + 2 +  2 )?N=2 d

for large N. For k 6= 0 and  = jkhj  C0 we replace the integrand in (5.18) by (1 + 2 +  2 )?(N ?s)=2. The integral is O(1), and the term is of order hq?s?1 = h?jmj+jpj+1 ; that is (5.19) jP^mp (k)j  Ch?jmj+jpj+1 ; jkhj  C0 ; k 6= 0 : If the exponent of h is  0, we can majorize by C jkjjmj?jpj?1, as in the estimate (5.14) for F^mp . On the other hand, for jkhj  1, we drop the 1 in the second factor in (5.18), and nd that the integral is of order ?n for large n, so that (5.20) jP^mp (k)j  Ch?jmj+jpj+1 jkhj?n ; jkhj  1 : Next we consider the cut-o function . We rst show that the dierence in the transform due to  is negligible for jkj  1, i.e., we show that the transform of (5.21) Qmp () = p Dxm Gh (J)
(1 ? (J)) is rapidly decreasing for large k. If we apply Ds to Qmp with s of large degree, we can estimate using the chain rule and the fact that jDx` Gh (x)j  C jxj?`?1 for jxj  h. We nd jDs Qmp ()j  C jjjpj?jmj?1?jsj, and of course Qmp () = 0 for  near 0 and for  large. Thus Ds Qmp is bounded in L1 independent of h for jsj large. This implies that jQ^ mp (k)j  C jkj?N for jkj  1 with N large. This estimate is small enough that it will not aect the estimates already obtained for jkj  1 when Fmp , Pmp are multipled by (J). There are additional terms where  is dierentiated, but these terms vanish near 0; consequently their derivatives are bounded independent of h, and the transform is again rapidly decreasing. It remains to estimate the transforms of Fmp ()(J) and Pmp ()(J) for jkj  1. In the rst case the transform is Z ? 1 (5.22) p Dm (Gh (J))e?ik(J) d2 : (2) We can integrate by parts so that the m derivatives apply to p e?ik(J). The result is an integrand bounded by jGh(J)j, for jkj  1, over the support of . Since jGh (x)j  h?1 and jGh (x)j  jxj?1 for jxj  h, we can easily show that this last integral is independent of h. For the case of Pmp  we cannot integrate by parts because we have x-derivatives; however, we can use the scaling of h in h to show that the transform is bounded by Ch?jmj+jpj+1 . Finally, combining the estimates (5.14),(5.15),(5.19),(5.20) and the last two paragraphs, we obtain the estimates stated for K^ mp . It can be checked that derivatives in 0 do not seriously aect the estimates, so that the same bounds can be obtained in this case; the dependence on 0 is through J = J(0 ) in Kmp and therefore through M in the transform. Proof of Theorem 5.1. We rst consider the operator (5.2). As explained before, we reduce consideration to Gh0 and form the Taylor expansion in terms (5.8). The operator with kernel (5.8) has the form (4.20), (4.22), with symbol given by (4.24), and it follows from the estimates of Lemma 5.3 and from Lemma 4.2 that this operator is bounded of order jmj ? jpj ? 1. The rst term, with m = p = 0, is

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

33

of order ?1 and the others ?2 or less, since jpj  2jmj. It remains to show that the remainder is of order ?2. Suppose we carry the expansion for Gh0 to terms with jmj < N and the one for zj` to terms of degree  N in . Then we have two kinds of remainders: One has the form FN (j ; `)(zj` )N , where FN is an average of DxN Gh0 over the line from yj` to xj ? x`. The other type of remainder, from expanding zj` , is Dxm Gh0 (yj` )(j ? ` )p g(j ? ` ) where g is some smooth function and jpj  N + 2jmj ? 1. It is routine to check that derivatives up to order N ? 1 of either remainder with respect to ` , treated as a continuous variable, are uniformly bounded. We can use this fact to estimate the remainders crudely. We need to show that RE 2 is bounded as an operator on L2, where R is the operator corresponding to a remainder term. If R has kernel R(j ; `), we can regard (RE 2 f)(j ) as the inner product of R(j ;
) with E 2f. Since E 2 is self-adjoint, this will be uniformly bounded by jf jL2 provided E 2R(j ;
) is bounded uniformly in j ; `, where E 2 acts on the ` variable. It can be shown that the spectral derivative E 2R is bounded by derivatives Dm R for jmj  5 (cf. [31], pp. 250-4), and thus the operator RE 2 is bounded provided N  6. The decomposition A = ?A(1) ? A(2) essentially follows from separating out the leading term in the expansion, namely Gh0(J(j )(j ? ` )). This term contributes to A an operator of the form (4.22) whose symbol is given by (4.24) with K replaced by Gh0  J. Since G^ h (k) < 0, and since (Gh0  J)^(k) ? (Gh  J)^(k) decreases rapidly for large k, as seen below (5.21), we will take A(1) to be a similar operator with kernel K whose transform agrees with (?Gh  J)^(k) for large k. Let  : [0; 1) ! [0; 1] be a smooth function with (r) = 1 for r  1 and (r) = 0 for r  2. We de ne ^ = (jkj)+(1 ? (jkj))(?Gh  J)^(k). Then A(2) = ?A ? A(1) the kernel K by K(k) has order ?2. We proceed to verify the properties of the symbol a1 of A(1) , which is the sum (4.24). The estimates of Lemma 5.3 show that the sum converges; the sum is positive since the terms are. To show that a1 is bounded below by jkj?1 it is enough to show that ?(Gh  J)^(k)  c1 jkj?1 for 1  jkj  C0=h. This can be seen by rewriting (5.11) with (k:`) = jkj(; ) as (5.23)

?(Gh  J)^(jkj) = C jkj?1

Z

1

?1

jM(; )j?2 ^(M(; )jkjh) d

and recalling that ^ is positive and smooth. The estimate for Ds a1 follows from the lemma. With regard to the k-dierence in (4.10), we note that i@(Gh  J)^=@k is the transform of  Gh  J. This transform decays like jkj?2, according to the Lemma. The estimate for
(k)Ds a1 follows, bounding k-dierences by derivatives for jkj  2. Next we consider the other operators, rst ignoring the smooth factor  and assuming x~j = xj . For the operator with kernel Dx Gh (xj ? x` ) the argument above applies for the boundedness, except that the rst term in the expansion now has m = 1; p = 0, so that the leading operator has order 0. A similar remark applies to Dx2 Gh . If the factor (j ; `) is included, we expand it in the second variable in powers of j ? ` and incorporate this into the previous expansion, so that there is no qualitative change. For @Gh =@n(` ) = ?n(` )
rGh , we use a similar expansion for n(` ) in powers of j ? `. The leading term is ?n(j )
rGh(J(j )(j ? `)), and this is zero since Gh is radial and the range of J(j ) is tangent at j ; the remaining operator is of order ?1. In case  = 0 on the diagonal, the leading term in the expansion of  is zero. The other terms have at least one factor of (j ? ` ),

34

J. THOMAS BEALE

eectively raising the index p by 1 in each term, and the order of the operator is changed as stated. Finally we consider the eect of perturbing xj to x~j . First, using the pointwise estimates (2.4) for derivatives of Gh , we can show 8 > jmj = 0 ; C j loghj ; jmj = 1 ; (5.24) : ?jmj+1 `2I Ch ; jmj = 2 for jx~j ? xj j  c0h, and the same for the sum over j; see e.g. Lemma 5 of [16] or Lemma 3.2 of [7]. It follows from Young's inequality that the discrete integral operator with kernel Dxm Gh (~xj ? x~` ), regarded an operator from L2h to itself, has norm bounded by the right side of (5.24) (e.g., cf. [7], pp. 13-14). The same is true if the kernel is multiplied by a uniformly bounded factor such as . Now suppose the operator with kernel Dx Gh (xj ? x` ) is perturbed by replacing xj ; x` by x~j ; x~`. The dierence of the two kernels can be written as an average of Dx2 Gh , evaluated along the line from xj ? x` to x~j ? x~`, multiplied by (~xj ? xj ) ? (~x` ? x` ). If we assume that max jx~j ? xj j  C0h, then we nd, using the above with jmj = 2, that the operator on L2h corresponding to this dierence is bounded in norm by Ch?1
h = C. A similar argument works for the kernel Dx Gh (~xj ? x~`), and this veri es the statement in the theorem regarding this operator. It also improves the statement made above for bounds on L2h operators in the case jmj = 1, replacing O(j loghj) by O(1). That is, the norm of the operator on L2 with kernel Dxm Gh (xj ? x` ), is O(1) if jmj = 0; 1, and O(h?jmj+1 ) if m  2. We can now use this improved statement with jmj = 1 to verify the assertion that Gh (~xj ? x~`) gains one derivative: the error due to the change in the x's is found to be O(h) as an operator on L2 , and hE is a bounded operator on L2 . We can treat the operator hDx2 Gh similarly. In the special case with  = 0 on the diagonal, the pointwise estimates for the kernel are better, and we can show that the operator Dx2 Gh (~xj ? x~` ) is bounded on L2 by an argument very similar to that above. We can then use that fact to show that Dx Gh (~xj ? x~` ) gains one derivative. In turn, this last fact can be used to show that Gh (~xj ? x~`) gains two derivatives. Proof of Theorem 5.2. We prove that B1 is bounded and then remark on B2 . Let us say that Dh acts in the r-direction, so that e.g. Dh g has transform h?1(kr h)g(k) for k 2 I. We will write (k) = (1 + jkj)1=2 for the symbol of E. Regarding the sum as an inner product, we can rewrite the operator as X  ((j ; `)Ef(` )) h2 : B1 f(j ) = Gh (xj ? x` )Dh;r (5.25) `2I

We will compare this with X  Ef(` ) h2 : (5.26) B0 f(j ) = Gh (xj ? x` )(j ; `)Dh;r `2I

The previous theorem tells us that Gh  gains two derivatives, since  = 0 on the diagonal, and it follows that B0 is bounded on L2 , independent of h. (We think of  as E times a bounded operator.) Thus we need only treat the dierence of Dh;r the two operators.

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

35

To express (B1 ? B0 )f, we write f and  in Fourier sums: X  ik`h ; (j ; `) = X (j ; )ei`h : f(` ) = f(k)e (5.27) k2I

 2I

;`

k

Since  is assumed smooth, (j ; ) decays rapidly in . Then (5.28) X X  ik`h h2 (B1 ? B0 )f(j ) = (j ; )j jei`hGh (xj ? x`) b(; k)(k)f(k)e where we have set h?1 (((r + kr )h) ?  (kr h)) = j jb(; k). Since we assume  is

continuous with bounded derivative, we have j ((r + kr )h) ?  (kr h)j  Chjr j, so that b(; k) is uniformly bounded. We think of the k-sum as an operator on f and write this as X (5.29) (B1 ? B0 )f(j ) = (j ; )j jei`hGh (xj ? x`)(E f)(`h)h2 ;`

 The operator  is bounded on where  is de ned by ( f)(k) = b(; k)f(k). L2 uniformly in  as well as h. Finally, we can think of the above as a sum over  of operators applied to E f, with kernel (j ; )j jei`hGh (xj ? x` ). Each such operator gains one derivative according to Thm. 5.1, i.e., when composed with E, it gives a bounded operator on L2 for each , with norm depending on the smoothness in (j ; `) of the factor (j ; )j jei`. To sum over , we need a bound for some high norm of the factor ei` , and this bound grows with . However, assuming  is smooth enough, (j ; ) decays so rapidly in  that the sum is bounded, and we can conclude that B1 ? B0 is bounded on L2 independent of h. It suces to prove the second statement with E 2 replaced by Dh;r , r = 1; 2, since E 2 can be written in terms of these. For this we compare (5.5) with X (5.30) Dh Gh (xj ? x` )(j ; `)Dh f(` ) h2 : `2I

The latter is bounded according to (5.4), and the dierence can be treated as in the above argument. Finally we need one more fact about sums with Dm Gh : Lemma 5.4. Suppose (; 0) is a smooth, periodic function. Let x~ be a grid function approximating the coordinate mapping x so that jx~(j ) ? x(j )j  c0 h for each j . Then the sum X (5.31) Dm Gh (~xj ? x~` )(j ; `)h2 `2I

is bounded uniformly in j and h, where m is a multi-index with jmj = 1, i.e., Dm denotes any rst derivative in x. If  = 0 when  = 0, then the sum is also uniformly bounded for any m with jmj = 2. Proof. The essential point is that, since Gh is radial, DGh is odd in x and D2 Gh

is even. The proofs of the two cases are similar; we prove the second. Only a small neighborhood of the diagonal matters; we assume  is supported on such a neighborhood and replace Gh with Gh . At rst we assume x~ is identical with x. For simplicity we take j = 0 and assume x0 = x(0 ) = 0. Since (0; ) = 0 at  = 0, we can write (0; ) = 1() + 2 () where 1 is odd in  and 2 = O(2)

36

J. THOMAS BEALE

as  ! 0. Thus if we replace  by 1 in the sum and average with the sum where ` is changed to ?` we get 1 X[D2 G (x( )) ? D2 G (x(? ))] ( )h2 : (5.32) h ` h ` 1 ` 2 `2I

Since x is smooth, x(?) = ?x() + O(2), and D2 Gh (x(?` )) = D2 Gh (?x(` )) +R3, where the remainder R3 is O(2` ) times an average of D3 Gh at points near x(?`). When we substitute this in the last sum, the D2 Gh terms cancel because of the evenness property. We are left with the sum of R3 1 . This is of order j`j?1 for ` 6= 0, since D3 Gh is order j`j?4 and 1 of order j`j. This corresponds to an integrable singularity in 2-D, and the sum is bounded independent of h. For the sum with 2 , we arrive more directly at a similar estimate, using 2 = O(jj2) and D2 Gh = O(jj?3). To complete the proof, we estimate the change in the sum when x is perturbed to x~. Again with j = 0, x(0) = 0, the change in D2 Gh is O(h) times an average of D3 Gh at points near x(?`). For some C1 , jx~(`)j > h and jx(`)j > h when j`j > C1h. Then for such ` we again use D3 Gh = O(j`j?4) to bound the resulting sum by X X (5.33) h
j`j?4j`jh2 = h3 j`hj?3 = C : `6=0

`6=0

The remaining `'s are O(1) in number, and it can be checked that their contribution to the sum is O(1). 6. The convergence proof We now prove that the numerical solution converges to the presumed smooth solution with error O(h3 ). As usual we separate consistency and stability, beginning with the former. Assuming an exact solution given, we rst compare the velocity v, as computed by the scheme at one time from x(j ; t); (j ; t), with the exact velocity. First X1 ; X2 are found at j using the discrete operator Dh , and then N; n; X1 ; X2 . The tangential part of v is X rTh  = (6.1) (Dh;k )Xk : k=1;2

The full computed velocity is (6.2) v = rTh  + wn where w is found from the discrete integral equation (1.25), with the operator Kh de ned by (1.26) and fj by (1.27). The following lemma, proved in x7, concerns the solvability of (1.25). Lemma 6.1. Given x(; t) smooth for 0  t  T , the operator I+2Kh is invertible for h suciently small, and the norm of (I+2Kh )?1 as an operator on L2h is bounded uniformly in h and t. It is evident that rTh  and n dier by O(h3) from the exact values. To compare w with n, determined by the harmonic extension of , we note that n satis es the exact integral equation (1.3). The quadrature results in x3 show that n has truncation error O(h3) in the discrete version (1.25). Then Lemma 6.1 implies that w ? n = O(h3 ) at least in L2h . Thus v of (6.2) also diers by O(h3 ) from the actual

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

37

velocity. A similar statement follows for Bernoulli's equation, and in summary the exact solution satis es the following with p = 3: xt = v + O(hp ) ; t = 21 jvj2 ? gx3 + O(hp ) (6.3) in L2h : It seems that O(h3) accuracy is not quite enough for simple treatment of the nonlinear terms in the stability estimates. For this reason we use a version of Strang's method (see [26], x5.6). That is, we construct a modi cation of the exact solution which matches the scheme to higher order, and then compare the computed solution with this modi cation. The following lemma (proved in x7) asserts that this is possible. Lemma 6.2. Suppose x(0)(; t); (0)(; t) is a smooth, exact solution of the water wave equations (1.1), (1.2). Then there exist smooth functions x(1)(; t), (1) (; t) diering from x(0) (; t), (0)(; t) by O(h3), so that if v(1) is found from x(1); (1) according to (6.1),(6.2),(1.25)-(1.30), then (6.3) holds for x(1); (1); v(1) for any p < 4. Also w(1) is within O(h3 ) of (1) n in L2h . For the rest of this section we compare quantities such as v(1) , found by the scheme from x(1) ; (1) at each time, with the corresponding quantities in the computed solution. For simplicity we will drop the superscript (1) for the former. In contrast, we use tildes for the computed solution, determined by (6.4) x~t = v~ ; ~t = 21 jv~j2 ? g~x3 with v~ found from x~; ~ as in (6.1),(6.2),(1.25)-(1.30). We write x for x~ ? x and similarly for other dierences. Thus, e.g. X1 = Dh;1 x, X~ 1 = Dh;1 x~, X1 = X~ 1 ?X1 . Then v = v~ ? v is the error in the computed velocity due to x; . Comparing (6.3),(6.4) we have (6.5) in L2h ; p < 4 : (x)t = v + O(hp ) ; ()t = 21  jvj2 ? gx3 + O(hp ) We will show that x;  = O(hp?1=2) in L2h . With p near 4, it will then follow from the triangle inequality that x~; ~ dier from the original exact solution in L2h by O(h3 ), thereby completing the proof that the scheme converges. Our main task now is to express the error v in computed velocity in terms of x, . To do this we need to nd the errors in such quantities as n, rTh , and w. For the latter we vary the integral equation which determines w. Because we are proving that x;  are almost O(h7=2) in L2h , we will assume a lower order bound on the errors, (6.6) jxjL2h  h3 ; jjL2h  h3 ; jwjL2h  h2 ; 0  t  T : It follows that (6.7) jxjL1h  h2 ; jjL1h  h2 ; jwjL1h  h ; 0  t  T : as well as jDh xjL1h  C0h etc. This assumption is helpful in estimating nonlinear terms in the 's. It will be removed at the end of the argument. In order to make use of the gain of derivatives in certain operators related to Gh , seen in x5, we will use the discrete Sobolev space Hh?1 of grid functions with norm X  j2h2 : jf j2Hh?1 = (1 + jkj2)?1 jf(k) (6.8) k2I

38

J. THOMAS BEALE

We can write f = E
E ?1 f, so that jE ?1f jL2h = jf jHh?1 , with E as in (4.16). Thus if an operator A is of order ?1 on L2h , in the language of x5, then A is bounded from Hh?1 to L2h , independent of h. In the next lemma we note some facts which will be used. Lemma 6.3. We have jf jHh?1  jf jL2h and jf jL2h  Ch?1jf jHh?1 for any f . The  are bounded from L2 to H ?1 , uniformly in h. If is a grid operators Dh;` , Dh;` h h function which either is the restriction of a smooth function or satis es j jL1h  Ch, then j f jHh?1  C0jf jHh?1 , where C0 depends on but not h. Proof. The rst two statements follow from de nitions. For the last, we represent + be the rst-order forward dierence in elements of Hh?1 as dierences. Let Dh;` + is bounded above and direction `, ` = 1; 2. Since the Fourier symbol of Dh;` ? 1 below by jk`j for k 2 I, we can write any f 2 Hh as f = f0 + Dh;+ 1 f1 + Dh;+ 2 f2 , where f0 ; f1; f2 are bounded in L2h by the Hh?1 norm of f. Then, e.g., Dh;+ 1 f1 = Dh;+ 1 ( f1 ) ? (Dh;+ 1 )f1] , where f1] is f1 with a coordinate shift. The factor Dh;+ 1 is uniformly bounded under either hypothesis on , and it follows that Dh;+ 1 f1 has the appropriate bound in Hh?1, and in turn the same holds for f. Lemma 6.4. Assume (6.6). Then for the tangent vectors X` ; X~`, ` = 1; 2, we have jX~ ` ? X` jHh?1  C jx~ ? xjL2h . The same is true for the errors in N , jN j, n, and X` . Proof. The statement for X` is apparent from the de nition of the norm. For the other quantities we use the following general fact: Suppose u is a smooth function of  on the periodic square, possibly vector valued, and suppose u~j is an approximation to u(j ) such that ju~ ? ujHh?1  C jx~ ? xjL2h . Assume (6.6), and let F be a smooth function on an open set covering the range of u. Then (6.9) jF(~u) ? F(u)jHh?1  C 0 jx~ ? xjL2h : The Lemma follows from this fact, since the various quantities are functions of X1 ; X2. To verify this claim, we rst note that for each j , (6.10) ju~j ? u(j )j  Ch?1jujL2h  Ch?2jujHh?1  Ch?2
h3 = O(h) by hypothesis, so that u~ is uniformly close to u. Next with u = u(j ) we write F(~u) ? F(u) = DF(u)(u) + R(~u; u)(u)2, where R is an average of the second derivative of F between u and u~. The rst term is bounded in Hh?1 by C jujHh?1 according to Lemma 6.3, since DF(u) is a smooth factor. For the second term we estimate in L2h . The rst factor is uniformly bounded, and for the rest we have (6.11) j(u)2 jL2h  C jujL1h jujL2h  ChjujL2h  Ch
h?1 jujHh?1 where we have used (6.10). Thus both terms are bounded by jujHh?1  C jxjL2h , and (6.9) is proved. In order to simply expressions for errors, we will not be speci c about terms contributing to v which are bounded in L2h by jxjL2h , jjL2h , since they do not aect the numerical stability of the scheme. For this reason we use the generic notation B(f) for a grid function which is bounded in L2h , uniformly in h, by jf jL2h . Thus we write remainders such as B(x), B().

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

39

In the next lemma we nd the error in the unit normal, n = n~ ? n, due to x = x~ ? x, and also the error in the tangential gradient of an arbitrary function , de ned as in (6.1). Here and below we use the product rule (1.17) for Dh . Lemma 6.5. Assuming (6.6), the error in n due to x is (6.12) n = ?rTh (n
x) + B(x) : Similarly, if  is a smooth function of  and ~ is an approximation on the grid, then

?




?




 rTh  = rTh () ? rTh (x
rTh ) ? n
rTh  n + B(x) + B() ?
= rTh () ? rTh (x
rTh ) + rTh (n
x)
rTh  n + B(x) + B() : Proof. Thinking of n = X1  X2 as a pointwise function of Xk = Dh;k x, k = 1; 2, we introduce the in nitesimal change X n_ = (6.15) (@n=@Xk )
Xk

(6.13) (6.14)

k=1;2

where of course Xk = X~k ? Xk . Then (6.16) n(j ) = n( _ j ) + O(jXk (j )j2) : Dierentiating n
n = 1 shows n_
n = 0, i.e., n( _ j ) is tangential at x(j ), and similarly n_
Xk = ?n
Xk . Then X X n_ = (n_
Xk )Xk = ? (6.17) (n
Xk )Xk : k=1;2

k=1;2 Now since Xk = Dh;k (x), we can write n
Xk = Dh;k (x
n) + B(x), Then the above becomes n_ = ?rTh (x
n) + B(x). We combine this with (6.16) to obtain the result (6.12); the remainder in (6.16) can be estimated in L2h using (6.6) by (6.18) jDh (x)jL1h jDh (x)jL2h  Ch?1jDh (x)j2L2h  Ch?3 jxj2L2h  C jxjL2h :

Next we nd a similar expression for Xk = X~k ? Xk . As before, let X_ k be the in nitesimal change in Xk due to x. Again Xk ? X_ k is quadratic in Dh (x), and _ it can be estimated as in (6.18). We have X_ k
X` = ?Xk
X` and X_ k
n = ?Xk
n. Then X X (X_ k
X` )X` + (X_ k
n)n = ? (6.19) X_ k = (Xk
X` )X` ? (Xk
n)n _ `=1;2

`=1;2

and after moving the `-derivative and replacing n,_
? (6.20) Xk = ?rTh (x
Xk ) + rTh (n
x)
Xk n + B(x) Finally, X ?
(6.21) Dh;k () Xk + (Dh;k )Xk + Dh;k () Xk :  rTh  = k=1;2

The last term can be estimated by C jjL2h using (6.6). Combining (6.20),(6.21) we then obtain (6.14), and (6.13) follows. In the next lemma we nd an expression for v in terms of w. Lemma 6.6. With v as in (6.2) and similarly for v~, we have, assuming (6.6), ?
(6.22) v = rTh ( ? v
x) + (w)n ? n
rTh  n + B(x) + B() :

40

J. THOMAS BEALE

Proof. Varying (6.2) we have ?




?




(6.23) v =  rTh  + (wn) =  rTh  + (w)n + w(n) + (w)(n) : From (6.13) we have ?
?
?
(6.24)  rTh  = rTh () ? rTh rTh 
x ? n
rTh  n + B(x) + B() and from (6.12) (6.25) w(n) = ?wrTh (n
x) + B(x) = ?rTh (wn
x) + B(x) : For the last term in (6.23) we can estimate using (6.6) and Lemma 6.4 (6.26) j(w)(n)jL2h  C jwjL1h jnjL2h  Ch
h?1 jnjHh?1  C jxjL2h : Now we substitute (6.24)-(6.26) into (6.23). We can combine terms in (6.24),(6.25) using (6.2) to obtain (6.22). The calculation of v is now reduced to that of w. To nd w, we subtract the discrete integral equation (1.25) for w from that for w, ~ and obtain an expression with variations in the discrete integrals entering the equation. We state the variations below; they are derived in x7. In the sums, rGh = rGh (xj ? x` ) etc., while in the integrals rGh = rGh (xj ? x(0 )). (6.27) 

(6.28) 

X

`

+ (nj )


X

`

rGh
N` [(n`
nj )w` ? wj ]h2 = Z

X

rGh
N` [(n`
nj )w` ? wj ]h2

`  0 0 (rGh
N(x( ))n(x( ))n (0) d0

rTh Gh
nj w` jN` jh2 = +

X

`

X

+ hB(w) + B(x)

rTh Gh
nj (w` )jN` jh2

` T  rh Gh
rTh(w` n`
x`)jN` jh2 Z

+ (nj )
rT Gh n(0 )jN(0)jd0 + hB(w) + B(x) (6.29) " 

X

X nj
[Dh Gh ; X]`  (rTh ` ? rTh j )h2 ? nj
rGh
N` (rTh ` ? rTh j )h2 ` ` X T  T T = rh Gh
rh (` ? rh `
x` )jN` jh2 ` Z + (nj )
rGh  (N(0)  rT (0))d0 + B(x) + B()

#

When we substitute these into the integral equation for w, we can collect the terms. The integrals in (6.27) and (6.28) can be combined, and the second sum in (6.28) can be combined with the sum in (6.29) using (6.2). The rst terms in (6.27) and (6.28) have the form Kh (w), with Kh de ned by (1.26), except for the correction. It can be checked that the error in the correction terms is of type

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

41

hB(w)+B(x)+B(); the operator hPh Dh is bounded from Hh?1 to L2h because of the high k cut-o. Consequently we get (6.30) ( 12 I + Kh )(w) = L( ? v
x) + (n)
S + hB(w) + B(x) + B() where the operator L and the function S are de ned by X (Lf)j = rTh Gh (xj ? x`)
(rTh f` )jN` jh2 (6.31) (6.32)

Sj =

Z

` Z  rGh  (N(0 )  rT (0))d0 ?

rGh n(0 )jN(0)jd0 :

Lemma 6.1 assures us that (6.30) can be solved for w; the solvability is not aected by the small term hB(w). The following related lemma guarantees that w~ can be found from its discrete integral equation. Lemma 6.7. Under assumption (6.6) the operator I + 2K~ h , obtained by replacing x with x~ in (1.26), is invertible for h suciently small, and the operator norm of the inverse is uniformly bounded. The discrete integral equation can be solved by the simple iteration (1.34).

The principal part of the operator L is of central importance in estimating the growth of numerical error. The following lemma identi es this principal part and its properties. Lemma 6.8. The operator L on L2h de ned by (6.31) has the form L = 21  + B0 , where B0 is bounded on L2h , uniformly in h, and  is an operator with the following properties: (i)  is self-adjoint as an operator on L2h , f is real-valued for real f , and

(6.33)

(f; f)L2h  c0jE 1=2f j2L2h = c0

X

k2I

 j2 : (1 + jkj2)1=2 jf(k)

(ii) There are operators B1 ; : : :; B5 bounded on L2h , uniformly in h, so that

(6.34)  = B1 E ;  = EB2 ; 1 = 2 1 = 2 (6.35)  = E B3 E ; E = B 4  + B5 : Here E is the absolute derivative de ned in (4.16). Proof. The result depends on the assumptions (1.14)-(1.16) for Dh as well as the positivity of ?G^ h . If we write out rTh we have X Lf(j ) = ? Dh; Gh (xj ? x`)(Dh; f` )(X`
X` )jN` jh2 (6.36) `;;

  with ;  = 1; 2, or, after summing by parts and setting b ` = (X`
X` )jN` j =  (g )` jN` j, (6.37) X  [b Dh; f` ]h2 = ? X G b D Dh; f` h2 + (Bf)(j ) : Lf(j ) = ? Gh Dh; h ` h; ` `;;

`;;

 and note that the To see that the remainder is bounded, we use (1.17) for Dh;  operator with kernel Gh , timesa smooth factor, is of order ?1, according to Thm.  Dh; f 5.1. Next we replace b` with bj ; the dierence is an operator applied to Dh;r

42

J. THOMAS BEALE

with the kernel Gh times a function which vanishes at ` = j. This operator has order ?2, by Thm. 5.1, and so the error is again bounded in f. We now have X  Lf(j ) = ? b (6.38) j (ADh; Dh; f)(j ) + (Bf)(j ) ;

where A is the operator of Thm. 5.1, with kernel Gh . There we wrote A as A = ?A(1) ? A(2) , where A(2) is of order ?2, and thus contributes a bounded term to L, while A(1) is an operator of type (4.2) whose symbol a1 (; k) is of class S ?1 with a1(; k)  c1(1 + jkj2)?1=2 for k 2 I. Thus L = 0 + B, where 0 is the operator of type (4.2) with symbol 0 of class S 1 , X 0(; k) = h?2 a1(; k) b ()(k h)(k h) : (6.39) ;

We have 0 (?k) = 0 (k), which implies that 0 preserves real functions. Since the matrix b is positive de nite and symmetric, and h?1j(k h)j  C jkj for k 2 I, we have 0 (; k)  c2 jkj for k 2 I. Thus Lemma 4.5 applies to 0; the estimate (4.18) holds for 0 , and 0 ? 0 is bounded. Now we set  = 0 + 0 + C0, where C0 is a large enough constant. Then  ? 20 is bounded, and the properties of  follow from Lemmas 4.4-4.6. In view of this lemma, equation (6.30) reduces to (6.40) (I + 2Kh )(w) = ( ? v
x) + 2(n)
S + hB(w) + B(x) + B() : We now show that the term Kh (w) is negligible. The two main terms on the right above can be written as Ef + B(x) + B() for some f with jf jL2h  C(jxjL2h + jjL2h ); this follows from (6.12) and (6.34). Then w is (I + 2Kh )?1 Ef plus a remainder, and the main term in Kh (w) is (I +2Kh )?1 Kh Ef. By Lemma 6.1, the last is bounded in L2h norm by jKh Ef jL2h . Theorems 5.1, 5.2 imply that the operator Kh is of order ?1; for the second term in Kh we use the fact that Xk`
nj = 0 when ` = j. Thus jKh Ef jL2h  C jf jL2h , and it follows that Kh (w) = hB(w) + B(x) + B(). We can now absorb Kh (w) into the remainder in (6.40), and then invert (I + hB) acting on w, simplifying (6.40) to (6.41) (w) = ( ? v
x) + 2(n)
S + B(x) + B() : Next we substitute (6.41) into (6.22) to nd a de nitive expression for v. After combining we have (6.42) ?
v = rTh ( ? v
x) + ( ? v
x)n + n
2S ? rT  n + B(x) + B() : (We have replaced rTh  with the exact rT, at a cost of hB(n) = B(x).) However, the third term above is negligible, as a consequence of potential theory: If  is a periodic harmonic function below the surface x(), then (6.43) Z Z ?
@ dS(y) ? rT (x) = 0 2 rT (x)G  n(y) rT (y) dS(y) ? 2 rT (x)G @n(y) where G = G (x ? y) and rT (x) means the part of the gradient which is tangential at x. This follows from the same argument as for the integral equation (1.3), except that we take the limiting tangential derivative at the surface. If we replace G by Gh ,  by , and @=@n by w, the above holds with error O(h), because of the

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

43

regularization of G and the modi cation in Lemma 6.2. The resulting expression almost matches 2S ? rT , except that the latter has rG rather than rT (x)G . But n is tangential except for a B(x) term; see (6.12). Therefore n
(2S ?rT ) has the form hB(n) + B(x) and thus B(x). In summary, we have (6.44) v = rTh ( ? v
x) + ( ? v
x)n + B(x) + B() : We want to replace the discrete approximation v in (6.44) with xt , the actual velocity of the modi ed exact solution. From (6.3) we see that xt ? v is O(hp ) in p?2 L2h and thus O(hp?1) in L1 h . The change in v is then h (B(x)+B()). Since p > 2, this can be absorbed into the remainder. Thus we introduce (6.45) u1 =  ? xt
x and rewrite (6.44) as (6.46) v = rTh u1 + (u1 )n + B(x) + B() : We are now ready to estimate the growth in time of the errors in x; . Because of the special structure of v in (6.46), we can use an argument similar to that in [5], [6]. We will replace x;  by new state variables u1; u2; u3, where u1 is de ned by (6.45) and (6.47) u2 = x
n so that, recalling (6.5), we have in L2h (6.48) u2;t = (x)t
n + B(x) = v
n + B(x) + O(hp ) = u1 + B(x) + B() + O(hp ) : We will choose u3 to account for (x)T while taking advantage of the similarity of the two main terms in (6.46). We will need to know that B4 of Lemma 6.8, and later also B3 , have bounded time derivatives. This follows from their construction in x4 once we know that the t-derivative of the symbol a1 (; k; h) in Thm. 5.1 has bounds similar to those for a1 itself, i.e. as a symbol of class S ?1 . The essential point is that Gh (J) depends on t only through the Jacobian J, and its transform through the related matrix M (see (5.11)); dierentiation in t does not change the degree of homogeneity in  or k, and the estimates for Gh apply to @Gh =@t as well. To de ne u3 , we rst write rTh = B6  + B7 , using (6.35), with B6 ; B7 bounded uniformly in h, t, determined by B4 ; B5. We now de ne (6.49) u3 = (x)T ? B6 (x
n) : Applying @=@t, using the boundedness of B6;t , and dropping less important terms, we nd (6.50) u3;t = (v)T ? B6 (v
n) + B(x) + O(hp ) = rTh u1 ? B6 (u1) + B(x) + O(hp ) = B(x) + B() + O(hp ) with no important terms remaining. Next we derive the equation for u1;t. From (6.5) we obtain t = 21 (v
v) ? gx3 + O(hp ) = v
v + 12 jvj2 ? gx3 + O(hp ) : (6.51) For the other term in u1 we have (6.52) (xt
x)t = xt
v + xtt
x + O(hp ) :

44

J. THOMAS BEALE

When we subtract, we have a term (v ? xt)v. From (6.5), v ? xt = O(hp?1 ) ?1 in L1 h , and from (6.46) v has the form h (B(x) + B()). Then this term is h(B(x) + B()), since p > 3. Combining terms we have u1;t = ?(xtt + ge3 )
x + 12 jvj2 + hB(x) + hB() + O(hp ) : (6.53) For an exact solution of the Euler equations, ?(xtt + ge3 ) = rp, the pressure gradient. At the water surface, this is a normal vector, since p = 0 there. Moreover, we saw in (1.35) that ?rp
n  c0 > 0 at the surface. We have modi ed the exact solution by O(h3), and thus in our case we have (6.54) xtt + ge3 = ?cn + O(h3) ; c = c(; t)  c0 > 0 : The rst term in (6.53) is thus ?cx
n + h3B(x) = ?cu2 + h3 B(x). We need an estimate for the (v)2 term in (6.53). From (6.46) and (6.6) we have jvjL1h = O(h), and also jvjL2h  C jEu1jL2h  h?1=2jE 1=2u1jL2h , and combining these two facts, j(v)2 jL2h  h1=2jE 1=2u1jL2h . We can now estimate the growth of a norm of (u1 ; u2; u3). We de ne a norm so that the two principal terms will cancel in the growth estimates, U = 21 f(u1; u1)L2h + jc1=2u2j2L2h + jc1=2u3j2L2h g  U1 + U2 + U3 : (6.55) We have (6.56) jxj2L2h + jj2L2h  CU ; jE 1=2u1 j2L2h  CU and we can rewrite our equations as (6.57) u1;t = ?cu2 + h1=2r1 + hp 1 ; u2;t = u1 + r2 + hp 2 ; u3;t = r3 + hp 3 where the r's are bounded in L2h by U 1=2, and the 's by a constant. Then (6.58) U1;t = (u1;t; u1) + 12 (t u1; u1) = ?(cu2 ; u1) + h1=2(r1 ; u1) + 12 (B3;tE 1=2 u1; E 1=2u1) + hp (1 ; u1) The second term can be estimated by (6.59) Ch1=2jr1jL2h jEu1jL2h  Ch1=2jr1jL2h h?1=2jE 1=2u1jL2h  CU and the third term by C jE 1=2u1 j2L2h  CU. Recalling that  = B1 E, we estimate the last term by (6.60) Chp ju1jL2h  ChpjEu1jL2h  Chp?1=2jE 1=2u1jL2h  Chp?1=2U 1=2 and in summary (6.61) U1;t = (u1;t; u1) + R1 ; jR1j  C(U + hp?1=2 U 1=2) : For U2 we have U2;t = (cu2;t; u2) + 21 (ctu2 ; u2) = (cu1; u2) + R2 (6.62) with jR2j bounded by C(U + hp U 1=2), and similarly U3;t = R3. When we add, the two main terms in (6.61),(6.62) cancel, and we have Ut  C1(U + h2p?1). Then U(t)  C2 h2p?1 for 0  t  T, and from (6.56) (6.63) jxjL2h + jjL2h  C3hp?1=2 ; 0  t  T :

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

45

It follows from (6.46) and Lemma 6.6 that (6.64) jvjL2h ; jwjL2h  C4hp?3=2 ; 0  t  T : We can now remove the assumption (6.6) used to derive (6.63),(6.64). We can take p near 4, so that p ? 1=2 is almost 7=2. For h small enough, (6.63),(6.64) imply that the L2h norms of x; , hw remain below h3=2 as long as they are less than h3, and therefore they never reach h3 . Thus for small h (6.63),(6.64) hold without quali cation, and in view of Lemma 6.2 the computed solution stays within O(h3) of the actual solution in L2h . 7. Proofs of lemmas Here we present the proofs of Lemma 6.1, Lemma 6.2, estimates (6.27)-(6.29), and Lemma 6.7. Proof of Lemma 6.1. The invertibility of I + 2Kh will follow from that of the operator I + 2K of which it is an approximation, where K is the operator acting on periodic functions of , (7.1)

Z

Kf() = rG (x() ? x(0))
n(x())f(0 )jN(0)jd0 :

Thus Kf is the adjoint double layer potential due to f. With time t xed for the moment, the operator K is compact, as an operator on L2 of the periodic square, so that its nonzero spectrum consists only of eigenvalues. Classical arguments show that the operator 2K and its adjoint 2K have spectral radius strictly less than 1; that is, if  is an eigenvalue of 2K, then jj  (2K) < 1. This fact depends on the periodic geometry; for bounded domains there is an eigenvalue with jj = 1. A proof that (2K) = (2K ) < 1 was given in [3] for double layer potentials in periodic regions in two dimensions; essentially the same proof applies to the present case of 3-D regions with doubly periodic boundaries. Now let K(t) be the operator at time t for 0  t  T. We argue that for some constants " > 0 and C0 > 0, and for all  complex with jj  1 ? ", and all t with 0  t  T, the operator I ? 2K(t) is invertible, and k(I ? 2K(t))?1k  C0 . First, since (2K(t)) < 1, then also (2K(t) + A) < 1 if A is small enough in norm as an operator on L2 ; this is shown in [10], Thm. 1.37. It can be checked that K(t) depends continuously on t, as an operator on L2 , using the smoothness assumptions. Thus the above statement tells us that for each t0 there is a t-interval near t0 such that (2K(t)) is bounded away from 1 on the interval. Since [0; T] is a compact set, the same is true for this entire interval, i.e. there is some " > 0 so that (2K(t)) < 1 ? " for each t 2 [0; T]. This shows that I ? 2K(t) is invertible for each t and  with jj  1 ? ". The inverse depends continuously on t; , and so a uniform bound for the norm follows for  in a bounded set by compactness. Finally, we have a uniform bound for the norm of K(t), and from this we can easily obtain a bound for the inverse for large . This veri es the assertion. The above fact allows us to show that operators close to K(t) have similar properties. Suppose A is any operator on L2 with kAk < (2C0 )?1, where C0 is as above. Then for any t and any  with jj  1 ? ", we have I ? (2K ? A) = (I ? 2K)[I + (I ? 2K)?1A]. By the statement above, each of the last two factors is invertible, and the inverse has norm bounded by 2C0. Thus we have shown that (2K(t) ? A) < 1 ? ", and in particular k(I + 2K(t) ? A)?1 k  2C0. To complete

46

J. THOMAS BEALE

the proof of the lemma, it will be sucient to show that Kh can be regarded as a small perturbation of K when h is small. We rst treat the operator Kh0 , which is Kh with the two terms with wj omitted. To show that I + 2Kh0 is invertible on L2h , we identify it with an operator on L2 ([?; )2); i.e., we extend from the discrete variable j 2 Ih to  2 [?; )2. We partition the square by de ning Bj = f 2 [?; )2 : 0   ? (j ) < h;  = 1; 2g for j 2 I. If we write Kh0 in the form X Kh0 f(j ) = Kh (j ; `)f(` ) h2 (7.2) `2I

we can extend the kernel by de ning Kh (; 0) = Kh (j ; `) for  2 Bj , 0 2 B` . This kernel then gives an integral operator on L2 ([?; )2). Now if f is a function on Ih, we can extend it to a piecewise constant function on [?; )2 by de ning f() = f(` ) for  2 B` . The integral operator applied to this extended f produces a piecewise constant function with values given by (7.2). The L2 norms of the extended f and its image match the L2h norms of the discrete functions. (Cf. [6], Lemma 5; [7], p. 16.) Thus if we show that the extended operator I + 2Kh0 has an inverse, bounded independent of h, then the same will be true for the discrete operator. In view of the above remarks, it will be sucient to show that the discrete operator Kh0 , when extended, is close to the corresponding continuous version. The two terms of Kh0 have the form X J1 f(j ) = N`
rGh (xj ? x` )1 (j ; `)f(` ) h2 (7.3) (7.4)

`2I

J2 f(j ) =

X

Dh Gh (xj ? x` )2 (j ; `)f(` ) h2

`2I where 1 ; 2 are smooth, 2 is zero on the diagonal, Dh is the discrete derivative with respect to ` , and N` is computed from x() using Dh . Thus, in the rst case, we compare J1 extended to L2 with Z (7.5) N(0)
rGh (x() ? x(0))1 (; 0)f(0 ) d0 :

(We have already seen in Thm. 2.3 that the change in the integral due to regularizing G to Gh is O(h3 ).) The error to be accounted for is that from replacing x() by xj , x(0) by x` , and N(0) by N` in the kernel when  2 Bj , 0 2 B` . We use an estimate for the rst derivative of the kernel. We consider Gh instead of Gh , since the remainder is smooth. Since Gh is radial, we have h 0 0 (7.6) N(0 )
rGh (x() ? x(0 )) = jx() ? x(0)j?1 @G @r [N( )
(x() ? x( )] : The last factor is O(j ? 0j2 ) for  near 0 , since x() ? x(0) is almost tangent. From the estimates (2.4) for Gh , we have j(@=@r)m Gh j  C j ? 0 j?m?1 for j ? 0j  h and O(h?m?1 ) otherwise, where m = 1 or 2. Combining these facts, we nd for the kernel K(; 0) of (7.5) that jDK(; 0)j  C j ? 0j?2 for j ? 0j  h and O(h?2 ) otherwise; here D denotes a rst derivative in  or 0. Thus the error in discretizing the kernel is of order h
h?2 for ; 0 close, and order h
j ? 0j?2 when they are at least O(h) apart. We can estimate the norm of the operator resulting from this error just as in the last part of the proof of Theorem 5.1; the

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

47

norm is bounded by Chj loghj. Replacing N(` ) with N` leads to a smaller error since they dier by O(h3). We argue similarly for the second operator J2 . We need to compare the kernel 2DGh (x() ? x(0)) with discrete values of 2 Dh Gh . (Again the change in the integral from regularizing G to Gh is O(h3).) Assuming rst that Dh is a dierence operator, the error in replacing DGh (xj ? x` ) with Dh Gh (xj ? x` ) can be estimated by writing (Dh ? D)Gh as an integral of D2 Gh . This is controlled by a bound for h
D2 Gh on an O(h) neighborhood of xj ? x` . Since 2 (j ; `) = 0, we nd that the contribution to the error in the kernel is then of order h
jj ? ` j?2 when j ; ` are apart, and order h?1 when they are close. Just as before, this estimate leads to an error in the operator with norm of order hj loghj. Next we assess the error in discretizing 2 DGh . The error in the kernel has terms of the form h
2 D2 Gh and h
DGh . For either of these we obtain the same estimates as before, using again the fact that 2 = 0 on the diagonal. This completes the treatment of J2 if Dh is a dierence operator. If not, we can reduce to the special case by comparing the discrete 2Dh Gh with one where Dh is replaced by a dierence operator. For this we note that Dh(1) ? Dh(2) has the form hBE 2 for some bounded B, and by Thm. 5.2, (5.5), the two operators dier by O(h). In summary, we have shown that the error in discretizing each operator J1, J2 is O(hj loghj), and thus Kh0 ? K is of the same order. Finally we consider Kh ? Kh0 . Each of the two terms multiplies wj by a factor. For the correction term this factor is O(h); for the other term it is a sum which is within O(h) of the rst integral in (1.5), which is zero. Thus Kh ? Kh0 is O(h), and therefore Kh ? K is O(hj loghj), as operators on L2h . Proof of Lemma 6.2. The exact equations for x; , as functions of ; t, have the

form xt = v , t = jvj2=2 ? gx3 , with v determined by v = rT  + wn , and with ( 21 I + K(x))w = f(x; ). We think of this system in the form zt = F(z), where z = (x; ) and F is a nonlocal mapping. Let H s be the Sobolev space of scalar-valued periodic functions of  with s derivatives in L2 , X s the space of such vector-valued functions, and Z s the space of pairs x 2 X s ,  2 Y s. We choose s0 > 1 and check that F is C 2 from Z s1 to Z s0 when s1 is large enough relative to s0 . First, the single and double layer potentials de ne bounded operators from H s0 to H s0+1 , and they are at least C 2 as functions of x 2 X s1 ; the x-derivatives do not change the order of the singularity. Then f : Z s1 ! H s0 is C 2. Since I +2K is a Fredholm operator on H s0 and is 1-1, it is invertible, and the inverse is a C 2 function of x 2 X s1 . It now follows easily that F : Z s1 ! Z s0 is C 2 . We need to use DF(z 0), the derivative of F at the exact, smooth solution z 0 . We presume that z 0 exists in Z s2 for some s2 larger than s1 . Finding DF(z 0 ) amounts to linearizing the equations about z 0, and this can be done as for the discrete equations in x6 but more easily. Writing the in nitesimal variation in z in the u variables of (6.45), we have DF(z 0 )u = (?cu2; u1; 0)+(0; B2u; B3 u), where  is a positive pseudodierential operator of order 1 and B2 ; B3 are bounded in H s1 by u2; u3 2 H s1 . For solutions of ut = DF(z 0 )u, we can estimate as in x6 but in high norms. Let Um = (Dm u1; Dm u1)L2 + jDm (c1=2u2)j2L2 + jDm (c1=2 u3)j2L2 where m a multi-index for -derivative D, and let U = Um , jmj  s1 . Then Ut  C1U and thus U(t)  C2U(0) for 0  t  T. (Cf. [5], p. 1282.) The discrete equations have a similar structure zh;t = Fh (zh ), with vh determined by wh satisfying ( 12 I + Kh (xh ))wh = fh (zh ). The consistency argument shows that

48

J. THOMAS BEALE

for any z = (x; ) 2 Z s1 we have ( 12 I + Kh (x))n = fh (z) + h3r31(z) + r41(z; h), where r31(z) 2 H s0 , depending smoothly on z 2 Z s1 , and r41 is uniformly O(h4 ). To invert, we use the fact that ( 12 I + Kh (x))?1 diers as an operator on L2 from ( 21 I + K(x))?1 by O(hj loghj), as shown in the proof of Lemma 6.1. Thus we can write n = ( 21 I + Kh (x))?1 fh (z) + h3 r32(z) + r42(z; h), where r32(z) = ( 12 I + K(x))?1r31(z) depends smoothly in H s0 on z as above, and r42 is bounded in L2h by Ch4 j loghj. (The log h can be removed by a more careful argument.) Similarly, for any z 2 Z s1 , we nd F(z) ? Fh (z) = h3 r3(z) + r4 (z; h), where r3 : Z s1 ! Z s0 is smooth and r4(z; h) is O(h4j loghj) in L2h . We can also assume that r3(z 0 ) 2 Z s1 since z 0 2 Z s2 . Now we can combine the above to modify z 0 . First, we can solve the linear equation t = DF(z 0 ) ? r3(z 0 ),  = 0 at t = 0, with  2 Z s1 ; this can be done in a standard way from the linear estimate above. Now set z 1 = z 0 + h3  in Z s1 . Then F(z 1 ) = F(z 0 ) + h3DF(z 0 ) + O(h6 ) in Z s0 , since F is C 2. We check that z 1 has the desired property. We have zt1 = F(z 0) + h3DF(z 0 ) ? h3 r3(z 0 ) = F(z 1) ? h3 r3(z 0 ) + O(h6) = Fh (z 1 ) ? h3 r3(z 0 ) + O(h6) + h3 r3(z 1 ) + r4(z 1 ; h). Now since z 1 ? z 0 = O(h3) in Z s1 , r3(z 1 ) ? r3(z 0 ) = O(h3 ) in Z s0 . Thus the O(h3) part cancels, and and we conclude that zt1 = Fh (z 1 ) + O(h4j loghj) in L2h . The two components of this equation give (6.3) for p < 4. The last assertion follows from combining this with (6.2). Proof of (6.27). We need to assess the dierence between X (7.7) rGh (~xj ? x~` )
N~` [(~n`
n~ j )w~` ? w~j ] h2 and the same where the computed tilde quantities are replaced with the quantities without tildes. To do this we will consider various separate dierences by adding and subtracting. (All sums are over ` 2 I). We note rst that under the hypothesis (6.6) we have, using Lemmas 6.3, 6.4 (7.8) jnjL1h  Ch?1jnjL2h  Ch?2 jnjHh?1  Ch?2jxjL2h  Ch?2 h3 = Ch and similarly for w, N. We begin with the sum X (7.9) rGh (xj ? x` )
N~` [(~n`
n~ j )w~` ? w~j ] h2 and consider several dierences which, when added to the above, will result in the corresponding sum with no tildes. The rst such dierence will be the replacement of N~` by N` . We split this into two terms, X (7.10) rGh (xj ? x`)
(N` ) [(n`
nj )w` ? wj ]h2 ; X (7.11) rGh (xj ? x` )
(N` ) [((n`
nj )w` ) ? wj ] h2 : We view (7.10) as the action on N of the kernel rGh times (n`
nj )w` ? wj . The latter factor is the restriction of a smooth function vanishing on the diagonal, except for a small error, according to the last statement in Lemma 6.2. We ignore the error, since it can be handled as below for (7.11). By Thm. 5.1 the resulting operator has order ?1 on L2h , and according to the remark before Lemma 6.3 it is bounded from Hh?1 to L2h . Thus the sum (7.10) is bounded in L2h by C jN jHh?1  C jxjL2h . Now for (7.11), the quantity in brackets is uniformly O(h), according to (7.8) and the analogue for w. We can use the boundedness of rGh on L2h , again from Thm.

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

49

5.1, to estimate this term by ChjN jL2h  C jN jHh?1  C jxjL2h . Next, to replace n~ ` by n` , we consider the dierence X (7.12) rGh (xj ? x` )
N` [w~`(n` )]
n~ j h2 : Now the operator rGh
N(` ) is of order ?1 by Thm. 5.1, where N(` ) is the exact normal. Since N` ? N(` ) = O(h3) and rGh has order 0, it follows that rGh
N` is also of order ?1. Using this fact and Lemmas 6.3, 6.4 we can estimate (7.12) by C jw(n) ~ jHh?1  C jnjHh?1  C jxjL2h . We have now replaced (7.9) with X (7.13) rGh (xj ? x` )
N` [(n`
n~ j )w~` ? w~j ]h2 : If we replace w~` , w~j with w` , wj , the resulting dierence is the rst term on the right in (6.27) plus X (7.14) (nj )
rGh (xj ? x` )
N` n` (w` )h2 : Since nj is uniformly O(h), the last is bounded in L2h by ChjwjL2h . Then if we replace n~ j by nj , the dierence is X (7.15) (nj )
w` n` rGh (xj ? x` )
N` h2 : In the sum above, we can replace w` with n (`) etc. using Lemma 6.2 with small error. The resulting sum is then within O(h) of the integral term in (6.27), according to Thm. 3.7. The term from the quadrature error is bounded in L2h by ChjnjL2h  C jxjL2h . We have now replaced the sum (7.9) with the corresponding sum without tildes, and obtained the two principal terms in (6.27). Next we consider the dierence which results from replacing x~j ; x~` with xj ; x` in the sum X (7.16) rGh (~xj ? x~` )
N` [(n`
nj )w` ? wj ]h2 : The dierence is a similar sum with  rGh instead of rGh . We write  rGh = R2(j ; `)
(xj ? x`), where R2 denotes an average of second derivatives of Gh along the line from xj ? x` to x~j ? x~` . Since the bracket term vanishes on the diagonal, Thm. 5.1 tells us that the x` part of the error is bounded in L2h by C jxjL2h . For the second part, we factor xj out of the sum; the remaining sum is bounded, uniformly in j and h, according to Lemma 5.4. Thus the second part of the error has the same bound as the rst. We are now left with the dierence X (7.17)  rGh
[N` ((n`
nj )w` ? wj )]h2 : By considering cases we can see that this is of type B(x). From  rGh we have R2 times a factor of xj or x`, and each term in the second factor has at least one dierence in either j or `. We discuss two speci c terms. For the term X X (7.18) R2
N~` (xj )(wj )h2 = (xj )(wj ) R2
N~` h2 we note that jN~` j is uniformly bounded, and wj is uniformly O(h). The last sum is O(h?1 ), uniformly in j; this can be seen from the pointwise estimates (2.4) for

50

J. THOMAS BEALE

Dm Gh . Thus the term has the form xj times a uniformly bounded function of j. Next we consider the term X (7.19) R2
N~` (xj )(n` )w~` h2 : We can estimate j(n)w~jL2h by Ch2. Then from Thm. 5.1, hR2 is a bounded operator on L2h , so that the sum omitting the xj factor is O(h) in L2h , and therefore O(1) in L1 h . This is then multiplied by xj , and the resulting term is B(x). The other terms can be treated by the same methods. Proof of (6.28). We have to nd the dierence between the sum X (7.20) ? Dh;k Gh (~xj ? x~` )X~k`
n~ j w~` jN~` jh2 : and the corresponding sum without tildes. We have written out rTh Gh , recalling (1.23), with a sum over k = 1; 2 implicit. We begin with the sum X (7.21) ? Dh;k Gh (xj ? x`)X~k`
(~nj ? n~ ` ) w~` jN~` jh2 and add dierences successively to remove the tildes. We have inserted ?n~ ` ; this does not aect the sum, since n~ ` ? X~k` . We rst replace n~ j with nj , resulting in the dierence, after a summation by parts, X  [X  w` jN` j]h2 (7.22) ?(nj )
Gh (xj ? x` )Dh;k k` plus a remainder. The remainder has (Xk` w` jN` j), which is O(h2) is L2h , acted on by Dh;k Gh . Since Gh is an operator of order ?1 on L2h , by Thm. 5.1, Dh;k Gh is bounded, so that the image is O(h2 ) in L2h and thus O(h) in L1 h . This is multiplied by nj , which is bounded in L2h by h?1jxjL2h , so that remainder term is of type B(x). As for the main term (7.22), the sum on the right is close to the integral in (6.28), after replacing w` with n(` ) etc., with O(h) quadrature error. This last error gives a term bounded by hjnjL2h or jxjL2h . Next we replace n~ ` by n` . The dierence is X (7.23) Dh;k Gh (xj ? x` )Xk`
(n` ) w` jN` jh2 plus a remainder again with (Xk` w` jN` j). As before the remainder is B(x). We can substitute in (7.23) for n from (6.12) and obtain the second term on the right in (6.28), up to a remainder which is B(x). We now have reduced (7.20) to X (7.24) ? Dh;k Gh (xj ? x` )X~k`
(nj ? n` ) w~` jN~` jh2 : Now since nj ? n` vanishes at j = `, Thm. 5.2 tells us that the operator with kernel (nj ? n` )Dh;k Gh gains a derivative. Thus, to replace X~k` jN~` j, we estimate w(X ~ k jN j) in Hh?1 by jxjL2h . The corresponding sum is then bounded by jxjL2h . To nish with (7.20), we replace w~` with w`, obtaining as the dierence the rst term on the right of (6.28). Next we consider the dierence X (7.25)  [Dh;k Gh (xj ? x` )]Xk`
nj w` jN` jh2 X  [X  w` jN` j]h2 : = [Gh (xj ? x` )]
nj Dh;k k`  We write Gh as R1(j ; `)
(xj ? x`), where R1 is an average of DGh from xj ? x` to x~j ? x~` . For the x` part, we conclude from the boundedness of DGh on

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

51

L2h that the sum is bounded by jxjL2h . For the second part, we bring xj outside the sum. The remaining sum is uniformly bounded according to Lemma 5.4, and this term is B(x). We are left with more nonlinear terms, which can be handled as at the end of (6.27), after summing by parts to remove Dh;k from Gh as above. Proof of (6.29). The dierence has the form (7.26) (n
) = (n)
 + n
 + (n)
() Since nj is uniformly O(h), by (6.6), it will be evident that the last term can be neglected after we have treated the second term. The sums in the rst term are within O(h) of the corresponding integrals, according to Thm. 3.7. (We use summation by parts.) The resulting error term is B(hn) = B(x). The integrals can be converted to the integral in (6.29) by reversing the steps which led to (1.7). Now our task is to evaluate the second term in (7.26). We treat the rst sum in (6.29). As before we begin with X ~ `  (r~ Th ~` ? r~ Th ~j )h2 (7.27) nj
[Dh Gh ; X] P where Gh = Gh (xj ? x` ) and r~ Th f = k=1;2(Dh;k f)X~k . We remove the tildes in steps; we rst replace r~ Th ~` with rTh `. The dierence is X (7.28) nj
[Dh Gh ; X]`  (rTh ` )h2 plus a remainder with products of 's. This remainder is B(x); to see this, we note that Dh Gh is a bounded operator on L2h and products such as (X1` )(Dh;2 ` ) are bounded in L2h by jxjL2h , using (6.6). In (7.28) we can substitute for (rTh ` ) from (6.14) obtaining, up to a negligible remainder, X

nj
[Dh Gh ; X]`  rTh (` ? rTh `
x` )h2 X ?
+ ([Dh Gh ; X]`  n` )
nj rTh (n`
x` )
rTh ` h2 : In the rst sum we substitute nj = n` + (nj ? n` ) and split the sum into two parts. The n` -part can be converted to the rst term on the right in (6.29) using a vector identity as in (1.6). For the remaining part in the rst sum of (7.29), we have a smooth factor which vanishes at ` = j, and we can use Thm. 5.2 to show that the sum is B(x) + B(). As for the second sum in (7.29), we note that the cross product is tangential at x`, so that the scalar product with nj is zero when ` = j. Thus we can use Thm. 5.2 again to bound this term by jxjL2h . Next we replace r~ Th ~j in (7.27) with rTh j . The dierence is X (7.30) nj
(rTh j )  [Dh Gh ; X]` h2 except for an error nonlinear in the 's which is again negligible. The sum in (7.30) is within O(h) of an integral which is zero; see (1.5),(1.8),(1.28). The term with the quadrature error is of order hj rTh jL2h , which in view of (6.14) is B() + B(x). We have now reduced (7.27) to X ~ `  (rTh ` ? rTh j )h2 : (7.31) nj
[Dh Gh ; X] The last factor vanishes when ` = j, and we can therefore use Thm. 5.2 to replace X~` with X` , with a dierence B(x), thus completing the treatment of (7.27). (7.29)

52

J. THOMAS BEALE

The most important remaining term from the rst sum on the left in (6.29) is X (7.32) nj
[(Dh Gh ); X]`  (rTh ` ? rTh j )h2 : This can be treated similarly to the corresponding term in (6.28). The other terms are nonlinear and can be handled like earlier cases. The variation of the second sum can be treated by techniques as for (6.27) but is simpler. We need to use the fact that rTh j is multiplied by a discretization of an integral which is zero by (1.5). Proof of Lemma 6.7. The operator K~ h diers from Kh by an operator which is O(h). This can be shown by estimates related to (6.27), (6.28) but simpler. By the remarks in the proof of Lemma 6.1, 2K~ h then has spectral radius less than one, and it follows that the discrete integral equation can be solved by simple iteration. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22.

References C. Anderson and C. Greengard, On vortex methods, SIAM J. Numer. Anal. 22 (1985), 413-40. V. K. Andreev, Stability of unsteady motions of a uid with a free boundary, VO Nauka, Novosibirsk, 1992 (in Russian). G. Baker, D. Meiron, and S. Orszag, Generalized vortex methods for free-surface ow problems, J. Fluid Mech. 123 (1982), 477-501. G. Baker, D. Meiron, and S. Orszag, Boundary integral methods for axisymmetric and threedimensional Rayleigh-Taylor instability problems, Physica D 12 (1984), 19-31. J. T. Beale, T. Y. Hou and J. S. Lowengrub, Growth rates for the linearized motion of uid interfaces away from equilibrium, Comm. Pure Appl. Math. 46 (1993), 1269-1301. J. T. Beale, T. Y. Hou and J. S. Lowengrub, Convergence of a boundary integral method for water waves, SIAM J. Numer. Anal. 33 (1996), 1797-1843. J. T. Beale and A. Majda Vortex methods, I: Convergence in three dimensions, Math. Comp. 39 (1982), 1-27. J. T. Beale and A. Majda, High order accurate vortex methods with explicit velocity kernels, J. Comput. Phys. 58 (1985), 188-208. J. Broeze, E. F. G. Van Daalen, and P. J. Zandbergen, A three-dimensional panel method for nonlinear free surface waves on vector computers, Comput. Mech. 13 (1993), 12-28. D. Colton and R. Kress, Integral Equation Methods in Scattering Theory, Wiley, New York, 1983. G. H. Cottet and P. A. Raviart, Particle methods for the one-dimensional Vlasov-Poisson equations, SIAM J. Numer. Anal. 21 (1984), 52-76. J. Goodman, T. Y. Hou and J. Lowengrub, Convergence of the point vortex method for the 2-D Euler equations, Comm. Pure Appl. Math. 43 (1990), 415-30. I. S. Gradshteyn and I. M. Ryzhik, Table of Integrals, Series, and Products, Academic Press, New York, 1980. L. Greengard and V. Rokhlin, A new version of the fast multipole method for the Laplace equation in three dimensions, Acta Numer. 6 (1997), 229-69. B. Gustafsson, H. O, Kreiss, and J. Oliger, Time dependent problems and dierence methods, Wiley, New York, 1995. O. Hald, The convergence of vortex methods, II, SIAM J. Numer. Anal. 16 (1979), 726-755. D. Haroldsen and D. Meiron, Numerical calculation of three dimensional interfacial potential

ows using the point vortex method, SIAM J. Sci. Comput. 20 (1999), 648-83. T. Y. Hou, Z. Teng, and P. Zhang, Well-posedness of linearized motion for 3-D water waves far from equilibrium, Commun. P.D.E. 21 (1996), 1551-85. T. Y. Hou and P. Zhang, Stability of a boundary integral method for 3-D water waves, submitted to SIAM J. Numer. Anal. R. E. Kleinman and G. F. Roach, Boundary integral equations for the three-dimensional Helmholtz equation, SIAM Review 16 (1974), 214-36. M. S. Longuet-Higgins and E. D. Cokelet, The deformation of steep surface waves on water, I. A numerical method of computation, Proc. Roy. Soc. London A 350 (1976), 1-26. J. Lowengrub, M. Shelley, and B. Merriman, High-order and ecient methods for the vorticity formulation of the Euler equations, SIAM J. Sci. Comput. 14 (1993), 1107-42.

A BOUNDARY INTEGRAL METHOD FOR 3-D WATER WAVES

53

23. J. N. Lyness, An error functional expansion for N-dimensional quadrature with an integrand function singular at a point, Math. Comp. 30 (1976), 1-23. 24. J. N. Lyness, A survey of numerical cubature over triangles, in Proc. Symp. Appl. Math. 48, A.M.S., Providence, 1994. 25. A. Majda, J. McDonough and S. Osher, The Fourier method for nonsmooth initial data, Math. Comp. 32 (1978), 1041-81. 26. R. D. Richtmyer and K. W. Morton, Dierence Methods for Initial Value Problems, 2nd ed., Wiley, New York, 1967. 27. J. E. Romate, The numerical simulation of nonlinear gravity waves, Engrng. Analysis Bdry. Elts., 7 (1990), 156-66. 28. J. E. Romate and P. J. Zandbergen Boundary integral equation formulations for free-surface

ow problems in two and three dimensions, Comput. Mech. 4 (1989), 267-82. 29. J. Strain, Fast potential theory. II. Layer potentials and discrete sums, J. Comput Phys. 99 (1992), 251-70. 30. M. Taylor, Pseudodierential Operators, Princeton Univ. Press, Princeton, NJ, 1981. 31. M. Taylor, Partial Dierential Equations, Springer, New York, 1996. 32. W. Tsai and D. Yue, Computations of nonlinear free-surface ows, Ann. Rev. Fluid Mech. 28 (1996), 249-78. 33. T. Vinje and P. Brevig, Numerical simulation of breaking waves, Adv. Water Resources 4 (1981), 77-82. 34. Sijue Wu, Well-posedness in Sobolev spaces of the full water wave problem in 2-D, Invent. Math. 130 (1997), 39-72. 35. Sijue Wu, Well-posedness in Sobolev spaces of the full water wave problem in 3-D, J. Amer. Math. Soc. 12 (1999), 445-95. Department of Mathematics, Duke University, Durham, NC 27708-0320

E-mail address :

[email protected]