Mean exit time for diffusion on irregular domains

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,

$$\begin{equation}\frac{D}{r}\frac{\mathrm{d}}{\mathrm{d}r}\left[r\frac{\mathrm{d}T\left(r\right)}{\mathrm{d}r}\right]=-1,\quad 0{< }r{< }R.\end{equation} \tag{ 3 }$$

The solution of equation (3) with appropriate boundary conditions, T(R) = 0 and dT(0)/dr = 0, is given by

$$\begin{equation}T\left(r\right)=\frac{{R}^{2}-{r}^{2}}{4D}.\end{equation} \tag{ 4 }$$

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.

Zoom In Zoom Out Reset image size

For the special case where ∂Ω is an ellipse centered at the origin it is convenient to work in Cartesian coordinates. The ellipse, given by

$$\begin{equation}\frac{{x}^{2}}{{a}^{2}}+\frac{{y}^{2}}{{b}^{2}}=1,\end{equation} \tag{ 5 }$$

has width 2a > 0 and height 2b > 0. The solution of equation (1) on this domain is given by [24]

$$\begin{equation}T\left(x,y\right)=\frac{{a}^{2}{b}^{2}}{2D\left({a}^{2}+{b}^{2}\right)}\left[1-\frac{{x}^{2}}{{a}^{2}}-\frac{{y}^{2}}{{b}^{2}}\right].\end{equation} \tag{ 6 }$$

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.

Zoom In Zoom Out Reset image size

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=\mathcal{R}\left(\theta \right)$, subject to the condition that $\mathcal{R}\left(\theta \right)$ 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

$$\begin{equation}\mathcal{R}\left(\theta \right)=R\left(1+\varepsilon g\left(\theta \right)\right),\end{equation} \tag{ 7 }$$

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 $\mathcal{O}\left(1\right)$ smooth periodic function with period 2π. We assume that the solution can be written as

$$\begin{equation}T\left(r,\theta \right)={T}_{0}\left(r,\theta \right)+\varepsilon {T}_{1}\left(r,\theta \right)+{\varepsilon }^{2}{T}_{2}\left(r,\theta \right)+\cdots +{\varepsilon }^{n}{T}_{n}\left(r,\theta \right)+\mathcal{O}\left({\varepsilon }^{n+1}\right).\end{equation} \tag{ 8 }$$

When the $\mathcal{O}\left({{\epsilon}}^{n+1}\right)$ term is neglected in equation (8) the solution is referred to as the $\mathcal{O}\left({{\epsilon}}^{n}\right)$ perturbation solution. To proceed we evaluate T(r, θ) on $\mathcal{R}\left(\theta \right)$ and expand about ɛ = 0 to give,

$$\begin{align}\hfill 0& =T\left(R+\varepsilon Rg\left(\theta \right),\theta \right),\hfill \\ \hfill 0& =T\left(R,\theta \right)+\varepsilon Rg\left(\theta \right){\left.\frac{\partial T}{\partial r}\right\vert }_{r=R}+\cdots +\frac{{\left(\varepsilon Rg\left(\theta \right)\right)}^{n}}{n!}{\left.\frac{{\partial }^{n}T}{\partial {r}^{n}}\right\vert }_{r=R}+\mathcal{O}\left({\varepsilon }^{n+1}\right).\hfill \end{align} \tag{ 9 }$$

Substituting equation (8) into equation (9) and equating powers of ɛ leads to,

$$\begin{equation}\mathcal{O}\left(1\right):\quad {T}_{0}\left(R,\theta \right)=0,\end{equation} \tag{ 10 }$$

$$\begin{equation}\mathcal{O}\left({\varepsilon }^{\ell }\right):\quad {T}_{\ell }\left(R,\theta \right)+\sum\limits _{k=1}^{\ell }\frac{{\left(Rg\left(\theta \right)\right)}^{k}}{k!}{\left.\frac{{\partial }^{k}{T}_{\ell -k}}{\partial {r}^{k}}\right\vert }_{r=R}=0,\end{equation} \tag{ 11 }$$

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

$$\begin{equation}\mathcal{O}\left(1\right):\quad D{\nabla }^{2}{T}_{0}=-1,\qquad {T}_{0}\left(R,\theta \right)=0,\end{equation} \tag{ 12 }$$

$$\begin{equation}\mathcal{O}\left({\varepsilon }^{\ell }\right):\quad {\nabla }^{2}{T}_{\ell }=0,\qquad {T}_{\ell }\left(R,\theta \right)=-\sum\limits _{k=1}^{\ell }\frac{{\left(Rg\left(\theta \right)\right)}^{k}}{k!}{\left.\frac{{\partial }^{k}{T}_{\ell -k}}{\partial {r}^{k}}\right\vert }_{r=R},\end{equation} \tag{ 13 }$$

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

$$\begin{equation}{T}_{\ell }\left(r,\theta \right)=\frac{{a}_{0}}{2}+\sum\limits _{k=1}^{\infty }\left[{a}_{k}{r}^{k}\enspace \mathrm{cos}\left(k\theta \right)+{b}_{k}{r}^{k}\enspace \mathrm{sin}\left(k\theta \right)\right],\end{equation} \tag{ 14 }$$

where

$$\begin{equation}{a}_{k}=\frac{1}{\pi {R}^{k}}{\int }_{0}^{2\pi }{T}_{\ell }\left(R,\theta \right)\mathrm{cos}\left(k\theta \right)\enspace \mathrm{d}\theta ,\qquad {b}_{k}=\frac{1}{\pi {R}^{k}}{\int }_{0}^{2\pi }{T}_{\ell }\left(R,\theta \right)\mathrm{sin}\left(k\theta \right)\enspace \mathrm{d}\theta ,\end{equation} \tag{ 15 }$$

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.

Zoom In Zoom Out Reset image size

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.

We proceed by deriving expressions for the exit time on a perturbed ellipse by considering ∂Ω as

$$\begin{equation}x=a\left(1+\varepsilon g\left(\theta \right)\right)\mathrm{cos}\enspace \theta ,\end{equation} \tag{ 16 }$$

$$\begin{equation}y=b\left(1+\varepsilon h\left(\theta \right)\right)\mathrm{sin}\enspace \theta ,\end{equation} \tag{ 17 }$$

where 2a > 0 and 2b > 0 are the width and height of the unperturbed ellipse, respectively, with a > b and g(θ) and h(θ) are $\mathcal{O}\left(1\right)$ smooth periodic functions with period 2π. Working in Cartesian coordinates, we suppose the solution takes the form

$$\begin{equation}T\left(x,y\right)={T}_{0}\left(x,y\right)+{\epsilon}{T}_{1}\left(x,y\right)+{{\epsilon}}^{2}{T}_{2}\left(x,y\right)+\cdots +{{\epsilon}}^{n}{T}_{n}\left(x,y\right)+\mathcal{O}\left({{\epsilon}}^{n+1}\right).\end{equation} \tag{ 18 }$$

To proceed, we impose equation (2) at the boundary specified in equations (16) and (17) and expand about ${\epsilon}=0$,

$$\begin{align}\hfill 0& =T\left(a\enspace \mathrm{cos}\theta +a{\epsilon}g\left(\theta \right)\mathrm{cos}\enspace \theta ,b\enspace \mathrm{sin}\theta +b{\epsilon}h\left(\theta \right)\mathrm{sin}\enspace \theta \right)\hfill \\ \hfill 0& =T\left(a\enspace \mathrm{cos}\enspace \theta ,b\enspace \mathrm{sin}\enspace \theta \right)+\sum\limits _{k=1}^{n}\frac{1}{k!}\sum\limits _{i=1}^{k}\left(\genfrac{}{}{0.0pt}{}{k}{i}\right)\frac{{\partial }^{k}T}{\partial {x}^{i}\partial {y}^{k-i}}{\left(a{\epsilon}g\left(\theta \right)\mathrm{cos}\enspace \theta \right)}^{i}{\left(b{\epsilon}h\left(\theta \right)\mathrm{sin}\enspace \theta \right)}^{k-i}+\mathcal{O}\left({{\epsilon}}^{n+1}\right),\hfill \end{align} \tag{ 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 partial derivatives like this on the boundary of the unperturbed ellipse. For brevity it is convenient to define

$$\begin{align}\hfill {H}_{\ell }\left(\theta \right)& =-\frac{\partial {T}_{\ell -1}}{\partial y}bh\left(\theta \right)\mathrm{sin}\enspace \theta -\frac{\partial {T}_{\ell -1}}{\partial x}ag\left(\theta \right)\mathrm{cos}\enspace \theta -\frac{1}{2}\frac{{\partial }^{2}{T}_{\ell -2}}{\partial {y}^{2}}{\left(bh\left(\theta \right)\mathrm{sin}\enspace \theta \right)}^{2}\hfill \\ \hfill & \quad -\frac{{\partial }^{2}{T}_{\ell -2}}{\partial x\partial y}\left(ag\left(\theta \right)\mathrm{cos}\enspace \theta \right)\left(bh\left(\theta \right)\mathrm{sin}\enspace \theta \right)-\cdots -\frac{1}{\ell !}\frac{{\partial }^{\ell }{T}_{0}}{\partial {x}^{\ell }}{\left(ag\left(\theta \right)\mathrm{cos}\enspace \theta \right)}^{\ell }\hfill \\ \hfill & =-\sum\limits _{k=1}^{\ell }\sum\limits _{i=0}^{k}\frac{1}{k!}\left(\genfrac{}{}{0.0pt}{}{k}{i}\right)\frac{{\partial }^{k}{T}_{\ell -k}}{\partial {x}^{i}\partial {y}^{k-i}}{\left(ag\left(\theta \right)\mathrm{cos}\enspace \theta \right)}^{i}{\left(bh\left(\theta \right)\mathrm{sin}\enspace \theta \right)}^{k-i}.\hfill \end{align} \tag{ 20 }$$

Substituting equation (18) into equation (19), and equating powers of gives

$$\begin{equation}\mathcal{O}\left(1\right):\quad {T}_{0}\left(a\enspace \mathrm{cos}\enspace \theta ,b\enspace \mathrm{sin}\enspace \theta \right)=0,\end{equation} \tag{ 21 }$$

$$\begin{equation}\mathcal{O}\left({{\epsilon}}^{\ell }\right):\quad {T}_{\ell }\left(a\enspace \mathrm{cos}\enspace \theta ,b\enspace \mathrm{sin}\enspace \theta \right)={H}_{\ell }\left(\theta \right),\end{equation} \tag{ 22 }$$

for ℓ = 1, ..., n. Substituting equation (18) into equation (1) leads to the following family of boundary value problems

$$\begin{equation}\mathcal{O}\left(1\right):\quad D{\nabla }^{2}{T}_{0}=-1,\qquad {T}_{0}\left(a\enspace \mathrm{cos}\enspace \theta ,b\enspace \mathrm{sin}\enspace \theta \right)=0,\end{equation} \tag{ 23 }$$

$$\begin{equation}\mathcal{O}\left({{\epsilon}}^{\ell }\right):\quad {\nabla }^{2}{T}_{\ell }=0,\qquad {T}_{\ell }\left(a\enspace \mathrm{cos}\enspace \theta ,b\enspace \mathrm{cos}\enspace \theta \right)={H}_{\ell }\left(\theta \right),\end{equation} \tag{ 24 }$$

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 $\mathcal{O}\left({{\epsilon}}^{\ell }\right)$ boundary value problem given in equation (24) involves the boundary data H_ℓ (θ), which depends on T_ℓ−1, T_ℓ−2, ..., T₀ 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,

$$\begin{equation}{H}_{\ell }\left(\theta \right)=\frac{{a}_{0}}{2}+\sum\limits _{k=1}^{\infty }\left({a}_{k}\enspace \mathrm{cos}\left(k\theta \right)+{b}_{k}\enspace \mathrm{sin}\left(k\theta \right)\right),\quad 0{\leqslant}\theta {< }2\pi ,\quad \mathrm{w}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e},\quad \end{equation} \tag{ 25 }$$

$$\begin{equation}{a}_{k}=\frac{1}{\pi }{\int }_{0}^{2\pi }{H}_{\ell }\left(\theta \right)\mathrm{cos}\left(k\theta \right)\enspace \mathrm{d}\theta ,\qquad {b}_{k}=\frac{1}{\pi }{\int }_{0}^{2\pi }{H}_{\ell }\left(\theta \right)\mathrm{sin}\left(k\theta \right)\enspace \mathrm{d}\theta .\end{equation} \tag{ 26 }$$

We then compute ${u}_{\ell }\left(x,y\right)=\mathfrak{R}\left({\left(x+iy\right)}^{\ell }\right)$ and ${v}_{\ell }\left(x,y\right)=\mathfrak{I}\left({\left(x+iy\right)}^{\ell }\right)$ for ℓ = 1, 2, ..., and also define u₀(x, y) = 1 and v₀(x, y) = 0, given explicitly as

$$\begin{equation}{u}_{\ell }\left(x,y\right)=\sum\limits _{j=0}^{\lfloor \frac{\ell }{2}\rfloor }\left(\genfrac{}{}{0.0pt}{}{\ell }{2j}\right){x}^{\ell -2j}{\left(-1\right)}^{j}{y}^{2j},\end{equation} \tag{ 27 }$$

$$\begin{equation}{v}_{\ell }\left(x,y\right)=\sum\limits _{j=0}^{\lfloor \frac{\ell +1}{2}\rfloor }\left(\genfrac{}{}{0.0pt}{}{\ell }{2j+1}\right){x}^{\ell -2j-1}{\left(-1\right)}^{j}{y}^{2j+1}.\end{equation} \tag{ 28 }$$

The next step is to compute $2{\mathcal{T}}_{k}\left(\mathrm{cos}\enspace \theta \right)$ for k = 1, 2, ..., where ${\mathcal{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

$$\begin{equation}{U}_{k}\left(x,y\right)=\frac{1}{{\left(a+b\right)}^{k}+{\left(a-b\right)}^{k}}\sum\limits _{r=0}^{\frac{k}{2}}{C}_{2r}{\left({a}^{2}-{b}^{2}\right)}^{\frac{k}{2}-r}{u}_{2r}\left(x,y\right),\end{equation} \tag{ 29 }$$

$$\begin{equation}{V}_{k}\left(x,y\right)=\frac{1}{{\left(a+b\right)}^{k}-{\left(a-b\right)}^{k}}\sum\limits _{r=0}^{\frac{k}{2}}{C}_{2r}{\left({a}^{2}-{b}^{2}\right)}^{\frac{k}{2}-r}{v}_{2r}\left(x,y\right),\end{equation} \tag{ 30 }$$

and for odd k ⩾ 1,

$$\begin{equation}{U}_{k}\left(x,y\right)=\frac{1}{{\left(a+b\right)}^{k}+{\left(a-b\right)}^{k}}\sum\limits _{r=0}^{\frac{k-1}{2}}{C}_{2r+1}{\left({a}^{2}-{b}^{2}\right)}^{\frac{k-1}{2}-r}{u}_{2r+1}\left(x,y\right),\end{equation} \tag{ 31 }$$

$$\begin{equation}{V}_{k}\left(x,y\right)=\frac{1}{{\left(a+b\right)}^{k}-{\left(a-b\right)}^{k}}\sum\limits _{r=0}^{\frac{k-1}{2}}{C}_{2r+1}{\left({a}^{2}-{b}^{2}\right)}^{\frac{k-1}{2}-r}{v}_{2r+1}\left(x,y\right).\end{equation} \tag{ 32 }$$

This gives us the solution of our $\mathcal{O}\left({{\epsilon}}^{\ell }\right)$ problem,

$$\begin{equation}{T}_{\ell }\left(x,y\right)=\frac{{a}_{0}}{2}+\sum\limits _{k=1}^{\infty }\left({a}_{k}{U}_{k}\left(x,y\right)+{b}_{k}{V}_{k}\left(x,y\right)\right),\end{equation} \tag{ 33 }$$

where a₀, 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.

Zoom In Zoom Out Reset image size

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.

Zoom In Zoom Out Reset image size

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 ${\left(x-{x}_{c}\right)}^{2}+{\left(y-{y}_{c}\right)}^{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 $\mathcal{R}\left(\theta \right)=R\left(1+{\epsilon}g\left(\theta \right)\right)$, where we approximate g(θ) by

$$\begin{equation}g\left(\theta \right)={A}_{0}+\sum\limits _{n=1}^{G}\left({A}_{n}\enspace \mathrm{cos}\left(n\theta \right)+{B}_{n}\enspace \mathrm{sin}\left(n\theta \right)\right).\end{equation} \tag{ 34 }$$

If the points ${\left\{\left({x}_{i},{y}_{i}\right)\right\}}_{i=1}^{N}$ represent the given boundary, we compute the polar angle for each point θ_i . To represent our boundary in this way we require

$$\begin{equation}\frac{1}{R}\sqrt{{x}_{i}^{2}+{y}_{i}^{2}}-1={\epsilon}{A}_{0}+{\epsilon}\sum\limits _{n=1}^{G}{A}_{n}\enspace \mathrm{cos}\left(n{\theta }_{i}\right)+{\epsilon}\sum\limits _{n=1}^{G}{B}_{n}\enspace \mathrm{sin}\left(n{\theta }_{i}\right),\quad i=1,2,\dots ,N.\end{equation} \tag{ 35 }$$

We estimate the coefficients A₀, 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

$$\begin{equation}\sum\limits _{i=1}^{N}{\left(\frac{1}{R}\sqrt{{x}_{i}^{2}+{y}_{i}^{2}}-1-{\epsilon}{A}_{0}-{\epsilon}\sum\limits _{n=1}^{G}{A}_{n}\enspace \mathrm{cos}\left(n{\theta }_{i}\right)-{\epsilon}\sum\limits _{n=1}^{G}{B}_{n}\enspace \mathrm{sin}\left(n{\theta }_{i}\right)\right)}^{2}.\end{equation} \tag{ 36 }$$

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 $\mathcal{R}\left(\theta \right)$ 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 $\mathcal{O}\left({\varepsilon }^{2}\right)$ 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.

Zoom In Zoom Out Reset image size

To quantify the discrepancy we introduce a measure of the difference between the two solutions,

$$\begin{equation}e\left(x,y\right)=100\frac{\left\vert {T}_{\mathrm{n}}\left(x,y\right)-{T}_{\mathrm{p}}\left(x,y\right)\right\vert }{{\mathrm{max}}_{\left(x,y\right)\in {\Omega}}\vert {T}_{\mathrm{n}}\left(x,y\right)\vert },\end{equation} \tag{ 37 }$$

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

$$\begin{equation}g\left(\theta \right)={A}_{0}+\sum\limits _{n=1}^{G}\left({A}_{n}\enspace \mathrm{cos}\left(n\theta \right)+{B}_{n}\enspace \mathrm{sin}\left(n\theta \right)\right),\end{equation} \tag{ 38 }$$

$$\begin{equation}h\left(\theta \right)={C}_{0}+\sum\limits _{n=1}^{H}\left({C}_{n}\enspace \mathrm{cos}\left(n\theta \right)+{D}_{n}\enspace \mathrm{sin}\left(n\theta \right)\right).\end{equation} \tag{ 39 }$$

As before, the boundary is given by a set of points, ${\left\{\left({x}_{i},{y}_{i}\right)\right\}}_{i=1}^{N}$, and we compute the polar angle θ_i for each point. Representing the boundary using this approach leads to two systems of linear equations

$$\begin{equation}\frac{{x}_{i}}{a\enspace \mathrm{cos}{\theta }_{i}}-1={\epsilon}{A}_{0}+{\epsilon}\sum\limits _{n=1}^{G}{A}_{n}\enspace \mathrm{cos}\left(n{\theta }_{i}\right)+{\epsilon}\sum\limits _{n=1}^{H}{B}_{n}\enspace \mathrm{sin}\left(n{\theta }_{i}\right),\quad i=1,2,\dots ,N,\end{equation} \tag{ 40 }$$

$$\begin{equation}\frac{{y}_{i}}{b\enspace \mathrm{sin}{\theta }_{i}}-1={\epsilon}{C}_{0}+{\epsilon}\sum\limits _{n=1}^{H}{C}_{n}\enspace \mathrm{cos}\left(n{\theta }_{i}\right)+{\epsilon}\sum\limits _{n=1}^{H}{D}_{n}\enspace \mathrm{sin}\left(n{\theta }_{i}\right),\quad i=1,2,\dots ,N.\end{equation} \tag{ 41 }$$

As for Tasmania, we estimate the coefficients, A₀, A_n , B_n for n = 1, 2, ..., G in equation (40) and C₀, 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(b). All MATLAB files required to replicate the boundary extraction and fitting are available on GitHub. The solution in figure 7(b) is the $\mathcal{O}\left({\varepsilon }^{2}\right)$ perturbation solution. Visual comparison of the numerical and perturbation 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 very fine-scale features of the coastline of Taiwan that would never be accurately represented by a perturbed ellipse.

Zoom In Zoom Out Reset image size

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 $\mathcal{P}\in \left[0,1\right]$. Here, θ is sampled from a uniform distribution, $\theta \sim \mathcal{U}\left[0,2\pi \right]$. This discrete process corresponds to a random walk with diffusivity $D=\mathcal{P}{\delta }^{2}/\left(4\tau \right)$ [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, $\mathcal{P}=1$ and δ = 1 × 10⁻², giving D = 2.5 × 10⁻⁵.

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}_{\text{sim}}\left(x,y\right)=\left(1/N\right){\sum }_{i=1}^{N}{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³ 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 provide a visual comparison of the impact of varying N. The solution in figure 8(a) is the exact solution of equation (1) on the unit disc, equation (4). Data in figures 8(b)–(f) show estimates of the mean exit time on the unit disc generated using N = 60, 125, 250, 500 and 1000 particles per node in the finite volume mesh. Comparing data in figures 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 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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.