A fast solver for the narrow capture and narrow escape problems in the sphere

by   Jason Kaye, et al.
NYU college

We present an efficient method to solve the narrow capture and narrow escape problems for the sphere. The narrow capture problem models the equilibrium behavior of a Brownian particle in the exterior of a sphere whose surface is reflective, except for a collection of small absorbing patches. The narrow escape problem is the dual problem: it models the behavior of a Brownian particle confined to the interior of a sphere whose surface is reflective, except for a collection of small patches through which it can escape. Mathematically, these give rise to mixed Dirichlet/Neumann boundary value problems of the Poisson equation. They are numerically challenging for two main reasons: (1) the solutions are non-smooth at Dirichlet-Neumann interfaces, and (2) they involve adaptive mesh refinement and the solution of large, ill-conditioned linear systems when the number of small patches is large. By using the Neumann Green's functions for the sphere, we recast each boundary value problem as a system of first-kind integral equations on the collection of patches. A block-diagonal preconditioner together with a multiple scattering formalism leads to a well-conditioned system of second-kind integral equations and a very efficient approach to discretization. This system is solved iteratively using GMRES. We develop a hierarchical, fast multipole method-like algorithm to accelerate each matrix-vector product. Our method is insensitive to the patch size, and the total cost scales with the number N of patches as O(N log N), after a precomputation whose cost depends only on the patch size and not on the number or arrangement of patches. We demonstrate the method with several numerical examples, and are able to achieve highly accurate solutions with 100,000 patches in one hour on a 60-core workstation.



There are no comments yet.


page 5

page 7

page 9

page 27


On the Numerical Solution of Fourth-Order Linear Two-Point Boundary Value Problems

This paper introduces a fast and numerically stable algorithm for the so...

A spectrally accurate method for the dielectric obstacle scattering problem and applications to the inverse problem

We analyze the inverse problem to reconstruct the shape of a three dimen...

A domain decomposition solution of the Stokes-Darcy system in 3D based on boundary integrals

A framework is developed for a robust and highly accurate numerical solu...

A fast direct solver for two dimensional quasi-periodic multilayered medium scattering problems

This manuscript presents a fast direct solution technique for solving tw...

Towards optimal boundary integral formulations of the Poisson-Boltzmann equation for molecular electrostatics

The Poisson-Boltzmann equation offers an efficient way to study electros...

Modeling 3D Surface Manifolds with a Locally Conditioned Atlas

Recently proposed 3D object reconstruction methods represent a mesh with...

Quadrature by Parity Asymptotic eXpansions (QPAX) for scattering by high aspect ratio particles

We study scattering by a high aspect ratio particle using boundary integ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

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 disk-shaped and comment briefly on more general shapes in the conclusion.

Figure 1: A sphere partially covered by disk-shaped patches. We assume each patch is of radius . We also assume that distinct patches are separated by a distance of at least . In the figure, this means that the regions bounded by the dashed lines do not overlap.

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 Dirichlet-Neumann boundary conditions [8, 9]. More precisely, for the capture problem, if the far-field particle concentration is set to be , then satisfies the exterior Laplace equation:


A scalar quantity of interest is the total equilibrium flux of particles through :


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:


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 :


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 :


Thus the average MFPT over may be obtained from the average MFPT on .

Given an arrangement of patches, we present here a fast, high-order 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 Dirichlet-Neumann interfaces [14, 15]. Direct discretization, using either PDE-based 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 ill-conditioned, especially in cases involving large numbers of small patches.

The formulation presented here is well-conditioned, 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 large-scale 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:

  • We make use of the Neumann Green’s functions for the interior and exterior of the sphere to recast (1) and (3) as first-kind integral equations for a density on .

  • Given a patch radius , we precompute the solution operator for the corresponding one-patch 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 many-patch integral equation, we use the solution operator for the one-patch integral equation as a block-diagonal “right preconditioner”. This yields a second-kind Fredholm system of equations which, upon discretization, is well-conditioned and has a small number of degrees of freedom per patch.

  • We solve the resulting linear system by iteration, using GMRES, and accelerate each matrix-vector 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 first-kind 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 first-kind integral equations into Fredholm equations of the second-kind, using the solution operator for the one-patch 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 one-patch integral equation. We demonstrate the performance of the method with numerical experiments in Section 10.

Figure 2: MFPT plotted just inside the unit sphere for an example with random well-separated patches of radius . The integral equation associated with this problem was solved in minutes on a 60-core workstation, to an residual error of approximately . Further details are given in Section 10.2.

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:


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


with as for fixed . It can be shown, using Green’s second identity, that


solves the exterior Neumann problem (6). When , is given explicitly by


If, in addition, , then


The interior Neumann problem is given by


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


As before,


solves the interior Neumann problem (11). When , is given by


If, in addition, , this reduces to


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 ,

Figure 3: MFPT plotted just inside the unit sphere for an example with uniformly distributed patches of radius . The integral equation associated with this problem was solved in seconds on a 60-core workstation, and in minutes on a four-core, eight-thread laptop, to an residual error of approximately . Further details are given in Section 10.2.

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


Let us denote the unknown Neumann data on by . Then (8) implies that for , we have


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 free-space 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 first-kind integral equation


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


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 free-space Green’s function for the Laplace equation (up to a constant scaling factor), we obtain the first-kind integral equation for the Dirichlet problem on an open surface, which we can denote in operator form by

This is a well-studied 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:


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 non-uniqueness 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


for some non-zero constant . Then given by


solves (3). We will therefore seek a method to produce a solution of (21) for some .

Figure 4: MFPT plotted just inside the unit sphere for an example with random, clustered patches of radius . The integral equation associated with this problem was solved in seconds on a 60-core workstation, and in minutes on a four-core, eight-thread laptop, to an residual error of approximately . Further details are given in Section 10.2.
Lemma 2



where satisfies the first-kind integral equation


for . Then solves (21) with , for defined as in (19), and .

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.

To calculate the average MFPT directly from , we plug (22) into (5) to obtain


To calculate , we use the representation (23):

A calculation using the explicit form (15) of gives

for any . We therefore have

Plugging this into (25) and replacing by gives


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


For the sake of simplicity, we assume that each patch has the same radius . We also assume that the patches are well-separated, 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 many-patch integral equation (27) in the form


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 one-patch integral equation with data is defined by


where is an unknown density on .

Remark 3

Writing (28) in the form

and observing that is a smooth function for well-separated from , we see that each satisfies a one-patch integral equation with smooth data. Conversely, if satisfy (28), then each is smooth on .

It is convenient to make use of an orthonormal basis of smooth functions on each patch, so that for smooth on we have


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 one-patch integral equation with data given by the basis element :


Thus, if a smooth function on is expanded as , then the solution of the one-patch integral equation with data is given by . This motivates the following definition.

Definition 3

We denote the solution operator of the one-patch integral equation in the basis by

For and , satisfies

We denote the solution operator of the one-patch integral equation in the truncated basis by

For and , satisfies

Note that the construction of requires solving the one-patch 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 well-conditioned 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


Since the patches and are well-separated, is a compact operator for , so that (32) is a Fredholm system of the second kind. The corresponding truncated system takes the form


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 well-conditioned Fredholm equation of the second kind.

Remark 5

The original unknowns may be recovered from the solution of (32) or (33) using the formula


Thus, we may think of the unknowns as a representation of the unknown density in the basis .

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 well-known 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


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 Gauss-Legendre 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


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 one-patch case, the summation term vanishes, and is sufficient. For multiple patches, the choice of depends largely on how well-separated 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 one-patch solver

While the details of our solver for the one-patch 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 half-space) 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 right-hand side and by

We assume that is discretized to high-order accuracy by a matrix with


so that the discretized system takes the form


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


for any smooth , with high-order accuracy. In the next section, we will use this quadrature to discretize the operators .

The solutions of the one-patch 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 one-patch 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


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 trapezoidal-Legendre 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,


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


Since we may precompute the row vector of length , the cost to compute is .

When the system (41) is solved iteratively, each matrix-vector product is dominated by the computation of the “multiple scattering events”


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 matrix-vector products, the cost would be