On the discretization of Laplace's equation with Neumann boundary conditions on polygonal domains

In the present paper we describe a class of algorithms for the solution of Laplace's equation on polygonal domains with Neumann boundary conditions. It is well known that in such cases the solutions have singularities near the corners which poses a challenge for many existing methods. If the boundary data is smooth on each edge of the polygon, then in the vicinity of each corner the solution to the corresponding boundary integral equation has an expansion in terms of certain (analytically available) singular powers. Using the known behavior of the solution, universal discretizations have been constructed for the solution of the Dirichlet problem. However, the leading order behavior of solutions to the Neumann problem is $O(t^{\mu})$ for $\mu \in (-1/2,0)$ depending on the angle at the corner (compared to $O(C+t^{\mu})$ with $\mu>1/2$ for the Dirichlet problem); this presents a significant challenge in the design of universal discretizations. Our approach is based on using the discretization for the Dirichlet problem in order to compute a solution in the"weak sense"by solving an adjoint linear system; namely, it can be used to compute inner products with smooth functions accurately, but it cannot be interpolated. Furthermore we present a procedure to obtain accurate solutions arbitrarily close to the corner, by solving a sequence of small local subproblems in the vicinity of that corner. The results are illustrated with several numerical examples.


Introduction
Laplace's equation arises in a vast array of contexts (electrostatics, harmonic functions, low-frequency 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 second-kind 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 better-conditioned 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 compression-based 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 two-dimensions, 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 high-order 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 4 and 5 we present our numerical algorithm and the associated analysis. Finally, in section 6 we illustrate its application with several numerical experiments.

Mathematical preliminaries 2.1 Boundary value problems
Given a polygonal domain Ω ⊂ R 2 with boundary Γ and outward-pointing unit normal ν, as well as a function f : Γ → R, we consider the following four boundary value problems.
1. The interior Dirichlet problem for Laplace's equation: 2. The exterior Dirichlet problem for Laplace's equation: 3. The interior Neumann problem for Laplace's equation: 4. The exterior Neumann problem for Laplace's equation: Remark 2.1. The existence and uniqueness of the solutions to the above equations is a well-known result (see [15] for example).

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.

Layer potentials
Definition 3.1. Given a function σ : Γ → R, the single-layer potential is defined by where Similarly, the double-layer potential is defined via the formula In the following we will often refer to the function σ as the density which generates the corresponding potential.
where ν(x) is the inward-pointing normal to Γ at x. It will often be convenient to work instead with a parametrization of K. In particular, if γ : [0, L] → Γ is a counterclockwise arclength parametrization of Γ, we denote by k : 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 x approaches a point x 0 = γ(t 0 ) (where x 0 is not a corner vertex) from the inside along a path such that for some α > 0. Then for any continuous function σ : Γ → R, Similarly, if x approaches a point x 0 = γ(t 0 ) from the outside then for any continuous function σ : Γ → R, Next we define the following operator which arises in the study of Neumann boundary value problems. Definition 3.3. Let S be the single-layer potential operator and ν · ∇S denote its normal derivative restricted to Γ. In particular, for x 0 ∈ Γ, where γ(t 0 ) = x 0 .
The following proposition relates the normal derivative of the single-layer operator to the double-layer operator. Its proof follows directly from Definitions 3.1 and 3.3.

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 second-kind integral equations.
Moreover, the solution to the interior Dirichlet problem for Laplace's equation with boundary data f is given by u(y) = D[σ](y) for all y ∈ Ω.
Moreover, the solution to the interior Neumann problem for Laplace's equation with boundary data f is given by u(y) = S[σ](y) for all y ∈ Ω.
Moreover, the solution to the interior Neumann problem for Laplace's equation with boundary data f is given by u(y) = S[σ](y) for all y ∈ Ω.

Corner expansions
In the remainder of this section, we assume Γ is an open wedge with sides of length one and interior angle πα with 0 < α < 2. Let γ : [−1, 1] → Γ be an arc length parametrization of Γ and ν : [−1, 1] → R 2 be the inward-pointing 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 0 < α < 2 and that N is a positive integer. Let · and · denote the ceiling and floor functions, respectively, and define L, L, M , and M by the following formulas Suppose further that σ is defined via the formula where b 0 , b 1 , . . . , b N and c 1 , c 2 , . . . , c N are arbitrary real numbers and the functions σ α,N (i) and ν α,N (i) are defined as follows otherwise.

(38)
If f is defined by and σ is defined by (36) then there exist two sequences of real numbers β 0 , β 1 , . . . and γ 0 , γ 1 , . . . such that  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 α ∈ (0, 2). Let γ : (−δ, δ) → R 2 be an arclength parametrization of Γ in the vicinity of the corner, with γ(0) coinciding with the corner. If the boundary data, f, is analytic on either side of the corner then there exist unique real numbers b 0 , b 1 , . . . , b N and c 0 , c 1 , . . . , c N such that the density, ρ, defined by (36) satisfies the interior Dirichlet boundary integral equation to within an error O(t N +1 ) for t within δ of the corner. For the Neumann problems the representation is the same with the powers in the expansion reduced by one.

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.

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 s 1 , . . . , s N , and weights w 1 , . . . , w N , which enable interpolation of the left-and right-hand sides of the boundary integral equation with precision . In other words, given f ( 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 We note that scaling by the square root of the weights in the above equation is equivalent to solving the problem in the L 2 sense, and results in discretized operators with condition numbers which are close to those of the original physical systems [2]. The new unknowns are resulting in the linear system

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 2δ centered about each corner. The former are discretized using a standard smooth quadrature rule such as nested Gauss-Legendre quadrature while the latter are discretized using a custom set of interpolation nodes constructed in the following way. First, all functions of the form x µ , µ ∈ {0} ∪ [1/2, 50], x ∈ [0, 1] are discretized using nested Gauss-Legendre panels in x and a single Gauss-Legendre panel in µ. This creates a N × M matrix where N denotes the number of spatial discretization points r i and M denotes the number of µ j chosen. M and N are increased until it is guaranteed that using Lagrange interpolation from the nested discretization the function x µ can be interpolated to within an L 2 error less than on the interval [0, 1] for any µ in the specified range. A singular value decomposition is then performed on the N × M matrix. Let K denote the number of singular values greater than . The right singular vectors correspond to discretizations of an orthonormal set of functions φ 1 , . . . , φ K such that x µ is in the span of φ 1 , . . . , φ K to within an accuracy of .
Finally, a set of interpolation points x j , j = 1, . . . , K and quadrature weights w j , j = 1, . . . , K are chosen for φ 1 , . . . , φ K such that the matrix U ij = φ i (x j ) √ w j is wellconditioned. In practice suitable interpolation points can be obtained by using the roots of φ K+1 and calculating the corresponding weights by solving a linear system. The corresponding discretization nodes and weights for the corner-containing intervals of Γ are obtained by suitable translations and scalings of {x j } and {w j }.

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 s i and t 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 s i 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 whereφ j 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 t = 0 corresponding to the corner itself. Moreover, in light of symmetry between the two legs of the wedge it suffices to design quadratures assuming s j lies in the half of a corner panel parametrized by (−δ, 0).
where x j was one of the original discretization nodes generated on the interval [0, 1].
Remark 4.2. By interpolating from the discretization nodes to these quadrature nodes we obtain a set of weightsW i,j such that if s 1 , . . . , s 2K correspond to the discretization of a corner parametrized by (−δ, δ) with 0 corresponding to the corner then for all i = 1, . . . , 2K and m = 1, . . . , K.
After all the quadratures have been constructed the result is an N ×N linear system the solution of which gives an approximation to σ sampled at the discretization nodes. Theorem 7. Let A be the N × N matrix obtained by discretizing the interior Dirichlet problem in the preceding manner. In particular if f ∈ S is piecewise analytic and f = ( can be interpolated to a functionσ which is within of the true density σ in an L 2 -sense.

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 (−1/2, 0)) 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 10 −14 ). 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 Dirichlet-type boundary integral equation with the diagonal term scaled and the other a Neumann-type 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.

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.
is a bounded invertible operator and that A is an operator such that for all f and g in some subspace S ⊂ L 2 ([0, L]). Here ·, · denotes the inner product on L 2 ([0, L]) and · denotes the norm for L 2 ([0, L]). Then, for all functions f and g in S where * denotes the adjoint.
The following corollary follows immediately from the previous result.
where f , g are the discretizations of f and g scaled by the square roots of the discretization weights, and σ is the solution to the exterior Neumann problem with boundary data f .
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 S 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 y ∈ R 2 \ Ω 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 K(y, γ(t)) ∈ S . Thus accurate values of the solution in the far-field can be obtained almost immediately.
Remark 5.2. Similarly, if the point y ∈ R 2 \ Ω 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 u(y) 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.

Weak corner re-solving
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 re-discretize 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 x 1 , . . . , x N and corresponding weights w 1 , . . . , w N . 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 Gauss-Legendre or Chebyshev nodes). In the following we assume that an M -point Gauss-Legendre quadrature rule is used and the corner panels are discretized using P nodes (together with a collection of orthonormal functions on that interval φ 1 , . . . , φ P ).
Additionally, we denote the discretization of the interior Dirichlet operator (using the custom quadratures described in Section 4.1) by A.
√ w i and f : ∂Ω → R is the right-hand 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 For notational convenience we let γ : [−δ, L − δ] → ∂Ω be a counterclockwise arclength parametrization of ∂Ω such that γ(0) corresponds to a vertex and γ[−δ, δ] corresponds to a corner panel.
For a panel γ([s 1 , s 2 ]) with discretization nodes x i , . . . , x i+M corresponding to a Gauss-Legendre panel the density is smooth and thus it is expected to be well-represented in the basis of Legendre polynomials (shifted and scaled to the interval [s 1 , s 2 ]). Hence standard interpolation techniques can be used to obtain an accurate approximation to the density σ on the interval s 1 ≤ s ≤ s 2 . Typically we use 16th order Gauss-Legendre 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 > 0 there exists an M such that if the Gauss-Legendre panels are discretized using an M -point Gauss-Legendre rule then the density on that panel can be interpolated to relative precision in an L 2 -sense. (We discuss a sketch of a proof in appendix B) For corner panels the nodes were constructed to enable stable interpolation of densities s µ , µ ∈ 0 ∪ [1/2, 50], on the interval s ∈ (−δ, δ) -assuming for simplicity that the corner is at 0 and the panel is of length 2δ. As mentioned above, the density is expected to contain terms of the form s µ for some finite collection of µ in the interval (−1/2, 1/2), 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 r is the smallest distance of an evaluation point from the corner then only log 2 r/δ levels are required. Each problem involves the solution of a small linear system (typically less than 100 × 100) 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 k and is omitted.
Proposition 5.1. Suppose that f be a piecewise-analytic function in S and σ = (A T ) −1 f 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 s 1 , . . . , s P correspond to the corner panel associated with the interval (−δ, δ), s P +1 , . . . , s P +M correspond to the Gauss-Legendre panel immediately to the left associated with the interval (−2δ, −δ), and s P +M +1 , . . . , s P +2M to the Gauss-Legendre panel immediately to the right associated with the interval (δ, 2δ). Then is an analytic function of t for all t ∈ (−2δ, 2δ).
In light of this we consider the following integral equation We note that the solution to (54) is equal to the solution of the original boundary integral equation (31) restricted to the interval [−2δ, 2δ]. Taking the adjoint of (54) we obtain which is a Dirichlet boundary integral equation for a wedge with a piecewise analytic right-hand side. In particular, we can discretize the operator using the method summarized in the previous section. Specifically, we subdivide the interval [−δ, δ] into three subintervals I 0 = [−δ, δ/2], L 0 = [−δ/2, δ/2], and J 0 = [δ/2, δ]. On I 0 and J 0 we place standard Gauss-Legendre discretization nodes, while on L 0 we use the custom discretization scheme for corners, outlined in Section 4.1 (see [] for a detailed description of the method). On the intervals [−2δ, −δ] and [δ, 2δ] we use the same discretization nodes and weights as in the original system for those intervals (we call these panels K 0 and Q 0 respectively). Let f 0 denote the right-hand side of (54) evaluated at these discretization nodes and scaled by the square roots of the corresponding weights. Let A 0 be the discretization of the interior Dirichlet problem operator (ie. the operator acting on σ on the left-hand side of (55)). We note that due to the scale invariance of Laplace's equation for polygonal domains the portion of A 0 corresponding to the self-interaction of L 0 is a submatrix of the original matrix A. All other blocks can be generated using the discretization nodes as quadrature nodes. The analysis of the previous section then shows that if σ 0 is the solution of the equation then σ 0 gives a weak solution to the integral equation (54), i.e. for any function g which is analytic on [−2δ, 0] and [0, 2δ] the inner product g, σ can be calculated to precision using the solution σ 0 . Moreover, since the true density σ is smooth on [δ, 2δ] and [−2δ, −δ] the Gauss-Legendre 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 non-local -in order to compute the right-hand side for the sub-problem one must evaluate contributions from the rest of the domain.
The following theorem shows that the right-hand side f 0 can be computed only using local data (i.e. values of the weak solution in the vicinity of the corner).
Then δ). In particular, iff 0 is the vector of length N 0 with entries defined by Proof. We begin by observing that both f and h are analytic on the interval (−δ, δ). In particular, they can be accurately interpolated using φ 1 , . . . , φ P on the interval (−δ, δ). Hence, Let A 0 be the (P + 2M ) × (P + 2M ) submatrix of A corresponding to the first P + 2M rows and columns of A, and A red be the P + 2M × (N − P − 2M ) submatrix of A corresponding to selecting the first P + 2M rows of A and all but the first P + 2M columns of A. Using this notation, the first P + 2M rows of (60) can be re-written as The first term on the right-hand side is (h(s 1 ), . . . , h(s P +2M )) T . Substituting this into the previous equation, we obtain The result follows by substituting the above equality into (59).
This can be iterated to obtain an interpolable approximation to the density on L 0 = [−δ, δ]. In particular, we consider the restriction of the exterior Neumann integral equation, as well as its, adjoint to the interval I 0 and J 0 . For the right-hand side we use the original right-hand side f minus the contribution from the remainder of the domain. In particular, if we define The corresponding adjoint equation is given by Once again, we divide L 0 into three intervals I 1 , L 1 , and J 1 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 I 0 ∪ L 0 ∪ J 0 which can be interpolated on I 0 , and J 0 to within precision .
This process can be repeated an arbitrary number of times to yield a sequence of solutions σ j , j = 0, 1, 2, . . . together with corresponding intervals I 0 , I 1 , . . . and J 0 , J 1 , . . . on which it can be interpolated.
Note that if x is a point a distance r away from the corner then after J = 1+log 2 r/d such subdivisions x will be at least twice the corner panel length away from the corner. Thus K(x, ·) will be smooth when restricted to the corner panel [−δ/2 J , δ/2 J ] and hence will be integrated accurately using the corner panel discretization nodes and weights.
6 Numerical results

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 re-solves) 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 2 −200 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 2 −200 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. Figure 1: Problem domain and panel discretization of the boundary. The discretization on the left is based on using the Dirichlet discretization at corner panels (indicated by blue panels) discussed in section 4.1, while the discretization on the right is a sample discretization with 2 levels of refinement in the vicinity of the corner. All the panels in black are discretized using scaled Gauss-Legendre nodes. The square ticks indicate location of the charges x j for defining the boundary data for the scattering problem.

Remark 2.
Simple arguments from complex analysis show that when using graded meshes, in order to obtain full machine precision (∼ 1.11 × 10 −16 ) for solutions of the Neumann problem at any point in the interior at least 10 −16 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 2 −100 . 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 2 −100 away from the corner, we use a smallest panel size of 2 −200 .
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 re-solving (see remark 5.1). For a target location y, we compute the potential via the formula, 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 where the locations x j 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 f 1 = ν 1 or f 2 = ν 2 . Let σ 1 and σ 2 denote the corresponding solutions. The polarization tensor can be expressed in terms of the solutions σ 1 and σ 2 as 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 re-solving 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 re-solving 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 near-singular 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.

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 where there is one x j located in each of the spokes, one of the x j is in the central disc, and the remaining 50 x j are chosen randomly in the exterior of the bounding disc containing the domain. The strengths c j are chosen such that they average to 0. The domain contains 108 corners, was discretized using 22240 nodes and required 105 iterations to converge to a residue of 10 −15 . The matrix at each iteration was applied using an FMM whose tolerance was also set to 10 −15 . The solution was computed in 15 secs, and plotted at a 500 × 500 grid of targets in 6.5 secs. All of the results have been computed on a single core on a Macintosh machine with Intel core i5 2.3GHz processors. In fig. 5, we plot the scattered field, the boundary data, and the computed density.

Conclusion and future work
In this paper we described a method for obtaining solutions to Laplace's equation with Neumann boundary conditions on polygonal domains given an accurate discretization of a corresponding Dirichlet problem. The resulting solutions are accurate in a "weak sense", allowing evaluation of the solution at points which are located sufficiently far from the boundary of the domain. We then presented a method for using these "weak solutions" to obtain accurate solutions to the Neumann problem in an L ∞ -sense arbitrarily close to the corner in a computationally efficient manner. Though the present paper treats only Laplace's equation for polygonal domains, the method shown here extends much more broadly. In particular, the approach easily extends to accommodate curved boundaries. Moreover, in addition to Laplace's equation, this approach can be easily adapted to solve the Helmholtz equation and the biharmonic equation with analogous boundary conditions for which the nature of singularities of corresponding integral equations have been analyzed [17,20]. A manuscript detailing this extension is currently in preparation.

Acknowledgments
J. Hoskins was supported in part by AFOSR FA9550-16-1-0175 and by the ONR (award no. N00014-14-1-0797). The authors would like to thank Alex Barnett, Leslie Greengard, Michael O'Neil, and Vladimir Rokhlin for many useful discussions, and Jeremy Magland for providing sector plotting tools.
A Approximation of data on corner panels for re-solve Here we give explicit bounds for the rate of convergence of the contribution of the rest of the boundary to a corner panel. In particular, given a polygonal domain with boundary Γ, let x denote a vertex of Γ and C = Γ ∩ B r (x), where B r (x) is the ball of radius of r centered at x. We choose r so that Γ ∩ B 2r (x) corresponds to a wedge with internal angle πα and side lengths 2r.
Theorem 9. Let Γ be the boundary of a polygon and x be a vertex. Let r > 0 be a real number such that Γ ∩ B 2r (x) corresponds to a wedge with internal angle πα and side lengths 2r, where B R (x) denotes a ball of radius R centered at x. Let L denote the length of Γ and γ : [−L/2, L/2] → Γ be an arclength counterclockwise parameterization of Γ such that γ(0) = x. Finally, for any f ∈ L 2 (Γ), let H : [0, r] → R be the function defined by Then H is analytic in a neighborhood of 0 with Taylor series coefficients {a n } satisfying Proof. Without loss of generality we can assume that Γ is shifted, oriented and parameterized so that x = γ(0) = 0 and the leg of the wedge corresponding to positive t is oriented along the positive x axis. Then Since γ(s) − γ(t) > 2r it follows that In particular, H has a Taylor series about t = 0, B Strong approximation of density away from corner panels In this section, we demonstrate that for a panel which is sufficiently far from the corner and discretized using M Gauss-Legendre nodes, the density computed using the adjoint of a Dirichlet discretization can be interpolated accurately at any point on the panel. As before let Ω denote a polygonal domain with boundary Γ. Let L denote the length of the boundary, and let γ : [0, L] → R 2 denote an arc-length parameterization of the boundary. Assume that the discretization satisfies the following assumptions: 1. All panels which are not at a corner, are separated from the closest corner by at least their panel length.
2. Let E ρ (Γ i ) denote the Bernstein ρ−ellipse (see [21]) corresponding to the panel Γ i , and let ρ i be such that E ρ i (Γ i ) does not intersect Γ \ S Γ i , where S Γ i is the edge containing Γ i (see fig. 6). Let ρ 0 = min i ρ i > 1.
3. All panels which are not adjacent to a panel at the corner, are separated from the corner by 2r c where r c denotes the length of the panel at the corner.
Under these assumptions, it can be shown that the accuracy of computing the Legendre coefficients of the density (at panels which are not adjacent to a corner panel) for the Neumann problem using the adjoint discretization is related to the accuracy in the computation of the solution to an associated Dirichlet problem. Let A denote the operator corresponding to interior Dirichlet problem using a double layer potential. Let g denote the right hand side for the Neumann problem, let σ denote the corresponding solution. Let f be a Legendre polynomial of degree n scaled to the panel γ([s 1 , s 2 ]) and 0 everywhere else. Then where σ f is the solution of the interior Dirichlet problem with boundary data f using a double layer potential. Using lemma 1, the above statement implies that the error in computing the M Legendre coefficients of the density for the Neumann problem is the same as the error in computing the solution of a Dirichlet problem with data given by a Legendre polynomial on the same panel.
There are two concerns which must be addressed. First, the accuracy of computing f (s) for any point s ∈ [0, L] using an M point Gauss-Legendre quadrature on [s 1 , s 2 ], and secondly, the resolution of the functionf (s) on the given discretization of the boundary. For any s which is contained on the same segment as γ([s 1 , s 2 ]), the kernel k(s, t) is identically 0. Thus the boundary dataf (s) = 0 on the same edge as the panel γ([s 1 , s 2 ]). We further observe that f (t) is an entire function when extended to the complex plane, since it is a Legendre polynomial. Moreover, the nearest singularity of the function k(s, t) in the complex plane as a function of t is at γ(s). From assumption 2, it follows that the error in computingf (s) using an M point Gauss-Legendre rule is bounded by Cρ −M 0 , where the constant C is related to the smoothness of k(s, t) as a function of t ∈ C. Thus the functionf (s) can be computed to any desired precision by increasing the order of quadrature nodes used to compute the integrals.
With regards to the resolution of the of the functionf (s) on the given discretization of the boundary, we note that the closest singularity of the functionf (s) when restricted to a panel away from the corner and not on the same edge as γ([s 1 , s 2 ]) is the closest point on the panel γ([s 1 , s 2 ]). However, by assumption 2, the error in resolving the functionf (s) using an M point Gauss-Legendre or Chebyshev panel is bounded by Cρ −M 0 . Note that the behavior off (s) in the complex plane is related to the behavior of k(s, t) in the complex s plane and hence the constant C is O(1). For the panels, at the corner, based on the proof in Appendix A, the error in resolving the functionf (s) when truncated to a Taylor series of order N is less than C2 −N , since all points on γ([s 1 , s 2 ]) are well-separated from corners by twice the panel length r c . Thus, by making the panels small enough, ρ can be increased arbitrarily to obtain desired tolerances on the boundary dataf (s) on the corresponding discretization of the boundary.
Thus, the boundary dataf (s) is piecewise analytic , which can be approximated to any desired tolerance by appropriately reducing the panel sizes. This is the precise setup for which the discretization of the Dirichlet problem is designed to obtain accurate solutions to the densityσ.