An Exploration in Diffusion and the Mean First Passage
Time
Henry Jones
Advised by Dr. Ike Agbanusi
April 2023
Abstract
The subject of interest in a variety of applied and theoretical fields, the Mean
First Passage Time (MFPT) is a solution to a particular Poisson equation deeply
connected to the classical diffusion equation. Within the scope of our exploration,
the MFPT is interpreted as the average time for a diffusing particle or random walker
to arrive or ‘react’ at an absorbing boundary. This thesis aims to contextualize the
study of diffusion and random walks and examine comparisons between the discrete
and continuous. We begin with several foundational derivations before delving into
some analytical and numerical results for a particularly interesting MFPT problem.
1 Introduction
Consider a particle in a two dimensional, diffusive medium bounded within some containing do-
main. Given an initial distribution of particles and the particular geometry of the domain, one
could ask a variety of interesting questions about the effect of simple and non-simple geometries
on the time-evolution of the particle distribution. Such a scenario is easily made more complicated
(and realistic) by allowing the boundary to ‘leak’ particles, or for particles to react in some way
along certain portions of the boundary. While such questions are mathematical abstractions, there
are obvious applications to a multitude of biological and chemical processes [1,2,3].
Interestingly, problems such as these are similarly understood by random walks on domains of
interest. As will be seen in the following section, the evolution of an n-dimensional random walk
in the continuum limit is exactly described by the long-time behavior of the classical diffusion
equation. As such, the computational ease of random walks provide a powerful computational
approximation of diffusion in a variety of domains. To understand this fundamental connection
between random walks and diffusion, we turn to a derivation of the diffusion equation from a
random-walk.
1
1.1 A Derivation of the Diffusion Equation
We here show the derivation of the diffusion equation using a discrete random-walk process in
one dimension. To begin, consider the partition of the real line with intervals of length ∆xand
endpoints x0, x1, . . . , xnin increasing order. For each discrete time step ∆t, we define a function
Ui,j (ti, xj) = Pr{Xj=xj}to be the probability the random variable Xjassumes the state (in this
case location on the real line) xjfor j= 0,1, . . . , n at time tifor ti=i∆twhere i∈N.
Next, we let pand qbe the probability of increasing or decreasing Xjby ∆xfor a single
increment in time ∆t. If we suppose our random variable must change state every time iteration,
we must have p+q= 1, and thus
U(ti+1, xj) = pU (ti, xj−1) + qU (ti, xj+1).
If we assume Uis at least once differentiable with respect to tand at least twice with respect
to x, we can expand our equation to give
U(ti, xj) + Ut(ti, xj)∆t+O(∆t)2
=p"U(ti, xj)−Ux(ti, xj)∆x+1
2Uxx(ti, xj)(∆x)2+O(∆x)3#
+q"U(ti, xj) + Ux(ti, xj)∆x+1
2Uxx(ti, xj)(∆x)2+O(∆x)3#.
Since p+q= 1, we can cancel U(ti, xj) from both sides of our equation to give
Ut(ti, xj) = (q−p)∆x
∆tUx(ti, xj) + (∆x)2
2∆tUxx(ti, xj) + O(∆x)3
∆t−O(∆t)2
∆t
In the particular case of unbiased diffusion, we let p=q=1
2as diffusion to the right or left in one
dimension occurs with equal probability. Thus our equation simplifies to
Ut(ti, xj) = (∆x)2
2∆tUxx(ti, xj) + O(∆x)3
∆t−O(∆t)2
∆t.
Omitting the details of the continuum limit, we refine our time and space discretization by taking
∆t→0 and ∆x→0 in such a way as to ensure that (∆x)2
2∆t→k > 0. In this way our terms of
higher order numerators go to zero and the classical heat equation emerges:
Ut(t, x) = kUxx(t, x).
This result illustrates that the behavior of a diffusive medium can equivalently be described as
an unbiased stochastic process [4,5]. It is this connection between the long-time behavior of the
diffusion equation and random motion which allows for discrete random walks to be powerful
approximations of diffusion. Having established this deep connection, we now introduce the Mean
First Passage Time (MFPT).
1.2 First Passage Time Problems
For a particular particle in random motion in a domain Ω, the time taken for this particle to
arrive at a point of interest in Ω, say at an absorbing portion of the boundary ∂Ωa,is known as the
2
particle’s First Passage Time. Note that this first passage time only has relevance to the particles
starting position x0∈Ω. Given the stochastic nature of any such walk, of greater interest to the
applied sciences is the expected value of this random variable, i.e. the MFPT. That is, given a
particle’s particular starting position x0, the MFPT is the average time elapsed before the particle
arrives at ∂Ωa.
Depending on the particular field and physical process under investigation, a huge variety
of domains and boundary conditions can be considered. For instance, in the investigation of
chemo-receptor dynamics and membrane process, reactions have particular activation energies and
chemical species must be appropriately oriented such that an enzyme-substrate complex(or any
analogous intermediary complex) can form (we note all chemical reactions are constrained by
reactant orientation). Such considerations require more realistic assumptions.
To this end, many different models have been constructed in various fields, including the well-
known Robin boundary condition. The Robin boundary condition allows for particles to have
specific ‘reaction’ probabilities at ∂Ωawhich allows for more realistic physical modelling. In the
context of our investigation, the simplest boundary conditions are the classical Dirichlet-Neumann
boundary conditions. Application of the Dirichlet boundary condition on ∂Ωaencodes the first-
contact reaction or absorption and the Neumann condition on ∂Ω\∂Ωaencodes the inert nature
or impermeability of the remainder of the boundary. The following section and derivation of the
MFPT will involve this simpler case of Dirichlet-Neumann boundary conditions.
1.3 A Diffusive Boundary Problem
Let Ω be an interval in Ror a domain in R2or R3. Consider some initial distribution of particles
in Ω which diffuse over time until they exit some particular boundary of Ω which we denote by
∂Ωa. We will assume that initially there is a single particle at location x0∈Ω.
We endeavor to determine the MFPT over the domain, i.e. the average time a particle will take
to exit through ∂Ωa.As we are interested in the long-term diffusive behavior, our investigation is
rooted in the classical diffusion equation.
We interpret U(t, ¯x) as the probability of locating a particle at ¯xat time twhich satisfies
the classical diffusion equation. For now, we specify the boundary condition on the entire domain
boundary ∂Ω as the Dirichlet type which defines the probability of finding particles on the boundary
to be zero. This boundary condition encodes the instantaneous ‘absorption’ expected of a Dirichlet
boundary condition. Our initial distribution is the “Dirac mass” which places all probability for
time t= 0 at x0, the initial location.
Our boundary value problem is thus
Ut−k∆U= 0; t > 0x∈Ω
U(t, x) = 0; t > 0x∈∂Ω (1)
U(0, x) = δ(x−x0).
The “Dirac Mass” initial condition is a generalized function called a distribution and by defi-
3
nition satisfies
ZR
δ(x−x0)ϕ(x)dx =ϕ(x0),(2)
for all smooth functions ϕthat vanish outside some compact set.
One way to solve this boundary value problem is to separate variables and write
U(t, x) = T(t)X(x).
Our diffusion equation then becomes
T′(t)X(x)−kT (t)∆X(x) = 0,
or
T′
kT =∆X
X=−λ.
Thus we have an eigenvalue problem for both Tand X:
T′=−λkT and −∆X=λX for X∈Ω, X|∂Ω= 0.
In one dimension on the interval (0, L), this eigenvalue problem for our spatial variable Xis a
simple second-order ordinary differential equation,
−X′′(x) = λX(x) for 0 < x < L, (3)
X(0) = X(L) = 0.(4)
We know that this ordinary differential equation has a discrete set of eigenvalues λn= (nπ
L)2for
n∈Nwith corresponding eigenfunctions ϕn(x) = sin(nπx
L).With these eigenvalues, we can solve
the first order ODE for Tto get
Tn(t) = An(e−λnkt) = Ane−(nπ
L)2kt.
Taking the the linear combination of our solutions, we see that
U(t, x) = ∞
X
n=1
Ane−kλntϕn(x) (5)
=∞
X
n=1
Ane−(nπ
L)2ktsin nπx
L.(6)
All that remains is to determine our coefficients Anby our initial condition. By our knowledge
of Fourier coefficients,
An=2
LZL
0
U(0, x) sin nπx
Ldx (7)
=2
LZL
0
δ(x−x0) sin nπx
Ldx (8)
=2
Lsin nπx0
L.(9)
Thus for the one-dimensional case with absorbing endpoints on (0, L),
U(t, x) = 2
L
∞
X
n=1
e−(nπ
L)2kt sin nπx0
Lsin nπx
L.(10)
4
Having dealt with a simple example, let us now consider a more general situation. We now
assume Ω is simpply a bounded domain with a smooth or piece-wise smooth boundary.
We again look to solve the eigenvalue problem in the spatial variable (using Dirichlet boundary
conditions),
−∆X(x) = λX(x) for x∈Ω (11)
X|∂Ωa= 0.(12)
For such domains it is a somewhat advanced fact that there is a discrete set of positive eigen-
values 0 < λ1< λ2≤λ3≤λ4≤. . . ≤λn≤. . . going to infinity eigenfunctions ϕn(x). Except for
simple domains, it is difficult to obtain exact expressions for such eigenvalues and eigenfunctions.
Despite these complications, there are a few useful properties we will take for granted which
allow us to gain insight into this problem. First, the eigenfunctions are smooth in Ω; that is, they
and all their derivatives are continuous. Secondly, such eigenfunctions can be chosen to satisfy the
“normalization”
ZΩ
(ϕn(x))2dx = 1.
And finally it is known that the n-th eigenvalue satisfies
λn≈cdV ol(Ω)n2
d,
where dis the dimension and V ol(Ω) is the volume of the domain Ω, and cdis some constant
dependent only on d.
Taking these advanced properties for granted, it follows that the solution of (1) is
U(t, x) = ∞
X
n=1
Ane−kλntϕn(x),
which is analogous to the one-dimensional solution. The coefficients are given by
An=ZΩ
U(0, x)ϕn(x)dx (13)
=ZΩ
δ(x−x0)ϕn(x)dx (14)
=ϕn(x0).(15)
Our general solution is therefore
U(t, x, x0) = ∞
X
n=1
e−kλntϕn(x)ϕn(x0),(16)
where we note that now the solution depends upon the starting point x0. Having established this
diffusion solution, the reader will recall that our investigation is interested in the MFPT.
1.4 The Mean First Passage Time
First, we can let the random variable τdenote the first time a diffusing particle hits the ‘target’
or domain boundary of interest. We denote the probability of finding our particle somewhere in Ω
5
without exiting through the domain boundary ∂Ωaas Pr(τ > t). By our result in the preceding
section we have
Pr(τ > t) = ZΩ
U(t, x, x0)dx.
Then for any positive random variable R, we can say that
E[R] = Z∞
0
Pr(τ > t)dt. (17)
It then follows that if we let V(x0) = E[τ],then
V(x0) = Z∞
0ZΩ
U(t, x, x0)dxdt. (18)
Returning to our general formulation for U(t, x, x0),we see that
V(x0) = Z∞
0ZΩ
∞
X
n=1
e−kλntϕn(x)ϕn(x0)dxdt. (19)
With the appropriate assumptions of smoothness and integrability, (in part granted by the conver-
gence of the exponential and the earlier assumption of the ‘normalization’ or L2condition of our
eigenfunctions) the order of integration and summation can be interchanged to write
V(x0) = ∞
X
n=1
ϕn(x0)ZΩ
ϕn(x)dx Z∞
0
e−kλntdt. (20)
The right-most integral can be evaluated as
Z∞
0
e−kλnt=e−kλnt
−kλn
∞
0
=1
kλn
.
Therefore we have
V(x0) = ∞
X
n=1
ϕn(x0)
kλnZΩ
ϕn(x)dx =∞
X
n=1 1
kλnZΩ
ϕn(x)dxϕn(x0).(21)
The final observation, which is certainly not obvious, is the fact that this expression is a solution
to a certain Poisson PDE, namely V(y) solves
−∆yV=1
kfor y∈Ω (22)
V|∂Ωa= 0.(23)
To see this, consider seeking a solution of the general form −∆V(y) = g(y) with the same mixed
Dirichlet-Neumann boundary conditions in our direct solution. Given the assumptions made earlier
about the eigenvalues and eigenfunctions for such a problem on this domain, we can write gas a
linear combination in terms of these eigenfunctions:
g(y) = ∞
X
n=0
fnϕn(y),
where
fn=ZΩ
g(t)ϕn(t)dt.
Similarly we can write V(y) as ∞
X
n=0
Bnϕn(y)
6
for some coefficients Bnto be determined.
If we then insert this form of V(y) into the Poisson equation, we are left with
−∆V(y) = −∆∞
X
n=0
Bnϕn(y) = ∞
X
n=0
Bn(−∆ϕn) (24)
=∞
X
n=0
Bnλnϕn(y).(25)
Hence
∞
X
n=0
Bnλnϕn(y) = ∞
X
n=0 ZΩ
g(t)ϕn(t)dtϕn(y).(26)
If we let An=RΩg(t)ϕn(t)dt and equate coefficients, we see that
Bnλn=Anor Bn=An
λn
.
Comparing this result to (21) we see that we must have
1
λnZΩ
g(t)ϕn(t)dt =1
kλnZΩ
ϕn(t)dt, (27)
which means we must have g(t) = 1
k,thus completing this derivation of the Poisson equation for
the MFPT.
2 Analytical Investigations
Having laid the foundations for our analytical investigation, several relevant MFPT solutions
will be presented in one, two, and three dimensions before introducing the MFPT for the narrow
escape problem, a large motivation for this thesis.
2.1 One Dimensional Solutions
2.1.1 Dirichlet-Neumann Boundary Value Problem on the Line
We begin with our Poisson problem on the line [0, b] with a Neumann boundary condition (BC) at
x= 0 and a Dirichlet BC at x=L > 0.Thus our problem is given by
∆V=−1
k
Vx(0) = 0
and V(b)=0.
Integrating twice we find
V(x) = ZZ −1
kdxdx =−x2
2k+c1x+c2.
Using our mixed Dirichlet-Neumann BCs we see that
V(x) = −x2
2k+b2
2k.(28)
7
We note too that following the direct method detailed in Section 1.3, the same problem’s
solution is given by
2
kb
∞
X
n=0 2b
(1 + 2n)π3sin (1 + 2n)π
2cos (1 + 2n)π
2bx (29)
which can be seen graphically.
2.1.2 Robin-Neumann Boundary Value Problem on the Line
Prior to this example, we introduce the Robin boundary condition as a ‘reactive’ BC with behav-
ior between the absolute absorption and reflection of Dirichlet and Neumann boundary conditions
respectively. The Robin BC at the end our the line at x=bis
V(b) = −σV ′(b),
where σis the constant parameter which determines the behavior. Note that when σ= 0,we have
a Dirichlet BC, and as σ→ ∞,we have a Neumann BC.
Similar to the preceding result the Robin-Neumann problem is the solution of
∆V=−1
k
Vx(0) = 0
and V(b) = −σV ′(b).
Integrating twice we find
V(x) = ZZ −1
kdxdx =−x2
2k+c1x+c2.
Using our mixed Neumann-Robin boundary conditions we see that
V(x) = −x2
2k+b2
2k+σb
k.
Thus when σ > 0, we have greater MFPTs than the simpler Neumann-Dirichlet problem in the
previous section. This increase is constant with respect to the starting position x.
2.2 Two and Three Dimensional Results
In considering simple two and three dimensional domains, we have to consider the distinct
form of the Laplacian operator in polar and spherical coordinates. As will be seen in section 3,
non-radial symmetry leads to singularities and great complexity, so our analysis in this section
pertains to radially symmetric solutions.
2.2.1 A Reflective Perimeter Disk
Consider the MFPT on the unit disk with an absorbing inner circle of radius b < 1.The
boundary conditions for such a problem are Vr(1) = 0 and V(b) = 0 (also V(1), Vr(b)<∞). Using
the radially symmetric, polar Laplacian, we can write
V′′ +1
rV′=−1
k.
8
Letting U=V′, we have
U′+1
rU=−1
k.
Thus we have integrating factor eR1
rdr =eln r=rand thus
d
dr rU=−r
k.
Therefore
rU =−Zr
k=−r2
2k+c1.
Reverting back to Vwe have
V′=−r
2k+c1
r.
Integrating once more we have
V=−r2
4k+c1ln(r) + c2.
Imposing our boundary conditions we have c1=1
2kand c2=b2
4k−ln(b)
2kThus we have
V(r) = −r2
4k+ln(r)
2k+b2
4k−ln(b)
2k.
We can easily consider the case of our inner disk with a Robin boundary condition, which leaves
gives
V(r) = −r2
4k+ln(r)
2k−σ −b
2k+1
2kb !+b2
4k−ln(b)
2k.
2.2.2 A Reflective Outer Sphere with Reactive Inner Sphere
By a very similar process using the spherically symmetric Laplacian ∆V=Vrr +2
rVr, the solution
with an absorbing inner sphere of radius b < 1 is given by
V(r) = −r2
6k−1
3kr +b2
6k+1
3kb ,
which holds on the interval r∈[b, 1].
In the case of a reactive (Robin condition) center, we can see that
V(r) = −r2
6k−1
3kr +b2
6k+1
3kb −σ −b
3k+1
3kb2!.
3 The Narrow Escape Problem
Of particular interest to the field of mathematical biology is the so-called Narrow Escape Problem
(actually a class of problems), which investigate the asymptotic behavior of the MFPT in domains
which small absorbing ‘windows.’ The most studied of such problems is that of a disk which has a
small absorbing window of arc length 2ϵin its perimeter. In this section we show two approaches
to the problem: we first derive and match the setup for the narrow escape problem in [1] in the real
variables before following an approach by Chen and Friedman [6] which makes use of conformal
mappings in the complex plane to produce the desired asymptotics.
9
3.1 Real Variables
The dimensionless (lacking diffusion constant) MFPT is given by the mixed Neumann-Dirichlet
boundary value problem
∆v(r, θ) = −1, r < 1,for 0 ≤θ < 2π,
v(r, θ)|r=1 = 0,for |θ−π|< ϵ, (30)
∂v(r, θ)
∂r r=1
= 0,for |θ−π|> ϵ.
As demonstrated in Singer et al., the substitution
u=v−1−r2
4
reduces (30) to a Laplace equation problem
∆u(r, θ)=0, r < 1,for 0 ≤θ < 2π,
u(r, θ)|r=1 = 0,for |θ−π|< ϵ, (31)
∂u(r, θ)
∂r r=1
=1
2,for |θ−π|> ϵ.
Separating variables, we have
u(r, θ) = R(r)Θ(θ)
and thus
R′′Θ + R′Θ
r+RΘ′′
r2= 0.
Multiplying by r2
RΘwe are left with the eigenvalue problem
r2R′′
R+rR′
R=−Θ′′
Θ
Treating the Θ equation first, we have λΘ+Θ′′ = 0.
We know consider the three possible cases for these two eigenvalue problems.
λ=0: Then for some constants aand b,
Θ(θ) = aθ +b.
But because we are working on the disk, Θ must be 2πperiodic, so this this cannot be our solution.
λ<0: In this case, letting λ=−µwe have
Θ′′ =µΘ.(32)
As this is a second order, linear, homogeneous ODE, we look for solutions of the form Θ(θ) = eβθ .
Assuming a solution of this form, it can be seen that
eβθ(β2−µ) = 0,
and therefore the characteristic equation has two real roots β=±√µwhich suggests the ansatz
Θ(θ) = ae√µθ +be−√µθ +c
10
for some constants a, b, and c. But again our solution must be periodic in θ, which this solution is
not, so we reject the λ < 0 case.
λ>0: As before, we let λ=µ2and look for solutions of the form Θ(θ) = eβθ. Assuming a
solution of this form, it can be seen that
eβθ(β2+µ2)=0,
and therefore the characteristic equation has two real roots β=±iµ which suggests the ansatz
Θn(θ) = Ancos(µnθ) + Bnsin(µnθ).
The natural boundary conditions for polar coordinates require that Θ(0) = Θ(2π) and Θ′(0) =
Θ′(2π). Imposing these conditions we see that
An=Ancos(µn2π) + Bnsin(µn2π)
and
Bn=Bncos(µn2π)−Ansin(µn2π).
Thus An=Bnsin(µn2π)
1−cos(µn2π)so we have
1 = cos(µn2π)−sin2(µn2π)
1−cos(µn2π),
and thus
1−cos(µn2π) = cos(µn2π)(1 −cos(µn2π)−sin2(µn2π).
Finally this gives us 1 = cos(µn2π) so we have µn=n.
Now considering the radial function Rn(r) for λ=n2>0 our differential equation becomes
r2R′′
n+rR′
n−n2Rn(33)
which is an Euler Differential equation whose ansatz is derived by the change of variables t= ln(r)
and a function φsuch that R(r) = φ(ln(r)) = φ(t).Consequently,
R′=1
r
dφ
dt
and
R′′ =1
r2 d2φ
dt2−dφ
dt !.
Then in our new coordinates (33) becomes
d2φ
dt2−dφ
dt +dφ
dt −nφ(t) = d2φ
dt2−n2φ(t) = 0
which is the same eigenvalue problem as (32) so we know we have ansatz of the form
φ(t) = aent +be−nt (34)
for some constants aand b. Then converting back to the original radial coordinates we see that
R(r) = arn+br−n.
11
Although it may seem as though the radial function has only one boundary condition, the natural
boundary condition is that
lim
r→0R(r)=∞,
so we discard the r−ncomponent of the solution. Therefore particular solutions are given by
un(r, θ) = rn(Ancos(nθ) + Bnsin(nθ)),(35)
and the general solution is a linear combination of all particular solutions
u(r, θ) = 1
2A0+∞
X
n=1
rn(Ancos(nθ) + Bnsin(nθ)).(36)
By the construction of our problem, the symmetry requires an even solution in θ, so we know
that Bn= 0.Then imposing boundary conditions we are left with the deal series relations
1
2A0+∞
X
n=1
Ancos(nθ) = 0 for |θ−π|< ϵ (37)
∞
X
n=1
nAncos(nθ) = 1
2for |θ−π|> ϵ. (38)
This is essentially the start of the method presented in [1] which relies heavily on the theory of
special functions. In part due to the lack of progress with this dual series, a method using conformal
mappings is given next.
3.2 An Approach in the Complex Plane
After a first course in complex analysis, it seemed worthwhile to pursue a method developed by
Chen and Friedman [6] which makes use of conformal mappings to map our problem into a simpler,
unbounded domain where a variety of approaches are possible.
Consider the same problem in the complex plane, where our boundary value problem is defined
on the domain Ω = {reiθ |r < 1,−π≤θ≤π}with boundary ∂Ω
∂Ωa={eiθ | − ϵ≤θ≤ϵ},0< ϵ ≪1.
We also set a=eiϵ,¯a=e−iϵ.
Just as before, we define our boundary value problem for our solution v:
∆v=−1 in Ω, vn= 0 on ∂Ω\∂Ωa, v = 0 on ∂Ωa.(39)
Similar to the substitution used in Singer et al., consider writing vin the form
v(z) = 1− |z|2
4+ ln |1−z| − ϕ(z).
Then ϕsatisfies
∆ϕ= 0 in Ω, ϕn= 0 on ∂Ω\∂Ωa, ϕ(z) = g(z) on ∂Ωa,(40)
where g(z) = ln |1−z|.Following the method in [6], this Dirichlet-Neumann problem will be solved
for a general gby way of conformal mappings.
12
We begin by considering the conformal mapping
ζ:= a−z
1−az .(41)
All fractional linear transformations such as this can be written as a composition of four transfor-
mations, i.e. ζ=T4◦T3◦T2◦T1(z),where in our case
T1(z) = z+1
−a=z−¯a, T2(z) = 1
z
T3(z) = −a2+ 1
a2(z),and T4(z) = z+1
a=z+ ¯a.
Note that we have used that fact since ¯a=1
a.The motivation for the use of this particular
transformation is more easily understood through an examination of these four transformations.
First, the unit circle in Cis translated to have center at −¯aby T1(z).Upon composition with
the inversion transformation T2(z),¯a→ ∞, the unit circle is mapped to a half plane in the second,
third, and fourth quadrants, with boundary marked by a line in the plane as illustrated below.
Figure 1
Unit Circle Image under T1(z) Image under T2◦T1(z)
Upon further composition with T3, this half plane is rotated and scaled such that ∂Ω is mapped
a line parallel to the real axis. Finally the mapping T4(z) translates our half plane by ¯asuch that
ζ(a) = 0, ∂Ωais mapped to the negative real axis, and ∂Ω\∂Ωais mapped to the positive real
axis.
Next we look to extend ϕ(ζ(z)) by ‘even’ reflection from Im(ζ)>0 to Im(ζ)<0 such that
ϕ(¯
ζ) = ϕ(ζ).Since ϕsatisfies Laplace’s equation in Ω, is it harmonic under the mapping ζfor
Im(ζ)>0 and is continuous on the positive real axis with zero normal derivative (along the
imaginary ζaxis). While time did not permit an exploration of this advanced concept, the ‘even’
extension of ϕ(ζ) to the negative half plane can be shown to be harmonic first in the lower half plane,
and then along the positive real axis. A variety of sources suggest proving these two results related
to the Schwarz Reflection Principle using the mean value and half-ball mean value property[7,8].
The Neumann condition essentially ensures the extended ϕto be have continuous first derivative
along the positive real axis which allows ϕto be harmonic on this line. After this extension, ϕis
harmonic on the whole ζ-plane except the negative real axis, as we much have the same boundary
value g(ζ).We then consider the conformal mapping
η=pζ=ra−z
1−az =: η(z)
where we have specified our branch of the complex logarithm to be from −πto π. With this choice
of branch, the even extension of the ζplane is mapped onto the half plane ℜ(η)>0 and the
13
two sides of the negative real ζaxis are mapped onto the imaginary axis of the η-plane. Thus by
the even extension of ϕ(ζ) and the ηmapping, we are essentially left with a half-plane Dirichlet
problem instead of a mixed boundary value problem. Note too that η(ζ(z)) maps points from the
original domain ¯
Ω to the η-plane, and thus because ηhas an inverse which is one-to-one given by
η−1(z) = a−z2
1−az2,(42)
η−1(z) maps points from the η-plane to ¯
Ω.Thus our construction of ηand determination of
its inverse has allowed us to rewrite our boundary value problem as a more familiar half-plane
boundary value problem which we will approach using the classical Fourier Transform.
Our boundary value problem in the complex η-plane is
∆ϕ= 0 for ℜ(η)>0, ϕ(η−1(it)) = g(η−1(it)) for t∈R.
First, we make note of the use and motivation for the Fourier Transform. The Fourier Transform
is a continuous extension of the Fourier series which involves a continuous set of eigenfunctions.
The classical transform is invertable and the Fourier Transform ˆ
f(ξ) of a real-valued function f(x)
is defined by
ˆ
f(ξ) = Z∞
−∞
e−ixξf(x)dx.
The transform is said to be invertable because if a functions transform ˆ
fis known, one can easily
recover fwith
f(x) = 1
2πZ∞
−∞
eixξ ˆ
f(ξ)dξ.
Because the integral transformation is of another variable, one can easily compute derivatives
(where Fdenotes the operator notation of the transform):
F[f′(ξ)] = iξ ˆ
f(ξ) and
f′(x) = 1
2πZ∞
−∞
iξ ˆ
f(ξ)eixξdξ.
Note that in two-dimensions
ˆ
f(ξ1, ξ2) =
∞
Z
−∞
∞
Z
−∞
eixξ1+iyξ2f(x, y)dxdy.
The Fourier Tranform transforms a function into a new function which describes the frequencies
of the original function. The transform ˆ
fis a complex-valued function of these frequency. Historical
motivation for the FT is in part due to the difficulty of imposing boundary conditions on unbounded
regions.
Returning to the problem at hand, if we consider the real and imaginary axes of the η-plane be
uand v, we can see that ∆ϕ=ϕuu +ϕvv = 0 noting that u, v ∈Rand u > 0.If we then take the
Fourier transform of this equation in the vvariable (as required by our boundary condition along
the v-axis [9]), we have
Fv[ϕuu +ϕvv = 0] (43)
ˆ
ϕuu(u, ξ2)−ξ2
2ˆ
ϕ(u, ξ2) = 0 (44)
14
with the boundary condition ˆ
ϕ(0, ξ2) = F[g(η−1(iξ2))].With ξ2fixed, the general solution to this
homogeneous, second order ordinary differential equation is
ˆ
ϕ(u, ξ2) = A(ξ2)e|ξ2|u+B(ξ2)e−|ξ2|u.
The exponential growth of the first term requires its omission, so setting u= 0 to determine B(ξ2),
we have
ˆ
ϕ(u, ξ2) = F[g(η−1(iξ2))]e−|ξ2|u.
We can then invert our formulation to give
ϕ(u, v) = 1
2πZ∞
−∞ F[g(η−1(iξ2))]e−|ξ2|ueivξ2dξ2.(45)
Similarly we can invert the transform of the boundary condition within the integral to find
ϕ(u, v) = 1
2πZ∞
−∞ Z∞
−∞
g(η−1(it))e−iξ2tdte−|ξ2|ueivξ2dξ2.(46)
Assuming the functions meet the advanced criteria for smoothness and are integrable, an inter-
change of integration yields
ϕ(u, v) = 1
2πZ∞
−∞
g(η−1(it))Z∞
−∞
e−|ξ2|u+iξ2(v−t)dξ2dt (47)
for u > 0.To compute this inner integral, we split the domain of integration into two regions from
−∞ to 0 and from 0 to ∞.Thus we have
Z0
−∞
eξ2u+iξ2(v−t)dξ2+Z∞
0
e−ξ2u+iξ2(v−t)dξ2
=Z0
−∞
eξ2(u+i(v−t))dξ2+Z∞
0
e−ξ2(u−i(v−t))dξ2
=eξ2(u+i(v−t)
u+i(v−t)0
−∞
+e−ξ2(u−i(v−t)
−(u−i(v−t))∞
0
=1
u+i(v−t)+1
u−i(v−t)
=u−i(v−t) + u+i(v−t)
u2+ (v−t)2=2u
u2+ (v−t)2.
Writing η−1(it) explicitly as in (42), we then have
ϕ(u, v) = 1
2πZ∞
−∞
ga+t2
1 + at22u
u2+ (v−t)2dt. (48)
Since our boundary condition holds on the imaginary axis under η, we know η−1(it) for t∈Rgives
{z|z∈∂Ωa}, so from the definition of g(z) = ln |1−z|we can rewrite our boundary condition
g(η−1(it)) as
g(a+t2
1 + at2) = ln
1 + at2−a−t2
1 + at2
= ln
(1 −a)(1 −t2)
1 + at2
= ln
1−a
2a·2(1 −t2)
1 + t2·1 + t2
¯a−1 + 1 + t2
.
Using this formulation (48) can be written as
ϕ(u, v) = f1+f2−f3,(49)
15
where
f1=1
πZ∞
−∞
u
u2+ (v−t)2ln
1−a
2a
dt, (50)
f2=1
πZ∞
−∞
u
u2+ (v−t)2ln
2(1 −t2)
1 + t2
dt, (51)
f3=1
πZ∞
−∞
u
u2+ (v−t)2ln
1 + ¯a−1
1 + t2
dt. (52)
We first address f1by simplifying to give the form of a familiar trigonometric integral after the
substitution s=v−t
u, ds =−1
u,
f1=1
πu ln
1−a
2aZ∞
−∞
1
1+(v−t
u)2dt =1
πln
1−a
2aZ∞
−∞ −1
1 + s2ds (53)
=1
πln
1−a
2a−arctan v−t
u∞
t=−∞
=1
πln
1−a
2aπ
2+π
2.(54)
Thus we have
f1= ln
1−a
2a
= ln
1
2e−ϵ−1
2
= ln
1
2(cos(−ϵ) + isin(−ϵ)) −1
2
(55)
= ln 1
2(cos(−ϵ)−1
2+i
2sin(−ϵ))
(56)
= ln s1
2(cos(−ϵ)−1
22
+1
2sin(−ϵ)2
(57)
= ln r1
4cos2(−ϵ)−1
2cos(−ϵ) + 1
4+1
4sin2(−ϵ) (58)
= ln r1−cos(ϵ)
2= ln sin ϵ
2.(59)
To handle f2and f3, we will only consider the asymptotics from the center of the circle, at the
point η(0) = √a= cos( ϵ
2) + isin( ϵ
2).Thus u= cos( ϵ
2) and v= sin( ϵ
2).Thus our kernel can be
written as
u
u2+ (v−t)2=cos ϵ
2
cos2(ϵ
2) + (sin( ϵ
2)−t)2(60)
=cos( ϵ
2)
cos2(ϵ
2) + sin2(ϵ
2)−2tsin( ϵ
2) + t2(61)
=cos( ϵ
2)
1−2tsin( ϵ
2) + t2.(62)
Note both numerator and denominator are positive for all ϵ < π, and the denominator is minimized
when at t= sin ϵ
2which approaches 0 as ϵ→0.Thus (62) is positive for all twith ϵ < π and is
bounded above by
cos( ϵ
2)
1−sin2(ϵ
2)(63)
which is bounded by some constant for all ϵ≤π
2and approaches 1 as ϵ→0.Therefore for any
ϵ < π
2, the Poisson kernel for the image of the origin under ηhas asymptotics O(1).
Thus the asympotics for f2at the origin of the circle are
f2=1
πZ∞
−∞ O(1) ln
2(1 −t2)
1 + t2
dt =1
πZ∞
−∞ O(1) ln
1−t2
1 + t2
dt (64)
=O(1)
πZ−1
−∞
ln t2−1
1 + t2dt +Z1
−1
ln 1−t2
1 + t2dt +Z∞
1
ln t2−1
1 + t2dt.(65)
16
Due to the weakness of the logarithmic singularities, all three of these integrals converge, though
time has not permitted a full explanation. From this we can conclude that f2=O(1).
For f3we can then write the asympotics for the logarithmic component to give
|f3|=1
πZ∞
−∞
(1)O(1)|¯a−1|
1 + t2dt. (66)
Notice that
|¯a−1|=q(cos(−ϵ)−1)2+ sin2(−ϵ) = p2−2 cos(ϵ) = 2 sin( ϵ
2) = O(ϵ)
Therefore we have
|f3| ≤ 1
πZ∞
−∞
O(ϵ)
1 + t2dt =1
πO(ϵ) arctan(t)
∞
−∞
=O(ϵ).(67)
Thus from the center of the unit circle, we have matched the asymptotics in [6]to show that
ϕ(u, v) = ln sin ϵ
2+O(ϵ) + O(1).(68)
Therefore from the preceding substitution, we have shown
v(z) = 1− |z|2
4+ ln |1−z|
sin ϵ
2
+O(ϵ) + O(1).
4 Development and Results of the Random Walk Numerical
Method
As a way to compare the different paradigms of the random walk and the analytic results
in the literature, a vectorized random walk model was build using the python library NumPy.
With its core written in the compiled language C, the library provides incredible computational
power to Python and as such is the backbone for most modern computational science in Python.
The NumPy based, vectorized random walk was motivated out of a desire to gain familiarity with
the core data structures of NumPy, namely n-dimensional arrays, and as an exercise creating a
potentially more powerful and elegant model. While such random walks would be quite easy to
implement with many nested iterables and conditionals, the computational inefficiency of such an
approach makes a more complicated, vectorized model attractive.
The core of our model is the use of a 3D NumPy array and framework for a vectorized random
walk based on ‘folding’ our array in the plane back into the unit square [10]. First, we initialize a
2D array Awith shape (N, d), i.e. Nrows and dcolumns. The Nrows hold the Nstarting points
of our random walk and the dcolumns hold the dcomponents of the d-dimensional points. Clearly
in a two dimensional random walk, d= 2 and the two columns hold the xand yvalues of these
starting positions. We then generate a 3D array Fof shape (N, d, S), where Sis the number of
steps we have chosen to take in our random walk. We generate this array by randomly assigning one
of the 2d(the quantity, not a dimension), d-dimensional orthogonal basis vectors to every location
in an array of shape (N, S), and then reassigning these values to their cumulative sum along the S
17
axis. Thus we are left with our array Fof shape (N, d, S). Then, via array broadcasting, we add
our starting array A‘slice-wise’ along the S-axis of F. From here, essentially modulus our values
in Fby the particular boundary constraints. In the case of a two dimensional walk on the unit
square in the first quadrant, we could simply take the modulus of all values in Fby 1 to ‘fold’ our
random walk inside the unit square. We essentially impose reflective boundary conditions on the
entire domain which allows us to search through the array to identify first passage times.
Figure 2
xT1,1xT1,2
xT2,1xT2,2
xT3,1xT3,2
.
.
..
.
.
xTN,1xTN,2
.
.
..
.
.
.
.
..
.
.
.
.
..
.
.
.
.
..
.
.
.
.
..
.
.
x31,1x31,2
x32,1x32,2
x33,1x33,2
.
.
..
.
.
x3N,1x3N,2
x21,1x21,2
x22,1x22,2
x23,1x23,2
.
.
..
.
.
x2N,1x2N,2
x11,1x11,2
x12,1x12,2
x13,1x13,2
.
.
..
.
.
x1N,1x1N,2
Time
Structure of the 3D array where the ‘Time’ axis corresponds to the Saxis and d= 2.
In seeking to compare analytical and numerical results, some attention must be paid to differ-
ence in scaling between the discrete and continuous. We note that in many random walk models
such as those in [11], chemical reaction or random walks firstly based in time, and the particular
reaction of movement of a particle does not occur in a predictable or fixed way. For instance in
[11] (8), the position of a particle in a one dimensional random walk is computed by
X(t+ ∆t) = X(t) + √2D∆tζ,
where ∆tis the fixed time increment, Dis the diffusion constant, and ζis a normally distributed
variable (mean of zero and unit variance). From this definition and the work in [12], we can think
of the mean displacement in one dimension as ∆¯x=√2D∆t.
In our model, our particle ‘hops’ along a fixed, regular mesh-grid with a particle always moving
one of four directions in the Cartesian plane, so as in [12], the mean displacement in two dimensions
is given by ∆¯x=√4D∆t. With these two relationships established, to compare analytical solutions
with units of the MFPT dictated by the diffusion constant being equal to 1 (as we have used the
dimensionless MFPT) we must multiply analytical solutions by
2
(∆x)2and 4
(∆x)2
in one and two dimensions respectively to compare our analytical solutions to our numerical results
which have units ∆t. Of particular importance is the coefficient R2for radius Rin the asymptotics
in [1], which in in our case is ( 1
2)2because the natural choice of a conformal map from a circle to
the unit square would have the circle to have radius 1
2as the unit square has side length of one.
18
4.1 One Dimensional Results for Confirmation
To first demonstrate the functionality of the model and its agreement with our analytical solutions,
the one dimensional random walks with Dirichlet and Robin boundary conditions were compared
to the results in section 2.1.1 and 2.1.2 (in the case that b= 1).
Figure 3
In the above numerical results, our partition of the subset of the real line [0,1] made use of
∆x= 0.02, so by referencing the calculations in sections 2.1.1 and 2.1.2, we found the scaled
Dirichlet analytical solutions to be
2
0.022−x2
2+1
2,
and the scaled Robin analytical solution to be
2
0.022−x2
2+1
2+σ,
where σis our Robin parameter whose relationship to ‘absorption’ probability 1 −qis given in [13]
as
σ=(∆x)q
1−q.
19
Here qis the ‘splitting probability’ or the chance that a particle is not absorbed at the boundary.
Note that in the figure above we have used q= 0.9 (i.e. σ= 0.18), which means particles are
absorbed at x= 1 only 10% of the time.
As seen in the above figure, numerical results for both the Robin and Dirichlet boundary
condition at x= 1 show good agreement with the analytical results developed in sections 2.1.1 and
2.1.2. The above numerical results were computed with a sample size of n= 30000,and the small
variance seen in the results would diminish with a larger sample size.
4.2 The Classical Dirichlet Narrow Escape Problem
Of great interest to this thesis was the numerical approximation of asympotics for solution to the
Narrow Escape Problem. In analytical settings our interest is in smooth domains like the unit disk,
however such domains are not nearly as conducive to random walks as are structured, rectangular
domain like the unit square. We can rely on the fact that the Laplacian operator is conformal-
invariant, and thus a mapping from the unit disk to the unit square yields the same asymptotic
results. Thus for the ease and practicality of computations, our numerical model for this problem
conducted a random walk on the unit square with a narrow opening in the middle of the x= 1
edge, with equation ∂Ωa={(x, y)|x= 1,|y−0.5|< ϵ}.
With the appropriate scaling procedure detailed in the preceding section, the following figures
compare our numerical results to the asymptotics (with the error removed for simplicity) developed
in section 3.2.
Figure 4
(a) First attempt with ∆x= 0.003125 (b) Refined data with ∆x= 0.001
20
(c) Compilation with refined data for ∆x= 0.001
The above figures suggest our model provides a useful approximation of the asymptotic behavior
of the MFPT but is limited in its accuracy. In the first attempt shown in (a), the model shows
poor agreement for larger ϵwhich is to be expected as the asymptotics derived in the literature
and the preceding section only hold for ϵ≪1.Though we would expect smaller ϵto yield more
accurate results, for ϵ < 0.04, our model exhibits increasingly poor estimates of the MFPT as best
seen in (b). The compilation of the data seen in (c) shows our model is capturing the rough trend
in behavior, but the decreasing agreement between model and analytical solution for very small ϵ
seen in (b) is a curious and interesting phenomena.
With our model running a finite number of steps, there is always a possibility of a partic-
ular “particle” never reaching the absorption window. Without more advanced statistical anal-
ysis to account for such random walks, we chose to o,it such walks with our calculation of the
MFPT and instead report the distribution of absorbed particles for the four smallest epsilons:
0.001,0.002,0.003,and 0.004.The data is shown below in Figure 5.
21
Figure 5.
Distribution of Absorbed Particles by Epsilon
As can be seen, the choice of 80,000,000 steps in our model (which was decided up after many
trials) was enough to absorb the vast majority of random walks such that the decreasing effect of
omitting walks never absorbed on the MFPT is evidently minimal.
4.3 A Comparison with a Robin Narrow Escape Problem
While literature to serve as a basis for comparison has not been established, the model was
extended to include the case of a reactive, Robin BC absorbing window.
22
Figure 6.
Narrow escape with small Robin condition window
Note that in all four panels the black dots are the results for the Dirichlet problem, while the
colored crosses are the corresponding Robin results. Reading left to right, top to bottom, we can
see the first two figures show fairly minimal upward translations or magnification of the Dirichlet
result, which is to be expected for these fairly high reaction probabilities as we have low q(splitting
probability) values. In the third figure, we see a notable magnification of the the Dirichlet results,
as the Robin data points are not simply translated by a constant to higher MFPT values as was
seen in the one-dimensional case. The fourth figures serves to give a sense of scale by compiling
all data.
5 A Comparison of Numerical and Analytic Solutions
In general, the figures demonstrate moderate agreement with the analytical results derived by
our earlier work and the more involved results in the literature. While our model continues to be
refined, it was been shown to be a reliable tool for demonstrating both the accuracy of analytical
solutions and also the deep connections between random walks and diffusion-based processes.
Of particular interest is the apparent difference in behavior of the 2D Robin results for the
23
Narrow escape problem and the 1D Robin results on the interval [0, b].It appears that the additional
dimension in the narrow escape problem has a magnifying effect on the MFPT. In many ways it
is no surprise that we observe non-constant change in the MFPT with respect to σof qin this
scenario. When a particle is deflected at the reactive window, it has half the chance of returning for
another attempt in the following step compared to the one-dimensional case. This means particle
deflections are in a sense more ‘consequential’ and delay absorption by high orders than in one-
dimension. Particularly noticeable for high qvalues, it would seem that the the Robin Narrow
Escape problems features a more-pronounced singularity as ϵ→0.
In our comparison between the discrete approximations and continuous asymptotics to the
Narrow Escape Problem, of great concern is the divergence of agreement as ϵ→0 in Figure 4
(b). This disagreement could have arisen from a variety of places, including (a) methodological
error (despite the many reviews of the computations), (b) a simple lack of computational power
to narrow our window to an even greater extent, and (c) a more fundamental difference between
diffusion and random walks in higher dimensions as mentioned in [2].
There is a chance that a minimum ϵof 0.001 as was computed in our case is simply not a small
enough window to create behavior in agreement with the asymptotics (Indicating case (b)). In
this case a significant increase in computational power would be required for model agreement.
This would likely mean that the slopes of analytical and numerical results would have to differ by
enough to allow for their convergence for ϵ < 0.001.
Case (c) is a phenomena mentioned in [2] (which is an undergraduate paper), but despite efforts
to follow the resources they cite, I have yet to find evidence for this vague side-note. If there is
such an in-congruence between the diffusion equation and random walks in higher dimensions, this
would certainly limit the two’s comparison, but we have already shown that regardless, there are
undeniable connections and fruitful comparisons to make. As an efficient computational tool, the
numerical model has proven to be a useful proxy for the interesting diffusion-based boundary value
problems considered in this thesis.
5.1 Future Directions
Over the course of this project spanning most of my senior year, there are a host of problems
which arose and a host of deep facts which have been taken for granted. For much of the process,
coming to an understanding of a solution of the dual series (37) was a large analytical goal.
Despite many attempts at our own solution, and attempts to understand the referenced work,
the background required for their analysis proved too involved. After a first course in complex
analysis, the approach using conformal mappings in [6] was of great interest so the last few weeks
have been spent working through their solution. In working through this solution, again I was met
with a variety of fundamental concepts which I did not have time to fully explore. Of these include
the conformal invariance of the Laplacian Operator, a variety of deep facts about generalized
eigenfunctions in ‘smooth’ domains, and the Schwarz reflection principle for harmonic and analytic
functions. The work in this paper relies in several crucial places on these advanced techniques,
and gaining more familiarity with these concepts would be of interest to me in the future.
24
In considering extensions and different approaches to numerical schemes, it would be a great
exercise to be become introduced to the finite volume method, as used in [2]. More advanced
papers with greater computational experience and resources have implemented extremely powerful
fast-solvers and finite element methods for domains of arbitrary complexity. While such approaches
are entirely out of reach for this thesis, it is interesting to note the variety of methods which have
been employed to solve different versions of these problems.
6 Conclusion
This thesis aims to have introduced the reader to the stochastic world of random walks and its
remarkable connections to diffusion. After deriving this connection explicitly, notions of first pas-
sage times and mean first passage times were expanded upon and illustrated mathematically. The
direct, first principles approach to a MFPT from the heat equation has been shown to be equivalent
to the much simpler Poisson equation ∆V=−1
kwith appropriate boundary conditions. From this
foundation, a few interesting examples of MFPT problems in one, two, and three dimensions were
illustrated, before deriving the dual series relations for the Narrow escape problem. Despite our
lack of an analytical solution to this dual series, and alternate method in the complex plane using
conformal mappings has been presented. In conjunction with this analytical work, a vectorized
random walk model using 3D NumPy arrays served as a basis for comparison between discrete
random walks and the asymptotics found in [1] and seen in section 3.2. After this comparison,
we conducted an exploration of a Robin boundary condition narrow escape problem which showed
behavior not seen in the one-dimensional Robin MFPT problem.
While we initially underestimated the complexities of the narrow escape problem, this thesis
has spurred a great amount of problem solving and learning since its inception in the summer of
last year.
25
7 References
[1] Singer, A., Schuss, Z.; Holcman, D. (2004, December 15). Narrow escape, part II: The circular
disk. arXiv.org. https://arxiv.org/abs/math-ph/0412050. Accessed 21 April 2023
[2] Maurer, N. (2020). Solving 2D Unconditional MFPT Problems with Arbitrary Confining Ge-
ometries. (T. Roberts, Ed.). AMSI Vacation Research Scholarships. Accessed 21 April 2023
[3] Polizzi, N. F., Therien, M. J.; Beratan, D. N. (2016, August 5). Mean first-passage times in
biology. Israel Journal of Chemistry. U.S. National Library of Medicine. https://www.ncbi.nlm.
nih.gov/pmc/articles/PMC5656268/#:~:text=The\%20mean\%20first\%2Dpassage\%20time\%20(MFPT)
\%20defines\%20an\%20average,started\%20at\%20some\%20initial\%20state. Accessed 21
April 2023
[4] Bell, J. An Alternative Heat Equation Derivation. http://www.math.umbc.edu/~jbell/pde_
notes/K_Random\%20Walk\%20Heat\%20Equation\%20Derivation.pdf. Accessed 21 April 2023
[5] Xinwei , Y. 1.1. Heat equation. 1.1.1. 1D Random walk. - ualberta.ca. https://www.math.
ualberta.ca/~xinweiyu/436.A1.12f/PDE_Sources_PDE.pdf. Accessed 21 April 2023
[6] Chen, X. C.; Friedman, A. (2011). Asymptotic analysis for the narrow escape problem. SIAM
Publications Library. https://epubs.siam.org/doi/10.1137/090775257. Accessed 21 April
2023
[7] Fok, A. (2015, May 11). Schwarz reflection principle for harmonic functions. Mathematics
Stack Exchange. https://math.stackexchange.com/questions/1276993/schwarz-reflection
-principle-for-harmonic-functions. Accessed 21 April 2023
[8] Sargera; Negro, G. (2013, February 14). Schwarz reflection principle for (real) harmonic func-
tions. Mathematics Stack Exchange. https://math.stackexchange.com/questions/302509/
schwarz-reflection-prinnciple-for-real-harmonic-functions. Accessed 21 April 2023
[9] Stein. The Fourier Transform on R. Radboud Universiteit. https://www.math.ru.nl/~mueger/
Stein4.pdf. Accessed 21 April 2023
[10] Panzer, P. (2018, February 14). Vectorized random walk in python with boundaries. Stack
Overflow. https://stackoverflow.com/questions/48777345/vectorized-random-walk-in
-python-with-boundaries. Accessed 21 April 2023
[11] Erban, R.; Chapman, S. J. (2007, February 15). Reactive boundary conditions for stochas-
tic simulations of reaction ... IOP Science Physical Biology. https://iopscience.iop.org/
26
article/10.1088/1478-3975/4/1/003. Accessed 21 April 2023
[12] Geller, B.; Redish, J. (2011, December 1). Diffusion and Random Walks. Diffusion and
random walks - Nexus Wiki. https://www.compadre.org/nexusph/course/Diffusion_and_
random_walks. Accessed 21 April 2023
[13] Bressloff, P. C. (2022, January 10). The narrow capture problem: An Encounter-based ap-
proach to partially reactive targets. The narrow capture problem: an encounter-based approach to
partially reactive targets. https://arxiv-export1.library.cornell.edu/abs/2201.01675v2.
Accessed 21 April 2023
[14] Grebenkov, D. S. (2016, December 23). Universal formula for the mean first passage time in
planar domains. Physical Review Letters. American Physical Society. https://journals.aps.
org/prl/abstract/10.1103/PhysRevLett.117.260201. Accessed 21 April 2023
8 Appendix
All finalized Python code written for this thesis is available on GitHub at https://github.com/
henrynjones/MFPT-Thesis.git.
27
