A numerical technique for linear elliptic partial differential equations in polygonal domains

Integral representations for the solution of linear elliptic partial differential equations (PDEs) can be obtained using Green's theorem. However, these representations involve both the Dirichlet and the Neumann values on the boundary, and for a well-posed boundary-value problem (BVPs) one of these functions is unknown. A new transform method for solving BVPs for linear and integrable nonlinear PDEs usually referred to as the unified transform (or the Fokas transform) was introduced by the second author in the late Nineties. For linear elliptic PDEs, this method can be considered as the analogue of Green's function approach but now it is formulated in the complex Fourier plane instead of the physical plane. It employs two global relations also formulated in the Fourier plane which couple the Dirichlet and the Neumann boundary values. These relations can be used to characterize the unknown boundary values in terms of the given boundary data, yielding an elegant approach for determining the Dirichlet to Neumann map. The numerical implementation of the unified transform can be considered as the counterpart in the Fourier plane of the well-known boundary integral method which is formulated in the physical plane. For this implementation, one must choose (i) a suitable basis for expanding the unknown functions and (ii) an appropriate set of complex values, which we refer to as collocation points, at which to evaluate the global relations. Here, by employing a variety of examples we present simple guidelines of how the above choices can be made. Furthermore, we provide concrete rules for choosing the collocation points so that the condition number of the matrix of the associated linear system remains low.

Integral representations for the solution of linear elliptic partial differential equations (PDEs) can be obtained using Green's theorem. However, these representations involve both the Dirichlet and the Neumann values on the boundary, and for a well-posed boundary-value problem (BVPs) one of these functions is unknown. A new transform method for solving BVPs for linear and integrable nonlinear PDEs usually referred to as the unified transform (or the Fokas transform) was introduced by the second author in the late Nineties. For linear elliptic PDEs, this method can be considered as the analogue of Green's function approach but now it is formulated in the complex Fourier plane instead of the physical plane. It employs two global relations also formulated in the Fourier plane which couple the Dirichlet and the Neumann boundary values. These relations can be used to characterize the unknown boundary values in terms of the given boundary data, yielding an elegant approach for determining the Dirichlet to Neumann map. The numerical implementation of the unified transform can be considered as the counterpart in the Fourier plane of the well-known boundary integral method which is formulated in the physical plane. For this implementation, one must choose (i) a suitable basis for expanding the unknown functions and (ii) an appropriate set of complex values, which we refer to as collocation points, at which to evaluate the global relations. Here, by employing a variety of 2015 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/ by/4.0/, which permits unrestricted use, provided the original author and source are credited.

Introduction
A new method for analysing boundary-value problems (BVP) for linear and for integrable nonlinear partial differential equations (PDEs) was introduced by the second author in the late Nineties [1][2][3]. This method, which is usually referred to as the unified transform (or the Fokas transform), has been applied to a variety of linear elliptic PDEs formulated in the interior of a polygon. Important results in this direction include the following: (i) for the Laplace, modified Helmholtz and Helmholtz equations, it is possible to express the solution in terms of integrals in the complex λ-plane (complex Fourier plane). These integrals contain certain integral transforms of the Dirichlet and of the Neumann values on the boundary of the polygon. Hence, these integral formulae provide the analogue of the classical Green's representations, but now the formulation takes place in the complex Fourier plane instead of the physical plane. (ii) The above transforms of the Dirichlet and of the Neumann boundary values are related via two simple algebraic equations called global relations. These relations provide a characterization of the generalized Dirichlet to Neumann map. (iii) By employing the integral representation and global relations mentioned in (i) and (ii), respectively, it has been possible to obtain exact solutions for a variety of problems for which apparently the usual approaches fail, see e.g. [4,5]. (iv) Ashton [6,7] has developed a rigorous approach for deriving well posedness results for linear elliptic PDEs using the new formalism. This includes the analysis of BVPs with distributional data and with corner singularities [8]. (v) The new method can be applied to linear PDEs with nonlinear boundary conditions, see e.g. [9][10][11]. (vi) The first steps have been taken towards extending the unified transform to three dimensions, see e.g. [9,12].
The analysis of the global relations yields a novel numerical technique for the numerical solution of the generalized Dirichlet to Neumann map, i.e. for the determination of the unknown boundary values in terms of the prescribed boundary data [13][14][15][16][17][18][19][20][21]. Substantial progress in this direction was made by Fornberg and co-worker [14,15]. The global relations couple the finite Fourier transform of the given boundary data with the finite Fourier transform of the unknown boundary values. For the determination of these boundary values one has to (a) choose appropriate basis functions and (b) suitable collocation points in the Fourier plane. For the numerical computation of the finite Fourier transforms of the basis functions, Fornberg and co-worker have used Legendre polynomials, and have employed the fact that the finite Fourier transform of the Legendre polynomials can be expressed in terms of the modified Bessel function of order half integer. In the above-mentioned papers, the authors have also used the so-called Halton nodes for collocation points, and have used the crucial observation that the conditioning of the associated linear system improves if the linear system is overdetermined. Here, following Fornberg and co-worker, we also use Legendre polynomials and also overdetermine the relevant system. However, motivated by the results of [21], we introduce a new choice of collocation points. In this way, we are able to improve dramatically the relevant condition number. For example, for a particular BVP considered in [15], the condition number improves from approximately 10 16 to 4.9.
There exist several different numerical methods available for solving linear elliptic PDEs, which include finite-difference methods, finite-element methods, boundary-integral equations and the method of particular solutions. The main novelty of the method developed in this paper, compared with standard methods, is that it is a boundary-based discretization that does not involve the computation of singular integrals (as opposed to the discretizations of boundary rspa.royalsocietypublishing.org Proc. R. Soc. integral equations). Similar ideas have been proposed before in the literature, and the relationship of the method here to these earlier approaches is discussed in detail in §6.
This paper is organized as follows: in §2, we review the concept of global relations and obtain these relations for the Laplace, modified Helmholtz and Helmholtz equations in the interior of a polygonal domain. In §3, we approximate the known and unknown boundary values in terms of Legendre polynomials, and derive the approximate global relations. In §4, we discuss a convenient choice for collocation points for the particular case of a convex polygon. In §5, we present several numerical calculations comparing the unified transform solution with an exact solution or a numerical solution obtained via the finite-element method. Finally, in §6, we discuss further these results.

Global relations
where ∂Ω denotes the boundary of the domain Ω, v is a solution of the Laplace equation, and ∂/∂N denotes the derivative along the outward normal direction of the boundary. Letting For the modified Helmholtz equation where k > 0 is the wavenumber, we choose the particular solution v = e (ik/2)(z/λ−λz) , λ ∈ C \ {0}, and then equation (2.2) yields the global relation For the Helmholtz equation

(a) A polygonal domain
In what follows, we consider the particular case where Ω is the interior of a polygon with n sides and corners {z j } n 1 . The jth side of this polygon, which is the side between the edges z j and z j+1 , can be parametrized by the expression where m j and h j are defined by (2.11) Using the above parametrization in equation (2.5) and employing ds = |h| dt, we find that the global relation for the Laplace equation in a polygonal domain becomes Typical boundary conditions for the jth side are given below: and Robin: In the above equations, g j is a given function and γ j is a given constant. If any of the above boundary conditions is prescribed on each side of the polygon, it follows that the global relation (2.12) involves n unknown functions. For the particular case that a Dirichlet boundary condition is prescribed on every side, it is rigorously shown in [6,7], that equation (2.12) together with the Schwartz conjugate of (2.12) uniquely determine the unknown values {∂u j /∂N } n 1 . In what follows, we present a simple procedure for computing numerically the n unknown boundary values for the Laplace, modified Helmholtz and Helmholtz equations in the case that any of the boundary conditions (2.13)-(2.15) is prescribed on the jth side of the polygon. We emphasize that it is possible to prescribe different types of boundary conditions on different sides.

The approximate global relation
We expand the functions {u j } n 1 and {∂u j /∂N } n 1 in terms of N ∈ Z + basis functions denoted by Substituting the expressions (3.1) into the global relation (2.12) and using equation (3.2), we find the following approximate global relation for the Laplace equation: .
Finally, the approximate global relation for Helmholtz equation is given by Recalling that equation (  Following Fornberg and co-worker [15], we let S l (t) = P l (t), where P l (t) denotes the Legendre polynomials of order l. It is well known that the Fourier transform of P l (t) can be expressed in terms of the spherical Bessel function J l [22], which in turn can be expressed in terms of the modified Bessel function I l+1/2 . Also, an analytical expression for the Fourier transform of P l (t) is given in [23]. Thus, the expressionsŜ l which appear in (3.3)-(3.5) can be computed by

Collocation points
In order to motivate a convenient choice for the collocation points, we concentrate on the particular side with index p. The global relation (2.12) can be written in the form where F j (λ, t) is given by Thus, for the particular side with index p, a natural choice for λ is h p λ = real number. However, it is also desirable to choose this real number in such a way that This condition guarantees that as λ → ∞, the contribution from the other sides tends to zero. It is shown in [21] that for a convex polygon the condition (4.3) can indeed be satisfied provided that the above number is negative. Thus, for the side p we choose λ by Thus, this can be written in the form of the second of equations (4.4) provided thatf > k|h p |. On the other hand, if k|h p | >f , then The above analysis shows that for a convex polygon an appropriate choice of the collocation points associated with the side p is λ = −h p f , f > 0. We find it convenient to write f in the form below. Thus, associated with the side p we choose the following M collocation points

Numerical results
In this section, we present the results of a parametric study of the key parameters in the unified transform, namely (R, n, N, M). Our chosen convex domains are the octagon shown in figure 1 and the trapezoid shown in figure 2. We compare the solution obtained via the unified transform with an analytical solution and the solution computed via the finite-element method.

(a) The octagon geometry
For the octagon geometry shown in figure 1, we compare an analytical solution to the solution computed via the above approach. We denote the analytic solution by u A (x, y) and the solution obtained via the unified transform by u UT (x, y).
We choose the same analytical solution as the one used in [15]. The parametrization error for this particular solution has been studied in some detail by [15]. The condition number K versus  To validate the unified transform solution, we compare it with the solution obtained using the finite-element method [24]. We have used linear finite elements and the finite-element mesh of the trapezoidal geometry involves a total of 1954 nodes, 3715 triangles and the boundary has 191 edges. This comparison is shown in figure 6. The results show an excellent agreement between the finite element method and the unified transform.

Conclusion
The so-called global relations play a crucial role in the construction of analytical solutions for both evolution and elliptic PDEs. For the Laplace, modified Helmholtz and Helmholtz equations, the global relations are equations (2.5), (2.7) and (2.9), respectively, as well as the equations obtained from these equations by taking the complex conjugate and then replacingλ by λ. By employing the fact that λ is arbitrary, the global relations provide an elegant characterization of the Dirichlet to Neumann map.
It has been realized since the work of [16] that the global relations can be solved numerically. In this direction, different numerical techniques have been derived by several authors, see e.g. [13][14][15][16][17][18][19][20][21]. Here, following Fornberg and co-worker, we use the Legendre basis and also overdetermine the relevant system; furthermore motivated by the results of [21], we introduce a simple choice of 'collocation' points, see equation (4.8). This choice involves the positive number R and the positive integer M; M is a measure of overdeterminancy and R/M is the distance between two consecutive points.
We provide strong numerical evidence that if M and R satisfy the constraints given by equation (4.9) then the condition number of the associated system is of order 1. For example, for the trapezoidal domain analysed in [14,15] the condition number reduces from O(10 8 ) to 4.7.
We next discuss the relation of the method presented in this paper with two other methods for solving linear elliptic PDEs in the literature which are also based on the relation (2.2): (i) The null field method for the solution of the Helmholtz equation in the exterior of a bounded obstacle (originally introduced by Waterman in [25,26]) and (ii) a method for the solution of the Helmholtz equation above a periodic rough surface introduced by DeSanto [27] and further developed by DeSanto and co-workers [28][29][30][31]. The advantage of these two methods, as well as of the method described in this paper, is that they are boundary-based discretizations that do not involve the computation of singular integrals (as opposed to the discretizations of boundary integral equations). Relations between these methods are discussed in detail in [32, §10] and [33, §4]. In the null-field method, u in (2.2) is the solution of the Helmholtz equation in the exterior of a bounded obstacle, and v is one of a countable family of separable solutions of the Helmholtz equation in polar coordinates that satisfies the appropriate radiation condition (see, e.g. [34, §7.5]). The main difference between the null-field method and the method in this paper are the following: (i) the former method is used for an exterior of an obstacle, whereas the current method is used for the interior of a polygon and (ii) in the null-field method the unknown boundary value is expanded in a global basis (i.e. one in which the support of the basis functions is the whole of ∂Ω), whereas the method in this paper uses local bases (where each basis function is supported only on one side of the polygon). Regarding the method of DeSanto, u in (2.2) is the solution of the Helmholtz equation above a periodic rough surface, and v is chosen to be one of a countable family of separable solutions of the Helmholtz equation in Cartesian coordinates (so-called generalized plane waves) that satisfies the appropriate radiation condition. This method has been implemented with a variety of different bases chosen to express the unknown boundary value. Indeed, the authors of [30] proposed two different bases for the unknown boundary value to be expanded in: a local basis consisting of piecewise-constant basis functions, and a global basis consisting of the traces of the functions v in (2.2), and the authors of [31]  Ethics statement. In this study, we have not conducted any experiments on human beings or animals. Furthermore, we are not using any data collected from human beings or animals by other research groups or organizations.
Data accessibility. We have not used any data (measurements) collected by any organization including the University of Cambridge. This work involves a clear and rigorous description of a numerical implementation of the Fokas transform for solution of linear elliptic PDEs for convex polygons. The described method does not require or rely on any data.
Author contributions. A.S.F. invented and proposed the unified transform (Fokas transform) in the Nineties and derived the global relations presented in this manuscript. S.A.S. had worked on an appropriate choice of the collocation points. P.H. continued the investigation of the optimal choice of the collocation points and finalized through a rigorous numerical study the optimal degree of overdetermination for the system of linear equations. By means of extensive numerical tests, he obtained the empirical rules for the optimal choice of the key parameters in the Fokas transform. These rules guarantee a low condition number of the system of linear equations. Furthermore, he provided the initial draft of this manuscript as well as all the numerical results presented in this manuscript.