## Abstract

Motivated by problems involving diffusion through small gaps, we revisit two-dimensional eigenvalue problems with localized perturbations to Neumann boundary conditions. We recover the known result that the gravest eigenvalue is *O*(|ln *ϵ*|^{−1}), where *ϵ* is the ratio of the size of the hole to the length-scale of the domain, and provide a simple and constructive approach for summing the inverse logarithm terms and obtaining further corrections. Comparisons with numerical solutions obtained for special geometries, both for the Dirichlet ‘patch problem’ where the perturbation to the boundary consists of a different boundary condition and for the gap problem, confirm that this approach is a simple way of obtaining an accurate value for the gravest eigenvalue and hence the long-term outcome of the underlying diffusion problem.

## 1. Introduction

Consider Brownian motion in a two-dimensional domain separated into two halves by a barrier. The domain and the barrier are impermeable, but the barrier contains a small hole. (Extensions to multiple barriers and multiple holes are straightforward but lengthy.) Then, the density inside the domain is obtained by solving the diffusion equation(1.1)with the Neumann condition *u*_{n}=0 on the boundary. Separation of variables shows that for a bounded domain the solution can be written as(1.2)where *λ*_{j} are positive decay rates obtained by solving the eigenproblem(1.3)again with Neumann boundary conditions. Note that the lowest eigenvalue is always *λ*_{0}=0, corresponding to steady state. The non-trivial large-time behaviour of the solution is determined by the eigenvalue *λ*_{1}.

The value of the flux, *u*_{n}, in the gap is unknown and is obtained as part of the solution. Physically one can see that if the hole is small enough, the density will become homogenized inside each subdomain much faster than it diffuses through the gap. Hence, if one starts with the density concentrated in one subdomain and the gap is small, the slow diffusion process through the gap will take place with nearly spatially uniform but time-dependent values of *u* in the gap. This motivates the Dirichlet ‘patch problem’, in which one considers subdomains in isolation and replaces the gap by a patch of boundary on which the density satisfies *u*=*C*, where *C* is a constant. Without loss of generality, we can take *C*=0 (for multiple gaps, things are not so straightforward); this is exact for symmetrical geometries (see §5).

Ward & Keller (1993; hereafter WK) obtained, by matched asymptotic expansions, the leading-order behaviour of perturbed eigenvalue problems when the boundary condition is changed over a small region of the boundary of size, *ϵ*. For two-dimensional cases, the change in the eigenvalue can be *O*(|ln *ϵ*|^{−1}) and the leading-order solution can be a poor approximation to the actual eigenvalue. Ward *et al*. (1993; hereafter WHK) presented a hybrid technique that sums all the terms in the expansion in (ln *ϵ*)^{−1} and hence gives much more accurate results.

Similar problems have been considered by a variety of authors. Extensions to the problem of low Reynolds number flow past a cylinder are given in Kropinski *et al*. (1995) and Keller & Ward (1996), to an oxygen transport problem in Titcombe & Ward (2000), to optimization of eigenvalues in Kolokolnikov *et al*. (2005). These works all draw from the original WK work. There has been a considerable further amount of more abstract work, for example, the recent results in Denzler (1999). The wide-ranging book of Maz'ya *et al*. (2000) presents a great many results on this and related problems, as well as further references. 9.1.3 and 9.1.5 of this book touch on the current problem, but the development is essentially formal and no constructive results are given.

For gaps in two dimensions, we give a new constructive method that is based on known solutions to fluid flow and wave scattering problems involving bodies with edges. Several explicit calculations are given and the general procedure deduced. The directness of the method is evident and we recover some results of WHK, as we must, without requiring any asymptotic matching. In addition, we show how to develop a small-*k* approximation to the Green's function, for any domain, that furnishes a more straightforward means of estimating the *O*(1) terms. We describe the approach by applying it to specific geometries (half disc and rectangle for the patch problem in §2, annulus disc, disc with barrier and multiple rectangles for the gap problem in §3). Then in §4, we look at a time-dependent problem involving diffusion through multiple connected rectangles. In §5, we apply the direct matched expansions technique to the half-disc problem to derive perturbations in the higher Fourier modes and show the usefulness of an undetermined scale factor in the WK method. We conclude in §6.

## 2. The patch problem

This is(2.1)with Neumann boundary condition *ϕ*_{n}=0 over most of the boundary, except for a small Dirichlet patch with dimension *O*(*ϵ*) where *ϕ*=0. The unperturbed problem always has the eigenvalue 0, corresponding to the steady-state solution to the original diffusion equation.

The leading-order value of the lowest eigenvalue obtained using the WK method is the same for all patch problems: it depends only on the value of the unperturbed eigenfunction which is essentially the inverse area for the gravest mode, and on the solution of an inner problem which is the same for all geometries. The result is(2.2)where *A* is the area of the domain. This estimate is handicapped by the very slow decay of the inverse logarithm. Often, the determination of further terms soon becomes unmanageable, putting an efficient estimate out of reach. Previous experience leads us to expect an estimate of the form *π*[*A* ln(2/*ϵ*)+*c*]^{−1}, where *c* is a constant that the first correction term in equation (2.2) suffices to determine. We present a method that directly gives an estimate for *c* that effectively eliminates the need to determine an asymptotic expansion involving inverse powers of ln(2/*ϵ*).

Our approach is to write the problem as an integral equation and expand it close to the patch. Then, we can extract the dependence of the eigenvalue on the patch location and width by keeping more and more terms in this expansion. The leading-order solution requires the solution of an associated Poisson problem, as does WHK.

### (a) The half disc

Suppose that *ϕ*(*r*, *θ*) satisfies(2.3)subject to(2.4)and(2.5)(We take *λ*=*k*^{2} here for convenience.) The Dirichlet patch corresponds to the region |*x*−*x*_{0}|<*ϵ* of the *x*-axis, on which *ϕ*=0 (the symmetry allows us to take *x*_{0}>0 in the calculations of this section and, as expected, the results depend only on ). Then, we may write(2.6)(2.7)Here, the Green's function *G*(*x*, *y*, *x*′, *k*^{2}) satisfying(2.8)and(2.9)is given by(2.10)where *ϵ*_{0}=1 and *ϵ*_{n}=2(*n*≥1) is the Neumann symbol. (The first term in equation (2.10) is the Green's function of the Helmholtz operator in an unbounded domain and Graf's addition theorem enables the Neumann boundary condition to be satisfied by adding separated solutions of the Helmholtz equation.) The appearance of the delta function ensures that the unit inward flux is confined to the point (*x*′, 0). It is assumed that , so that the denominators in *G* do not vanish.

The behaviour of *G* along the patch is then given by(2.11)(2.12)where we must require that *k* not be too large, to be precise *kϵ*=*o*(1). We have defined a new function *G*^{*} according to(2.13) and denote derivatives of *G*^{*}(*x*_{1}, *x*_{2}, *k*^{2}) with respect to *x*_{1} and *x*_{2}, respectively, and, like *G*^{*}, are evaluated here and below at (*x*_{0}, *x*_{0}). We have, on inserting the small-argument form of *Y*_{0} into equation (2.11),(2.14)

(2.15)

The new function *G*^{*} (a function of *x*_{0} only), also used by WHK, is the key to the approach and may be obtained numerically from the solution to the following problem, which is just the previous Helmholtz equation with the logarithmic term removed. Write(2.16)Then,(2.17)with appropriate Neumann conditions on the boundary. We obtain(2.18)

Another formulation is to write(2.19)This leads to(2.20)along with appropriate Neumann conditions on the boundary. Then,(2.21)

Returning to the main problem, the square-root edge singularity is anticipated by setting(2.22)where *T*_{n} is a Chebyshev polynomial of the first kind, chosen because the square-root factor is the weight function in the orthogonality property of {*T*_{n}:*n*≥0} (Davis & Scharstein 1993). Substitution in equation (2.6), with (*x*, *y*)=(*x*_{0}+*ϵs*, 0), i.e. on the patch, then gives(2.23)We note an additional advantage of the Chebyshev polynomials, namely(2.24)and substitute into equation (2.12) to obtain(2.25)The orthogonality properties of Chebyshev polynomials yield the equation(2.26)The left-hand side of equation (2.26) corresponds to the zero Dirichlet boundary condition in the patch. More complicated Dirichlet boundary conditions could be considered by expanding *ϕ*(*x*_{0}+*ϵs*, 0) as a Chebyshev series on the patch. The truncation procedure below would then have non-zero left-hand sides in the equations.

The *T*_{0}(*s*) and *T*_{1}(*s*) components of equation (2.26) give the coupled equations(2.27)(2.28)which together yield *a*_{1}/*a*_{0}=*O*(*ϵ*) and lead to the eigenvalue equation(2.29)Keeping two terms in the Chebyshev expansion suffices to obtain equation (2.29). Clearly this process can be continued to as many terms as needed in the series equation (2.22). However, more terms will be needed in the non-singular parts of the Green's function (2.11), both from expanding the *Y*_{0} function and from the series.

Equation (2.29) is valid for all eigenvalues *k* that are not too large. Neglecting the order term, we can solve for *k* numerically using standard root-finding procedures. Alternatively, we can obtain leading-order predictions by expanding for small *k*,(2.30)and for small (*k*−*j*_{1m}) (*m*>0),(2.31)These estimates apply to a perturbation of the azimuthal mode *j*_{0}(*kr*), which has eigenvalue *k*=*j*_{1m}, where *j*_{1m} is the *m*th root of . Thus, and *J*_{0}(*j*_{1m})*Y*_{1}(*j*_{1m})=2/(*πj*_{1m}).

To summarize, and in order of decreasing accuracy, we now have three predictions for the lowest eigenvalue: equation (2.29), which requires numerical root-finding, the small-*k* and small-(*k*−*j*_{1m}) approximations (2.30) and (2.31), and the basic WK result(2.32)

To gain an idea of the accuracy of the different results that have been presented, figure 1 shows the eigenvalue *λ* as a function of *ϵ* for the cases *x*_{0}=0 and 0.3. The curves show values obtained using the formulae above. The initial guess for the root-finding procedure was the WK prediction (2.2). Full numerical solutions to equation (2.3) with appropriate boundary conditions were used for comparison.

The full solutions were first computed using the Matlab finite-element toolbox. Four levels of grid refinement were used. Numerical tests showed that the eigenvalues had converged to two significant figures in the range of *ϵ* down to 10^{−4}. The finite-element toolbox used cannot deal with barriers of zero width. (This may not be a limitation for physical problems but is awkward when comparing with exact solutions.) A spectral approach was also used, following closely Weideman & Reddy (2000; see also Trefethen 2000), with 50 Chebyshev functions in *r* and 50 Fourier modes in *θ*. The spectral approach gave plausible results for *x*_{0}=0, but is so poor for *x*_{0}=0.3 that the results are not shown. The failure of the spectral method is to be expected, since the collocation points on the boundary grotesquely under-resolve the Dirichlet patch. For example, when *x*_{0}=0.3 the Dirichlet patch corresponds to a single collocation point on the boundary for *ϵ*=10^{−2}. The step-like structure of the circles is simply a result of the number of collocation points in the Dirichlet patch being quantized and staying fixed over a range of *ϵ*. This is not to say that an efficient spectral method cannot be found. Rather we see that for an equivalent amount of programming work (in fact quite a bit less), the asymptotic approach is a fast and accurate option. The numerical solution to the Poisson problem (2.17) both using finite elements and spectral methods was easy, fast and accurate.

Representative timings are given in table 1: the main conclusion to be drawn is that the time in the asymptotic method is independent of *ϵ*. The FE method gets more expensive as more and more grid refinement is required (the refinement is done automatically by Matlab).

We have not compared the lowest asymptotic eigenvalues for the half disc to exact results for *ϵ* down to 10^{−8}, but they are clearly getting more accurate as *ϵ* decreases. We also note that similar asymptotic expansions in cases where an exact solution exists were shown by WK to be very accurate.

### (b) The general case

The methods of the previous subsection are evidently independent, except for the construction of Green's function, of the half-disc geometry. Suppose that *ϕ*(*x*, *y*) satisfies ∇^{2}*ϕ*+*k*^{2}*ϕ*=0 in a domain *D* subject to the homogeneous Neumann condition ∂*ϕ*/∂*n*=0 on the piecewise smooth boundary *C* except on a small arc *C*_{p} where *ϕ*=0. The latter is assumed to be not close to a corner, which case needs separate discussion, as illustrated by the example in §3*b*. Dimensionless parameterization of the small arc is achieved using a length-scale *L* associated with *D* to denote distance along the arc by *L*(*s*_{0}+*ϵs*);|*s*|<1. Using the customary outward normal direction, we generalize equation (2.6) as(2.33)in which the source point that generates *G* is on *C*_{p} (in an abuse of notation we replace the second vector argument of *G* by the corresponding value of arc length). We observe that, if **r**,

**′ are two points separated by a distance**

*r**ϵL*(

*s*−

*w*) along an arc of curvature

*ρ*, then , and deduce that equation (2.12) remains valid, after replacing

*x*

_{0}by

*s*

_{0}. A similar error term arises from curvature variations, assumed to occur on the scale

*L*. Use of the outward normal requires us to replace equation (2.22) by(2.34)after which the earlier manipulations establish the general form of the eigenvalue equation (2.29),(2.35)Here, the function

*G*

^{*}is now defined by (more abuse of notation)(2.36)where

*r*_{ip}denote points on

*C*

_{p}with positions given by dimensionless arc length values as defined above. The result (2.35) is analogous to equation (2.7) of Kolokolnikov

*et al*. (2005) for the patch case.

We can recover the basic WK result by finding a crude approximation for *G*^{*}. We note that integrating equation (2.8) over the domain, using the divergence theorem for the first term and then substituting the point source in the boundary term, gives the exact relation . Hence, the leading term in *G*^{*} is *A*^{−1}*k*^{−2}, and not an *O*(1) term, as is well known (cf. e.g. Kolokolnikov *et al*. 2005). Substituting this into equation (2.29) gives equation (2.32), which is the same as equation (2.2). Hence, the basic WK estimate is dangerous because, while it comes from the dominant term in *G*^{*}, the eigenvalue relation (2.29) is approximated by neglecting *O*(1) terms compared with logarithmic terms.

With interest focused on the eigenvalue perturbed from zero, we aim to solve the problem for *G* as a series in *k*^{2}. The form of the Helmholtz operator suggests that we write(2.37)where *A* is the area of *D* and the regular functions *K*^{0} and *K*^{1} satisfy(2.38)The boundary conditions, derived from the low-order terms in the expansion of *Y*_{0}, are(2.39)(2.40)In the first instance, these conditions are applied on *C*−*C*_{p}; their extension to *C* as shown completes the definitions of *K*^{0} and *K*^{1}. Evidently, *K*^{0} is determined by a Poisson problem to within an arbitrary constant, whose value is given by the identity(2.41)obtained by applying the divergence theorem to the second of equations (2.38).

In the half-disc example, equations (2.39) and (2.40) yield(2.42)(2.43)(2.44)with *a*_{n} immaterial, whence(2.45)and equation (2.41) yields(2.46)Substitution in equation (2.35) then gives(2.47)and equation (2.30) is recovered.

### (c) The rectangle

We now apply this technique to the rectangle 0<*y*<*a*, 0<*x*<*b* with a Dirichlet patch at *x*=0, |*y*−*y*_{0}|<*ϵ*. We need the Green's function to the problem with a source within the patch as above. We then need to isolate the logarithmic singularity of the Green's function and the linear correction to it, which will not have a simple closed form as in equation (2.14).

Let *G*(*x*, *y*, *y*′, *k*^{2}) be the Green's function such that(2.48)subject to the Neumann boundary conditions(2.49)(2.50)(2.51)It is assumed that *k*^{2} does not coincide with an eigenvalue and noted that equation (2.50) requires a source of strength 2 at (0, *y*′). Thus, by construction,(2.52)with the constants included to ensure convergence.

Then, evidently,(2.53)Substitution of(2.54)and the definition,(2.55)in equation (2.53) then yields, as described above,(2.56)where the subscripts 1 and 2 denote differentiation with respect to *y* and *y*′, respectively. For the Dirichlet patch, we set *ϕ*(0, *y*_{0}+*ϵs*)=0 in equation (2.56) and deduce that(2.57)because *a*_{1}/*a*_{0}=*O*(*ϵ*) in general. The WK result equation (2.2) gives (*π*/*ab*)[ln(2/*ϵ*)]^{−1}.

The only part of the above procedure that is not entirely straightforward is the calculation of *G*^{*} and its various derivatives. This is analogous to the associated problem of WHK. The function *G*^{*}(*y*_{0}, *y*_{0}, *k*^{2}) is defined by equation (2.55) and can be obtained by solving a Poisson equation analogous to that in the previous section. Once *G*^{*} is known, equation (2.57) can then be solved by a root-finding procedure, since in general the dependence of *G*^{*} on *k* is known only numerically.

One procedure that can be used for the rectangle is a resummation of the series for *G* and then judicious use of the Ewald sum procedure. This is outlined in appendix A. By contrast, the small-*k* approximation method described by equations (2.37)–(2.41) is remarkably successful. The problem for *K*^{0} is first reduced to the elementary construction of a harmonic function that has zero normal derivatives on the three sides *x*=0, *y*=0 and *y*=*a*. The logarithmic terms with *n*=0 in equation (2.52) yieldwhich can be summed by evaluating its *x*- and *y*-derivatives. The additional function of *x*′, determined so that the sum is zero at the origin, is immaterial here since *K*^{0} contains an arbitrary constant. Thus, we have(2.58)in which the log term appears on the left-hand side because the same unwanted term is included in the series, *x*^{2}/2*ab* is a particular solution of ∇^{2}*K*^{0}=(*ab*)^{−1} in equation (2.38) and , subject to at *x*=0 and at *y*=0, *a*. Hence, since the left-hand side of equation (2.58) must have zero *x*-derivative at *x*=*b*,(2.59)The boundary conditions (2.40) on *K*^{1} yield the required jumps(2.60)(2.61)and subsequent substitution of these and equations (2.56) and (2.59) into equation (2.41) gives, after lengthy algebra,(2.62)With *G*^{*} defined by equation (2.55), substitution in equation (2.56) then givessince the ln singularities in *K*^{0} have to cancel. Thus, is revealed the simple result(2.63)Note that, as in equation (2.30), the algebraic form of the estimate provides warning of the method's failure when the patch approaches a corner.

Figure 2 compares the different results. The naive spectral method, which is entirely analogous to that in §2*a*, performs poorly yet again while the finite element method and the asymptotic result are close again. Once again, the WK result is clearly the correct asymptotic limit, but has slow inverse logarithmic decay. Computation times are summarized in table 2; they have the same characteristics as before.

## 3. The gap problem

We now consider problems with gaps between domains. The relevant geometries are shown in figure 3 for §§3*a*,*b* and 4.

### (a) Disc and annulus

Suppose that the disc 0≤*r*<*R*_{1} is connected to the annulus *R*_{1}<*r*<*R*_{2} by the small gap *r*=*R*_{1}, |*θ*|<*ϵ*≪1. The interior Green's function *G*_{i}(*r*, *θ*−*θ*_{0}, *k*^{2}) and the annular Green's function *G*_{a}(*r*, *θ*−*θ*_{0}, *k*^{2}) satisfying(3.1)(3.2)and such that(3.3)(3.4)are given by(3.5)(3.6)The appearance of the periodic delta function ensures that the unit inward (disc or annulus) flux is confined to the point *r*=*R*_{1}, *θ*=*θ*_{0}. The associated logarithmic singularity can be displayed by rewriting equation (3.5) and the first line of equation (3.6) in the forms(3.7)in whichThe second and the last terms in equation (3.7) cancel with the third term to leave equation (3.5), but the sum is now well behaved near the singular point. Similar comments apply to *G*_{a}. The assumption that *k*^{2} is not an eigenvalue of either the disc or the annulus ensures that the denominators in *G*_{i} and *G*_{a} are non-zero.

Suppose that *ϕ*(*r*,*θ*) satisfies equations (3.1), (3.2) and (3.4) and(3.8)Then, evidently,(3.9)so that(3.10)and(3.11)so that(3.12)Matching of equations (3.10) and (3.12) yields(3.13)where, from equations (3.7) and (3.6),(3.14)(3.15)since . The square-root edge singularity is anticipated by setting(3.16)and assuming that *a*_{0}≠0. Then the eigenvalue equation, to leading order in *ϵ* and obtained by substitution in equation (3.13), is(3.17)This method may be readily extended to nested annuli with arbitrarily placed gaps in a manner analogous to that given above for the rectangles. The extension to non-circular boundaries follows that for the single domain in §2*b* and the general form of the eigenvalue equation (2.35) is then replaced by(3.18)consistent with equation (3.17).

For the small-*k* approximation in the circular geometry considered above, the application of equations (2.37)–(2.41) to both the Green's functions yields , , which can be deduced from equation (3.17).

### (b) Circle with barrier

We now consider a circle of radius 1 with a barrier extending across its diameter from −1 to 1−2*ϵ*. Owing to the symmetry of this problem, the lowest eigenvalue of this problem is the same as that for the Dirichlet patch problem above. However, putting *x*_{0}=1−*ϵ* into equation (2.30) indicates a breakdown in the approximation procedure. To overcome this problem, we use the following approach.

Since *x*_{0}=1−*ϵ*, the singularity in *G* is close to the boundary *r*=1. Hence, we exploit the identity,by replacing equation (2.10) by(3.19)in which the additional terms sum to zero. Then,(3.20)Thus, at small *k*,(3.21)because . Substitution of equation (3.21) into equation (2.23), namely(3.22)now yields, at leading order,(3.23)The double integral isbecauseThe remaining integral iswhere(3.24)Since *I*′(*x*)=(arctan *x*)/*x* and *I*(0)=0, it follows that(3.25)where *C*_{G}=0.915… is Catalan's constant. Hence, equation (3.23) yields the estimate(3.26)

Unfortunately, simple finite-element methods, such as the Matlab toolbox, cannot handle domains that are divided into two by a barrier with zero width as in figure 3*a*,*b*. A barrier with small but non-zero width could be used, but then one would not strictly be comparing the same values. More sophisticated software such as PLTMG (Bank 1994) can handle zero-width barriers but we did not pursue this approach. Another option is to use the extension to the method of particular solutions due to Betcke & Trefethen (2005). This approach represents the solution by a sum of Fourier–Bessel terms that have the appropriate inverse square-root behaviour at the tip of the barrier and enforces the correct Neumann condition on the rest of the boundary. To ensure numerical well posedness, additional points in the interior of the domain are used to force the solution away from zero in the interior. However, the number of points on the boundary close to the gap, which is obviously crucial to the success of the method, depends on the total number of points used. Smaller gaps require more points to resolve the geometry of the boundary close to the gap. In practice, the method is beset by the same problems as the naive spectral approach mentioned above.

### (c) Multiple rectangles

Suppose that the *N*+1 rectangles {*R*_{i};0≤*i*≤*N*} lie at *ib*≤*x*≤(*i*+1)*b*, 0≤*y*≤*a* and are connected by gaps at *x*=*ib*, *Y*_{i}−*ϵ*<*y*<*Y*_{i}+*ϵ*(l≤*i*≤*N*). Then, as in equation (2.53),(3.27)(3.28)(3.29)The square-root edge singularity is again anticipated by setting(3.30)for the leading term, and then matching in the gaps yields the asymptotic eigenvalue problem in the form of a tridiagonal system of *N* homogeneous equations for {*a*_{i0};1≤*i*≤*N*}, namely(3.31)

Note that for *N*=1, equation (3.31) just recovers the small-*k* result (2.63) for the Dirichlet patch, as discussed in §1, because *x*=*b* is then a line of symmetry. Hence, since *G*^{*}(*Y*_{i}, *Y*_{i}, *k*^{2}) and *G*(*b*, *Y*_{i}, *Y*_{i±1}, *k*^{2}) have the term (*abk*^{2})^{−1}, the leading-order behaviour for small *k* of the *N* eigenvalues is *λ*_{WK} times the eigenvalues of the tridiagonal matrix with diagonal elements unity and super/subdiagonal elements 1/2. By setting up a recurrence relation, this matrix is found to have characteristic polynomial(3.32)where *U*_{N}(cos *ψ*)=sin(*N*+1)*ψ*/sin *ψ* is a Chebyshev polynomial of the second kind. Thus, the leading-order eigenvalues are just 1−cos(*jπ*/(*N*+1)) (1≤*j*≤*N*) times the basic WK result (2.32) for the Dirichlet patch problems, for any *Y*_{i} not close to 0 or *a*. Evidently, these depend on area but not on *Y*_{i}. However, inclusion of the *O*(1) terms in all elements of equation (3.31) gives *N* approximate eigenvalues that are split off from the set *λ*_{WK}[1−cos(*jπ*/(*N*+1))] (1≤*j*≤*N*). Figure 4 shows typical eigenvalues for multiple rectangles. Extensions to rectangles of different widths and heights are straightforward. Figure 5 shows some example results.

## 4. Diffusion through multiple rectangles

If we revert to the initial-value problem(4.1)and apply the Laplace transform(4.2)then(4.3)The solution is given by equations (3.27)–(3.29), with and −*p* instead of *ϕ* and *k*^{2}, respectively, and the additional term (*N*+1)/*p* on the right-hand side of the first equation. If we then mimic equation (3.30) with *A*_{i0}(*p*) instead of *a*_{i0}, the matching conditions analogous to equation (3.31) are(4.4)and the determination of the poles is the above eigenvalue problem. Since, either from equation (A 1) or its definition, *G* has the term −1/*pab*, the functions {*A*_{i0}(*p*);1≤*i*≤*N*} are regular at *p*=0 and hence the flux through each gap eventually vanishes. The matching conditions yield different equations whose solution is such that(4.5)and subsequent substitution in the integrals for verifies that(4.6)as expected.

The truncated system (4.4) was solved to produce the Laplace-transformed variables *A*_{i0}(*p*). These Laplace transforms were subsequently inverted by following the method of Talbot (1979), yielding the time evolution of the leading-order flux values in the *i*th gap as *πa*_{i0}(*t*), in the notation of §3*c*.

## 5. Matched asymptotic expansions for the half-disc with *x*_{0}=0

Here, we consider the lowest three Fourier modes and match more directly and informally than WK. On defining inner coordinates by(5.1)the leading-order inner field is Laplacian and, to be zero at *ξ*=0 (the gap) and have zero flux across the neighbouring boundary (*η*=0,*π*), must be of the form(5.2)To satisfy equation (2.3) with Neumann boundary conditions, the zeroth-mode outer field issince is small. But equation (5.1) yields(5.3)and so matching is achieved with only *c*_{0}≠0 in equation (5.2), whence the *x*_{0}=0 versions of equation (2.29) and hence equations (2.30) and (2.31) are recovered. These estimates are, for *k*^{2}≪1,and, for ,(5.4)

This result is obtained by the WK method with(5.5)(5.6)where *μ*(*ϵ*) is a scale factor to be determined. The right-hand side of equation (5.5) is the limit of the Laplacian field (5.2) away from the patch. Thus,(5.7)because inclusion of *J*_{0}(*j*_{1m}*r*) is a redundant rescaling of *ϕ*_{0}. Condition (2.4) determines *B*_{1}=*π*[*J*_{0}(*j*_{1m})]^{2}/4 and subsequent matching of the ln *r* terms and the constants in equations (5.5) and (5.7) yieldsand hence equation (5.4). The use of *μ*(*ϵ*) is advantageous and the method is evidently equivalent to the more structured matched asymptotic expansions approach.

Reverting to the direct method, the first-mode outer field that satisfies equation (2.3) with Neumann boundary conditions iswhere is the *m*th zero of and is small. But equation (5.1) yields(5.8)and so matching is achieved with only *c*_{1}≠0 in equation (5.2), whenceIf , the estimate reduces to(5.9)which is an *O*(*ϵ*^{2}) correction, consistent with the term in equation (2.11).

This pattern is continued in higher modes but less simply. The second-mode outer field is approximated bysince is small. Thus, the expansionshows that matching with the inner field requires thatwhence the eigenvalue estimate is(5.10)which is an *O*(*ϵ*^{4}) correction, consistent with the term in equation (2.11). Moreover, the restriction in equation (5.2) to odd functions of *ξ* mandates additional termsin {*ϕ*_{inner},*ϕ*}, respectively, whereThus, the eigenvalue correction is *O*(*ϵ*^{4}) but the inner field generates *O*[*ϵ*^{2}/ln(2/*ϵ*)] correction to the outer field. Higher-order terms in *ϕ*_{inner} involve Mathieu functions.

## 6. Conclusion

Our method gives a dominant role to *G*^{*}, which can be computed more easily than the eigenfunctions, for any piecewise smooth boundary. In particular, we have solved for *G*^{*} numerically for the rectangle. Moreover, we give a straightforward construction for small *k*, which furnishes the desired estimate of the lowest eigenvalue, i.e. the one that is a perturbation from zero.

The numerical solution of the full patch problem is computationally demanding. For values of *ϵ* smaller than can be attained for the patch problem, we can be reassured by the fact that exact solutions that exist for related problems such as that in Titcombe & Ward (2000) show excellent agreement between the exact solution and the eigenvalue obtained from resummed series corresponding to equation (2.29) here.

Bigg & Tuck (1982) use matched asymptotic expansions to calculate Helmholtz resonances in cavities with small openings. Our approach could be applied to that problem.

## Acknowledgments

The authors gratefully acknowledge helpful discussions with Professor William R. Young of the Scripps Institution of Oceanography and with Professor Joseph B. Keller.

## Footnotes

- Received September 22, 2006.
- Accepted November 16, 2006.

- © 2006 The Royal Society