Convergence analysis of a high-order Nyström integral-equation method for surface scattering problems

# Convergence analysis of a high-order Nyström integral-equation method for surface scattering problems

Oscar P. Bruno1, Víctor Domínguez2  & Francisco-Javier Sayas3
11Applied and Computational Mathematics, California Institute of Technology, Pasadena, CA 91125, USA – bruno@acm.caltech.edu
22 Dep. Ingeniería Matemática e Informática, E.T.S.I.I.T., Universidad Pública de Navarra. 31500 - Tudela, Spain – victor.dominguez@unavarra.es
33 Department of Mathematical Sciences, University of Delaware, Newark, DE 19716, USA – fjsayas@math.udel.edu
September 10, 2019
###### Abstract

In this paper we present a convergence analysis for the Nyström method proposed in [Jour. Comput. Phys. 169 pp. 2921-2934, 2001] for the solution of the combined boundary integral equation formulations of sound-soft acoustic scattering problems in three-dimensional space. This fast and efficient scheme combines FFT techniques and a polar change of variables that cancels out the kernel singularity. We establish the stability of the algorithms in the norm and we derive convergence estimates in both the and norms. In particular, our analysis establishes theoretically the previously observed super-algebraic convergence of the method in cases in which the right-hand side is smooth.

Key words: Acoustic Scattering, Boundary Integral Equations, Nyström methods, FFT
MSC: 65N38, 65N35, 65T40

## 1 Introduction

In this paper we present a full convergence analysis of the high-order Nyström method introduced by Bruno and Kunyansky in [6, 7] for boundary integral equations (BIE) related to the scattering of time-harmonic acoustic waves by arbitrary (smooth) surfaces in three-dimensional space. In particular, our analysis establishes theoretically the previously observed super-algebraic convergence of the method for smooth right-hand sides. To the best of the authors’ knowledge, this proof constitutes one of the few instances of analysis of a Nyström type method for a boundary integral equation in three dimensions.

As is well known, Galerkin methods for BIE of the first kind have enjoyed thorough theoretical analyses since their inception—on the basis of ellipticity properties and discrete Fredholm theory. Compactness arguments can also be used to establish convergence of Galerkin methods for equations of the second kind. Few results exist, on the other hand, concerning convergence for three-dimensional BIE collocation methods—in which finite-element bases are used for approximation, but testing relies on point sampling. We refer to the excellent text-book [3, Chapter 9] for a brief introduction on this topic.

This state of affairs has led to the widespread perception that, being even more “discrete” than collocation schemes, Nyström methods for BIEs of the second kind with weakly singular kernels could not be easily analyzed. This paper will hopefully help dispell this view and lend additional credibility to Nyström methods—whose qualities can be very attractive in practice [6].

One of the main difficulties in the design of three-dimensional integral solvers concerns development of high-quality quadrature rules for approximation of singular integral terms. Wienert [24] constructed a singular integration rule on the basis of spherical-harmonic transforms for surfaces for which a diffeomorphism to the sphere can be constructed, (see also [12, Section 3.6]). A Galerkin version of this approach was introduced and analyzed in [16] where, in particular, the superalgebraic convergence of the method was established. Reference [15] presents a collocation method, with corresponding analysis, for the Laplace equation, which shares the good convergence properties of the method by Wienert but which again is limited to surfaces for which a smooth mapping from the sphere is available. We emphasize that such limitations, which are highly restrictive in practice, are not imposed by the Nyström method studied in this paper.

We thus consider the numerical method [6] for the second-kind combined field integral equation associated to the problem of sound-soft scattering by a three dimensional obstacle with smooth boundary . The method relies on a series of geometric constructions: 1) Representation of the surface by means of a set of overlapping smooth charts (parametrizations); 2) A smooth partition of unity subordinated to these charts which decomposes both the overall integral operators as well as the solution of the BIE as the sum of contributions defined on the parametrized patches; 3) A family of floating smooth cut-off functions (see (2.10) below) that is used to isolate the singularity of the kernel function, and thus produces a splitting of the integral operator of the form —whereby the regular operator is an integral operator with a smooth kernel, and the singular part which enjoys a reduced domain of integration but contains the weak singularity of the integral kernel.

The Nyström method under consideration is based on approximation of the integral operator by a quadrature rule which treats the regular and singular parts and separately. The local charts mentioned above are used to push forward uniform grids from the unit square to the surface; these grids are used for both, approximation and integration. The quadrature rule used for the regular operator is based on application of the two-dimensional trapezoidal rule in the parameter space for each one of the parametric domains: since the corresponding integrands are smooth with compact support strictly contained in the unit square, these trapezoidal quadratures give rise to super-algebraically accurate approximations. For the singular part, a change of variables to polar coordinates around each integration point is applied. This procedure results in a smooth integrand to which, upon necessary Fourier-based interpolations, the trapezoidal rule is applied for radial and angular integration—once again yielding super-algebraically accurate approximations.

As a result of these constructions we obtain an algorithm that evaluates the action of an approximating discrete operator on a continuous function, using only the values of the function on the selected quadrature points.

Our theoretical treatment relies on use of both existing and new analytical tools. In a first key step of the analysis the problem is re-expressed as a system of integral equations in a space of periodic functions. This is accomplished by means of yet another set of local cut-off functions, whose presence does not affect the actual numerical scheme. The use of periodic Sobolev spaces allows us to take full advantage of numerous results for approximation of Fourier series and interpolation operators [22], as well as the theory of collective compact operators by Anselone [2]. Recasting the numerical scheme of quadrature type as a discretization method in -type (Sobolev) spaces gives rise to a number of difficulties. In particular, for the sake of the analysis we introduce Fourier projection operators, bi-periodic trigonometric interpolation operators, and a discrete operator [10] that produces a linear combination of Dirac delta distributions on a uniform grid from a smooth function input. The convergence analysis for the operators arising from the regular part of the original boundary integral operator follows the lines of the theory [10, 14] on periodic integral equations in one variable. The final (rather technical) element of our convergence proof is a detailed analysis of the integration error arising from the numerical polar coordinate integration algorithm for products of smooth functions, “shrinking” floating cut-off functions and bi-variate trigonometric polynomials, in terms of Sobolev norms of the latter polynomials.

We point out that, for efficiency, a variety of acceleration techniques were used in [6, 7] in conjunction with the Nyström algorithm we consider. Some of these algorithmic components have been taken into account in our analysis. In the formulation considered in the present paper, for example, the computation of the singular part requires the evaluation of bivariate trigonometric polynomials that approximate the unknown at points on a polar grid. A deliberate choice of the radial quadrature nodes for this integral makes it possible to reduce this process to 1–dimensional trigonometric interpolation problems (see section 2.3 and Figure 3 below) on the horizontal and vertical lines of the grid. In [6], such trigonometric polynomials are then approximated by means of piecewise Hermite interpolation on dyadic grids, which can be evaluated much more rapidly than the either of the underlying trigonometric polynomials. We have analyzed the effect of these additional approximations on the full convergence of the Nyström method. The corresponding results can be found in Appendix B; briefly, upon correct parameter selections, the resulting (more efficient) method retains the super-algebraic convergence of the original approach. Additional, more sophisticated acceleration techniques which, based on use of equivalent sources and Fast Fourier Transforms, provide a means to reduce the solution cost for high-frequency problems (but on which the Nyström method itself does not depend, and whose use is not advantageous for problems of lower frequencies) were introduced in [6]. The impact of such equivalent-source acceleration methods on convergence are not considered either in the present paper or in the Appendix B of this paper.

This paper is structured as follows. In Section 2 we describe, in a compact form, the Nyström method under consideration, and we state the main convergence results of this paper. In Section 3 we then recast both the continuous and discrete problems as systems of equations in spaces of biperiodic functions, we derive bounds, on various norms, of the main continuous integral operators in our periodic formulation, we establish unique solvability of the continuous system of periodic integral equations as well as the equivalence of this system to the original BIE, and, finally, we state the main approximation results in the biperiodic frame: norm convergence of the discrete operators to the continuous ones together with corresponding error estimates. In Section 4 we present key estimates on the singular operators that appear in the biperiodic formulation, and in Section 5, in turn, we provide the proofs of the main results stated in Section 3 and Section 2—in that order. Appendix A is devoted to the Sobolev error analysis mentioned above for the polar-integration of products of smooth functions, cut-offs and trigonometric polynomials. Finally, in Appendix B we describe and analyze a slight variant of the numerical method, where one of the interpolation processes is substituted by polynomial interpolation.

We conclude our introduction with two remarks concerning notation.

###### Remark 1.1

To make a clearly visible distinction between points on the surface and coordinates for their parametrization, we use boldface letters (e.g. ) for points on the scattering surface , and underlined letters (such as ) for points in . The coordinates of such points will be denoted according to .

###### Remark 1.2

Throughout this paper the letter denotes a positive constant independent of the parameters and and any other variable quantities appearing in the equation. When necessary a subscripted letter is used—either to avoid confusion or to explicitly signify dependence on parameters other than and .

## 2 The Nyström method

### 2.1 The Boundary Integral Equation

We consider the problem of time-harmonic acoustic scattering by a sound-soft obstacle with smooth boundary in three-dimensional space:

 ⎡⎢⎣ΔU+κ2U=0in ext(S)% ,U=−Uincon S,∂rU−iκU=o(r−1)as r→∞. (2.1)

Here is the wave number, is the incident wave, and denotes the radial derivative. Letting denote the fundamental solution of the Helmholtz equation,

 Φκ(r,r′):=exp(iκ|r−r′|)4π|r−r′|,

then the solution of (2.1) can be expressed as the combined (or Brakhage-Werner [4]) potential

 U(r)=∫S∂Φκ(r,r′)∂ν(r′)ψ(r′)dr′−iη∫SΦκ(r,r′)ψ(r′)dr′,

where denotes the outward normal derivative on and is a coupling parameter. The density is the unique solution of the integral equation

 12ψ(r)+Aψ(r)=−Uinc(r)∀r∈S, (2.2)

where

 Aψ(r)=∫SK(r,r′)ψ(r′)dr′:=∫S(∂Φκ(r,r′)∂ν(r′)−iηΦκ(r,r′))ψ(r′)dr′. (2.3)

Various choices of the coupling parameter have been proposed for accuracy and numerical efficiency; see e.g. [8, 11, 13, 6]. Note that the kernel can be expressed in the form , where

 K0(r,r′):=ηsin(κ|r−r′|)4π|r−r′|+i(κcos(κ|r−r′|)−sin(κ|r−r′|)|r−r′|)(r−r′)⋅ν(r′)4π|r−r′|2 (2.4)

and

 K1(r,r′):=−iηcos(κ|r−r′|)4π|r−r′|−κ|r−r′|sin(κ|r−r′|)+cos(κ|r−r′|)4π|r−r′|(r−r′)⋅ν(r′)|r−r′|2. (2.5)
###### Remark 2.1

It is easy to check that . The kernel , on the other hand, while not , can be integrated with super-algebraic accuracy by means of a combination of a polar change of variables and the trapezoidal rule; see [6] and equation (2.23).

In this paper we focus on the Brakhage-Werner formulation presented above; a similar analysis can be used to treat the closely related Burton-Miller formulation [9].

### 2.2 Geometry

The numerical method studied in this work relies heavily on the use of a system of local charts for description of the surface . We will thus use a set of a number of open overlapping coordinate patches that cover ,

 S=J⋃j=1Sj, (2.6)

each one of which is the image of a parametrization

 rj:Dj⟶Sju––:=(u1,u2)⟼rj(u1,u2),

where . We assume that can be extended to a bijective diffeomorphism between and so that, in particular, the Jacobians

 aj(u––):=∣∣∣∂rj(u1,u2)∂u1×∂rj(u1,u2)∂u2∣∣∣>0,j=1,…,J (2.7)

() are functions of .

The method requires explicit use of a smooth partition of unity on , subordinated to the covering (2.6), that is

 ωj≥0,suppωj⊂Sj% andJ∑j=1ωj≡1.

The hypotheses on availability of local charts and a partition of unity is not restrictive in practice: such parameterizations can be constructed for smooth arbitrary geometries (see e.g. [5]). We will also assume that the boundary of is the finite union of Lipschitz arcs, a restriction which, again, is easy to accommodate [5].

For any and we let

 B(r,δ):={r′∈S:|r′−r|<δ}

and, selecting parameters that will otherwise be fixed throughout this paper, we define

 Sjδ:=¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯⋃r∈suppωjB(r,2ϵ1δ)=⋃r∈suppωj¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯B(r,2ϵ1δ). (2.8)

Clearly there exists such that for all

 ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯B(r,ϵ1δ0)∩Sj2δ0=∅∀r∈S∖Sj (2.9)

—as it can be checked by considering a pair of points that realize the distance between the boundaries of and . In particular, this implies that . The final element in our geometric constructions is a function such that

Given we now define the functions

 ηδ(r,r′):=υ(r,δ0|r′−r|δ). (2.10)

Clearly

 suppηδ(r,⋅)⊂¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯B(r,ϵ1δ),∀r∈S. (2.11)

In view of the definition of the sets and (2.9) we also have

 suppηδ(r,⋅)∩Sj3δ/2=∅∀r∈S∖Sj2δ. (2.12)

Given and considering the decomposition in term of the kernel functions given in equations (2.4) and (2.5), we define the regular part of the kernel of the integral operator (2.3) as

 Kreg,δ(r,r′):=K0(r,r′)+(1−ηδ(r,r′))K1(r,r′). (2.13)

Clearly . The remainder is the singular part of the kernel,

 Ksing,δ:=K−Kreg,δ=ηδK1, (2.14)

which, like the kernel , can be integrated accurately by means of a polar change of variables; see Remark 2.1. The parameter , which controls the support of the kernel , plays an essential role in both, the performance of the algorithm and its theoretical analysis; see Remark 2.4 for details.

We next introduce

 Kij(u––,v–):=K(ri(u––),rj(v–))aj(v–),

(see equation (2.7)) with corresponding definitions of and (cf. equations (2.13) and (2.14)). Noting that , and , are defined in , we extend these functions by zero (possibly thereby introducing discontinuities on the boundary of ) to the full product of squares . Clearly,

 ∫SK(ri(u––),r′)ξ(r′)dr′=J∑j=1∫DjKij(u––,v–)ωj(rj(v–))ξ(rj(–v))dv–∀u––∈Di.

Therefore, if is the exact solution of (2.2), then (), is a solution of the system

 12ψi(u––)+J∑j=1∫I2Kijreg,δ(u––,v–)ωj(rj(v–))ψj(v–)dv– +J∑j=1∫I2Kijsing,δ(u––,v–)ωj(rj(v–))ψj(v–)dv– = −Uinc(ri(u––)) (2.15) ∀u––∈Di,i=1,…,J.
###### Remark 2.2

Note that the functions that appear in the integrals over in equation (2.15) are only defined in . However, they are multiplied by the cutoff function which vanishes outside and thus provides a natural extension for the product throughout .

The solution of this system is used in Section 3 to reconstruct the solution of the original problem (2.2).

We can clearly distinguish two different types of integral operators in (2.15), namely, integral operators with smooth and singular kernels. The discretizations of these operators are produced, accordingly, by means of two different strategies—as discussed in the following section.

### 2.3 Discretization of the integral operators in equation (2.15)

For each patch , we select a positive integer , we let , and we introduce the grid points

 x––jm––:=hjm––=hj(m1,m2),0≤m1,m2≤Nj−1.

We assume that these grids are quasiuniform: letting

 N:=minj=1,…,JNj,h:=1/N=maxj=1,…,Jhj,

we assume there exists such that

 h

This is not a restrictive assumption in view of the assumed smoothness of and , and, therefore, of the solution .

The algorithm under consideration is a Nyström method which produces point-wise values of the unknown at the discretization points . The quadrature rules that ultimately define the method are described in what follows; for notational simplicity our approximate quadrature formulae use the set

 ZNj:={m––:=(m1,m2)∈Z2:0≤m1,m2≤Nj−1} (2.17)

of two-dimensional summation indices .

The approximate integration method used by the algorithm to treat the regular portion of the kernel is, simply, based on the trapezoidal rule,

 ∫I2Kijreg,δ(⋅,v–)ξ(v–)dv–≈h2j∑m––∈ZNjKijreg,δ(⋅,x––jm––)ξ(x––jm––), (2.18)

with the convention that for . Notice that for a function compactly supported in

 Ωj:=(rj)−1(suppωj)⊂Dj, (2.19)

(such as with ) the rule embodied in equation (2.18) gives rise to super-algebraic convergence.

To approximate integrals of the form

 ∫I2Kijsing,δ(u––,v–)ξ(v–)dv–,ξ∈C∞(I2),suppξ⊆Ωj, (2.20)

that include the singular part of the kernel, the algorithm uses polar coordinates around the singularity. To properly account for contributions arising from various regions delineated by local charts and overlaps we let

 r–ji:=(rj)−1∘ri:(ri)−1(Si∩Sj)⊂Di→DjSijδ:=Si∩Sjδ⊂Si∩Sj,Ωijδ:=(ri)−1(Si∩Sjδ)=(ri)−1(Sijδ); (2.21)

see Figure 2 and equation (2.8). A close inspection of the definition of the sets shows that if . Since, additionally,

 suppKijsing,δ(u––,⋅)⊂suppηδ(ri(u––),rj(⋅))⊂(rj)−1(¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯B(ri(u––),ϵ1δ)∩Sj),

we conclude that the integral (2.20) with , vanishes for all such that , i.e., for . The algorithm therefore only produces approximations of (2.20) for . To do this consider the relations

 ∫R2Kijsing,δ(u––,v–)ξ(v–)dv– (2.22) =∫R2Kijsing,δ(u––,r–ji(u––)+w––)ξ(r–ji(u––)+w––)dw–– =12∫2π0(∫∞−∞|ρ|Kijsing,δ(u––,r–ji(u––)+ρe–(θ))ξ(r–ji(u––)+ρe–(θ))dρ)dθ

where and note that the function

 (ρ,θ)↦|ρ|Kijsing,δ(u––i,r–ji(u––)+ρe–(θ)) (2.23)

is smooth (as shown in Section 4, cf. [6, 7]), -periodic in and compactly supported as a function of the variable . In what follows we temporarily let , so that the dependences on the patch indices and parametric coordinate are not displayed.

The integral in the angular variable is approximated using a trapezoidal rule on a uniform partition of in subintervals of length . Therefore, we consider the angles

 θjp:=pkj,p=0,…,Θj−1

and approximate (2.22) by

 kj2Θj−1∑p=0(∫∞−∞Kijsing,δ(u––,z–+ρe–(θjp))|ρ|ξ(z–+ρe–(θjp))dρ). (2.24)

Recall that a uniform Cartesian grid in with mesh length has been introduced in the approximation of the regular part (2.18). We will now use points of this grid to approximate the radial integrals in (2.24).

If , that is, modulo , , we look for the intersections

 {x––=z–+ρe(θjp):ρ∈R}∩{x––=(qhj,ξ):ξ∈R},q=0,…,Nj−1,

i.e., the intersections of the double rays stemming from with angle with the vertical lines of the uniform grid: the points of intersections are given by

 (qhj,z2+(qhj−z1)tanθjp),q=0,1,…,Nj−1; (2.25)

clearly, the distance between two consecutive points is Alternatively we can express these points in the form

 z–+ρhjq(z–,θjp)e–(θjp),whereρhq(z–,θ):=qh−z1cosθ.

For these values of , the corresponding integral in (2.24) is approximated by

 hjc(θjp)Nj−1∑q=0Kijsing,δ(u––,z–+ρhjq(z–,θjp)e–(θjp))|ρhjq(z–,θjp)|ξ(z–+ρhjq(z–,θjp)e–(θjp)), (2.26)

where . For the remaining angles, intersections are computed with the corresponding horizontal lines:

 {x––=z–+ρe–(θjp):ρ∈R}∩{x––=(ξ,qhj):ξ∈R},q=0,…,Nj−1. (2.27)

The quadrature points are deliberately selected to lie on lines with one of the coordinates equal to . As discussed in the introduction, such a selection was introduced in [6] to enable fast interpolation of the function on the radial quadrature points (that is, to produce a fast algorithm for evaluation of the operator in (2.33)).

Using the angle-dependent weight

 c(θ):=min{1|cosθ|,1|sinθ|}, (2.28)

as well as the nodal points radii

 ρhq(z–,θ):=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩qh−z1|cosθ|,if |cosθ|≥√2/2,qh−z2|sinθ|,otherwise,

expression (2.26) provides approximation of all the integrals in (2.24).

In sum, we define our discrete singular operator (which depends on and , or, equivalently, on and ) by

 (2.29)

where

In what follows we use the scaling

 Θj≈N1+αj,α>0, (2.30)

i.e., we assume that there exist such that for all , . We point out that is discontinuous at but both side limits exist. It is immaterial which value (the right or left limit or even the average of both) is set as value to at these points.

### 2.4 The numerical method

We are now in a position to lay down a fully discrete version of equations (2.15). The unknowns are taken to be approximate values

 (2.31)

of the true function values . In what follows we denote by the array whose entry is for all indices . We consider the set of indices

 Z∗Nj:={m––:=(m1,m2)∈Z2:−Nj/2≤m1,m2

(note that the cardinality of is ), the space of trigonometric polynomials

 TNj:=span{exp(2πim––⋅u––):m––∈Z∗Nj},

the interpolation operator defined by the equations

 RNjφ(x––jℓ–)=⎧⎨⎩φjℓ–,ℓ–∈Ωjh,0,ℓ–∈ZNj∖Ωjh,

and the vectors with components . Using these notations, the discrete Nyström equations for the system (2.15) are given by

 (2.33) =−Uinc(ri(x––iℓ–))

for and , where the binary operator denotes componentwise product of arrays: for example, is the array whose -th entry is given by . Once these point values have been computed, the Nyström method provides the reconstruction formula

 ψih=−2J∑j=1(h2j∑m––∈ΩjhKijreg,δ(⋅,x––jm––)ωjm––φjm––+Lijδ,hRNj(ωj⊙φj))−2Uinc∘ri, (2.34)

(). With these definitions, we clearly have

 ψjh(x––jℓ–)=φjℓ–ℓ–∈Ωjh.

Finally, the parameter-space continuous functions can be assembled into a single continuous approximate solution defined on :

 ψh(r):=∑j∈I(r)ωj(r)ψjh((rj)−1(r)),where I(r):={j:r∈suppωj}. (2.35)

The main convergence result of this paper can now be stated; its proof is given in Sections 3 through 5. The regularity estimates given in the statement of the theorem are expressed in terms of standard Sobolev norms on the surface (see for instance [1, 21, 23] or any standard text on Sobolev spaces) .

###### Remark 2.3

The notation means that

 chβ≤δ≤Chβ,

for some constants and , independent of and . Hereafter, the parameter will be taken to be fixed. Dependence of constants on this parameter will not be shown explicitly.

###### Theorem 2.1

Let be the solution of (2.2), and assume with . Then for sufficiently large values of the Nyström equations (2.33) admit a unique solution. Letting denote the reconstructed function on as defined by (2.34) and (2.35), we have the error estimate

 ∥ψh−ψ∥L2(S)≤Ct(|logh|32ht(1−β)+ht−1)∥ψ∥Ht(S), (2.36)

for all . Finally, for all and

 ∥ψh−ψ∥L∞(S)≤Ct,ε(ht(1−β)−(1+ε)β+ht−1)∥ψ∥Ht(S). (2.37)

Theorem 2.1 tells us that for a smooth surface and a right-hand side (from which it follows that ) the Nyström algorithm under consideration converges super-algebraically fast.

###### Remark 2.4

As mentioned in Section 2.2, the parameter plays central roles in both the theory and the actual performance of the the algorithm under consideration. With regards to the latter we briefly mention here that use of the floating cut-off function (2.10), whose support is controlled by the parameter , helps restrict the use of the costly polar integration scheme to a small region around the singular point (thus reducing the overall computing time required by the algorithm), and, further, it enables acceleration via a sparse, parallel-face FFT-based equivalent-source technique; see [6] for details. (In particular we note that the value is used in [6, 7] for optimal speed of the equivalent-source accelerated Nyström method.) The parameter also has a significant impact on our theoretical treatment. Indeed, one of the most delicate points in our stability analysis concerns the convergence in norm of a discrete singular operator to a corresponding continuous singular operator. One of the terms in our estimate of the norm of the difference of these operators (equation (5.18) in Proposition 5.7) is bounded by . By taking we ensure that this terms also tends to zero, from which the desired convergence in norm results.

## 3 Biperiodic framework

### 3.1 Continuous equations

The analysis of the method will be carried out by recasting the system (2.15) as a system of periodic integral equations with all unknown functions defined on the unit square . To introduce our periodic formulation we utilize a second family of cut-off functions, , that depends on and satisfies the following assumptions:

 supp˜ωjδ⊂Ωjj3δ/2=(rj)−1(Sj3δ/2), (3.1a) ˜ωjδ ≡ 1 in Ωjjδ=(rj)−1(Sjδ)⊃Ωj, (3.1b) ∥∂α––˜ωjδ∥∞,I2≤Cα––δ−|α––|∀α––≥0–, (3.1c)

where for a given non-negative bi-index ,

 ∂α––ξ(u––)=∂α1u1∂α2u2ξ(u1,u2)

and is the norm. We will use the characteristic function

 ˜ωj0(u––):={1% if u––∈Ωj,0otherwise, (3.2)

which can be viewed as the limit of as .

Using these functions, we define the following periodic integral operators

 Aijreg,δξ := (ωi∘ri)∫I2Kijreg,δ(⋅,u––)˜ωjδ(u––)ξ(u––)du––, (3.3) Aijsing,δξ := (ωi∘ri)∫I2Kijsing,δ(⋅,u––)˜ωjδ(u––)ξ(u––)du––, (3.4) Aijδξ := (Aijreg,δ+Aijsing,δ)ξ, (3.5)

as well as the right-hand sides , properly extended by zero to the full unit square . If , then and

 supp(Kijsing,δ(u––,⋅)˜ω