1 Introduction
Laplace’s equation arises in a vast array of contexts (electrostatics, harmonic functions, lowfrequency acoustics, percolation theory, homogenization theory, and the study field enhancements in vacuum insulators for example) and serves as a useful model problem for the study of general elliptic partial differential equations (PDEs). As such, effective numerical methods for quickly and robustly solving Laplace’s equation with high accuracy are desirable. Approaches based on potential theory proceed by reducing PDEs to secondkind boundary integral equations (BIEs), where the solution to the boundary value problem is represented by layer potentials on the boundary of the domain. Once these boundary integral equations are discretized the resulting linear systems are betterconditioned than those obtained by directly discretizing the PDE. When the boundary of the domain is smooth there are numerous methods for solving BIEs quickly and accurately (see
[7], for example).Near corners, however, the solutions to both the partial differential equations and corresponding boundary integral equations may have singularities, preventing the application of many traditional methods. Fortunately, a number of approaches have been developed to obviate this difficulty. One class of methods proceeds by introducing many additional degrees of freedom in the vicinity of the corners. In order to prevent the resulting linear systems from becoming intractably large one can use a variety of methods for
compressing the linear system, effectively eliminating the extra degrees of freedom added in the vicinity of the corners. Moreover, the corner refinement and compression can be done in tandem resulting in fast and accurate solvers for elliptic PDEs (see [9], [11], [16], [8] and [10] for one approach called recursive compressed preconditioning, and [5], [4],[1], and [2] for other compressionbased methods for solving Laplace’s equation). Unfortunately, this approach becomes considerably more expensive in three dimensions limiting its application in that context.Another class of methods is based on approximating the solution to the twodimensional problem by rational functions [6] with poles exponentially clustered near the corners. While this approach allows for fast evaluation of the solution near the boundary of the domain, current implementations are specialized to twodimensions, and do not scale well for large problems.
Finally, a recent approach is based on leveraging explicit representations of the solutions to the BIEs in the vicinity of the corner as sums of fractional powers depending on the angle [18, 19]. Using these representations one can construct highorder discretizations which introduce relatively few extra degrees of freedom near the corners (i.e. an amount which is comparable to the number required for smooth portions of the boundary). This approach has been used to generate efficient discretizations for Dirichlet problems for Laplace’s equation on polygonal domains[14].
In this paper we describe a method for solving Laplace’s equation on polygonal domains with Neumann boundary conditions given only a discretization of a corresponding Dirichlet problem. Our approach is based on using the discretization of a suitable adjoint problem. In particular, we show that if the transpose of the discretization of a suitable Dirichlet BIE is used, then the resulting solution will be accurate in a “weak sense”; namely, it can be used to compute inner products with smooth functions accurately, though it cannot be interpolated. We then show how this solution can be used to obtain accurate solutions to the Neumann problem arbitrarily close to a corner by solving a set of local subproblems in the vicinity of that corner.
The paper is organized as follows. In section 2 we review relevant mathematical results associated with Laplace’s equation. Section 3 describes the reduction of boundary value problems to boundary integral equations via potential theory, and reviews the analytic behavior of solutions near a corner. In sections 5 and 4 we present our numerical algorithm and the associated analysis. Finally, in section 6 we illustrate its application with several numerical experiments.
2 Mathematical preliminaries
2.1 Boundary value problems
Given a polygonal domain with boundary and outwardpointing unit normal , as well as a function we consider the following four boundary value problems.

The interior Dirichlet problem for Laplace’s equation:
(1) (2) 
The exterior Dirichlet problem for Laplace’s equation:
(3) (4) (5) 
The interior Neumann problem for Laplace’s equation:
(6) (7) (8) 
The exterior Neumann problem for Laplace’s equation:
(9) (10) (11)
Remark 2.1.
The existence and uniqueness of the solutions to the above equations is a wellknown result (see [15] for example).
3 Boundary integral equations
A classical technique for solving the four Laplace boundary value problems given above is to reduce them to boundary integral equations. Before describing this procedure we first define the single and double layer potential operators and summarize their relevant properties.
3.1 Layer potentials
Definition 3.1.
Given a function , the singlelayer potential is defined by
(12) 
where
(13) 
Similarly, the doublelayer potential is defined via the formula
(14) 
In the following we will often refer to the function as the density which generates the corresponding potential.
Definition 3.2.
For we define the kernel by
(15) 
where is the inwardpointing normal to at It will often be convenient to work instead with a parametrization of In particular, if is a counterclockwise arclength parametrization of we denote by the function defined by
(16) 
The following theorems describe the behavior of the single and double layer potentials in the vicinity of the boundary curve
Theorem 1.
Suppose the point approaches a point (where is not a corner vertex) from the inside along a path such that
(17) 
for some Then for any continuous function
(18)  
(19)  
(20) 
Similarly, if approaches a point from the outside then for any continuous function
(21)  
(22)  
(23) 
Next we define the following operator which arises in the study of Neumann boundary value problems.
Definition 3.3.
Let be the singlelayer potential operator and denote its normal derivative restricted to In particular, for
(24) 
where
The following proposition relates the normal derivative of the singlelayer operator to the doublelayer operator. Its proof follows directly from Definitions 3.1 and 3.3.
Proposition 3.1.
Let be defined as above. Let denote the normal derivative of in the sense of the previous definition. Then where denotes the adjoint operator with respect to the inner product
(25) 
where is a counterclockwise arclength parametrization of In particular, for all
(26) 
and
(27) 
3.2 Reduction of boundary value problems
In this section we describe the conversion of the Laplace boundary value problems (interior Dirichlet, exterior Dirichlet, interior Neumann, and exterior Neumann) to secondkind integral equations.
Theorem 2 (Interior Dirichlet problem for Laplace’s equation).
For every in there exists a unique which satisfies
(28) 
Moreover, the solution to the interior Dirichlet problem for Laplace’s equation with boundary data is given by for all
Theorem 3 (Exterior Dirichlet problem for Laplace’s equation).
For every in there exists a unique which satisfies
(29) 
for all Moreover, the solution to the exterior Dirichlet problem for Laplace’s equation with boundary data is given by for all
Theorem 4 (Interior Neumann problem for Laplace’s equation).
For every in such that there exists a unique which satisfies
(30) 
Moreover, the solution to the interior Neumann problem for Laplace’s equation with boundary data is given by for all
Theorem 5 (Exterior Neumann problem for Laplace’s equation).
For every in there exists a unique which satisfies
(31) 
Moreover, the solution to the interior Neumann problem for Laplace’s equation with boundary data is given by for all
3.3 Corner expansions
In the remainder of this section, we assume is an open wedge with sides of length one and interior angle with Let be an arc length parametrization of and be the inwardpointing normal to The following theorem gives an explicit representation of the solutions of the boundary integral equation (28) in this geometry.
Theorem 6 ([18]).
Suppose that and that is a positive integer. Let and denote the ceiling and floor functions, respectively, and define and by the following formulas
(32)  
(33)  
(34)  
(35) 
Suppose further that is defined via the formula
(36) 
where and are arbitrary real numbers and the functions and are defined as follows
(37) 
(38) 
If is defined by
(39) 
and is defined by (36) then there exist two sequences of real numbers , and such that
(40) 
for all Conversely, suppose that has the form (40). Suppose further that is an arbitrary positive integer. Then, for all angles there exist unique real numbers and such that defined by (36), solves equation (39) to within an error
Remark 3.1.
A similar result holds for the case where the identity term in (39) is replaced by its negative; the change in sign corresponds to replacing the boundary integral equation for the interior Dirichlet problem with the boundary integral equation corresponding to exterior Dirichlet problem. Similar expansions also hold for both the exterior and interior Neumann problems, in which case the singular powers are obtained by subtracting one from the singular powers arising in the Dirichlet problem.
The following corollary, proved in [18] gives a characterization of the solutions to the Dirichlet and Neumann boundary integral equations in the vicinity of a corner.
Corollary 3.1.
Let be the boundary of a polygonal region and suppose one of its corners has interior angle where Let be an arclength parametrization of in the vicinity of the corner, with coinciding with the corner. If the boundary data, is analytic on either side of the corner then there exist unique real numbers and such that the density, defined by (36) satisfies the interior Dirichlet boundary integral equation to within an error for within of the corner. For the Neumann problems the representation is the same with the powers in the expansion reduced by one.
4 Numerical preliminaries
In this section we summarize the numerical tools which are necessary for the main result. In particular we summarize the method for discretizing the boundary integral equation for the Dirichlet problem described in [14], which uses the expansion in Theorem 6.
4.1 Discretization of the Dirichlet problem
In this section we sketch an algorithm for solving the interior Dirichlet boundary integral equation using a Nyström method; the exterior Dirichlet boundary integral equation can be discretized in a similar way. See [14] for a thorough description of the method.
The Nyström method proceeds as follows. We begin by constructing a discretization of the boundary with nodes and weights which enable interpolation of the left and righthand sides of the boundary integral equation
(41) 
with precision In other words, given for the values and can be obtained for all to within
Once these nodes and weights have been generated we proceed by enforcing equality of (41) at the discretization nodes, which yields the system of equations
(42) 
We note that scaling by the square root of the weights in the above equation is equivalent to solving the problem in the sense, and results in discretized operators with condition numbers which are close to those of the original physical systems [2]. The new unknowns are Next, for each interpolation node we find a collection of weights such that
(43) 
resulting in the linear system
(44) 
4.1.1 Obtaining interpolation nodes
The boundary is separated into a collection of intervals which are at least a fixed distance (measured in terms of arclength) away from a corner and the collection of intervals of length centered about each corner. The former are discretized using a standard smooth quadrature rule such as nested GaussLegendre quadrature while the latter are discretized using a custom set of interpolation nodes constructed in the following way.
First, all functions of the form are discretized using nested GaussLegendre panels in and a single GaussLegendre panel in . This creates a matrix where denotes the number of spatial discretization points and denotes the number of chosen. and are increased until it is guaranteed that using Lagrange interpolation from the nested discretization the function can be interpolated to within an error less than on the interval for any
in the specified range. A singular value decomposition is then performed on the
matrix. Let denote the number of singular values greater thanThe right singular vectors correspond to discretizations of an orthonormal set of functions
such that is in the span of to within an accuracy ofFinally, a set of interpolation points and quadrature weights are chosen for such that the matrix is wellconditioned. In practice suitable interpolation points can be obtained by using the roots of and calculating the corresponding weights by solving a linear system. The corresponding discretization nodes and weights for the cornercontaining intervals of are obtained by suitable translations and scalings of and
4.1.2 Construction of quadrature rules
Once the discretization has been constructed it is necessary to construct appropriate quadrature for the integrals appearing in equation (42). When and do not belong to the same corner panel (in particular when either is not itself contained in a corner panel) then the weights and nodes associated with the discretization can be used as the quadrature rule. When corresponds to a corner panel special care must be taken. Instead, using an algorithm for generating generalized Gaussian quadratures [3], quadrature nodes are chosen which integrate
(45) 
where is a suitably scaled and translated copy of the singular function obtained in the discretization step, and for ease of exposition we assume that the corner panel corresponds to in the parametrization with corresponding to the corner itself. Moreover, in light of symmetry between the two legs of the wedge it suffices to design quadratures assuming lies in the half of a corner panel parametrized by
Remark 4.1.
Due to scale invariance, it suffices to compute quadratures for
(46) 
where was one of the original discretization nodes generated on the interval
Remark 4.2.
By interpolating from the discretization nodes to these quadrature nodes we obtain a set of weights such that if correspond to the discretization of a corner parametrized by with corresponding to the corner then
(47) 
for all and
After all the quadratures have been constructed the result is an linear system the solution of which gives an approximation to sampled at the discretization nodes.
Definition 4.1.
Let denote the set of functions which can be interpolated from their values at the discretization nodes to any point in with a relative accuracy of That is to say that for if denotes the function obtained by interpolating using the values then
The results of this algorithm are summarized in the following theorem (see ).
Theorem 7.
Let be the matrix obtained by discretizing the interior Dirichlet problem in the preceding manner. In particular if is piecewise analytic and then
(48) 
can be interpolated to a function which is within of the true density in an sense.
4.2 Discretization of the Neumann problem
In principle a similar method could be employed to discretize the Neumann boundary integral equations. Unfortunately, the singular nature of the powers (the smallest in the expansion given in Theorem 6 lies in the range ) makes it difficult to produce universal discretizations and quadratures which work for large ranges of angles. When the above method is run on these problems, discretization nodes tend to accumulate close to the corner (within ). Apart from posing certain numerical challenges, it also makes the task of finding suitable quadrature formulae difficult. Instead, a different set of discretization nodes and a different set of quadrature nodes can be constructed for each angle, though this would significantly increase the precomputation cost.
Finally, in many applications one already has a discretization of the Dirichlet problem. For example, when considering Laplace transmission problems or triple junction problems one has to solve two decoupled boundary integral equations: one of them a Dirichlettype boundary integral equation with the diagonal term scaled and the other a Neumanntype boundary integral equation with the identity term scaled (see [12] and [13] for example). In such cases it is convenient to reuse the Dirichlet discretization for the Neumann problem.
5 Numerical apparatus
5.1 Adjoint discretization
The following lemma relates the discretization of the inverse of an operator to the adjoint of the discretization of its inverse. Its proof follows directly from the definition of the adjoint and is omitted.
Lemma 1.
Suppose is a bounded invertible operator and that is an operator such that
(49) 
for all and in some subspace Here denotes the inner product on and denotes the norm for Then, for all functions and in
(50) 
where denotes the adjoint.
The following corollary follows immediately from the previous result.
Corollary 5.1.
Let be the matrix obtained by discretizing the interior Dirichlet problem and be the collection of functions given by Definition 4.1. Then for all functions
(51) 
where are the discretizations of and scaled by the square roots of the discretization weights, and is the solution to the exterior Neumann problem with boundary data .
Hence a discretization of the Neumann problem can be obtained simply by taking the adjoint of the Dirichlet problem. The resulting density obtained is accurate in a weak sense, ie. its inner products against functions in are accurate to within an error of
We conclude this section with a few remarks.
Remark 5.1.
We observe that if the solution to the boundary value problem is being calculated at a point more than one panel length away from the boundary curve then the Neumann density obtained using the above result will give an accuracy of ie. the function Thus accurate values of the solution in the farfield can be obtained almost immediately.
Remark 5.2.
Similarly, if the point at which the solution to the Neumann boundary value problem is to be calculated lies close to a smooth panel then the density near that point can be interpolated to a finer set of quadrature points and the value of can once again be obtained to precision We note, however, that in general the density in the vicinity of a corner cannot be interpolated accurately. This follows from the fact that the interpolation scheme constructed is only guaranteed to interpolate the powers arising in the Dirichlet problem accurately near the corner. The collection of singular powers arising in Neumann problems contain negative powers which are not contained in this set and hence are not interpolated accurately.
5.2 Weak corner resolving
In this section we address the problem highlighted in the previous one; namely, the accurate evaluation of the solution to the exterior Neumann problem in the vicinity of a corner. Our approach is based on the observation that the potential generated by the density on the boundary outside of a sufficiently small neighborhood of the corner is smooth when evaluated in the vicinity of the corner. This allows us to convert the problem of evaluating the potential near the corner (given the approximation to the density obtained using the adjoint approach described in the previous section) into a purely local one. In particular, we rediscretize only a small neighborhood of the corner which in turn allows us to evaluate the potential arbitrarily close to the corner to within a small factor of machine precision.
In the following we assume that we are given a discretization of the interior Dirichlet boundary integral equation (28) with nodes and corresponding weights In particular, we assume that the discretization nodes are obtained by subdividing the boundary into panels. Those panels which contain a vertex are discretized using a custom discretization scheme (see Section 4.1) while the remaining panels are discretized using a standard smooth quadrature rule (such as GaussLegendre or Chebyshev nodes). In the following we assume that an point GaussLegendre quadrature rule is used and the corner panels are discretized using nodes (together with a collection of orthonormal functions on that interval ).
Additionally, we denote the discretization of the interior Dirichlet operator (using the custom quadratures described in Section 4.1) by Let where and is the righthand side of the exterior Neumann problem. Finally, let be the approximation to the density (scaled by the square roots of the weights) obtained by solving the linear system
(52) 
For notational convenience we let be a counterclockwise arclength parametrization of such that corresponds to a vertex and corresponds to a corner panel.
For a panel with discretization nodes corresponding to a GaussLegendre panel the density is smooth and thus it is expected to be wellrepresented in the basis of Legendre polynomials (shifted and scaled to the interval ). Hence standard interpolation techniques can be used to obtain an accurate approximation to the density on the interval . Typically we use 16th order GaussLegendre panels and choose their sizes so that their length is no more than their distance to the nearest corner. This latter choice guarantees that for any there exists an such that if the GaussLegendre panels are discretized using an point GaussLegendre rule then the density on that panel can be interpolated to relative precision in an sense. (We discuss a sketch of a proof in appendix B)
For corner panels the nodes were constructed to enable stable interpolation of densities on the interval  assuming for simplicity that the corner is at and the panel is of length As mentioned above, the density is expected to contain terms of the form for some finite collection of in the interval and hence will not in general be stably interpolable on the interval However, it is possible to use the density obtained using (52) to construct a sequence of nested problems in the neighborhood of the corner, the solutions of which enable accurate interpolation of the density arbitrarily close to the vertex. The number of these problems depends only on the distance of the closest evaluation point to the corner. In particular, if is the smallest distance of an evaluation point from the corner then only levels are required. Each problem involves the solution of a small linear system (typically less than ) and as such can be performed quickly. Furthermore, we note that the algorithm can be easily parallelized to treat multiple corners concomitantly.
We begin with the following proposition, the proof of which follows immediately from the definition of the kernel and is omitted.
Proposition 5.1.
Suppose that be a piecewiseanalytic function in and is the approximation to the Neumann density obtained using the adjoint of the discretization for the interior Dirichlet boundary integral equation. Further suppose that the discretization nodes are ordered so that correspond to the corner panel associated with the interval correspond to the GaussLegendre panel immediately to the left associated with the interval and to the GaussLegendre panel immediately to the right associated with the interval Then
(53) 
is an analytic function of for all
In light of this we consider the following integral equation
(54) 
We note that the solution to (54) is equal to the solution of the original boundary integral equation (31) restricted to the interval Taking the adjoint of (54) we obtain
(55) 
which is a Dirichlet boundary integral equation for a wedge with a piecewise analytic righthand side. In particular, we can discretize the operator using the method summarized in the previous section. Specifically, we subdivide the interval into three subintervals and On and we place standard GaussLegendre discretization nodes, while on we use the custom discretization scheme for corners, outlined in Section 4.1 (see for a detailed description of the method). On the intervals and we use the same discretization nodes and weights as in the original system for those intervals (we call these panels and respectively). Let denote the righthand side of (54) evaluated at these discretization nodes and scaled by the square roots of the corresponding weights. Let be the discretization of the interior Dirichlet problem operator (ie. the operator acting on on the lefthand side of (55)). We note that due to the scale invariance of Laplace’s equation for polygonal domains the portion of corresponding to the selfinteraction of is a submatrix of the original matrix All other blocks can be generated using the discretization nodes as quadrature nodes.
The analysis of the previous section then shows that if is the solution of the equation
(56) 
then gives a weak solution to the integral equation (54), i.e. for any function which is analytic on and the inner product can be calculated to precision using the solution Moreover, since the true density is smooth on and the GaussLegendre discretization allows accurate interpolation of the density on those regions.
Remark 5.1.
Though the above method produces a viable method for reducing the problem, as written the reduction is nonlocal — in order to compute the righthand side for the subproblem one must evaluate contributions from the rest of the domain.
The following theorem shows that the righthand side can be computed only using local data (i.e. values of the weak solution in the vicinity of the corner).
Theorem 8.
Suppose that is the discretization of the righthand side of (55) corresponding to nodes Further suppose that is the matrix with entries where are the orthonormal functions on spanning on that interval. Let be the vectorvalued function defined by
(57) 
and be the approximation to the solution in the vicinity of the corner obtained by solving the original system (52). Then for all In particular, if is the vector of length with entries defined by
(58) 
then
Proof.
We begin by observing that both and are analytic on the interval In particular, they can be accurately interpolated using on the interval Hence,
(59) 
A similar argument shows that is interpolable on and On the other hand, by construction,
(60) 
Let be the submatrix of corresponding to the first rows and columns of and be the submatrix of corresponding to selecting the first rows of and all but the first columns of Using this notation, the first rows of (60) can be rewritten as
(61) 
The first term on the righthand side is Substituting this into the previous equation, we obtain
(62) 
The result follows by substituting the above equality into (59). ∎
This can be iterated to obtain an interpolable approximation to the density on In particular, we consider the restriction of the exterior Neumann integral equation, as well as its, adjoint to the interval and For the righthand side we use the original righthand side minus the contribution from the remainder of the domain. In particular, if we define
(63) 
then restricted to the interval satisfies
(64) 
The corresponding adjoint equation is given by
(65) 
Once again, we divide into three intervals and and discretize each interval as before. After solving the corresponding discretization of (64) using the adjoint of the discretization of the integral operator appearing in (65) we obtain a weak solution of on the interval which can be interpolated on and to within precision
This process can be repeated an arbitrary number of times to yield a sequence of solutions together with corresponding intervals and on which it can be interpolated.
Note that if is a point a distance away from the corner then after such subdivisions will be at least twice the corner panel length away from the corner. Thus will be smooth when restricted to the corner panel and hence will be integrated accurately using the corner panel discretization nodes and weights.
6 Numerical results
6.1 Accuracy
In this section, we demonstrate the accuracy of the proposed numerical method (both in the weak sense described above, as well as in the classical sense after sufficiently many resolves) on the triangular domain shown below. The reference solution for each of the examples is computed using a discretization with a graded mesh in the vicinity of the corners, where the smallest panel at the corner is times the length of the first macroscopic panel away from the corner (see fig. 1). In these examples, the solutions are computed via dense linear solves.
Remark 1.
Though is significantly smaller than machine precision, the matrix entries corresponding to the corner interactions can be computed accurately by translating the corners to the origin when computing interactions of nearby points.
Remark 2.
Simple arguments from complex analysis show that when using graded meshes, in order to obtain full machine precision () for solutions of the Neumann problem at any point in the interior at least away from a corner, it suffices to choose the smallest panel (i.e. the size of the panel closest to the corner) to be of size However, resulting values of the density will not be accurate to machine precision at all nodes. In fact the quality of the density deteriorates as one approaches the corner. Thus, in order to obtain accurate point values of the density to machine precision at all points which are at least away from the corner, we use a smallest panel size of .
The potential at target locations which are sufficiently far from the boundary (i.e. at least one panel length away from every panel) is the inner product of the density with a smooth function and hence can be computed accurately without resolving (see creftype 5.1). For a target location , we compute the potential via the formula,
(66) 
In fig. 2, we compute the error in the solution at target locations for a scattering problem whose right hand side is given by a collection of three interior charges
(67) 
where the locations are denoted by square dots in fig. 2. Note that the density plotted as a function of arclength goes to infinity at the corner vertices, indicating that the native Dirichlet discretization presented in section 4.1 wouldn’t have sufficed. However, the potential in the volume is accurate to 14 digits at target locations away from the boundary.
Another example of a “weak quantity” is the polarization tensor associated with a domain. This requires the solution of the exterior problems with boundary data
or . Let and denote the corresponding solutions. The polarization tensor can be expressed in terms of the solutions and as(68) 
The polarization tensor as computed by the reference solution, and the error in computation using the adjoint discretization are given by
(69) 
In order to demonstrate the accuracy of the corner resolving approach in obtaining the true density at the corner panels, we apply the procedure discussed in section 5.2 iteratively, and compare the obtained density with the reference density after 20,40,60, and 80 iterations of resolves in the vicinity of one of the corners. The reference density and the errors are shown in fig. 3. Furthermore, to highlight the need for special purpose discretizations in the vicinity of corners in the adjoint discretization, we also compare the solution computed using a graded mesh in the vicinity of corners, where the size of the smallest panels for both discretizations are equal.
After resolving the density, the solution is evaluated on a tensor product polar grid, where the grid is exponentially spaced in the radial direction and equispaced in the angular direction. For evaluation points (targets) close to panels which are not at the corner, we use adaptive integration in order to resolve the nearsingular behavior of the kernel for accurate computation of the integrals. For target locations close to the corner panel, since we do not have the capability to interpolate the density, we use the underlying smooth quadrature rules for computing their contribution. The reference solution and the errors are demonstrated in fig. 4.
6.2 Performance
In this section, we demonstrate the performance of the solver by solving a scattering problem in the exterior of a “broken wheel” region. The boundary data is given by
Comments
There are no comments yet.