Mean exit time for diffusion on irregular domains

Many problems in physics, biology, and economics depend upon the duration of time required for a diffusing particle to cross a boundary. As such, calculations of the distribution of first passage time, and in particular the mean first passage time, is an active area of research relevant to many disciplines. Exact results for the mean first passage time for diffusion on simple geometries, such as lines, discs and spheres, are well-known. In contrast, computational methods are often used to study 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 time on irregular domains formed by perturbing the boundary of a disc or an ellipse. Classical perturbation expansion solutions are then constructed using the exact solutions available on a disc and an ellipse. We apply the perturbation solutions to compute the mean first exit time on two naturally-occurring irregular domains: a map of Tasmania, an island state of Australia, and a map of Taiwan. Comparing the perturbation solutions with numerical solutions of the elliptic boundary value problem on these irregular domains confirms that we obtain a very accurate solution with a few terms in the series only. MATLAB software to implement all calculations is available at https://github.com/ProfMJSimpson/Exit_time.


Introduction
Mathematical models describing diffusive transport phenomena are fundamental tools with broad applications in many areas of physics [1][2][3], engineering [4,5] and biology [6][7][8]. Traditional analysis of mathematical models of diffusion often focus on idealised problems with relatively simple geometries, whereas practical applications routinely involve more complicated geometries that are normally dealt with using computational approaches. While computational approaches are essential in many circumstances [9,10], exact analytical insight is preferable, where possible, because it can lead to mathematical expressions that explicitly highlight key relationships that are not always revealed computationally.
A fundamental property of diffusion is the concept of particle lifetime, which is a particular application of the more general concept of the first passage time [1][2][3]. Estimates of particle lifetime provide insight into the timescale required for a diffusing particle to reach a certain target, such as an absorbing boundary. Many results are known about the particle lifetime for diffusion in simple geometries, such as lines and discs [1,2]. Generalising these results to deal with other geometrical features, such as wedges [11,12], symmetric domains [13][14][15], growing domains [16,17], slender domains [18][19][20], small targets [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 irregular two-dimensional (2D) geometries. Our approach is to use classical results for diffusion on a disc and an ellipse, and then construct perturbation solutions for more general geometries [24,25]. The new perturbation solutions are evaluated using symbolic computation, and the results compare very well with numerical solutions of the  (1) and (2). (c) Numerical solution of equations (1) and (2). Parameters are τ = 1, P = 1, δ = 1 × 10 −2 and D = 2.5 × 10 −5 . The triangular mesh used to construct the solution in (a) and (c) has 632 nodes and 1183 triangular elements. For (a), 1000 random walks starting from each node were generated.
governing boundary value problem, and with averaged data from repeated stochastic simulations. Once the new perturbation solutions are derived and validated, we consider two case studies to show how these tools can be used to take a particular 2D region and represent it as a perturbed ellipse or disc. This allows us to use the new perturbation solutions to analyze the mean particle lifetime in irregular 2D geometries. MATLAB code provided on https://github.com/ProfMJSimpson/Exit_time is used to compute: (i) the symbolic evaluation of the perturbation solution; (ii) the numerical finite volume solution of the boundary value problem; and, (iii) the stochastic random walk algorithm.

Results and discussion
Standard arguments that relate unbiased random walk models on domains with absorbing boundary conditions can be used to show that the mean exit time is given by the solution of a linear ellipse partial differential equation [1][2][3]. In this work we seek solutions of that equation, where D > 0 is the diffusivity and T (x) is the mean exit time for a particle released at location x. Here the diffusivity is related to the random walk model, D = Pδ 2 /(4τ ), where δ > 0 is the step length, τ > 0 is the duration between steps and P ∈ [0, 1] is the probability that the particle undergoing Brownian motion attempts to undergo a step of length δ in the time duration τ [1][2][3]. The continuum description is valid in the constrained limit δ → 0, τ → 0 and δ 2 /τ held constant [1][2][3]. The stochastic random walk model is Pearson's walk in R 2 , and full details of how simulations are performed are given in appendix A. Key results are presented and discussed in the following format. In sections 2.1 and 2.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, and in section 2.5 we apply the new approximate solutions to study the exit time for diffusion in two naturally-occurring geometries.

Mean exit time on a disc
For the special case where Ω is a disc of radius R > 0 centered at the origin it is convenient to work in polar coordinates, where the mean first exit time depends upon the radial coordinate only, The solution of equation (3) with appropriate boundary conditions, T(R) = 0 and dT(0)/dr = 0, is given by Results in figures 1(a)-(c) provide a visual comparison of exit time data from repeated stochastic simulations, the exact solution and the numerical solution of equations (1) and (2) for a unit disc, R = 1. A complete description of the continuous space, discrete time stochastic algorithm and the finite volume numerical method used with an unstructured triangular mesh to solve equation (1) is given in appendix A.  (1) and (2). Parameters are τ = 1, P = 1, δ = 1 × 10 −2 and D = 2.5 × 10 −5 . The triangular mesh used to construct the solution in (a) and (c) has 1240 nodes and 2356 triangular elements. For (a), 1000 random walks starting from each node were generated.

Mean exit time on an ellipse
For the special case where ∂Ω is an ellipse centered at the origin it is convenient to work in Cartesian coordinates. The ellipse, given by has width 2a > 0 and height 2b > 0. The solution of equation (1) on this domain is given by [24] T(x, y) = Results in figures 2(a)-(c) provide a visual comparison of exit time data from repeated stochastic simulations, the exact solution and the numerical solution of equation (1) for an ellipse with a = 2 and b = 1.

Mean exit time on a perturbed disc
We begin working on irregular domains by calculating expressions for the exit time on a perturbed disc. Using plane polar coordinates (r, θ), we consider a region Ω with boundary r = R(θ), subject to the condition that R(θ) is a single-valued function of θ to ensure that any ray drawn from the origin intersects precisely one point of the boundary ∂Ω. If our region is conceived as a modest perturbation of a circular disc of radius R we can write where R > 0 is the radius of the unperturbed disc, θ ∈ [0, 2π) is the polar angle, ε 1 is a small dimensionless parameter and g(θ) is an O(1) smooth periodic function with period 2π. We assume that the solution can be written as T(r, θ) = T 0 (r, θ) + εT 1 (r, θ) + ε 2 T 2 (r, θ) + · · · + ε n T n (r, θ) + O(ε n+1 ).
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, Substituting equation (8) into equation (9) and equating powers of ε leads to, for = 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, with the boundary conditions given by equations (10) and (11). This family of boundary value problems can be summarised as for = 1, . . . , n. The solution of equation (12) is equation (4), and each higher order term, T (r, θ), = 1, . . . , n, is the solution of Laplace's equation on a disc with prescribed boundary data associated with the previous term. Each higher order term can be obtained using separation of variables, giving where for = 1, . . . , n. These solutions are straightforward to evaluate symbolically, and an MATLAB implementation of the symbolic evaluation of them is available on GitHub. Results in figures 3(a)-(c) provide a visual comparison of exit time data from repeated 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/20. These results show that the truncated perturbation solution is very accurate, despite using only the first three terms in equation (8) and just the first 25 terms to approximate the infinite sum in equation (14). Perturbation solutions with different choices of truncation, different choices of ε, and different choices of g(θ) can be evaluated using the MATLAB algorithms available on GitHub. Additional information about how the choice of truncation in equation (8) affects the accuracy of the perturbation solution is given in the appendix B.
It is worth explicitly pointing out some interesting features that arise when we solve equation (1) on a perturbed disc. In sectors where g(θ) > 0, there are points of Ω that lie beyond the circle r = R, whereas in sectors where g(θ) < 0, there are points that lie within the circle r = R but are not within Ω. This situation often arises when solving boundary value problems where the location of the boundary is perturbed [24,26]. The key point being that the domain of r in the functions T(r, θ) and T i (r, θ) for i = 1, 2, 3, . . . , is the same, r < R(1 + εg(θ)). However, the boundary value problems that define the terms in the perturbation solution, T i (r, θ) for i = 1, 2, 3, . . . , are defined on the original unperturbed domain, r < R. Accordingly, our finite volume calculations and stochastic simulations are performed directly on the perturbed domain, r < R(1 + εg(θ)), whereas the perturbation solution is calculated on the unperturbed domain r < R with the 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 construct such approximate solutions. The same situation arises when we solve equation (1) on a perturbed ellipse, as we shall now explain.

Mean exit time on a perturbed ellipse
We proceed by deriving expressions for the exit time on a perturbed ellipse by considering ∂Ω as where 2a > 0 and 2b > 0 are the width and height of the unperturbed ellipse, respectively, with a > b and g(θ) and h(θ) are O(1) smooth periodic functions with period 2π. Working in Cartesian coordinates, we suppose the solution takes the form To proceed, we impose equation (2) at the boundary specified in equations (16) and (17) and expand about = 0, 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 partial derivatives like this on the boundary of the unperturbed ellipse. For brevity it is convenient to define Substituting equation (18) into equation (19), and equating powers of gives for = 1, . . . , n. Substituting equation (18) into equation (1) leads to the following family of boundary value problems for = 1, . . . , n. The solution of equation (23) is equation (6), and each higher order term is the solution of Laplace's equation on the ellipse with prescribed boundary data associated with the previous terms. For example, the O( ) boundary value problem given in equation (24) involves the boundary data H (θ), which depends on T −1 , T −2 , . . . , T 0 as evident in equation (20). Following Jackson [27] we construct the solution in terms of harmonic polynomials by expanding the boundary data H (θ) as a Fourier series, We then compute u (x, y) = R (x + iy) and v (x, y) = I (x + iy) for = 1, 2, . . . , and also define u 0 (x, y) = 1 and v 0 (x, y) = 0, given explicitly as  (1) and (2). Parameters are τ = 1, P = 1, δ = 1 × 10 −2 and D = 2.5 × 10 −5 . The triangular mesh used to construct the solution in (a) and (c) has 1243 nodes and 2342 triangular elements. For (a), 1000 random walks starting from each node were generated.
The next step is to compute 2T k (cos θ) for k = 1, 2, . . . , where T k is the Chebyshev polynomial of the first kind of degree k, and extract the coefficients of cos r θ, C r , in this polynomial. Using these expressions we compute, for even k 2 and for odd k 1, This gives us the solution of our O( ) problem, where a 0 , a k and b k are the Fourier coefficients from the boundary data, equations (25) and (26). This solution is straightforward to evaluate symbolically, and an MATLAB implementation to evaluate it is available on GitHub. Results in figures 4(a)-(c) provide a visual comparison of exit time data from repeated stochastic simulations, the perturbation solution 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. The solutions 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 produce a 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 using the MATLAB algorithms on GitHub.

Case studies
To showcase how the perturbation solutions can be applied to naturally-occurring geometries, we now turn to the problem of taking an image of a region, approximating that region as a perturbed disc or perturbed ellipse, and then using our approach to estimate the exit time from that region. For this exercise we consider the islands of Tasmania and Taiwan. In particular we represent the boundary of Tasmania as a perturbed disc and the boundary of Taiwan as a perturbed ellipse, based on the maps in figure 5. Note that we have deliberately omitted any scale on the maps in figure 5 since we wish to focus on the role of shape rather than size. This allows us to represent the boundary of Tasmania as a perturbed unit disc, and to compare the exit time from this realistic geometry with the exit times from the more idealised shapes in figures 1 and 3. Similarly, we represent the boundary of Taiwan as a perturbed ellipse on a scale that allows us to compare the exit time distribution from this realistic geometry with the results in figures 2 and 4. We now explain how we process the images in figure 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, we use MATLAB image processing tools to describe the boundary of Tasmania as a set of points as described in the appendix C, and we fit the disc (x − x c ) 2 + (y − y c ) 2 = R 2 to those points using a spline approximation in MATLAB's cscvn [30] function. Second, we shift this disc so that it is centered at the origin, and assume that the boundary takes the form R(θ) = R(1 + g(θ)), where we approximate g(θ) by If the points {(x i , y i )} N i=1 represent the given boundary, we compute the polar angle for each point θ i . To represent our boundary in this way we require We estimate the coefficients A 0 , A n , B n for n = 1, 2, . . . , G by computing the least-squares solution of equation (35) that minimises the sum of squared residuals This least-squares solution is computed using MATLAB's backslash operator. For simplicity we set = 1/10 in this case study. Given our estimates of the coefficients in equation (35), our approximation of the boundary is shown in figure 6, where we note that the approximation amounts to neglecting fine-scale features of the Tasmanian coastline. 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 drawn from the origin intersects precisely one point of the boundary ∂Ω. This condition can only be met by neglecting the fine-scale structures of the boundary, especially at the South-East part of Tasmania where there is an obvious peninsula. Given this approximation, we apply the perturbation analysis in section 2.3 to give the results in figure 6. Figure 6(a) shows the numerical solution of equation (1) within the region enclosed by the boundary obtained by truncating equation (35) with G = 3 with boundary coordinate data obtained from the map in figure 5(a). All MATLAB files required to replicate the boundary extraction and fitting are available on GitHub. Figure 6(b) shows the O(ε 2 ) perturbation solution, where infinite sums are approximated using the first 25 terms in equation (14). Visual comparison of the numerical and perturbation solutions indicates that the perturbation solution is remarkably accurate given that the domain in figures 6(a) and (b) is quite far from a unit disc. In fact, the maximum difference between the boundary in figures 6(a) and (b) and the underlying unit disc is, approximately, 0.64, confirming that this domain is a reasonably large perturbation of a unit disc. Careful comparison of the numerical and perturbation solutions show some discrepancy, particularly near the southern and north eastern portions of the boundary.
To quantify the discrepancy we introduce a measure of the difference between the two solutions, where T n (x, y) is the numerical finite volume solution and T p (x, y) is the truncated perturbation solution, such that e(x, y) is a convenient measure of the percentage relative error as a function of position. The plot of e(x, y) in figure 6(c) confirms that the perturbation solution is remarkably accurate in the interior of the domain, with small discrepancies along some of the boundary. While the small discrepancy along some parts of the boundary are visually discernable, these differences are not overwhelming, and the basic features of the numerical solution is captured by the perturbation solution. More accurate perturbation solutions can be constructed by including more terms in the truncated perturbation solution, including more terms in the infinite sums, or both. These options may be explored using the MATLAB algorithms provided on GitHub. More details about the implications of approximating Tasmania by pseudo-Tasmania are given in appendix D.
To represent Taiwan as a perturbed ellipse we first follow image pre-processing outlined in the appendix C. The method developed by Szpak et al [31] allows us to approximate the boundary as an ellipse with a particular orientation and center. We then rotate and shift this identified ellipse so that it is centered at the origin with semi-major axis along the x-axis. To approximate the boundary of Taiwan as a perturbed ellipse, equations (16) and (17), we represent the functions g(θ) and h(θ) as (C n cos(nθ) + D n sin(nθ)) .
As before, the boundary is given by a set of points, , and we compute the polar angle θ i for each point. Representing the boundary using this approach leads to two systems of linear equations As for Tasmania, we estimate the coefficients, A 0 , A n , B n for n = 1, 2, . . . , G in equation (40) and C 0 , C n , D n for n = 1, 2, . . . , 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 equations (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 by the boundary obtained by truncating equations (38) and (39) with G = H = 9 that is based on boundary data obtained from the map in figure 5

Conclusions and outlook
In this work we consider the canonical problem of determining the mean first passage time for diffusion, which requires the solution of an elliptic partial differential equation on the domain of interest. This problem has been studied, in detail, both analytically and computationally, with many known exact solutions for relatively simple geometries, such as lines, discs and spheres [1][2][3]. The calculation of exact expressions for the mean first passage time for more complicated geometries is an active, and ongoing field of research. In this work we present new solutions for the mean first passage time for diffusion on irregular two-dimensional domains, where these solutions are obtained in terms of a perturbation of the classical results for the mean first passage time on a disc or an ellipse. The expressions we derive for perturbed discs and perturbed ellipses are tested using numerical solutions of the governing partial differential equation. We show that the perturbation solutions rapidly converge to the numerical solution with just a small number of terms that are straightforward to evaluate using MATLAB code supplied on GitHub. Finally, we show how to estimate the exit time in naturally-occurring shapes by 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 consider the most fundamental transport mechanism: unbiased diffusion, but it is also possible to consider generalisations of equation (1) such as drift-diffusion, diffusion-decay [32,33] or more complicated discrete mechanisms including Lévy flights [34,35]. A further extension is to consider calculating both the mean first passage time and higher moments of exit time [36]. All problems in this work consider exit times by specifying absorbing boundary conditions in the random walk model, which correspond to homogeneous Dirichlet conditions in the partial differential equation model. These can be extended to mixed boundary conditions where some parts of ∂Ω are absorbing, and other parts of ∂Ω are reflecting on the original, unperturbed boundary. It would then be very interesting to consider perturbations with such mixed boundary conditions. A more substantial extension of this work would be to consider the solution of equation (1) on more complicated shapes 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 a smaller circle, we are unable to directly apply the techniques developed in this work, and so a different approach is required to construct approximate solutions of equation (1). In addition to the more mathematical extensions described here, it is also possible to consider extensions of the present work that are more computational. For example, further consideration could be given to the way in which the perturbation solutions presented in this work are evaluated. In this work we evaluate the perturbation solutions using symbolic tools in MATLAB since it is convenient for us to provide a single software to perform stochastic random walk simulations, finite volume numerical calculations and to evaluate the perturbation solution in a single programming language. We note, however, that working with a different symbolic language could be more efficient, especially if additional terms in the perturbation solution are to be evaluated.

Acknowledgments
This work is supported by the Australian Research Council (DP200100177) and Queensland University of Technology for providing a summer research scholarship to DJV. We thank two referees for helpful suggestions.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// github.com/ProfMJSimpson/Exit_time.

A.1. Stochastic simulations
We simulate particle lifetime distributions using continuous space, discrete time random walk 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 δ > 0, to x(t + τ ) = x(t) + δ(cos θ, sin θ) with probability P ∈ [0, 1]. Here, θ is sampled from a uniform distribution, θ ∼ U[0, 2π]. This discrete process corresponds to a random walk with diffusivity D = Pδ 2 /(4τ ) [3]. Simulations continue until the particle steps across the boundary of the domain, and the exit time is recorded. All simulations in this work use τ = 1, P = 1 and δ = 1 × 10 −2 , giving D = 2.5 × 10 −5 .
To estimate the mean exit time we consider N identically-prepared realisations of the discrete process with starting position x(0) = (x, y). This gives us N estimates of the exit time from which we calculate the mean, T sim (x, y) = (1/N) N i=1 t i , where t i is the exit time from the ith identically prepared stochastic realisation. In practice we use the stochastic model to estimate T sim (x, y) with N = 1 × 10 3 simulations, and these estimates are obtained at a number of spatial points in Ω. We consider an unstructured triangular meshing of Ω, and we estimate T sim (x, y) at each node. This means that for a meshing with M nodes, we perform a total of M × N stochastic simulations. For example, the results in figure 1(a) are generated with N = 1000 identically-prepared realisations at each of the M = 632 nodes, giving a total of 632 000 stochastic simulations. The resulting point estimates of T sim (x, y) at each node are interpolated across Ω using the interp option in MATLAB's shading function [37] to provide a continuous estimate of the distribution of the mean exit time.
Results in figure 8 figure 9 compares the exact solution and simulation-based estimates as a function of N in terms of equation (37). Again, we see that e(x, y) approaches zero as N increases. These results justify our choice of N = 1000 in the main document since we roughly have e(x, y) < 5% for 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 to longer simulation times.

A.2. Finite volume calculations
We solve equation (1) numerically using a finite volume approximation [38] to discretize the governing equation over an unstructured triangular meshing of Ω. To perform these calculations we use mesh generation software, GMSH [39]. The finite volume method is implemented using a vertex centered strategy with nodes located at the vertices in the mesh. Control volumes are constructed around each node by connecting the centroid of each triangular element to the midpoint of its edges [40]. Linear finite element shape functions are used to approximate gradients in each element. Assembling the finite volume equations  yields a sparse linear system that can be stored and solved efficiently. For each numerical solution reported in this work we report the number of nodes and elements in the finite volume mesh, and in each case use a prescribed mesh element size of 0.08 in GMSH [39] to generate these meshes. A MATLAB implementation of the numerical algorithm is available on GitHub.  We now quantify this comparison using equation (37). Additional results in figure 10 show plots of e(x, y) for the same problem examined in figure 3, except that here we use different truncations in the perturbation solution. Comparing results in figures 10(a)-(c) indicate that the perturbation solution rapidly approaches the finite volume solution using only a modest number of terms. Note that the perturbation solution in figure 3 is even more accurate then those in figure 8. A similar comparison in figure 11 confirms that the perturbation solution for the ellipse also rapidly approaches the numerical solution as with only a small number of terms. Results in figure 11 correspond to the problem previously explored in figure 4.

Appendix C. Image processing
To apply our analysis to the boundary of Tasmania we use MATLAB to produce an array representation of the boundary, and we smooth some of the boundary by refining certain jagged edges and the removal of some small peninsulas. After smoothing, we use the Sobel edge detection method in MATLAB's edge [41] function and the imclose [42] function to detect and refine the boundaries. Boundary points on the detected edges are obtained with bwboundaries [43]. Given a relatively dense set of points along the boundary, we retain 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 center of the region is at the origin. We then scale the data so that it is comparable to a unit disc [31]. After numerical and perturbation solutions are obtained on the region contained in this boundary we rotate the resulting solutions to match the shape of the original region.
To apply our analysis to the boundary of Taiwan we start by manually smoothing some jagged portions of the boundary, and then use the Sobel method with imclose and bwboundaries [42,43] to represent the boundaries in terms of a dense set of points. We work with every 40th point, and shift the region so that the centroid is at the origin, followed by a counterclockwise rotation of π/3 so that the semi-major axis co-insides with the x-axis, allowing us to apply the results in section 2.4. The data is then scaled so 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-fitting ellipse, and we then use the least-squares procedure described in section 2.5 to calculate trigonometric polynomial representations of g(θ) and h(θ).

Appendix D. Comparing Tasmania and pseudo-Tasmania
To quantitatively examine the implications of smoothing the boundaries of the natural coastlines we show additional results in figure 12 where we compare exit times on Tasmania and pseudo-Tasmania. Results in figures 12(a) and (b) show the exit time on Tasmania using repeated stochastic simulations and the finite volume numerical solution of equations (1)-(2), respectively. Figure 12(c) repeats the perturbation solution on pseudo-Tasmania shown previously in figure 6(b). Visual inspection of the solutions on Tasmania and pseudo-Tasmania in figure 12 indicate that the solutions compare very well 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 the inland location of Cradle Mountain, shown in figure 12 as a purple disc. Our results give T sim = 9427 using repeated stochastic simulations on Tasmania, whereas we obtain T p = 9673 on pseudo-Tasmania, giving a discrepancy of just 2.6%.