## Abstract

Consider the boundary-value problem for the field *ψ*(* x*) which satisfies the linear partial differential equation in an arbitrary domain with data given on the boundary . It is generally believed that, unless is the union of constant coordinate lines in a separable coordinate system for the operator , the problem cannot be solved by the classical method of eigenfunction expansions. We show how this apparent limitation can be overcome. The key idea is to embed in a larger embedding domain , which is endowed with a complete set of eigenfunctions of the operator , where the

*λ*

_{n}are the eigenvalues. We can now expand

*ψ*(

*) in terms of this set, i.e. . Although the unknown scalars {*

**x***a*

_{n}} can no longer be determined by the use of an inner product, a least-squares procedure which minimizes the error in the boundary data yields the scalars to as high a precision, in principle, as needed. Examples are given of steady heat conduction in two and three dimensions, governed by Laplace's equation, and of Stokes flow in a container, governed by the biharmonic equation, all in non-simple domains. The scope of a powerful classical method has, by this extension, been enlarged very considerably. It is believed that it will be of great use in solving practical, linear boundary-value problems, which until now had to be solved by brute force numerical methods.

## 1. Introduction

It is generally perceived that the method of eigenfunction expansions (MEE), while powerful for simple geometries admitting complete sets of eigenfunctions, is useless for complex geometries. It is also believed that, even when suitable eigenfunctions exist, they are not of use in solving linear boundary-value problems, except in those special cases where the boundary data are prescribed along the coordinate lines. In all other cases, it is believed that there is no recourse except to numerical methods and perhaps to some special methods such as the boundary integral technique. In the classical treatises, such as Morse & Feshbach (1953), Titchmarsh (1958) and Carslaw & Jaeger (1959), a very large number of boundary-value problems are solved by the MEE, but all of them involve only simple geometries in the coordinate systems such as the Cartesian, polar, spherical polar, bipolar, etc. which lead to separable solutions. All the authors cited sought solutions that were exact, in the sense that the expansion coefficients could be determined exactly in terms of the parameters of the problem. There exist, of course, methods which are not purely computational to tackle problems in complex geometries. One can refer to Finlayson (1972) and the references cited therein for the literature on the weighted residual methods, such as collocation, Galerkin and least-squares, which have been used to obtain approximate solutions to such problems.

Almost 50 years on, the situation does not seem to have changed. A glance through some contemporary texts in mathematical physics and engineering show that they seem to hold the same perception—that the MEE is only for simple special geometries. After discussing the method and giving examples, Greenberg (1998) concludes on p. 972 that it is essential that the boundary and initial conditions be given on ‘constant coordinate curves’. In a similar manner, Hassani (1999) concludes (p. 526) ‘Problems most suitable for Cartesian coordinates have boundaries with rectangular symmetry such as boxes or planes’. In dealing with two-dimensional, steady-state diffusion in plates, Gustafson & Wilcox (1999) state ‘For most plate shapes—the problem can only be solved by numerical methods’ (p. 598). The whole treatise of Lurie & Vasiliev (1995) is an application of biharmonic eigenfunctions to the two-dimensional problems of elasticity. Once again, all the problems are restricted to those of simple geometry; in fact, here the authors further restrict themselves to those problems where biorthogonality leads to the expansion coefficients. Thus, it is still perceived that geometry is a serious limitation for the MEE.

The purpose of this communication is to suggest a very general method, the method of embedding, to overcome this limitation and to show by examples that it works for a variety of problems, thus enormously extending the scope of the MEE. This method should not be confused with the ‘fictitious domain’ method used primarily in finite-element numerical solutions of complicated flow or flow–structure interaction problems. It is, of course, possible that there have been sporadic, individual attempts, maybe even successful ones, to overcome the geometrical limitations of the MEE. But if there have been such attempts, they do not seem to have been sufficiently general or clear enough to have made an impression on the scientific and academic community. We have not so far been able to find, nor have we heard of, the specific use of such a method of extending the scope of the MEE, either in the Russian or any other literature.

## 2. The extension to arbitrary domains

Let us suppose that in a bounded domain in *R*^{2} or *R*^{3}, with boundary , we wish to determine a field *ψ*(* x*) satisfying the linear equation(2.1)with suitable conditions prescribed on the boundary . Here, is a differential operator. Although we will treat

*ψ*(

*) as a scalar field, all the arguments go through even if it is a vector field, in which case would be a matrix or vector differential operator. Our aim is to solve for the field using an eigenfunction expansion procedure. Let us suppose that there is no standard coordinate system in which is the union of constant coordinate lines or that no suitable set of eigenfunctions of , complete on ,1 can be found.*

**x**We now make the following assumptions: (i) that the problem in is well posed and has a unique solution, (ii) that is convex and simply connected, (iii) the boundary data are continuous on and (iv) there exists a domain with boundary , in which can be embedded, which is endowed with a complete set of eigenfunctions of . It is possible that some of these assumptions are unnecessarily restrictive, but at this stage, where a new idea is to be tested, caution may be in order. In fact, the restriction to convex domains may be unnecessary, and we will indeed consider an example where the method works even though is not convex.

The basic idea is illustrated in figure 1. Embed in as shown. Now it immediately becomes clear that we can attempt to use the eigenfunctions in to solve the problem (2.1) in . Let be the set of eigenfunctions of , which is complete on ; here, we deliberately show the explicit dependence on the eigenvalues *λ*_{n}. Initially, we assume that the eigenfunctions are real; this restriction will be lifted in the following. The ordering of the eigenvalues is in increasing values of their magnitude. It is often the case that {*ϕ*_{n}} is itself the union of a finite number of infinite subsets, each of which is a set of eigenfunctions, i.e.(2.2)say, where {*ϕ*_{n}} is made up of three independent subsets , and . This structure arises because the different sets either satisfy different symmetry conditions or homogeneous conditions on different subsets of . This will become very clear from the examples in §3.

To keep the notation simple, we will deal with the set of eigenfunctions {*ϕ*_{n}} alone, with the clear understanding that it may be the union of a number of subsets. We now write down a formal eigenfunction expansion in for the field *ψ*(* x*),(2.3)where the unknown scalars {

*a*

_{n},

*n*=1, 2, …} have to be determined from the boundary data on . We note at this stage the considerable progress that has been made. The formal expansion (2.3) satisfies the field equation (2.1) exactly. Moreover, since the {

*ϕ*

_{n}} are complete on the larger, embedding domain , we expect the set to be complete on also.2 The only remaining task is to determine the scalars

*a*

_{n}. Any clinging to the notions of orthogonality or biorthogonality at this stage will lead to unnecessary complications, if not failure. Here we determine the coefficients by the least-squares procedure, which has already proven to be reliable and robust (Shankar 1993, 1997, for example).

First truncate the series in (2.3) to the first *N* terms and choose *M* (*M*>*N*) points on the boundary , not even necessarily evenly spaced. Let *ψ*_{0}(* x*) be the given boundary data on . Then the square of the error (the difference between the value given by (2.3) and

*ψ*

_{0}(

*)) at a given point*

**x**

**x**_{k}on the boundary is just(2.4)while the total error squared over all the

*M*points on the boundary is just(2.5)Now the scalars

*a*

_{n}have to be chosen so that

*E*

^{2}is a minimum, i.e. so that ∂

*E*

^{2}/∂

*a*

_{i}=0,

*i*=1, 2,…,

*N*. This leads to a linear system of

*N*equations for these unknowns,(2.6)Naturally one expects that as

*N*,

*M*→∞, the

*a*

_{n}will converge to their asymptotic values and calculations fully justify this expectation.

It was assumed above that there was a single condition, i.e. that *ψ*(* x*)=

*ψ*

_{0}(

*) on ; in general, depending on the nature of*

**x***ψ*and there could be more. If there are

*I*such conditions, we define

*I*pointwise error squared functions on such as equation (2.4) and write for the total error squared . The rest of the procedure is then as before.

In general, the eigenfunctions {*ϕ*_{n}} of in will be complex. However, the procedure is exactly as above, except that now the scalars {*a*_{n}} will also have to be complex. Instead of equation (2.3) we now write(2.7)Truncation to *N* terms followed by the least-squares procedure leads to 2*N* equations for the 2*N* real unknowns .

Thus, whether the eigenfunctions are real or complex, whether the domain is two- or three-dimensional, or whether the field is a scalar or a vector, the same uniform procedure leads to a solution of the boundary-value problem for equation (2.1) in the arbitrary domain . In §3, this procedure will be illustrated by a number of examples which will show how well it works.

## 3. Illustrative examples of eigenfunction expansions in arbitrary domains

### (a) Laplace's equation in the plane: planar, steady heat conduction

Consider steady, planar heat conduction in the solid shown in figure 2; the domain here is the union of a semicircular disc of radius *r*=1 and a rectangle of dimensions *h*×2. The temperature *ψ*(* x*) satisfies Laplace's equation, i.e. ∇

^{2}

*ψ*=0, and we will here assume Dirichlet data,

*ψ*(

*)=*

**x***ψ*

_{0}(

*), on the boundary. Suitable eigenfunctions do not seem to exist for this geometry, where the boundaries are not all coordinate lines in some suitable system. Using the methods of §2, we will solve this problem in three different ways.*

**x**First, choose as a disc of sufficiently large radius centred at the origin O (see figure 2*a*), such that the boundary encloses . This domain has the complete set of real eigenfunctions . Thus, on we have the expansion(3.1)where the real coefficients *a*_{n} and *b*_{n} can be determined by truncation to *N*_{1} terms and using a least-squares procedure with *M* points on the boundary. Figure 2*b* shows the field for a particular choice, , of the boundary data. Note that there is no symmetry in the field as the boundary data are asymmetrical. With *M*=3*N*_{1}, we find that the maximum error on the boundary is less than 0.7×10^{−2} for *N*_{1}=10, less than 0.29×10^{−3} for *N*_{1}=25 and less than 0.88×10^{−4} with *N*_{1}=50.

A natural and obvious question is: if other embedding domains exist, will they all give the same answer? Of course, if the theory presented here is correct they, indeed, must. We illustrate this by considering the rectangular embedding domain with boundary , part of which overlaps and choosing as origin O′, the lower, left corner of the rectangle. Now the eigenfunctions are of the form ; and ; note that it is *essential* that both sets be used, since it is the union of these which is complete on . Let 3 and ; the field can then be written as,(3.2)where the *a*_{1n}, *a*_{2n}, *b*_{1n} and *b*_{2n} are the four sets of real scalars to be determined. We once again emphasize the importance of taking eigenfunctions from both boundaries in order to get a complete set. The results for the same boundary-value problem as above, with truncation to *N*_{2} terms, and with *N*_{2}=25 and *M*=150 (i.e with the same number of unknown coefficients and boundary points as in the first case) are shown in figure 2*b* as dashed lines. The two results are indistinguishable. The maximum error on the boundary is now less than 0.16×10^{−3}, i.e. a little greater than in the earlier case.

Similarly, one may ask: what if one uses embedding eigenfunctions which satisfy Neumann boundary conditions on , i.e. if ∂*ϕ*/∂*n*=0 on ; will they work? Again, the theory of §2 would imply that they would. In figure 2*a*, is the larger rectangle which embeds ; it is of size (*h*+1+2*b*)×(2+2*b*). The eigenvalue sets are just as before, but the sines in (3.2) get replaced by cosines. We once again have a representation similar to (3.2) and the procedure goes through. The results of this calculation are shown as dot-dashes and these too are indistinguishable from the earlier two; however, the maximum error is now less than 0.17×10^{−2} for *N*_{2}=25 and *M*=150. In all three cases, *the error over most of the boundary is an order of magnitude smaller than the maximum error*. It is impressive and reassuring that three different embedding procedures agree so well on a problem that could not have been solved by the classical method.

Although in the problem considered above the boundary condition was of the Dirichlet type, the same method will, of course, work with the general homogeneous boundary condition of the form on . In fact, it appears that as long as the field satisfies a linear equation, even nonlinear boundary conditions can be handled. This fact has implications for the study of capillary–gravity waves on the free surfaces of liquids in containers.

### (b) Laplace's equation in R^{3}: three-dimensional, steady heat conduction

Figure 3*a* shows an axisymmetric body of complex shape subjected to a non-axisymmetric boundary temperature distribution; thus, the temperature distribution in the solid is fully three-dimensional. The body section is made up of two circular arcs of different radii, 1 and *ρ*_{0}, with the other two segments AB and CD being straight lines; AE is the axis of symmetry. Once again ∇^{2}*ψ*=0 in the body. Traditional eigenfunction methods fail in this complex geometry. We choose the cylinder of half section AFGE as the embedding domain and use cylindrical polar coordinates with the origin at E. The temperature can be expanded as follows:(3.3)where the radial eigenvalues *λ*_{mn} satisfy *J*_{m}(*λ*_{mn})=0, while the axial ones are simply *nπ*/*h*. Note that the terms in the first curly braces constitute the radial eigenfunctions satisfying homogeneous Dirichlet conditions on while the ones in the second constitute the axial ones satisfying those conditions on . Together, they form a complete set on . The three sets of scalars, {*a*_{mn}}, {*b*_{mn}} and {*c*_{mn}}, have to be determined by the usual procedure of trunction to *N*_{1} terms and using least-squares. Figure 3*b* shows the temperature distribution for a particular surface distribution having an azimuthal dependence corresponding to *m*=2. The form of the distribution is , where on BC and on DE when *θ*_{D}<*θ*<*θ*_{D}+2*δ*; *f*(*θ*)=0 elsewhere. Here *a*, *b*, *θ*_{1}, *θ*_{D} and *θ*_{3} are scalars. Computations show that with *N*_{1} equal to 10, 20 and 30 the maximum error on the boundary was less than 5×10^{−3}, 8.1×10^{−4} and 4.1×10^{−4}, respectively; over most of the boundary the error was an order of magnitude less. Note that, whereas the problem in the previous subsection could have been solved by (numerical) conformal mapping, this one cannot be so handled. Here, one would normally require a finite difference or finite element, fully numerical treatment.

### (c) The biharmonic equation: two-dimensional Stokes flow in a container

As a more complicated example involving complex eigenvalues, we now consider slow viscous flow in a container full of liquid, the steady fluid motion being generated by the motion of the lid. In figure 4*a*, the width is 1, the straight sidewalls are of height *h*, while the bottom is trench-like, with the bottom walls at an angle *θ*_{0} to the horizontal. Since the geometry is symmetric about AD and the lid speed *v*_{x}(*x*) will be assumed to be symmetric in *x*, the velocity field will be antisymmetric about AD. The origin is at the mid-point of AD and (*x*, *z*) are the Cartesian coordinates used. The field is represented by the stream function *ψ*(*x*, *z*) which satisfies the biharmonic equation, ∇^{4}*ψ*=0, in the domain. Standard eigenfunctions do not exist for this geometry. We now choose the rectangular domain FBGH shown in figure 4*a*, for which eigenfunctions do exist, as .

Recall the antisymmetric and symmetric (Papkovich–Fadle) reduced eigenfunctions on (Joseph & Sturges 1978; Shankar 1993):(3.4a)(3.4b)(3.5a)(3.5b)All the eigenvalues, the roots of (3.4*b*) and (3.5*b*), are in this case complex. The dominant eigenvalues are *λ*_{1}≈4.2124+2.507i and *μ*_{1}≈7.4977+2.7687i. Let *h*_{1}=0.5 tan *θ*_{0}+*h* be the maximum depth of the container and let *α*=*h*_{1}/2 and *β*=0.5/*h*_{1}. One can now write down an eigenfunction expansion for the stream function,(3.6)where the *a*_{n}, *b*_{n}, *c*_{n} and *d*_{n} are all complex scalars that are to be determined by the boundary data. It should be noted that the form (3.6) takes into account the antisymmetric nature of the field about AD. Also that the complete set on involves both sets of reduced eigenfunctions and, in general, acting on both sets of boundaries. If *θ*_{0} were 0, only the *a*_{n} and *b*_{n} would be needed. In all the calculations, the lid speed was taken to be for (1/2−*δ*)<*x*<1/2 and 1 for all other *x*≥0; i.e. uniform except for a smooth drop to the sidewall.

Figure 4*b* shows the results of a calculation for *h*=1 and *θ*_{0}=60°. The streamline pattern shows a single primary eddy, a large corner eddy at the bottom and a second, smaller, corner eddy which is but poorly resolved. The maximum error on the boundary with *N*_{1}=50 and *M*=12*N* was 4.4×10^{−4}. The calculation, based on a simple, unoptimized, almost profligate code, takes 12 s on a Pentium 4 PC running at 2 GHz including the calculation of the field at 100×100 points; the calculation of the coefficients takes only 3 s. This calculation was checked with *N*_{1}=100 and the field was indistinguishable for the most part from that in figure 4*b*; however, the second corner eddy was indeed resolved much better. Once again we point out that for such geometries it was until now believed that the eigenfunction method could not be applied. Interestingly, Ratkowsky & Rotem (1968) used expansions based on the real separable solutions to the biharmonic equation in polar coordinates with least-squares on the boundary to solve Stokes flow in the rectangular container. Although this is a very inefficient choice, especially for large aspect ratios, in the context of the embedding method it is one where the embedding domain can be considered to be a sufficiently large circle embedding the rectangle with the origin in the fluid. Although we have not tried this procedure on the geometry considered here, we believe it should work, though inefficiently.

### (d) The field in a domain which is not convex

It may have been noted that although the convexity of had been assumed, it had not actually been used in the analysis. Following the useful suggestion of a referee, we now consider as a final example a field in the non-convex domain shown in figure 5.

We wish to determine the field *ψ*(*r*, *θ*) satisfying ∇^{2}*ψ*=0 in the half-annular ring <*r*<1, 0<*θ*<*π* with *ψ*(*r*, *θ*)=sin *θ* on the boundary. This is a case where the exact solution to the problem,(3.7)can be written down by inspection. If the embedding method were to be applied to this or a similar domain, the natural embedding domain would be an annulus. However, following the suggestion of the referee, we consider using the *full disc r*<1, even though we know it will be unnecessarily larger than the given domain. This is interesting, and possibly disturbing, for the following reason. For this embedding domain, the natural eigenfunctions are of the type *r*^{n}(cos *nθ*, sin *nθ*) alone; and it may be feared that these will be unable to represent equation (3.7), since the *r*^{−n}(cos *nθ*, sin *nθ*) terms are now missing. It may even appear that this may be a case, in view of the term involving 1/*r* in equation (3.7), where the smooth solution in becomes singular in some sense when extended to . In spite of these apprehensions, the calculations were done using the disc eigenfunctions alone and the results are shown in figure 5. The solid lines represent the field generated by the exact solution (3.7), while the solution by the embedding method with *N*=50, *M*=3*N* is shown by the dotted lines. There is good agreement everywhere except near *θ*=*π*/2, where it is less good as the errors increase in this neighbourhood; these errors could, of course, be reduced with more terms in the series. This example also shows, as might be expected, that the convergence will generally be slower the more unnecessarily large the embedding domain is. We should also caution, in spite of the success in this particular case, that much work remains to be done to establish more definitively when the method will always work and when it may fail.

## 4. Discussion and conclusion

It has been shown here that the eigenfunction method can be applied, contrary to conventional belief, to many bounded domains of complex shape which in themselves are not endowed with complete sets of eigenfunctions. Although the method has been shown to work in a variety of cases, there are many issues that will have to be investigated in the future. It was assumed here, to avoid possible difficulties, that the domains were convex, but nowhere was the assumption explicitly used. In fact, in the example of §3*d*, a non-convex domain was considered and the method appeared to work, even if somewhat inefficiently because of the large size of ; convexity does not seem necessary. Another practically important issue is that of the connectedness of the domain. Here, we had assumed the domains to be simply connected. A simple example, heat conduction in a plate with a central hole, seems to suggest that at least for some multiply connected domains and operators, suitable embedding domains exist; but the method needs to be extended. The general issue is open. Given a number of embedding domains, which one is likely to be ‘optimal’ in some sense? Intuitively, one would guess that optimum would be ones which are as small as possible. However, it is not clear whether the coordinate systems used will play a role in determining the optimum. The field equation (2.1) was homogeneous and one can enquire about the inhomogeneous case. There seems to be no problem here: we only need to find any inhomogeneous solution *ψ*_{inh} and write *ψ*=*ψ*_{inh}+*ψ*_{h}; then *ψ*_{h} can be solved for by the method suggested here.

Similarly, there are issues regarding the implementation of the least-squares procedure on the boundary which need further investigation. One can, for example, ask whether the *M* points on the boundary, for prescribed *N* and *M*, can be distributed ‘optimally’. Would clustering help near corners? For all the problems considered here, the points were distributed equidistant or almost equidistant along the boundary. Our position is that as *N*, *M*→∞ the computed field must converge to the unique solution. It is hard to say more without detailed investigation. The situation is similar to numerical integration using, for example, the trapezoidal rule. As long as the number of intervals is large, the exact nature of the partitioning is irrelevant; however, for particular integrands there may well be ‘optimal partitions’. There are many practically important questions that need to be investigated. In addition, there are, no doubt, important mathematical issues underlying the method that need to be addressed, but these are beyond the scope of this paper and the competence of the author.

Finally, it must be admitted that no attempt has been made in this brief paper to compare the eigenfunction method with other possible methods, since the only purpose here was to extend its scope considerably. However, it seems that where applicable, the eigenfunction method would be more convenient, faster and more accurate than any numerical method. The only methods that we can possibly conceive of as being comparable are the boundary integral or surface-singularity methods (Higdon 1985; Pozrikidis 1992), which also yield fields that satisfy the field equations exactly. Each method has its advantages and its disadvantages. The principal advantages of the boundary integral methods are their generality and their capacity to handle unbounded domains. However, these methods use singular solutions, which need desingularization techniques; further, the boundary needs to be panelled, which is a nuisance. The eigenfunction method appears to be conceptually much simpler and easier to handle. Its main strengths are: (i) its conceptual simplicity, (ii) its high accuracy, (iii) the exact satisfaction of the field equations, (iv) the potential to obtain the asymptotic nature of the field, (v) the absence of the need for any gridding or paneling, (vi) the ability to store the details of the field in a relatively small number of coefficients and (vii) the potential to examine any part of the field in great detail.

The main weakness of the eigenfunction expansion method was its restriction to the simplest of geometries. We believe that this weakness has been eliminated to a considerable extent in this paper.

## Acknowledgments

I would like to thank the referees for their useful comments and Mrs D. Shobha for help in preparing the manuscript.

## Footnotes

↵

^{†}Present address: 33/1 Kasturba Road Cross, Bangalore 560001, India.↵By a complete set of functions on , we mean a set of functions that is rich enough to represent, modulo certain smoothness restrictions, essentially arbitrary data on . The example in §3 will make this notion clearer.

↵Equivalently, we assume here that the solution in has a nice extension to . A referee has pointed out that this may not always be so, in which case the method would fail. While this is true, so far every case that has been tried has been successful.

↵Since the data may not be homogenous at the corners of the rectangle, one needs to add the terms

*c*_{1}*x*+*c*_{2}*y*, for example, with appropriate values for*c*_{1}and*c*_{2}; see Shankar (2004).- Received July 22, 2004.
- Accepted April 5, 2005.

- © 2005 The Royal Society