## Abstract

A recently discovered transform approach allows a large class of PDEs to be solved in terms of boundary and/or contour integrals. We introduce here a spectrally accurate numerical discretization of this approach for the case of Laplace's equation on a polygonal domain, and compare it against an also spectrally accurate implementation of the traditional boundary integral formulation.

## 1. Introduction

About a decade ago, Fokas and collaborators introduced a new analytical transform technique that provides closed-form solutions (in terms of contour integrals) for a variety of PDE boundary value (BV) and initial-boundary value (IBV) problems (Fokas 1997, 2001, 2002, Fokas & Pelloni 2005; Spence 2010). For several dispersive IBV problems, numerical implementations of these integrals have been found to offer opportunities that in some respects improve on previously available techniques (Flyer & Fokas 2008; Fokas 2008; Fokas *et al.* 2009). In the case of the Dirichlet-to-Neumann map problem for elliptic equations, the Fokas approach provides a direct alternative to regular boundary integral (BI) methods in the sense that it also requires discretization only along the domain boundary, and not throughout the interior of the domain. The new formulation provides a spectral (Fourier) space analogue of the usual BI method, which is formulated in the physical space.

This present study will focus on simple polygonal domains and, in particular, on the trapezoidal domain displayed in figure 1. We introduce here a spectrally accurate implementation (assuming the solutions themselves are sufficiently smooth), i.e. the relevant errors decrease as *O*(e^{−cN}) where *N* is the number of parameters used for discretization along each side, and *c*>0 is some constant. Thus, low values for *N* often suffice for high accuracy. Most previous numerical implementations of this new analytical formulation have been limited to a relatively low *O*(*N*^{−2}) algebraic rate of convergence (Fulton *et al.* 2004; Fokas 2008; Fokas *et al.* 2009), thus using relatively large *N*-values and iterative solution procedures. The possibility of spectral convergence (by discretizing with Chebyshev rather than Fourier expansions along the boundaries) was noted in Sifalakis *et al.* (2008) and Smitheman *et al.* (2010). However, this approach encountered conditioning and convergence issues in cases of irregular polygons.

The goals of the present study are to show:

— Legendre (rather than Chebyshev) expansions allow key integrals to be evaluated in closed form,

— a least-squares approach provides good numerical conditioning,

— no iterative linear algebra algorithms are needed,

— corner singularity corrections can be successfully incorporated in the procedure.

Brief numerical comparisons are provided between the new and the traditional BI methods, denoted by BI-F and BI-T respectively, and also with a direct collocation method based on sets of Laplace solutions.

We recall that BI-T implementations (as surveyed in Pozrikidis (2002)), often encounter difficulties with certain quadrature evaluations. As a result, it is common to discretize the boundary into a large number of *boundary elements*, generally leading to algebraic convergence rates. We use here for our comparisons an implementation of the BI-T formulation which also features spectral accuracy. It should be noted that conformal mapping based methods can be very effective for the Laplace equation on polygonal domains (Driscoll & Trefethen 2002). However, we do not consider this approach since, in contrast to BI methods, it cannot be generalized to other types of PDEs.

## 2. Governing equation and three analytical formulations

We consider Laplace's equation
2.1over the bounded two-dimensional domain *Ω* shown in figure 1. Either *q*(*x*,*y*) or *q*_{n}(*x*,*y*) (normal derivative in the outward direction) is specified along each section of the domain boundary ∂*Ω*. This paper focuses on the Dirichlet-to-Neumann map problem, which amounts to finding, along each section, the boundary information that was not specified. We describe next three different analytical formulations which all lead to practical numerical algorithms. None of them will require any discretization within *Ω* to solve the problem.

### (a) New boundary integral formulation

The first analytical formulation that we consider is based on the works of Fokas cited above (hence the acronym BI-F). The Dirichlet-to-Neumann map is described by the following pair of integral relations that must hold when integrating around the domain boundary
2.2
2.3The integrals are taken in the positive direction, and *s* is the arc length. The complex variable *z*=*x*+*iy* denotes locations in the (*x*,*y*)-plane. Complex conjugation is defined in the usual way, i.e. . The equations (2.2) and (2.3) are called the *global relations*. They hold *for any value* of the complex variable *k*. By choosing different values for *k*, they provide an unlimited number of independent relations between *q*(*s*) and *q*_{n}(*s*).

We are in this study not concerned with calculating *q*(*x*,*y*) at locations in the domain interior. This can however be achieved by evaluating a succession of integrals, as described in Fokas (2008, eqns (11.3)–(11.5)).

### (b) Traditional boundary integral formulation

The second analytical formulation we consider is the BI-T method. In this case, the key relation linking *q*(** ξ**) and

*q*

_{n}(

**) around the boundary is 2.4where 2.5is the infinite domain Green's function and**

*ξ**G*

_{n}(

**,**

*x***) is its outward pointing normal derivative. This function**

*ξ**G*(

**,**

*x***) satisfies equation (2.1) for**

*ξ***≠**

*x***where**

*ξ***is an abbreviated notation for (**

*x**x*,

*y*). The variable

*s*represents again the arc length, and

**=**

*ξ***(**

*ξ**s*) traces out the boundary ∂

*Ω*. Since

**can be any point along the boundary ∂**

*x**Ω*, we again obtain an unlimited number of independent relations between

*q*(

**) and**

*ξ**q*

_{n}(

**).**

*ξ*If, in addition to solving the Dirichlet-to-Neumann problem, one wants to obtain the solution *q*(** x**) throughout the inside of

*Ω*, this can, in the BI-T formulation, be achieved very simply by replacing (1/2)

*q*(

**)|**

*x*_{x∈∂Ω}with

*q*(

**)|**

*x*_{x∈Ω}in equation (2.4), i.e. 2.6The difference by a factor of two stems from the fact that equation (2.4), with

**∈∂**

*x**Ω*, should be interpreted as a principal value integral. It therefore picks up only half of the contribution from the singularity in

*G*

_{n}(

**,**

*x***)**

*ξ**q*(

**) at the location**

*ξ***=**

*ξ***compared to the case when**

*x***is inside**

*x**Ω*.

If one adheres to the (unusual) convention of assigning the value of a principal value integral to the limit when reaching the boundary from the inside, it will hold that (cf. §3*b*), and the relation (2.6) can then, for our four-sided polygonal domain, be written as:
2.7where
This formulation is very convenient because, as soon as *q*^{(ν)}(** ξ**) and are expanded in Legendre polynomials, all the integrals in equation (2.7) are available in closed form (appendix B). Also, it directly leads to a linear system for the unknown Legendre coefficients, as illustrated in figure 4.

The present BI-T implementation is not the only version that provides spectral accuracy. We prefer it here because of it similarity to the BI-F discretization (in terms of Legendre expansions), its ease of implementation and the freedom from numerical quadrature.

### (c) Collocation with sets of Laplace solutions

It is straightforward to find infinite sets of solutions to Laplace's equation. The collocation with sets of Laplace solutions (the CLS approach) amounts to combining such solutions, so that the specified boundary conditions (BCs) also become satisfied. Examples of Laplace solutions include the real and/or the imaginary parts of any set of analytical functions, such as *z*^{μ} or *e*^{iμz}, *μ* integer. It is often advantageous to consider solution sets that include functions with certain corner properties. The analytical forms of the corner singularities are described in Lehman (1959) and Wigley (1964), with some numerical approaches for handling them noted in Li & Lu (2000). With regard to our specific test domain (using the notation in figure 1), we note that the functions
2.8satisfy Laplace's equation
for *r*>0 as well as the BCs
2.9Hence, the same will hold for all linear combinations of these functions
2.10Assuming that we are only interested in bounded solutions, we include in the sum only terms with *μ*≥1. The functions in this sum are known as stress (or flux) intensity factors, with the coefficient for the *μ*=1 function often of particular interest.

For the remaining three corners (with interior angles 45^{°} and 90^{°}), singularities will arise only if the BC are inconsistent with the PDE (a situation encountered in Test case 3) or are inconsistent with differentiated versions of the PDE. For example, the expressions
2.11satisfy Laplace's equation Δ*h*_{μ}=0 for *r*>0 together with *h*_{μ}(*r*,0)=0 and *h*_{μ}(*r*,*π*/2)=*r*^{2μ}, i.e. they violate Δ^{μ}*h*_{μ}≡(∂^{2}/∂*x*^{2}+∂^{2}/∂*y*^{2})^{μ}*h*_{μ}=0, *μ*>1, at *r*=0. When such violations occur, there exist logarithmic singularities at the corner, as indicated by equation (2.11). In the worst locally bounded case (*μ*=1), this type of singularity will be of the form , i.e. much milder than the singularities that are typical for when sides with Dirichlet and Neumann conditions meet at a 135^{°} corner. We recall that corner singularities can also arise for IBV problems, often with more complex singularity structures (Flyer & Fornberg 2003).

## 3. Numerical implementation strategies

### (a) BI-F formulation

If the unknown parts of *q*(*s*) and *q*_{n}(*s*) in equations (2.2) and (2.3) are discretized by *N* parameters along each of the four sides, we need to use *K*≥2*N* different complex values for *k*, in order to get a sufficient number of equations for the unknowns parameters. Using Legendre expansions for these discretizations leads to spectral accuracy and also to key integrals of a form that can be evaluated analytically. The first step in this process is to parametrize each side of the domain using a new variable *t*∈[−1,1]. The equations (2.2) and (2.3) then become
3.1and
3.2For the present domain geometry, we obtain
3.3For each of these *k*-values, chosen as discussed in appendix A, we evaluate (analytically if possible, else numerically, for example, by adaptive Gaussian quadrature or via conversion to Legendre expansions), the integrals in equations (3.1) and (3.2) for which either *q* or *q*_{n} have been specified. The unknown function along each of the four sides *ν*=1,…,4 will be represented as Legendre expansions with unknown coefficients
3.4Substituting the expansion (3.4) into equations (3.1) and (3.2) leads to integrals of the general form
3.5Here, *α* is the coefficient that appears in front of the *t*-term in the exponent −*ikz* in case of equation (3.1) and in case of equation (3.2), when using the appropriate formula for *z* (dependent on the side number according to the leftmost column in equation (3.3)). The integral stated in equation (3.5) is an entire function of *α*. We have written the factor in front of the Bessel *I* function as instead of , in order for the branch cuts along the negative real axis for the square root and for the Bessel *I* function to cancel out when using the standard branch cut conventions for these two functions.

Each of the four domain sides contains *N* unknowns (, *m*=0,1,…,*N*−1). Hence equations (3.1) and (3.2) take the form of 2*K* linear equations for these 4*N* expansion coefficients . As noted above, one needs to choose *K*≥2*N*. Numerical tests show that larger values of *K* generally gives more robust and accurate results, with *K*=5*N* a good routine choice. The general structure of the overdetermined linear system is illustrated in figure 2. A standard linear systems solver (such as the ‘\’-operator in Matlab) provides the unknown expansion coefficients. The corresponding boundary functions then follow from equation (3.4).

### (b) BI-T formulation

The relation (2.7) will hold for every value ** x** that we choose on ∂

*Ω*. As a first step, we again change the parametrization separately for each side from arc length

*s*to

*s*=

*s*(

*t*) with

*t*∈[−1,1]. Relative to this interval in a

*t*-plane, the sample point

**=(**

*x**x*,

*y*) will be at some location (

*α*,

*β*)—illustrated in figure 3 for the case that the sample point is located on a different side than the one we just now are integrating along. Along the

*t*-axis, equation (2.5) will take the form where

*c*is some constant that follows from the variable change that brought the side to

*t*∈[−1,1]. Then Hence, the integrals that arise from equation (2.7) will become of the general forms 3.6and 3.7respectively, where

*f*(

*t*) will either be known or unknown functions. If

**is located on the side that is here parametrized by**

*x**t*∈[−1,1], the factor

*β*/((

*α*−

*t*)

^{2}+

*β*

^{2}) in the integral (3.7) approaches

*πδ*(

*t*−

*α*) as

*β*↘0. Most numerical quadrature methods face difficulties both for

*β*=0 (when the integral (3.6) features a log-type singularity at

*t*=

*α*), and for small non-zero values of

*β*.

Legendre expansions are again well suited for discretization of the *q* and the *q*_{n}-functions along sides *ν*=1,…,4:
3.8and
3.9Substituting such expansions for *f*(*t*) in the integrals (3.6) and (3.7) leads thus to integrals of the general type
3.10and
3.11Just as with the Legendre integrals (3.5) arising in the context of the BI-F, integrals of the general forms (3.10) and (3.11) can be quickly evaluated in closed form, as shown in appendix B. The algorithm for *B*_{m}(*α*,*β*) is designed so that the value *B*_{m}(*α*,0) coincides with the limit value for *β*↘0. In our context, the algorithm therefore provides the convenience of making equation (2.7) valid when ** x**∈∂

*Ω*.

Substituting equations (3.8) and (3.9) into equation (2.7) and carrying out the changes of variables outlined above gives, for every ** x**-value that we choose on the boundary ∂

*Ω*, a linear relation between all the and entries, of the general form 3.12Here all the 8

*N*entries

*σ*

_{m,ν}and

*τ*

_{m,ν}(

*m*=0,1,…,

*N*−1,

*ν*=1,…,4) are explicitly known numbers (obtained by evaluating the integrals

*A*

_{m}(

*α*,

*β*) and

*B*

_{m}(

*α*,

*β*) described in equations (3.10) and (3.11)). Figure 4 illustrates the ‘structure’ of these 4

*N*relations (3.12) (one relation for each of the 4

*N*sample points along ∂

*Ω*) when they are written out in matrix/vector form. Having arrived at this linear algebra relationship, the numerical solution proceeds as follows: Along each side on which either

*q*

^{(ν)}(

*t*) or is given, use equations (3.8) or (3.9) to compute the known entries in the long column vector shown in figure 4. After these are multiplied with their respective columns in the coefficient matrix and moved to the system's right-hand side, we are left with a square linear system for the unknown Legendre expansion coefficients.

Numerical tests showed that there were no benefits to be gained by using a larger number of sample points than unknown coefficients. They also showed it to be a good strategy to select the sample points ** x** along each of the four sides as the zeros of the degree

*N*Chebyshev polynomial. This clusters nodes towards the ends of each side, but avoids evaluations exactly at the corners.

### (c) Collocation with sets of Laplace solutions formulation

We will not attempt here to present a systematic discussion of how to best use this approach for general Dirichlet-to-Neumann problems. Our present ad hoc implementation of choosing some set of Laplace solutions and then forcing their combination to satisfy all the BCs, requires only a few lines of code and can yield excellent results. A different CLS version (enforcing the BCs weakly via Lagrange multipliers) was introduced in Georgiou *et al.* (1996) and later used in Elliotis *et al.* (2002). For our implementation, in the last two test cases of this paper, the function set we use already satisfies the BCs along two of the four sides, and it only remains to collocate along the remaining two sides. We do this by selecting *N* nodes along each of these two sides (again distributed as the Chebyshev zeros), and then impose that the function set (2.10), truncated to 2*N* terms, obeys the boundary data at these points. The structure of the resulting linear system is seen in figure 5.

## 4. Statements of the three test cases

We will consider three different choices of BCs for the domain in figure 1

### (a) Test case 1 (no singularity)

This test case is chosen so that the provided boundary data will be satisfied by
4.1The above function is the real part of the entire function *e*^{z+1+2i}. We specify *q* along the sides 1, 2, 4 and *q*_{n} along side 3 (all according to equation (4.1)), and the task is to compute the unknown boundary values along each of the four sides, i.e. *q* along side 3 and *q*_{n} along sides 1, 2, 4.

### (b) Test case 2 (one corner singularity)

Following Elliotis *et al.* (2005) and Igarashi & Honma (1996), we specify:
This problem features one singular point located at the 135^{°} corner (*x*,*y*)=(0,1).

### (c) Test case 3 (all corners singular)

This case is taken from Elliotis *et al.* (2002). A reformulation of Elliotis *et al.* (2002, eqn (4)) in terms of *q*(*x*,*y*) along sides 1 and 2 gives the following BCs:
In this case, there are singularities at all four corners. The strongest singularity is located at (0,1), and is locally of the type given by equation (2.8).

## 5. Numerical results for the three test cases

### (a) Test case 1

In terms of the variable *t*, the values for *q*(*t*) and *q*_{n}(*t*) along the four sides follow directly from equation (4.1)
5.1Since we have the exact answer available, it is straightforward to calculate the largest absolute error in any of the four unknown functions. This yields the result displayed in figure 6. Several observations can be made:

— The BI-F method requires the selection of the parameter

*R*—the radius of the circle within which the*k*-values are to be taken. Numerical experimentation suggests that the choice of*R*is not critical and, in the present test, values within the range 8≤*R*≤20 generally worked well. A good general strategy is to choose*R*so as to keep the system condition number low.— For the CLS method, we used as basis functions {1,{Re(

*z*^{k}),Im(*z*^{k}),*k*=1,2,…}}, i.e. {1,*x*,*y*,*x*^{2}−*y*^{2},2*xy*,…}.— The excellent performance of the CLS approach seen in figure 6 is somewhat misleading, as the solution to the present test problem is the real part of an entire function. Typical Laplace solutions will feature singularities in domain corners or, at the very least, in the near vicinity of the domain, ruining spectral convergence of this approach (or convergence altogether), unless care has been taken in finding a specially designed set of expansion functions. The other two methods (BI-F and BI-T) do not assume continuity at corners or any regularity when solutions are continued to outside the domain. This makes the BI approaches applicable to more general situations.

— If the condition number for a linear system approaches 10

^{16}, there is a serious danger that rounding errors have compromised many of the digits (when using standard 16 digit arithmetic precision). Although not an ideal measure of numerical sensitivity (especially when very high), the low condition numbers seen for both of the BI implementations in figure 7 indicate good robustness for these methods against numerical ‘noise’ (due to either machine rounding errors or data inaccuracies).

### (b) Test case 2

The functions *g*_{μ}(*r*,*θ*) defined in equation (2.8) satisfy the BCs on sides 3 and 4. We base our CLS implementation on these functions (*μ*=1,2,3,…), collocating with the BCs along sides 1 and 2. This gives (in less than 10 lines of Matlab code) the coefficient set shown in table 1. These coefficient values are consistent with those presented in Elliotis *et al.* (2005), but are probably more accurate since increasing the number of coefficients leaves the values shown in table 1 invariant to all displayed decimal places (rather than causing higher decimal places to ‘drift’, leaving no reliable significant digits in Elliotis *et al.* (2005) for coefficients past *μ*=20). The coefficients decay at a steady exponential rate that suggests a radius of convergence *R*_{c} of about *R*_{c}=2.4. Thus, the region of convergence extends across the full problem domain (for which no point is located further from the expansion corner than ). Such a large value for *R*_{c} should not be expected for more general BCs, but is a great convenience for this particular test problem, since we therefore obtain a highly accurate reference solution to compare the numerical BI approaches against.

Since for this problem when (recalling (2.8)), we cannot again display the max norm of the errors around the boundary. Instead, we display in figure 8 the values computed by the BI-F method and their errors as functions of *t* along the four sides. In this computation, we have used a very coarse discretization: *N*=5, i.e. only Legendre modes *P*_{0}(*t*) to *P*_{4}(*t*) along each side, and a total of 28 *k*-values (chosen within a radius *R*=15—as noted before not a particularly critical parameter setting). The errors are clearly dominated by the singularity at the 135^{°} corner, corresponding to *t*=1 for *q*^{(3)}(*t*) and to *t*=−1 for . We next supplement the 20 expansion functions (the 5 Legendre polynomials along the 4 sides) with the most singular corner function (along sides 3 and 4 only). For the same parameter settings, we now obtain the result shown in figure 9. The corner error has dropped in size by several orders of magnitude, and it is no longer the dominant error. Also the errors are now smaller around all sides. In this calculation, the coefficient for the *g*_{1}(*r*,*θ*)-term becomes 1.1277—a reasonable approximation to 1.1280. Increasing *N* (above 5) improves this approximation, and reduces the errors further along the sides. The inclusion of one (or more) corner singularity function leads again to integrals that can be evaluated in closed form, as is shown in appendix C.

Since the primary goal of the present study is to verify the viability of the new BI-F formulation, we do not present test results here with the BI-T formulation for Test cases 2 and 3. We only note that the same type of corner correction procedure is again available, although algebraically more complex. The counterparts to the integrals *I*_{1}, *I*_{2}, *J*_{1}, *J*_{2} in appendix C can then be obtained from the closed-form expression for , which combines three _{2}*F*_{1} hypergeometric functions.

### (c) Test case 3

The CLS approach—immediate collocation along sides 1 and 2 using the singularity functions *g*_{μ}(*r*,*θ*) from equation (2.8)—gives for the leading coefficients the values shown in table 2. These values (limited to the number of digits that can be believed to be reliable) are in excellent agreement with the result obtained in Elliotis *et al.* (2002, table III), where they in turn were found to compare favourably against still other calculations in the literature. The coefficients decay to zero much slower than in Test case 2. The present expansion's radius of convergence must satisfy *R*_{c}≤1 because of the singularity at the origin (*x*,*y*)=(0,0). In spite of the fact that the leading coefficients, as seen in table 2, may give an appearance of decaying reasonably fast, the expansion cannot converge over the full domain. We nevertheless take it as the reference solution to compare against. The BI-F method (with *N*=5, *R*=8, no corner correction) produces the results shown in figure 10. The largest discrepancy is clearly at the 135^{°} corner. The wobbly fine structure of the ‘error’ along side 1 corresponds to ‘noise’ in the badly determined reference solution and not to errors in the BI-F solution. We leave it as a future project to pursue corrections at all the four corners and to compare against more accurate reference solutions (which for example could be obtained via a conformal mapping approach).

## 6. Conclusions

BI methods are routinely used in a wide variety of applications, and the literature surrounding them is immense—including journals dedicated to the topic. It was therefore exciting when an altogether different BI formulation was proposed recently, based on the new ‘global equation’ approach arising from analytical function theory. This spectral space BI formulation is algebraically no more complex than earlier BI versions, and it leads to genuinely different numerical implementation opportunities. We have shown here that it readily provides spectral accuracy and that it compares well against an also spectral implementation of the usual BI-T formulation. In particular, the new method gives rise to integrals that are either available in closed form or, when involving boundary data, have integrands that are just as smooth as the data. This is in distinct contrast to standard BI formulation, for which the required integrals can feature logarithmic singularities and sharp gradients.

The scope of this study has been limited to Laplace's equation on a single polygonal domain. Its goal has only been to highlight a novel opportunity that requires future studies before its place in the arsenal of numerical PDE methods can be properly assessed. For a generalization of the present BI-F implementation to the Helmholtz and the modified Helmholtz equations, see Davis (2011).

## Acknowledgements

The work of B.F. was supported by the NSF Grants DMS-0611681 and DMS-0914647. N.F. was supported by the NSF Grants ATM-0602100 and DMS-0934317. National Center for Atmospheric Research is supported by the National Science Foundation. Numerous discussions with Prof. Athanassios S. Fokas and with Christopher-Ian Davis are gratefully acknowledged.

## Appendix A. Selection of *k*-values: Halton nodes

The *k*-values we use should be well separated from each other, so that all the relations that are obtained by substituting them into equations (2.2) and (2.3) will be relatively independent of each other. On the other hand, *k*-values of large magnitude will make the integrals strongly dominated by a very narrow range of *z*-values. Numerical experiments, described in Davis (2011, ch. 5), suggest it to be a good strategy to distribute the *k*-values somewhat randomly within a domain |*k*|≤*R*, where *R* has to be appropriately chosen. Direct use of a random number generator is a bad strategy since, by the nature of randomness, dense clusters will occasionally form. Several options are available for creating scattered numbers that lack any apparent regularities, but which nevertheless will ‘avoid’ each other in a systematic manner. Halton nodes have this feature, and are particularly easy to generate. For example, the routine *haltonset* in Matlab's Statistics Toolbox or *haltonseq* from Matlab Central can be used to generate node sets in a *d*-dimensional unit cube. We use *d*=2, and then ignore those falling outside the inscribed circle to the unit square. Translation and scaling will then provide good node sets satisfying |*k*|≤*R*.

## Appendix B. Closed-form expressions for the integrals (3.10) and (3.11)

The algorithm for evaluating the integrals *A*_{m}(*α*,*β*) and *B*_{m}(*α*,*β*) defined by the integrals (3.10) and (3.11) (for *α*≠±1), starts by defining *z*=*α*+i*β*, , . One calculates
B 1and then discards the imaginary part of the obtained results. Similarly,
B 2and again one discards the imaginary parts. The complex conjugates in the definitions of *L*_{1} and *L*_{2} ensure that, when using the standard convention for the branch cut of the logarithm function, the result for *B*_{m}(*α*,0) will agree with the limit *B*_{m}(*α*,*β*) for *β*↘0.

The proofs for equations (B 1) and (B 2) are similar, so we give it only for the first case. The first three lines can be verified directly. Since
the task is to show that the recursion holds for , *m*=3,4,… The standard three term recursion for Legendre polynomials can be written
Subtracting (2*m*−1)*zP*_{m−1}(*t*) from both sides, multiplying by and integrating in *t* from −1 to 1 gives
We next replace (2*m*−1)*P*_{m−1}(*t*) according to the identity and integrate by parts. The integral above then becomes
and the recursion follows.

For given smooth functions *f*(*t*), integrals of the types (3.6) and (3.7) can be quickly and accurately evaluated by (i) evaluating *f*(*t*) at *Chebyshev* node locations, in order to obtain its Chebyshev expansion coefficients by means of a fast Fourier transform (FFT), (ii) converting these to Legendre expansion coefficients (for example, by the *cheb2leg* routine from Matlab Central), and (iii) applying equations (B 1) and (B 2) to the resulting integrals.

## Appendix C. Integrals arising when using corner singularity functions with the BI-F approach

To include the corner basis function in the representations of the boundary functions along sides 3 and 4, the key relations that need to be considered are:

side 3:

side 4:
Here, . Although the integral itself is well known as the incomplete gamma function, it is better to work with *H*(*a*,*z*) because this is an entire function (and not a multivalued function with a singularity at the origin). From a numerical perspective, the continued fraction expansion for *H*(*a*,*z*) is particularly convenient to use (cf. Olver *et al.* 2010, eqn (8.9.1))
For effective numerical evaluation of such expansions (e.g. the algorithms by Steed & Lentz; see for example, Press *et al.* (1992)). This same approach—reformulating boundary integrals to *H*-functions—will apply to every one of the corner singularity functions given by equation (2.8).

- Received January 12, 2011.
- Accepted April 27, 2011.

- This journal is © 2011 The Royal Society