1 Introduction
We consider the numerical solution of two related problems which arise in the study of Brownian diffusion by a particle in the exterior or interior of a porous sphere. We denote the open unit ball centered at the origin in by , and assume that the sphere is partially covered by small patches of radius , measured in arclength (Fig. 1). For the sake of simplicity, we assume that the patches are diskshaped and comment briefly on more general shapes in the conclusion.
The union of the patches is referred to as the absorbing boundary and denoted by . The remainder of the boundary, , is referred to as the reflecting boundary. The first problem, called the narrow capture problem, is to calculate the concentration , at equilibrium, of Brownian particles at with a given fixed concentration far from the origin, assuming that particles are absorbed (removed) at . The second problem, called the narrow escape problem, is to calculate the mean first passage time (MFPT) in , namely the expected time for a Brownian particle released at to first reach . In both settings, particles are reflected from . In this paper, we sometimes refer to the narrow capture problem as the exterior problem, and the narrow escape problem as the interior problem.
These problems have received quite a lot of attention in the mathematics and biophysics communities since the seminal work of Berg and Purcell [1]. We do not seek to review the biophysical background here, but note that the absorbing patches serve as a simplified model for either surface receptors (the capture mechanism) or pores (the escape mechanism) in an otherwise impermeable membrane. We refer the reader to [1, 2, 3, 4, 5, 6, 7] for more detailed discussions of applications and a selection of work on related biophysical models.
Standard arguments from stochastic analysis show that both and satisfy a Poisson equation with mixed DirichletNeumann boundary conditions [8, 9]. More precisely, for the capture problem, if the farfield particle concentration is set to be , then satisfies the exterior Laplace equation:
(1) 
A scalar quantity of interest is the total equilibrium flux of particles through :
(2) 
This is sometimes referred to as the capacitance of the system (see Remark 1). For the escape problem, the MFPT satisfies the interior Poisson equation:
(3) 
Here, the quantity of interest is the average MFPT  that is the average, over all possible initial particle positions, of the expected time to escape from through :
(4) 
Here, and in the remainder of the paper, refers to the derivative in the outward normal direction; points towards the interior of for the exterior problem, and towards the exterior of for the interior problem. In order to understand how the distribution of absorbing patches on the surface affects , and the associated quantities and , a variety of asymptotic and numerical methods have been developed (see [1, 10, 11, 12, 13, 4, 5] and the references therein).
Remark 1
The total flux defined in (2) is sometimes referred to as the capacitance because of a connection to electrostatics. Imagine that the ball is a dielectric with low permittivity, and that is a collection of perfectly conducting patches on its surface, connected by infinitesimally thin wires so that they act as a single conductor. Suppose also that this object is surrounded by a dielectric with high permittivity and that the outer dielectric is enclosed by an infinitely large perfectly conducting sphere, with a unit voltage drop from the outer conductor to the conducting patches. Then, letting the ratio of the permittivity of the outer dielectric to that of the inner dielectric approach , the electrostatic potential outside satisfies (1), and the electrostatic capacitance of the system is given by .
Remark 2
The total flux is computed directly from the Neumann data on , as seen from (2). Likewise, the average MFPT can be computed directly from the Dirichlet data on . For this, we use Green’s second identity,
with and . Using that , , and that for , and , we obtain
Applying the divergence theorem to the second term, dividing by , and using that , gives an alternative expression for :
(5) 
Thus the average MFPT over may be obtained from the average MFPT on .
Given an arrangement of patches, we present here a fast, highorder accurate numerical scheme for the evaluation of , , , and , of particular use when is large and is small. Such computations are numerically challenging, partly because solutions of elliptic boundary value problems of mixed type are singular near DirichletNeumann interfaces [14, 15]. Direct discretization, using either PDEbased methods or integral equation methods, would require many degrees of freedom to resolve the singularities in and . Further, the resulting linear systems would be large and illconditioned, especially in cases involving large numbers of small patches.
The formulation presented here is wellconditioned, is nearly identical for the capture and escape problems, and suffers no loss in accuracy or increase in computational cost as is decreased. To make largescale problems practical, we have developed a fast algorithm, so that the cost per GMRES iteration [16] is of the order , rather than . Our method involves the following ingredients:

Given a patch radius , we precompute the solution operator for the corresponding onepatch integral equation, assuming smooth Dirichlet data which is expanded in a rapidly converging series of Zernike polynomials. We analytically incorporate a square root singularity in the induced density at the Dirichlet/Neumann interface.

To solve the manypatch integral equation, we use the solution operator for the onepatch integral equation as a blockdiagonal “right preconditioner”. This yields a secondkind Fredholm system of equations which, upon discretization, is wellconditioned and has a small number of degrees of freedom per patch.

We solve the resulting linear system by iteration, using GMRES, and accelerate each matrixvector product by means of a fast algorithm modeled after the fast multipole method (FMM). The fast algorithm uses the interpolative decomposition [17] to derive a compressed representation of the outgoing field induced by the density on a patch, a hierarchical organization of patches into groups at different length scales, and a spectral representation of the smooth incoming field due to densities on distant patches.
Though most of the past work on the narrow capture and narrow escape problems is based on asymptotics, we wish to highlight the numerical work of Bernoff and Lindsay, who also proposed an integral equation method for the narrow capture problem for the sphere and the plane based on the Neumann Green’s function [12]. Our approach to discretization shares several characteristics with theirs: both methods incorporate a square root singularity into the density on each patch analytically, and both use a representation in terms of Zernike polynomials for smooth Dirichlet data on each patch.
The paper is organized as follows. In Section 2, we introduce the analytical framework for our method, reformulate the boundary value problems as firstkind integral equations using single layer potentials, and explain how to calculate the scalar quantities and directly as functionals of the layer potential densities. In Section 3, we show how to transform the firstkind integral equations into Fredholm equations of the secondkind, using the solution operator for the onepatch integral equation as a preconditioner. In Sections 4, 5, and 6 we describe our discretization approach for the full system of equations, and in Section 7 we introduce the technical tools involved in our fast algorithm. In Section 8 we describe the full method, including our fast algorithm to accelerate the application of the system matrix. In Section 9, we provide a detailed description of the solver for the onepatch integral equation. We demonstrate the performance of the method with numerical experiments in Section 10.
2 Analytical setup
Our approach to solving the exterior and interior problems (1) and (3) uses a representation of each solution as an integral involving the corresponding Neumann Green’s function. This representation leads to an integral equation, and the scalar quantity of interest  or  can be calculated directly from its solution.
2.1 Neumann Green’s functions for the sphere
Let us first consider the exterior Neumann problem:
(6) 
Here is a bounded domain, and a given continuous function on . This problem has a unique solution, and if is the unit ball in , it may be obtained using the exterior Neumann Green’s function , which is known analytically [18, 19]. is symmetric, and satisfies
(7) 
with as for fixed . It can be shown, using Green’s second identity, that
(8) 
solves the exterior Neumann problem (6). When , is given explicitly by
(9) 
If, in addition, , then
(10) 
The interior Neumann problem is given by
(11) 
where is a bounded domain and is a continuous function defined on the boundary, with the additional constraint that must satisfy the consistency condition
This problem has a solution which is unique up to an additive constant. The consistency condition precludes the existence of an interior Green’s function with zero Neumann data. Rather, for the unit ball in , we have an interior Neumann Green’s function , also known analytically [18, 19]. It is again symmetric and satisfies
(12) 
As before,
(13) 
solves the interior Neumann problem (11). When , is given by
(14) 
If, in addition, , this reduces to
(15) 
This is the same as (10) except for the sign of the second term. In other words, the restrictions of the interior and exterior Green’s functions to the boundary are nearly identical.
The following lemma, which we will require in the next section, follows from the second property in (12) and the symmetry of .
Lemma 1
Let be an open subset of and let be continuous on . Then for ,
2.2 The narrow capture problem
We turn now to the narrow capture problem, which is the simpler of the two. We first modify the BVP (1) by defining , so that solutions decay as . The function satisfies the modified equations
(16) 
Let us denote the unknown Neumann data on by . Then (8) implies that for , we have
(17) 
By analogy with classical potential theory, we refer to this as a single layer potential representation with density supported on . Since the dominant singularity of the kernel is that of the freespace Green’s function for the Laplace equation, this single layer potential is continuous up to . Taking the limit as and using the second condition in (16), we obtain the firstkind integral equation
(18) 
where , with the weakly singular kernel . Assuming that we can solve (18) for , it follows that , given by (17), is the solution to (16), and that solves (1). Furthermore, since on , the total flux from (2) will be given by
where we have introduced the shorthand
(19) 
We will not prove the existence of a solution to (18), but sketch a possible approach. If we replace the kernel in (18) with its first term , which is the freespace Green’s function for the Laplace equation (up to a constant scaling factor), we obtain the firstkind integral equation for the Dirichlet problem on an open surface, which we can denote in operator form by
This is a wellstudied problem, which has a unique solution in the Sobolev space given data in [20]. Writing the full single layer potential operator in the form , where is a compact pseudodifferential operator of order , we may rewrite (18) in the form of a Fredholm integral equation of the second kind:
(20) 
Thus, to prove existence and uniqueness for the single patch equation, one can apply the Fredholm alternative to (20). That is, one need only show that the homogenous version of the single patch equation has no nontrivial solutions. This is straightforward to prove when is sufficiently small, since the norm of goes to zero as goes to zero and the corresponding Neumann series converges. We conjecture that the result holds for any .
2.3 The narrow escape problem
The analytical formulation of the narrow escape problem is somewhat more complicated than that of the narrow capture problem, largely because of the nonuniqueness of the interior Neumann problem, but it leads to a similar integral equation. We first recast the Poisson problem (3) as a Laplace problem with inhomogeneous boundary conditions. Assume that satisfies
(21) 
for some nonzero constant . Then given by
(22) 
solves (3). We will therefore seek a method to produce a solution of (21) for some .
Lemma 2
Proof: The function is harmonic in , and by Lemma 1, it satisfies the third condition of (21) with , as long as . Taking to and using the continuity of the single layer potential up to , we find that will satisfy the second condition of (21) as long as satisfies (24).
It remains only to show that if satisfies (24), then . If not, then given by (23) satisfies (21) with , as does the constant function . It follows from Green’s identity that solutions to (21) with the same value of are unique, so we must have . The formula (14) for shows that if , then , so if we have
a contradiction.
The question of the existence of a solution to (24) is
analogous to that for (18), which was discussed in
Section 2.2.
3 A multiple scattering formalism
We have shown that the solutions of the two boundary value problems of interest, as well the associated scalars and , may be obtained by solving (18) and (24), respectively, on the collection of absorbing patches. These integral equations differ only by the sign of one term in their respective kernels, as seen in Section 2.1. Since our treatment of the two cases is the same, we drop the subscripts on and , and discuss the solution of
where is an unknown density on . Letting , where is the th patch, and letting be the restriction of to , we write this equation in the form
(27) 
For the sake of simplicity, we assume that each patch has the same radius . We also assume that the patches are wellseparated, in the sense that the distance between the centers of any two patches in arc length along the surface of the sphere is at least . That is, any two patches are separated by a distance greater than or equal to their own radius. For , we define by
More specifically, we define each such operator in a coordinate system fixed about the center of . Since all the patches have the same radius, the operators are therefore identical, and we denote by . Thus we may rewrite the manypatch integral equation (27) in the form
(28) 
The aim of this section is to reformulate (28) as a Fredholm system of the second kind in an efficient basis.
Definition 1
Let be a smooth function on some patch . The onepatch integral equation with data is defined by
(29) 
where is an unknown density on .
Remark 3
It is convenient to make use of an orthonormal basis of smooth functions on each patch, so that for smooth on we have
(30) 
in the usual sense, with
We postpone until Section 4 a discussion of our particular choice of the basis , which will be constructed using Zernike polynomials. We will denoted by the vector of the first coefficients:
Definition 2
Let be a smooth function on defined by (30), with , computed as above. The projection operators and are defined by
with defined in the same manner for . The synthesis operators and are defined by
and are left inverses of and , respectively.
Finally, we define to be the solution of the onepatch integral equation with data given by the basis element :
(31) 
Thus, if a smooth function on is expanded as , then the solution of the onepatch integral equation with data is given by . This motivates the following definition.
Definition 3
We denote the solution operator of the onepatch integral equation in the basis by
For and , satisfies
We denote the solution operator of the onepatch integral equation in the truncated basis by
For and , satisfies
Note that the construction of requires solving the onepatch integral equations with data to obtain , and that the construction of requires solving the first of these equations. For a fixed patch radius , these solutions are universal and do not depend on the number or arrangement of patches in the full problem.
Given , we are now able to rewrite the integral equation (28) as a wellconditioned Fredholm system of the second kind in the basis . On , we define a function by
Substituting into (28), we have
To transform to the basis , we write in the form and multiply on the left by to obtain
(32) 
Since the patches and are wellseparated, is a compact operator for , so that (32) is a Fredholm system of the second kind. The corresponding truncated system takes the form
(33) 
where we have used the approximation .
Remark 4
We refer to the approach described above as a multiple scattering formalism by analogy with the problem of wave scattering from multiple particles in a homogeneous medium. In the language of scattering theory, one would say that for the th patch, the boundary data is the known data (), perturbed by the potential “scattered” from all other patches, namely . Solving the system (28) corresponds to determining how the collection of uncoupled single patch solutions needs to be perturbed to account for the “multiple scattering” effects.
The approach developed above, where are the unknowns, has many advantages over solving (28) directly, even with as a left preconditioner. By working in the spectral basis, we avoid the need to discretize on each patch, the number of degrees of freedom per patch is significantly reduced, and the linear system is a wellconditioned Fredholm equation of the second kind.
Remark 5
We turn now to the construction of an orthonormal basis for smooth functions on a patch, the construction of the singular solutions , and the efficient solution of the discretized multiple scattering system (33).
4 A basis for smooth functions on a patch
It is wellknown that the Zernike polynomials are a spectrally accurate, orthogonal basis for smooth functions on the disk. For a thorough discussion of these functions, we refer the reader to [21]. Here, we simply summarize their relevant properties.
The Zernike polynomials on the unit disk , are given by
with , , and
where is a Jacobi polynomial on . The Jacobi polynomials are orthogonal on with respect to the weight function . Thus, for fixed , the functions are orthogonal on with respect to the weight function . This gives the orthogonality relation
(35) 
The natural truncation of this basis is to fix a cutoff mode in both the radial and angular variables, and to let . This yields basis functions. To use this basis on a generic patch , we define a polar coordinate system about the patch center, for which is the distance in arc length along the sphere from the center, and is the polar angle. We rescale the radial variable from to , transforming the Zernike polynomials to functions on . Finally, the basis functions discussed in Section 3 can be defined as the scaled Zernike polynomials up to mode .
¿From the orthogonality relation (35), the projection operators and are obtained as normalized inner products against Zernike polynomials in polar coordinates. This Zernike transform
can be implemented numerically using a tensor product quadrature with a GaussLegendre rule in the radial variable and a trapezoidal rule in the angular variable. The number of grid points required to obtain the exact Zernike coefficients of a function in the space spanned by
is ; we denote this number by . We refer to these points as the Zernike sampling nodes (see [21] for further details).Remark 6
Rewriting (33) in the form
(36) 
we see that the truncation error compared with (32) depends on how well the smooth function
is represented in the space spanned by . In the onepatch case, the summation term vanishes, and is sufficient. For multiple patches, the choice of depends largely on how wellseparated the patches are. Since the Zernike basis is spectrally accurate, grows only logarithmically with the desired precision. In practice, a posterioriestimates are easily obtained for any fixed configuration by inspection of the decay of the Zernike coefficients in the computed solution.
5 Informal description of the onepatch solver
While the details of our solver for the onepatch integral equation
are deferred to Section 9, we outline the general approach here. First, we note that in the absence of curvature (i.e. a flat disk on a halfspace) and with the associated terms of the Green’s function removed, the solution is known to have a square root singularity at the disk edge [12, 14, 15, 20, 22]. In our case, we will explicitly include this square root singularity in the representation of , but also allow for weaker singularities  which we have observed and will demonstrate in Section 9.3  by using a discretization that is adaptively refined toward the edge .
Assume then that we have discretized the patch using a suitable polar mesh with fine grid points, denoted by . The fine grid points for different patches are identical relative to the coordinate systems of their own patches. We denote the corresponding samples of the righthand side and by
We assume that is discretized to highorder accuracy by a matrix with
(37) 
so that the discretized system takes the form
(38) 
We will also require a set of quadrature weights, denoted by and identical for each patch, that permit the accurate integration over of the product of an arbitrary smooth function with the discretized density , taking into account the fact that has an edge singularity. That is, we assume that
(39) 
for any smooth , with highorder accuracy. In the next section, we will use this quadrature to discretize the operators .
The solutions of the onepatch integral equations (31) may be obtained in a precomputation, after which we have access to the functions sampled on the fine grid. We assemble these functions into an matrix with
is then the discretization of the operator , mapping the first Zernike coefficients of a smooth function to the solution of the corresponding onepatch integral equation sampled on the fine grid. If we denote by the discretization of the synthesis operator as an matrix,
then we have, as in Definition 3,
In short, the precomputation amounts to solving this matrix system for .
6 Discretization of the multiple scattering system
We return now to the multiple scattering system (33). The unknowns on are defined in the truncated Zernike basis as . We will need as intermediate variables the fine grid samples of . From Remark 5, we define the sampling vector by
In order to discretize the integral operators for , we note that is smooth for , , and use the quadrature (39). This yields
(40) 
Setting to be the th Zernike sampling node on , we define the matrix by
Thus, maps a density sampled on the fine grid on to the smooth field it induces at the Zernike sampling nodes on . Lastly, we discretize the truncated Zernike transform as a matrix using the trapezoidalLegendre scheme described in Section 4.
Definition 4
The discrete Zernike transform is defined to be the mapping of a smooth function sampled on the Zernike sampling nodes to its Zernike coefficients.
We can now write the multiple scattering system (33) in a fully discrete form,
(41) 
where is the vector of length with all entries equal to . Since , , and , this is a linear system of dimensions , with degrees of freedom per patch. As a discretization of a Fredholm system of the second kind, it is amenable to rapid solution using an iterative method such as GMRES [16].
We now describe how to calculate the constants and from the solution of (41). We saw in Sections 2.2 and 2.3 that these can be computed directly from . Using the fine grid quadrature (39), we have
(42) 
Since we may precompute the row vector of length , the cost to compute is .
When the system (41) is solved iteratively, each matrixvector product is dominated by the computation of the “multiple scattering events”
(43) 
for . That is, for each patch , we must compute the Zernike coefficients of the field induced on that patch by the densities on all other patches. Note that if we were to calculate the above sums by simple matrixvector products, the cost would be
Comments
There are no comments yet.