Mean exit time for diffusion on irregular domains
Matthew J Simpson, Daniel J Vandenheuvel, Joshua M Wilson, Scott W McCue, Elliot J Carr
MMean exit time for diffusion on irregular domains
Matthew J. Simpson ∗ , Daniel J. Vandenheuvel , Joshua M. Wilson ,Scott W. McCue , and Elliot J. Carr School of Mathematical Sciences, Queensland University of Technology,Brisbane, Queensland 4001, AustraliaFebruary 9, 2021
Abstract
Many problems in physics, biology, and economics depend upon the durationof time required for a diffusing particle to cross a boundary. As such, calculationsof the distribution of first passage time, and in particular the mean first passagetime, is an active area of research relevant to many disciplines. Exact results forthe mean first passage time for diffusion on simple geometries, such as lines, discsand spheres, are well–known. In contrast, computational methods are often used tostudy the first passage time for diffusion on more realistic geometries where closed–form solutions of the governing elliptic boundary value problem are not available.Here, we develop a perturbation solution to calculate the mean first passage timeon irregular domains formed by perturbing the boundary of a disc or an ellipse.Classical perturbation expansion solutions are then constructed using the exactsolutions available on a disc and an ellipse. We apply the perturbation solutions to ∗ To whom correspondence should be addressed. E-mail: [email protected] a r X i v : . [ phy s i c s . b i o - ph ] F e b ompute the mean first exit time on two naturally–occurring irregular domains: amap of Tasmania, an island state of Australia, and a map of Taiwan. Comparingthe perturbation solutions with numerical solutions of the elliptic boundary valueproblem on these irregular domains confirms that we obtain a very accurate solutionwith a few terms in the series only. MATLAB software to implement all calculationsis available at https://github.com/ProfMJSimpson/Exit time. Keywords: Random walk; Hitting time; Passage time; Perturbation.2
Introduction
Mathematical models describing diffusive transport phenomena are fundamental toolswith broad applications in many areas of physics [1–3], engineering [4,5] and biology [6–8].Traditional analysis of mathematical models of diffusion often focus on idealised problemswith relatively simple geometries, whereas practical applications routinely involve morecomplicated geometries that are normally dealt with using computational approaches.While computational approaches are essential in many circumstances [9,10], exact analyt-ical insight is preferable, where possible, because it can lead to mathematical expressionsthat explicitly highlight key relationships that are not always revealed computationally.A fundamental property of diffusion is the concept of particle lifetime, which is a par-ticular application of the more general concept of the first passage time [1–3]. Estimatesof particle lifetime provide insight into the timescale required for a diffusing particle toreach a certain target, such as an absorbing boundary. Many results are known aboutthe particle lifetime for diffusion in simple geometries, such as lines and discs [1, 2]. Gen-eralising these results to deal with other geometrical features, such as wedges [11, 12],symmetric domains [13–15], growing domains [16, 17], slender domains [18–20], smalltargets [21, 22] or arbitrary initial conditions [23] is an active area of research.In this work, we develop new expressions for particle lifetime for diffusion in irregulartwo–dimensional (2D) geometries. Our approach is to use classical results for diffusionon a disc and an ellipse, and then construct perturbation solutions for more generalgeometries [24, 25]. The new perturbation solutions are evaluated using symbolic com-putation, and the results compare very well with numerical solutions of the governingboundary value problem, and with averaged data from repeated stochastic simulations.Once the new perturbation solutions are derived and validated, we consider two casestudies to show how these tools can be used to take a particular 2D region and representit as a perturbed ellipse or disc. This allows us to use the new perturbation solutions to3nalyse the mean particle lifetime in irregular 2D geometries. MATLAB code providedon https://github.com/ProfMJSimpson/Exit time is used to compute: (i) the symbolicevaluation of the perturbation solution; (ii) the numerical finite volume solution of theboundary value problem; and, (iii) the stochastic random walk algorithm.
Standard arguments that relate unbiased random walk models on domains with absorbingboundary conditions can be used to show that the mean exit time is given by the solutionof a linear ellipse partial differential equation [1–3]. In this work we seek solutions of thatequation, D ∇ T ( x ) = − , x ∈ Ω , with (1) T = 0 , on ∂ Ω , (2)where D > T ( x ) is the mean exit time for a particle released atlocation x . Here the diffusivity is related to the random walk model, D = P δ / (4 τ ),where δ > τ > P ∈ [0 , δ in the time duration τ [1–3]. The continuum description is valid in theconstrained limit δ → τ → δ /τ held constant [1–3]. The stochastic randomwalk model is Pearson’s walk in R , and full details of how simulations are performed aregiven in Appendix A.Key results are presented and discussed in the following format. In Sections 2.1 and2.2 we review known exact solutions to Equation (1) on a disc and an ellipse, respectively.In Sections 2.3 and 2.4 we develop new approximate solutions on irregular domains, andin Section 2.5 we apply the new approximate solutions to study the exit time for diffusion4n two naturally–occurring geometries. For the special case where Ω is a disc of radius
R > Dr dd r (cid:20) r d T ( r )d r (cid:21) = − , < r < R. (3)The solution of Equation (3) with appropriate boundary conditions, T ( R ) = 0 andd T (0) / d r = 0, is given by T ( r ) = R − r D . (4)Results in Figure 1(a)–(c) provide a visual comparison of exit time data from repeatedstochastic simulations, the exact solution and the numerical solution of Equations (1)–(2) for a unit disc, R = 1. A complete description of the continuous space, discrete timestochastic algorithm and the finite volume numerical method used with an unstructuredtriangular mesh to solve Equation (1) is given in Appendix A.Figure 1: Mean exit time on the unit disc. (a) Averaged data from repeated stochas-tic simulations. (b) Exact solution of Equations (1)–(2). (c) Numerical solution of Equa-tions (1)–(2). Parameters are τ = 1, P = 1, δ = 1 × − and D = 2 . × − . Thetriangular mesh used to construct the solution in (a) and (c) has 632 nodes and 1183triangular elements. For (a), 1000 random walks starting from each node were generated.5 .2 Mean exit time on an ellipse For the special case where ∂ Ω is an ellipse centered at the origin it is convenient to workin Cartesian coordinates. The ellipse, given by x a + y b = 1 , (5)has width 2 a > b >
0. The solution of Equation (1) on this domain isgiven by [24] T ( x, y ) = a b D ( a + b ) (cid:20) − x a − y b (cid:21) . (6)Results in Figure 2(a)–(c) provide a visual comparison of exit time data from repeatedstochastic simulations, the exact solution and the numerical solution of Equation (1) foran ellipse with a = 2 and b = 1.Figure 2: Mean exit time on an ellipse with a = 2 and b = 1 . (a) Averaged datafrom repeated stochastic simulations. (b) Exact solution of Equation (1). (c) Numericalsolution of Equations (1)–(2). Parameters are τ = 1, P = 1, δ = 1 × − and D =2 . × − . The triangular mesh used to construct the solution in (a) and (c) has 1240nodes and 2356 triangular elements. For (a), 1000 random walks starting from each nodewere generated. We begin working on irregular domains by calculating expressions for the exit time ona perturbed disc. Using plane polar coordinates ( r, θ ), we consider a region Ω withboundary r = R ( θ ), subject to the condition that R ( θ ) is a single-valued function of θ toensure that any ray drawn from the origin intersects precisely one point of the boundary6 Ω. If our region is conceived as a modest perturbation of a circular disc of radius R wecan write R ( θ ) = R (1 + εg ( θ )) , (7)where R > θ ∈ [0 , π ) is the polar angle, ε (cid:28) g ( θ ) is an O (1) smooth periodic function withperiod 2 π . We assume that the solution can be written as T ( r, θ ) = T ( r, θ ) + εT ( r, θ ) + ε T ( r, θ ) + · · · + ε n T n ( r, θ ) + O ( ε n +1 ) . (8)When the O ( ε n +1 ) term is neglected in Equation (8) the solution is referred to as the O ( ε n ) perturbation solution. To proceed we evaluate T ( r, θ ) on R ( θ ) and expand about ε = 0 to give,0 = T ( R + εRg ( θ ) , θ ) , T ( R, θ ) + εRg ( θ ) ∂T∂θ (cid:12)(cid:12)(cid:12)(cid:12) r = R + · · · + ( εRg ( θ )) n n ! ∂ n T∂θ n (cid:12)(cid:12)(cid:12)(cid:12) r = R + O ( ε n +1 ) . (9)Substituting Equation (8) into Equation (9) and equating powers of ε leads to, O (1) : T ( R, θ ) = 0 , (10) O ( ε (cid:96) ) : T (cid:96) ( R, θ ) + (cid:96) (cid:88) k =1 ( Rg ( θ )) k k ! ∂ k T (cid:96) − k ∂r k (cid:12)(cid:12)(cid:12)(cid:12) r = R = 0 , (11)for (cid:96) = 1 , . . . , n . With this information we substitute Equation (8) into Equation (1)to give a family of boundary value problems for each term in the power series, withthe boundary conditions given by Equations (10)–(11). This family of boundary value7roblems can be summarised as O (1) : D ∇ T = − , T ( R, θ ) = 0 , (12) O ( ε (cid:96) ) : ∇ T (cid:96) = 0 , T (cid:96) ( R, θ ) = − (cid:96) (cid:88) k =1 ( Rg ( θ )) k k ! ∂ k T (cid:96) − k ∂r k (cid:12)(cid:12)(cid:12)(cid:12) r = R , (13)for (cid:96) = 1 , . . . , n . The solution of Equation (12) is Equation (4), and each higher orderterm, T (cid:96) ( r, θ ), (cid:96) = 1 , . . . , n , is the solution of Laplace’s equation on a disc with prescribedboundary data associated with the previous term. Each higher order term can be obtainedusing separation of variables, giving T (cid:96) ( r, θ ) = a ∞ (cid:88) k =1 (cid:2) a k r k cos( kθ ) + b k r k sin( kθ ) (cid:3) , (14)where a k = 1 πR k (cid:90) π T (cid:96) ( R, θ ) cos( kθ ) d θ, b k = 1 πR k (cid:90) π T (cid:96) ( R, θ ) sin( kθ ) d θ, (15)for (cid:96) = 1 , . . . , n . These solutions are straightforward to evaluate symbolically, and aMATLAB implementation of the symbolic evaluation of them is available on GitHub.Results in Figure 3(a)–(c) provide a visual comparison of exit time data from re-peated stochastic simulations, the perturbation and numerical solution of Equation (1)for a perturbed unit disc with R = 1, g ( θ ) = sin(3 θ ) + cos(5 θ ) − sin θ , and ε = 1 / ε , and different choices of g ( θ ) can be evaluated using theMATLAB algorithms available on GitHub. Additional information about how the choiceof truncation in Equation (8) affects the accuracy of the perturbation solution is given inthe Appendix B. 8igure 3: Mean exit time on a perturbed disc. (a) Averaged data from repeatedstochastic simulations. (b) O ( ε ) perturbation solution. (c) Numerical solution of Equa-tions (1)–(2). Parameters are τ = 1, P = 1, δ = 1 × − and D = 2 . × − . Thetriangular mesh used to construct the solution in (a) and (c) has 636 nodes and 1189triangular elements. For (a), 1000 random walks starting from each node were generated.It is worth explicitly pointing out some interesting features that arise when we solveEquation (1) on a perturbed disc. In sectors where g ( θ ) >
0, there are points of Ω thatlie beyond the circle r = R , whereas in sectors where g ( θ ) <
0, there are points that liewithin the circle r = R but are not within Ω. This situation often arises when solvingboundary value problems where the location of the boundary is perturbed [24, 26]. Thekey point being that the domain of r in the functions T ( r, θ ) and T i ( r, θ ) for i = 1 , , , . . . , is the same, r < R (1 + εg ( θ )). However, the boundary value problems that definethe terms in the perturbation solution, T i ( r, θ ) for i = 1 , , , . . . , are defined on theoriginal unperturbed domain, r < R . Accordingly, our finite volume calculations andstochastic simulations are performed directly on the perturbed domain, r < R (1 + εg ( θ )),whereas the perturbation solution is calculated on the unperturbed domain r < R withthe boundary conditions on r = R (1 + εg ( θ )) projected onto the unperturbed circle r = R . It is precisely this feature of the perturbation approach that allows us to constructsuch approximate solutions. The same situation arises when we solve Equation (1) on aperturbed ellipse, as we shall now explain. 9 .4 Mean exit time on a perturbed ellipse We proceed by deriving expressions for the exit time on a perturbed ellipse by considering ∂ Ω as x = a (1 + εg ( θ )) cos θ, (16) y = b (1 + εh ( θ )) sin θ, (17)where 2 a > b > a > b and g ( θ ) and h ( θ ) are O (1) smooth periodic functions with period 2 π . Workingin Cartesian coordinates, we suppose the solution takes the form T ( x, y ) = T ( x, y ) + εT ( x, y ) + ε T ( x, y ) + · · · + ε n T n ( x, y ) + O ( ε n +1 ) . (18)To proceed, we impose Equation (2) at the boundary specified in Equation (16) and (17)and expand about ε = 0,0 = T ( a cos θ + aεg ( θ ) cos θ, b sin θ + bεh ( θ ) sin θ )0 = T ( a cos θ, b sin θ ) + n (cid:88) k =1 k ! k (cid:88) i =1 (cid:18) ki (cid:19) ∂ k T∂x i ∂y k − i ( aεg ( θ ) cos θ ) i ( bεh ( θ ) sin θ ) k − i + O ( ε n +1 ) , (19)where we evaluate the partial derivatives on the boundary of the unperturbed ellipse: x = a cos θ , and y = b sin θ . From this point on in this Section we evaluate all par-tial derivatives like this on the boundary of the unperturbed ellipse. For brevity it is10onvenient to define H (cid:96) ( θ ) = − ∂T (cid:96) − ∂y bh ( θ ) sin θ − ∂T (cid:96) − ∂x ag ( θ ) cos θ − ∂ T (cid:96) − ∂y ( bh ( θ ) sin θ ) − ∂ T (cid:96) − ∂x∂y ( ag ( θ ) cos θ )( bh ( θ ) sin θ ) − · · · − (cid:96) ! ∂ (cid:96) T ∂x (cid:96) ( ag ( θ ) cos θ ) (cid:96) = − (cid:96) (cid:88) k =1 k (cid:88) i =0 k ! (cid:18) ki (cid:19) ∂ k T (cid:96) − k ∂x i ∂y k − i ( ag ( θ ) cos θ ) i ( bh ( θ ) sin θ ) k − i . (20)Substituting Equation (18) into Equation (19), and equating powers of ε gives O (1) : T ( a cos θ, b sin θ ) = 0 , (21) O ( ε (cid:96) ) : T (cid:96) ( a cos θ, b sin θ ) = H (cid:96) ( θ ) , (22)for (cid:96) = 1 , . . . , n . Substituting Equation (18) into Equation (1) leads to the followingfamily of boundary value problems O (1) : D ∇ T = − , T ( a cos θ, b sin θ ) = 0 , (23) O ( ε (cid:96) ) : ∇ T (cid:96) = 0 , T (cid:96) ( a cos θ, b cos θ ) = H (cid:96) ( θ ) , (24)for (cid:96) = 1 , . . . , n . The solution of Equation (23) is Equation (6), and each higher orderterm is the solution of Laplace’s equation on the ellipse with prescribed boundary dataassociated with the previous terms. For example, the O ( ε (cid:96) ) boundary value problem givenin Equation (24) involves the boundary data H (cid:96) ( θ ), which depends on T (cid:96) − , T (cid:96) − , . . . , T as evident in Equation (20). Following Jackson [27] we construct the solution in terms ofharmonic polynomials by expanding the boundary data H (cid:96) ( θ ) as a Fourier series, H (cid:96) ( θ ) = a ∞ (cid:88) k =1 ( a k cos( kθ ) + b k sin( kθ )) , ≤ θ < π, where , (25) a k = 1 π (cid:90) π H (cid:96) ( θ ) cos( kθ ) d θ, b k = 1 π (cid:90) π H (cid:96) ( θ ) sin( kθ ) d θ. (26)11e then compute u (cid:96) ( x, y ) = (cid:60) (cid:0) ( x + iy ) (cid:96) (cid:1) and v (cid:96) ( x, y ) = (cid:61) (cid:0) ( x + iy ) (cid:96) (cid:1) for (cid:96) = 1 , , . . . ,and also define u ( x, y ) = 1 and v ( x, y ) = 0, given explicitly as u (cid:96) ( x, y ) = (cid:98) (cid:96) (cid:99) (cid:88) j =0 (cid:18) (cid:96) j (cid:19) x (cid:96) − j ( − j y j , (27) v (cid:96) ( x, y ) = (cid:98) (cid:96) +12 (cid:99) (cid:88) j =0 (cid:18) (cid:96) j + 1 (cid:19) x (cid:96) − j − ( − j y j +1 . (28)The next step is to compute 2 T k (cos θ ) for k = 1 , , . . . , where T k is the Chebyshevpolynomial of the first kind of degree k , and extract the coefficients of cos r θ , C r , in thispolynomial. Using these expressions we compute, for even k ≥ U k ( x, y ) = 1( a + b ) k + ( a − b ) k k (cid:88) r =0 C r (cid:0) a − b (cid:1) k − r u r ( x, y ) , (29) V k ( x, y ) = 1( a + b ) k − ( a − b ) k k (cid:88) r =0 C r (cid:0) a − b (cid:1) k − r v r ( x, y ) , (30)and for odd k ≥ U k ( x, y ) = 1( a + b ) k + ( a − b ) k k − (cid:88) r =0 C r +1 (cid:0) a − b (cid:1) k − − r u r +1 ( x, y ) , (31) V k ( x, y ) = 1( a + b ) k − ( a − b ) k k − (cid:88) r =0 C r +1 (cid:0) a − b (cid:1) k − − r v r +1 ( x, y ) . (32)This gives us the solution of our O ( ε (cid:96) ) problem, T (cid:96) ( x, y ) = a ∞ (cid:88) k =1 ( a k U k ( x, y ) + b k V k ( x, y )) , (33)where a , a k and b k are the Fourier coefficients from the boundary data, Equation (25)–Equation (26).This solution is straightforward to evaluate symbolically, and a MATLAB implemen-tation to evaluate it is available on GitHub. Results in Figure 4(a)–(c) provide a visual12omparison of exit time data from repeated stochastic simulations, the perturbation solu-tion and the numerical solution of Equation (1) for a perturbed ellipse with a = 2, b = 1, g ( θ ) = sin(3 θ ) + cos(5 θ ) − sin θ , h ( θ ) = cos(3 θ ) + sin(5 θ ) − cos θ , and ε = 1 /
20. Thesolutions in Figure 4 show that the truncated perturbation solution can be very accurate,with just three terms in Equation (18), and 25 terms in Equation (33) required to producea good match with the numerical solution. Additional terms in the perturbation solution,or perturbation solutions for different choices of ε , g ( θ ), or h ( θ ) can be evaluated usingthe MATLAB algorithms on GitHub.Figure 4: Exit time on a perturbed ellipse. (a) Averaged data from repeated stochas-tic simulations. (b) O ( ε ) perturbation solution. (c) Numerical solution of Equations(1)–(2). Parameters are τ = 1, P = 1, δ = 1 × − and D = 2 . × − . The triangularmesh used to construct the solution in (a) and (c) has 1243 nodes and 2342 triangularelements. For (a), 1000 random walks starting from each node were generated. To showcase how the perturbation solutions can be applied to naturally–occurring ge-ometries, we now turn to the problem of taking an image of a region, approximating thatregion as a perturbed disc or perturbed ellipse, and then using our approach to estimatethe exit time from that region. For this exercise we consider the islands of Tasmania andTaiwan. In particular we represent the boundary of Tasmania as a perturbed disc andthe boundary of Taiwan as a perturbed ellipse, based on the maps in Figure 5.13 a) (b)N N
Figure 5:
Case studies. (a) Tasmania [28]. (b) Taiwan [29].Note that we have deliberately omitted any scale on the maps in Figure 5 since wewish to focus on the role of shape rather than size. This allows us to represent theboundary of Tasmania as a perturbed unit disc, and to compare the exit time from thisrealistic geometry with the exit times from the more idealised shapes in Figure 1 andFigure 3. Similarly, we represent the boundary of Taiwan as a perturbed ellipse on ascale that allows us to compare the exit time distribution from this realistic geometrywith the results in Figure 2 and Figure 4. We now explain how we process the images inFigure 5 to extract data that allows us to compute the exit times.To represent Tasmania as a perturbed disc we follow a two-step procedure. First, weuse MATLAB image processing tools to describe the boundary of Tasmania as a set ofpoints as described in the Appendix C, and we fit the disc ( x − x c ) + ( y − y c ) = R tothose points using a spline approximation in MATLAB’s cscvn [30] function. Second, weshift this disc so that it is centered at the origin, and assume that the boundary takes14he form R ( θ ) = R (1 + εg ( θ )), where we approximate g ( θ ) by g ( θ ) = A + G (cid:88) n =1 ( A n cos( nθ ) + B n sin( nθ )) . (34)If the points { ( x i , y i ) } Ni =1 represent the given boundary, we compute the polar angle foreach point θ i . To represent our boundary in this way we require1 R (cid:113) x i + y i − εA + ε G (cid:88) n =1 A n cos( nθ i ) + ε G (cid:88) n =1 B n sin( nθ i ) , i = 1 , , . . . , N. (35)We estimate the coefficients A , A n , B n for n = 1 , , . . . , G by computing the least–squaressolution of Equation (35) that minimises the sum of squared residuals N (cid:88) i =1 (cid:32) R (cid:113) x i + y i − − εA − ε G (cid:88) n =1 A n cos( nθ i ) − ε G (cid:88) n =1 B n sin( nθ i ) (cid:33) . (36)This least–squares solution is computed using MATLAB’s backslash operator. For sim-plicity we set ε = 1 /
10 in this case study. Given our estimates of the coefficients inEquation (35), our approximation of the boundary is shown in Figure 6, where we notethat the approximation amounts to neglecting fine-scale features of the Tasmanian coast-line. For clarity we refer to this region as pseudo-Tasmania . As we pointed out previously,our approach requires that R ( θ ) is a single-valued function of θ such that any ray drawnfrom the origin intersects precisely one point of the boundary ∂ Ω. This condition canonly be met by neglecting the fine-scale structures of the boundary, especially at theSouth-East part of Tasmania where there is an obvious peninsula. Given this approxi-mation, we apply the perturbation analysis in Section 2.3 to give the results in Figure6. 15igure 6:
Mean exit time on pseudo-Tasmania. (a) Numerical solution of Equations(1)–(2). (b) O ( ε ) perturbation solution with G = 3. (c) Discrepancy between thenumerical and perturbation solution. All results correspond to D = 2 . × − . Thetriangular mesh used to construct the solution in (a) has 1188 nodes and 2250 triangularelements.Figure 6(a) shows the numerical solution of Equation (1) within the region enclosedby the boundary obtained by truncating Equation (35) with G = 3 with boundary coordi-nate data obtained from the map in Figure 5(a). All MATLAB files required to replicatethe boundary extraction and fitting are available on GitHub. Figure 6(b) shows the O ( ε )perturbation solution, where infinite sums are approximated using the first 25 terms inEquation (14). Visual comparison of the numerical and perturbation solutions indicatesthat the perturbation solution is remarkably accurate given that the domain in Figure6(a)–(b) is quite far from a unit disc. In fact, the maximum difference between the bound-ary in Figure 6(a)–(b) and the underlying unit disc is, approximately, 0 .
64, confirmingthat this domain is a reasonably large perturbation of a unit disc. Careful comparison ofthe numerical and perturbation solutions show some discrepancy, particularly near thesouthern and north eastern portions of the boundary.To quantify the discrepancy we introduce a measure of the difference between the twosolutions, e ( x, y ) = 100 | T n ( x, y ) − T p ( x, y ) | max ( x,y ) ∈ Ω | T n ( x, y ) | , (37)where T n ( x, y ) is the numerical finite volume solution and T p ( x, y ) is the truncated per-turbation solution, such that e ( x, y ) is a convenient measure of the percentage relative16rror as a function of position. The plot of e ( x, y ) in Figure 6(c) confirms that theperturbation solution is remarkably accurate in the interior of the domain, with smalldiscrepancies along some of the boundary. While the small discrepancy along some partsof the boundary are visually discernable, these differences are not overwhelming, and thebasic features of the numerical solution is captured by the perturbation solution. Moreaccurate perturbation solutions can be constructed by including more terms in the trun-cated perturbation solution, including more terms in the infinite sums, or both. Theseoptions may be explored using the MATLAB algorithms provided on GitHub. More de-tails about the implications of approximating Tasmania by pseudo-Tasmania are givenin Appendix D.To represent Taiwan as a perturbed ellipse we first follow image pre-processing out-lined in the Appendix C. The method developed by Szpak et al. [31] allows us to ap-proximate the boundary as an ellipse with a particular orientation and centre. We thenrotate and shift this identified ellipse so that it is centered at the origin with semi–majoraxis along the x –axis. To approximate the boundary of Taiwan as a perturbed ellipse,Equation (16) and (17), we represent the functions g ( θ ) and h ( θ ) as g ( θ ) = A + G (cid:88) n =1 ( A n cos( nθ ) + B n sin( nθ )) , (38) h ( θ ) = C + H (cid:88) n =1 ( C n cos( nθ ) + D n sin( nθ )) . (39)As before, the boundary is given by a set of points, { ( x i , y i ) } Ni =1 , and we compute thepolar angle θ i for each point. Representing the boundary using this approach leads totwo systems of linear equations x i a cos θ i − εA + ε G (cid:88) n =1 A n cos( nθ i ) + ε H (cid:88) n =1 B n sin( nθ i ) , i = 1 , , . . . , N, (40) y i b sin θ i − εC + ε H (cid:88) n =1 C n cos( nθ i ) + ε H (cid:88) n =1 D n sin( nθ i ) , i = 1 , , . . . , N. (41)17s for Tasmania, we estimate the coefficients, A , A n , B n for n = 1 , . . . , G in Equa-tion (40) and C , C n , D n for n = 1 , . . . , H in Equation (41) by computing the least–squares solution of each linear system using MATLAB’s backslash operator. In summary,we can represent the boundary of Taiwan using g ( θ ) and h ( θ ) given by Equation (16)and 17 for some choice of ε , which we again take to be ε = 1 /
10 in this case study.Figure 7(a) shows the numerical solution of Equation (1) within the region enclosed bythe boundary obtained by truncating Equation (38)–(39) with G = H = 9 that is basedon boundary data obtained from the map in Figure 5(b). All MATLAB files required toreplicate the boundary extraction and fitting are available on GitHub. The solution inFigure 7(b) is the O ( ε ) perturbation solution. Visual comparison of the numerical andperturbation solutions indicates that the perturbation solution is remarkably accurate,and the plot of e ( x, y ), given by Equation (37) in Figure 7(c) confirms this accuracy,even along the boundaries. Again, this accuracy is obtained through neglecting the veryfine-scale features of the coastline of Taiwan that would never be accurately representedby a perturbed ellipse.Figure 7: Mean exit time on pseudo-Taiwan. (a) Numerical solution of Equations(1)–(2). (b) O ( ε ) perturbation solution. (c) Discrepancy between the numerical andperturbation solution. Parameters are τ = 1, P = 1, δ = 1 × − , D = 2 . × − , G = H = 9. The triangular mesh used to construct the solution in (a) has 961 nodes and1803 triangular elements. 18 Conclusions and outlook
In this work we consider the canonical problem of determining the mean first passagetime for diffusion, which requires the solution of an elliptic partial differential equationon the domain of interest. This problem has been studied, in detail, both analyticallyand computationally, with many known exact solutions for relatively simple geometries,such as lines, discs and spheres [1–3]. The calculation of exact expressions for the meanfirst passage time for more complicated geometries is an active, and ongoing field ofresearch. In this work we present new solutions for the mean first passage time fordiffusion on irregular two–dimensional domains, where these solutions are obtained interms of a perturbation of the classical results for the mean first passage time on a disc oran ellipse. The expressions we derive for perturbed discs and perturbed ellipses are testedusing numerical solutions of the governing partial differential equation. We show thatthe perturbation solutions rapidly converge to the numerical solution with just a smallnumber of terms that are straightforward to evaluate using MATLAB code supplied onGitHub. Finally, we show how to estimate the exit time in naturally–occurring shapesby representing the boundaries of Tasmania and Taiwan as a perturbed disc and ellipse,respectively, and then evaluating the exit time on these shapes.There are many avenues available to extend the work presented here. Here we considerthe most fundamental transport mechanism: unbiased diffusion, but it is also possible toconsider generalisations of Equation (1) such as drift–diffusion, diffusion–decay [32, 33] ormore complicated discrete mechanisms including L´evy flights [34,35]. A further extensionis to consider calculating both the mean first passage time and higher moments of exittime [36]. All problems in this work consider exit times by specifying absorbing bound-ary conditions in the random walk model, which correspond to homogeneous Dirichletconditions in the partial differential equation model. These can be extended to mixedboundary conditions where some parts of ∂ Ω are absorbing, and other parts of ∂ Ω are19eflecting on the original, unperturbed boundary. It would then be very interesting toconsider perturbations with such mixed boundary conditions. A more substantial exten-sion of this work would be to consider the solution of Equation (1) on more complicatedshapes that are not small, smooth perturbations of a disc or an ellipse. If, for example,we consider the case where Ω corresponds to a larger circle whose boundary just touches asmaller circle, we are unable to directly apply the techniques developed in this work, andso a different approach is required to construct approximate solutions of Equation (1). Inaddition to the more mathematical extensions described here, it is also possible to con-sider extensions of the present work that are more computational. For example, furtherconsideration could be given to the way in which the perturbation solutions presentedin this work are evaluated. In this work we evaluate the perturbation solutions usingsymbolic tools in MATLAB since it is convenient for us to provide a single software toperform stochastic random walk simulations, finite volume numerical calculations andto evaluate the perturbation solution in a single programming language. We note, how-ever, that working with a different symbolic language could be more efficient, especiallyif additional terms in the perturbation solution are to be evaluated.
Acknowledgements:
This work is supported by the Australian Research Council (DP200100177)and Queensland University of Technology for providing a summer research scholarship toDJV. We thank two referees for helpful suggestions.20 ppendix A: Numerical methods
Stochastic simulations
We simulate particle lifetime distributions using continuous space, discrete time randomwalk stochastic simulations. Time is discretized with constant time steps of duration τ >
0. In each time step, a particle, at location x ( t ), attempts to step a distance δ > x ( t + τ ) = x ( t ) + δ (cos θ, sin θ ) with probability P ∈ [0 , θ is sampled from auniform distribution, θ ∼ U [0 , π ]. This discrete process corresponds to a random walkwith diffusivity D = P δ / (4 τ ) [3]. Simulations continue until the particle steps acrossthe boundary of the domain, and the exit time is recorded. All simulations in this workuse τ = 1, P = 1 and δ = 1 × − , giving D = 2 . × − .To estimate the mean exit time we consider N identically–prepared realisations ofthe discrete process with starting position x (0) = ( x, y ). This gives us N estimates ofthe exit time from which we calculate the mean, T sim ( x, y ) = (1 /N ) N (cid:88) i =1 t i , where t i is theexit time from the i th identically prepared stochastic realisation. In practice we use thestochastic model to estimate T sim ( x, y ) with N = 1 × simulations, and these estimatesare obtained at a number of spatial points in Ω. We consider an unstructured triangularmeshing of Ω, and we estimate T sim ( x, y ) at each node. This means that for a meshingwith M nodes, we perform a total of M × N stochastic simulations. For example, theresults in Figure 1(a) are generated with N = 1000 identically–prepared realisations ateach of the M = 632 nodes, giving a total of 632000 stochastic simulations. The resultingpoint estimates of T sim ( x, y ) at each node are interpolated across Ω using the interp optionin MATLAB’s shading function [37] to provide a continuous estimate of the distributionof the mean exit time.Results in Figure 8 provide a visual comparison of the impact of varying N . The so-lution in Figure 8(a) is the exact solution of Equation (1) on the unit disc, Equation (4).21ata in Figure 8(b)–(f) show estimates of the mean exit time on the unit disc generatedusing N = 60 , , ,
500 and 1000 particles per node in the finite volume mesh. Com-paring data in Figure 8(b)–(f) we see that the fluctuations in the estimates of T sim ( x, y )appear to decrease, as expected, as N increases. Additional results in Figure 9 com-pares the exact solution and simulation-based estimates as a function of N in terms ofEquation (37). Again, we see that e ( x, y ) approaches zero as N increases. These resultsjustify our choice of N = 1000 in the main document since we roughly have e ( x, y ) < N = 1000. Of course the algorithms available on GitHub can be used to generate T sim ( x, y ) for larger N , but this comes with the drawback that increasing N leads tolonger simulation times.Figure 8: Stochastic simulations as N increases. (a) Exact solution for the meanexit time on the unit disc, Equation (4). (b)–(f) averaged simulation data for the meanexit time generated using N = 60 , , ,
500 and 1000 particles released per node inthe finite volume mesh. 22igure 9:
Comparison of stochastic simulations and exact solution as N in-creases. (a)–(f) Comparison of the exact solution for the mean exit time on the unitdist with simulation data for N = 60 , , , ,
750 and 1000, respectively. All resultsare presented in terms of e ( x, y ), given by Equation (37). Finite volume calculations
We solve Equation (1) numerically using a finite volume approximation [38] to discretizethe governing equation over an unstructured triangular meshing of Ω. To perform thesecalculations we use mesh generation software, GMSH [39]. The finite volume method isimplemented using a vertex centered strategy with nodes located at the vertices in themesh. Control volumes are constructed around each node by connecting the centroidof each triangular element to the midpoint of its edges [40]. Linear finite element shapefunctions are used to approximate gradients in each element. Assembling the finite volumeequations yields a sparse linear system that can be stored and solved efficiently. For eachnumerical solution reported in this work we report the number of nodes and elementsin the finite volume mesh, and in each case use a prescribed mesh element size of 0 . ppendix B: Truncation effects Results in Figures 3–4 compare estimates of mean exit time, T , for a perturbed disc andsphere, respectively. In these figures we compare T from repeated stochastic simulations,a truncated perturbation solution and a fine-mesh finite volume solution of the governingboundary value problem. In these comparisons we choose a particular truncation ofthe perturbation solution to ensure that the perturbation solution and the finite volumesolutions compare reasonably well. Here, we explore how the choice of truncation impactsthe accuracy of the perturbation solutions. The comparison between the perturbationand finite volume solutions in Figures 3–4 correspond to O ( ε ) perturbation solutionswith 25 terms retained in the truncated infinite sum. The visual comparison between theperturbation and finite volume solutions in Figures 3–4 indicates that the perturbationsolution is very accurate.Figure 10: Perturbation solution on a perturbed disc.
Comparison of perturbationsolution to the finite volume approximation of the exit time on the perturbed disc, using:(a) O ( ε ) perturbation solution with 1 term in the truncated infinite sum; (b) O ( ε ) pertur-bation solution with 10 terms in the truncated infinite sum; and, (c) O ( ε ) perturbationsolution with 10 terms in the truncated infinite sum. All results are presented in termsof e ( x, y ), given by Equation (37). 24igure 11: Perturbation solution on a perturbed ellipse.
Comparison of perturba-tion solution to the finite volume approximation of the exit time on the perturbed ellipse,using: (a) O ( ε ) perturbation solution with 1 term in the truncated infinite sum; (b) O ( ε )perturbation solution with 10 terms in the truncated infinite sum; and, (c) O ( ε ) pertur-bation solution with 10 terms in the truncated infinite sum. All results are presented interms of e ( x, y ), given by Equation (37).We now quantify this comparison using Equation (37). Additional results in Figure 10show plots of e ( x, y ) for the same problem examined in Figure 3, except that here we usedifferent truncations in the perturbation solution. Comparing results in Figure 10(a)–(c)indicate that the perturbation solution rapidly approaches the finite volume solution usingonly a modest number of terms. Note that the perturbation solution in Figure 3 is evenmore accurate then those in Figure 8. A similar comparison in Figure 11 confirms thatthe perturbation solution for the ellipse also rapidly approaches the numerical solutionas with only a small number of terms. Results in Figure 11 correspond to the problempreviously explored in Figure 4. 25 ppendix C: Image processing To apply our analysis to the boundary of Tasmania we use MATLAB to produce an arrayrepresentation of the boundary, and we smooth some of the boundary by refining certainjagged edges and the removal of some small peninsulas. After smoothing, we use the Sobeledge detection method in MATLAB’s cscvn [30] function and the imclose [41] functionto detect and refine the boundaries. Boundary points on the detected edges are obtainedwith bwboundaries [42]. Given a relatively dense set of points along the boundary, weretain each 30th point to give the boundary in Figure 6. We compute the mean of the x and y coordinates and shift the data so that the centre of the region is at the origin.We then scale the data so that it is comparable to a unit disc [31]. After numerical andperturbation solutions are obtained on the region contained in this boundary we rotatethe resulting solutions to match the shape of the original region.To apply our analysis to the boundary of Taiwan we start by manually smoothingsome jagged portions of the boundary, and then use the Sobel method with imclose and bwboundaries [41, 42] to represent the boundaries in terms of a dense set of points. Wework with every 40th point, and shift the region so that the centroid is at the origin,followed by a counterclockwise rotation of π/ x –axis, allowing us to apply the results in Section 2.4. The data is then scaledso that it is comparable to an ellipse with semi-major axis 2 and semi-minor axis 1 [31].We now apply the procedure described in Section 2.5 [31] to identify the best–fittingellipse, and we then use the least–squares procedure described in Section 2.5 to calculatetrigonometric polynomial representations of g ( θ ) and h ( θ ).26 ppendix D: Comparing Tasmania and pseudo-Tasmania To quantitatively examine the implications of smoothing the boundaries of the natu-ral coastlines we show additional results in Figure 12 where we compare exit times onTasmania and pseudo-Tasmania. Results in Figure 12(a)-(b) show the exit time on Tas-mania using repeated stochastic simulations and the finite volume numerical solutionof Equations (1)–(2), respectively. Figure 12(c) repeats the perturbation solution onpseudo-Tasmania shown previously in Figure 6(b). Visual inspection of the solutions onTasmania and pseudo-Tasmania in Figure 12 indicate that the solutions compare verywell in the interior of the domain, with the key difference being along the coastlines,as one might anticipate. To provide a more quantitative comparison we consider theinland location of Cradle Mountain, shown in Figure 12 as a purple disc. Our resultsgive T sim = 9427 using repeated stochastic simulations on Tasmania, whereas we obtain T p = 9673 on pseudo-Tasmania, giving a discrepancy of just 2 . Comparing Tasmania and pseudo-Tasmania. (a) Averaged data fromrepeated stochastic simulations on Tasmania. (b) Numerical solution of Equations (1)–(2) on Tasmania. (c) Exit time on pseudo-Tasmania using the perturbation solution fromFigure 6(b). On each map we show the approximate location of Cradle mountain (purpledisc). The triangular mesh used to construct the solution in (a) and (b) has 3356 nodesand 6369 triangular elements. For (a), 1000 random walks starting from each note weregenerated. 27 eferences [1] Redner S (2001) A Guide to First Passage Processes. Cambridge University Press.[2] Krapivsky PL, Redner S, Ben-Naim E (2010) A Kinetic View of Statistical Physics.Cambridge University Press.[3] Hughes BD (1995) Random Walks in Random Environments. Oxford UniversityPress.[4] Bear J (1972) Dynamics of Fluids in Porous Media. Elsevier.[5] Bird RB, Stewart WR, Lightfoot EN (2002) Transport Phenomena. John Wiley andSons.[6] Murray JD (2002) Mathematical Biology I: An Introduction. Third edition. Springer,New York.[7] Kot M (2003) Elements of Mathematical Ecology. Cambridge University Press, Cam-bridge.[8] Edelstein-Keshet L (2005) Mathematical Models in Biology. SIAM, Philadelphia.[9] L¨otstedt P, Meinecke L (2015) Simulation of stochastic diffusion via first exit times.Journal of Computational Physics. 300: 862–886.[10] Meinecke L, L¨otstedt P (2016) Stochastic diffusion processes on Cartesian meshes.Journal of Computational and Applied Mathematics. 294: 1–11.[11] Dy DLL, Esguerra JP (2008) First-passage-time distribution for diffusion through aplanar wedge. Physical Review E. 78: 062101.[12] Chupeau M, B´enichou O, Majumda SN (2015) Survival probability of a Brownianmotion in a planar wedge of arbitrary angle. Physical Review E. 91: 032106.2813] Vaccario G, Antoine C, Talbot J (2015) First–passage times in d –dimensional het-erogeneous media. Physical Review Letters. 115: 240601.[14] Rupprecht J-F, B´enichou O, Grebenkov DS, Voituriez R (2015) Exit time distri-bution in spherically symmetric two–dimensional domains. Journal of StatisticalPhysics. 158: 192–230.[15] Carr EJ, Ryan JM, Simpson MJ (2020) Diffusion in heterogeneous discs and spheres:new closed-form expressions for exit times and homogenization formulae. Journal ofChemical Physics. 153: 074115.[16] Simpson MJ, Baker RE (2015) Exact calculations of survival probability for diffusionon growing lines, disks and spheres: the role of dimension. Journal of ChemicalPhysics. 143: 094109.[17] Simpson MJ, Sharp JA, Baker RE (2015) Survival probability for a diffusive processon a growing domain. Physical Review E. 91: 042701.[18] Kurella V, Tzou JC, Coombs D, Ward MJ (2014) Asymptotic analysis of first passagetime problems inspired by ecology. Bulletin of Mathematical Biology. 77: 83-125.[19] Lindsay AE, Kolokolnikov T, Tzou JC (2015) Narrow escape problem with a mixedtrap and the effect of orientation. Physical Review E. 91: 032111.[20] Grebenkov DS, Metzler R, Oshanin G (2019) Full distribution of first exit times inthe narrow escape problem. New Journal of Physics. 21: 122001.[21] Lindsay AE, Tzou JC, Kolokolnikov T (2017) Optimization of first passage times bymultiple cooperating mobile traps. Multiscale Modeling and Simulation. 15: 915–947.[22] Grebenkov DS, Skvortsov AT (2020) Mean first–passage time to a small absorbingtarget in an elongated planar domain. New Journal of Physics. 22: 113024.2923] Nyberg M, Ambj¨ernsson T, Lizana L (2016) A simple method to calculate first-passage time densities with arbitrary initial conditions. New Journal of Physics. 18:063019.[24] McCue SW, King JR (2011) Contracting bubbles in Hele–Shaw cells with a power–law fluid. Nonlinearity. 24: 613–641.[25] McCollum WM (2014) Laplace’s equation on perturbed domains. Colorado Schoolof Mines. Master of Science Thesis.[26] Farlow SJ (1982) Partial Differential Equations for Scientists and Engineers. DoverBooks, New York.[27] Jackson D (1944) The harmonic boundary value problem for an ellipse or an ellipsoid.The American Mathematical Monthly. 51: 555–563.[28] Blank map of Tasmania. Retrieved February 2021 Tasmania.[29] Blank map of Taiwan. Retrieved February 2021 Taiwan.[30] Mathworks cscvn . Retrieved February 2021 cscvn.[31] Szpak ZL, Chojnacki W, van den Hengel A (2015) Guaranteed ellipse fitting witha confidence region and an uncertainty measure for centre, axes, and orientation.Journal of Mathematical Imaging and Vision. 52: 173–199.[32] Ellery AJ, Simpson MJ, McCue SW, Baker RE (2012). Critical time scales foradvection-diffusion-reaction processes. Physical Review E. 85: 041135.[33] Ellery AJ, Simpson MJ, McCue SW, Baker RE (2012) Moments of action provideinsight into critical times for advection-diffusion-reaction processes. Physical ReviewE. 86: 031136. 3034] Wardak A (2020) First passage leapovers of L´evy flights and the proper formulationof absorbing boundary conditions. Journal of Physics A: Mathematical and Theo-retical. 53: 375001.[35] Padash A, Chechkin AV, Dybiec B, Magdziarz M, Shokri B, Metzler R (2020) Firstpassage time moments of asymmetric L´evy flights. Journal of Physics A: Mathemat-ical and Theoretical. 53: 275002.[36] Carr EJ, Simpson MJ (2018) Rapid calculation of maximum particle lifetimes fordiffusing particles in complex geometries. Journal of Chemical Physics. 148: 094113.[37] Mathworks shading . Retrieved February 2021 shading.[38] Eymard R, Gallou¨et T, Herbin R (2000) Finite Volume Methods, Handbook ofNumerical Analysis, Volume 7. North-Holland, Amsterdam.[39] Geuzaine C, Remacle F-F (2009) Gmsh: A 3–D finite element mesh generatorwith built–in pre– and post–processing facilities. International Journal for NumericalMethods in Engineering. 79: 1309–1331.[40] Carr EJ, Perr´e P, Turner IW (2016) The extended distributed microstructure modelfor gradient–driven transport: A two–scale model for bypassing effective parameters.Journal of Computational Physics. 327: 810–829.[41] Mathworks imclose . Retrieved February 2021 imclose.[42] Mathworks bwboundariesbwboundaries