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.

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

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 , satisfies Green's second identity where ∂Ω denotes the boundary of the domain Ω, v is a solution of the Laplace equation, and denotes the derivative along the outward normal direction of the boundary. Letting noting that , is a particular solution of the Laplace equation, and using the identity equation (2.2) yields the global relation For the modified Helmholtz equation where k>0 is the wavenumber, we choose the particular solution , and then equation (2.2) yields the global relation

For the Helmholtz equation we choose the particular solution , and then equation (2.2) yields the global relation For the Laplace, modified Helmholtz and Helmholtz equations with 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 with λ.

A polygonal domain

In what follows, we consider the particular case where Ω is the interior of a polygon with n sides and corners . The jth side of this polygon, which is the side between the edges z and z, can be parametrized by the expression where m and h are defined by 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:

In the above equations, g is a given function and γ 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 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 and in terms of basis functions denoted by : where and are constants. Equations (2.13)–(2.15) imply that for each j there exists a linear relation between and .

Let denote the Fourier transform of , i.e. 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: Similarly the approximate global relation for the modified Helmholtz equation is given by Finally, the approximate global relation for Helmholtz equation is given by

Recalling that equation (2.12) implies a linear relation between and , it follows that equation (3.3) together with its Schwartz conjugate are two equations with 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 points where we can determine the unknown constants. Similar considerations are valid for the modified Helmholtz and Helmholtz equations.

The choice of the basis functions

Following Fornberg and co-worker [15], we let , where P(t) denotes the Legendre polynomials of order l. It is well known that the Fourier transform of P(t) can be expressed in terms of the spherical Bessel function J [22], which in turn can be expressed in terms of the modified Bessel function I. Also, an analytical expression for the Fourier transform of P(t) is given in [23]. Thus, the expressions 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(λ,t) is given by

Thus, for the particular side with index p, a natural choice for λ is hλ=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 where is an arbitrary positive constant, and denotes the complex conjugate of h.

Similarly for the modified Helmholtz equation, the analogue of the first of equations (4.4) is the equation which can also be written in the form of the second of equations (4.4). For the Helmholtz equation, the analogue of equation (4.5) is the equation Thus, this can be written in the form of the second of equations (4.4) provided that . On the other hand, if , then

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 2nM 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 of the relevant linear system depends on the quadruple (R,n,N,M). By scrutinizing a variety of numerical experiments, we have found the following necessary conditions for low :

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.

Figure 1.

An octagon (typical polygon). The nodes are labelled counter clockwise. (Online version in colour.)

Figure 2.

The trapezoidal geometry. (Online version in colour.)

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 uA(x,y) and the solution obtained via the unified transform by uUT(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 versus the number of basis functions 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 . It can be observed in all three cases, namely the Laplace, the modified Helmholtz and the Helmholtz equations, the choice of M=nN, R/M=2 results in a low condition number.

Figure 3.

The condition number for the octagon geometry. The blue solid line is M=nN/2, R=2M; the red solid line with dots is M=nN,R=2M and the solid black line with dots is M=2nN,R=2M. For the modified Helmholtz and Helmholtz equations . N, number of basis functions on each side. (a) Condition number for Laplace equation; (b) condition number for modified Helmholtz equation and (c) condition number for Helmholtz equation.

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 .

Figure 4.

The solution and the corresponding absolute error on the first side (a typical side) of the octagon geometry for the modified Helmholtz equation. Subplot (a) shows the solution, where the blue solid line with dots is the analytical solution. The red dashed line with the circle is the unified transform solution. Subplot (b) shows the absolute error i.e. . The wavenumber is . The parameters associated with the unified transform solution are n=8, N=20, M=160 and R=320. The condition number of the system matrix is .

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: The results of the parametric study for the trapezoid is shown in figure 5. The results reconfirm the parametric study of the octagon and reveal that the choice M=nN, R/M=2 ensures a low condition number.

Figure 5.

The condition number for the trapezoidal geometry. The blue solid line is M=nN/2, R=2M; the red solid line with dots is M=nN, R=2M and the solid black line with dots is M=2nN, R=2M. For the modified Helmholtz and Helmholtz equations N, number of basis functions on each side. (a) Condition number for Laplace equation; (b) condition number for modified Helmholtz equation and (c) condition number for Helmholtz equation.

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.

Figure 6.

The solution u(t) (s=1,…,4) on each side of the trapezoid for the modified Helmholtz equation. On the plots in the first column, the solid blue line with dots is the solution using the unified transform, whereas the red line with circles is the solution using the finite-element method. The figures on the second column show the corresponding absolute difference. The parameters associated with the unified transform are . Using these parameters the condition number of the system matrix is .

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-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(108) 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.