Corrected trapezoidal rules for boundary integral equations in three dimensions

The manuscript describes a quadrature rule that is designed for the high order discretization of boundary integral equations (BIEs) using the Nyström method. The technique is designed for surfaces that can naturally be parameterized using a uniform grid on a rectangle, such as deformed tori, or channels with periodic boundary conditions. When a BIE on such a geometry is discretized using the Nyström method based on the Trapezoidal quadrature rule, the resulting scheme tends to converge only slowly, due to the singularity in the kernel function. The key finding of the manuscript is that the convergence order can be greatly improved by modifying only a very small number of elements in the coefficient matrix. Specifically, it is demonstrated that by correcting only the diagonal entries in the coefficient matrix, \(O(h^{3})\) convergence can be attained for the single and double layer potentials associated with both the Laplace and the Helmholtz kernels. A nine-point correction stencil leads to an \(O(h^5)\) scheme. The method proposed can be viewed as a generalization of the quadrature rule of Duan and Rokhlin, which was designed for the 2D Lippmann–Schwinger equation in the plane. The techniques proposed are supported by a rigorous error analysis that relies on Wigner-type limits involving the Epstein zeta function and its parametric derivatives.

The authors would like to thank Alex Barnett, Mike O’Neil, Vladimir Rokhlin, and Johan Helsing for sharing valuable perspectives and insights. The work reported was supported by the Office of Naval Research (Grant N00014-18-1-2354) and by the National Science Foundation (Grants DMS-1620472, DMS-2012606, and DMS-1952735).

Authors and Affiliations

1. Oden Institute for Computational Engineering and Sciences, University of Texas at Austin, Austin, TX, 78712, USA

   Bowei Wu

2. Department of Mathematics, University of Texas at Austin, Austin, TX, 78712, USA

   Per-Gunnar Martinsson

Corresponding author

Correspondence to Bowei Wu.

Wigner-type limit formulae for the Laplace SLP

To complete the expressions (71), we give formulae for \({\mathcal {C}}_1,{\mathcal {C}}_2\) and \({\mathcal {C}}_{01}\). Firstly, \({\mathcal {C}}_1\) and \({\mathcal {C}}_2\) are the Wigner-type limits associated with \({\mathcal {W}}^{(2)}\) defined by (69), which we restate here

where the numerator is a homogeneous polynomial:

such that

where all the involved derivatives are evaluated at \({\mathbf {0}}.\) Thus we have the following Wigner-type limits

where

is a vector-valued second-derivative operator, and where the constants in the vectors \(L^{(A)}_{3,1}\) and \(L^{(A)}_{3,2}\) correspond to the terms \(u\,{\mathcal {W}}^{(2)}({\mathbf {u}})\) and \(v\,{\mathcal {W}}^{(2)}({\mathbf {u}})\), respectively, therefore

where “\(*\)” represents convolution (given that these are coefficients of the product of polynomials). Next, \({\mathcal {C}}_{01}\) is the Wigner-type limit associated with the function \({\mathcal {W}}^{(3)}\), given by

then, expanding the denominator as above, one derives that

where

is a vector-valued third-derivative operator, and where the vectors \(L^{(A)}_4\) and \(L^{(A)}_6\) correspond to the expansions of \((2{\mathbf {d}}^1\cdot {\mathbf {d}}^3+{\mathbf {d}}^2\cdot {\mathbf {d}}^2)\) and \((2{\mathbf {d}}^1\cdot {\mathbf {d}}^2)^2\), respectively, such that

where

\(O(h^5)\) corrected trapezoidal rule for the Laplace DLP

Consider the DLP kernel \({\mathcal {K}}({\mathbf {u}})\) given by

then the \(O(h^5)\) quadrature for the DLP is given by

(89)

where \(\sigma \) is a smooth density function and \(U_2\) is the 9-point stencil as defined in (72). The correction weights satisfy the system (80) but with the quantities \(D_0\)–\(D_5\) replaced by (92) below, which we will derive next.

To derive the formulae for the correction weights, we follow an analysis similar to Sect. 5.1. First the DLP kernel \({\mathcal {K}}({\mathbf {u}})\) is expanded using the binomial series

thus

Then the numerators \(O((r^2-Q)^m)\) are further expanded as in (44); the term \(-{\mathbf {r}}\cdot ({\mathbf {r}}_u\times {\mathbf {r}}_v)\) can be expanded as

where

where the involved constants are defined as follows (all involved derivatives are evaluated at \({\mathbf {0}}\)),

Similar to Table 1, we summarize in Table 3 the singularity information associated with the DLP kernel expansion (91).

Similar to (71), we define the following quantities based on Table 3 which are useful for constructing the \(O(h^5)\) quadrature (89) for the DLP:

where

The involved vector-valued differential operators, besides (86, 88), are defined as

and the involved constant vectors are

and where \(L^{(A)}_3,L^{(A)}_4\) and \(L^{(A)}_6\) are as appeared in “Appendix A” in the formulae for the SLP.

Substituting the above constants \(D_0\)–\(D_5\) into the system (80) results in the correction weights for the DLP in (89).

\(O(h^5)\) corrected trapezoidal rule for the normal gradient of Laplace SLP

We will omit the derivations in this section and just present the quadrature formula for the normal gradient of the Laplace SLP on a surface.

The kernel of the SLP normal gradient, centered at \({\mathbf {0}}\), is

then the associated \(O(h^5)\) quadrature is given by

(94)

where \(P({\mathbf {u}})=\sigma ({\mathbf {u}})\,|{\mathbf {r}}_u({\mathbf {u}})\times {\mathbf {r}}_v({\mathbf {u}})|\) is smooth and \(\sigma \) is a smooth density function, and where \(U_2\) is the 9-point stencil as defined in (72). The correction weights satisfy the system (80) but with the quantities \(D_0\)–\(D_5\) replaced by the following:

where

such that

and where \(L_A^{(3)}, L_A^{(4)}, L_A^{(6)}\) are as defined in Sect. A, and \(\square ^{(k)}, k=1,2,3,4\), are as defined in (86, 88, 93).

\(O(h^5)\) corrected trapezoidal rules for the Helmholtz potentials

The Helmholtz SLP, DLP, and the normal gradient of SLP (denoted SLPn) are related to the Laplace potentials by the following expansions:

where the first term in each expansion is the corresponding Laplace kernel. Therefore based on the \(O(h^5)\) quadratures for the Laplace potentials, one only needs to include the errors up to \(\varTheta (h^3)\) from the additional terms to construct the corresponding quadratures for the Helmholtz kernels. Specifically, these relevant additional terms are

where the term \(i\kappa \) is regular, all other terms associate to Wigner-type limits that are \(O(h^3)\), hence only the leading errors need correction. By analyses similar to Sect. 5, quadratures for the Helmholtz potentials are constructed as follows:

• The \(O(h^5)\) quadrature for the Helmholtz SLP uses the same correction weights as the Laplace SLP constructed by (71) and (81), except that \(D_0\) is modified as

  $$\begin{aligned} D_0\leftarrow D_0+\Big (i(\kappa h) + \frac{(\kappa h)^2}{2}Z(-1)\Big )\,h. \end{aligned}$$

• The \(O(h^5)\) quadrature for the Helmholtz DLP uses the same correction weights as the Laplace DLP constructed by (92) and (81), except that \(D_0\) is modified as

  $$\begin{aligned} D_0\leftarrow D_0-\frac{(\kappa h)^2}{2}(L^{(B)}_2\cdot \square ^{(1)})Z(-1)\,h. \end{aligned}$$

• The \(O(h^5)\) quadrature for the normal gradient of the Helmholtz SLP uses the same correction weights as the corresponding Laplace potential constructed by (95) and (81), except that \(D_0\) is modified as

  $$\begin{aligned} D_0\leftarrow D_0-\frac{(\kappa h)^2}{2}(L^{(C)}_2\cdot \square ^{(1)})Z(-1)\,h. \end{aligned}$$

Formulae for Z(s) and its parametric derivatives

To simplify notations, first define \(s_1 :=\frac{s}{2}, s_2 := 1-s_1\) and the customed incomplete gamma function

then

where

1.1 Scalar derivatives of Z(s)

We denote the scalar parametric k-th derivative of Z(s) as

To derive the derivative formulae based on (96), first note the following properties of g(s, x)

These recurrence relations will allow one to evaluate all the derivatives of Z(s) with only one gamma function evaluation at each point, greatly reducing the cost.

The formulae \(\square ^kZ(s)\) for \(k\le 4\) are then given by

To further reduce costs, recall for example the \(O(h^3)\) errors used \(Z(-1)\) while the O(h) errors used Z(1), thus one can compute \(Z(s+2)\) along the way of computing Z(s), with

In all the formulae for the Epstein zeta function, the various involved quantities are defined as follows.

Firstly, the involved quadratic forms and their derivatives are

where

and the other involved constants such as \(H_D\) and their parametric derivatives are given by

Note that if L, M, N are the second fundamental form coefficients, then \(H_D\) and \(K_D\) become the mean and Gaussian curvatures.

1.2 Vector-valued derivatives \(\square ^{(k)}Z(s)\) via the scalar derivatives \(\square ^{k}Z(s)\)

We show that the vector-valued k-th derivative \(\square ^{(k)}Z(s)\), as defined in (86, 88, 93), can be expressed as a combination of scalar derivatives \(\square ^{k}Z(s)\), for \(k=2,3,4\).

Firstly, the 5-component vector \(\square ^{(2)}Z(s)\) can be computed via 5 scalar second derivatives:

so

Secondly, the 7-component vector \(\square ^{(3)}Z(s)\) can be computed via 7 scalar third derivatives:

so

Lastly, to compute \(\square ^{(4)}\), we first define

and denote the 9 components of \(\square ^{(4)}\) as \(\square ^{(4)}_i\), \(i=1,\dots ,9\), then

where

Proof of the convergence rate of Wigner-type limits

Remark 9

The proofs of Theorem 8 and Lemma 2 in this section are inspired by [30, Theorem 3.1, Lemma 3.3], which are the one-dimensional analogs. As a result, our proofs complete the convergence analysis in [30] for their two-dimensional quadrature.

Throughout this section, \(\eta ({\mathbf {u}})\in C_c^{\beta }([-a,a]^2)\) is a \(\beta \)-times continuously differentiable function that is compactly supported in \([-a,a]^2\), and satisfies

for all non-negative integers i and j such that \(0<i+j\le 2K.\)

We restate the definition (53) as follows:

with

Theorem 7

Using the notations in Theorem 4, then for any integer \(K \ge 1\) and \(m, d \ge 0\) such that

the estimate

holds provided that \(\eta ({\mathbf {u}})\) satisfies (101) and is 4K-times contiuously differentiable.

Proof

Let \(s = 2(m-d)+1\) then \(3-2K \le s \le 1\), then using Theorem 8 (to be proved below) with \(\beta = 4K\) yields

Taking appropriate parametric derivatives on both sides according to Theorem 4 gives

where the order of the error term remains unchanged because of the following reason.

Parametric differentiation in the form of Theorem 4 only acts on the component \(Q({\mathbf {u}})^{-s/2}\) and is in such a way that does not change its homogeneous order in \({\mathbf {u}}\), which is \(\varTheta (|{\mathbf {u}}|^{-s})\). Therefore the analysis in Theorem 8 holds under parametric differentiation. \(\square \)

Theorem 8

If \(s\in {\mathbb {R}}, s\le 1<2K\) and \(\eta \) satisfies (101) and is at least \(\beta \) times continuously differentiable such that \(\beta + s \ge 2K+3\), \(K\ge 1\), then

Proof

To analyze the convergence property of

we define a smooth cutoff function \(\phi ({\mathbf {u}})\), such that \(\phi ({\mathbf {u}}) \equiv \phi (-{\mathbf {u}})\) and

and denote \(f({\mathbf {u}},h) := \eta ({\mathbf {u}}h)\,Q({\mathbf {u}})^{-\frac{s}{2}}\), then a partition of unity using \(\phi ({\mathbf {u}})\) gives

where the fact that \(f({\mathbf {0}},h)\phi ({\mathbf {0}}) = 0\) and \(1-\phi ({\mathbf {i}})=0\) if \({\mathbf {i}}\ne {\mathbf {0}}\) are used. We then estimate the quantities \(I_1\) and \(I_2\).

• We first bound \(I_2\). Use \(|\eta ({\mathbf {u}}h)-1|\le C' |{\mathbf {u}}h|^{2K}\) for some constant \(C'\), then

  $$\begin{aligned} \begin{aligned}&\left| \int _{|{\mathbf {u}}|\le 1}(\eta ({\mathbf {u}}h)-1)\,Q({\mathbf {u}})^{-\frac{s}{2}}(1-\phi ({\mathbf {u}}))\,\mathrm {d}{\mathbf {u}}\right| \\&\quad \le C' h^{2K} \int _{|{\mathbf {u}}|\le 1}|{\mathbf {u}}|^{2K}|Q({\mathbf {u}})|^{-\frac{s}{2}} |1-\phi ({\mathbf {u}})|\mathrm {d}{\mathbf {u}}\le Ch^{2K}, \end{aligned} \end{aligned}$$

  where we used the fact \(|{\mathbf {u}}|^{2K}|Q({\mathbf {u}})|^{-s/2} = O(|{\mathbf {u}}|^{2K-s})\) is integrable since \(2K-s>0\). So, if we define

  $$\begin{aligned} Z_2(s) = \int _{|{\mathbf {u}}|\le 1}Q({\mathbf {u}})^{-\frac{s}{2}}(1-\phi ({\mathbf {u}}))\,\mathrm {d}{\mathbf {u}}, \end{aligned}$$

  which is independent of h and integrable for \(s\le 1\), then

  $$\begin{aligned} \begin{aligned} |I_2&- Z_2(s)| = \left| \int _{|{\mathbf {u}}|\le 1}(\eta ({\mathbf {u}}h)-1)\,Q({\mathbf {u}})^{-\frac{s}{2}}(1-\phi ({\mathbf {u}}))\,\mathrm {d}{\mathbf {u}}\right| \le Ch^{2K}. \end{aligned} \end{aligned}$$

• For the summation and integral in \(I_1\), since the function \(f({\mathbf {u}},h)\) is compactly supported, the Poisson summation formula (cf. [36, Chapter 4]) gives

  $$\begin{aligned} \sum _{|{\mathbf {i}}|<\infty }f({\mathbf {i}},h)\phi ({\mathbf {i}}) - \int _{|{\mathbf {u}}|<\infty }f({\mathbf {u}},h)\phi ({\mathbf {u}})\,\mathrm {d}{\mathbf {u}}= \sum _{|{\mathbf {k}}|\ne 0} {\hat{f}}_\phi ({\mathbf {k}},h), \end{aligned}$$

  (104)

  where \({\mathbf {k}}= (k_1,k_2)\in {\mathbb {Z}}^2\), and where

  $$\begin{aligned} {\hat{f}}_\phi ({\mathbf {k}},h) := \int _{|{\mathbf {u}}|<\infty }f({\mathbf {u}},h)\phi ({\mathbf {u}})e^{-2\pi i {\mathbf {k}}\cdot {\mathbf {u}}}\,\mathrm {d}{\mathbf {u}}. \end{aligned}$$

  Using Lemma 2 which we will prove next, there is a smooth function \(Z_1(s)\) independent of h, such that

  $$\begin{aligned} \left| \sum _{|{\mathbf {k}}|\ne 0} {\hat{f}}_\phi ({\mathbf {k}},h)-Z_1(s)\right| = O(h^{2K+1}). \end{aligned}$$

Combining all of above, we have

On the other hand, by Theorem 4,

is analytic. Since both \(Z_1\) and \(Z_2\) are smooth and neither of them depend on h, one must have

\(\square \)

Lemma 2

Assume that \(s\le 1\) and \(\eta \) satisfies (101) and is at least \(\beta \) times continuously differentiable such that \(\beta + s \ge 2K+3\), \(K\ge 1\). Then for \({\hat{f}}_\phi ({\mathbf {k}},h)\) as defined in the Poisson summation formula (104), there is a smooth function \(Z_1(s)\) that does not depend on h, such that

Before proving the lemma, we define the following simplified notations to be used throughout the proof. Denote the function

and for \({\mathbf {k}}=(k_1,k_2)\) and \({\mathbf {u}}=(u,v)\) define the following

We also use the shorthand notation

to denote the derivatives of any smooth function g when the associated \({\mathbf {k}}\) is clear from context.

Proof

Since \(\eta ({\mathbf {u}}h)\phi ({\mathbf {u}})S({\mathbf {u}},s)\) is compactly supported in \(\frac{1}{2}\le {\mathbf {u}}\le \frac{a}{h}\) and at least \(\beta \)-times continuously differentiable, then after integration-by-parts \(\beta \) times we have

Define

then, since \(\phi ^{(j)}({\mathbf {u}})\) is compactly supported for \(j\ge 1\) and \(|S^{(\beta )}({\mathbf {u}},s)| = O(|{\mathbf {u}}|^{-s-\beta })\) integrable for \(\beta +s\ge 2K+3\ge 5\), we have

for some constant C, where \(d_{\beta ,j}\) are binomial coefficients. So \(W({\mathbf {k}})\) is independent of h and uniformly bounded with the appropriate \({\mathbf {k}}\) and \(D^j\) correspondence as defined before the proof.

Next, define \(c({\mathbf {u}}):=\phi ({\mathbf {u}})\eta ({\mathbf {u}}h)\) and consider

Using the fact that \(\phi ^{(j)}({\mathbf {u}}) \equiv 0\) for \(j > 0\) and \(|{\mathbf {u}}| > 1\), we have

therefore, using the fact that \(\phi ({\mathbf {u}})\equiv 0\) for \(|{\mathbf {u}}|\le \frac{1}{2}\), \(\phi ({\mathbf {u}})\equiv 1\) for \(|{\mathbf {u}}|\ge 1\) and \(\eta ({\mathbf {u}}h)\equiv 0\) for \(|{\mathbf {u}}|\ge \frac{a}{L}\), the integrals in (105) can be rewritten as follows

where \(\delta _{i,j}\) is the Kronecker delta. Note that \(\eta ^{(\ell )}({\mathbf {0}}) = 0\) for \(\ell \le 2K\), so by Taylor’s theorem, and using the fact \(\eta ({\mathbf {0}})=1\), we have

for some constant \(C'\). Then we can estimate each of \(I_1^{(\ell )},I_2,I_3\) as follows.

• For the first sum of integrals

  $$\begin{aligned} |I_1^{(\ell )}| \le h^\ell \, (C' \,h^{2K+1-\ell })\int _{\frac{1}{2}\le |{\mathbf {u}}|\le 1}|{\mathbf {u}}|^{2K+1-\ell }\,|S^{(\beta -j)}({\mathbf {u}},s)|\mathrm {d}{\mathbf {u}}\le C\,h^{2K+1}, \end{aligned}$$

  so

  $$\begin{aligned} \left| \sum _{\ell =0}^jd_{j,\ell }I_1^{(\ell )}\right| \le \sum _{\ell =0}^jd_{j,\ell }|I_1^{(\ell )}|\le C\,h^{2K+1}. \end{aligned}$$

• For the second integral

  $$\begin{aligned} \begin{aligned} |I_2|&\le h^j\,(C'\,h^{2K+1-j})\int _{1<|{\mathbf {u}}|\le \frac{a}{h}}|{\mathbf {u}}|^{(2K+1-j)+(-s-\beta +j)}\mathrm {d}{\mathbf {u}}\\&= C'\,h^{2K+1}\int _{1<|{\mathbf {u}}|\le \frac{a}{h}}|{\mathbf {u}}|^{2K+1-s-\beta }\mathrm {d}{\mathbf {u}}\\&\le C (h^{2K+1} + h^{s+\beta -2}). \end{aligned} \end{aligned}$$

• For the third integral

  $$\begin{aligned} |I_3| \le \int _{\frac{a}{h}<|{\mathbf {u}}|<\infty }|{\mathbf {u}}|^{-s-\beta }\mathrm {d}{\mathbf {u}}\le C\,h^{s+\beta -2}. \end{aligned}$$

Therefore, substituting all above estimates back into (105) yields

By assumption, \(\beta +s \ge 2K+3\) and \(s\le 1\), so substitute \({\mathbf {k}}\mapsto 2\pi {\mathbf {k}}\) in the above formula yields

Therefore, if we define

which converges since \(\beta \ge 2K+3-s \ge 4\) given \(s\le 1\) and \(K\ge 1\), then

\(\square \)

Reprints and permissions