Multiresolution representation of operators with boundary conditions ...

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS ON SIMPLE DOMAINS GREGORY BEYLKIN* , GEORGE FANN+ , ROBERT J. HARRISON+, CHRISTOPHER KURCZ* , LUCAS MONZÓN*

Abstract. We develop a multiresolution representation of a class of integral operators satisfying boundary conditions on simple domains in order to construct fast algorithms for their application. We also elucidate some delicate theoretical issues related to the construction of periodic Green’s functions for Poisson’s equation. By applying the method of images to the non-standard form of the free space operator, we obtain lattice sums that converge absolutely on all scales, except possibly on the coarsest scale. On the coarsest scale the lattice sums may be only conditionally convergent and, thus, allow for some freedom in their definition. We use the limit of square partial sums as a definition of the limit and obtain a systematic, simple approach to the construction (in any dimension) of periodized operators with sparse non-standard forms. We illustrate the results on several examples in dimensions one and three: the Hilbert transform, the projector on divergence free functions, the non-oscillatory Helmholtz Green’s function and the Poisson operator. Remarkably, the limit of square partial sums yields a periodic Poisson Green’s function which is not a convolution. Using a short sum of decaying Gaussians to approximate periodic Green’s functions, we arrive at fast algorithms for their application. We further show that the results obtained for operators with periodic boundary conditions extend to operators with Dirichlet, Neumann, or mixed boundary conditions.

1. Introduction The primary goal of this paper is to develop a multiresolution representation of a class of integral operators satisfying boundary conditions on simple domains and construct fast algorithms for their application. As a practical consequence of our approach, we show that a minor modification of the fast algorithms for free space operators in [24, 9, 6], yields a fast algorithm (of the same complexity) for the operator satisfying boundary conditions. Another goal of this paper is to elucidate some delicate theoretical issues related to the method of images for the construction of periodic Green’s functions for Poisson’s equation. Indeed, due to the slow decay of the Poisson’s kernel, the solution of the periodic problem is not unique and, in fact, several physically meaningful periodic Green’s functions have been discussed in the literature (over a long period of time). Within our approach, these Green’s functions are easy to describe as particular choices of just a few parameters in the construction. Date: October 26, 2011. Key words and phrases. multiresolution, non-standard form, projector on divergence free functions, Poisson Green’s function, non-oscillatory Helmholtz Green’s function, Hilbert transform, periodic boundary conditions, separated representations. * This research was partially supported by DOE/ORNL grant 4000038129, and NSF grants DMS-0612358 and DMS-100995. + RJH gratefully acknowledges support from NSF OCI 0904972. GF gratefully acknowledges support from Applied Mathematics Program of the Office of Advanced Scientific Computing Research, U.S. Department of Energy. His work was performed at the Oak Ridge National Laboratory, which is managed by UT-Battelle, LLC under Contract No. De-AC05-00OR22725.

1

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

2

In our approach we apply the method of images not to the free space operator itself but to its nonstandard form. The non-standard form splits the action of the operator to an infinite set of scales and, for appropriate classes of operators, yields a sparse representation [7]. For operators with kernels whose partial derivatives decay faster than the kernel itself (e.g., the Calderon-Zygmund operators), the non-standard form is sparse on all scales, except for the coarsest scale. We use the rapid decay of the coefficients of the non-standard form to construct its periodized version and to show that, on all scales except possibly the coarsest scale, the lattice sums converge absolutely and require no further analysis. On the coarsest scale, for some of the coefficients, the lattice sums may be only conditionally convergent and, thus, allow for some freedom in their definition. For such coefficients a summation convention needs to be specified and we choose the limit of square partial sums for that purpose. In this way, we obtain a systematic, simple approach to the construction (in any dimension) of periodized operators with sparse non-standard forms. We illustrate the results on several examples in dimensions one and three: the Hilbert transform, the projector on divergence free functions (the so-called Leray projector), the non-oscillatory Helmholtz Green’s function and the Poisson operator. The Poisson Green’s function appears in many fields including electrostatics, material sciences, and molecular dynamics (see e.g. [18, 27, 31]). The standard method of images when applied directly to the free space kernel yields only a formal result that requires interpretation, a key topic in lattice sum literature. As it turns out, the periodic Poisson Green’s function is non-unique which explains the appearance of several versions in the literature (see e.g. [21, 26] for a review). An early seminal contribution was made by P. Ewald [20], although the history of lattice sums starts earlier and we refer to [21] for a historical overview and results prior to 1980. Due to the slow decay of the Poisson kernel, kx − yk−1 , the analysis of its periodization turns out to be more delicate than for (even slightly) faster decaying kernels. Similar difficulties arise in other periodic problems with operators exhibiting the same rate of decay, e.g., Stokes operator recently considered in [28]. For the Poisson kernel, our approach identifies several specific components of the periodized non-standard form which converge only conditionally and, moreover, are not limits of the corresponding components of, e.g., kx − yk−1 e−µkx−yk as µ → 0 or other possible operator limits. As a peculiar consequence, the limit of square partial sums yields a Green’s function which is not a convolution, even though it may be natural to expect the method of images, according to its formal form, to always produce a convolution kernel. As a consequence, using such Green’s function to solve the Poisson equation yields solutions which are not necessarily mean-free. Our algorithms approximate the operator kernel via a separated representation given by a short linear combination of decaying Gaussians with positive exponents and coefficients, which immediately reduces the computational cost and yields a non-standard form with elements given by one dimensional sums. As a result, for any finite accuracy, we obtain an efficient separated representation in any dimension d ≥ 2 and associated fast algorithms. This type of approximation via Gaussians (see e.g. [11, 12, 13, 14]) has been successfully used in [24, 9, 6] to construct fast and accurate algorithms for applying free space convolution kernels for any user supplied finite accuracy. Using the non-standard form of free space operators, we show that, on simple domains, the periodized non-standard form also yields fast and accurate algorithms for applying periodic operators as well as for applying operators satisfying Dirichlet or Neumann boundary conditions. We also note that the Fast Multipole Method (FMM) [23, 16] may also be used to apply such periodic operators. We limit our presentation to the non-standard forms of weakly singular or singular operators. We note that non-standard forms may also be constructed for hyper-singular operators [8]. However, periodization of such operators does not present a challenge due to the rapid decay of their kernels away from singularities and we do not discuss them in this paper. In order to limit the size of the paper, we do not present numerical results. We note, however, that the speed of algorithms

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

3

for Green’s functions with boundary conditions is essentially the same as that for the free space case. Indeed, we show that the operators effectively coincide on the wavelet scales which are those dominating the computational cost. We start in Section 2 by introducing the non-standard form for convolution operators in dimension d = 1 using multiwavelet bases [1, 2, 3]. In this case only one term may require an appropriate interpretation and we illustrate this using the Hilbert transform as an example. In Section 3 we construct the non-standard form in dimension d = 3 for operators with periodic boundary conditions. As examples, we then analyze the projector on divergence free functions, the non-oscillatory Helmholtz Green’s function and, in Section 4, the Poisson Green’s function. In Section 5 we describe a fast algorithm for applying these operators using separated representations. In Section 6, we construct Green’s functions which incorporate Dirichlet, Neumann, or mixed boundary conditions on simple domains. Finally, we provide some closing remarks in Section 7 and collect most proofs in the appendix. 2. Preliminaries 2.1. Multiresolution and non-standard form. In this section we review a multiresolution approach for representing and applying operators in one dimension. Since we use multiwavelets as the underlying basis for the multiresolution representation, we briefly describe their properties (see also [1, 3, 9, 6]). We then turn to the non-standard form of operators in multiwavelet bases and describe a class of operators which becomes effectively sparse in this representation (see also [7, 6]). We then construct an operator with periodic boundary conditions by applying the method of images to the components of the non-standard form and illustrate the result with the Hilbert transform. The notation used below deviates slightly from usual wavelet notation, however, its introduction facilitates the higher dimensional description presented in later sections. m 2.1.1. Multiwavelets. Let P[a,b] denote the space of polynomials of degree less than m restricted to the interval [a, b]. Let us define subspaces [ 2 m Vj = P[2 −j l,2−j (l+1)] ⊂ L (R) l∈Z

for j ∈ N, where N denotes all non-negative integers. These subspaces are nested V0 ⊂ V1 ⊂ · · · ⊂ Vj ⊂ · · ·

S j;l 2 and ∞ j=0 Vj = L (R). We select scaling functions to form an orthonormal basis of Vj , ψi;0 (x) = 2j/2 ψi;0 (2j x − l), j ∈ N, l ∈ Z, where (√ 2i + 1Pi (2x − 1), x ∈ [0, 1] (1) ψi;0 (x) = , i ∈ {0, . . . , m − 1}, 0, otherwise and Pi are the i-th order Legendre polynomials. We will need the cross-correlation functions of the scaling functions, Z ψi;0 (x + y)ψi′ ;0 (y)dy, (2) Φii′ (x) = R

i, i′

∈ {0, . . . , m − 1}. Due to orthogonality of the scaling functions where supp(Φ ) ⊂ [−1, 1] for in (1), these functions have vanishing moments (see [9, §2.2]), Z Φii′ (x)xk dx = 0 for i + i′ ≥ 1, and 0 ≤ k ≤ i + i′ − 1. (3) ii′

R

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

4

We define the wavelet subspaces Wj as Wj ⊕ Vj = Vj+1 , so that Vj+1 = V0 ⊕ W0 ⊕ · · · ⊕ Wj . j;l We denote the multiwavelets, an orthonormal basis of Wj , as ψi;1 for i ∈ {0, . . . , m − 1} and l ∈ Z. We do not need an explicit expression for the multiwavelets and only use their vanishing moments property, Z j,l ψi;1 (x)xk dx = 0 for i, k = 0, . . . , m − 1, l ∈ Z, and j ∈ N, (4) R

which follows from orthogonality of the subspaces Wj and Vj . Also we need the cross-correlation functions of the wavelets, Z ψi;s (x + y)ψi′ ;s′ (y)dy, (5) Φii′ ;ss′ (x) = R

ss′

i, i′

where = 11, 10, 01 and ∈ {0, . . . , m − 1}. In this notation Φii′ ;00 = Φii′ in (2) are the cross-correlations of the scaling functions. In L2 (Rd ) we use the tensor-product basis formed by products of multiwavelets and scaling functions ′ j;l from the same scale. For example, in dimension d = 2 the basis for Wj is given by ψi;1 (x)ψij;l ′ ;1 (y), ′

′

′

j;l j;l j;l j;l j;l ψi;1 (x)ψij;l ′ ;0 (y), and ψi;0 (x)ψi′ ;1 (y), whereas the basis for Vj is given by ψi;0 (x)ψi′ ;0 (y), where (x, y) ∈ R2 , i, i′ ∈ {0, . . . , m − 1}, and l, l′ ∈ Z.

Since multiwavelet bases include discontinuous basis functions (like those of the Haar basis), using them as the underlying basis for the non-standard form (see below) limits our discussion to, at most, singular operators. If we were to extend our approach to include hyper-singular operators, it would be necessary to use sufficiently smooth wavelets as the underlying basis [8].

2.1.2. Non-standard form. The non-standard form of an operator T is based on the telescopic series representation T

= P0 T P0 +

∞ X j=1

(6)

= P0 T P0 +

∞ X

(Pj T Pj − Pj−1 T Pj−1 ) (Qj T Qj + Qj T Pj + Pj T Qj ) ,

j=0

where Qj and Pj are the orthogonal projectors, Qj : L2 (Rd ) → Wj and Pj : L2 (Rd ) → Vj and Pj+1 = Pj + Qj . For the purposes of this paper we distinguish the wavelet part of the non-standard form, namely, (7)

Twavelet = {Qj T Qj , Qj T Pj , Pj T Qj }j∈N

from the scaling part of the non-standard form, (8)

Tscaling = P0 T P0 .

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

5

2.1.3. Example in one dimension. Let K be the kernel of the convolution operator Z K(x − y)f (y)dy. (9) (T f ) (x) = R

The elements of the scaling part of the non-standard form are computed as the projection of the kernel onto the scaling functions, Z Z Z ′ 0;l 0;l′ (10) K(x)Φii′ (x + l′ − l)dx (y) dydx = K(x − y)ψ (x)ψ = Tii0;l−l ′ ;00 i;0 i′ ;0 R ZR R ′ K(x + l − l )Φii′ (x)dx, = R

l, l′

i, i′

for ∈ Z and ∈ {0, . . . , m − 1}. Similarly, the elements of the wavelet part of the non-standard form are computed as the projection of the kernel onto a multiwavelet functions, Z Z Z ′ j;l j;l′ (11) K(x)Φii′ ;ss′ (2j x + l′ − l)dx (y) dydx = K(x − y)ψ (x)ψ = Tiij;l−l ′ ;ss′ i;s i′ ;s′ R R R Z = 2−j K(2−j (x + l − l′ ))Φii′ ;ss′ (x)dx, R

for j ∈ N, moments.

l, l′

∈ Z,

ss′

= 11, 10, 01, and i, i′ ∈ {0, . . . , m − 1}, where m is the number of vanishing

Remark. Limiting our analysis to singular operators assures existence of the elements of the nonstandard form in (10) and (11). Definition 1. We say that an operator T is integral-defined if the elements of its non-standard form (10) and (11) are given by either absolutely or conditionally convergent integrals. Examples of integral-defined operators include weakly singular and singular Calderon-Zygmund operators, and various classes of pseudo-differential operators, see [8]. Proposition 2. Let T be an integral-defined operator with a convolution kernel K ∈ C m (R\{0}), m ≥ 1, satisfying, cα (12) |∂xα K(x)| ≤ for cα > 0, 0 ≤ α ≤ m and β ≥ 1. |x|α+β

Then, represented in a multiwavelet basis with m vanishing moments, j ∈ N, l, l′ ∈ Z, and i, i′ = 0, . . . , m − 1, the elements of the wavelet part of the non-standard form satisfy
− min{m,m}−β j;l−l′ , Tii′ ,ss′ ≤ Cj 1 + l − l′ ′

′ where Tiij;l−l ′ ;ss′ is given by (11), ss 6= 00, with constants Cj > 0 that depend on the scale j but not on l or l′ . The elements of its scaling part satisfy
− min{i+i′ ,m}−β 0;l−l′ , Tii′ ,00 ≤ C0 1 + l − l′ ′

′ where Tii0;l−l ′ ;00 is given by (10) and the constant C0 does not depend on l or l .

See Appendix 8.2 for the proof. Proposition 2 states that the non-standard form of operators satisfying (12) are effectively sparse. Indeed, the operator norm of the difference between the infinite m × m block-Toeplitz matrices n o ′ n ol,l′ ∈Z n ol,l′ ∈Z ′ l,l ∈Z j;l−l′ j;l−l′ Qj T Qj = Tiij;l−l , Q T P = T , and P T Q = T ′ ;11 j j j j ii′ ;10 ii′ ;01 ′ ′ ′ i,i ∈{0,...,m−1}

i,i ∈{0,...,m−1}

i,i ∈{0,...,m−1}

and their banded versions (obtained by setting to zero blocks with |l − l′ | ≥ b) decay rapidly at least as b− min{m,m}−β , where b is the width of the band. Hence, for any finite but arbitrary accuracy,

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

6

the entries outside the band may be discarded resulting in a representation of the operator in terms of banded matrices and, therefore, yielding a fast algorithm for its application (see e.g. [7]). 2.2. Operators with periodic boundary conditions. Given a convolution operator T of the form (9), the method of images is the standard approach to construct an associated operator T satisfying a periodic boundary condition. Specifically, # Z 1 "X K(x − y + n) f (y)dy, (13) T f (x) = 0

n∈Z

where (T f ) (x) = (T f ) (x + 1) for x ∈ [0, 1]. However, the sum in (13) may require further analysis since it may diverge or converge only conditionally. Instead of considering (13), we first construct the non-standard form of the free space operator (9) and then apply the method of images to the elements of the non-standard form. By linearity, given the elements (10)-(11) of the non-standard form of the free space operator (7)- (8), we may construct lattice sums on each scale separately. As the method of images for the elements of the wavelet part of the non-standard form, we define the periodized operator on scale j ∈ N as X j;l−l′ +2j n ′ Tii′ ;ss′ , (14) Tiij;l−l ′ ;ss′ = n∈Z

for l, l′ ∈ {0, . . . , 2j −1}, ss′ = 11, 10, 01, and i, i′ ∈ {0, . . . , m−1}. In this way, restricting indices l, l′ to the set {0, . . . , 2j −1} in (14), limits the integration in (11) to a unit interval while the summation over index n achieves the periodicity. If the kernel K satisfies the assumptions of Proposition 2, we have Cj j;l−l′ X j;l−l′ +2j n X , Tii′ ;ss′ ≤ Tii′ ;ss′ ≤ ′ + 2j n|)min{m,m}+β (1 + |l − l n∈Z n∈Z and, hence, the sum in (14) converges absolutely for any choice of multiwavelet basis with m ≥ 1 vanishing moments. We note that the sum in (14) formally corresponds to the projection of the periodized kernel on the wavelet subspaces, Z 1Z 1X X j;l−l′ +2j n ′ j;l K(x − y + n)ψij;l Tii′ ;ss′ = ′ ;s′ (y)ψi;s (x)dydx 0

n∈Z

=

Z

0 n∈Z

1

X

−1 n∈Z

K(x + n)Φii′ ;ss′ (2j x + l′ − l)dx,

but while the series on the left hand side is absolutely convergent for ss′ 6= 00, the sum on the right hand side may not converge. For ss′ = 00, the elements of the scaling part of the non-standard form satisfy C0 0;n , Tii′ ;00 ≤ min{i+i′ ,m}+β (1 + |n|)

and, thus, unless i + i′ = 0, (15)

Tii0;0 ′ ;00 =

X

Tii0;n ′ ;00 ,

n∈Z

is absolutely convergent for any choice of multiwavelet basis with m ≥ 1 vanishing moments. However, for i = i′ = 0, the absolute convergence is not guaranteed and we choose a symmetric summation convention, namely, (16)

0;0 T00;00

= lim

N →∞

N X

n=−N

0;n T00;00 .

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

7

Remark 3. For kernels satisfying the assumptions of Proposition 2, on fine scales (j ≫ 1) and for fixed l, l′ , we have X j;l−l′ +2j n ′ Tii′ ;ss′ ∼ Tiij;l−l ′ ;ss′ n∈Z

due to the rapid decay of the terms in the sum. Thus, for any given accuracy ǫ > 0, on a sufficiently fine scale the norm of the difference between the non-standard form of the free space and the periodized operators is less than ǫ.

2.2.1. Example. Assuming that the non-standard form of the Hilbert transform (a singular operator with kernel K(x) = 1/π p.v. 1/x) is available in the multiwavelet basis with m ≥ 1, we consider its periodic version, Z 1X Z 1 1 1 p.v. f (y)dy = p.v. cot π(x − y) f (y)dy, π x−y+n 0 0 n∈Z

where f ∈

L2 [0, 1].

Since the kernel K satisfies Proposition 2, the Hilbert transform is effectively sparse in the nonstandard form. Furthermore, all elements of the wavelet part of the non-standard form of the Hilbert transform with periodic boundary conditions converge absolutely. In fact, due to rapid convergence of the series, we may compute them directly via (14). For the scaling part of the non-standard form, all elements in (15) converge absolutely except for 0;0 0;0 T00;00 . Let us show that T00;00 = 0 according to the definition (16). Indeed, we have 0;0 T00;00

=

0;0 T00;00

where 0;n T00;00

N   X 0;n 0;−n T00;00 + T00;00 , + lim N →∞

1 = p.v. π

n=1

Z

1 −1

Φ00 (x) dx. x+n

0;0 Seeing that Φ00 (x) = 1 − |x| is an even function, it follows that T00;00 = 0 due to parity. Also, for n 6= 0, we have  Z  Z 1 1 1 2x 1 1 1 0;n 0;−n T00;00 + T00;00 = + Φ00 (x)dx = 0, Φ00 (x)dx = 2 π −1 x + n x − n π −1 x − n2

where the integrals are well defined since Φ00 (±1) = 0. In this example the elements of the non-standard form coincide with those obtained using the kernel p.v. cot(πx) [6]. 3. Periodization of the non-standard form in three dimensions In this section we develop the non-standard form for operators in dimension d = 3. As in dimension d = 1, we construct the operator with periodic boundary conditions by applying the method of images to the non-standard form of the free space operator. We demonstrate that, as in dimension d = 1, all elements of the wavelet part of the non-standard form and nearly all elements of its scaling part converge absolutely. With several representative examples, we illustrate how to analyze the remaining elements of the scaling part of the non-standard form. In what follows, we denote the standard vector p-norm by kxkp .

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

8

3.1. Non-standard form in dimension three. Let us consider integral operators given by a convolution kernel in dimension d = 3, Z K(x − y)f (y)dy (17) (T f ) (x) = R3

R3 .

for x, y ∈ The basis functions (both scaling and multiwavelet) are the tensor product of the one-dimensional basis functions described in Section 2.1.1 and are denoted as j;l1 j;l2 j;l3 Ψj;l i;s (x) = ψi1 ;s1 (x1 )ψi2 ;s2 (x2 )ψi3 ;s3 (x3 ),

(18)

where x = (x1 , x2 , x3 ), j ∈ N, l = (l1 , l2 , l3 ) ∈ Z3 , i = (i1 , i2 , i3 ) ∈ {0, . . . , m−1}3 and s = (s1 , s2 , s3 ) j;l with s1 , s2 , s3 = 0 or 1. Thus, in this notation, the scaling functions correspond to Ψi;0 . We also use the cross-correlation functions of the wavelets,

(19)

Φii′ ;ss′ (x) =

Z

R3

0;0 Ψ0;0 i;s (x + y)Ψi′ ;s′ (y)dy,

for ss′ 6= 00. Since most of our analysis deals with the cross-correlations of the scaling functions, instead of denoting them as Φii′ ;00 , we simplify their notation as Φii′ , (20)

Φii′ (x) = Φi1 i′1 (x1 )Φi2 i′2 (x2 )Φi3 i′3 (x3 ),

where Φii′ are one dimensional cross-correlations of the scaling functions in (2). The elements of the wavelet part of the non-standard form of the operator in (17) are given by Z Z Z j;l−l′ j;l′ j;l K(x)Φii′ ;ss′ (2j x + l′ − l)dx, (21) Tii′ ;ss′ = K(x − y)Ψi′ ;s′ (y)Ψi;s (x)dydx = R3

ss′

for 6= 00, j ∈ N, given by (22)

′ Tii0;l−l ′ ;00

= =

R3 l, l′ ∈

Z

Z

R3

R3

R3

Z3 ,

Z

R3

and

i, i′

K(x −

= {0, . . . , m −

1}3 ,

while the elements of the scaling part are

′ 0;l y)Ψ0;l i′ ;0 (y)Ψi;0 (x)dydx

=

Z

R3

K(x)Φii′ (x + l′ − l)dx

K(x + l − l′ )Φii′ (x)dx,

for l, l′ ∈ Z3 and i, i′ = {0, . . . , m − 1}3 . We have an extension of Proposition 2: Proposition 4. Let T be an integral-defined operator (i.e., (21) and (22) are either absolutely or conditionally convergent) with convolution kernel K ∈ C m (R3 \{0}), m ≥ 3, satisfying, −|α|−β

|D α K(x)| ≤ cα kxk2

(23) where

Dα

∂ |α| /∂xα1 1 ∂xα2 2 ∂xα3 3 ,

=

for

cα > 0, 0 ≤ |α| ≤ m and

α = (α1 , α2 , α3 ) ∈

N3

β ≥ 1,

and |α| = α1 + α2 + α3 .

Then, represented in a multiwavelet basis with m vanishing moments, j ∈ N, l, l′ ∈ Z3 and i, i′ = {0, . . . , m − 1}3 , the elements of the wavelet part of the non-standard form satisfy
− min{m,m}−β j;l−l′ , Tii′ ;ss′ ≤ Cj 1 + kl − l′ k2 ′

′ where Tiij;l−l ′ ;ss′ is given by (21), ss 6= 00, with constants Cj ≥ 0 that depend on the scale j but not on l, l′ . The elements of the scaling part satisfy
− min{|i+i′ |,m}−β 0;l−l′ , Tii′ ;00 ≤ C0 1 + kl − l′ k2 ′

where Tii0;l−l ′ ;00 is given by (22).

The proof of Proposition 4 is similar to that for Proposition 2 and is presented in Appendix 8.3.

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

9

Remark. We note that the estimates of the rate of decay in Proposition 4 allow us to set to zero all blocks with kl − l′ k2 ≥ b, where b is some parameter chosen according to the desired accuracy. We show in Section 5 that separated representation of operators allows us to use ordinary banded matrices to take advantage of this property. Thus, the bandwidth parameter b should not be confused with the bandwidth of a matrix organized in a lexicographical order to represent a multidimensional operator. 3.2. Operators with periodic boundary conditions in dimension three. Using the same approach as in dimension d = 1, we apply the method of images to the non-standard form of the free space operator to construct the operator satisfying the periodic boundary condition. As before, the wavelet part elements of the non-standard form are given by X j;l−l′+2j n ′ Tii′ ;ss′ , (24) Tiij;l−l ′ ;ss′ = n∈Z3

for j ∈ N, l, l′ ∈ {0, . . . , 2j − 1}3 , ss′ 6= 00, and i, i′ ∈ {0, . . . , m − 1}3 , and we assume that the kernel satisfies the assumptions on Proposition 4. Since X j;l−l′+2j n X Cj j;l−l′ , ≤ Tii′ ;ss′ Tii′ ;ss′ ≤ ′ j nk )min{m,m}+β (1 + kl − l + 2 3 3 2 n∈Z n∈Z

the sum in (24) converges absolutely for any choice of multiwavelet basis with m > 3 − β vanishing moments. From now on, we assume m > 3 − β and with this condition all elements of the wavelet part of the non-standard form are well defined. For ss′ = 00, the elements of the scaling part of the non-standard form satisfy C0 0;n . Tii′ ;00 ≤ min{|i+i′ |,m}+β (1 + knk2 ) and, thus, for |i + i′ | > 3 − β,

(25)

Tii0;0 ′ ;00 =

X

Tii0;n ′ ;00

n∈Z3

is absolutely convergent for any choice of multiwavelet basis with m > 3 − β vanishing moments. For |i + i′ | ≤ 3 − β, we select a particular summation convention, the so-called square partial sums, X (26) Tii0;0 Tii0;n ′ ;00 = lim ′ ;00 . N →∞

knk∞ ≤N

Next we construct the non-standard form for several operators with periodic boundary conditions. We start with the projector on divergence free vector functions since its kernel decays relatively fast making the analysis simpler. Later, we consider the Poisson operator in free space and construct all possible operators with periodic boundary conditions consistent with its free space version. Within our approach it is immediate how to identify the few elements of the non-standard form responsible for this lack of uniqueness. This non-uniqueness is due to the slow decay of the free space Poisson kernel and, for the particular example we choose to present, it leads to a periodic operator which is not a convolution. 3.3. Projector on divergence free functions with periodic boundary conditions. The projector on divergence free vector functions (the so-called Leray projector) is given by the matrix of convolution kernels,   3xι xι′ διι′ 1 − , (27) Pιι′ (x) = διι′ δ(x) − 4π kxk32 kxk52

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

10

where ι, ι′ = 1, 2, 3 and διι′ denotes the Kronecker delta function (see e.g. [17] for more details). This operator may be obtained using the Riesz transform, see the derivation in, e.g., [25]. Observing that the first term in (27) is the identity operator (if ι = ι′ ), it is sufficient to consider the non-standard forms of the free space operators [9],  Z  διι′ 1 3xι xι′ ′ p.v. (28) Tιι′ f (x) = 3 − kx − yk5 f (y)dy ι, ι = 1, 2, 3, 4π 3 kx − yk R 2 2 and use them to construct the periodized non-standard form. Since operators in (28) satisfy Proposition 4 with β = 3, all elements of the wavelet part of the periodized non-standard form converge absolutely for any multiwavelet basis (m ≥ 1). We have Proposition 5. Let us consider the non-standard form of operators Tιι′ (28) in a multiwavelet basis. Then (i) The elements (24) of the wavelet part of the periodized non-standard form,   Z διι′ 3(xι + nι )(xι′ + nι′ ) 1 X j;l−l′ p.v. Φii′ ;ss′ (2j x + l′ − l)dx, Tii′ ;ss′ ;ιι′ = 3 − 5 4π 3 kx + nk kx + nk [−1,1] 2 2 3 n∈Z

converge absolutely on all scales j ∈ N.

(ii) For |i + i′ | ≥ 1, the elements (25) of the scaling part of the non-standard form ,   Z 3(xι + nι )(xι′ + nι′ ) διι′ 1 X 0;0 Φii′ (x)dx, Tii′ ;00;ιι′ = p.v. 3 − 4π kx + nk52 [−1,1]3 kx + nk2 3 n∈Z

converge absolutely. (iii) For i = i′ = 0, the elements (26) of the scaling part of the non-standard form vanish,   Z X 1 διι′ 3(xι + nι )(xι′ + nι′ ) 0;0 (29) T00;00;ιι′ = lim p.v. − Φ00 (x)dx = 0. 4π N →∞ kx + nk32 kx + nk52 [−1,1]3 knk∞ ≤N

Proof. The absolute convergence of (i) and (ii) follows directly from Proposition 4. To demonstrate (iii), we show that the sum in (29) is zero for any fixed N . Since the result does not depend on the choice of indices, if ι 6= ι′ , we set ι = 1 and ι′ = 2. Thus, we consider Z N X 3(x1 + n1 )(x2 + n2 ) Φ00 (x1 )Φ00 (x2 )Φ00 (x3 ) dx1 dx2 dx3 . p.v. 3 kx + nk52 [−1,1] n1 ,n2 ,n3 =−N

Since the function Φ00 is even, parity considerations and symmetry of summation with respect to zero, imply that each individual term obtained by expanding (x1 + n1 )(x2 + n2 ) = x1 x2 + n1 x2 + n2 x1 + n1 n2 vanishes. For ι = ι′ , we set ι = ι′ = 1 and consider Z N X −2(x1 + n1 )2 + (x2 + n2 )2 + (x3 + n3 )2 p.v. Φ00 (x1 )Φ00 (x2 )Φ00 (x3 ) dx. kx + nk52 [−1,1]3 n1 ,n2 ,n3 =−N

The three terms in the numerator cancel each other, since the sum for each term is independent of the choice of indices.  Remark 6. The same approach applies to the periodization of the Riesz transform itself.

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

11

3.4. Non-oscillatory Helmholtz Green’s function with periodic boundary conditions. Let us consider the problem
(30) −∆ + µ2 u(x) = f (x) (31) u(x + n) = u(x) for x ∈ [0, 1]3 , µ > 0, n ∈ Z3 , and f ∈ L2 ([0, 1]3 ). Although this problem is easily handled by the standard method of images, we apply our approach in order to show that the limit as µ → 0 does not cover all possible constructions available for the case µ = 0. We consider the solution to (30) and (31) Z u(x) =

where GµH satisfies

[0,1]3

GµH (x − y)f (y)dy,


−∆x + µ2 GµH (x − y) = δ(x − y)

GµH (x − y + n) = GµH (x − y).

We obtain GµH by applying the method of images to the free space Green’s function, Gµfree (x) =

1 e−µkxk2 , 4π kxk2

this time yielding (for x ∈ / Z3 ) the absolutely convergent sum, 1 X e−µkx+nk2 (32) GµH (x) = , x ∈ [0, 1]3 and µ > 0. 4π kx + nk 2 3 n∈Z

By Proposition 4, all elements of the non-standard form for (32) are given by absolutely convergent sums and the usual method of images and that applied to the non-standard form coincide. In the next proposition, we explicitly obtain values of the elements of the scaling part of GµH as functions of µ which, for µ = 0, are given by conditionally convergent sums. Later, in Section 4, we compare these elements with those for the Poisson kernel. ′

′

′

′

Proposition 7. Let Tii0;0 ′ ;00 (µ) with i = (i1 , i2 , i3 ) and i = (i1 , i2 , i3 ) denote an element of the scaling part of the periodized non-standard form of the operator of kernel GµH (x) in (32). It holds that (i) If i = i′ = 0, we have

1 . µ2 ′ ′ (ii) If for any j, 1 ≤ j ≤ 3, ij + ij is odd, then Tii0;0 ′ ;00 (µ) = 0. In particular, if |i + i | = 1 , then 0;0 T00;00 (µ) =

Tii0;0 ′ ;00 (µ) = 0.

(iii) If |i + i′ | = 2,

   0 2 eµ − 12+6µ+µ2 0;0 0;0 ) ( Tii′ ;00 (µ) = Ti′ i;00 (µ) = (12−6µ+µµ4)(e , µ −1)   0,

for i ∈ {(1, 1, 0), (1, 0, 1), (0, 1, 1)} and i′ = (0, 0, 0), for i = i′ ∈ {(1, 0, 0), (0, 1, 0), (0, 0, 1)} ,

for i ∈ {(2, 0, 0), (0, 2, 0), (0, 0, 2)} and i′ = (0, 0, 0).

For µ → 0, we have (33)

0;0 lim Tii0;0 ′ ;00 (µ) = lim Ti′ i;00 (µ) =

µ→0

µ→0

1 , for i = i′ ∈ {(1, 0, 0), (0, 1, 0), (0, 0, 1)} . 60

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

12

See Appendix 8.6 for the proof. The formulas derived in the proof may be used to explicitly compute other elements of the non-standard form. 4. Poisson Green’s function with periodic boundary conditions In this section we consider the problem (34) (35)

−∆u(x) = f (x) u(x + n) = u(x)

for x ∈ [0, 1]3 , n ∈ Z3 and f ∈ L2 ([0, 1]3 ) satisfying the mean-free condition Z f (x)dx = 0. (36) [0,1]3

Due to the slow decay of the free space Green’s function 1 (37) Gf ree (x) = , 4πkxk2 the usual method of images produces a divergent series, X 1 1 X . (38) Gf ree (x + n) = 4π kx + nk2 3 3 n∈Z

n∈Z

However, in our approach to obtain the periodized non-standard form, no “interpretation” of (38) is required, since we apply the method of images not to Gf ree but to its non-standard form. Moreover, using that Gf ree satisfies Proposition 4, all elements of the wavelet part of the non-standard form Gp converge absolutely for any multiwavelet basis with number of vanishing moments m ≥ 3. Thus, to construct Gp , we only need to examine the elements of the scaling part of the non-standard form. By selecting square partial sums as a method of summation (see (26)), these elements are computed as X Z X Φii′ (x) 1 0;0 0;n lim dx, |i + i′ | > 0, (39) Tii′ ;00 = lim Tii′ ;00 = N →∞ 4π N →∞ kx + nk 3 2 [−1,1] knk∞ ≤N

knk∞ ≤N

yielding a particular Green’s function Gp . We remark that other conventions may lead to different variants of the periodized Green’s function consistent with the free space operator Gf ree (see also Remark 9 below). Note that we do not use the sum in (26) to define the element, i = i′ = 0, since in this case it would lead to a divergent sum (see Theorem 8). Instead, we set the value of this element to zero which effectively restricts the domain of the operator Gp to mean-free functions f . Surprisingly, the Green’s function Gp resulting from our summation convention (26), is not a convolution. This is consistent with the fact that the sums in (39) with indices 1 ≤ |i + i′ | ≤ 2 are conditionally convergent and, some of the resulting elements, are not limits, as µ → 0, of the corresponding elements of the convolution operator in (3.4) (see Proposition 7). Theorem 8. Let us consider the non-standard form of the operator Gf ree (37) in a multiwavelet basis with m ≥ 3 vanishing moments. Then (i) The lattice sums in (24) defining wavelet part elements of the periodized non-standard form ′ Tiij;l−l ′ ;ss′

converge absolutely.

=

X

n∈Z3

′ +2j n Tiij;l−l ′ ;ss′

Z Φii′ ;ss′ (2j x + l′ − l) 1 X dx = 4π kx + nk2 3 3 [−1,1] n∈Z

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

13

(ii) For |i + i′ | ≥ 3, the lattice sums defining the scaling part elements of the periodized non-standard form (25), Z Φii′ (x) 1 X 0;0 dx, (40) Tii′ ;00 = 4π 3 kx + nk2 3 [−1,1] n∈Z

converge absolutely.

(iii) For 1 ≤ |i + i′ | ≤ 2, the lattice sums in (39) for the scaling part of the periodized non-standard form X Z 1 Φii′ (x) 0;0 (41) Tii′ ;00 = lim dx, 4π N →∞ [−1,1]3 kx + nk2 knk∞ ≤N

converge conditionally.

0;0 (iv) For |i + i′ | = 0, with the summation convention (iii), the lattice sum for the element T00;00 0;0 diverges. By setting it to zero, T00;00 = 0, we effectively restrict the domain of the periodized R operator to the class of functions with zero mean [0,1]3 f (x)dx = 0.

See Appendix 8.4 for the proof.

Remark 9. The fact that only a few elements of the non-standard form are given by conditionally convergent sums permits a characterization of all possible versions of the periodic Poisson Green’s function. Our approach offers a unified way of constructing such Green’s functions and, perhaps, explains difficulties encountered in their usual interpretation. Some of these different periodic Green’s functions may be found in the literature [19, 15, 30]. The fact that in computing the periodic Poisson Green’s function one encounters conditionally convergent sums is well known. Assigning different values to such sums explains the differences in e.g., [20] and [32] approaches to lattice summation. A particular choice is made in the context of the Fast Multipole Method [23, Section 4]. For a discussion of this issue see [22, Section 3]. In the next proposition we obtain the values of several elements Tii0;0 ′ ;00 of the non-standard form. In particular, we obtain all values of the elements given by conditionally convergent series. Recall that those elements correspond to indexes satisfying |i + i′ | ≤ 2. Proposition 10. Let Tii0;0 ′ ;00 with i = (i1 , i2 , i3 ) and i = (i1 , i2 , i3 ) denote an element of the scaling part of the periodized non-standard form of the operator of kernel Gfree in (37). It holds that ′

′

′

′

′ (i) If for any j, 1 ≤ j ≤ 3, ij + ij is odd, then Tii0;0 ′ ;00 = 0. In particular, if |i + i | is odd, then ′

Tii0;0 ′ ;00 = 0.

(ii) If |i + i′ | = 2, Tii0;0 ′ ;00

=

Ti0;0 ′ i;00

=

  0

2 ,  45

−

1√ , 36 5

for i ∈ {(1, 1, 0), (1, 0, 1), (0, 1, 1)} and i′ = (0, 0, 0), for i = i′ ∈ {(1, 0, 0), (0, 1, 0), (0, 0, 1)} , for i ∈ {(2, 0, 0), (0, 2, 0), (0, 0, 2)} and i′ = (0, 0, 0).

See Appendix 8.7 for the proof. Remark 11. If in the Proposition above |i + i′ | ≥ 3 and one of the coordinates of the multi-indices ′ ′ is zero, e.g. i = (i1 , i2 , i3 ), i′ = (0, i2 , i3 ) but i1 > 0 then See Appendix 8.8 for the proof.

Tii0;0 ′ ;00 = 0.

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

14

0;0 0;0 We observe that the non-zero values of Ti0;00 and T0i;00 for |i + i′ | = 2 are due to the slow decay of the kernel. Indeed, comparison of Propositions 7 and 10 shows that, for some indices ii′ , the limits 0;0 of Tii0;0 ′ ;00 (µ) as µ → 0 do not match the values of Tii′ ;00 in Proposition 10. Thus, this mismatch may be attributed to the slow decay of the Poisson kernel. On the other hand, there is no mismatch for all terms defined by absolutely convergent sums since in that case the order of summation and integration may be exchanged.

Moreover, if we were to modify Gf ree outside of an arbitrarily large ball of radius R as to increase the rate of decay from 1/R to 1/R1+δ , δ > 0, then no mismatch will occur in e.g. the elements with indices i ∈ {(2, 0, 0), (0, 2, 0), (0, 0, 2)} and i′ = (0, 0, 0). In fact, a much stronger result is true. Proposition 12. Consider a kernel G(x1 , x2 , x3 ), locally integrable, even on each coordinate and such that, for some positive δ and M , its partial derivatives satisfy (42)

kGxj (x)k2 ≤

C , for kxk2 ≥ M, kxk2+δ 2

where C is a constant. If ϕj , 1 ≤ j ≤ 3, denote three bounded functions on [0, 1] and one of them is odd about 1/2 then X Z ϕ1 (x1 )ϕ2 (x2 )ϕ3 (x3 )G(x1 + n1 , x2 + n2 , x3 + n3 )dx = 0. lim N →∞

knk∞ ≤N

[0,1]3

In particular, the scaling elements of the periodized non-standard form of G, Tii0;0 ′ ;00 , vanish for ′ i ∈ {(2, 0, 0), (0, 2, 0), (0, 0, 2)} and i = (0, 0, 0). See Appendix 8.9 for the proof. Remark 13. We may also consider the limit using the oscillatory Helmholtz kernel eiκr /r. Sending κ → 0 (as in [10]) yields a particular periodic Green’s function for the Poisson’s kernel also obtainable by Ewald’s method [20]. Thus, in practical applications, the selection of the Green’s function of the periodized Poisson kernel may depend on physical considerations that either emphasize the long range behavior of this kernel or use its properties only in a finite region. The effect of such choice on the solutions and their behavior on the boundary of a periodic cell is further discussed in the next section. 4.1. On mean-free and weak solutions of the periodic Poisson equation. Let us show that our construction yields a solution of the periodic problem (46)-(47) that is not mean free which, in turn, implies that the periodized operator is not a convolution. Note that if u is a solution of (46)-(47), then Z u(y)dy, (43) u(x) − [0,1]3

R is a mean-free solution of the same problem. However, in our construction [0,1]3 u(y)dy may not be zero as we demonstrate below. Since Z X 0;0 Z f (x1 , x2 , x3 )ψi′ ;0 (x1 )ψi′ ;0 (x2 )ψi′ ;0 (x3 )dx1 dx2 dx3 , u(x)dx = T0i′ ;00 (44) [0,1]3

i′

[0,1]3

1

2

3

where ψi′ ;0 are the one-dimensional scaling functions defined in (1), from Theorem 8 part (iv) and j

Proposition 10, we conclude that the only non-zero terms of the sum in (44) correspond to the three

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

15

multi-indices i′ = (0, 0, 2), i′ = (2, 0, 0) and i′ = (0, 0, 2). Hence, we obtain Z Z 1 u(x)dx = − f (x1 , x2 , x3 ) [P2 (2x1 − 1) + P2 (2x2 − 1) + P2 (2x3 − 1)] dx1 dx2 dx3 . 36 [0,1]3 [0,1]3 Expanding P2 (2t − 1) = 1 − 6t + 6t2 and using that f is mean-free, the last equation is equivalent to Z Z
1 f (x1 , x2 , x3 ) x1 + x2 + x3 − x21 − x22 − x23 dx1 dx2 dx3 . u(x)dx = (45) 6 [0,1]3 [0,1]3 This last condition is also derived in the literature (but with more restrictive assumptions on the function f ). See, e.g., [5, Eq. 29], [26, Eq. 38] or [29, Eq. 8]. Further analysis of (45) leads us to consider the weak formulation of the problem (46)-(47), Z Z (46) f (x)ϕ(x) dx ∇u(x) · ∇ϕ(x) dx = [0,1]3

[0,1]3

(47)

u(x + n) = u(x)

where the test functions ϕ ∈ C ∞ ([0, 1]3 ) are supported on [0, 1]3 and also satisfy (47). Defining p0 (x) = x1 + x2 + x3 − x21 − x22 − x23 ,

(48)

and integrating by parts on the left hand side of (46), yields Z Z Z u(x)dx − f (x)p0 (x)dx = 6

u dS,

∂([0,1]3 )

[0,1]3

[0,1]3



where ∂ [0, 1]3 denotes the boundary of the unit box and dS is the measure on ∂ [0, 1]3 . Combining the last equation with (45), we obtain that our construction produces a solution with the additional property Z u dS = 0, (49) ∂([0,1]3 )

i.e., the integral of the solution over the boundary vanish. 4.2. An analytic expression for the periodized Green’s function. The non-standard form approach for the construction of the periodic Poisson kernel provides the coefficients in the multiwavelet basis of a solution u for the problem (34)-(35) under the assumption (36). The solution u so obtained is not mean-free and satisfies the boundary condition (49). Let us now describe analytically the Green’s function that yields this solution for the problem. Let us consider X e2πin·x ∈ L2 ([0, 1]3 ). G0 (x) = 4π 2 knk2 n6=0

Formally, G0 solves the problem (50) (51)

−∆G0 (x) = −1+ P

G0 (x + n)

=

X

n∈Z3

δ(x − n)

G0 (x),

where n∈Z3 δ(x − n) is the periodic delta function (for the box [0, 1]3 ). The mean-free condition (36) on the function f yields a solution u for the problem (34)-(35) as Z G0 (x − y)f (y) dy. u(x) = [0,1]3

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

Since the periodicity of G0 yields Z

[0,1]3

G0 (x − y) dx =

Z

[0,1]3

16

G0 (x) dx

we also have that the solution u is mean-free. We now modify G0 as to obtain a Green’s function G yielding the boundary condition (49). Note that for y = (y1 , y2 , y3 ) ∈ [0, 1]3 , Z 3 X −2πinj yj 3 X X e 1 (yj2 − yj + 1/6) = −p0 (y) + , G0 (x − y) dx = 2 (52) = 2 2 2 4π nj ∂([0,1]3 ) j=1 j=1 nj 6=0

where p0 is the polynomial in (48). Let’s define for x, y ∈ [0, 1]3

G(x, y) = G0 (x − y) + G1 (x, y),

where

  3 X 1 xj y j  kxk2 + kyk2 − 2 G1 (x, y) = − 6 j=1

3 which we extend periodically as G1 (x+n, y) P= G1 (x, y) and G1 (x, y+n) = G1 (x, y) for x, y ∈ [0, 1] 3 and any n ∈ Z . Although −∆G(x, y) = n∈Z3 δ(x − n), we observe that G is not a convolution, since k (x − y) (mod 1)k2 6= kx (mod 1) − y (mod 1)k2 , where mod 1 indicates periodization on the unit box. For the Green’s function G, the corresponding solution u satisfies Z G(x, y)f (y) dy, u(x) = [0,1]3

yielding Z Z u(x) dx = [0,1]3

[0,1]3

Z

!

1 G1 (x, y) dx f (y) dy = 6 [0,1]3

Z

1 (p0 (y) − 1) f (y) dy = 6 [0,1]3

Z

[0,1]3

p0 (y)f (y) dy,

which coincides with (45). On the other hand, combining Z 7 G1 (x, y) dx = − + p0 (y) 6 ∂([0,1]3 ) with (52) we obtain Z

∂([0,1]3 )

u(x) dx =

Z

[0,1]3

Z

∂([0,1]3 )

(G0 (x − y) + G1 (x, y)) dx f (y) dy = 0.

We refer to e.g. [4] for a different construction of G0 . 5. Separated representations We use approximation via Gaussians as a tool for constructing separated representations of operator kernels to obtain fast algorithms for their application. Such approximation separates along each coordinate direction, thus simplifying the computation of the non-standard form and yielding a fast algorithm to apply the operator. Approximation via Gaussians (see e.g. [11, 12, 13, 14]) has been successfully used in [24, 9, 6] to construct fast and accurate algorithms for applying free space convolution kernels for any user supplied finite accuracy. Our goal in this section is to extend this approach to periodized kernels constructed in Sections 3 and 4. As an example, we consider convolutions with kernels of the form (53)

−µkxk2 K(x) = p1 (x1 )p2 (x2 )p3 (x3 )kxk−β , 2 e

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

17

where β and µ are non-negative parameters, both not simultaneously zero, and pγ is a polynomial, −µkxk2 , may be efficiently approxiγ = 1, 2, 3. We note that both, kxk−β and e−µkxk2 , or kxk−β 2 2 e mated by short sums of Gaussians for any user selected accuracy ǫ and distance from the origin δ (see Theorem 6 and Proposition 8 of [14]). In fact, the number of terms is shown to be proportional
2 to log δ−1 and log ǫ−1 (although in practice we observe essentially log ǫ−1 dependence). Substi−µkxk2 yields a separated representation tuting in (53), the approximation by Gaussians of kxk−β 2 e of the free space kernel K. In this section we show that the periodized operator has a separated representation as well. Once equipped with the separated representation of the non-standard form, we may use the algorithms described in [24, 9, 6] (with minor modifications) to apply the periodized non-standard form. Such algorithms have the same complexity as those for the free space operators. We may write (54) where

−β r − Gβ (r) ≤ ǫr −β , for all δ ≤ r ≤ R,

(55)

Gβ (r) =

Nβ X

an e−αn r

2

n=1

with positive an and αn , and

Gβ (r) < (ǫ + 1) r −β , for all r > 0.

(56)

The bound (56) may be obtained following the derivation of [14, Lemma 4]. We may also write −µr ǫ e , for δ ≤ r ≤ R, (57) − Gµ (r) ≤ ǫ+1 where (58)

Gµ (r) =

Nµ X

dn e−ηn r

2

n=1

and dn and ηn are positive. Hence, combining (54)-(58), we obtain a sum of Gaussians approximation for r −β e−µr in the range r ∈ [δ, R],  h i  −β −µr (59) − Gβ (r)Gµ (r) ≤ r −β − Gβ (r) e−µr + e−µr − Gµ (r) Gβ (r) ≤ 2ǫr −β . r e

The number of terms in the sub-optimal approximation Gβ (r)Gµ (r) may be reduced further by using the reduction algorithms in [13, 14]. As a consequence, we obtain an approximation of the kernel (53) as ˜ K(x) = p1 (x1 )p2 (x2 )p3 (x3 )

M X

−τm kxk22

wm e

=

M X

2

2

2

wm p1 (x1 )e−τm x1 p2 (x2 )e−τm x2 p3 (x3 )e−τm x3 ,

m=1

m=1

where the number of terms, M , depends logarithmically on ǫ and δ, and the parameters τm and wm ˜ the non-standard form inherits the separation along are positive. Due to the functional form of K, each coordinate direction and we obtain M X j;l −l′ j;l −l′ j;l −l′ j;l−l′ ˜ wm t˜i1 i1′ ;s11s′ ;m;1 t˜i2 i2′ ;s22s′ ;m;2 t˜i3 i3′ ;s33s′ ;m;3 , (60) Tii′ ;ss′ = m=1

where (61)

′ t˜j;l−l ii′ ;ss′ ;m;γ

=

Z Z R

1

1

2

2

2

R

3

3

′

j,l (x)ψij,l pγ (x − y)e−τm (x−y) ψi,s ′ ,s′ (y)dxdy.

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

18

′ Thus, in order to compute T˜iij;l−l ′ ;ss′ , it is sufficient to evaluate one dimensional integrals with the cross-correlations of the scaling functions (see (20)), Z ′ 2 pγ (x)e−τm x Φii′ (2j+1 x + l − l′ )dx = t˜j+1;l−l ii′ ;00;m;γ

R

and then apply the quadrature mirror filters for the multiwavelets (see [3, eq 3.25a 3.25b 3.25c ′ 3.25d]) to construct all the coefficients t˜j;l−l ii′ ;ss′ ;m;γ for s = 11, 10, 01. We note that to apply the operator we may also use the modified non-standard form [6] which only requires the projection of the operator onto cross-correlation functions of the scaling functions. Applying the method of images to (60), we obtain the coefficients of the non-standard form of the operator with periodic boundary conditions, ′ T˜iij;l−l ′ ;ss′ =

(62) where in each direction

M X

′

′

1

m=1

′

1

˜tj;l−l ii′ ;ss′ ;m;γ =

(63)

′

j;l −l j;l −l j;l −l wm˜ti1 i1′ ;s11s′ ;m;1˜ti2 i2′ ;s22s′ ;m;2˜ti3 i3′ ;s33s′ ;m;3 , 2

2

X

′

3

3

j

j;l−l +2 n t˜ii ′ ;ss′ ;m;γ ,

n∈Z j;l−l′ +2j n

with t˜ii′ ;ss′ ;m;γ defined in (61). Clearly (62) is in separated form with the same separation rank as its free space counterpart (60) and, moreover, (63) provides a simple recipe for computing its components. ′ +2j n Remark 14. By first computing the blocks T˜iij;l−l of the non-standard form of the free space ′ ;ss′ ′ ˜ we have a simple way to evaluate via (63) the corresponding blocks T˜ j;l−l approximation K, ii′ ;ss′ for ′ +2j n the approximation of the periodized operator. Since the norm of the blocks t˜j;l−l ii′ ;ss′ ;m;γ in (63) decays rapidly with n, only a few terms participate in the sum for a given accuracy. In fact, on finer scales (large j) only the ′term with n = 0 needs to be kept. We may estimate the error j;l−l′ ˜ j;l−l′ j;l−l Tii′ ;ss′ − Tii′ ;ss′ , where Tii′ ;ss′ are the blocks of the non-standard form of the original operator ′ +2j n ′ +2j n K, by using Proposition 4 together with the estimates for Tiij;l−l − T˜iij;l−l given in [9, ′ ;ss′ ′ ;ss′ Theorem 4]. However, an exception to using (63) for computing operator blocks has to be made for conditionally convergent elements on the coarsest scale whose definition reqiures special attention (see Proposition 10).

Remark 15. Our approach applies to any Bravais lattice. We note that for a non-rectangular lattice the non-standard form does not separate along each coordinate and further approximations are required. 6. Dirichlet, Neumann and mixed boundary conditions Using the results for the periodic case, we now have the necessary tools to efficiently apply operators with Dirichlet, Neumann or mixed boundary conditions on simple domains. We note that although the resulting integral operators are no longer convolutions, they have a simple algebraic structure and, as a result, the algorithm for applying these operators is similar to those described in the previous section.

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

19

As an example, let us consider the problem
−∆ + µ2 u(x) = f (x) for x ∈ D (64) (65) u(x) = 0 for x ∈ ∂D,

where µ ≥ 0 and D = [−1/2, 1/2]3 . A solution to (64) which satisfies (65) is given by Z Gµ (x, y)f (y)dy, u(x) = D

where

Gµ

satisfies


−∆x + µ2 Gµ (x, y) = δ(x − y) Gµ (x, y) = 0 for x ∈ ∂D

(66) (67)

and ∆x denotes the Laplacian with respect to x. Let us first consider the case where µ > 0. Even though the integral operator Gµ is not a convolution, it may be written as     x1 − y 1 x2 − y 2 x3 − y 3 x1 − y 1 x2 − y 2 x3 + y 3 + 1 µ µ µ G (x, y) = GH , , , , − GH 2 2 2 2 2 2     x1 − y 1 x2 + y 2 + 1 x3 − y 3 x1 − y 1 x2 + y 2 + 1 x3 + y 3 + 1 µ µ , , , , − GH + GH 2 2 2 2 2 2     x1 + y 1 + 1 x2 − y 2 x3 + y 3 + 1 x1 + y 1 + 1 x2 − y 2 x3 − y 3 µ µ + GH , , , , − GH 2 2 2 2 2 2     x1 + y 1 + 1 x2 + y 2 + 1 x3 + y 3 + 1 x1 + y 1 + 1 x2 + y 2 + 1 x3 − y 3 µ µ , , , , (68) + GH − GH , 2 2 2 2 2 2

where the periodic Green’s function GµH is constructed as in Section 3.4 to satisfy
1 (69) −∆x + 4µ2 GµH (x − y) = δ(x − y). 2 The changes in the equation relative to (66) are due to the way variables appear in (68) and to the dimension of the space, d = 3. Since GµH has period one and is even in each variable, for x ∈ ∂D the terms in (68) cancel each other so
that Gµ satisfies the Dirichlet boundary condition (67). For x 6= y inside D, we have −∆x + µ2 Gµ (x, y) = 0 since each of the eight terms in (68) vanishes. The only singularity is at x = y, in which case the first term in (68) yields (66).

The non-standard form of Gµ is then constructed by using Propositions 4 and 7 for each term in (68). However, in contrast with Proposition 7, part (i), in the next proposition we show that the 0;0 element T00;00 of the non-standard form of Gµ is finite as µ → 0. This allows us to obtain Gµ for µ = 0. We note that, unlike the periodic Green’s function, the Green’s function for the Dirichlet problem is unique. 0;0 Proposition 16. The element T00;00 of the non-standard form of Gµ is given by
−2
−2
−2 n2 + 21 n3 + 12 n1 + 12 1 X 0;0 T00;00 = 4
2 1 2 π ) + (n2 + 1 )2 + (n3 + 1 )2 + µ 3 (n1 + n∈Z

2

2

2

2π

which converges to a positive constant as µ → 0.

0;0 Proof. Using (22), we write the element T00;00 of Gµ as

Z

D3

Z

D3

µ

G (x, y) dx dy =

Z

D3

Z

8 X (−1)l+1 GµH (al (x, y)) dx dy,

D 3 l=1

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

20

where al (x, y) denotes the argument of the lth term in (68). As in the proof in Section 8.6, we compute the Fourier coefficients of GµH (x) and obtain X 1 X e−2µkx+nk2 = cn e2πin·x , 2π kx + nk 2 3 3

GµH (x) =

n∈Z

with cn =

1 1 . 2π 2 n2 +n2 +n2 +( µ )2 1 2 3 π 0;0 T00;00

(70)

n∈Z

Integrating to obtain =

X

Z

cn

D3

n∈Z3

P8

Z

D3

8 X

0;0 T00;00 ,

we write

(−1)l+1 e2πin·al (x,y) dx dy.

l=1

Note that the integrand l=1 (−1)l+1 e2πin·al (x,y) may be expressed in separated form as                x +y +1 x +y +1 x +y +1 x −y x −y x −y 2πin1 1 2 1 2πin2 2 2 2 2πin3 3 2 3 2πin2 2 2 2 2πin3 3 2 3 2πin1 1 2 1 −e −e −e e e , e so that

R

R

D3 D3

P8

l+1 e2πin·al (x,y) dx dy l=1 (−1) 3 Z Z  Y

j=1 D

=

j −yj 2



j −yj 2



2πinj

x

2πinj

x

e

equals to

D

3 Z Z  Y

j=1 D

e

D

=

3 Z Z Y

j=1 D

2πinj

e

2πinj

x

2πinj

x

−e −e

x

j −yj 2



D

j +yj +1 2

j −yj +1 2

 

dxj dyj

dxj dyj

(1 − (−1)nj ) dxj dyj

where, for each j, we changed variables yj 7→ −yj . Therefore, we may rewrite the series (70) using n = (n1 , n2 , n3 ) with only odd indices nj . Thus, we compute Z 1 2  x −y  Z Z 2 1 1 2πi(2nj +1) j 2 j dxj dyj = e e2πi(nj + 2 )xj dxj =
2 . −1 D D π 2 nj + 1 2

2

Combining integrals in each coordinate, we obtain the result.



6.0.1. Separated representation of Gµ . The number of terms in the construction of the Green’s function satisfying boundary conditions in (68) grows exponentially with the dimension, i.e. if d = 2 we have four terms, of d = 4 we have 16 terms, etc. On the other hand, the number of terms ˜ µ , of Gµ is independent of the dimension. Indeed, using (59), we in the separated approximation, G µ approximate G by (71)

˜ µ (x, y) = G

M X

m=1

wm

X

Sm,n1 (x1 , y1 )

n1 ∈Z

X

Sm,n2 (x2 , y2 )

X

n2 ∈Z

n3 ∈Z

2

2

Sm,n3 (x3 , y3 )

where Sm,n (x, y) = e−τm (x−y+n) − e−τm (x+y+n+1) .

In other words, the approximations of individual terms in (68) combine to yield (71). Thus, to ˜ µ , we only need to compute the integrals, construct the non-standard form of G Z 1Z 1 ′ 2 j,l − ˜j;l−l′+2j n tii′ ;ss′ ;m = e−τm (x−y+n) ψij,l ′ ,s′ (y)ψi,s (x)dxdy 0

0

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

and ′ j + ˜j;l+l +2 (n+1) tii′ ;ss′ ;m

=

Z

0

1Z 1 0

21

′

j,l e−τm (x+y+1+n) ψij,l ′ ,s′ (y)ψi,s (x)dxdy 2

′

j

+2 n and for j ∈ N, n ∈ Z, l, l′ ∈ {0, . . . , 2j − 1}, i, i′ ∈ {0, . . . , m − 1}. The integrals − t˜j;l−l ii′ ;ss′ ;m ′ +2j n j;l−l + t˜ are simplified further and reduce to one dimensional integrals using cross and autoii′ ;ss′ ;m correlations of wavelet and scaling functions (see Section 5). As a result, the non-standard form is given by ′ T˜iij;ll ′ ;ss′ =

where ′

M X

1

m=1

t˜j;ll ii′ ;ss′ ;m =

′

′

′

j;l l j;l l j;l l wm t˜i1 i1′ ;s1 1 s′ ;m t˜i2 i2′ ;s2 2 s′ ;m t˜i3 i3′ ;s3 3s′ ;m ,

X

1

− ˜j;l−l′ +2j n tii′ ;ss′ ;m

n∈Z

2

2

3

3

′

jn

j;l+l +2 −+ t˜ii ′ ;ss′ ;m



.

Remark 17. Although we discussed the Poisson Green’s function with Dirichlet boundary conditions, this approach extends to any operator which is effectively sparse in the non-standard form and whose kernel may be approximated by a separated representation. Also we may use the same approach for operators with Neumann or mixed boundary conditions.

Remark 18. We note that the Fast Multipole method provides an alternative approach to the treatment of boundary conditions, see [23, Section 4].

7. Conclusions and remarks We have described an approach to construct and apply a class of operators with periodic boundary conditions. The non-standard form of the corresponding free space operator provides the foundation for our approach and allows us to analyze these operators on a hierarchy of scales. This analysis is operator independent and reveals that the wavelet part of the non-standard form is always well defined. Depending on the properties of the kernel for large arguments, we have shown that the scaling part of the non-standard form may have elements which require special attention. With the use of separated representations via Gaussians, we obtain fast algorithms for application of these operators that are minor modification of their free space versions. For simple domains, we construct Green’s functions satisfying Dirichlet, Neumann or mixed boundary conditions also yielding separated representations of operators and fast algorithms for their application. We would like to note that it may be possible to use our approach as a tool for constructing Green’s functions for finite size lattices. While interior cells may be well approximated by a periodic construction, cells near the boundary usually require a different approximation. In this scenario Fourier methods are not available since the Poisson summation formula no longer applies, while the direct summation is not computationally effective. In contrast, the multiresolution approach (for a given accuracy) only requires modifications in the vicinity of the boundary on all scales except for the coarsest. Indeed, due to the rapid decay of the lattice sums on wavelet subspaces, only a few elements are affected by their neighbors.

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

22

8. Appendix 8.1. Properties of cross-correlation functions. We use the following properties of Φii′ (x) which follow from [9, Proposition 3 ], ′

(72)

Φii′ (x) = (−1)i+i Φi′ i (x),

(73)

Φii′ (−x) = (−1)i+i Φii′ (x),

′

(74) (75)

Z

(76)

Φi0 (1 − x) = (−1)i+1 Φi0 (x) for i > 0, x ∈ [0, 1], Φ00 (1 − x) = 1 − Φ00 (x) for x ∈ [0, 1],

1

Φ00 (x)dx = 1.

−1

8.1.1. Moments on the interval [0, 1]. From the definition (2), we have Z 1 √ (77) Φi0 (x) = 2i + 1 Pi (2t − 1)dt, x ≥ 0. x

In particular, (78)

which, together with (74), implies (79)

Φi0 (1) = 0, i ≥ 0 Φi0 (0) = 0, i > 0.

Integrating by parts Z

1

Φi0 (x)

0



xk+1 k+1

′

√

2i + 1 dx = k+1

Z

1 0

Pi (2x − 1)xk+1 dx

and using the orthogonality of the Legendre polynomials, we obtain Z 1 Φi0 (x)xk dx = 0 for i > k + 1 and k ≥ 0. (80) 0

8.2. Proof of Proposition 2. It is enough to prove the proposition for |l − l′ | ≥ 3 since the ′ general resultfollows by modifying C j to ′include the case |l − l | ≤ 2. We first prove the estimate for j;l−l′ R R j;l j;l Tii′ ,10 = R R K(x − y)ψi;1 (x)dx ψi′ ;0 (y) dy . Denoting Il = [2−j l, 2−j (l + 1)], we note that the j;l j;l multiwavelet ψi;1 is supported on Il and the scaling function ψi;0 on Il′ . For each y ∈ Il′ , consider −j−1 the Taylor expansion of the function K(· − y) centered at x0 = 2 (2l + 1), the mid-point of Il ,

K(x − y) = K(x0 − y) + · · · +

K (ν−1) (x0 − y) K (ν) (ξ − y) (x − x0 )ν−1 + (x − x0 )ν (ν − 1)! ν!

where ν = min {m, m} and ξ is between x and x0 and, hence, ξ ∈ Il . Since the multiwavelets have vanishing moments (4), for 0 ≤ n ≤ ν − 1  Z Z ′ 1 j;l (x − x0 )n ψi;1 (x)dx K (n) (x0 − y)ψij;l ′ ;0 (y)dy = 0. n! R R For the remainder term in the Taylor expansion, using (12), we obtain Z Z Z 1 cν j;l′ (ν) ν j;l ν j;l (y)dydx ≤ K (ξ − y)ψ − x ) ψ (x) (x − x ) ψ (x) (x ′ 0 0 i;1 i;1 i ;0 ν! ν! Il R R Z 1 j;l′ (y) (81) dydx. ν+β ψi′ ;0 Il′ |ξ − y|

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS j;l Using Hölder’s inequality and kψi;0 kL2 (R) = 1 we estimate

Z

I l′

Z

1 j;l′ ψ (y) dy ≤ ′ |ξ − y|ν+β i ;0

I l′

1 dy |ξ − y|2(ν+β) Z 1

≤ 2j(ν+β−1/2)

0

!1/2

′ ′

j ;l

ψi′ ;0

23

L2 (R)

1 du |l − l′ − (u − η)|2(ν+β)

1/2

,

where we changed variables u = 2j y − l′ and used that ξ ∈ Il to write ξ = 2−j (η + l) for η ∈ [0, 1]. Since |(u − η)| ≤ 1 and |l − l′ | ≥ 3 we obtain |l − l′ − (u − η)| ≥ |l − l′ | − |(u − η)| ≥ |l − l′ | − 1 ≥ (1 + |l − l′ |)/2.

(82)

The other term in (87) is estimated as 1/2 Z Z



j;l 2ν ν j;l (x − x0 ) dx

ψi;1 (x − x0 ) ψi;1 (x) dx ≤

L2 (R)

Il

Il

=

s

2−j−2ν−2jν , 2ν + 1

j;l since kψi;1 kL2 (R) = 1.

Combining these estimates we obtain the result with cν Cj = √ 2j(β−1)+β . ν! 2ν + 1 ′ j;l−l′ The proof for Tiij;l−l ′ ,01 and Tii′ ,11 follows in a similar fashion because on each of these terms at least one multiwavelet is present. 0;l−l′ R 1 ′ ′ It remains to prove the estimate for Tii′ ,00 = −1 K(x + l − l )Φii (x)dx . First assume i + i′ ≥ 1.

This time we use the Taylor expansion of K(· − (l′ − l)) centered at x0 = 0, so that ′

K(x + l − l ) =

ν−1 X K (n) (l − l′ )

n=0

n!

xn +

K (ν) (ξ − (l′ − l) ν x , ν!

where ν = min {i + i′ , m} ≥ 1 and ξ is between 0 and x ∈ [−1, 1] , and thus |ξ| ≤ 1. Due to the vanishing moments of Φii′ (3) the first ν terms in the Taylor expansion vanish. Using (12) Z 1 ν Z 1 |xν Φii′ (x)| x (ν) ′ cν ′ (x)dx K (ξ + l − l )Φ dx ≤ ii ν! ′ ν+β −1 |ξ + l − l | −1 ν! ≤

where aν = maxi,i′ (83)

nR

1 ν ′ −1 |x Φii (x)| dx

o

cν aν 2ν+β

ν! (1 + |l − l′ |)ν+β

and we estimated

|ξ + l − l′ | ≥ |l − l′ | − |ξ| ≥ |l − l′ | − 1 ≥ (1 + |l − l′ |)/2

using |ξ| ≤ 1 and |l − l′ | ≥ 3. The result follows with C0 =

cν aν 2ν+β . ν!

For the case i = i′ = 0, we first use (12) to bound the kernel and then apply an estimate equivalent to (83).

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

24

8.3. Proof of Proposition 4. √ Proof. It is enough to prove the result for kl − l′ k2 ≥ 2 d + 1 since the general result follows by √ modifying Cj to include the case kl − l′ k2 < 2 d + 1. We first prove the estimate for ′

Tiij;l−l ′ ;ss′ =

(84)

Z

Rd

Z

′

Rd

j;l K(x − y)Ψj;l i′ ;s′ (y)Ψi;s (x)dydx,

with ss′ 6= 00. Let’s assume that s 6= 0 and denote Il = [2−j l1 , 2−j (l1 + 1)]× · · · × [2−j ld , 2−j (ld + 1)]. j;l′ Thus, Ψj;l i;s is a multiwavelet supported on Il while the function Ψi′ ;s′ is supported on Il′ . For each y ∈ Il′ let us consider the Taylor expansion
of the function K(· − y) centered at x0 = 2−j−1 (2l1 + 1), . . . , 2−j−1 (2ld + 1) , (85)

K(x − y) =

X 1 1 α D K(x0 − y)(x − x0 )α + D α K(ξ − y)(x − x0 )α , α! α!

X

|α|=ν

|α|≤ν−1

where ν = min {m, m} and ξ = (1 − θ)x0 + θx, with θ ∈ [0, 1] and, hence, ξ ∈ Il . We write ξ = 2−j (η + l) , η ∈ [0, 1]d .

(86)

Substituting (85) into (84) and using that the multiwavelets have vanishing moments (4), we observe that all terms with |α| ≤ ν − 1 do vanish. For the remainder term in the Taylor expansion, using (23), we obtain Z X cα Z 1 j;l′ j;l−l′ α j;l (y) (87) dydx. Ψ (x − x0 ) Ψi;s (x) Tii′ ;ss′ ≤ ′ ′ i ;s |α|+β α! Il I ′ kξ − yk l 2 |α|=ν ′

Hölder’s inequality and kΨj;l i′ ;s′ kL2 (Rd ) = 1 yield Z

1 Il′

|α|+β

kξ − yk2

j;l′ (y) Ψ i′ ;s′ dy ≤

Z

1

Il′

2(|α|+β)

kξ − yk2

dy

!1/2

.

By changing variables u = 2j y − l′ in the last integral and using that y ∈ Il′ and (86) we obtain Z Z 1 1 2j(|α|+β)−jd dy = 2 du. 2(|α|+β) 2(|α|+β) Il′ kξ − yk2 [0,1]d kη − u + l − l′ k2 √ √ Since kη − uk2 ≤ d and kl − l′ k2 ≥ 2 d + 1, we estimate kη − u + l − l′ k2 ≥ kl − l′ k2 −

√

d≥

1 + kl − l′ k2 , 2

and therefore Z

Il′

1 kξ −

2(|α|+β) yk2

dy

!1/2

≤ 2−dj/2 2(j+1)(|α|+β) 1 + kl − l′ k2

Substituting the last inequality into (87), we now bound the integral 1/2 Z Z



j;l 2α α j;l (x − x0 ) dx

Ψi;s 2 (x − x0 ) Ψi;s (x) dx ≤ Il

Il

L (Rd )


−|α|−β

.

2−dj/2 2−(j+1)|α| = Qd √ , r=1 2αr + 1

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS



j;l where we used Ψi;s Z

Il

L2 (Rd )

25

= 1 and 2−j (lr +1) 

 1 2αr t − 2 (lr + ) dt 2 −j r=1 2 lr  Z 1 d d Y Y 2−2αr (j+1) 1 2αr −dj −2αr j−j du = 2 . 2 = u− 2 2αr + 1 0 r=1 r=1

d Z Y (x − x0 )2α dx =

−j

Combining these estimates we obtain the result with X c √ α Cj = 2β 2−j(d−β) . α! 2α + 1 |α|=ν

It remains to prove the estimate for ′

Tii0;l−l ′ ;00 =

(88)

Z

[−1,1]d

K(x + l − l′ )Φii′ (x)dx.

First assume |i + i′ | ≥ 1. This time we use the Taylor expansion of K(· + l − l′ ) centered at the origin, so that X 1 X 1 D α K(l − l′ )xα + D α K(l − l′ + θx)xα , K(x + l − l′ ) = α! α! |α|=ν

|α|≤ν−1

where ν = min {|i + i′ |, m} ≥ 1 and θ ∈ [0, 1]. Substituting into (88) and using that Φii′ have vanishing moments (3), we observe that all terms with |α| ≤ ν − 1 do vanish. For the remainder term in the Taylor expansion, using (23), X cα aα X cα Z 2|α|+β |xα Φii′ (x)| 0;l−l′ dx ≤ , Tii′ ;00 ≤ α! [−1,1]d kl − l′ + θxk|α|+β α! (1 + kl − l′ k2 )|α|+β 2 |α|=ν |α|=ν nR o where aν = maxi,i′ ,|α|=ν [−1,1]d |xα Φii′ (x)| dx and we estimated (89)

kl − l′ + θxk2 ≥ kl − l′ k2 − kθxk2 ≥ kl − l′ k2 −

The result follows with C0 = 2ν+β

√

d≥

1 + kl − l′ k2 . 2

X cα aα . α!

|α|=ν

For the case i = i′ = 0, we first bound the kernel in (88) and then apply an estimate equivalent to (89).  8.4. Proof of Theorem 8. Proof. The absolute convergence in (i) and (ii) follows directly from Proposition 4 with β = 1. For (iii-iv) we use the Taylor expansion   1 kxk22 3(x · n)2 1 1 x·n − + + O (90) = − kx + nk2 knk2 knk32 2knk32 2knk52 knk42

for x ∈ [−1, 1]3 and n 6= 0 . Note that the case n = 0 corresponds to the the elements of the non-standard form of the free-space Green’s function, which we assume well defined and that the

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

26

remainder term in (90) leads to an absolutely convergent sum. We proceed by substituting the first four terms of the expansion (90) into the integrand of (41), and consider   X Z kxk22 3(x · n)2 1 x·n − + Φii′ (x)dx. (91) lim − 3 3 5 N →∞ knk knk 2knk 2knk 3 2 [−1,1] 2 2 2 knk ≤N ∞ n6=0

By symmetry considerations and using that Φii′ (x) = Φi1 i′1 (x1 )Φi2 i′2 (x2 )Φi3 i′3 (x3 ) is a separable function, we now check that the last three terms of the Taylor expansion in (91) lead to a zero sum. In fact, X nk = 0, k = 1, 2, 3, knk32 knk ≤N ∞ n6=0

because the sum contains indexes of the form n and −n which cancel each other. For the other two terms, we write


kxk22 3(x · n)2 1 − x21 2n21 − n22 − n23 + x22 2n22 − n21 − n23 + x23 2n23 − n21 − n22 + = 3 5 5 2knk2 2knk2 2knk2 (92) +6n1 n2 x1 x2 + 6n1 n3 x1 x3 + 6n2 n3 x2 x3 ) and note that

X

knk∞ ≤N n6=0

2n2k − n2k′ − n2k′′ = 0, k, k ′ , k′′ = 1, 2, 3, knk52 X

knk∞ ≤N n6=0

nk nk ′ = 0, k, k ′ = 1, 2, 3, 5 knk2

where we used, for the first sum, an appropriate change of indexes and, for the second sum, we added first over terms of the form nk, − nk . The convergence is conditional since X

knk∞ ≤N n6=0

|nk | , knk32

X

knk∞ ≤N n6=0

n2k , knk52

X

knk∞ ≤N n6=0

|nk nk′ | → ∞ as N → ∞, k, k ′ = 1, 2, 3. knk52

It remains to consider the term 1/knk2 in (91). Due to vanishing moments of Φii′ , Z 1 (93) Φii′ (x)dx = 0, |i + i′ | ≥ 1, knk 2 2 [−1,1]

which finishes the proof of (iii).

For i = i′ = 0, we first use that Φ00 is even and that the sum over n is the same as the sum over −n to rewrite X Z X Z Φ00 (x) Φ00 (x) dx = 8 dx. [0,1]3 kx + nk2 [−1,1]3 kx + nk2 For x ∈

(94)

[0, 1]3

knk∞ ≤N

knk∞ ≤N

we have

Φ00 (x) = (1 − x1 )(1 − x2 )(1 − x3 ) =



1 ϕ(x1 ) + 2

   1 1 ϕ(x2 ) + ϕ(x3 ) + , 2 2

where ϕ(t) = 1/2− t. By expanding the product in (94), we observe that, with the only exception of the term corresponding to the product of the three constants, all other terms satisfy the assumptions 0;0 of Lemma 20 part (2). Hence, the value of T00;00 is Z X 1 1 lim dx1 dx2 dx3 , 4π N →∞ [0,1]3 kx + nk2 knk∞ ≤N

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

27

which, changing variables xj 7→ xj − nj on each j = 1, 2, 3 yields Z Z 1 1 1 1 lim dx = dx = ∞. 4π N →∞ [−N,N +1]3 kxk2 4π R3 kxk2 0;0 Thus, the summation convention (41) yields a non-finite element T00;00 . To deal with this situation, we simply set the value of this element to zero which is equivalent to restrict the domain of the operator to mean-free functions. 

8.5. Auxiliary results for the computation of non-standard form elements. The vanishing moments and symmetries of the cross-correlation functions (20) allow us to explicitly compute elements of the periodized non-standard forms. The relevant properties and how we use them to compute these elements are captured on the following results. Lemma 19. Let ϕ be a bounded function with odd symmetry about 1/2 (95)

ϕ(1 − t) = −ϕ(t), 0 ≤ t ≤ 1.

Then

N Z X

n=−N

1

ϕ(t)h(t + n)dt =

0

Z

1

ϕ(t)h(t + N )dt, 0

for any even function h such that the integrals exist. Proof. Let I be I=

N Z X

n=−N

1

ϕ(t)h(t + n)dt.

0

Splitting the sum in non-negative and negative values of n and changing variables t 7→ 1 − t on the latter, the assumption (95) yields N Z 1 N Z 1 X X I = ϕ(t)h(t + n)dt − ϕ(t)h(1 − t − n)dt =

n=0 0 N Z 1 X n=0 0

ϕ(t)h(t + n)dt −

because h is even.

n=1 0 N Z 1 X n=1 0

ϕ(t)h(t + n − 1)dt =

Z

1

ϕ(t)h(t + N )dt,

0



Lemma 20. Let ϕj , 1 ≤ j ≤ 3, be three bounded functions on [−1, 1] such that one of them, e.g. ϕ1 , is odd and let G(x1 , x2 , x3 ) be a locally integrable function, even on each variable. Then X Z ϕ1 (x1 )ϕ2 (x2 )ϕ3 (x3 )G(x1 + n1 , x2 + n2 , x3 + n3 )dx = 0 (96) lim N →∞

knk∞ ≤N

[−1,1]3

Proof. Let C denote a constant whose value may change along the derivation. Observe that X Z ϕ1 (x1 )ϕ2 (x2 )ϕ3 (x3 )G(x1 + n1 , x2 + n2 , x3 + n3 )dx knk∞ ≤N

[−1,1]3

is always well defined because ϕj are bounded and G is locally integrable. Let g(t) = G(t, x2 + n2 , x3 + n3 ) and isolate the sum over the index n1 and the integral over x1 . Splitting the integral

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

28

over [0, 1] and [−1, 0] and changing variables x1 → −x1 in the latter yields X Z 1 X Z 1 ϕ1 (x1 )g(x1 + n1 )dx1 ϕ1 (x1 )g(x1 + n1 )dx1 = −1

|n1 |≤N

|n1 |≤N

−

0

X Z

|n1 |≤N

1

ϕ1 (x1 )g(−x1 + n1 )dx1 . 0

Since in the last term the sum over n1 is the same as the sum over −n1 and g is an even function, the two terms in the previous equation cancel each other and we obtain the result.  Proposition 21. Let ϕj , 1 ≤ j ≤ 3, denote three bounded functions on [0, 1]. It holds that A: If ϕ1 is odd about 1/2, then (97) lim

N →∞

X

knk∞ ≤N

Z

[0,1]3

ϕ1 (x1 )ϕ2 (x2 )ϕ3 (x3 ) 2π dx = − kx + nk2 3

Z

1

tϕ1 (t) dt

 Z

1

ϕ2 (t) dt



1

ϕ3 (t) dt .

0

0

0

 Z

B: If ϕ1 is even about 1/2 and mean free, then Z X Z ϕ1 (x1 ) 4π 1 2 (98) lim dx = t ϕ1 (t) dt. N →∞ 3 0 [0,1]3 kx + nk2 knk∞ ≤N

C: If ϕ1 is mean free, then X Z (99) lim N →∞

[0,1]3

knk∞ ≤N

ϕ1 (x1 ) dx = −2π kx + nk2

Z

1 0

4π tϕ1 (t) dt + 3

Z

1

t2 ϕ1 (t) dt.

0

For simplicity, the proposition is stated for the Poisson kernel G(x) = kxk−1 , but similar results hold for any radially symmetric kernel with enough decay at infinity and, thus, to linear combination of such kernels. However, due to the slow decay of the Poisson kernel, the proof of Proposition 21 is more challenging than the one for kernels with faster decay at infinity. Proof. We use the same notation as in the proof of Lemma 20. Note that, the same argument given in that proof shows that X Z ϕ1 (x1 )ϕ2 (x2 )ϕ3 (x3 ) + p (100) SN = dx. 3 (x1 + n1 )2 + (x2 + n2 )2 + (x3 + n3 )2 knk ≤N [0,1] ∞

is well-defined for all N . To prove part A, we use Lemma 19 with h(t) = G(t, x2 + n2 , x3 + n3 ) to write (101) Z 1 X Z 1 ϕ1 (x1 )G(x1 + N, x2 + n2 , x3 + n3 )dx1 . ϕ1 (x1 )G(x1 + n1 , x2 + n2 , x3 + n3 )dx1 = IN = |n1 |≤N

0

0

The assumption of odd symmetry for ϕ1 implies vanishes at the endpoints of [0, 1], [1]

(102)

R1 0

[1]

ϕ1 (t)dt = 0 and, thus, ϕ1 (x) =

[1]

ϕ1 (0) = ϕ1 (1) = 0.

Integrating by parts the last term of the identity (101) and using (102), we obtain Z 1 [1] (x1 + N )ϕ1 (x1 )G(x1 + N, x2 + n2 , x3 + n3 )3 dx1 . (103) IN = 0

Rx 0

ϕ1 (t)dt,

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

29

Hence, substituting (103) into (100) yields Z [1] X (x1 + N )ϕ1 (x1 )ϕ2 (x2 )ϕ3 (x3 ) + dx. (104) SN = 3 ((x + N )2 + (x + n )2 + (x + n )2 )3/2 [0,1] 1 2 2 3 3 |n2 |≤N,|n3 |≤N Since (x1 + N )2 + (x2 + n2 )2 + (x3 + n3 )2

(105)

[1]


−3/2

≤ N −3 ,

in the limit for N → ∞, only the numerator term N ϕ1 (x1 )ϕ2 (x2 )ϕ3 (x3 ) provides a non-zero contribution in (100) and hence Z [1] X ϕ1 (x1 )ϕ2 (x2 )ϕ3 (x3 ) + + dx. (106) S∞ = lim SN = lim N 2 2 2 3/2 N →∞ N →∞ 3 |n2 |≤N,|n3 |≤N [0,1] ((x1 + N ) + (x2 + n2 ) + (x3 + n3 ) ) Isolating again the integral with respect to x1 and integrating by parts on that variable, we obtain Z 1 1 [1] [2] ϕ1 (x1 )G(x1 + N, x2 + n2 , x3 + n3 )3 dx1 = ϕ1 (x1 )G(x1 + N, x2 + n2 , x3 + n3 )3 0

0

(107)

+

Z

1

0

[2]

where ϕ1 (x) =

Rx

((x1 +

[2] 3(x1 + N )ϕ1 (x1 ) N )2 + (x2 + n2 )2 + (x3

+ n3 )2 )5/2

dx1 ,

[1]

ϕ1 (t)dt. Since the integrand of the right hand side is bounded by
−5/2 C(N + 1) (x1 + N )2 + (x2 + n2 )2 + (x3 + n3 )2 ≤ C(N + 1)N −5 , 0

+ is bounded by the contribution of this term to the value of S∞ X N (N + 1) (N + 1) (2N + 1)2 C 1 = C N5 N4 |n2 |≤N,|n3 |≤N

which vanishes as N → ∞. It follows that Z X [2] + S∞ = ϕ1 (1) lim N N →∞

|n2 |≤N,|n3 |≤N

ϕ2 (x2 )ϕ3 (x3 )

[0,1]2

((1 +

N )2

+ (x2 + n2 )2 + (x3 + n3 )2 )3/2

dx.

By the same argument, this time integrating by parts first respect to x2 and then respect to x3 , we obtain X 1 [2] [1] [1] + S∞ = ϕ1 (1)ϕ2 (1)ϕ3 (1) lim N 2 2 2 3/2 N →∞ |n2 |≤N,|n3 |≤N ((1 + N ) + (1 + n2 ) + (1 + n3 ) ) [1]

where ϕ2 (x) =

Rx

[1]

ϕ2 (t)dt and ϕ3 (x) = X I = lim N 0

N →∞

=

lim N

N →∞

Rx 0

|n2 |≤N,|n3 |≤N N +1 X

ϕ3 (t)dt. Let us denote 1 ((1 + N )2 + (1 + n2 )2 + (1 + n3 )2 )3/2

N +1 X

n2 =−N +1 n2 =−N +1

1 (1

+ N )2

+ n22 + n23


3/2 .

Observe that we may consider the sums in the range |n2 | ≤ N + 1, |n3 | ≤ N + 1 because the limit is zero when n2 or n3 are set to −N or −N − 1. Therefore, X 1 1 I = lim 3/2  2 N →∞ (N + 1) n n |n2 |≤N +1,|n3 |≤N +1 1 + ( 2 )2 + ( 3 )2 N +1 N +1

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

30

and hence I is the limit of Riemann sums for the continuous function (1+x2 +y 2 )−3/2 in the interval [−1, 1]2 , yielding Z Z 1 2 p dy = π. (1 + x2 + y 2 )−3/2 dxdy = 2 I= 2 2 3 2 2+y [−1,1] (1 + y ) [−1,1]

1 R R1 [1] [2] [1] 1 To finish the proof, note that ϕ1 (1) = tϕ1 (t) − 0 tϕ1 (t)dt = − 0 tϕ1 (t)dt and ϕj (1) = 0 R1 0 ϕ(t)dt. For part B, let us denote ϕ = ϕ1 . Since t − 1/2 is odd about 1/2 we have Z 1 Z 1 1 (t − )ϕ(t)dt = (108) 0= tϕ(t)dt, 2 0 0

since, by assumption, ϕ is even about 1/2 and mean-free. Denoting the successive anti-derivatives of ϕ by Z t ϕ[j−1] (s)ds, ϕ[j] (t) = 0

where

ϕ[0] (t)

= ϕ(t), we observe that the mean free property of ϕ yields ϕ[1] (0) = ϕ[1] (1) = 0.

(109)

Also, integration by parts and (108) yields [2]

(110)

ϕ (1) = −

Z

1

tϕ(t)dt = 0 = ϕ[2] (0).

0

Thus, ϕ[1] and ϕ(2) vanish at the endpoints of [0, 1] and hence  Z  Z 1 1 2 1 1 d2 2 [2] s ϕ (s)ds = s ϕ(s)ds. (111) ϕ[3] (1) = 2 0 ds2 2 0

Similarly to the proof of part A, we use these properties of the anti-derivatives of ϕ to show that, in order to compute X Z ϕ(x1 ) + + S∞ = lim SN = lim dx, N →∞ N →∞ [−1,1]3 kx + nk2 knk∞ ≤N

it is enough to consider the sums over n2 and n3 within the range −N, N − 1 instead of −N, N . In fact, let’s consider the term n3 = N and integrate by parts. Using (109) we obtain Z Z ϕ(x1 ) ϕ[1] (x1 )(x1 + n1 ) dx = dx, 1/2 3/2 [0,1]3 ((x1 + n1 )2 + (x2 + n2 )2 + (x3 + N )2 ) [0,1]3 ((x1 + n1 )2 + (x2 + n2 )2 + (x3 + N )2 ) which, now using (110), equals Z ϕ[2] (x1 ) [0,1]3

((x1 + n1 )2 + (x2 + n2 )2 + (x3 + N )2 )3/2

Since

and

X

+

ϕ[2] (x1 )(x1 + n1 )2 ((x1 + N )2 + (x2 + n2 )2 + (x3 + N )2 )5/2

(x1 + n1 )2 + (x2 + n2 )2 + (x3 + N )2

|n1 |≤N,|n2 |≤N

X

|n1 |≤N,|n2 |≤N

(x1 +n1 )2 (x1 + n1 )2 + (x2 + n2 )2 + (x3 + N )2


−3/2


−5/2

dx.

≤ (2N + 1)2 N −3

≤ (2N +1)N −5

X

|n1 |≤N

(1+n1 )2 ≤ cN −1 ,

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

31

+ leads to a sequence which tends to 0 as N → ∞. Setting the term corresponding to n3 = N in SN n2 = N leads to a similar estimate yielding + S∞ = lim SN , N →∞

where SN =

N −1 X

N X

N −1 X

n1 =−N n2 =−N n3 =−N

Z

[−1,1]3

ϕ(x1 ) dx. kx + nk2

Changing variables xj 7→ xj − nj for j = 2, 3 and combining the sums with the integrals, we obtain Z NZ N X Z n1 +1 1 p dx2 dx3 dx1 . SN = 4 ϕ(x1 − n1 ) 2 x1 + x22 + x23 0 0 |n |≤N n1 1

We now explicitly compute the integrals over x2 and x3 and denote the result by ! Z N Z NZ N 1 N p dx3 arcsinh p dx2 dx3 = aN (x) = x2 + x22 + x23 x2 + x23 0 0 0     N N2 √ = 2N arcsinh √ . − xarctan x2 + N 2 x x2 + 2N 2

Observe that aN has particularly simple derivatives,

√ x x2 + 2N 2 d aN (x) = −arccot dx N2 d2 2N 2 √ a (x) = N dx2 (x2 + N 2 ) x2 + 2N 2
2 5N 4 x + 3N 2 x3 d3 aN (x) = − 3 . dx3 (x2 + N 2 )2 (x2 + 2N 2 ) 2

Hence, using (109) and (110), integration by parts yields, X Z n1 +1 d2 X Z n1 +1 d2 [2] [2] SN = 4 ϕ (x − n )a (x )dx = 4 ϕ (x − n ) aN (x1 )dx1 1 1 N 1 1 1 1 dx2 dx2 |n1 |≤N n1 |n1 |≤N n1 X X Z n1 +1 d2 d3 =4 ϕ[3] (1) 2 aN (n1 + 1) − 4 (112) ϕ[3] (x1 − n1 ) 3 aN (x1 )dx1 , dx dx n1 |n1 |≤N

|n1 |≤N

because ϕ[3] (1) = 0. The last term in (112), vanishes as N → ∞ since Z n1 +1 X d3 3 X d [3] [3] ϕ (x1 − n1 ) 3 aN (x1 )dx1 ≤ max ϕ (t) aN (x1 ) 3 dx dx t∈[0,1] |n1 |≤N n1 |n1 |≤N X ≤C N −3 = C(2N + 1)N −3 , |n1 |≤N

for some constant C. Therefore, using (111), lim SN = 2

N →∞

Z

1

2

t ϕ(t)dt lim 0

N →∞

N +1 X

n=−N +1

d2 aN (n) dx2

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

32

and, since limN →∞ aN (N + 1) = limN →∞ aN (−N ) = 0, Z 1 N 1 X d2 2 t ϕ(t)dt lim lim SN = 4 aN (n). N →∞ 2 N →∞ dx2 0 n=−N

The result follows observing that N N 1 X 1 X d2 aN (n) = 2 dx2 N

(113)

n=−N

n=−N

1 q
n 2 ) 2+ (1 + N


n 2 N

it is a Riemann Sum in the interval [−1, 1] for the continuous function the sum (113) converges to

Z

1 −1

(1 +

1 √

x2 )

2

+ x2

dx =

1√ . (1+x2 ) 2+x2

As N → ∞,

π . 3

For part C, given a mean free function ϕ1 we write it as ϕ1 (t) = ϕodd (t) + ϕeven (t), where ϕ1 (t) − ϕ1 (1 − t) ϕ1 (t) + ϕ1 (1 − t) and ϕeven (t) = . 2 2 Since both ϕ1 and ϕodd are mean free, the same holds for ϕeven . Using parts A and B and the definitions of ϕodd and ϕeven , the result follows adding  Z 1 Z 1 Z 1 2 2 1 2 tϕodd (t)dt = − t ϕ1 (t)dt ϕ1 (t)dt = − t− − 3π 0 3π 0 2 3π 0

(114)

ϕodd (t) =

and 4 3π

Z

1

0

4 t ϕeven (t)dt = 3π 2

Z

0

1

1 t −t+ 2 2



ϕ1 (t)dt =

4 3π

Z

1

0


t2 − t ϕ1 (t)dt.



8.6. Proof of Proposition 7. ′

Proof. Since from (72) we have that Φi′ i (x) = (−1)i+i Φii′ (x), it is enough to show the result for Tii0;0 ′ ;00 . By (26) and (22), Tii0;0 ′ ;00 (µ)

= lim

N →∞

X

Tii0;n ′ ;00 (µ)

knk∞ ≤N

= lim

N →∞

X

knk∞ ≤N

Z

[−1,1]3

Gµfree (x+n)Φi1 i′1 (x1 )Φi2 i′2 (x2 )Φi3 i′3 (x3 )dx.

Therefore, Lemma 20 implies that Tii0;0 ′ ;00 (µ) vanishes whenever any of the functions Φij i′ , j = 1, 2, 3 j is odd, which, by (73), is the case if ij and ij ′ have different parity. We have proved part (ii). Next consider the case of ij and ij ′ having the same parity for all j. In this case all the functions Φij i′j , j = 1, 2, 3 are even and X Z 0;0 Gµfree (x + n)Φi1 i′1 (x1 )Φi2 i′2 (x2 )Φi3 i′3 (x3 )dx Tii′ ;00 (µ) = 8 lim N →∞

=8

Z

knk∞ ≤N

[0,1]3

[0,1]3

Φi1 i′1 (x1 )Φi2 i′2 (x2 )Φi3 i′3 (x3 )GµH (x)dx.

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

33

For part (iii), by symmetry of the kernel, it is sufficient to consider only one of the elements listed on each of the three cases. The case i = (1, 1, 0) and i′ = (0, 0, 0) follows from part (ii). For the other two cases, we write Φii′ (x) = Φ(x1 )Φ00 (x2 )Φ00 (x3 ) = Φ(x1 )(1 − |x2 |)(1 − |x3 |), where Φ is either Φ11 or Φ20 . Part (i) is also covered considering Φ(x1 ) = 1 − |x1 |. We write    Z 1 1 0;0 Φ(x1 ) (115) Tii′ ;00 (µ) = 8 + ϕ(x2 ) + ϕ(x3 ) GµH (x)dx, 2 2 [0,1]3

where ϕ(t) = 1/2 − t. Expanding the product in (115) and using that GµH is even about 1/2 on each variable, GµH (1 − x1 , 1 − x2 , 1 − x3 ) = GµH (x1 , x2 , x3 ), we obtain that all terms vanish, with the only exception of the term corresponding to 1/4 Φ(x1 ). Hence, Z 0;0 Φ(x1 )GµH (x)dx. (116) Tii′ ;00 = 2 [0,1]3

√
For Φ(x1 ) = Φ20 (x1 ) = − 5x1 1 − 3x1 + 2x21 , Tii0;0 ′ ;00 (µ) vanishes because Φ is odd about 1/2. For Φ(x1 ) = 1 − x1 = 1/2 + ϕ(x1 ), we obtain part (i) since Z Z ∞ Z 1 1 e−µkxk µ 0;0 GH (x)dx = T00;00(µ) = e−µr r dr = 2 . dx = 4π R3 kxk µ [0,1]3 0

It only remains to consider Φ(x1 ) = Φ11 (x1 ) = 1 − 3x1 + 2x31 in (116). This case is more delicate and to obtain the answer we describe a more general approach for the computation of integrals of the form Z I=

B

ϕ1 (x1 )ϕ2 (x2 )ϕ3 (x3 )GµH (x)dx

where ϕj ∈ L2 ([0, 1]3 ). Since Gµfree ∈ L1 (R3 ), it follows from [33, Thm. 2.4] that GµH ∈ L1 ([0, 1]3 ) ˆ 
and that its Fourier coefficients are given by gn = Gµfree (n) = 1/ 4π 2 n21 + 4π 2 n22 + 4π 2 n23 + µ2 . Hence {gn }n∈Z3 ∈ l2 (Z3 ), and, therefore GµH ∈ L2 ([0, 1]3 ). As a result, using Parseval’s identity we obtain ϕ b1,n1 ϕ b2,n2 ϕ b3,n3 1 X I = 2 µ
2 2 2 2 4π 3 n +n + n + n∈Z

1

1 = 2 lim 4π N →∞

(117)

2

3

X

knk∞ ≤N

2π

ϕ b1,n1 ϕ b2,n2 ϕ b3,n3

n21 + n22 + n23 +

µ
2 2π

,

where ϕ bj,nj is the Fourier coefficient nj of the function ϕj . In particular, if ϕ2 ≡ ϕ3 ≡ 1, using Parseval’s identity but now for functions of one variable, we obtain Z 1   1 X 1 ϕ b1,n1 1 X \ −µ|x| I= 2 ϕ b1,n1 e = ϕ1 (t)Aµ (t)dt,
2 = 4π 2µ 2µ 0 n n2 + µ n∈Z

n∈Z

2π

where

Aµ (t) =

X

n∈Z

e−µ|t+n| =

eµt + eµ(1−t) , t ∈ [0, 1], eµ − 1

is odd about 1/2. Writing ϕ1 (t) = ϕodd + ϕeven (t) as in (114), we have Z 1 Z 1 1 1 I= ϕeven (t)Aµ (t)dt = ϕeven (t)eµt dt. 2µ 0 µ (eµ − 1) 0

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

34

Therefore, for i = i′ = (1, 0, 0) we obtain

Z 1 12 − 6µ + µ2 eµ − 12 + 6µ + µ2 1 − 6t + 6t2 µt 2 0;0 e dt = . Tii′ ;00 (µ) = µ (eµ − 1) 0 2 µ4 (eµ − 1)  8.7. Proof of Proposition 10. Proof. For part (i), observe that the separable function Φii′ (x) contains a function of the form Φij i′

j

′

such that ij + ij is odd. By (73), such a function is odd. Hence, the result follows from Lemma 20. ′

For part (ii), since (72) asserts that Φi′ i (x) = (−1)i+i Φii′ (x) , it is enough to show the result for Tii0;0 ′ ;00 . By symmetry of the kernel, it is sufficient to consider only one of the elements listed on each of the three cases. The case i = (1, 1, 0) and i′ = (0, 0, 0) follows from part (i). For the other two cases, we write Φii′ (x) = Φ(x1 )Φ00 (x2 )Φ00 (x3 ) = Φ(x1 )(1 − |x2 |)(1 − |x3 |),

where Φ is either Φ11 or Φ20 . By (73), Φii′ (x) is even in all of its arguments and hence

X Z Φ(x1 ) ϕ(x2 ) + 12 ϕ(x3 ) + 21 2 0;0 (118) Tii′ ;00 = lim dx1 dx2 dx3 , π N →∞ kx + nk2 [0,1]3 knk∞ ≤N

where ϕ(t) = 1/2 − t. By expanding the product in (118), we observe that, with the only exception of the term corresponding to 1/4 Φ(x1 ), all other terms satisfy the assumptions of Proposition 21 part A. These terms do not contribute to the value of Tii0;0 ′ ;00 because all of them contain the factor R1 0 Φ(t)dt, which is zero due to (3) and that Φ is an even function, Z Z 1 1 1 Φ(t)dt = 0. Φ(t)dt = 2 −1 0 Hence, Tii0;0 ′ ;00 equals

1 lim 2π N →∞

X

knk∞ ≤N

Z

[0,1]3

Φ(x1 ) dx1 dx2 dx3 . kx + nk2

The result follows using √ Proposition 21 part
C applied to either Φ(x1 ) = Φ11 (x1 ) = 1 − 3x1 + 2x31  or Φ(x1 ) = Φ20 (x1 ) = − 5x1 1 − 3x1 + 2x21 with x1 ∈ [0, 1]. 8.8. Proof of Remark 11.

′

Proof. Taking into account part (i) of Proposition 10, it is enough to consider ij + ij to be even, for all j = 1, 2, 3 and hence (73) yields Φi1 0 (x1 )Φi2 i′ (x1 )Φi3 i′ (x1 ) X Z 2 0;0 2 3 Tii′ ;00 = lim dx1 dx2 dx3 π N →∞ kx + nk2 [0,1]3 knk∞ ≤N

Since i1 is even and greater than one, (74) implies that Φi1 0 (x1 ) is odd about 1/2. Thus, using Proposition 21 part A, we have   Z 1  Z 1 Z 1 4 0;0 Φi2 i′ (t) dt , Φi2 i′ (t) dt tΦi1 0 (t) dt Tii′ ;00 = − 2 2 3 0 0 0

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

35

which, by (80), vanishes if i1 ≥ 4. It remains to consider the case i1 = 2. Since |i + i′ | > 2, either ′ ′ i2 + i2 or i3 + i3 is positive; thus, due to (3) and that the functions Φij i′ are even, at least one of j the two functions Φi2 i′ or Φi2 i′ is mean-free.  2

2

8.9. Proof of Proposition 12. Proof. Let SN =

X

knk∞ ≤N

Z

[0,1]3

ϕ1 (x1 )ϕ2 (x2 )ϕ3 (x3 )G(x1 + n1 , x2 + n2 , x3 + n3 )dx,

and assume that ϕ1 is odd about 1/2. Repeating the steps performed at the beginning of the proof of Proposition 21, we obtain Z 1 X Z 1 [1] ϕ1 (x1 )Gx1 (x1 + N, x2 + n2 , x3 + n3 )dx1 , ϕ1 (x1 )G(x1 + n1 , x2 + n2 , x3 + n3 )dx1 = − |n1 |≤N

0

0

[1]

where ϕ1 (x) =

Rx

ϕ1 (t)dt. Hence, by the assumption (42), we have
−(2+δ)/2 |SN | ≤ C(2N + 1)2 (x1 + N )2 + (x2 + n2 )2 + (x3 + n3 )2 ≤ C(2N + 1)2 N −2−δ , 0

where N ≥ M and C is a constant. The result follows.



References [1] B. Alpert. A class of bases in L2 for the sparse representation of integral operators. SIAM J. Math. Anal, 24(1):246–262, 1993. [2] B. Alpert, G. Beylkin, R. Coifman, and V. Rokhlin. Wavelet-like bases for the fast solution of second-kind integral equations. SIAM J. Sci. Comput., 14(1):159–184, 1993. [3] B. Alpert, G. Beylkin, D. Gines, and L. Vozovoi. Adaptive solution of partial differential equations in multiwavelet bases. J. Comput. Phys., 182(1):149–190, 2002. [4] J. Batt and G. Rein. A rigorous stability result for the Vlasov-Poisson system in three dimensions. Ann. Mat. Pura Appl., 164(1):133–154, 1993. [5] P. Becker and P. Coppens. About the Coulomb potential in crystals. Acta Crystallographica Section A, 46(4):254– 258, Apr 1990. [6] G. Beylkin, V. Cheruvu, and F. Pérez. Fast adaptive algorithms in the non-standard form for multidimensional problems. Appl. Comput. Harmon. Anal., 24(3):354–377, 2008. [7] G. Beylkin, R. Coifman, and V. Rokhlin. Fast wavelet transforms and numerical algorithms, I. Comm. Pure Appl. Math., 44(2):141–183, 1991. Yale Univ. Technical Report YALEU/DCS/RR-696, August 1989. [8] G. Beylkin and R. Cramer. A multiresolution approach to regularization of singular operators and fast summation. SIAM J. Sci. Comput., 24(1):81–117, 2002. [9] G. Beylkin, R. Cramer, G.I. Fann, and R.J. Harrison. Multiresolution separated representations of singular and weakly singular operators. Appl. Comput. Harmon. Anal., 23(2):235–253, 2007. [10] G. Beylkin, C. Kurcz, and L. Monzón. Fast algorithms for Helmholtz Green’s functions. Proceedings of the Royal Society A, 464(2100):3301–3326, 2008. doi:10.1098/rspa.2008.0161. [11] G. Beylkin and M. J. Mohlenkamp. Numerical operator calculus in higher dimensions. Proc. Natl. Acad. Sci. USA, 99(16):10246–10251, August 2002. [12] G. Beylkin and M. J. Mohlenkamp. Algorithms for numerical analysis in high dimensions. SIAM J. Sci. Comput., 26(6):2133–2159, July 2005. [13] G. Beylkin and L. Monzón. On approximation of functions by exponential sums. Appl. Comput. Harmon. Anal., 19(1):17–48, 2005. [14] G. Beylkin and L. Monzón. Approximation of functions by exponential sums revisited. Appl. Comput. Harmon. Anal., 28(2):131–149, 2010. [15] J. Callaway and M. Glasser. Fourier coefficients of crystal potentials. Phys. Rev., 112(1):73–77, Oct 1958. [16] M. Challacombe, C. White, and M. Head-Gordon. Periodic boundary conditions and the fast multipole method. The Journal of Chemical Physics, 107(23):10131–10140, 1997. [17] A. Chorin and J. Marsden. A mathematical introduction to fluid mechanics. Springer, third edition, 1993.

MULTIRESOLUTION REPRESENTATION OF OPERATORS WITH BOUNDARY CONDITIONS

36

[18] S. de Leeuw, J. Perram, and E. Smith. Simulation of electrostatic systems in periodic boundary conditions I: Lattice sums and dielectric constants. Proceedings of the Royal Society of London Series A, Mathematical and Physical Sciences, 373(1752):27–56, 1980. doi:10.1098/rspa.1980.0135. [19] H. M. Evjen. On the stability of certain heteropolar crystals. Phys. Rev., 39(4):675–687, Feb 1932. [20] P. Ewald. Die berechnung optischer und elektrostatischer gitterpotentiale. Ann. Phys., 64:253–287, 1921. [21] M.L. Glasser and I.J. Zucker. Lattice sums. In H. Eyring and D. Henderson, editors, Theoretical Chemistry: Advances and Perspectives, volume 5, pages 67–139. Academic Press, 1980. [22] L. Greengard and M. C. Kropinski. Integral equation methods for Stokes flow in doubly-periodic domains. J. Engrg. Math., 48(2):157–170, 2004. [23] L. Greengard and V. Rokhlin. A fast algorithm for particle simulations. J. Comp. Phys., 73(1):325–348, 1987. [24] R.J. Harrison, G.I. Fann, T. Yanai, Z. Gan, and G. Beylkin. Multiresolution quantum chemistry: basic theory and initial applications. J. Chem. Phys., 121(23):11587–11598, 2004. [25] T. Kato and G. Ponce. Commutator estimates and the Euler and Navier-Stokes equations. Comm. Pure Appl. Math, XLI:891–907, 1988. [26] E. V. Kholopov. Convergence problems of coulomb and multipole sums in crystals. Physics-Uspekhi, 47(10):965– 990, 2004. [27] C. Kittel. Introduction to Solid State Physics. John Wiley and Sons Inc., sixth edition, 1986. [28] D. Lindbo and A-K. Tornberg. Spectrally accurate fast summation for periodic Stokes potentials. J. Comput. Phys., 229(23):8994–9010, 2010. [29] M. O’Keeffe and J. C. H. Spence. On the average Coulomb potential and constraints on the electron density in crystals. Acta Crystallographica Section A, 50(1):33–45, Jan 1994. [30] C. G. Poulton, L. C. Botten, R. C. McPhedran, and A. B. Movchan. Source-neutral Green’s functions for periodic problems in electrostatics, and their equivalents in electromagnetism. R. Soc. Lond. Proc. Ser. A Math. Phys. Eng. Sci., 455:1107–1123, 1999. [31] D. Rapaport. The art of molecular dynamics simulation. Cambridge University Press, 2004. [32] Lord Rayleigh. On the influence of obstacles arranged in a rectangular order upon the properties of a medium. Philos. Mag., 34:481, 1892. [33] E. Stein and G. Weiss. Fourier Analysis on Euclidean Spaces. Princeton Univ. Press, Princeton NJ, 1971. Current address: * Department of Applied Mathematics, University of Colorado at Boulder, 526 UCB, Boulder, CO 80309-0526, USA., + Oak Ridge National Laboratory, P.O. Box 2008, MS6367, Oak Ridge, TN 37831, USA