## Abstract

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.

## 1. 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–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–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–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 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.

## 2. Global relations

For completeness, we first recall the global relations for the Laplace, modified Helmholtz and Helmholtz equations. Solutions of the Laplace equation in the closed domain *Ω* denotes the boundary of the domain *Ω*, *v* is a solution of the Laplace equation, and *k*>0 is the wavenumber, we choose the particular solution

For the Helmholtz equation
*u* real, a second global relation can be obtained from the global equations (2.5), (2.6), and (2.9) via Schwartz conjugation, i.e. by taking the complex conjugate of these equations and then replacing

### (a) A polygonal domain

In what follows, we consider the particular case where *Ω* is the interior of a polygon with *n* sides and corners *j*th side of this polygon, which is the side between the edges *z*_{j} and *z*_{j+1}, can be parametrized by the expression
*m*_{j} and *h*_{j} are defined by
*s*=|*h*| d*t*, we find that the global relation for the Laplace equation in a polygonal domain becomes
*jth* side are given below:

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

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 *j*th side of the polygon. We emphasize that it is possible to prescribe different types of boundary conditions on different sides.

## 3. The approximate global relation

We expand the functions *j* there exists a linear relation between

Let

Recalling that equation (2.12) implies a linear relation between *nN* unknown constants. However, it is crucial to observe that λ is an *arbitrary* complex number. Thus by evaluating equation (3.3) and its Schwartz conjugate at the

### (a) The choice of the basis functions

Following Fornberg and co-worker [15], we let *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

## 4. 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
*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
*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
*h*_{p}.

Similarly for the modified Helmholtz equation, the analogue of the first of equations (4.4) is the equation

The above analysis shows that for a *convex* polygon an appropriate choice of the collocation points associated with the side *p* is *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

In this way we choose *nM* collocation points, and hence by evaluating the two global relations at these points we obtain 2*nM* equations for *nN* unknowns. Hence, *M* must satisfy the constraint *M*≥*N*/2. Clearly, *M* specifies the ‘overdeterminedness’ of the linear system (the larger the *M* the larger the overdeterminedness). On the other hand, the parameter *R*/*M* defines the distance between two consecutive collocation points.

The error of approximating the Dirichlet and Neumann boundary values with the expressions (3.1) depends on *N*, whereas the condition number *quadruple* (*R*,*n*,*N*,*M*). By scrutinizing a variety of numerical experiments, we have found the following necessary conditions for low

## 5. 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 *N* on each side of the octagon for different choices of (*M*,*R*) is shown in figure 3. In the case of the Helmholtz and modified Helmholtz equations, we have fixed the wave number at *M*=*nN*, *R*/*M*=2 results in a low condition number.

The comparison of the analytical solution and the solution obtained via the unified transform method for one side of the octagon is shown in figure 4. It can be observed that the unified transform provides an accurate solution. Furthermore, the choice of *M*=*nN*, *R*/*M*=2 results in the low condition number of

### (b) The trapezoidal geometry

We now consider a case for which we cannot obtain an analytical solution. Consider the trapezoidal domain shown in figure 2. We choose a Robin boundary condition on side 1 and homogeneous Neumann boundary condition on the remaining sides:
*M*=*nN*, *R*/*M*=2 ensures a low condition number.

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.

## 6. 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

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–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–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] proposed a variant of the second basis, where the complex-conjugates of the functions *v*, instead of the functions themselves, are used. An important point to note is that in both the null-field method and the method of DeSanto there is *no* flexibility in the choice of *v*, and thus there is *no* analogue of the question of how to choose λ in the global relations.

## Disclaimer

The manuscript has been authored by employees of the Department of Applied Mathematics and Theoretical Physics at the University of Cambridge.

## 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.

## Funding statement

This research has been funded by the EPSRC grant, EP/H04261X/1: ‘Integrability in Multidimensions and Boundary Value Problems.’

## Conflict of interests

We have no conflicts of interest.

- Received October 20, 2014.
- Accepted January 14, 2015.

© 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.