Ubiquitous evaluation of layer potentials using Quadrature by Kernel-Independent Expansion
UUbiquitous evaluation of layer potentials using Quadrature byKernel-Independent Expansion
Abtin Rahimian · Alex Barnett · Denis ZorinAbstract
We introduce a quadrature scheme —
QBKIX — for the high-order accurate evaluation of layerpotentials associated with general elliptic
PDE s near to and on the domain boundary. Relying solely onpoint evaluations of the underlying kernel, our scheme is essentially
PDE -independent; in particular, noanalytic expansion nor addition theorem is required. Moreover, it applies to boundary integrals withsingular, weakly singular, and hypersingular kernels.Our work builds upon Quadrature by Expansion (
QBX ), which approximates the potential by an ana-lytic expansion in the neighborhood of each expansion center. In contrast, we use a sum of fundamentalsolutions lying on a ring enclosing the neighborhood, and solve a small dense linear system for theircoefficients to match the potential on a smaller concentric ring.We test the new method with Laplace, Helmholtz, Yukawa, Stokes, and Navier (elastostatic) ker-nels in two dimensions ( ) using adaptive, panel-based boundary quadratures on smooth and cornerdomains. Advantages of the algorithm include its relative simplicity of implementation, immediate exten-sion to new kernels, dimension-independence (allowing simple generalization to ), and compatibilitywith fast algorithms such as the kernel-independent FMM . The boundary integral method is a powerful tool for solving linear partial differential equations (
PDE s)of classical physics with piecewise constant material coefficients, with applications including electro-magnetic scattering, molecular electrostatics, viscous fluid flow, and acoustics. It involves exploitingGreen’s theorems to express the solution in terms of an unknown “density” function defined on the do-main boundaries or material interfaces, using the physical boundary condition to formulate an integralequation for this density, and finally obtaining a linear algebraic system via Galerkin, Nyström, or otherdiscretization. Compared to commonly used differential formulations, boundary integral methods havea number of advantages: decreasing the dimension of the problem that needs to be discretized, avoidingmeshing the volume, and improving conditioning. For instance, the integral equation can often be cho-sen to be a Fredholm equation of the second kind, resulting in a well-conditioned linear system whichcan be solved by a Krylov subspace methods in a few iterations. All these considerations are particularlyimportant for problems with complicated and moving geometries [ ] .The main difficulty in using boundary integral methods is the need to evaluate singular and nearly-singular integrals: (i) Evaluating system matrix entries requires evaluation of the potential on the sur- Abtin Rahimian · Denis ZorinCourant Institute of Mathematical Sciences, New York University, New York, NY 10003 E-mail: [email protected], E-mail:[email protected] BarnettDepartment of Mathematics, Dartmouth College, Hanover, NH 03755 E-mail: [email protected] a r X i v : . [ m a t h . NA ] D ec Abtin Rahimian et al. − − − − − − (a) − − − − − − (b) Fig. 1.1:
Evaluation error plotted in the solution domain due to approximating the Laplace double-layer potential Eq. (1.1) usinga quadrature designed for smooth functions. Logarithm of absolute error, log | ˜ u ( x ) − u ( x ) | , where u is the true solution and ˜ u is the discrete approximation using smooth quadrature is plotted for the case of constant density φ ≡ . (a) shows compositequadrature with M = (left) or M = (right) panels each with q = Gauss–Legendre nodes. (b) shows the global compositetrapezoid rule with N = (left) or N = (right) nodes. face, which involves a singular integral; (ii) Once the density is solved for, the desired solution must stillbe evaluated in the form of a potential. As an evaluation point approaches the boundary of the domain,the peak in the resulting integrand becomes taller and narrower, giving rise to what is referred to as a near-singular integral. The result is an arbitrarily high loss of accuracy, if the distance from points to thesurface is not bounded from below, when a quadrature scheme designed for smooth integrands is used [
2, Section 7.2.1 ] and [ ] .Figure 1.1 illustrates the near-singular evaluation of the solution u of the Dirichlet Laplace equationin a simple smooth domain, which is represented by the double-layer potential u ( x ) = π (cid:90) Γ ∂∂ n y log 1 (cid:107) x − y (cid:107) · φ ( y ) d s ( y ) , (1.1)where φ is the density defined on the boundary Γ . The growth in error as x approaches Γ is apparentin all four plots (showing panel-based and global quadratures with different numbers of nodes N ).Although the width of the high-error layer near the boundary shrinks like 1 / N [ ] , the error alwaysreaches (cid:79) ( ) at the boundary. The goal of this paper is to present a flexible scheme that handles bothtasks (singular and near-singular evaluation) to high-order accuracy in a kernel-independent (i.e., PDE -independent) manner.
Related work . Designing quadrature schemes for singular and near-singular integrals has a long andrich history [ ] . Until recently, the quadrature methods were designed specifically for either on-surface evaluation or near-surface evaluation. Many of the on-surface integration quadrature are specificto a certain type of kernel (singularity), e.g., log | r | in or 1 / | r | in [ ] ;the former case is reviewed in [ ] .A popular method for on-surface quadrature is the product integration (in , for the global trapezoidrule see [
2, Section 4.2 ] or [
38, Section 12.3 ] , and for panel-based rules see [ ] ). In this context, ananalytic convolution of the kernel with each function in some basis set is found, reducing evaluation ofthe integral to projection of the boundary density onto that basis set.Another approach for on-surface evaluation is singularity subtraction, where the integrand is modi-fied by subtracting an expression that eliminates its singularity [
16, Chapter 2 ] and [ ] . However,this leaves high-order singularities in the kernel which makes the higher derivatives of the kernels un-bounded, limiting the accuracy of the quadrature scheme. Alternatively, for weakly singular kernels, onecan use transformations to cancel the singularity by the decay of area element (e.g., in using Duffytransformation [ ] or polar coordinates) [ ] . To achieve a high con-vergence order, these methods need some form of partition of unity so that a high-order polar patch canbe constructed around each point [ ] . biquitous evaluation of layer potentials using QBKIX One can also regularize the kernel and then exploit quadrature schemes for smooth functions [ ] . However, to achieve higher accuracy, the effect of regularization needs to be corrected by usinganalytic expressions (e.g., asymptotic analysis) for the integrand [ ] . Finally, there exist special high-order quadrature schemes for domains with corners, either via reparametrization [ ] , panel-wisegeometric refinement [ ] , or by custom generalized Gaussian quadratures [ ] .We now turn to near-singular integrals (evaluation close to the surface), which has traditionallybeen handled as a distinct task [ ] . Beale and coauthors [ ] use regularizationmethods to remove the singularity of the integral. To correct the error introduced by the regularization,they perform asymptotic analysis and find correction expressions. Some authors used singularity cancel-lation (e.g., using local polar coordinates) in evaluating near-singular integrals [ ] . Interpolationalong carefully-chosen lines connecting distant points (where a smooth quadrature is accurate) to anon-surface point has also been successful [ ] .Recently, unified approaches to on-surface and close evaluation have been proposed, the first beingthe Laplace high-order global and panel-based quadratures of Helsing and Ojala [ ] . This approachhas been extended to near-singular Stokes single- and double-layer kernels with global [ ] and panel-based [ ] quadrature. The use of local expansions — analytic separation of variables to the PDE solutionsanalogous to a Taylor series in the complex plane — for the evaluation of integrals near the boundarywas introduced in [ ] .In this scheme, a refined smooth quadrature is needed to accurately evaluate the expansion coeffi-cients via the addition theorem. It was observed that the expansion can also be used to evaluate at targetpoints on the boundary of the domain, if certain conditions are satisfied [ ] ; this was used to constructa unified quadrature scheme — Quadrature by Expansion ( QBX ) — for near and on-surface evaluationof integrals [ ] . Racch [ ] recently showed how to efficiently combine QBX evaluations with the fastmultipole method.However, powerful as they are,
QBX schemes require both a local expansion and addition theoremparticular to each
PDE , which would be algebraically tedious especially for vector-valued
PDE s such asStokes and elastostatics. This motivates the need for a scheme that can handle multiple
PDE s withoutcode changes. The present work fills this gap.
Overview and model problems . As with
QBX , we construct an approximate representation for
PDE solu-tions in a small region abutting the boundary, then use it for near and on-surface evaluations. However, incontrast to
QBX , our representation is an equivalent density on a closed curve enclosing this region; whendiscretized, this gives a ring of “proxy” point sources (also known as the method of fundamental solutions [ ] ). Matching is done at a second smaller ring of “check” points where a refined smooth quadrature isaccurate, thus the only dependence on the PDE is via point-to-point kernel evaluations — the method iskernel-independent, and essentially
PDE -independent.We focus on Dirichlet boundary-value problems (cid:76) u = Ω , (1.2) u = f on Γ , (1.3)where Ω is a simply-connected interior domain with smooth boundary Γ , for the following partial dif-ferential operators: (cid:76) u = ∆ u Laplace, ( ∆ − λ ) u Yukawa, ( ∆ + ω ) u Helmholtz ( Im ω ≥ ) , ∆ u − ∇ p Stokes ( subject to ∇· u = ) , ∆ u + − ν ∇ ∇· u Elastostatic. (1.4)
Abtin Rahimian et al.
To obtain well-conditioned formulations of the problem, we represent the solution of Eqs. (1.2–1.4) for x ∈ Ω by the double-layer potentials u ( x ) = (cid:68) [ φ ]( x ) : = (cid:90) Γ ∂ Φ ( x , y ) ∂ n y φ ( y ) d s ( y ) , (1.5)where Φ is the fundamental solution for the operator (cid:76) , and φ is an unknown density. The fundamentalsolutions for the operators listed in Eq. (1.4) are given in Appendix A. A standard step (see, e.g., [ ] )is now to substitute Eq. (1.5) into the boundary condition and use the jump relation for the potential toobtain the second-kind integral equation − φ ( x ) + ( D φ )( x ) = f ( x ) , for x ∈ Γ , (1.6)where D is the restriction of (cid:68) to the curve. Here, the integral implicit in the integral operator D mustbe taken in the principal value sense. Discretization and overall approach . In general, a smooth quadrature is a set of nodes x i ∈ Γ withassociated weights w i , such that (cid:90) Γ f d s ≈ N (cid:88) i = w i f ( x i ) , (1.7)holds to high accuracy for smooth functions on Γ — including the density φ . In this work, we use q -node Gauss–Legendre quadrature scheme on panels, and for convergence tests, we increase the numberof panels while holding q fixed. Upon discretization, Eq. (1.6) will be approximated by the linear system N (cid:88) j = A i j φ j = f ( x i ) , i =
1, . . . , N , (1.8)whose solution φ = { φ j } Nj = approximates the density values at the collocation points. In practice, forlarge problems, the matrix A is not constructed explicitly, but instead the matrix-vector product A φ isevaluated using the fast multipole method. We test the QBKIX scheme both for applying matrix A (i.e.,on-surface evaluation) and evaluating the solution at arbitrary points, near-evaluation in particular.The system matrix elements are computed using the Nyström method [
38, Ch. 12 ] . If the operator D is smooth on Γ × Γ , we use a smooth Nyström formula; e.g., for Laplace, A i j = ∂ Φ ( x i , x j ) ∂ n x j w j , i (cid:54) = j , − − κ ( x j ) π w j , i = j , (1.9)where κ ( x ) is the curvature at x ∈ Γ . This discretization achieves super-algebraic convergence. However,for Yukawa and Helmholtz in , and all elliptic kernels, singular quadrature is needed.In contrast to established approaches using specialized singular quadratures, we follow the ideaunderlying the QBX method: applying A to a vector φ is equivalent to evaluating the interior limit of thedouble-layer potential due to a smooth density interpolated from φ . This observation leads to the QBKIX idea: use a fast algorithm combined with the smooth quadrature scheme, Eq. (1.7), for point evaluation away from the surface — at points we refer to as check points — and interpolate from these points to theon surface point, to compute A φ for the Krylov iteration. As this interpolation can be done using pointson one or both sides of the surface, in Section 4.2 we compare “one-sided” and “two-sided” variants of QBKIX with respect to their spectra and iterative convergence rates. biquitous evaluation of layer potentials using
QBKIX Although we are focusing on interior Dirichlet tests and Nyström-style sampled representation of thedensity in this work,
QBKIX is applicable for Neumann or other boundary conditions, and Galerkin andother discretization types. Moreover, while the approach presented in this paper is restricted to , thereis no fundamental obstacle to an extension to .The rest of the paper is structured as follows. In Section 2 we present the QBKIX algorithm for in-tegration. We present an error analysis in Section 3. In Section 4, and report the results of numericalexperiments quantifying the accuracy of the method for a number of representative problems.
Given a closed curve Γ ⊂ (cid:82) with interior Ω , and Dirichlet data f on Γ , our goal is to numerically solvethe integral equation (1.6) for density and evaluate the solution of the underlying PDE at an arbitrarytarget point x ∈ Ω . We assume that Γ is parametrized by a 2 π -periodic piecewise-smooth function X ( t ) ,so that the arc length element is d s = | X (cid:48) ( t ) | d t , | X (cid:48) ( t ) | is bounded from below, and that X ( t ) and thedata function f ( t ) may be evaluated at any t ∈ [
0, 2 π ) . The boundary is subdivided into panels , whichcan be of different lengths, on which the native quadrature rule is defined (we use Gauss–Legendrequadrature), at q nodes x j per panel. We assume that the density is available as a vector of samples φ ( x j ) at the quadrature nodes.2.1 Single-point evaluationWe describe our method in the simplest form for computing the solution accurately at a given point x .We assume that there is a single point on Γ closest to x , on a panel of length L . We assume that at adistance 2 δ along the normal to the panel at any point, the native quadrature meets the target accuracyof evaluation, so the distance from x to the surface is less than 2 δ . We discuss how δ is chosen and howto ensure that this condition holds after the algorithm formulation.The local geometric configuration of various types of points we are using in our algorithm is shownin Figure 2.2. The setup shown in the image is for computing the potential accurately for any point x inside a disk B c δ of radius δ centered at c , touching the surface at a point x on a panel of length L .The points we use in the algorithm are placed on two concentric circles with the same center as theevaluation disk B c δ . The proxy points on a circle ∂ B c R of a radius R > δ , where we compute equivalentdensity values, are used to approximate the solution inside B c δ . The check points z i are on a circle ∂ B c r c of a radius r c < δ . At these points, we evaluate the solution accurately by using a smooth quadrature onpanels refined by a factor β . The check points are used to compute the equivalent density values at theproxy points as described below.The algorithm depends on a number of parameters; these parameters need to be chosen appropri-ately to achieve an overall target accuracy. Specific choices are discussed in the next section. The keysteps in the algorithm are(1) Set-up of proxy and check points.
We choose a center c ∈ Ω at a distance δ from Γ , such that x is no further from c than δ . E.g., for x ∈ Γ , we set c = x − δ n , where n is the outward normal. n p proxy points y j are arranged equally on the circle of radius R with center c , where R > δ is of order L . Similarly n c check points z i are arranged on the concentric circle of radius r c < δ (Fig. 2.1).(2) Upsampling the density.
Each panel is split into β panels corresponding to equal ranges of t , togive a set of β N fine-scale nodes ˜ x l with weights ˜ w l . The global factor β is chosen so that thesolution can be evaluated accurately at the check points, i.e., at a distance δ − r c from the surface.The density is interpolated from its original samples φ ( x j ) on each panel, using q th order Lagrangeinterpolation to the fine-scale nodes, to give the refined vector of samples ˜ φ l , l =
1, . . . , β N . Abtin Rahimian et al. x δ B c δ checkcircle ∂ B c r c r c z i R proxycircle ∂ B c R y j ρ + " singularityquadraturepoints L Fig. 2.1: S CHEMATIC OF A KERNEL - INDEPENDENT EXPANSION . Geometry of
QBKIX , with proxy and check circles centered at c neara panel of length L of the boundary Γ discretized with q Gauss–Legendre sample points. The evaluation domain B c δ is a disccentered at c of radius δ (dashed circle abutting the boundary at x ). The points z i are the check points on the circle ∂ B c r c ofradius r c , and y j are the proxy points on the circle ∂ B c R of radius R. For error analysis, the singularities of the exact solutionare assumed to be at a distance farther than ρ from c . Note that, for clarity, the relative sizes of circles and distances betweensamples are different from the ones actually used. (3) Direct upsampled evaluation at check points.
The integral is evaluated at each check point z i using the fine-scale boundary native quadrature:˜ u ( z i ) = β N (cid:88) l = ∂ Φ ( z i , x l ) ∂ n x l ˜ φ l ˜ w l . (2.1)Denote by ˜ u : = { ˜ u ( z i ) } n c i = the column vector of these values at the check points.(4) Solving for the equivalent density values.
Next, we construct an n c × n p matrix Q with elements Q i j = Φ ( z i , y j ) . (2.2)Applying Q to a vector of density values at proxy points computes a periodic trapezoidal rule ap-proximation to the single-layer potential corresponding to this density evaluated at check points.Then we solve a small, dense, and ill-conditioned linear system Q α = ˜ u , (2.3)in the least-squares sense, to get the set of proxy density values α : = { α j } n p j = . The ill-conditioningarises from the exponential decay of singular values in the single-layer operator between concentriccircles (see Fig. 3.1). Despite this, if Eq. (2.3) is solved in a backward-stable manner, a high-accuracyresult is obtained (cf. [ ] , we explain the details below for completeness).(5) Evaluation of the proxy sources at the target.
Finally, the equivalent density is evaluated at thetarget x , ˆ u ( x ) = n p (cid:88) j = α j Φ ( x , y j ) , (2.4)We may view this as an approximation for the true solution u in the basis of fundamental solutionscentered at the proxy points, that holds to high accuracy in the disk B c δ .Figure 2.2 illustrates the stages of QBKIX evaluation for a set of target points lying in a single disk B c δ .The final evaluation of Eq. (2.4) over the disc of target points has around 12 digits of accuracy. biquitous evaluation of layer potentials using QBKIX Direct evaluation with the boundary consisting of 26 patches (a)
Error using the nativequadrature, stage (1).
Direct evaluation with the boundary consisting of 26 patches (b)
Error using the refined (native)quadrature, log | ˜ u − u | , stages(2) & (3). max ( Error ) = e + QBKIX error (the boundary consisting of 26 patches and kernels LaplaceDL and LaplaceDL) (c)
Error of the expansion withinB c δ , log | ˆ u − u | . Stages (4) & (5). − − − − − − − − Fig. 2.2: S TAGES OF
QBKIX
CONSTRUCTION . The stages given in Section 2 are illustrated using plots of the log of the evaluationerror near the boundary, for the double-layer density φ ≡ for Laplace’s equation. The evaluation disc B c δ (dashed circle), checkcircle ∂ B c r c (solid circle) are shown, and proxy points are not shown. Handling the ill-conditioned linear solves . The ill-conditioned system Eq. (2.3) is solved by applying aregularized pseudo-inverse, as follows. Let (cid:34) pinv be the desired relative accuracy for inversion; typicallywe set (cid:34) pinv = − . Then, taking the singular value decomposition ( SVD ) [ ] Q = U Σ V ∗ with Σ = diag { σ j } being the diagonal matrix of singular values, we write Σ † : = diag { σ † j } where σ † j = (cid:168) σ − j , σ j > (cid:34) pinv σ ,0, otherwise. (2.5)Then we use the solution α : = V ( Σ † U ∗ u ) . (2.6)Note that the matrices U ∗ and V must be applied in two separate steps (as indicated by the parenthesis)for backward stability [ ] , since a matrix-vector multiply with the single pseudo-inverse matrix Q † : = V Σ † U ∗ is unstable due to round-off error caused by its large entries. If k is the number of singular valuesgreater than (cid:34) pinv , i.e., the numerical (cid:34) pinv -rank of the matrix Q , the factors V and U ∗ have sizes n p × k and k × n c respectively. Parameter summary . The algorithm described above uses a number of parameters, which we summa-rize here.The following parameters are defined globally: – The quadrature order q , which determines the number of samples per panel, and both far-fieldevaluation accuracy and, together with β , the accuracy of evaluation at check points. This parameteris selected arbitrarily based on the desired overall accuracy. We use q =
16, which is sufficient for fulldouble precision of integration in the far field. – The panel refinement factor β which needs to be chosen to maintain desired accuracy for check pointevaluation. – The numbers of proxy points n p and check points n c ; the former determines how accurate the ap-proximation inside B c δ can be and the latter is chosen to have enough sampling.Three additional parameters, the accurate evaluation distance δ , the proxy point circle radius R andthe check point circle radius r c , are panel-dependent, and are chosen with respect to panel size L . Acareful choice of all of these, as fractions of L , is needed to achieve a target error without requiringexcessive refinement. We discuss the choice of these parameters in Section 3. Abtin Rahimian et al.
Defining panels . In our experiments, we consider two ways of defining panels. The first approach isprimarily needed to understand the convergence of the method with respect to the number of panels,i.e., for a given number of panels, we determine the error. In this case, we simply partition the parametricdomain of X ( t ) into M equal-sized intervals, with one panel corresponding to each interval. We assumethe parametrization to be sufficiently close to an arclength parametrization, so that the panel length haslittle variation, and choose M to be fine enough so that the geometric condition on the check points issatisfied.In a more practical scenario, when a target error is specified, we need to determine panel sizesadaptively. The key requirement that needs to be satisfied by panels is that the accuracy of check-pointevaluation at stage 2 matches the target accuracy in the far field (i.e., points farther than 2 δ from theboundary). The adaptive refinement starts with one panel covering the entire boundary, then recursivelysplitting panels into two equal pieces in parameter t , until all panels are deemed admissible or theirlength is less than a set tolerance (cid:34) .A panel is admissible if (i) the interpolation of X ( t ) and f ( t ) from a q -node panel at the collocationpoints of the two q -nodes Gauss–Legendre panels (obtained by splitting the coarse panel to two pieces)matches the direct evaluation of X and f on the finer nodes, to a maximum absolute tolerance (cid:34) a , whichwe choose as 10 − unless stated otherwise; (ii) it is no more than twice the parameter length of thatof its neighbors; (iii) the length of the panel does not exceed a given fraction of the minimal radius ofcurvature at a point of the panel, or is less than a minimal length proportional to the target error; and(iv) any check point corresponding to a point x is not closer than δ − r c to any point on the surface.The second criterion ensures that the panels are the leaves of a balanced binary tree, which is neededfor accurate evaluation of integrals at the check points. For domains with sharp corners, the forth andsecond conditions imply dyadic refinement of panel length bounded below by panel minimum length (cid:34) l .In both cases, the result is a set of N nodes x j = X ( t j ) , where t j are the parameter values of thenodes, with weights w j = | X (cid:48) ( t j ) | w (cid:48) j where w (cid:48) j are the Gauss–Legendre weights scaled by the panelparametric lengths. This native quadrature approximates the boundary condition f with target accuracy (cid:34) a . It follows from Eq. (1.6) that this also holds for the density φ , as φ to be no less smooth than f and X .2.2 On-surface evaluation for iterative solution of the linear systemAs discussed in the introduction, one context where singular quadratures are needed is for applying A ,the matrix discretization of the operator ( − I + D ) , to the current density vector φ during the iterativesolution of Eq. (1.8). This matrix-vector multiplication is equivalent to evaluation of the interior limit ofthe double-layer potential at the nodes due to the smooth interpolant of the density vector. As with QBX [
34, Sec. 3.5 ] , one may exploit this in two different ways. – One-sided
QBKIX : as stated above, we use the interior limit of the potential at the nodes for A φ . – Two-sided
QBKIX : we average the interior and exterior limits of the potential at the nodes, which,by canceling the jump relation terms, applies a matrix approximation to the operator D . We thenexplicitly add − φ to the answer.Although mathematically equivalent, these two variants smooth high-frequency components in the den-sity differently: one-sided QBKIX tends to dampen these components, leading to an accumulation ofeigenvalues of A around zero. This has a negative impact on convergence. In contrast, for two-sided QBKIX , since the approximation of D tends to damp high-frequency components, the explicit inclusionof − I ensures that these components end up being multiplied by a number very close to − , whichleads to better clustering of the spectrum and improved convergence rates. We present a numericalcomparison of these two alternatives in Section 4.2. biquitous evaluation of layer potentials using QBKIX x , the brute-force approach is to run the algorithm described above,including construction of check and proxy points, for each sample point separately. This is highly ineffi-cient, and the following obvious optimizations can be applied: – The upsampled density on the fine-scale nodes need be computed only once, and each expansioncenter may be chosen to cover several targets; this requires increasing evaluation disk radius δ ,adjusting other parameters accordingly. – The
SVD of matrices Q may be precomputed. For translation- and scale-invariant kernels, (i.e., allkernels we consider except Yukawa and Helmholtz) these matrices do not depend on the choice ofthe center and circle radii, as long as the ratio R / r c is fixed. – One may use the kernel-independent FMM method for evaluation of the solution at the check pointsfor all target points at once.We consider the complexity of using
QBKIX for the task of on-surface evaluation at all boundary nodes x ∈ Γ . For a boundary with M panels and q -node Gauss–Legendre quadrature on each, there are N = Mq nodes in total. We use a conservative assumption that a distinct set of check and proxy points is usedfor each of the targets. Then, using KIFMM , the evaluation of the boundary integral from the β -refinedboundary to the check points is (cid:79) (cid:0) ( β + n c ) N (cid:1) . We assume that the factorization of the pseudo-inversefor computing the equivalent densities α is precomputed. The cost of applying the factors V and U ∗ , ofsizes n p × k and k × n c , for targets point is (cid:79) (cid:128) k ( n c + n p ) N (cid:138) . The cost of evaluation of the approximationfrom proxy density values at target points is (cid:79) (cid:128) N n p (cid:138) .We conclude that the overall cost is (cid:79) (cid:128) ( β + n c + kn c + kn p + n p ) N (cid:138) , which for typical choices β = n c = n p reduces to (cid:79) (cid:128) kn p N (cid:138) . We see that the scheme is linear in N , but with a prefactor of order k (since, as discussed in the next section, n p is of order k ). The two-sided variant involves anotheroverall factor of 2.If the same check and proxy points are used for a number of targets, an additional, potentiallyvery large, constant-factor speedup can be obtained. The speedup factor is proportional to the averagenumber of targets handled by each set of check and proxy points. In this section, we present theoretical results, focusing on the cases of scalar u governed by the Laplaceequation ∆ u = ( ∆ + ω ) u = ω . We expect similar resultsfor other elliptic PDE s in Eq. (1.4).We split
QBKIX into two stages: (i) evaluation of u on the check points using a refined native quadra-ture, with the associated error e c ; (ii) solution of a small linear system to determine the equivalentdensity values α at the proxy points that best represent u at the check points. This is followed by evalu-ating the approximation of u at target points using these density values.At the first stage, the error e c is effectively the smooth quadrature error of the refined panels. Theprimary focus of our analysis is on the second stage. We analyze the error behavior in the idealizedsituation of exact arithmetic and infinitely many check points, obtaining the dependence of the second-stage error e on δ , R , ρ , and n p . We then describe a heuristic model for the effects of finite-precisioncomputations, which adds an extra term to e , depending on e c , δ , r c , and k .We use the overall error model, along with experiments, to provide a choice of the various parametersin the scheme resulting in the on- and near-surface evaluation errors of the same magnitude as the far-field integration errors. u on the check points is done by approximating the exact integral Eq. (1.5)by Eq. (2.1) using q -node Gauss–Legendre quadrature on panels (subdivided by factor β ). For a flatpanel, the error e c in this evaluation is bounded by standard quadrature estimates giving a term of theform C q ( L / ( β d )) q (cid:107) φ (cid:107) C q where d = δ − r c is the closest distance of check points to the panel, and φ denotes the density for which we evaluate the integrals. Our adaptive refinement procedure ensuresthat the formula still holds, as the radius of curvature of the panel is larger than its length, and hencelarger than δ .This estimate has the form of the second term in [
34, Theorem 1 ] , and for convergence as the panellength L going to zero, it requires L / d to converge to zero as well. Instead of following this route, we fixthe ratio L / d to a constant, by choosing δ and r c as fractions of L . If L / ( β d ) is sufficiently small, a high-order quadrature for sufficiently large q allows us to compute the integrals with any desired precision.For instance, when q =
16, it is sufficient to use L / ( β d ) = /
2, to obtain an error on the order of 10 − at distance d from the panel.3.2 Error of the proxy point representation in exact arithmeticNext, we analyze the dependence of the error (computed in exact arithmetic) of the second stage of QBKIX on the number of proxy points n p , the proxy circle radius R , and the distance r from the center c to the evaluation point. The distance r could be either smaller than δ if targets are away from thesurface, equal to δ if B c δ touches the surface at a single point, or exceed δ if there are several on-surfacetargets in B c δ ; we focus our attention to the case where r ≤ δ .Let ˆ u be given by the proxy representation, Eq. (2.4), with equivalent density values α j at proxypoints y j , j =
1, . . . , n p . We consider evaluation of the approximation ˆ u in B c r , the disc of radius r centered at c , given correct values for u at a very large number of check points n c , so that we can replacethe discrete least-squares problem we solve with a continuous one.Let the equivalent densities α j be chosen to minimize the L error on the check circle, i.e., α = arg min α ∈ (cid:67) np (cid:107) ˆ u − u (cid:107) L ( ∂ B c rc ) . (3.1)By convergence of the periodic trapezoidal quadrature on the check points, this corresponds to the n c → ∞ limit of the QBKIX scheme. Let e ( r ) : = sup x ∈ Ω ∩ B c r | ˆ u ( x ) − u ( x ) | , (3.2)be the upper bound on the pointwise error in the part of the disc lying inside the closure of the domain.We have the following bounds on e when u is sufficiently regular, meaning that any singularities in thecontinuation of u is further than some distance ρ > δ from the center of the expansion c . Theorem 3.1
Let u be continuable as a regular solution to the Laplace or Helmholtz equation in the closeddisc of radius ρ centered at c . Let R > δ in the Laplace case. Let the QBKIX equivalent density values at proxypoints be solved in exact arithmetic in the least-squares sense on the check circle as in Eq. (3.1) , and let e bedefined by Eq. (3.2) where ˆ u is the expansion in Eq. (2.4) . Then, in a disc of radius re ( r ) ≤ C (cid:0) r ρ (cid:1) n p / , ρ r < R , C (cid:112) n p (cid:0) rR (cid:1) n p , ρ r = R , C (cid:0) rR (cid:1) n p , ρ r > R , (3.3) where in each case, C indicates a constant that may depend on u (and ω in the Helmholtz case), c , r, andR but not on n p . biquitous evaluation of layer potentials using QBKIX Proof
Following the technique of Barnett and Betcke [
4, Theorem 3 ] , we only need to show that thereexists some choice of density values α for which the estimate holds; the least-squares solution cannot beworse than this. We choose density values α to cancel the Fourier coefficients with frequency | n | < n p / ˆ u − u on the check circle.By uniqueness of the local expansion for the regular PDE solution (in polar coordinates, (cid:80) n ≥ a n r n e in θ for Laplace or (cid:80) n ∈ (cid:90) a n J n ( ω r ) e in θ for Helmholtz) this choice of density values also cancels the sameFourier coefficients on any circle centered at c with radius less than R . Applying [
4, Theorem 3 ] for theHelmholtz case, the L -norm of the error on the circle of radius r obeys a bound of the form Eq. (3.3).Barnett and Betcke [
4, Section 2.1 ] produce the Laplace case as a limit of the Helmholtz case; however,one also needs the result that the constant single-layer density generates the constant potential log R / r c ,which excludes R = r c because it can only produce zero-mean data on the circle.Finally, we need to show that the sup norm of the error on the circle of radius r is bounded by the L -norm; this holds since the error ˆ u − u is a regular PDE solution in a disc with radius strictly larger than r , namely B c min ( R , ρ ) . Thus, its Fourier coefficients on the r -circle decay exponentially in | n | , and are thussummable with a bound controlled by the L norm. In the case where B c r lies partially outside Ω , onemay continue u as a regular PDE solution in the disc and apply the above. (cid:117)(cid:116)
Remark 3.1
The above derivation relies on analysis from the literature on the method of fundamentalsolutions (
MFS ). The original result for the Laplace equation is due to Katsurada [
32, Theorem 2.2 ] ,which considers the case n c = n p and restricted to r = r c . We extend this result to include extrapolation from the check radius r c out to larger radii r .Remarkably, r c does not appear in Eq. (3.3), because in exact arithmetic it does not matter at whatradius the Fourier coefficients are matched. In the next section we will see that in practice roundingerror strongly affects the choice of r c since the extrapolation is ill-conditioned.A surprising aspect of Theorem 3.1 is that u may have singularities closer to the center than the proxyradius R and yet exponential convergence still holds; this is closely related to the Runge approximationtheorem. Remark 3.2
The two regimes in Eq. (3.3) may be interpreted as follows: • r < R ρ : the solution u is relatively rough (has a nearby singularity), and error is controlled by thedecay of the local expansion coefficients a n of u for orders beyond n p / • r > R ρ : the solution u is smooth, and error is controlled instead by aliasing (in Fourier coefficientspace) due to the discreteness of the proxy point representation on the proxy circle.We observe in numerical experiments that when the boundary is adaptively refined based on the bound-ary data as in Section 2, L ≈ ρ and the expansion centers that dominate the error in a domain aretypically those that are near to a singularity of the solution. Such centers are typically in the roughregime.Note that the boundary Γ may intersect the closed disc, and still u may be continued as a PDE solutioninto the closed disc. This requires the boundary data f or density to be analytic — see [ ] for relatedanalysis of QBX in this case.
Remark 3.3 (Extension of analysis to other kernels)
It is clearly of interest to have a kernel-independentextension of Theorem 3.1 that would apply also to vector
PDE s such as Stokes. Initial attempts suggestthis requires significantly more complicated analysis, since to use the method of the above proof oneneeds to be able to write down a proxy coefficient vector α that produces a single Fourier mode on thecheck circle plus exponentially decaying amounts of aliased modes, which is challenging even in theStokes case. We leave this for future work. r c in Theorem 3.1 relies on exact arithmetic; since the extrapolation from r c to alarger r is ill-conditioned. Moreover, due to finite precision, there are possibly fewer than n p functionsavailable to cancel the Fourier coefficients. As a result, we need to study the effect of rounding error on ˆ u − u . Rather than attempting a rigorous analysis, we present a heuristic model and demonstrate that itagrees well with numerical observations.We first show that the n th singular value of the matrix Q in Eq. (2.2) decays as n ( r c / R ) n / , i.e.,marginally faster than exponentially. In the continuous limit ( n p , n c → ∞ ), this corresponds to the decayof the eigenvalues of the single-layer operator with kernel Φ , whose eigenfunctions are the Fouriermodes, since the operator is convolutional. For the Laplace equation, the potential defined in polarcoordinates centered at c as v ( r , θ ) = (cid:40) ( R / n )( r / R ) n e in θ , r ≤ R , ( R / n )( r / R ) − n e in θ , otherwise ,solves the PDE everywhere except at r = R , where the jump in radial derivative is e in θ . We conclude that v is the single-layer potential due to the n th Fourier mode density. Substituting r = r c , and recalling thatthe n th singular value is eigenvalue for the frequency n /
2, as the frequencies are in the range − n / n /
2, we conclude that σ n = n ( r c / R ) n / .The above argument also applies for the Stokes case except due to having two vector components, n th singular value of matrix Q corresponds to the eigenvalue for frequency n /
4. The Helmholtz case —although there are (cid:79) ( ω ) eigenvalues that do not decay — is asymptotically identical to Laplace [ ] . To verify this asymptotic behavior, in Fig. 3.1 we show the decay of singular values forseveral kernels.When the pseudoinverse of Q is computed based on Eq. (2.5), only k singular values lying above (cid:34) pinv σ are retained. The corresponding singular vectors approximate the lowest Fourier modes up tofrequency | n | < k / PDE cases). Thus, equating up to constants the k th singular valueabove to (cid:34) pinv , the ranks of the matrices in the pseudoinverse are k ≈ min (cid:128) k m , n p (cid:138) , k m = ( /(cid:34) pinv ) log ( R / r c ) , (3.4)and the highest (Nyquist) frequency they can represent is k / − − Index ( n ) S i n g u l a r V a l u e s n ) 0 100 200Index ( n ) R = R = R = R = R = R = ( T ) R = ( T ) R = ( T ) R = ( T ) R = ( T ) Laplace Helmholtz Stokes
Fig. 3.1: S INGULAR VALUES OF PROXY TO CHECK MATRIX . The solid lines are the singular values of Q for different R and differentsingle-layer kernels, and the dashed lines labeled ( T ) are the theoretical decay: n ( r c / R ) n / for Laplace or Helmholtz, and n ( r c / R ) n / for Stokes, where n denotes the index of the singular value. Other parameters are r c = , n p = , n c = .For the Helmholtz problem, the dashed lines show the asymptotic bound for the singular values and are not accurate for smallindices; the interested reader is referred to [
4, Eq. (14) ] . biquitous evaluation of layer potentials using QBKIX The values of ˜ u at the check points have error bounded by e c , so in this model we expect the errorsto be amplified (by considering the local expansion as above) to become e c ( r / r c ) k / at the evaluationradius r .3.4 Error bounds and optimal parameter choicesCombining the results from Sections 3.2 and 3.3 for a kernel-independent expansion, using n p proxypoints, the error is bounded by e ( r ) ≤ C (cid:18) r ρ (cid:19) k / + C e c (cid:18) rr c (cid:19) k / , ρ r < R , C (cid:129) rR (cid:139) n p + C e c (cid:18) rr c (cid:19) k / , ρ r > R , (3.5)where C represents possibly different constants in each case (omitting the case ρ r = R ).In Fig. 3.2, we show how this formula models the error growth for a single kernel-independentexpansion interpolating a Laplace solution in free space with a known nearest singularity at variousdistances ρ , for a typical choice of ratio R / r c =
8. The key observation is that, despite its simplicity, ourmodel Eq. (3.5) explains well the observed error behavior. Other salient features of the plots include: – As r increases beyond r c , errors grow rapidly dominated by the second term in the error estimate. – The error is mostly controlled by k and increasing n p beyond k m ≈
27 (defined in Eq. (3.4)) has notangible effect unless ρ r > R (i.e., right half of left plot).Figure 3.3 instead continuously varies R /ρ (the inverse scaled singularity distance), showing the sameeffect: a relatively distant singularity allows high accuracy expansion out to larger r / R . Choice of parameters . Using the model Eq. (3.5), one can make choices for R , r c , δ , and n p to achievea desired accuracy (cid:34) . An unknown in applying this in a practical setting is the singularity distance ρ .However, in any high-accuracy choice of boundary quadrature, such as the adaptive panel quadrature ofSection 2, panels are refined such that the data f and hence the density φ and the solution u are smooth − − r c R /ρ rR E rr o r r c rR r c ρ/ R rR n p = n p = n p = n p = n p = n p = ( T ) n p = ( T ) n p = ( T ) n p = ( T ) n p = ( T ) Fig. 3.2: E RROR BOUNDS FOR L APLACE
QBKIX
WITH KNOWN SINGULARITY . Errors e observed (solid lines) and predicted by Eq. (3.5) (dashed lines) for a single expansion with different singularity distances ρ = R , R , and R, and different numbers of proxypoints n p . The expansion is centered at c = [
0, 0 ] and the solution u ( x ) = − log | x − x | , x = ρ e i / is a harmonic functionwith a singularity at distance ρ . Laplace single-layer kernel is used for the expansion. The error is the maximum error over theB c r as defined in Eq. (3.2) . The proxy to check radius ratio is R / r c = , the number of checks is set to n c = n p , e c = − , andk m ≈ (given by Eq. (3.4) with (cid:34) pinv = − ). The constants C in Eq. (3.5) were chosen to qualitatively match the trend lines(all set to ). R / r r / R - . - . - . - . - . - . - . - . Fig. 3.3: E RROR AT DIFFERENT EVALUATION RADII . The error for evaluation of a single expansion with various R and r, but fixedr c = ρ/ and ρ . The expansion is interpolating a harmonic function (similar to the one used in Fig. 3.2) with singularity atdistance ρ = , using the Laplace double-layer kernel. The dotted lines are r = mr c for m =
1, 2, and . In practice, we have nodirect control on R ρ , and it is implied by the panel size. Here we chose n p = , and n c = n p ; the trends are the same for lowern p and n c . on the local panel scale L , thus we expect singularities to be at least of order L distant from the center.Indeed, we experimentally observe (in tests where we know the location of singularity, e.g., Fig. 3.4or Section 4.3) that when the panels are adaptively refined, L < ρ , and consequently the convergencebehavior is most like the left-hand plot of Fig. 3.2.Given the target accuracy of (cid:34) for the solution and the selected native quadrature order q , the adap-tive refinement of boundary sets the panel length L . We use the following steps to glean the value ofother parameters. Since the constants in the error estimates are problem dependent and unknown, weset them to unity. To have a concrete example, we pick (cid:34) = − and q = Setting δ : By construction, points farther than 2 δ from the boundary are evaluated using thenative quadrature. To meet the desired error (cid:34) at these points, L δ ≈ (cid:34) / q , which implies δ ≈ L / (cid:34) = − , q = Setting k m , R / r c , and n p : Requiring that the two terms in the error estimate (i.e., proxy pointrepresentation and extrapolation errors) have similar contribution at the on surface point ( r = δ )and assuming that L ≈ ρ we can estimate the minimum required k based on the proxy represen-tation error in the rough regime: (cid:18) δρ (cid:19) k / ≈ (cid:34) or k ≈ (cid:34) log ( δ/ L ) , (3.6)implying k ≈
32 for L /δ = (cid:34) = − . Since k is bounded by min ( k m , n p ) , knowing minimum k implies a lower bound for k m and n p . Therefore, reorganizing Eq. (3.4), we have R / r c = (cid:34) / k pinv ≈ (cid:34) pinv = − .(3) Setting r c /δ and β : Inspecting the extrapolation error at an on surface point, we have e e ( δ ) ≈ e c (cid:18) δ r c (cid:19) k / ≈ (cid:18) L β ( δ − r c ) (cid:19) q (cid:18) δ r c (cid:19) k / ≈ (cid:18) L βδ (cid:19) q ( − θ ) q θ k / , (3.7)where θ = r c /δ . This expression attains its minimum at θ = k q + k . For q =
16 and k =
32, wehave θ = /
3. As we require that two terms in the error estimate have similar contribution, weuse e e ( δ ) and estimate β : β ≈ L / δ ( − θ ) θ k / q (cid:34) / q , (3.8)implying β =
5, for the choices of parameter listed above. biquitous evaluation of layer potentials using
QBKIX
150 0.2 0.410 − − − − δ/ L e δ/ L ˜ ρ = ρ = ρ = ρ = ρ = Fig. 3.4: E RROR VS . CENTER AND SINGULARITY DISTANCES . The induced error for singularities and centers at various distances fromthe boundary for the Laplace Dirichlet interior
BVP , in the domain shown in Fig. 4.2. The boundary data is generated by puttinga Laplace singularity at distance ˜ ρ from the boundary — the singularity distance to the center of expansion is ρ ≥ ˜ ρ + δ . Thedensity is solved directly and QBKIX is used only for evaluation. The error is computed using the known solution correspondingto the boundary data. The left plot shows the errors for the case with fixed number of panels on the boundary (M = panels).In this plot, because L is fixed, L / ˜ ρ is decreasing by increasing ˜ ρ . The right plot shows the errors for adaptive refinement of theboundary with (cid:34) a = − . Here, since L is chosen adaptively due to the boundary data, it increases as the solution becomessmoother. Because, L is chosen proportional to ˜ ρ , the error curves almost collapse to one. We use n p = , n c = n p , r c = δ/ and R = r c . In both cases, the center of expansion is located based on the panel size at distance δ . Note that we have not analyzed the effect of finite n c , but find that the choice n c = n p behavesindistinguishably from the limit n c → ∞ ; we attribute this to the rapid convergence of the periodictrapezoid rule on the check points. In this section, we present the results of numerical tests demonstrating the accuracy and versatility of the
QBKIX algorithm for on-surface evaluation needed for the boundary integral equation solver and solutionevaluation close to the surface. In the following experiments, unless noted otherwise, we use
QBKIX forboth tasks.4.1 Convergence with respect to the number of panelsIn Table 4.1, we report the convergence of the solution evaluated at the interior points using non-adaptive boundary quadrature with increasing number of panels. The test solution is the potential dueto a set of singularities at the source points shown outside the domain. These source points are used togenerate the boundary data f and the reference solution to check the error. For all problems, the double-layer formulation is used, except for the Helmholtz for which a combined-field formulation u = (cid:68) [ φ ] + i ω (cid:83) [ φ ] , where (cid:83) is the single-layer potential [
14, Section 3.2 ] , is used. This representation addressesproblems associated with resonance of the complementary domain. The double-layer (or combined-field)density φ is solved using QBKIX to evaluate the matrix-vector product in each iteration of
GMRES . Theerror in the density is quantified by computing the solution from φ , Eq. (1.5), at a set of target points inthe interior of the domain. For the first three kernels, which are smooth, we also report the convergenceusing the Nyström ( direct ) evaluation, Eq. (1.9), which by comparison against one- or two-sided QBKIX shows how much of the error is due to
QBKIX .In all cases, it can be seen that
QBKIX gives high-order convergence rate that is independent of thetype of the kernel. We notice that the error performance of the two-sided variant is worse than one-
Geometry Kernel Quadrature Absolute error (Number of panels)Laplace Direct 2.90 e − ( ) e − ( ) e − ( ) e − ( ) QBKIX (one) 3.39 e − ( ) e − ( ) e − ( ) e − ( ) QBKIX (two) 2.25 e − ( ) e − ( ) e − ( ) e − ( ) Laplace Direct 5.80 e − ( ) e − ( ) e − ( ) e − ( ) QBKIX (one) 2.49 e + ( ) e − ( ) e − ( ) e − ( ) QBKIX (two) 4.29 e − ( ) e − ( ) e − ( ) e − ( ) Stokes Direct 1.48 e − ( ) e − ( ) e − ( ) e − ( ) QBKIX (one) 2.89 e − ( ) e − ( ) e − ( ) e − ( ) QBKIX (two) 6.95 e − ( ) e − ( ) e − ( ) e − ( ) Helmholtz ( ω = QBKIX (one) 2.12 e − ( ) e − ( ) e − ( ) e − ( ) QBKIX (two) 3.97 e − ( ) e − ( ) e − ( ) e − ( ) Yukawa ( λ = QBKIX (one) 1.60 e − ( ) e − ( ) e − ( ) e − ( ) QBKIX (two) 5.44 e − ( ) e − ( ) e − ( ) e − ( ) Elastostatic ( ν = QBKIX (one) 2.07 e − ( ) e − ( ) e − ( ) e − ( ) QBKIX (two) 3.17 e − ( ) e − ( ) e − ( ) e − ( ) Table 4.1: S OLUTION CONVERGENCE VS . NUMBER OF PANELS . Error in the solution to interior Dirichlet boundary value problemsusing non-adaptive M-panel quadrature and
QBKIX for solution. The subplots show Γ (solid) and the exterior sources used togenerate the solution, and interior test points. There are 40 source points outside the domain and error is measured on 40 pointsinside. The error is the maximum of absolute error over these interior points. The numerical parameters are n p = , n c = n p ,R = r c , and δ = r c . “Direct” indicates usage of the quadrature of Eq. (1.9) instead of QBKIX for the linear solve. “One” and“two” indicate one- or two-sided versions of on-surface
QBKIX discussed in Section 2.2. sided at the same number of panels (however, as we discuss below, it is valuable since it improves theconvergence rate of
GMRES ).4.2 Operator spectrum and
GMRES convergence rateWe now perform numerical tests of the one-sided and two-sided variants of on-surface evaluation of
QBKIX discussed in Section 2.2 and compare it to direct use of an accurate quadrature. To simplifycomparisons, we use an operator with a smooth kernel (Laplace). The spectra and convergence behaviorfor singular kernels is similar. In Fig. 4.1 we plot — for the domain shown in Fig. 4.2 and the Laplaceequation — the eigenvalues for four different approximations to the operator − + D : one-sided (interior) QBKIX , the one-sided (exterior)
QBKIX , two-sided
QBKIX , and the quadrature given by Eq. (1.9), to whichwe refer as direct . The exterior version of
QBKIX is constructed similarly to the interior variant discussedin Section 2. The only modification is that for each collocation point x on Γ , we place an expansioncenter at c = x + δ n . We see that the one-sided variants have clusters of eigenvalues near zero, whereasthe two-sided variant and the Nyström matrix have a cleaner spectrum with eigenvalue clustering onlyaround .A broader spread of the eigenvalues has a negative impact on GMRES convergence [ ] . Fig. 4.1,right, shows GMRES residual versus the iteration number for the interior, two-sided, and direct operatorswith two different right-hand sides (boundary data corresponding to a harmonic function and a randomright-hand side).The convergence of one-sided interior
QBKIX is identical to the Nyström method convergence up tothe residual magnitude on the order of numerical accuracy of
QBKIX , but it slows down once the residualdecreases below this value (near 10 − ). The two-sided variant has identical convergence behavior to thedirect method, and converges in a few iterations. We also show the residual for a random-right handside to expose the effect of near-zero eigenvalues: we see that convergence is very slow for the one-sidedscheme in this case, but for the two-sided scheme it is the same as for the true smooth data f . biquitous evaluation of layer potentials using QBKIX Index | λ | | Re ( λ ) | − − − − Iteration R e s i d u a l QBKIX interior
QBKIX two-sided Nyström
QBKIX exterior
QBKIX interior (random)
QBKIX two-sided (random) Nyström (random)interiorexteriortwo-sidedNyström
Fig. 4.1: T HE SPECTRA OF DISCRETIZATIONS OF THE L APLACE DOUBLE - LAYER OPERATOR . This figure shows eigenvalues, and the
GMRES convergence rate, for different discretizations of the Laplace double-layer operator in the domain shown in Fig. 4.2. Theleft plots show the real part and the magnitude of the eigenvalues corresponding the one-sided interior
QBKIX , one-sided exterior
QBKIX , two-sided
QBKIX , and the plain Nyström matrix. See Sections 2.2 and 4.2. The right plot shows the residual versus theiteration number for the three interior variants with two different right hand sides (boundary data corresponding to a harmonicfunction or random data). The residual of the two-sided and Nyström schemes are indistinguishable.
PDE sFor this set of tests, we use adaptive refinement as described in Section 2. We use
QBKIX both as theon-surface quadrature scheme when solving for the desired density as well as the evaluator for thenear-singular integrals. As before, we use boundary data sampled from a sum of fundamental solutionscentered at a set of points close to the boundary. Fig. 4.2 plots the error across the domain for all ofthe
PDE s listed in Eq. (1.4), on the points lying on a 600 ×
600 grid and interior to the domain. Whenan evaluation point is within 2 δ distance from the boundary, it is evaluated using the nearest QBKIX expansion. The remaining points are evaluated using Eq. (1.7) applied to Eq. (1.5).We observe that parameter choices which were selected for the Laplace equation perform well forthe other
PDE s. As expected, the highest error is due to expansions for panels adjacent to larger ones(e.g. Fig. 4.2(a)).4.4 Domain with a large number of cornersAs a final example, we use
QBKIX in a domain with 256 corners as shown in Fig. 4.3. A Laplace boundaryvalue problem is solved using
GMRES with tolerance for relative residual set to (cid:34) r = − . The boundarycondition is generated similar to the examples in Section 4.3 by placing 32 source points on a circle withradius 0.75 centered at [ ] (the domain’s bounding box is [
0, 1 ] × [
0, 1 ] ).The boundary of the domain is adaptively refined, with minimum panel length set to (cid:34) l = (cid:34) r / [ ] , solving for density in L sense. Considering this preconditioning and since the lastpanel in each side of the corner is of length smaller than (cid:34) r /
10, we set the density on those panels tozero (effectively deleting the last two panels). The
GMRES converges after 33 iterations; we use
KIFMM (with accuracy set to (cid:34) r /
10) for fast evaluation. k e k • = e (a) Laplace ( M = ) k e k • = e (b) Helmholtz ( M = ) k e k • = e − − − − − − − − (c) Yukawa ( M = ) k e k • = e (d) Stokes velocity ( M = ) k e k • = e − − − − − − − − (e) Elastostatic ( M = ) k e k • = e (f) Smooth stokes velocity ( M = ) k e k • = e − − − − − − − − (g) Smooth stokes pressure ( M = ) Fig. 4.2: T HE log OF POINTWISE ERROR . The interior Dirichlet boundary value problem is solved with known solution generatedby source points distributed over an exterior circle as shown in the lower figure in Table 4.1, apart from in (f) and (g) wherewe use the cubic flow with velocity u = [ y , x ] and pressure p = x y. Error is evaluated on the same fine grid used forvisualization ( × ). We use q = node Gauss–Legendre panels and set (cid:34) a = − in the adaptive panel quadratureset-up. M denotes the number of boundary panels. The expansion centers c are shown by black dots close to the boundary. In this paper we introduced a new quadrature scheme for the high-order accurate evaluation of layerpotentials associated with general elliptic
PDE on the domain boundary and close to it. The scheme —which builds local solution approximations using a refined evaluation and the solution of small linearsystems — relies solely on the evaluation of the underlying kernel, so is essentially
PDE -independent. Itis highly flexible, being agnostic as to the boundary condition type, the layer representation, and cru-cially, the dimension of the problem. We have analyzed the eror behavior of the scheme for Laplace and biquitous evaluation of layer potentials using
QBKIX k e k L = e − − − − − − − − (a) The log of pointwise error (b) The solution
Fig. 4.3:
QBKIX
IN A DOMAIN WITH
CORNERS . Helmholtz cases. It also fits naturally in the framework of existing fast kernel-independent algorithmsfor potential evaluation such as the
KIFMM , as it uses similar local approximations.We have tested its accuracy for three scalar- and two vector-valued Dirichlet boundary-valueproblems that are common in engineering problems. We have not attempted to optimize performance,and leave that for future work.There are several obvious extensions that have motivated this initial study that we plan to pursue:(1) Generalization to . High-order singular quadratures for surfaces are complicated, application de-pendent, and scarce. Since it requires only pointwise kernel evaluations, QBKIX is by design veryeasy to implement in using proxy and check surfaces, and would handle a wide class of PDE s. Theconstants will be larger, but the linear systems (anticipated to be of size around 10 ) would still bevery practical.(2) Generalization to other boundary conditions. QBX , and thus also
QBKIX , can apply without modifica-tion, for instance, the normal derivative of the double-layer operator, which is hypersingular.(3) Integration with
KIFMM . In this work, we only used kernel-independent
FMM for fast evaluation ofpotential on the check points. However, we expect performance gains by reusing the local expansionof
KIFMM as a
QBKIX expansion.(4) Local
QBKIX . The construction of local schemes which automatically handle general domains withthin features (i.e., with geodesically distant parts of the boundary in close proximity in space) with-out excessive refinement needed for the panel size to be on the order of feature size, is importantfor making the method practical. [ ] proposed the local version of QBX , in which only the contri-bution of the nearby panels to a target is evaluated using expansions, while contributions of moredistant panels is evaluated using standard quadrature. Implementing this idea is nontrivial however,as the end-points of the group of neighboring panels produce new singularities that can affect theconvergence rate.(5) Generalization of analysis to all kernels. As Remark 3.3 discusses, this is a nontrivial missing piecein the theoretical foundations.
Acknowledgements
We extend our thanks to Manas Rachh, Andreas Klöckner, Michael O’Neil, and Leslie Greengard for stim-ulating conversations about various aspects of this work. A.R. and D.Z. acknowledge the support of the US National ScienceFoundation (NSF) through grant DMS-1320621; A.B. acknowledges the support of the NSF through grant DMS-1216656.0 Abtin Rahimian et al.
Appendix A List of kernels
Here we list the kernels for the single- and double-layer potentials for the
PDE s considered text. In eachcase x and y are points in (cid:82) and r : = x − y . The single-layer kernel is the fundamental solution. Indouble-layer kernels, n is the unit vector denoting the dipole direction, which in the context of boundaryintegral formulation is the outward pointing normal to the surface. • Laplace: ∆ u =
0, (A.1) S ( x , y ) = − π log | r | , (A.2) D ( x , y ) = π r · n | r | , (A.3)lim y → x D ( x , y ) = − κ π , x , y ∈ Γ , (where κ is the signed curvature). (A.4) • Yukawa: ∆ u − λ u =
0, (A.5) S ( x , y ) = π K ( λ | r | ) , (A.6) D ( x , y ) = λ π r · n | r | K ( λ | r | ) , (A.7)where K , K are modified Bessel functions of the second kind of order zero and one, respectively. • Helmholtz: ∆ u + ω u =
0, (A.8) S ( x , y ) = i H ( ω | r | ) , (A.9) D ( x , y ) = i ω r · n | r | H ( ω | r | ) , (A.10)where H , H are respectively modified Hankel functions of the first kind of order zero and one. • Stokes: − ∆ u + ∇ p = ∇· u =
0, (A.11) S ( x , y ) = π (cid:18) − log | r | + r ⊗ r | r | (cid:19) , (A.12) D ( x , y ) = r · n π r ⊗ r | r | , (A.13)lim y → x D ( x , y ) = − κ π t ⊗ t , (A.14) P ( x , y ) = − π | r | (cid:18) − r ⊗ r | r | (cid:19) n . (A.15) • Navier: Linear elasticity for isotropic material with shear modulus µ and Poisson ratio ν , µ∆ u + µ − ν ∇ ∇· u =
0, (A.16) S ( x , y ) = − − ν π ( − ν ) log | r | + π ( − ν ) r ⊗ r | r | , (A.17) D ( x , y ) = − ν π ( − ν ) (cid:18) r · n + n ⊗ r − r ⊗ n | r | + − ν r · n r ⊗ r | r | (cid:19) . (A.18) biquitous evaluation of layer potentials using QBKIX References
1. Alpert, B.K.: Hybrid Gauss-trapezoidal quadrature rules. SIAM J. Sci. Comput. , 1551–1584 (1999)2. Atkinson, K.: The numerical solution of integral equations of the second kind. Cambridge University Press (1997)3. Barnett, A.H.: Evaluation of layer potentials close to the boundary for Laplace and Helmholtz problems on analytic planardomains. SIAM J. Sci. Comput. (2), A427–A451 (2014)4. Barnett, A.H., Betcke, T.: Stability and convergence of the Method of Fundamental Solutions for Helmholtz problems onanalytic domains. J. Comput. Phys. (14), 7003–7026 (2008)5. Barnett, A.H., Wu, B., Veerapaneni, S.: Spectrally-accurate quadratures for evaluation of layer potentials close to the boundaryfor the 2D Stokes and Laplace equations. SIAM J. Sci. Comput. (2014)6. Beale, J., Lai, M.C.: A method for computing nearly singular integrals. SIAM J. Numer. Anal. , 1902–1925 (2001)7. Beale, J.T., Ying, W., Wilson, J.R.: A Simple Method for Computing Singular or Nearly Singular Integrals on Closed Surfaces.Commun. Comput. Phys. pp. 1–21 (2015)8. Bogomolny, A.: Fundamental solutions method for elliptic boundary value problems. SIAM J. Numer. Anal. (4), 644–669(1985)9. Bremer, J.: On the nyström discretization of integral equations on planar curves with corners. Applied and ComputationalHarmonic Analysis (1), 45–64 (2012)10. Bremer, J., Gimbutas, Z.: A Nyström method for weakly singular integral operators on surfaces. J. Comput. Phys. ,4885–4903 (2012)11. Bremer, J., Rokhlin, V.: Efficient discretization of Laplace boundary integral equations on polygonal domains. J. Comput. Phys. , 2507–2525 (2010)12. Bremer, J., Rokhlin, V., Sammis, I.: Universal quadratures for boundary integral equations on two-dimensional domains withcorners. J. Comput. Phys. (22), 8259–8280 (2010)13. Bruno, O.P., Kunyansky, L.A.: A fast, high-order algorithm for the solution of surface scattering problems: basic implementa-tion, tests, and applications. J. Comput. Phys. , 80–110 (2001)14. Colton, D., Kress, R.: Inverse acoustic and electromagnetic scattering theory, Applied Mathematical Sciences , vol. 93, secondedn. Springer-Verlag, Berlin (1998)15. Corona, E., Rahimian, A., Zorin, D.: A tensor-train accelerated solver for integral equations in complex geometries (2015)16. Davis, P.J., Rabinowitz, P.: Methods of Numerical Integration. Academic Press, San Diego (1984)17. Duffy, M.G.: Quadrature over a pyramid or cube of integrands with a singularity at a vertex. SIAM journal on NumericalAnalysis (6), 1260–1262 (1982)18. Epstein, C.L., Greengard, L., Klöckner, A.: On the convergence of local expansions of layer potentials. SIAM J. Numer. Anal. , 2660–2679 (2013)19. Farina, L.: Evaluation of single layer potentials over curved surfaces. SIAM Journal on Scientific Computing (1), 81–91(2001)20. Ganesh, M., Graham, I.: A high-order algorithm for obstacle scattering in three dimensions. Journal of Computational Physics (1), 211–242 (2004)21. Graglia, R.D., Lombardi, G.: Machine precision evaluation of singular and nearly singular potential integrals by use of gaussquadrature formulas for rational functions. Antennas and Propagation, IEEE Transactions on (4), 981–998 (2008)22. Graham, I., Sloan, I.: Fully discrete spectral boundary integral methods for Helmholtz problems on smooth closed surfaces in (cid:82) . Numerische Mathematik (2), 289–323 (2002)23. Hackbusch, W., Sauter, S.A.: On numerical cubatures of nearly singular surface integrals arising in bem collocation. Computing (2), 139–159 (1994)24. Hao, S., Barnett, A.H., Martinsson, P.G., Young, P.: High-order accurate Nyström discretization of integral equations withweakly singular kernels on smooth curves in the plane. Adv. Comput. Math. (1), 245–272 (2014)25. Helsing, J.: Integral equation methods for elliptic problems with boundary conditions of mixed type. J. Comput. Phys. ,8892–8907 (2009)26. Helsing, J.: Solving integral equations on piecewise smooth boundaries using the RCIP method: a tutorial (2012). Preprint,34 pages, arXiv:1207.6737v3
27. Helsing, J., Ojala, R.: On the evaluation of layer potentials close to their sources. J. Comput. Phys. , 2899–2921 (2008)28. Hsiao, G., Wendland, W.L.: Boundary Integral Equations. Applied Mathematical Sciences, Vol. 164. Springer (2008)29. Järvenpää, S., Taskinen, M., Ylä-Oijala, P.: Singularity extraction technique for integral equation methods with higher orderbasis functions on plane triangles and tetrahedra. International journal for numerical methods in engineering (8), 1149–1165 (2003)30. Johnson, C.G., Scott, L.R.: An analysis of quadrature errors in second-kind boundary integral methods. SIAM Journal onNumerical Analysis (6), 1356–1382 (1989)31. Kapur, S., Rokhlin, V.: High-order corrected trapezoidal quadrature rules for singular functions. SIAM J. Numer. Anal. ,1331–1356 (1997)32. Katsurada, M.: A mathematical study of the charge simulation method. II. J. Fac. Sci. Univ. Tokyo Sect. IA Math. (1),135–162 (1989)33. Khayat, M.A., Wilton, D.R.: Numerical evaluation of singular and near-singular potential integrals. Antennas and Propagation,IEEE Transactions on (10), 3180–3190 (2005)34. Klöckner, A., Barnett, A.H., Greengard, L., O’Neil, M.: Quadrature by expansion: a new method for the evaluation of layerpotentials. J. Comput. Phys. (1), 332–349 (2013)2 Abtin Rahimian et al.35. Kolm, P., Rokhlin, V.: Numerical quadratures for singular and hypersingular integrals. Computers & Mathematics with Appli-cations (3), 327–352 (2001)36. Kress, R.: Boundary integral equations in time-harmonic acoustic scattering. Mathl. Comput. Modelling , 229–243 (1991)37. Kress, R.: On the numerical solution of a hypersingular integral equation in scattering theory. J. Comput. Appl. Math. ,345–360 (1995)38. Kress, R.: Linear Integral Equations, Appl. Math. Sci. , vol. 82, second edn. Springer (1999)39. Lowengrub, J., Shelley, M., Merriman, B.: High-order and efficient methods for the vorticity formulation of the euler equations.SIAM Journal on Scientific Computing (5), 1107–1142 (1993)40. Nachtigal, N.M., Reddy, S.C., Trefethen, L.N.: How Fast are Nonsymmetric Matrix Iterations? SIAM Journal on Matrix Analysisand Applications (3), 778–795 (1992)41. Ojala, R., Tornberg, A.K.: An accurate integral equation method for simulating multi-phase Stokes flow. J. Comput. Phys. ,145–160 (2015)42. Pozrikidis, C.: Boundary Integral and Singularity Methods for Linearized Viscous Flow. Cambridge Tests in Applied Mathe-matics. Cambridge University Press (1992)43. Quaife, B., Biros, G.: High-volume fraction simulations of two-dimensional vesicle suspensions. J. Comput. Phys. , 245–267 (2014)44. Rachh, M., Klöckner, A., O’Neil, M.: Fast algorithms for quadrature by expansion I: Globally valid expansions. arXiv preprintarXiv:1602.05301 (2016)45. Rahimian, A., Lashuk, I., Veerapaneni, S.K., Chandramowlishwaran, A., Malhotra, D., Moon, L., Sampath, R., Shringarpure,A., Vetter, J., Vuduc, R., Zorin, D., Biros, G.: Petascale Direct Numerical Simulation of Blood Flow on 200K Cores and Hetero-geneous Architectures. In: 2010 ACM / IEEE International Conference for High Performance Computing, Networking, Storageand Analysis, November, pp. 1–11. IEEE (2010)46. Schwab, C., Wendland, W.L.: On numerical cubatures of singular surface integrals in boundary element methods. NumerischeMathematik (1), 343–369 (1992)47. Sidi, A., Israeli, M.: Quadrature methods for periodic singular and weakly singular fredholm integral equations. Journal ofScientific Computing (2), 201–231 (1988)48. Strain, J.: Locally corrected multidimensional quadrature rules for singular functions. SIAM Journal on Scientific Computing (4), 992–1017 (1995)49. Tlupova, S., Beale, J.T.: Nearly singular integrals in 3d Stokes flow. Commun. Comput. Phys. (5), 1207–1227 (2013)50. Tornberg, A.K., Shelley, M.J.: Simulating the dynamics and interactions of flexible fibers in Stokes flows. Journal of Compu-tational Physics (1), 8–40 (2004)51. Trefethen, L.N., Bau III, D.: Numerical Linear Algebra. SIAM (1997)52. Veerapaneni, S.K., Rahimian, A., Biros, G., Zorin, D.: A fast algorithm for simulating vesicle flows in three dimensions. Journalof Computational Physics (14), 5610–5634 (2011)53. Yarvin, N., Rokhlin, V.: Generalized gaussian quadratures and singular value decompositions of integral operators. SIAMJournal on Scientific Computing (2), 699–718 (1998)54. Ying, L., Biros, G., Zorin, D.: A high-order 3D boundary integral equation solver for elliptic PDEs in smooth domains. J.Comput. Phys. , 247–275 (2006)55. Ying, W., Beale, J.T.: A fast accurate boundary integral method for potentials on closely packed cells. Commun. Comput. Phys.14