High-order corrected trapezoidal rules for a class of singular integrals

We present a family of high-order trapezoidal rule-based quadratures for a class of singular integrals, where the integrand has a point singularity. The singular part of the integrand is expanded in a Taylor series involving terms of increasing smoothness. The quadratures are based on the trapezoidal rule, with the quadrature weights for Cartesian nodes close to the singularity judiciously corrected based on the expansion. High-order accuracy can be achieved by utilizing a sufficient number of correction nodes around the singularity to approximate the terms in the series expansion. The derived quadratures are applied to the implicit boundary integral formulation of surface integrals involving the Laplace layer kernels.


Introduction
The trapezoidal rule is a simple and robust algorithm for approximating integrals.In general, it has second order accuracy, but when applied to compactly supported smooth functions, the accuracy is much higher.However, as the integrand becomes less smooth, the accuracy deteriorates, which makes the method unsuitable for singular integrals such as those found in boundary integral equations.This paper develops a systematic approach to derive high order corrected trapezoidal rules for integrals involving a class of integrands that are singular at one point.
Let f : R n \ {0} → R be a compactly supported function with an integrable singularity at 0. A crude way to approximate its integral is with the "punctured" trapezoidal rule T 0 h , where h denotes the discretization parameter.It equals the standard trapezoidal rule, but sets f ≡ 0 in a in a small h-dependent region N h surrounding the singular point This gives a low order accurate approximation.With R h [f ] denoting the error in the quadrature rule, we write (1.1) Figure 1: Singular behavior of a layer kernel.Two plots related to the double-layer kernel ∂G ∂ny (x, ȳ) = ∂G ∂ny (v(0, 0), v(θ, φ)) from boundary integral formulations, where v(θ, φ) is a surface parametrization centered around x = v(0, 0).The kernel is in the form (1.2) using (θ, φ) near its singular point.On the left, we plot the kernel ∂G ∂ny (v(0, 0), v(θ, φ)) on a uniform grid around (0, 0).On the right we plot (|w|, w/|w|), w := (θ, φ), obtained by multiplying the same kernel by |w|.We can clearly see that (0, u) is not constant in u.
One direction to improve the accuracy is to add back the function values at the excluded points in N h with judiciously chosen weights, such that they well approximate R h [f ].This can be seen as a correction of the standard trapezoidal rule, locally around the singularity.The overall simplicity of the method is therefore maintained.The approach has been used for example in [6,12,21,20] for the trapezoidal rule and in [18,10] for other quadrature methods.
In this article, we consider two-dimensional singular integrands f = s v, where v ∈ C ∞ c (R 2 ) and s is of the following form: for some smooth function : R × S 1 → R. Nevertheless, (|x|, x/|x|) is not necessarily a smooth function of x in 0 if (0, u) is non-constant in u.Singular functions of the type (1.2) are found in many applications.For instance, if g : R 2 → R 3 has a simple zero at the origin, then s(x) = 1/|g(x)| is of this type, as proven in Lemma 3.4.They also characterize the singular behavior of the kernels found in the boundary integral equations for elliptic problems.In three dimensions, the kernel is a function of two spatial variables, x, ȳ, but the integral typically involves the product of the kernel and a smooth function over a smooth and compact surface.In this setup, the singularity in the integrand depends on 1/|x − ȳ| and also on "the angle of approach", which corresponds to the way ȳ approaches x along a two dimensional surface.See Figure 1 for an illustration of the singular behavior of boundary integrals layer kernels.Note that we limit ourselves to compactly supported integrands.These can also be seen as the restrictions of periodic functions which are smooth away from the point singularity.If the integrand is not zero at the boundary of the integration domain, additional boundary corrections must also be introduced; see discussions in [6,1].
In [4], we derived a second order accurate method.In this paper, we generalize our approach systematically to derive higher order methods.Our approach is to Taylor expand the function in its first argument, and recognize that the smoothness of the remainder term increases with order, and can eventually be integrated accurately with the standard trapezoidal rule.We therefore only need to derive corrections for the leading Taylor terms, which are all of the form |x| j φ(x/|x|), for some j and φ : S 1 → R. The details are presented in Section 2.
One of the main motivations for the proposed approach is to provide the Implicit Boundary Integral Methods (IBIMs, see [7]) with high order convergent quadratures.IBIMs are volumetric integral formulations of classical boundary integrals and do not rely on explicit parameterization of surfaces (the "boundary" in the boundary integrals).The IBIM approach gives a way to compute accurate surface integrals, integral equations and variational problems on surfaces for other non-parametric methods, including the level set methods, e.g.[14,17,13,2] , and the closest point methods, e.g.[16,11].
In Section 3, we apply our new high order corrected trapezoidal rules to the singular integrals derived from IBIMs.In that formulation, the integrand is singular along a line and for each fixed plane it has a point singularity.To compute the volumetric integrals from IBIMs, our quadrature rules for integration in two dimensions are therefore applied plane by plane; see Section 3.1.We show in Theorem 3.2 that the resulting singularity on each plane is of the type in (1.2), and derive explicit expressions for the required terms in the expansion.Some efforts are needed to extract the needed geometrical information of the surface.Specifically, intrinsic information about the surface (principal directions and curvatures, and third derivatives of its local representation) together with extrinsic information (signed distance function to the surface) are needed to apply the quadrature rule in addition to the information needed for the IBIM formulation.
Finally, numerical simulations for selected problems in two and three dimensions are presented in Section 4.

The corrected trapezoidal rules
The standard trapezoidal rule has a low order accuracy when applied to singular integrals.In this section we show how one can raise the order of accuracy for integrands that are singular at a point, by correcting the computations at a few grid points close to the singularity.This type of corrections have been applied successfully in a few settings earlier.
We begin by defining the trapezoidal rules that we will work with.Let f be an integrable, compactly supported, function on R n .We are interested in approximating the integral R n f (x)dx by summation of the values of f on the uniform grid hZ n .Since f is supported in a compact set, the standard trapezoidal rule becomes the following simple Riemann sum: In this case the order of accuracy of the approximation is only limited by the regularity of f .If f ∈ C p c (R n ), the error is at worst O(h p ). See e.g. the discussion and proofs in [15].In particular, the trapezoidal rule enjoys spectral accuracy Here C p c (R n ) denotes the space of compactly supported functions on R n whose partial derivatives up to order p are continuous (of all orders, if p = ∞).
If f is smooth in R n \ {x 0 }, singular at x 0 , and R n f (x)dx exists as a Cauchy principal value, it is natural to modify the trapezoidal rule by excluding the summation over some grid nodes close to x 0 .We define the punctured trapezoidal rule with respect to N h as where N h (x 0 ) defines a small neighborhood around x 0 , the region being "punctured" from R n .When x 0 lies on a grid node, one typically sets N h (x 0 ) = {x 0 }; i.e.only the singularity point is removed from the standard trapezoidal rule.If x 0 does not lie on a grid node, one option is to remove the grid node x h that is closest to it.In this case, N h (x 0 ) = {x h }.
In general, N h (x 0 ) may contain several grid nodes, although the number is typically finite and independent of h.We will write N h,m to indicate that the set contains m nodes.
Here we consider an integrand that is the product of a smooth factor v and a singular factor s, which takes the form (1.2) near the origin.The punctured trapezoidal rule converges for such singular functions, albeit with a lower rate.In the case N h (x 0 ) = {x h }, we have the following theorem.
where the constant C is independent of h, but depends on j, and v.
The proof is given in Appendix A.1, where without loss of generality we consider x 0 = 0 and N h (0) = {0}.
We now give a brief summary of the steps that we shall take in Sections 2.1-2.5 to correct the trapezoidal rule for the two-dimensional case n = 2 and j = −1 in Theorem 2.1.
In Section 2.1 we expand in its first argument to derive a series of the form for some functions s k and q s.Theorem 2.1 states that the error R h , as defined in (1.1), in applying In Section 2.2 we derive a weight ω for approximating the error

Multiplication of the weight by any smooth function
where x h is the grid node in hZ 2 closest to x 0 .In addition, the weight depends on s k but not on h and v.With this weight, we define Consequently, We then generalize this approach systematically in Sections 2.3 and 2.4.Eventually, we obtain quadratures formally for any k ≥ 0.Here p ≥ p is a constant and x h,i are grid points near x 0 .These will be described more carefully later.
Remark 2.2.The quadrature rule Q p h [s k ] depends on the value of k in the subscript of s k , in addition to the function s k , but for simplicity of notation we will not make this distinction.
Finally, in Section 2.5 we combine the quadratures Q p h for s k v to define a quadrature U p h of order p ≥ 2 for the function s v (recall that s is expanded into a sum of s k for k = 0, 1, • • • , q, and q s).The order p specifies how many expansions terms are needed (q = p − 2) and which quadratures derived from correcting the punctured trapezoidal rules are needed for each term (Q p−1−k h for s k , and

Expansion of the singular function
To integrate the function s 2) with high order accuracy, we use a divide et impera strategy.For any u ∈ S 1 , we expand (r, u) with respect to the first variable, and approach each of the expansion components separately: becomes and we write s formally as the series where If we use q terms in this expansion, we expect that these terms s k strip away the singularity in s at x = 0 so that what is left behind from the expansion, i.e. the remainder term can be approximated directly with the (unmodified) trapezoidal rule and achieve the order of accuracy desired without needing special quadrature.This property is expressed in the following lemma.
For any integer q ≥ 0, there exist σ : R × S 1 → R such that σ ∈ C ∞ ((−r 0 , r 0 ) × S 1 ) and q s(x) = |x| q σ(|x|, x/|x|). (2.5) The proof of this lemma can be found in Appendix A.2. From this result and the previous Theorem 2.1 we can express the following lemma.
Lemma 2.4.The term q s in (2.4) is integrated by the punctured trapezoidal rule (2.2) with order q + 2: Hence, to get high order, it is sufficient to derive corrected trapezoidal rules for the expansion terms We start from first order correction for s k in Section 2.2 and end with a general description for an arbitrarily high order method for s in Section 2.5.

First order correction
Our goal in this section is to derive a first order in h correction for the punctured trapezoidal rule applied to for s k of the form (2.6).We assume ).This correction will yield an error with its largest part proportional to h k+2 .

The singular point rests on a grid node
Without loss of generality, we assume that x 0 = 0 and lies on a grid node.For such cases, the set N h typically contains only the grid node where the singularity is.In our case, N h (0) = {0}.The smoothness of s k in (2.6) increases with k.Theorem 2.1 tells us that for a function of this kind (two-dimensional, j = k − 1) the error behaves as: Following [12], one can show that the error has the form where ω[s k ] is a constant independent of v and h.In [12] this is proven for s 0 (x) = 1/|x| (k = 0 and φ 0 ≡ 1).
Hence we define the first order correction Q 1 h to the punctured trapezoidal rule as This quadrature rule thus corrects the trapezoidal rule in one node, the origin.It will then have an error of size O(h k+2 ).
Remark 2.5.In the special case when φ k ≡ 1, due to symmetry with respect to the grid node at 0, the O(h k+2 ) terms cancel out, and Q 1 h achieves an accuracy of O(h k+3 ).
To find the weight ω[s k ] we exploit the fact that it is independent of the smooth part, v, of the integrand.Therefore one may judiciously pick a smooth test function, g, which facilitates the computation of the weight.We choose a test function g ∈ C ∞ c (R 2 ) which is radially symmetric and g(0) = 1.We construct a family of weights {ω h } h such that the corrected rule with grid size h integrates exactly our test function g: We define ω[s k ] by the limit Note that since g is chosen to have compact support, T 0 h, N h [s k g] is a summation of a finite number of terms.By choosing g radially symmetric, g(x) = g(|x|), the two-dimensional Cauchy integral R 2 s k (x)g(x)dx can be efficiently approximated to machine precision, e.g. using a Gaussian quadrature, by passing to polar coordinates x = r(cos θ, sin θ)

The case of singular points lying off the grid
In most existing works, see [6,1,12,21], one assumes that the singularity lies in the origin x 0 = 0, or equivalently falls in one of the grid nodes.However, for integrals arising from the IBIM, one must consider the more general case with s k as in (2.6) and x 0 / ∈ hZ 2 .We let x h be the grid node closest to x 0 , satisfying as shown in the left plot of Figure 2. Correspondingly, we define N h,1 (x 0 ) = {x h }, and the punctured trapezoidal rule becomes One can observe from numerical simulations that (2.10) Hence, the weight ω is still independent of v and h, but now depends on the relative position of the singularity with respect to the grid, (α, β).Moreover, the function v is evaluated in x h rather than in the singular point x 0 .Following this observation, we define the first order correction Q 1 h to the punctured trapezoidal rule when x 0 does not fall on the grid as Again, this gives an overall error of size O(h k+2 ).
The weight only depends on the relative position of the singularity with respect to the grid.We therefore set x 0 = 0 and, fixed h and (α, β), shift the grid by (α, β)h.Hence the weight is defined as the limit of the sequence: where The test function g is chosen as in the previous case.The advantages of this choice are going to be the same, e.g. the integral )dx can be computed fast and accurately by passing to polar coordinates.

Second order correction
The goal now is to build a quadrature rule with error O(h k+3 ) for the integrand s k (x − x 0 )v(x) by approximating the quadrature error R h of the punctured trapezoidal rule T 0 h, N h,1 to higher order.We have Since both the singular integral and the punctured trapezoidal rule are linear in s k v, we generalize the ansatz for R h in (2.10) to achieve higher order accuracy by expanding v at a grid node xh close to x 0 : This ansatz requires three weights ω[ We further replace the partial derivatives of v at xh by finite differences of v on {x h,i } p i=1 , x h,1 = xh : where {µ x,i } p i=1 are the finite difference weights for the derivative ∂ ∂x v(x h,1 ).Here, the finite differences involve four grid nodes (p = 4) closest to to the singular point x 0 .They are shown in Figure 2 and given by where xh ∈ hZ 2 is the node such that We remark that α, β are different from the ones for first order correction.
Based on the ansatz above, we define the second order correction Q 2 h to the punctured trapezoidal rule by: (2.15) where We note that as long as the finite differences are first order accurate, using them will not change the formal accuracy of the quadrature rule.The next task is to find a suitable set of weights (2.16) As before, the weight only depends on the relative position of the singularity with respect to the grid.We therefore set x 0 = 0 and, fixed h and (α, β), shift the grid by (α, β)h.The four closest nodes are then Formula (2.16) suggests that we can set up four equations, involving four suitable functions {g j } 4 j=1 : for j = 1, 2, 3, 4 Fixed h and (α, β), this corresponds to imposing that the rule (2.15) integrates exactly the functions Then the weights are found as: We choose the test function , radially symmetric, such that g(0) = 1 and ∇g(0) = 0.This function behaves like the constant function one near 0 and decays to zero smoothly so that the integrand is compactly supported.These properties facilitate efficient and highly accurate numerical approximation of s k g j .We then use Out of the four conditions, the first three translate to the weights correctly integrating any function of the type s k q, q ∈ P 1 , i.e. q two-dimensional polynomial of degree at most one.Three is also the minimum number of points needed in a stencil to have first order accurate ∇v; the fourth node (and consequently the fourth condition) is unnecessary to reach the desired order.It is however useful for it allows us to consider the square fourpoint stencil (2.14) instead of four different three-point stencils necessary to describe the nodes closest to x 0 .
Computing the right-hand side of the linear system involves evaluating with high accuracy integrals with singular integrands By choosing g(x) = g(|x|) radially symmetric, we can write the integral in polar coordinates x = r(cos θ, sin θ): We compute the two factors with high accuracy using Gaussian quadrature.We also reuse the computed values for different parameters (α, β).

Higher order corrections
We now generalize the approach to construct higher order corrections to the punctured trapezoidal rule for (2.7).We expand the ansatz (2.16) used in the previous section to achieve higher order accuracy: where ν ∈ N 2 0 and the weights u ν ∈ R are independent of h and v.This ansatz requires p(p + 1)/2 =: p min weights.We replace the partial derivatives of v at x h,1 by sufficiently high order finite differences.Given p ≥ p min , let be a stencil of p nodes close to x 0 , where x h,1 is such that We approximate the derivatives of v using this stencil: where {µ ν,i } p i=1 are the finite difference weights for the derivative ∂ ν ∂x ν v(x h,1 ).We finally define the p-th order correction Q p h to the punctured trapezoidal rule as where As long as the finite differences for ∂ ν ∂x ν v(x h,1 ) have error ∼ O(h p−|ν| ) they will not affect the formal accuracy of the quadrature rule.
We now have to find a suitable set of weights {ω i [s k ; α, β]} p i=1 for the given s k and (α, β) so that for any smooth function v (2.18) Analogously to Section 2.3, the weights only depend on the relative position of the singularity with respect to the grid.We therefore set x 0 = 0 and, fixed h and (α, β), shift the grid by (α, β)h.Formula (2.18) suggests that we may set up p equations, involving p suitable test functions {g j } p j=1 , to uniquely define the weights {ω i } p i=1 .We proceed as in the previous Section and define the family of weights {ω i,h } p i=1 solution to for j = 1, . . ., p. Fixed h and (α, β), this corresponds to imposing that the rule (2.17) integrates exactly the functions s k (x − x 0 )g j (x − x 0 ), j = 1, . . ., p. Then the weights are found as: Table 1: Correction order and corresponding correction nodes.To increase the order of accuracy by p, the minimum number of nodes to correct is p min = p(p + 1)/2 but more nodes can be used.p ≥ p min is the number of nodes we used in our tests, corresponding to the stencils showed in Figure 3.
We use the function g similar to the one considered in Section 2.3, with the additional conditions that ∂ ν ∂x ν g(0) = 0 for all |ν| ≤ p − 1.This ensures that g is similar enough to the constant function g ≡ 1 near 0.
By choosing the p min functions {g j } p min j=1 equal to g multiplied by the p min monomials of degree at most p − 1 (x i y j , i, j ∈ N 0 , i + j ≤ p − 1) we impose that the method (2.17) integrates exactly all integrands of the type s k q, q ∈ P p−1 , i.e. two-dimensional polynomials of degree at most p − 1.The additional p − p min functions can be chosen for example as g multiplied by two-dimensional monomials of degree higher than p.
We use p ≥ p min because p min may not fit well with standard stencils.Thus it is possible to use more nodes than p min and impose additional conditions.For example in our implementations for first, second, third, and fourth order corrections we used p ≥ p min as shown in Table 1.A visualization of these stencils can be seen in Figure 3.

High order quadratures
In the previous sections we have shown how to deal with integrands of the kind (2.6) wherever the singularity point x 0 may lie, which means we can correct the trapezoidal rule for all terms in the expansion (2.3) 2).If we know these terms explicitly we can build a high order corrected trapezoidal rule for the integral We demonstrate the idea of successive corrections by deriving a second and then a third order accurate quadrature rule.
We first write Lemma 2.4 states that the punctured trapezoidal rule is second order accurate for integrating 0 s.If we apply the first order correction (2.11) to the punctured trapezoidal rule for the first term, we get a second order approximation.The explicit formula, with N h,1 = {x h } as in Section 2.2.2 and relative grid shifts (α 1 , β 1 ), is: This was the approach used in [4], although there (2.11) was used also on the second term instead of the punctured trapezoidal rule.
To achieve third order, we expand s further: We then use the second order correction (2.15) for integrating the first term, first order correction (2.11) for integrating the second, and the (uncorrected) punctured trapezoidal rule for 1 s; by Lemma 2.4 it is third order accurate for 1 s.
We use the set of correction nodes N h,4 (x 0 ) = {x h,i } 4 i=1 and define the corresponding relative grid shift (α 2 , β 2 ).Then the third order accurate rule In general, given the singular function s(x − x 0 )v(x), in order to build a quadrature rule U p h of order p ≥ 2 we need explicitly the first p − 1 (k = 0, . . ., p − 2) terms of the expansion (2.3) and apply to the term to the trapezoidal rule.The punctured trapezoidal rule is used for p−2 s.
x h,3 x h,4 x h = x h,5 x h,6 x h,7 x h,8 x h,9 x h,10 x h,11 x h, 12 Figure 3: Example of correction stencils.The stencils we tested for corrections p = 1 (p = 1 node: yellow square), p = 2 (p = 4 nodes: yellow and red squares), p = 3 (p = 6 nodes: yellow and red squares, and green circles), and p = 4 (p = 12 nodes: yellow and red squares, green circles, and cyan stars).The singularity node is x 0 (red circle).The nodes are {x h,i } 12 i=1 , and We can find an explicit expression for the quadrature rule U p h by specifying the stencils we use for the correction nodes.We denote by N n,p(p) the stencil of p(p) correction nodes to increase the order by p.We assume that the stencils are increasing: N n,p(p) ⊂ N n,p(p+1) .For example in our tests we took p(1) = 1, p(2) = 4, p(3) = 6, p(4) = 12, and 12 .This is shown in Table 1 and Figure 3.We call α p , β p the parameters describing the shift with respect to x 0 of the stencil of p(p) nodes: In Section 4.1 we show tests for the quadrature method (2.24) by combining first, second, third, and fourth order corrections.

Approximation and tabulation of the weights
Given functions of the kind (2.6), Fixed k ≥ 0 and (α, β), we write the function s k (specifically its factor φ k ) using its Fourier series: where {a j } ∞ j=0 and {b j } ∞ j=1 are the Fourier coefficients of φ k .Then, by linearity of the weights with respect to s k , we can write them as with i = 1, . . ., p.We can then approximate and tabulate the weights ω i [s k ; α, β] in the following way.We fix a stencil of parameters {(α m , β n )} m,n around (α, β) and basis functions {c m,n (α, β)} m,n such that we can approximate a function f : We let N be the number of Fourier modes used to approximate the weights.Then, given φ k , we first find the 2N + 1 coefficients a 0 , {a j , b j } N j=1 by using the Fast Fourier Transform.Then, and we can approximate the weight for (α, β) via So, for all expansion terms k = 0, 1, . . ., p − 2 used in (2.23), and the corresponding corrections , we need to compute and store the weights for the following constant and trigonometric functions, and all m, n in the stencil for (α, β).

Evaluating layer potentials in the implicit boundary integral formulation
We apply the high order quadrature methods from Section 2 to layer potentials used in Implicit Boundary Integral Methods (IBIM).To make the exposition clear we adopt the following convention.
Notation 3.1.We distinguish between variables in R 2 and R 3 by using boldface variables for vectors in R 2 and boldface variables with a bar for vectors in R 3 .For example, x ∈ R 2 and x ∈ R 3 .Moreover, for a vector y = (y 1 , y 2 ) ∈ R 2 and scalar y 3 ∈ R we frequently write ȳ = (y, y 3 ) to mean the vector (y 1 , y 2 , y 3 ) ∈ R 3 .For example, when f : R 3 → R m , we use the notations f (ȳ) ≡ f (y, y 3 ) ≡ f (y 1 , y 2 , y 3 ).
We consider the general form of a layer potential on a smooth, closed and bounded surface with K defined by one of the following kernels: (single-layer, SL) : (double-layer conjugate, DLC) : In (3.2), the vector nx is the normal vector to Γ at x * , pointing into the unbounded region R 3 \ Ω, where Ω is the bounded region enclosed by Γ. Analogously ny is the normal vector to Γ at ȳ.In preparation for the formulation of the implicit boundary integral methods we first define d Γ : R 3 → R to be the signed distance to the surface such that d Γ is negative inside Ω.Moreover, we let P Γ : R 3 → Γ be the closest point mapping that takes ȳ to a closest point on Γ: Let C Γ ⊂ R 3 be the set containing the all the points that have non-unique closest points on Γ.The reach τ Γ of Γ is defined as It depends on the local geometry (the curvatures) and the global structure of Γ (the Euclidean and geodesic distances between any two points on Γ).The reach is positive for some α > 0. The closest point mapping P Γ is invertible in the tubular neighborhood In this paper, we will assume that Γ is a closed bounded C 2 surface so that the mean and Gaussian curvatures are defined everywhere on the surface.Consequently, when ȳ lies within the reach of Γ, we have the explicit formula The surface integral (3.1) can then be reformulated into an equivalent volume integral using the Implicit Boundary Integral Methods [7,8], The "delta" function is defined as δ Γ,ε (ȳ where with H(ȳ) and G(ȳ) denoting respectively the mean and Gaussian curvatures of Γ η := {x ∈ R 3 : d Γ (x) = η} (see e.g.[8]).For |η| < τ Γ , J η is bounded away from zero.Moreover, unit mass.This is achieved by using δ ∈ C ∞ c (R) compactly supported in (−1, 1) with R δ(η)dη = 1.It turns out that for any positive ε, smaller than the reach of Γ, the IBIM is equal to the original layer potential for all x ∈ R 3 , If the surface Γ is smooth, the closest point mapping is also smooth ([3], Ch. 7, §8, Thm 8.4).As a consequence, if ρ is a smooth function on Γ and Γ is smooth, then ρ(P Γ (ȳ))δ Γ,ε (ȳ) is a smooth function R 3 , compactly supported in T ε .

Singular integrand in three dimensions and correction plane by plane
In this section, we will construct high order quadratures for 3) from the two-dimensional corrected trapezoidal rules (2.22).The three dimensional quadrature rules will be defined as the sum of integration over different coordinate planes, where the two dimensional corrected trapezoidal rule from the previous Section is applied to approximate the integration over each plane.The particular selection of the coordinate planes depends on the normal vector of Γ.
Without loss of generality, we consider a target point, x * = (x * , y * , z * ) ∈ Γ, at which the surface normal is n = (n 1 , n 2 , 1).The way to treat other cases are explained in Section 3.1.1.We denote the integrand in (3.3) by f , To approximate (3.3) the standard trapezoidal rule is first applied in the z-direction.With the grid points z k = kh, we get where we used the fact that f is compactly supported in T ε .We note that f is singular along the line since P Γ (ȳ 0 (z)) = x * for all z.Therefore, f (•, •, z) is singular at one point for each fixed z, by the assumption on n.Below we will derive the form of this singularity and we will show that it is of the same type (1.2) as considered in Section 2. See Figure 4 for an illustration of the line singularity, and an example of the singular behavior.Hence, we can use the corrected trapezoidal rules to approximate each integral in the sum in (3.6).
To connect back to the notation in Section 2 we write ȳ = (y, z) and ȳ0 (z) = (y 0 (z), z).Then we factorize f , for a fixed z, as where Note that the type of singularity for s depends on the properties of Γ at the target point (such as principal curvatures, principal directions, normal).Moreover, s depends smoothly on z.
We then use the corrected trapezoidal rule ) to compute the integrals on each plane, We denote by y ∆ (z) and (α 1 (z), β 1 (z)) the closest grid node to y 0 (z) and the relative grid shift parameters respectively, as defined in Section 2.2.2, and define N z h,1 (y 0 ) := {y ∆ (z)}.We also denote by {y ∆,i (z)} 4 i=1 and (α 2 (z), β 2 (z)) the four grid nodes surrounding y 0 (z) and the relative grid shift parameters respectively, as defined in Section 2.3, and define N z h,4 (y 0 ) := {y ∆,i (z)} 4 i=1 .From the definition in (3.9) we can compute the expansion (2.3) with q = 1 and find with s k (y; z) = |y| k−1 φ k (y/|y|; z), k = 0, 1.The expressions for s k are given in Theorem 3.2 below.We can then apply the additive splitting (2.21): {s(y − y 0 (z); z) − s 0 (y − y 0 (z); z)} v(y, z).Then the three-dimensional third order method V 3,z h , obtained by applying U 3 h plane-byplane along the z-direction, is given by If we apply the two-dimensional rule U 3 h plane-by-plane along the x-or y-direction we obtain the corresponding rules V 3,x h and V 3,y h respectively.These cases are discussed in the following Section 3.1.1.

Plane-by-plane correction for different normal directions
The normal direction n directly affects the decomposition of a three dimensional Cartesian grid into union of planes, on which we apply the new correction.We identify the dominant direction of n = (n x , n y , n z ), and discretize the volumetric integral along that direction.If the dominant direction is n z , the setup is the one described above.If it is n y , we discretize along the y-direction, and if it is n x , we discretize along the x-direction.
We shall use V p h for the general three-dimensional p-order corrected trapezoidal rule, and using the division presented for the change of coordinates, we define it as: (3.12)

Expansions of layer kernels
In this section we will analyze and expand the singular functions defined in (3.9) when K are the Laplace kernels (3.2): (SL) : (3.13) The approach developed in Section 2 requires analytic formulae of the expansions.This means that in order to adopt the third order quadrature rule (3.12) for the implicit boundary integral defined in (3.3), one needs explicit analytical expressions for the first two expansion functions in (3.10) related to the singular functions above.Through a third order approximation of the surface near the target point x * we find these functions, which are given in the following theorem.Theorem 3.2.Let x * ∈ Γ be the target point.Suppose that the normal n at x * satisfies nT ēz = 0, and that ȳ0 (z) ∈ T ε .Then, there is an r > 0, depending on z, such that all the singular functions defined in (3.13) can be written in the form where X ∈ C ∞ ((−r, r) × S 1 ).Moreover, the functions s X 0 (y; z) and s X 1 (y; z) in the expansion (3.10) are where ŷ = y/|y| and ψ j , ξ j and ξ1 are given explicitly in Section 3.2.4.
In order to use the expansions in the theorem in our quadrature method, we need to be able to evaluate the functions ψ j , ξ j and ξ1 .They depend on the local behavior of Γ at the target point x * , more precisely on the principal directions and curvatures, and the third derivatives of the function whose graph locally describes Γ.In Appendix B it is described how those quantities can be computed numerically using the closest point mapping.
In the subsequent subsections we will prove Theorem 3.2.First, in Section 3.2.1,we rotate the frame of reference and look at Γ locally as the graph of a two-dimensional function.Second, we expand the expressions we obtained around y = 0 in Section 3.2.2 and apply a general lemma to show (3.14) in Section 3.2.3.Finally, we use the expansions to derive expressions for s X 0 (y; z) and s X 1 (y; z) in Section 3.2.4.

Expressions of the layer kernels via the projection mapping
Let x * ∈ Γ be the target point.At x * we denote the surface principal directions τ1 , τ2 , the normal n, and the principal curvatures κ 1 , κ 2 .We introduce the principal basis B = ( τ1 , τ2 , n) and the notation The basis vectors used here are assumed to be normalized.If x are the coordinates in the B-basis for the point x in the canonical basis (ē x , ēy , ēz ), we denote by Q the (orthogonal) change of basis matrix, satisfying The surface Γ can now be parameterized locally in the B-coordinates.More precisely, in a neighborhood of the origin, The constant L depends on the maximum curvature of Γ and can be taken to be independent of x * .Moreover, since P Γ is smooth in the tubular neighborhood T ε of Γ, the mapping (y , z ) → P Γ (x * + (y , z ) B ) is smooth for (y , z ) ∈ T ε := {(y , z ) ∈ R 3 : x * + (y , z ) B ∈ T ε }.Therefore, for (y , z ) ∈ M L = T ε ∩ (I L × R), with L possibly smaller than L , we can use the B-basis and f , to write the closest point mapping as for some smooth function Clearly h(0, z ) = 0, which guarantees that L > 0.
We now write (y + y 0 (z), z) as a point in the B-basis centered in the target point x * , For (y , z ) ∈ M L we can then write the numerators and denominators of the layer kernels (3.13) using (3.17) and the orthogonality of Q: ) We next have to find how (y , z ) depends on y and z.From the definitions above we have Since ȳ0 (z) − x * is parallel to the normal n by definition, we can express this as .

Expansion of f and h
By the definition of the B-basis, the function f introduced above in Section 3.2.1 satisfies The Taylor expansions up to second order for f and ∇f are then given by where B is the third order trilinear term, and C is its bilinear gradient.With y = (x, y), they are given by B(y, y, y) := We next need to expand h.It is given by the following lemma, the proof of which can be found in the Appendix A.3.
with M given in (3.24).For (y , z ) ∈ M L the matrix is well-defined.The function and the Taylor expansion of h can be written in the form where C is defined in (3.26).

General form of the kernels
We now have expressions (3.23) of the kernels and expansions around y = y p = 0 of h and f .The next step is to prove (3.14), i.e. that the three kernels in (3.23) can all be written in the form |y| −1 (|y|, y/|y|).To do this we use the following lemma, a proof of which can be found in Appendix A.4.

Kernel expansions
The expansion of the kernels is based on the expansions of f in (3.25) and h in (3.28).We will skip most tedious intermediate calculations and focus on the end results.We recall that In the first step we expand the functions f (h(y , z )), D(z ) and h(y , z ) as functions of y, instead of y p and y as before.We get In this step we used the fact that χ j and ξ j are homogeneous of degree j + 1 and j + 2 respectively, so that χ j (y) = χ j (ŷ)|y| j+1 and ξ j (y) = ξ j (ŷ)|y| j+2 .From these expansions for f and h we obtain furthermore that Here ξ1 is homogeneous of degree three.Using (3.29) one can finally deduce the expansions of the kernels in (3.23).This concludes the proof of Theorem 3.2.We note that the matrix A and vector d contain elements of the principal directions and normal at the target point; see (3.16) and (3.21).The matrices D 0 and M are built from the principal curvatures of Γ, and the functions B and C contain the third derivatives of f ; see (3.25) and (3.26).In Appendix B we show how to numerically compute the information about the surface in the target point (κ 1 , κ 2 , τ1 , τ2 , and the third derivatives of f , f xxx , f xxy , f xyy , f yyy ) using the projection mapping P Γ and its derivatives.

Requirements for order higher quadratures for the singular IBIM integrals
Given the class of singular integrands, the main obstruction to obtaining higher order quadratures using the proposed approach is the smoothness of the surface.When applying the proposed method in the IBIM formulation using uniform Cartesian grids, one needs firstly a sufficiently accurate approximation of the distance function to the surface, d Γ , or the projection, P Γ , on the grid nodes.The construction of these functions are application dependent, but general methodologies do exist, see e.g.[13].If the surfaces are reconstructed on a grid by a level set method, then typically one does not expect that d Γ be more than 4th order accurate in the grid spacing due to the limitation imposed by commonly used level set reinitialization algorithms [13].This may cause a main bottleneck in practice.Then one needs to extract the surface's geometrical information from finite differences of d Γ or P Γ -in this paper, the related quantities to be approximated are the partial derivatives of f defined in (3.26),where f is defined in (3.5).In Appendix B the reader will find more details.
When the surface is sufficiently smooth, it has a non-zero reach; i.e. d Γ is smooth within T τ Γ for some τ Γ > 0. The Cartesian grid inside T τ Γ should be sufficiently dense to support the finite difference stencil around any node inside T ε , where ε ≤ τ Γ .Thus higher order approximations require denser grids around the surface to support the wider finite difference stencils used in high order finite differences.For example the second order corrected rule V 2 h needs curvature information which is obtained through a centered 5 point three-dimensional stencil.This implies that one needs accurate d Γ or P Γ within the distance of ε + 2h to the surface, and that ε + 2h should be smaller than the reach τ Γ .Analogously, the third order information about the surface needed for V 3 h is obtained using a 5 × 5 × 5 stencil around each node, which leads to the bound ε + 2 √ 2h < τ Γ .If the surface geometry varies "wildly", we envision that the proposed method should/could be generalized to multi-resolution gridding for efficiency.
We present an example supporting the above discussion in Section 4. We also refer the interested readers to the results and discussion in the recent paper [5], for an application of the proposed quadratures in computing the electrostatic potentials of large molecules in a solvent.

Numerical tests
In this Section we test the corrected trapezoidal rules derived in Section 2 and Section 3. In Section 4.1 we test the rules Q p h for integrating functions of the kind s k v from Section 2.4, and then the general rules U p h for integrating s v from Section 2.5.In Section 4.2 we test the third-order accurate quadrature rule V 3 h derived for the three-dimensional layer potentials discussed in Section 3.
Table 2: Maximum of the weights.Largest weight max p i=1 ω i in the stencil N h,p for different correction orders p = 1, 2, 3, 4, and different singularity order k = 0, 1, 2. The weights correspond to the ones used in the tests shown in Figure 5.All weights are nonnegative.

Corrections to the punctured trapezoidal rules in two dimensions
The quadrature rules discussed in Section 2 have been developed to correct any function of the kind where v is a smooth function, and then composite rules have been constructed to correct functions which can be expanded as where s(x) = s 0 (x) + s 1 (x) + s 2 (x) + . . . .

The function H
(1) α is the Hankel function of the first kind of degree α, and indicates the real part of a complex number.Although formally v is not compactly supported, it is smaller than the numerical machine precision outside [−2, 2] 2 , which we use as integration domain.
In Figure 5 we plot the difference between approximation values for grid sizes h and h/1.5, obtained for the four different quadratures Q p h , p = 1, 2, 3, 4, and the punctured trapezoidal rule T 0 h .The order of accuracy shown for integrating s k v, k = 0, 1, 2, is k+1 for the punctured trapezoidal rule and k + p + 1 for the quadrature Q p h , as expected.The error constant is determined by the value of (α, β), and in our tests we fixed (α, β) = (0.81, 0.46).The stencils used for the different quadratures are represented in Figure 3.The weights for the quadratures are all non-negative.Their maximum values are shown in Table 2.They are of moderate size also for the high order corrections.) with porder correction Q p h .For k = 0, 1, 2 (left, center, and right figures respectively) we present the difference between values obtained from grid sizes h and h/1.5, with the different methods.As expected, the order of accuracy is k + p + 1 where p is the order of the correction.
In order to test the general quadrature rule (2.24) we used the function where x =|x|(cos(ψ(x)), sin(ψ(x))), In Figure 6 we plot the difference between values obtained with grid sizes h and h/1.5 for the four different quadratures U p h (2.24) and the punctured trapezoidal rule T 0 h .The order of accuracy shown is 1 for the punctured trapezoidal rule and p for the quadrature U p h , which is what was expected.The error constant is determined by the value of (α, β).In all our tests we fixed (α, β) = (0.81, 0.46).The stencils used for the different quadratures Q k h needed to compose U p h are the same as the previous test, represented in Figure 3.

Evaluating the layer potentials in the IBIM formulation
We demonstrate the convergence and accuracy of the proposed quadrature rules by evaluating the single-layer, double-layer, and double-layer conjugate potentials with some smooth The integrals are first extended to the tubular neighborhood of T ε , as in (3.1) using the compactly supported C ∞ averaging function here a ≈ 7.51393 normalizes the integral R δ(η)dη to 1.

A numerical study on a smooth surface
The surface chosen for the tests is a torus, centered at a randomly chosen point in 3D and rotated with randomly chosen angles along the x-, y-and z-axes.This is to avoid any symmetry of the uniform Cartesian grid which can influence the convergence behavior.This setup includes all the essential difficulties one may encounter when applying the proposed method to a smooth surface: non-convexity, finite reach from the geometry, and asymmetry in the discretized system.The torus is described by the following parametrization where R 1 = 0.7, R 2 = 0.2, C imposes a translation, and The known density function ρ used in the test is defined using the parametrization of the torus: 196 sin θ − 0.29837 cos φ sin θ + 1.128 sin φ cos θ .
We present the errors computed for a sequence of grid size values {h i } i , where we used as reference value half of the smallest grid size h min = 1 2 min i h i .We our third order rule V 3 h (3.12).Moreover we compared with the previously developed second order rule, denoted by V 2 h , from [4].In the presented simulations, we take the component n z of n to be dominant if | tan θ| < √ 2, where n/|n| = (sin θ cos φ, sin θ sin φ, cos θ).If instead | tan θ| ≥ √ 2 and | tan φ| ≥ 1, we take n y to be dominant, and if | tan θ| ≥ √ 2 and | tan φ| < 1 we take n x to be dominant.We used θ and φ to determine the dominant direction because of their extensive use in the rest of the code.
At each target point x * , the total error is the sum of the errors of the two-dimensional rule applied on each plane.Recall that under the IBIM formulation, the kernel is singular along the surface's normal line passing through x * , and the singularity of the kernel on each plane lies at the intersection of the surface normal line and that plane.Since the normal lines of the surface generally do not align with the grid, the position of the singular point relative to the grid tends not to lie on any grid node.Recall further that the parameters α, β are used to described the position of the singular point relative to the closest grid node on the plane, and the error constants depend on them.Those parameters may change abruptly between planes, depending on which grid node in the plane is closest to the singular point.The closest grid nodes to each surface normal line certainly are expected to exhibit jumps as one refines the grids (decreases h).Thus, as noted in [4], the errors (4.5) as functions of h are generally not smooth.Consequently we cannot see a clear slope.
To show the overall convergence behavior we average the errors, defined in (4.5), over 20 target points, randomly chosen.The results can be seen in Figure 8.In the left column we present the averaged errors.In the right column we present a scatter plot of the errors at all the target points.We additionally highlight the errors corresponding to two specific target points to showcase an "average" error behavior (green line) and a "bad" error behavior (magenta line).By construction of the quadrature rule (3.12) we expect it to be third order accurate in h.However from the plots we observe order of accuracy ≥ 3.5.We conjecture that an additional cancellation of errors occurs when adding the results from each plane (see [4] for a related discussion regarding V 2 h ).A rigorous analysis of this behavior is beyond the scope of this article.
Of course to test our algorithms, we retain no information about the parametrizations.The test torus is represented only by d Γ and P Γ on the given grid.Figure 7 shows the torus that we use and the points used in the quadrature rule for a given grid configuration.We use fourth-order centered differencing of P Γ on the grid to approximate the Jacobian J Γ (see [8]).We also use fourth-order centered differencing of P Γ to find the third derivatives of f needed for the functions B and C in (3.26), as they are related via a linear system (see Appendix B).

A numerical study on a more complicated surface
We present a test of the quadrature rule applied to IBIM for a more complicated surface, shown in Figure 9.The surface represents the solvent-molecule interface of a complex biomolecular system immersed in a solvent [22].A level set representation of the surface is generated using VISM [19] by the authors of [22] on a 512 3 Cartesian grid.We compute the relative error in the double-layer identity using the proposed method.The relative error is defined as: with p = 0 denoting the punctured trapezoidal rule.In our setup, inherited from the shared data set, the signed distance function is accurate up to distance ≈ 9h, which is minimally adequate for the application of V 2 h .The computed values are E 0 (h) = 0.0159 for the punctured trapezoidal rule and E 2 (h) = 0.00142 for the second order corrected rule.Furthermore, when applying V 3 h , we notice that the resulting pointwise errors oscillate across grid nodes, xj , and do not appear to be smaller than those computed by V 2 h .On some xj , the error even appear to be larger that those computed by the punctured trapezoidal rule.This is expected because the grid does not yet resolve the fine geometry in this surface (notice in particular the narrow separation of the two connected components).
where we cut out the singularity point by multiplying by 1 − ψ around 0; the scaling by h ensures that, for fixed h, only the node in the singularity point is cut out.This allows us to split the error of the punctured trapezoidal rule as .
We will consider the two terms (I), (II) separately, and prove that both can be bounded by Ch j+n .
(I): Given the compact support of ψ, the integral is reduced to an integral over {|x| ≤ h}: since |x| j is integrable as j ≥ 1 − n.We have proven the estimate for the first term.
(II): For the second term, knowing that the volume of the fundamental parallelepiped of the lattice V := (hZ) n is h n and that the dual lattice is V * = (h −1 Z) n , we use the Poisson summation formula: Then the error in (II) is: Using integration by parts separately on each of the variables, we find For the Laplacian operator applied q times we therefore have We use this result to find an expression we can bound using Lemma A.2; given an integer q, we find Then the series of Fourier coefficients is The series converges if 2q > n, and the leading order is h j+n if 2q ≥ j + n, so by taking q ≥ max(1 + n/2, (n + j)/2), we find the result sought.Combining the results for (I) and (II), we find the bound This proves the theorem.
We use the notation x = (x 1 , x 2 , . . ., x n ) = n l=1 x l e l , and indicate with e l the l-th element of the standard R n basis.
, where ψ is such that Proof.Given β ∈ N n 0 , we first prove that there exist functions f β : R × We prove this by induction.The induction base β = 0 is true because . For the induction step we assume that (A.4) is true for β and prove it for β + e l : . By computing the derivative we find ) the same is also true for f β+e l .The next step is to expand the derivative in (A.3) and use (A.4), and then bound it: We use the properties of ψ, and the compact support of f ν .Let L > 0 be such that ∀ν ≤ β, supp f ν is contained in the ball B L (0).Note furthermore that the derivatives of ψ are compactly supported in the annulus {x ∈ R n : 1 2 ≤ |x| ≤ 1}.From this we can say that We use these bounds in the evaluation of the integral, and after passing to polar coordinates we arrive at (A.3) via The lemma is proven.

A.2 Proof of Lemma 2.3
For any u ∈ S 1 , we expand around r = 0 and write the remainder in integral form: (r, u) = q j=0 1 j! ∂ j r (0, u)r j + r q+1 q! and the result follows upon evaluating at y = y p = 0 and using (3.24).Since h is smooth on M L the matrix D must thus be well-defined.
For the second order term in the Taylor expansion, we write y p = (y 1 , y 2 ), h = (h 1 , h 2 ) T and F = (F 1 , F 2 ) T .We then get for j = 1, 2, From the expressions above we have that ∂F(0) ∂yp = D −1 (z).Therefore, evaluating at y = y p = 0, yields For the second function form, let r(x) := p(x) T ḡ(x); then ∇r(x) = ḡ(x) T Dp(x) + p(x) T Dḡ(x).Using the hypothesis p(0) T Dḡ(0) = 0, we write the expansion of r around x = 0 using the integral form of the remainder: From the hypotheses on the smoothness of ḡ and p, 2 is C ∞ ((−r 1 , r 1 ) × S m−1 ) and the result is proven.

B Computation of the derivatives of the local surface function
In this Section we will show how to find numerically the derivatives of f in the Implicit Boundary Integral Methods setting of Section 3. The derivatives are needed to evaluate the functions B and C of (3.26), which are used in the approximated kernels (3.15).
The first derivatives and the mixed second derivatives are zero by construction, so we will show how to find the pure second derivatives and all the third derivatives.
The pure second derivatives of f at P Γ (z), f xx , f yy , are the principal directions κ 1 , κ 2 of Γ at P Γ (z).We find the principal curvatures g 1 , g 2 of Γ η in z via the Hessian of d Γ at z: T where τ1 , τ2 are the principal directions and n is the normal to Γ in P Γ (z).In practice, the values of either P Γ or d Γ are given on the grid nodes.The principal directions and curvatures are computed from eigendecomposition of third order numerical approximations of the Hessian, H d Γ .Alternatively, one can obtain this information from the derivative matrix of P Γ , see [8].Then the following relation lets us find the principal curvatures κ 1 , κ 2 from g 1 , g 2 and η: The third derivatives of f can be found by computing the second derivatives with respect to y of h(y , z ) from Section 3.2.1.By differentiating twice (A.5) with respect to y = (x, y) with h(y , z ) = y p = (h 1 , h 2 ) and evaluating in y = 0, we find the following two linear systems: We find the first and second derivatives of h(y , z ) by computing the derivatives of P Γ in z and applying a change of basis transformation.
In the B basis, these points are expressed as vijk = x * +( wijk ) B , where wijk = Q −1 (v ijk − x * ).We apply finite differences (central differences of 4th order in this case) to the component of the nodes wijk = (X ijk , Y ijk , Z ijk ) to compute

Figure 2 :
Figure 2: Singularity unaligned to the grid.The parameters α, β are used to characterize the position of the singularity point x 0 (red circle) relative to the grid in two different settings.Left plot (first order correction): position of the singularity point relative to the closest grid node x h (yellow square).Right plot (second order correction): position of the singularity point relative to the four surrounding grid nodes x h,i , i = 1, 2, 3, 4 (red squares except x h,3 = x h which is yellow).

Figure 4 :
Figure 4: IBIM kernel singular behavior.The kernel K(x * , P Γ (ȳ)), for fixed x * ∈ Γ and ȳ ∈ T ε , is singular along the normal n to x * .The left figure illustrates how the kernel becomes singular for ȳ approaching any point of the line passing through x * with direction n.The center plot shows the double-layer conjugate kernel K = ∂G ∂nx plotted on a plane with fixed z, ȳ ∈ {(x, y, z) : x, y ∈ R}.The function will then have a point singularity in y 0 (z), and we plot the kernel for (x, y) close to y 0 (z).The right most plot shows the same kernel, multiplied by |y − y 0 (z)|.

Figure 5 :
Figure 5: Correction of s k in two dimensions.Error from integrating s k (4.1) with porder correction Q ph .For k = 0, 1, 2 (left, center, and right figures respectively) we present the difference between values obtained from grid sizes h and h/1.5, with the different methods.As expected, the order of accuracy is k + p + 1 where p is the order of the correction.

Figure 6 :
Figure 6: Corrected trapezoidal rules for a general function s in two dimensions.Corrected trapezoidal rules U p h for p = 2, 3, 4, 5 using additive splitting (2.24) for the function f = s v with singular integrand s (4.2).The first p − 1 terms of the expansion (2.3) (s k , k = 0, 1, . . ., p − 2) are needed to use U p h .In the plot we see that the punctured trapezoidal rule T 0 h has first order accuracy, and the corrections p h have order of accuracy p as predicted.
is the composition of three rotation matrices; Q x (a), Q y (a), and Q z (a) are the matrices corresponding to a rotation by an angle a around the x, y, and z axes respectively.The parameters used for the translation and the rotations were:C = 0.5475547095598521, 0.6864792402110276, 0.3502726366462485 • 10 −1 , a = 0.199487 • 10 1 , b = 0.2540979476510170 • 10 1 , c = 0.4219760487439292 • 10 1 .

Figure 7 :
Figure 7: Torus test surface.Left: the torus used in the tests.Right: the torus and the projections of the Cartesian grid nodes inside the tubular neighborhood T ε .The projected nodes serve as the quadrature nodes.

Figure 8 :
Figure 8: Errors in the evaluation of the three Laplace layer potentials.The errors (4.5) are computed for 20 randomly chosen target points on a tilted torus.The plots in the left column show the mean of the 20 errors.The plots in the right column show the scatter plot of the 20 target points.In the right plots we additionally highlight the behavior of two specific target points, to showcase a "bad" error (magenta line) and an "average" error (green line).

Figure 9 :
Figure 9: A solvent-molecule interface.The surface is computed by the VISM method for biomolecular system p53-MDM2 (PDB ID 1YCR) [9] from the Protein Data Bank (PDB).
+ 2f xxy xy + f xyy y 2 f yyy y 2 + 2f xyy xy + f xxy x 2 xxx , f xxy , f xyy , f yyy are the third order derivatives of f evaluated in 0.