Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The narrow escape problem : a matched asymptotic expansion approach Pillay, Samara 2008

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata

Download

Media
[if-you-see-this-DO-NOT-CLICK]
ubc_2008_fall_pillay_samara.pdf [ 907.98kB ]
Metadata
JSON: 1.0066535.json
JSON-LD: 1.0066535+ld.json
RDF/XML (Pretty): 1.0066535.xml
RDF/JSON: 1.0066535+rdf.json
Turtle: 1.0066535+rdf-turtle.txt
N-Triples: 1.0066535+rdf-ntriples.txt
Original Record: 1.0066535 +original-record.json
Full Text
1.0066535.txt
Citation
1.0066535.ris

Full Text

The Narrow Escape Problem: AMatched Asymptotic ExpansionApproachbySamara PillayB.Sc.(Hons.), The University of the Witwatersrand, 2005A THESIS SUBMITTED IN PARTIAL FULFILMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMaster of ScienceinThe Faculty of Graduate Studies(Mathematics)The University Of British Columbia(Vancouver)August, 2008c Samara Pillay 2008AbstractWe consider the motion of a Brownian particle trapped in an arbitrary boundedtwo or three-dimensional domain, whose boundary is re?ecting except for a smallabsorbing window through which the particle can escape. We use the method ofmatched asymptotic expansions to calculate the mean ?rst passage time, de?nedas the time taken for the Brownian particle to escape from the domain throughthe absorbing window. This is known as the narrow escape problem. Since themean escape time diverges as the window shrinks, the calculation is a singularperturbation problem. We extend our results to include N absorbing windowsof varying length in two dimensions and varying radius in three dimensions. Wepresent ?ndings in two dimensions for the unit disk, unit square and ellipse andin three dimensions for the unit sphere. The narrow escape problem has variousapplications in many ?elds including ?nance, biology, and statistical mechanics.iiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ixDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Brownian Motion . . . . . . . . . . . . . . . . . . . . . . . . . . 11.2 The Narrow Escape Problem . . . . . . . . . . . . . . . . . . . . 41.2.1 Statement of the Problem . . . . . . . . . . . . . . . . . 41.2.2 Derivation of the Model Equation . . . . . . . . . . . . . 51.3 Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . 81.4 Literature Review - Biological Context . . . . . . . . . . . . . . 101.4.1 Two Dimensions . . . . . . . . . . . . . . . . . . . . . . . 111.4.2 The Narrow Escape Problem in Three Dimensions . . . . 151.5 Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162 The Narrow Escape Problem in Two Dimensions . . . . . . . 182.1 Formulation of the Problem . . . . . . . . . . . . . . . . . . . . 192.2 Two-Dimensional Domain with One Hole on the Boundary . . . 212.2.1 Calculation of ?0 and ?0 . . . . . . . . . . . . . . . . . . 212.2.2 Calculation of the Mean First Passage Time . . . . . . . 282.3 Two-Dimensional Domain with N Holes on the Boundary . . . . 292.3.1 Calculation of ?0 and ?0 . . . . . . . . . . . . . . . . . . 302.3.2 Calculation of the Mean First Passage Time . . . . . . . 38iiiTable of Contents2.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 403 Numerical Realizations: The Neumann Green's Function . . 413.1 The Unit Disk . . . . . . . . . . . . . . . . . . . . . . . . . . . . 413.1.1 The Neumann Green's Function in a Unit Disk . . . . . . 413.1.2 The Mean First Passage Time in a Unit Disk with OneHole on the Boundary . . . . . . . . . . . . . . . . . . . . 433.1.3 The Mean First Passage Time in a Unit Disk with TwoHoles on the Boundary . . . . . . . . . . . . . . . . . . . 463.1.4 Equally spaced Points on a Unit Disk . . . . . . . . . . . 483.2 The Unit Square . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.2.1 The Neumann Green's Function in a Unit Square . . . . 513.2.2 Solution Method . . . . . . . . . . . . . . . . . . . . . . . 523.2.3 The Mean First Passage Time in a Unit Square with OneHole on the Boundary . . . . . . . . . . . . . . . . . . . . 573.2.4 The Mean First Passage Time in a Unit Square with TwoHoles on the Boundary . . . . . . . . . . . . . . . . . . . 623.3 Arbitrary Domains - A Numerical Approach . . . . . . . . . . . 663.3.1 The Boundary Element Method . . . . . . . . . . . . . . 673.3.2 The Unit Disk . . . . . . . . . . . . . . . . . . . . . . . . 713.3.3 The Ellipse . . . . . . . . . . . . . . . . . . . . . . . . . . 733.4 Perturbed Circular Domains - Analytical Approach . . . . . . . 793.4.1 Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . 793.4.2 Example . . . . . . . . . . . . . . . . . . . . . . . . . . . 863.5 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 894 The Narrow Escape Problem in Three Dimensions - The Sphere924.1 Derivation of the Neumann Green's function for a Sphere . . . . 924.2 Three-Dimensional Sphere with N Absorbing Patches on theBoundary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 974.2.1 Calculation of ?0 and ?0 . . . . . . . . . . . . . . . . . . 984.2.2 Calculation of the Mean First Passage Time . . . . . . . 1084.3 Three-Dimensional Sphere with N Absorbing Patches on theBoundary - An Alternate Derivation . . . . . . . . . . . . . . . . 1114.3.1 Calculation of the MFPT Directly . . . . . . . . . . . . . 1114.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118ivTable of Contents5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124AppendicesA Derivation of d . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128B Far-Field Behaviour of ?c . . . . . . . . . . . . . . . . . . . . . . . 130C Matlab Code . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132vList of Tables3.1 R(x;x0) for the unit disk with N mesh points on the boundary . 713.2 R(x0;x0) for the unit disk with N mesh points on the boundary 723.3 R(x;x0) for the ellipse with N mesh points on the boundary . . 743.4 R(x0;x0) for the ellipse with N mesh points on the boundary . . 764.1 Comparison of ?2, ?3 and ?n for various values of ", N = 1;2 and 4.117viList of Figures1.1 Structure of a dendritic spine and the trajectory of a receptor . . 121.2 A circular con?nement domain . . . . . . . . . . . . . . . . . . . 133.1 MFPT for the unit disk as the initial position moves from (?1;0)towards the absorbing arc at (1;0) with D = 1, d = 12 and " = 10?5. 453.2 MFPT from the centre of the unit disk with two holes on theboundary with D = 1, d = 12 and " = 10?5. . . . . . . . . . . . . 473.3 MFPT, v(x), from the centre of the unit disk with N symmet-rically located holes on the boundary with D = 1, d = 12 and" = 10?5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 513.4 MFPT for the unit square as the initial position moves from(0;0:5) towards the exit window at (1;0:5) with D = 1, d = 12and " = 10?5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 593.5 MFPT from the centre of the unit square as the exit windowmoves from (1;0:0125) to (1;0:9875) with D = 1, d = 12 and" = 10?5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 613.6 MFPT from the centre of a unit square with two exit windows,one ?xed at (0,0.5) and the other moving from (0,0.48) to (1,0.5),sides 1 and 2 with D = 1, d = 12 and " = 10?5. . . . . . . . . . . 633.7 MFPT from the centre of a unit square with two exit windows,one ?xed at (0,0.5) and the other moving from (0,0.48) to (1,0.5),side 3 and all sides with D = 1, d = 12 and " = 10?5. . . . . . . 643.8 Plot of the Neumann Green's function for the ellipse, Gm ? C,versus ?0 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 753.9 Plot of the MFPT, v(x), from the centre of the ellipse versus ?0with D = 1, d = 12 and " = 10?5 . . . . . . . . . . . . . . . . . . 783.10 Plot of the MFPT, v(x), from x = (1;0) for the ellipse versus ?0with D = 1, d = 12 and " = 10?5 . . . . . . . . . . . . . . . . . . 79viiList of Figures3.11 Plot of the perturbed unit disk, the curvature and ? with " = 0:1and a = 0:2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 883.12 Plot of Rm(x0;x0) for N = 600 (red) and N = 2400 (blue) with" = 0:1 and a = 0:2. . . . . . . . . . . . . . . . . . . . . . . . . . 89viiiAcknowledgementsI would like to thank my supervisors, Michael Ward and Anthony Peirce, forintroducing me to the ?eld of asymptotics and allowing me to work on thisproject. I appreciate all the help they have provided.I would like to thank Sydney Pachmann, a Mathematics graduate student,for her support during my degree. To all the graduate students, thank you formaking my time at UBC unforgettable. I also express my sincere appreciationto all the helpful sta? of the Mathematics Department and the IAM.Financial support was provided by my supervisors and a University GraduateFellowship.ixDedicationTo my family and Ryan.xChapter 1IntroductionBrownian motion describes the perpetual irregular motion of small particles,such as the random motion of smoke particles. This phenomenon was ?rststudied by the British botanist Robert Brown, who noticed the chaotic mo-tions of small grains of pollen, immersed in water under a microscope in 1827.Though Brown was never able to explain what he observed, he was able to dis-pel the notion that the random movements were exclusive to pollen particles byobserving the similar behaviour of dust particles. A mathematical descriptionof Brownian motion was ?rst formulated by Thorvald N. Thiele in 1800 in hispaper on the method of least squares. This mathematical formalism was com-pleted independently by Albert Einstein in 1905 and Marian Smoluchowski in1906.Initially, the term Brownian motion was reserved for the description of therandom movement of particles immersed in a liquid or gas. However, sinceBrown's discovery, the study of Brownian motion has been extensive, ?ndingapplications in various ?elds such as ?nance, biology, and statistical mechanics.This thesis is devoted to a speci?c aspect of Brownian motion, the mean ?rstpassage time, abbreviated as MFPT. We employ both analytic and numericaltechniques to solve this problem for the mean ?rst passage time and compareour ?ndings to relevant empirical results found in the literature.1.1 Brownian MotionWe will give a brief mathematical description of Brownian particles immersedin a ?uid in order to introduce certain key aspects of the motion. This is thesimplest physical model of Brownian motion.We assume that the size of the colloidal particles are much larger than themolecules of the surrounding ?uid. With this assumption it is clear that eachcollision has a negligible e?ect, however, the collective e?ect of multiple collisionswith the surrounding ?uid alters the path of the colloidal particle. We expect1Chapter 1. Introductionthat these collisions happen in rapid succession and in very high numbers. Forexample there are 1021 collisions per second for gold particles of radius 50?mimmersed in a ?uid under standard conditions. Since it is the collective e?ectof many collisions that cause an observable e?ect, we must describe the path ofa Brownian particle statistically.In the case under consideration, there are two main forces that act on theBrownian particle. Firstly, by Stokes law, which assumes low Reynolds number,there is a drag force exerted on a spherical colloidal particle by the ?uid. Thedrag force per unit force per unit mass acting on a spherical colloidal particleis ??v, where ? = 6?a?=m. Here, v is the particle's velocity, a is the radius ofthe particle, ? is the coe?cient of dynamic viscosity of the ?uid and m is themass of the particle.The second force acting on the Brownian particle is due to the individualcollisions with the molecules of the surrounding ?uid. These individual collisionsproduce instantaneous changes to the Brownian particle's acceleration. Thesechanges are random both in magnitude and direction. We denote this ?uctuatingforce by f(t) which satis?es the following properties [20]:? f(t) is statistically independent of v(t),? the variations in f(t) are more frequent than the variations in v(t),? the average of f(t) is zero.Newton's equations of motion lead us todv(t)dt = ??v(t) + f(t): (1.1)This equation is also known as Langevin's equation. This is a stochastic di?er-ential equation, which motivates the fact that Brownian motion is among thesimplest continuous-time stochastic processes. Typically, a stochastic di?eren-tial equation describing a stochastic or Brownian process will be driven by twoterms, a locally deterministic velocity or drift term and a volatility or Gaus-sian disturbance term. The Gaussian disturbance term or white noise term inequation (1.1) is the ?uctuating force f(t).The solution of a stochastic di?erential equation determines the transitionprobability density of the random process it describes. In our case, the solutionto equation (1.1) will determine the transition probability density p(v(t);t j v0)of the random process v(t), where v(0) = v0. We denote the probability density2Chapter 1. Introductionfunction from an initial position v0 to v at time t as p(v(t);t j v0). We will showlater that a di?usion equation describes the time evolution of the probabilitydensity function. The transition probability density function, or pdf as it usuallycalled, can be used to determine the probability. For example we can ?nd theprobability that v(t) 2 A given that the v(t) at t = 0 is v0 by integrating overthe pdf. That is,P(v(t) 2 A j v0(0) = v0) =ZAp(v(t);t;v0)dV: (1.2)We assume that v0 is given, and thusp(v(t);t j v0) ! ?(v ? v0) as t ! 0: (1.3)The situation under consideration can also be described by using statisticalmechanics. In fact, we can deduce the statistical properties of f(t) by comparingthe solution to equation (1.1) to known physical laws. From statistical physicswe know that p(v(t);t j v0) must tend towards the Maxwellian density for thetemperature T of the surrounding ?uid independently of v0 as t ! 1. Hence,p(v(t);t j v0) !? m2?kT?3=2exp??mjvj22kT!as t ! 1: (1.4)The conditions on p(v(t);t j v0) impose conditions on f(t). The solution to(1.1), derived upon using an integrating factor, isv(t) = v0e??t +Z t0e??(t?s)f(s)ds: (1.5)Since, for large t, v(t) ? Rt0 e??(t?s)f(s)ds, we can conclude that the integralin this equation must have the same properties as v(t) and satisfy the Gaussianproperties prescribed by the pdf (1.4). It it worth noting that the integral in(1.5) is not a standard Riemann integral. The integral is a stochastic integral.The white noise term f(t) is of locally unbounded variation, it is everywherecontinuous but no-where di?erentiable and thus the integral cannot be de?nedin straightforward manner. In fact, a new framework has to be used to dealwith this integral, namely the It^ calculus.3Chapter 1. IntroductionWe can describe the displacement of the Brownian particle, x(t), byx(t) = x0 +Z t0v(s)ds: (1.6)This is also a stochastic process driven by a locally deterministic velocity alongwith a Gaussian disturbance term. We state the probability density of x(t) forlarge t,p(x(t);t j x(0) = x0) ? (4?Dt)?3=2 exp?? jx ? x0j24Dt!; D = kTm? = kT6?a?:(1.7)The interested reader is directed to [20] for further details. It can be shown thatp satis?es the di?usion equation@p(x(t);t j x(0) = x0)@t = D?xp(x(t);t j x(0) = x0); (1.8)where D is the di?usion coe?cient. We have found the pdf for the process x.In addition, we can state the following properties of the process x [20]:? The increments x(t+s)?x(t) and x(t)?x(t?u) are independent of eachother and are independent of t for s ? 0 and u ? 0.? The paths of x(t) are continuous.? For s < t, x(t) ? x(s) has a Gaussian distribution with mean zero andvariance t ? s.These are properties of a standard stochastic process known as a Wiener process.It is equation (1.8), which describes the time evolution of the probability densityfunction for the process x that is of particular importance to this thesis.1.2 The Narrow Escape ProblemWe have described certain key aspects of Brownian motion. We now state thespeci?c research problem considered in this thesis.1.2.1 Statement of the ProblemWe consider the motion of a Brownian particle trapped in an arbitrary boundeddomain,   2 Rd d = 2;3, whose boundary is re?ecting, @ r, except for a small4Chapter 1. Introductionabsorbing window, @ a. We assume that @  = @ a+@ r is a d?1 dimensionalanalytic surface and that @  is su?ciently smooth. Furthermore, we assumethat the size of the absorbing window, centred at x0, is small in comparison tothe re?ecting portion of the boundary. We de?ne the small parameter " as" = j@ ajj@ rj? 1: (1.9)Alternatively, we may de?ne " = j@ aj. However, it must be understood thatthe absorbing part of the boundary is asymptotically small in comparison tothe re?ecting part.The trajectory of the Brownian particle is denoted by x(t). The mean ?rstpassage time or exit time is de?ned as the time taken for the Brownian particleto escape from the domain,  , through the absorbing arc, @ a, centred at x0,from some initial position x(0) 2  . The mean ?rst passage time, from a ?xedstarting position x(0) = x is de?ned asv(x) = E [? j x(0) = x] : (1.10)The focus of this thesis is to ?nd the mean ?rst passage time for a variety of twoand three dimensional domains,  , with N absorbing windows on the bound-ary. This is known as the narrow escape problem. Since the mean escape timediverges as the window shrinks, or as " ! 0, the calculation is a singular per-turbation problem. A plethora of work has already been done on this problem.In this thesis we use an alternative method based on matched asymptotic ex-pansions, a singular perturbation technique, to calculate the mean ?rst passagetime.1.2.2 Derivation of the Model EquationThe narrow escape problem is equivalent to solving an inhomogeneous mixedNeumann-Dirichlet boundary value problem for the Poisson equation. We willderive this equation. Our starting point is equation (1.8), which is a di?u-sion equation describing the time evolution of the probability density functionassociated with the process y. The probability density function satis?es the5Chapter 1. IntroductionFokker-Planck equation@p(y;t j x)@t = D?yp; y;x 2  ;@p(y;t;j x)@n(y) = 0; y 2 @ r;x 2  ;p(y;t;j x) = 0; y 2 @ a;x 2  ;p(y;0 j x) = ?(y ? x): (1.11)where @n = ^?r and ^ is the unit outward normal. The boundary conditions areeasily explainable in the context of the physical problem at hand. The Neumannboundary condition on the re?ecting part of the boundary represents the no ?uxboundary condition. The Dirichlet boundary condition on the absorbing partof the boundary represents the fact that the Brownian particle is absorbed atthe boundary. The ?rst passage time to the absorbing boundary is de?ned as? = infft > 0 : y(t) 2 @ ag: (1.12)We are interested in the mean ?rst passage time to @ a given that theparticle starts at some initial position x, that is x(0) = x. We de?ne thesurvival probability asS(x;t) =Z p(y;t j x)dy: (1.13)The nth moments of the ?rst passage time are given by?x = ?Z 10tn @@tS(x;t)dt (1.14)= ?tnS(x;t) j10 +nZ 10tn?1S(x;t)dt: (1.15)The ?rst line is obtained using integration by parts. In the second line, it isimportant to note that we can set the ?rst term to zero in both the upper andlower limit. This is a consequence of the fact that the pdf p(y;t j x) given by(1.11) tends to zero as t ! 1. Thus, in the limit as t ! 1, S(x;t) ! 0 muchfaster than tn ! 1 [18]. We are interested in the mean ?rst passage time,which is the ?rst moment of the equation above. Thus, we ?nd that the mean?rst passage time is the integral of the survival probability over t. Thus, the6Chapter 1. Introductionmean ?rst passage time satis?es the conditional expectation?x = E [? j y(0) = x] =Z 10Z p(y;t j x)dydt: (1.16)We return to the pdf, (1.11). We integrate this equation over t from 0 to 1and make use of the boundary conditions to ?nd??(y ? x) = D?yZ 10p(y;t j x)dt: (1.17)We let g0(y j x) = R10 p(y;t j x)dt. Using this we ?nd that?yg0(y j x) = ??(y ? x)D ; y;x 2  ; (1.18)@g0(y j x)@n(y) = 0; y 2 @ r;x 2  ; (1.19)g0(y j x) = 0; y 2 @ a;x 2  : (1.20)Thus, we can ?nd the MFPT by solving for g0 from the equation above andintegrating the ?nal solution over  . That is?x =Z g0(y j x)dy =Z 10Z p(y;t j x)dydt: (1.21)There is a simpler way to ?nd the MFPT. The pde given by (1.11) is theFokker-Planck equation. That means that the di?erentiation with respect to thespatial variable and the boundary and the initial conditions are given in termsof the forward variable, x, not the initial position, y. One can state a similarequation for the pdf in terms of backward variables, known as the Kolmogorovbackward equation@p(y;t j x)@t = ?D?xp; y;x 2  ; (1.22)@p(y;t j x)@n(x) = 0; x 2 @ r;y 2  ; (1.23)p(y;t j x) = 0; x 2 @ a;y 2  ; (1.24)p(y;t j x) = ?(y ? x): (1.25)There are several di?erences between the Kolmogorov backward equation andthe Fokker-Planck equation. Firstly, the spatial operator in the Fokker-Planck7Chapter 1. Introductionequation is the adjoint of the operator in the Kolmogorov equation, which actson the forward variable, y. In the case under consideration the operator is theLaplacian, which is a self-adjoint operator. However, there is a sign di?erencein the coe?cient of the Laplacian. In the Fokker-Planck equation, (1.11), itis a +1 while in the Kolmogorov equation, (1.22), it is a ?1. Furthermore,if one compares the boundary conditions (1.23) and (1.24) to the boundaryconditions of (1.11), we see that the initial position x starts on the boundaryin the Kolmogorov backward equation, while the ?nal position y is within  .Lastly, there is no initial condition but a ?nal condition in the Kolmogorovbackward equation. Thus, in the limit as t ! 1, the pdf p(y;t j x) tends to?(y ? x) and the pdf p(y;t j x) ! 0 as t ! 0. This is also di?erent to the pdfgiven by (1.11)We know the de?nition of the MFPT, given by (1.16). Thus, we mustintegrate (1.22) over   with respect to y and from 0 to 1 with respect to t.We let v(x) = R10 R  p(y;t j x)dydt. We ?nd that v(x) = ?x satis?es?v(x) = ? 1D; x 2  ; (1.26)v(x) = 0; x 2 @ a; (1.27)@v(x)@n(x) = 0; x 2 @ r: (1.28)The reader is directed to [20] for further reading and an alternate derivation ofthe MFPT, v(x) satisfying (1.26)-(1.28).We also de?ne the average MFPT, ?, based on an assumed uniform distri-bution of starting points,? = 1j jZ v(x)dx: (1.29)1.3 Literature ReviewThe mean ?rst passage time is a ?rst-passage probability, which is the proba-bility that a di?using particle or random walker ?rst hits a speci?ed site. The?rst instances of such a study came with Lord Rayleigh, who considered the ?uxthrough a small hole in the context of acoustics [17]. Today, there are manyapplications of ?rst-passage probabilities, for example ?uorescence quenching,integrate-and-?re neurons and the execution of buy/sell orders when a stockprice reaches a certain threshold [18]. Moreover, the narrow escape problem8Chapter 1. Introductionoften appears in the context of electrostatics and biology.As a result of the far reaching applications of the MFPT, the approachto solving the problem has been varied. Statistical, numerical and analyticaltechniques have all been employed to solve the narrow escape problem. Wepresent a brief overview of a few subject areas where the narrow escape problemarises.As we have shown, the narrow escape problem reduces to a mixed boundaryvalue problem for the Poisson equation. This scenario comes up in a varietyof contexts. In particular these problems arise in classical electrostatics, forexample, the electri?ed disk problem [6]. Furthermore, the problem arises inelasticity theory, di?usion and conductance theory and acoustics [17]. Theseproblems are generally solved using a separation of variables technique, with asubsequent summing of the resulting series. This approach to obtain analyti-cal results has proved to be challenging, especially when considering complexdomains, and ultimately, numerical techniques are employed.The ?rst passage probability and mean ?rst passage time play fundamentalroles in electrostatics as a result of their physical interpretation. In fact, theprobability of a particle, initially at a position x(0), escaping at the point x0 onthe boundary is equal to the electric ?eld at the point x0, when a point charge islocated at x(0) and the boundaries are grounded conductors [18]. This analogyis powerful, and can be applied to calculations in many areas. For example, onecan determine the `break-even' probabilities in the stock market more easily byconsidering this analogy.One can obtain digitized representations of composite materials via Brow-nian motion. The ?rst passage equations are adapted to digitized media andare solved numerically. The ?rst passage time plays a role in determining theconductivity, dielectric constant, and di?usion coe?cient of digitized compositemedia [32].In biology, the motion of ions, molecules or receptors are modelled as freeBrownian particles. In this context, the mean ?rst passage time represents themean time for a receptor to hit a target binding site [25], [23], [21], [24], [5],[2], [1]. The applications are far reaching including, receptor tra?cking in asynaptic membrane, calcium di?usion in dendritic spines and vesicle tra?ckingin cells to name a few.Bresslo? et. al. in [2] calculated the MFPT in a rectangular planar domainwith absorbing portions within the domain. The statement of the problem isslightly di?erent to equation (1.11) in that the pdf has periodic boundary condi-9Chapter 1. Introductiontions in addition to the mixed boundary conditions. The method of solution isa combination of statistical properties and an asymptotic approach. This modeldescribes protein receptor tra?cking within the membrane of a cylindrical den-drite.If one transforms the problem for the MFPT into an eigenvalue problem we?nd that the MFPT is inversely proportional to the principal eigenvalue, in thelimit of small windows. We shall explore this connection in detail in followingsections. In fact, it is the key idea used in this thesis to estimate the MFPT.The problem for the principal eigenvalue and eigenfunction has been extensivelystudied, for example in [8], and the references therein. The Neumann Green'sfunction for the domain under consideration is required to solve this problem.As a result, the form of the mean ?rst passage time is di?erent in two and threedimensions because of the di?erence in the singular behaviour in the NeumannGreen's function. In two dimensions it is a logarithmic singularity and in threedimensions it is a simple pole. We will devote a portion of this thesis to ?ndingthe Neumann Green's function in various domains.There are other closely related problems to the narrow escape problem, suchas Kolomogorov's exit problem and the problem for the Dwell time, which willbe discussed later. However, an extensive listing of all instances where thenarrow escape problem arises is not relevant to this thesis. Instead, we reviewthe relevant work of a few authors in the area who have produced results forescape from two and three-dimensional domains. In this thesis we intend toexpand upon these results, and compare our ?ndings with the work of theseauthors.1.4 Literature Review - Biological ContextWe review the work of certain authors who have determined results for theMFPT in the context of biology. We intend to compare our ?ndings to theresults of these authors.The calculation of the mean ?rst passage time is a typical problem in cellularbiochemistry. In particular, the function of neurobiological microstructures,such as dendrites, is largely unknown since experimental methods have failedto produce a complete understanding of the neuron. An understanding of themechanism by which newly synthesized proteins from the soma reach distantlocations on axons or dendrites is still elusive. Hence, a mathematical framework10Chapter 1. Introductionhas been developed to shed light where experimental methods have failed.1.4.1 Two DimensionsWe consider an example to emphasize the importance of the two-dimensionalnarrow escape problem in a neurobiological context as outlined in [5].The axon contains presynaptic active zones for neurotransmitter release.Dendrites, on the other hand, contain postsynaptic densities where receptorsthat bind neurotransmitters cluster. At most excitatory synapses in the brain, itis the dendritic spines that contain the postsynaptic densities (PSD). Dendriticspines are tiny membranous protrusions, on the order of micrometers, on thesurface of a dendrite. Generally, the spine head is bulbous in shape and isconnected to the stalk of the dendrite via a thin neck. Holcman and Schussin [5] consider the motion of a receptor inserted into a dendritic spine. Thephysical question they address is, how long does it take the receptor to reach its?nal destination on the postsynaptic density (PSD) from its point of insertion?The motion of the di?using receptor is considered to be free Brownian motionin the plane.Firstly, we describe the geometry of the surface of the dendritic spine asin [5]. We consider the surface of the domain to be planar, thus neglectingthe curvature. The surface contains many con?nement domains known as cor-rals. These corrals are smooth two-dimensional domains,  , whose boundaryis re?ecting, @ r, except for a small absorbing portion, @ a. The re?ectingboundary may represent a physical barrier or a potential barrier. We expectthat the random walker, the receptor in this case, gets trapped in these con-?nement domains during its journey towards its ?nal anchoring position on thePSD. Holcman and Schuss aim to ?nd the mean con?nement time, the time areceptor spends in a corral. Figure 1.1 depicts the geometry of the dendriticspine.11Chapter 1. IntroductionFigure 1.1: Structure of a dendritic spine and the trajectory of a receptorThe placement of the corrals, or con?nement domains, are also illustratedin Figure 1.1. The con?nement domain may be an arbitrary two-dimensionalshape. Figure 1.2 depicts a circular con?nement domain along with the trajec-tory of a random walker trapped within it.12Chapter 1. IntroductionFigure 1.2: A circular con?nement domainTypically, the trajectory of the Brownian particle occupies a much largerarea within the domain than what is depicted in Figure 1.2. Figures 1.1 and1.2 are based respectively on Figures 1 and 2 in [5].Physically, the con?nement domain represents sites where the di?using re-ceptor can bind to sca?olding proteins. Thus, the mean time a protein spendsin a con?nement domain provides information on the type of bonds a receptorcan make with sca?olding proteins. Binding increases the mean ?rst passagetime. In fact, a di?using receptor may bind to a sca?olding protein and neverleave the con?nement domain. This is in alignment with the fact that a receptorinserted into a dendritic spine far from the PSD may remain in the spine with-out reaching the PSD. Holcman and Schuss, in [5], are concerned with the meantime it takes a receptor to leave a con?nement domain. In other words, they didnot consider the possibility that receptors may be anchored in the con?nementdomain.13Chapter 1. IntroductionHolcman and Schuss considered a circular con?nement domain, depicted inFigure 1.2. The narrow escape problem for a circular domain with an absorbingarc on the boundary was studied in both [23] and [5].In [23] and [5], the mean ?rst passage time from the centre of a disk of radiusR with a small absorbing arc of length 2" on the boundary was shown to beE [? j x(0) = 0] = R2D?log 1" + log 2 + 14 + O(")?: (1.30)The mean ?rst passage time averaged with respect to an initial uniform distri-bution of starting points in the disk, from [23], isE? = R2D?log 1" + log 2 + 18 + O(")?: (1.31)The maximum mean ?rst passage time is attained when the initial position isantipodal to the centre of the absorbing arc. Singer et al., in [25], placed theinitial position at x = (r = 1;? = 0) and the centre of the absorbing arc atx0 = (r = 1;? = ?). For this situation Singer et al. calculated thatE [? j x = (r = 1;? = 0)] = R2D?log 1" + 2log 2 + O(")?: (1.32)Holcman and Schuss in [5] were able to calculate the mean time to anchoring inthe PSD using the con?nement time. This process underlies synaptic plasticity,the change in the number of receptors at a synapse. Synaptic plasticity isconjectured to underlie processes as complex as learning and memory.The results above conform to a general result derived by Singer et al. in[25], in which they showed that the leading order term and error estimate for atwo-dimensional Riemannian manifold with metric g isE? = j jg?D?log 1" + O(1)?; (1.33)where j jg is the Riemannian area of the domain with the assumption that" = j@ ajgj jg? 1. We see that the leading order term here is the same as thatobtained for the circle of radius R. Notice here, that to this order in ", that themean escape time E? is not a function of the initial position, x(0). Hence, it isdenoted E?, not E[? j x(0) = x].In [24], Singer et al. consider two-dimensional domains with cusps and cor-14Chapter 1. Introductionners. In the other words, they modify their domains to include non-smoothboundaries. If the absorbing window is located at a corner of angle ?, the meanlifetime isE? = j jg?D?log 1" + O(1)?: (1.34)Singer et al., in [24], ?nd that the MFPT in a rectangle with sides a and b withthe absorbing hole of size " at the corner isE? = 2j jD??log a" + log 2? + ?6 ba + 2?2 + O?"a;?4??; (1.35)where ? = e??b=a.They also consider an annulus and the domain bounded between two tangentcircles in [24]. The annulus has applications to a closely related problem for theDwell time, described in [30]. The Dwell time is the time a Brownian particlespends within a microdomain, a con?nement domain, including binding time,before it escapes through @ a. The con?nement domain is a circular annulusof outer radius R and inner radius ?. The boundary of the outer disk, D(R) isre?ecting except for a small absorbing portion, @ a. The boundary of the innerdisk, D(?) is absorbing. A Brownian particle is free to di?use in D(R) ? D(?).If it is absorbed at D(?), it is eventually released within the annulus. If it isabsorbed at @ a, it does not return to the annulus. The inner disk representsan immobile trap, a domain of chemical reactions, where it reacts with bindingmolecules.The results for a two-dimensional domain, given by (1.30)-(1.35), were de-rived using a separation of variables technique and the solution of certain dualintegral equations. We intend to reproduce the results obtained for the circu-lar disk with radius one. Furthermore, we intend to determine an asymptoticexpression for the MFPT in an arbitrary two-dimensional domain.1.4.2 The Narrow Escape Problem in Three DimensionsThe narrow escape problem has been studied for an arbitrary three-dimensionaldomain with an arbitrary exit window in [25]. The leading order term for theMFPT is expressed in terms of an integral equation, which can be solved forspeci?c geometries. Singer et al. were able to explicitly determine the leadingorder term for an absorbing ellipse on the boundary. The additional assumptionwas that the semi-major axis of the elliptic window, a, was to be much smaller15Chapter 1. Introductionthan the cube root of the volume, j j1=3. The MFPT was found to beE? ? j j2?DaK(e); (1.36)where e is the eccentricity of the elliptic hole, and K(:) is the complete ellipticintegral of the ?rst kind. For a circular hole, the result above reduces toE? ? j j4aD: (1.37)These results were obtained by solving for the MFPT using the properties of theNeumann Green's function in three-dimensions. It is clear from these resultsthat the MFPT depends on the geometry of the exit window, as well as on thedomain  . These results do not contain any error estimates, and since these areonly leading order approximations, the dependence of the MFPT on the initialposition and location of the absorbing window is lost. These results were knownto Lord Rayleigh in the context of acoustics [17]. For a spherical domain ofradius R with a circular hole of radius a on the boundary Singer et al., in [25],derived the second term and error estimate of the MFPT,E? = j j4aD?1 + "log 1" + O (")?: (1.38)This result was derived using Collins' method for solving dual series of integralequations and a subsequent expansion of these results in orders of " = aR ? 1.We intend to expand on the results for the spherical domain with a circular exitwindow by deriving the next term in the asymptotic expansion for the MFPT,and to extend the results to include N small windows on the boundary of thesphere.1.5 OutlineIn Chapter 2 we consider the narrow escape problem in two dimensions. Wederive an asymptotic expression for the MFPT in an arbitrary two-dimensionaldomain with one absorbing arc on the boundary. We extend this frameworkto include N absorbing arcs on the boundary. We use a matched asymptoticexpansions approach, a singular perturbation technique, to ?nd the MFPT.In Chapter 3, we use the ?ndings from Chapter 2 to obtain speci?c resultsfor the unit disk and unit square. We also consider methods to solve for the16Chapter 1. IntroductionNeumann Green's function in arbitrary two-dimensional domains. We considerboth an analytical approach and a numerical approach. In Chapter 4 we solvethe narrow escape problem for a three-dimensional sphere. This involves ?ndingthe Neumann Green's function for a sphere with a singularity on the boundaryof the sphere. Wherever possible, we compare our results to relevant ?ndingsin the literature.17Chapter 2The Narrow EscapeProblem in TwoDimensionsIn this chapter, we consider the motion of a Brownian particle trapped in anarbitrary bounded two-dimensional domain,  , whose boundary is re?ecting,@ r, except for a small absorbing window, @ a, through which the particlecan escape. We use the method of matched asymptotic expansions to calculatethe mean escape time. The concepts and method of solution introduced inthis chapter are important as we use similar techniques in Chapter 4, wherewe consider three-dimensional domains. Initially, we begin our analysis withan arbitrary two-dimensional domain with one absorbing arc or hole on theboundary. We then extend our results to include N absorbing arcs of varyinglength. In addition, we derive results for N identical holes.To recap, the length of the absorbing arc is asymptotically small in compar-ison to the length of the re?ecting boundary. In two-dimensions we de?ne " asl" = j@ aj ? 1; (2.1)where l is O(1).The MFPT v(x), from a ?xed starting position x is given by the probabilisticformula [5], [23] and [25]v(x) = E[? j x(0) = x]: (2.2)As the absorbing window shrinks, the mean free path time diverges. We expectthe behaviour v(x) ! 1 as " ! 0 and v(x) ! 0 as x ! x0. Note that x(0) = xis the initial position of the Brownian particle and x0 denotes the centre of theabsorbing arc. These are not to be confused. The MFPT v(x) satis?es the18Chapter 2. The Narrow Escape Problem in Two Dimensionsmixed boundary value problem [5], [23] and [25]?v(x) = ? 1D; x 2  ; (2.3)v(x) = 0; x 2 @ a; (2.4)@v(x)@n(x) = 0; x 2 @ r; (2.5)where D is the di?usion coe?cient. Our task is to calculate an asymptoticsolution for v(x) in the limit " ! 0 that satis?es (2.3)-(2.5).2.1 Formulation of the ProblemWe know from the literature, [23] and [8], and we will show that the mean ?rstpassage time is asymptotically inversely proportional to the principal eigenvalue.Motivated by this fact, we transform our problem into an eigenvalue problemfor the principal eigenvalue.We solve equations (2.3) to (2.5) by writing v(x) in terms of a complete setof eigenfunctions ?j(x), that isv(x) =1Xj=0cj?j(x); (2.6)where ?j(x) satis?es the eigenvalue problem??j + ?j?j = 0; x 2  ; (2.7)?j(x) = 0; x 2 @ a; (2.8)@?j(x)@n(x) = 0; x 2 @ r: (2.9)Recall that the eigenfunctions ?j satisfy the following properties(?j(x);?i(x)) =? 0 if i 6= j;1 if i = j; (2.10)where (u;v) denotes R  uvdx.To ?nd cj in (2.6) we use the orthogonality properties of the eigenfunctions?j. We substitute the expression for v(x) from (2.6) into (2.3) and use (2.7) to19Chapter 2. The Narrow Escape Problem in Two Dimensionsobtain1D =1Xj=0cj?j?j(x): (2.11)By multiplying both sides of (2.11) by ?i and integrating over the domain, andusing (2.10), we ?ndcj = (1;?j(x))D?j(?j(x);?j(x)): (2.12)Substituting the result (2.12) for cj into (2.6) we ?ndv(x) =1Xj=0(1;?j(x))D?j(?j(x);?j(x))?j(x): (2.13)where (?j(x);?j(x)) = 1.Note that we have an expression for the mean escape time, v(x), in terms ofan in?nite series, which is not simple to work with. Luckily, we can make someapproximations that simplify the expression for the escape time. Notice thatas " ! 0, or in other words as @ a ! x0, where x0 denotes the centre of theabsorbing arc, the eigenvalue problem given by equations (2.7) to (2.9) reducesto the unperturbed problem??j + ?j?j = 0; x 2  ; (2.14)@?j(x)@n(x) = 0; x 2 @ : (2.15)We know that the ?rst eigenvalue and eigenfunction for this problem are ?0 = 0and ?0 = j j?12 . Thus, in the limit as " ! 0, we know ?0 ! 0, while the other?j for j ? 1 remain bounded. We can expect ?0 ? O?1log "?as " ! 0, as willbecome more apparent in the calculations that follow. This implies that 1?0 ! 1and that 1?j << 1?0 for j ? 1 as " ! 0. Furthermore, by integrating (2.7)over the domain and using the Divergence Theorem along with the boundaryconditions (2.8) and (2.9) we ?ndZ ??jdx = ??jZ ?jdx;=Z@ r@n?jdS +Z@ a@n?jdS;=Z@ a@n?jdS: (2.16)20Chapter 2. The Narrow Escape Problem in Two DimensionsWe know ?j ? O (1) as " ! 0 for j ? 1 and we can assume that j@ aj ? O (").By considering the ?rst and third line of (2.16), we can deduce that R  ?j ? O (")for j ? 1. Thus, we can see that the most signi?cant contribution to (2.13) comesfrom the ?rst term, j = 0.Thus to leading orderv(x) ? (1;?0(x))D?0?0(x); (2.17)with (?0(x);?0(x)) = 1. Therefore, the expected lifetime or exit time of aBrownian particle is proportional to 1?0(") in accordance with theory, [8] and[23].Our task now is to ?nd ?0 and ?0 using strong localized perturbation theory.The absorbing arc @ a provides a perturbation that is large in magnitudebut small in extent, hence it is known as a strong localized perturbation. Gen-erally strong localized perturbations produce large changes to the solution in alocalized region. We use the method of matched asymptotic expansions to con-struct the solution in a similar way as outlined in [8], [33] and [34]. In particularthere will be expansions for the region in the vicinity of the absorbing arc, knownas the inner region, and expansions in the outer region, the region away fromthe hole. By matching these expansions we can accurately describe the largelocal changes that occur and the relatively small changes to the solution thatoccur away from the absorbing arc. This is how we will proceed.2.2 Two-Dimensional Domain with One Holeon the BoundaryWe start our analysis with one hole or absorbing arc on the boundary. In thenext section we will extend this framework to include N holes on the boundaryof varying length. We must ?rst calculate the principal eigenvalue ?0 and theprincipal eigenfunction ?0 before we can calculate the mean escape time v(x)given by equation (2.17).2.2.1 Calculation of ?0 and ?0We need to solve the system of (2.7)-(2.9) for the principal eigenvalue ?0(") andthe corresponding eigenfunction ?0. We start by expanding the eigenvalue?0(") = ?00 + ?(")?01 + (?("))2?02 + ? ? ? ; (2.18)21Chapter 2. The Narrow Escape Problem in Two Dimensionswith the condition that ?(") ! 0 as " ! 0, where ?(") = ? 1log("d). Here d isthe logarithmic capacitance, a constant, which is dependent on the length ofthe perturbing arc. This expansion was used for a similar problem, where theholes were located within the domain instead of on the boundary. It was shownby [34], that ?01 is independent of the position of the hole, x0. Though, theproblem studied in [34] is a slight variation of the problem at hand, the resultis analogous to our problem. Thus, we must go to higher orders to ?nd termsthat do depend on the location of the hole. This is important so that we obtaina formula for the MFPT that depends on the location of the hole as well ason the initial position. We can write ?0(") in terms of an in?nite logarithmicexpansion?0(") = ?currency1(?) + O? "log "?where ?(") = ? 1log("d): (2.19)In [33] it was shown how to formulate this expansion. We expand?0(") = ?currency1(?) + ??1 + ? ? ? ; (2.20)where ? ? ?m for m > 0.In the outer region, away from the absorbing arc, we expand the globalsolution as?0(x;") = ucurrency1(x;?) + ?u1(x;?) + ? ? ? : (2.21)Substituting (2.20) and (2.21) into (2.7) and (2.9) we obtain for O(?0) that?ucurrency1 + ?currency1ucurrency1 = 0; x 2   n fx0g;@nucurrency1 = 0; x 2 @ r;Z (ucurrency1)2 dx = 1; (2.22)with some singularity condition as x ! x0 to be found from matching innerand outer solutions. For order O(?1) we obtain?u1 + ?currency1u1 = ??1ucurrency1; x 2   n fx0g;@nu1 = 0; x 2 @ r;Z ucurrency1u1dx = 0: (2.23)The integral conditions on ucurrency1 and u1 in the last line of the system of (2.22) and22Chapter 2. The Narrow Escape Problem in Two Dimensions(2.23) are an implementation of the normalization condition on the eigenfunc-tions.Now in the inner region near the absorbing portion of the boundary @ a welety = (x ? x0)" ; (2.24)so that the arc is rescaled as @ 0 = @ a" ? O(1). In the inner region we letv(y;") = ?0(x0 + "y;"):Then we expandv(y;") = ?(")v0(y) + ? ? ? : (2.25)Substituting (2.25) into (2.7) and (2.8), we ?nd to order O(?(")0) that?yv0 = 0; y = @ 0;v0 = 0; y 2 @ 0: (2.26)We writev0(y) = A(?)vc(y); (2.27)with A(?) ? O(1) as " ! 0 and vc de?ned later.We will need to match v(y;") to ?0(x;") in the limit as x ! x0. The limitx ! x0 is equivalent to " ! 0, which implies that jyj ! 1. We will need anO(1) term at leading order. Notice that v0 is multiplied by a factor of ?(") inthe expansion (2.25). To ensure that we obtain an O(1) term in the expansionfor v(y;") in the limit as jyj ! 1 we would need v0(y) ? log jyj as jyj ! 1.We ?nd that vc(y) satis?es?yvc = 0; y = @ 0;vc = 0; y 2 @ 0;vc ! log jyj as jyj ! 1: (2.28)This ensures that the inner (local) solution has logarithmic behaviour at in?nity.The problem for vc(y), (2.28), has a unique solution and the behaviour at 1 is[8]vc(y) ? log jyj ? log d + p:yjyj2 : (2.29)23Chapter 2. The Narrow Escape Problem in Two Dimensionswhere d, the logarithmic capacitance and p = (p1;p2), the dipole vector, aredetermined from the length of the hole, @ a. In our case, where the length ofthe hole is j@ aj = 2", d = 12. In general, for an arc of length l, d = l4. Thisis obtained by a special solution in terms of elliptic cylindrical coordinates. Wepresent a derivation in Appendix A.The inner and outer solutions must match to obtain a consistent solution.Thus, from matching in the vicinity of the arc, in the limit as x ! x0 or jyj ! 1we must haveucurrency1(x;?) + ?u1(x;?) + ? ? ? = ?(")A(?)"log jyj ? log d + p:yjyj2#+ ? ? ? :Writing the inner variables in terms of outer variables using (2.24)ucurrency1(x;?) + ?u1(x;?) + ? ? ? = ?(")A(?)log jx ? x0j ? ?(")A(?)log("d) + ? ? ? :We recall that ?(") = ? 1log("d). Thus, from matching we ?nd thatucurrency1 ! A + A?(")log jx ? x0j as x ! x0: (2.30)This gives us the missing singularity condition as x ! x0. Thus ucurrency1 must satisfythe system (2.22)-(2.30). From the matching, we see that the next order in theinner expansion is O("?). To ensure matching of the inner and outer solutions? = O("?).To ?nd ucurrency1 and ?currency1 we introduce the Helmholtz Green's function G(x;x0;?currency1)and its regular part R(x;x0;?currency1). This Green's function satis?es?G + ?currency1G = 0; x 2  ;@nG = 0; x 2 @  n fx0g; (2.31)G(x;x0;?currency1) = ? 1? log jx ? x0j + R(x;x0;?currency1): (2.32)Notice that if we did not have an absorbing arc, but instead a hole in theinterior of the domain, so that the singular point x0 is in the interior of  , thebehaviour of the Green's function would beG(x;x0;?currency1) = ? 12? log jx ? x0j + R(x;x0;?currency1):However the singularity of the Green's function at the boundary is twice as large24Chapter 2. The Narrow Escape Problem in Two Dimensionsas it is in the interior of the domain, thus we have a factor of ? 1? in front oflog jx ? x0j.The solution to (2.22) with the correct factor for the logarithmic singularityisucurrency1 = ??A?G(x;x0;?currency1): (2.33)The normalization condition in (2.12) determines A as?2A2?2Z G2(x;x0;?currency1)dx = 1: (2.34)By using (2.32) and (2.33) we calculate ucurrency1 in the limit as x ! x0 asucurrency1(x;") = A? log jx ? x0j ? ?A?R(x0;x0;?currency1): (2.35)Comparing this with (2.30), we ?nd that ucurrency1 has the desired behaviour as x ! x0provided thatR(x0;x0;?currency1) = ? 1?? : (2.36)Equation (2.36) is a transcendental equation for ?currency1. To proceed we expand Gin powers of ?currency1 for ?currency1 << 1 asG(x;x0;?currency1) = 1?currency1 G0(x;x0) + G1(x;x0) + ?currency1G2(x;x0) + ? ? ? : (2.37)We substitute (2.37) into (2.31) to obtain for O ? 1?currency1 ? that?G0 = 0; x 2  ;@nG0 = 0; x 2 @  n fx0g: (2.38)For O ?(?currency1)0? we have?G1 = ?G0; x 2  ;@nG1 = 0; x 2 @  n fx0g;Z G1dx = 0: (2.39)25Chapter 2. The Narrow Escape Problem in Two DimensionsFor O ?(?currency1)j?1? with j = 1;2;3;? ? ??Gj = ?Gj?1; x 2  ;@nGj = 0; x 2 @  n fx0g;Z Gjdx = 0: (2.40)We see that G0 is a constant from (2.38). The condition that the integral of Gjover the domain must be zero for j ? 1 is a consequence of the normalizationcondition. It also ensures that we obtain a unique solution, as the Green'sfunction is only unique up to a constant as a result of the Neumann boundarycondition.We can solve for the G0js recursively. We can solve for G0 from (2.38) and(2.39). We make a semi-circular cut out of radius ?, ? << 1, near x0. Weintegrate (2.39) over   and using the Divergence Theorem and the boundarycondition for G1 we ?nd thatlim"!0Z n ??G1dx = lim"!0Z ?0??@G1@? j?="?"d';= ? lim"!0Z n ?G0dx; (2.41)where ? = jx ? x0j. We know that G1 ? ? 1? log ? as x ! x0, thus @G1@? = ? 1??.Substituting this into (2.41) we ?nd thatG0 = ? 1j j: (2.42)Thus G1 satis?es?G1 = 1j j; x 2  ;@nG1 = 0; x 2 @  n fx0g;Z G1dx = 0: (2.43)We call the Green's function satisfying (2.43) the Neumann Green's functionor modi?ed Green's function Gm(x;x0) with regular part Rm(x;x0). We willrefer to G1 as Gm. We know that as x ! x0Gm(x;x0) = ? 1? log jx ? x0j + Rm(x0;x0): (2.44)26Chapter 2. The Narrow Escape Problem in Two DimensionsUsing our expressions for G0 and G1 given by (2.42) and (2.44) respectively,we can rewrite expansion (2.37) asG(x;x0;?currency1) = ? 1?currency1 j j + Gm(x;x0) + O(?currency1): (2.45)Comparing the behaviour of G(x;x0:?currency1) from (2.35) and the behaviour ofG(x;x0) from (2.44) we see thatR(x;x0;?currency1) = ? 1?currency1 j j + Rm(x;x0) + O(?currency1): (2.46)Substituting the expression for R(x;x0;?currency1), given by (2.46), in the limit asx ! x0 into (2.36), we ?nd that?currency1 = ??j j (1 + ??Rm(x0;x0))+ O(?3): (2.47)We have found ?0 and ?0 to ?rst order in the outer region, from (2.20) and(2.21) we see that?0(?) = ??j j ? ?2?2j j Rm(x0;x0) + O(?3); (2.48)?0(x;?) = ?A?? 1?0 j j ? Gm(x;x0)?+ O(?currency1); (2.49)where Gm(x;x0) satis?es (2.43).It remains to ?nd the constant A = A(?). To ?nd the constant we mustimpose the normalization condition, (?0;?0) = 1, which givesA(?) = 1j j12(1 ? ??Rm(x0;x0)) + O(?2); (2.50)where we have used the fact that R  Gmdx = 0 and the expression for ?0 givenby (2.48). Substituting this into the expression for ?0(x;?), (2.49), and using(2.48) we have the following main result:Proposition 2.1: (One Hole) For " ! 0, the ?rst eigenvalue, ?0, and the?rst eigenfunction, ?0, of the system (2.7)-(2.9), have the two-term asymptotic27Chapter 2. The Narrow Escape Problem in Two Dimensionsbehaviour?0(?) = ??j j ? ?2?2j j Rm(x0;x0) + O(?3);?0(x;?) = j j?12 ? ?? j j?12 Gm(x;x0) + O(?2): (2.51)2.2.2 Calculation of the Mean First Passage TimeUsing (2.17) we can ?nd an expression for v(x) that holds in a general two-dimensional domain with one hole on the boundary. Substituting the leadingorder expressions for ?0 and ?0, given by (2.48) and (2.51) respectively, into(2.17) givesv(x) = j j?D ???1 + ? (Rm(x0;x0) ? Gm(x;x0))? + O(?); (2.52)where ? = ? log("d), with d the logarithmic capacitance. Hence we have thefollowing result:Proposition 2.2: (One Hole) For " ! 0, the mean ?rst passage time,v(x), given by (2.17), has the two-term asymptotic behaviourv(x) = j j?D (? log("d) + ? (Rm(x0;x0) ? Gm(x;x0))) + O(?): (2.53)The average MFPT, ?, is? = j j?D (? log("d) + ?Rm(x0;x0)) + O(?): (2.54)Recall that ? = 1j j R  v(x)dx.This expression for the mean ?rst passage time v(x) holds in a general two-dimensional domain. The modi?ed Green's function remains to be found for thegeometry under consideration. As " ! 0, we see that v(x) ! 1, in accordancewith the behaviour that we had anticipated. Also note that we have used theouter expansion for the eigenfunction to construct the mean ?rst passage time.In other words (2.53) holds in the outer region. To ?nd the behaviour in thelocal region, we must use the local expansion for the eigenfunction v(y). In thisregion, we do obtain the behaviour v(x) ! 0 as " ! 0 for jx ? x0j = O(").28Chapter 2. The Narrow Escape Problem in Two DimensionsWe compare our result to that of Singer et al., [25], in which they obtainedthe one-term expansion for the behaviour of v(x) in a two-dimensional domainwith one hole given byE? = j j?D?log 1" + O(1)?:Our result, (2.53), signi?cantly improves upon theirs. We have a two-termexpansion for the mean ?rst passage time. Consequently, our result depends onthe location of the hole, x0, the initial position, x, as well as on the length ofthe hole, through d, and the shape of the domain, through the Green's function.Note that if j@ aj = "l, then d = l4.2.3 Two-Dimensional Domain with N Holes onthe BoundaryWe can extend the expression for v(x) to include N absorbing arcs. The bulkof the analysis in the previous section remains the same, except that we replacethe single arc with N arcs, @ i for i = 1;2;::;N, of length O("), where xi canbe interpreted as the centre of the arc. Here the length of each arc isj@ ij = "li; (2.55)for some li = O(1). We assume that the arcs are non-overlapping. We restatethe problem for the mean ?rst passage time v(x) with N absorbing arcs on theboundary as?v(x) = ? 1D; x 2  ; (2.56)v(x) = 0; x 2 @ i; i = 1;2;:::;N; (2.57)@v(x)@n(x) = 0; x 2 @ r; (2.58)where @ i denotes the ith absorbing arc on the boundary.We proceed as before, by writing v(x) in terms of a complete set of eigen-29Chapter 2. The Narrow Escape Problem in Two Dimensionsfunctions ?j(x) which satisfy??j + ?j?j = 0; x 2  ; (2.59)?j(x) = 0; x 2 @ i; i = 1;2;:::;N; (2.60)@?j(x)@n(x) = 0; x 2 @ r: (2.61)After imposing the orthogonality and orthonormal properties, we arrive at thesame expression for the MFPT v(x) as in equation (2.17). The di?erence isthat our eigenfunctions satisfy a slightly di?erent problem to that for one hole.We must solve solve the system (2.59)-(2.61) for the leading order eigenfunctionand eigenvalue.2.3.1 Calculation of ?0 and ?0For " ! 0 we expand the eigenvalue ?0 and the eigenfunction ?0(x) as givenby (2.20) and (2.21). The subtlety to notice here, is that we no longer have oneparameter ?, instead, we have N parameters ?i for i = 1;2;:::N where?i = ? 1log("di); (2.62)where we incorporate the logarithmic capacitance, di of each arc. In this mannerwe capture the length of each arc. In particular, if j@ ij = "li, then d = li4 .Now we proceed as in the case for one hole. The problem in the outer regionaway from the hole at order O(?0) reads as?ucurrency1 + ?currency1ucurrency1 = 0; x 2   n fx1;x2;:::;xNg;@nucurrency1 = 0; x 2 @ r: (2.63)Now in the inner region near the absorbing portions of the boundary @ ifor i = 1;:::N we letyi = (x ? xi)" ;so that the arc is rescaled as j@ 0ij = j@ ij" = O(1). In the inner region we letv(yi;") = ?0(xi + "y;"):30Chapter 2. The Narrow Escape Problem in Two Dimensionsand we expandv(yi;") = ?(")v0(yi) + ? ? ? : (2.64)Proceeding as before we arrive at the order O(?(")0) equation for each hole?yv0 = 0; yi = @ 0i;v0 = 0; yi 2 @ 0i: (2.65)According to the same reasoning outlined for one hole, we want v0(yi) ?log jyij as jyij ! 1. This ensures that we can eventually match the inner andouter solutions. We writev0(y) = Ai(?i)vc(yi); (2.66)with Ai(?i) ? O(1) as " ! 0. So we ?nd that vc(yi) satis?es4yvc = 0; yi = @ 0i;vc = 0; y 2 @ 0i;vc ! log jyij as jyij ! 1: (2.67)Recall that if j@ ij = "li, then di = li4 . Matching the inner and outersolutions as before we ?nd that ucurrency1 must satisfy?ucurrency1 + ?currency1ucurrency1 = 0; x 2   n fx1;x2;:::;xNg;@nucurrency1 = 0; x 2 @ r;Z (ucurrency1)2d  = 1; (2.68)ucurrency1 ! Ai + Ai?i(")log jx ? xij as x ! xi: (2.69)This is similar to (2.22) and (2.30) for ucurrency1 for one hole. However, notice thatwe have N unknowns Ai for i = 1;:::;N, with only one normalization condi-tion. Thus, we will have a single relation between each of the Ai's set by thenormalization condition.We write the solution for ucurrency1 in terms of the Helmholtz Green's functionucurrency1 = ???Nk=1Ak?kG(x;xk;?currency1) (2.70)where k is the index of the N absorbing arcs. The Green's function G(x;xk;?currency1)31Chapter 2. The Narrow Escape Problem in Two Dimensionsmust satisfy a system similar to (2.31). The Green's function satis?es?G + ?currency1G = 0; x 2  ;@nG = 0; x 2 @  n fx1;x2;:::;xNg; (2.71)G(x;xk;?currency1) ? ? 1? log jx ? xkj + R(xk;xk;?currency1) as x ! xk:(2.72)Substituting (2.72) into (2.70) will give an expression for ucurrency1 in the limit asx ! xi. This expression must match to (2.69). After some simpli?cation, inthe limit as x ! xi, we require thatAi(1 + ??iR(xi;xi;?currency1)) + ?NXk=1;k6=iAk?kG(xi;xk;?currency1) = 0; i = 1;2;:::;N;(2.73)to satisfy the matching condition to the inner solution. This is di?erent fromthe result for one hole, and is a consequence of the fact that there are N holes.System (2.73) is an NXN homogeneous linear system for the N unknownsAi, i = 1;:::;N, which we write in matrix form asMA = 0;where A = (A1;A2;:::;AN)T is a vector and M is the matrix shown below,whose determinant must be zero in order for nontrivial solution to exist. Thisallows us to solve for the Ai numerically. Here we have denoted R(xi;xi;?currency1)by Rii(?currency1) and G(xi;xj;?currency1) by Gij(?currency1) for i 6= j. The matrix M is0BBBBBBB@1 + ??1R11(?currency1) ??2G12(?currency1) ??3G13(?currency1) ? ? ? ??NG1N(?currency1)??1G21(?currency1) 1 + ??2R22(?currency1) ??3G23(?currency1) ? ? ? ??NG2N(?currency1)??1G31(?currency1) ??2G32(?currency1) 1 + ??3R33(?currency1) ? ? ? ??NG3N(?currency1)... ... ... ... ...??1GN1(?currency1) ??1GN1(?currency1) ? ? ? ? ? ? 1 + ??NRNN(?currency1)1CCCCCCCAWe expand each of the Gij in orders of ?currency1 asG(xi;xj;?currency1) = 1?currency1 G0(xi;xj) + G1(xi;xj) + ?currency1G2(xi;xj) + ? ? ? : (2.74)Consequently, each of the Gij must satisfy equivalent forms of (2.38), (2.39)32Chapter 2. The Narrow Escape Problem in Two Dimensionsand (2.40). Applying the Divergence Theorem as we previously did in the caseof one hole, we ?nd that for each Gij we haveG0(xi;xj) = ? 1j j; (2.75)and thus G1 satis?es?G1(xi;xj) = 1j j; x 2  ;@nG1(xi;xj) = 0; x 2 @  n fx1;x2;:::;xNg; (2.76)Z G1(xi;xj)dx = 0;G1 ? ? 1? log jx ? xjj + R1(xj;xj) as x ! xj:(2.77)This is the problem for the modi?ed Green's function or Neumann's Greenfunction. We shall henceforth refer to G1 as Gm. Now we can rewrite (2.74) for?currency1 << 1 using Gm asG(xi;xj;?currency1) = ? 1?currency1 j j + Gm(xi;xj) + O(?currency1): (2.78)Comparing the behaviour of this expression as xi ! xj with (2.72), wherewe set x = xi, we ?nd thatR(xi;xi;?currency1) = ? 1?currency1 j j + Rm(xi;xi): (2.79)We substitute (2.78) and (2.79) into (2.73) to obtainAi?1 + ??iRm(xi;xi) ? ??ij j ?currency1?+ ?NXk=1;k6=iAk?k?Gm(xi;xk) ? 1?currency1 j j?= 0;(2.80)for i = 1;:::;N.It is convenient to rewrite this in matrix form as an eigenvalue problem.33Chapter 2. The Narrow Escape Problem in Two DimensionsFirstly we do some rearrangingAi(1 + ??iRm(xi;xi)) + ?NXk=1;k6=iAk?kGm(xi;xk)= ?j j ?currency10@Ai?i +NXk=1;k6=iAk?k1A; (2.81)for i = 1;:::;N.We can write this equation in matrix form asCA = ?j j ?currency1 BV A; (2.82)whereC = I + ??V; (2.83)andV =0BBBBBBB@?1 0 ? ? ? ? ? ? 00 ?2 0 ? ? ? 00 0 ?3 ? ? ? 0... ... ... ... ...0 ? ? ? ? ? ? 0 ?N1CCCCCCCA; B =0BBBBBBB@1 1 ? ? ? ? ? ? 11 1 ? ? ? ? ? ? 11 1 ? ? ? ? ? ? 1... ... ... ... ...1 1 ? ? ? ? ? ? 11CCCCCCCA;A =?A1 A2 ? ? ? ? ? ? AN?T;? =0BBBBBBB@Rm(1;1) Gm(1;2) Gm(1;3) ? ? ? Gm(1;N)Gm(2;1) Rm(2;2) Gm(2;3) ? ? ? Gm(2;N)Gm(3;1) Gm(3;2) Rm(3;3) ? ? ? Gm(3;N)... ... ... ... ...Gm(N;1) Gm(N;2) ? ? ? ? ? ? Rm(N;N)1CCCCCCCA:The Green's function is symmetric since Gm(i;j) = Gm(j;i) where Gm(i;j) de-notes Gm(xi;xj) and Rm(i;i) denotes Rm(xi;xi). As a result, ? is a symmetricmatrix. Let ?m = max ?k for k = 1;:::;N. For ?m su?ciently small, we caninvert C so that we obtain a matrix eigenvalue problem for the eigenvalue ?currency1?j jC?1BV A = ?currency1A: (2.84)34Chapter 2. The Narrow Escape Problem in Two DimensionsWe let ? = ?j jC?1BV , this gives us?A = ?currency1A; ? = ?j jC?1BV: (2.85)We can ?nd ?currency1 by ?nding the eigenvalue of the system above. We followthe procedure as outlined in [8]. Notice that each row of BV is the same[?1;?2;:::;?n], thus there is only one linearly independent row. Hence, BV hasrank one. As a result, we can conclude that ? has rank one. Thus, we canconclude that ?currency1 is the only unique non-zero eigenvalue of ?. Hence, ?currency1 =Trace(?). By using the structure of ? as de?ned in (2.85), we ?nd that?currency1 = Trace(?) = ?j jNXj=1?j? NXk=1cjk!; cjk = ?C?1?jk : (2.86)We are left to ?nd C?1. We have assumed that ?m << 1, thus the asymp-totic inverse of C isC?1 = I ? ??V + ? ? ? :Substituting this into equation (2.86) we ?nd that?currency1 = ?j jNXj=1?j? NXk=1Ijk ? ?(?V )jk!; (2.87)where ?V = PNj=1 PNk=1 ?k?jk and PNk=1 Ijk = 1 for each j. Thus?currency1 = ?j j0@NXj=1?j ? ?NXj=1?jNXk=1?k?jk1A + O(?3m): (2.88)To ?nd ucurrency1 we use (2.70) and (2.74). We have found a two-term expansionfor ?0 and ?0. We summarize our result as follows:Proposition 2.3: (N Holes) For " ! 0, the ?rst eigenvalue, ?0, and ?rsteigenfunction, ?0, of the system (2.59)-(2.61), have the two-term asymptotic35Chapter 2. The Narrow Escape Problem in Two Dimensionsexpansions?0(") = ?j j0@NXj=1?j ? ?NXj=1NXk=1?j?k?jk1A + O(?3m); (2.89)?0(x;") ? ??currency1 j jNXk=1Ak?k ? ?NXk=1Ak?kGm(x;xk): (2.90)This gives us a general expression for a two-dimensional domain with N holesof di?ering length on the boundary. It must be noted that we have not explicitlyimposed the normalization condition here to ?nd the relationship between theAi. The Ai can be found by ?nding the eigenvector of the matrix system (2.85).Notice that in the case of N identical holes, there is simpli?cation in theexpression for the eigenvalue ?0 and the eigenvalue ?0. Suppose that we haveN identical holes. Thus, each of the ?i are identical. The expression for ?0given by (2.90) reduces to?0(x;") ? ???0 j jNXk=1Ak ? ??NXk=1AkGm(x;xk): (2.91)Note that PNj=1 PNk=1 ?jk = PNk=1?Rm(xk;xk) + PNj=1;j6=k Gm(xj;xk)?.To simplify the notation, we de?ne p(x1;x2;:::;xN) byp(x1;x2;:::;xN) =NXk=10@Rm(xk;xk) +NXj=1;j6=kGm(xj;xk)1A (2.92)In the case of N identical holes we can impose the normalization conditionexplicitly to ?nd the relationship between the Ai. We calculate (?0;?0) using(2.91) and (2.89) to ?nd thatNXk=1Ak = Nj j12?1 ? ??N p(x1;x2;:::;xN)?+ O(?2); (2.93)where p is given by (2.92).As a Corollary to Proposition 2.3, we obtain the following result for N iden-tical holes:Corollary 2.4 (N Identical Holes) Suppose that the N holes are identical,36Chapter 2. The Narrow Escape Problem in Two Dimensionsthat is ?j = ?. The expressions for ?0 and ?0, given by (2.89) and (2.90),simplify to?0(x;") = ??j j (N ? ??p(x1;x2;:::;xN)) + O(?3); (2.94)?0(x;") = j j?12 ? ???Nk=1AkGm(x;xk) + O(?2): (2.95)For ? ? 1, we observe that the eigenvalue, ?0, given by (2.94) is largest whenthe hole locations x1;x2;:::;xN are chosen so as to minimize p(x1;x2;:::;xN).The mean escape time, which is inversely proportional to the principal eigen-value, will be a minimum in this case.There is in fact one more special case that we can consider. In the case of Nidentical holes, there arises the possibility of attaining a cyclic matrix. Recallfrom the matrix equation (2.85) where we have ?A = ?currency1A, we can ?nd A asthe eigenvector associated with the eigenvalue ?currency1. If we have an even number ofequally spaced identical points in certain geometries, we ?nd that the matrix ?is cyclic. For example, if we have two holes placed at the corners of a rectangle orsquare, we obtain a cyclic matrix. In addition, we ?nd that this property holdsfor any two points located on the circumference of the unit circle. Furthermore,we ?nd that the property holds for four points, located along lines of symmetry.These will be illustrated in following sections. The cause of this special propertyof the matrix ? is that the matrices ? and C and C?1 are cyclic. As a result, we?nd that the eigenvector associated with the eigenvalue ?currency1 has identical entries.Thus, A1 = A2 = ::: = AN. Using this unique property, we can ?nd a simpleexpression for v(x) analogous to the expression we found for one hole in (2.53).Notice that, for identical A0s, ?0 becomes?0(x;") ? ?NA??currency1 j j ? ??ANXk=1Gm(x;xk): (2.96)We can simplify the expression for PNk=1 Ak given by (2.93) to ?ndA = 1j j12?1 ? ??N p(x1;x2;:::;xN)?+ O(?2): (2.97)We can now state the new expression for ?0 in the case of a cyclic matrix37Chapter 2. The Narrow Escape Problem in Two DimensionsCorollary 2.5: (N Identical Holes, ? Cyclic) Suppose that the N holesare identical, that is ?j = ? and the holes are placed in such a manner that thematrix ? in system (2.85) is cyclic. The expression for ?0, given by (2.90),simpli?es to?0(x;") = j j?12 ? ?? j j?12NXk=1Gm(x;xk) + O(?2): (2.98)2.3.2 Calculation of the Mean First Passage TimeNext we use (2.17) to ?nd the MFPT in an arbitrary two-dimensional domainwith N holes on the boundary. We will have three di?erent expressions for theMFPT for the following cases, N holes of di?ering length on the boundary, Nidentical holes on the boundary, and N identical holes where we have a cyclicmatrix.For N holes of di?ering length we substitute the expressions for ?0 and ?0from (2.89) and (2.90) into (2.17). In terms of ?0, we ?nd a two-term expressionfor v(x):Proposition 2.6: (N Holes) For " ! 0, the two-term asymptotic be-haviour of the mean ?rst passage time, v(x), isv(x) ? ?PNk=1 Ak?kD(?0)2???0 j jNXk=1Ak?k ? ?NXk=1Ak?kGm(x;xk)!: (2.99)The average MFPT, ?, is? ? ?PNk=1 Ak?kD(?0)2???0 j jNXk=1Ak?k!: (2.100)Notice that the normalization condition has not been imposed in equation(2.99) for the MFPT. Therefore, we must divide expression (2.99) by (?0;?0),once we have found the Ai numerically by solving the eigenvalue problem (2.85).For N identical holes, we use expressions (2.94) and (2.95) for the eigen-function and eigenvalue. We have already imposed the normalization conditionin this case. We obtain that the mean ?rst passage time for a two-dimensional38Chapter 2. The Narrow Escape Problem in Two Dimensionsdomain with N identical holes on the boundary is given by the following:Corollary 2.7: (N Identical Holes) For " ! 0, the two-term asymptoticbehaviour of the mean ?rst passage time, v(x), isv(x) = j j?ND?? log("d) ? j j12 ?NXk=1AkGm(x;xk) + ?N p(x1;x2;:::;xN)!+O(?);(2.101)where p(x1;x2;:::;xN) is de?ned by (2.92). The average MFPT, ?, is? = j j?ND?? log("d) + ?NNXk=1p(x1;x2;:::;xN)!+ O(?): (2.102)In the case of the matrix ? being cyclic we ?nd:Corollary 2.8: (N Identical Holes, ? Cyclic) For " ! 0, the two-termasymptotic behaviour of the mean ?rst passage time, v(x), isv(x) = j j?ND?? log("d) + ??1N p(x1;x2;:::;xN) ?NXk=1Gm(x;xk)!!+ O(?): (2.103)The average MFPT, ?, is? = j j?ND?? log("d) + ?N p(x1;x2;:::;xN)?+ O(?): (2.104)Notice that by setting N = 1 in equations (2.98) and (2.103), we recover(2.51) and (2.53). Also notice that in (2.102) and (2.104) that the averageMFPT is minimized for a con?guration of arcs that minimize p(x1;x2;:::;xN).Holcman and Schuss, in [5], consider the case of N absorbing arcs of di?eringlength. We cannot directly compare our results with those of Holcman andSchuss, since their results are in terms of an in?nite series. In the case of Nidentical arcs, Holcman and Schuss ?nd that the leading order term for the mean?rst passage time is proportional to 1N . The leading order term in our results,(2.101) and (2.103), are also proportional to 1N . Our results for the mean ?rstpassage time depend on x, xk, for k = 1;:::;N, the length of the arcs through d39Chapter 2. The Narrow Escape Problem in Two Dimensionsand the shape of the domain,  . This is an improvement on the results derivedin [5].2.4 DiscussionIn this chapter we solved for the MFPT in an arbitrary two-dimensional domain.We transformed the problem at hand to an eigenvalue problem. We used themethod of matched asymptotic expansions to solve for the principal eigenvalue,?0, and eigenfunction, ?0. In the case of one hole on the boundary, our expres-sion for the MFPT, (2.53), agrees with the leading order term derived by Singeret al. in [25]. Furthermore, our result is an improvement on the result derivedby Singer et al. since we have the next term in the asymptotic expansion for theMFPT, which depends on the position of the hole at x0 and the initial positionat x. We were able to extend this framework to ?nd general expressions forthe MFPT with N holes on the boundary. We found expressions for the MFPTwith N holes of di?ering size on the boundary, (2.99) and N identical holes onthe boundary, (2.101). Furthermore, we found the MFPT, equation (2.103), inthe special case of a cyclic matrix in the system (2.85).In the next chapter we will apply these results to calculate the MFPT in afew domains.40Chapter 3Numerical Realizations:The Neumann Green'sFunctionIn this chapter, we apply the results derived in Chapter 2 to a few specialdomains. It is clear from our expressions for the MFPT in Chapter 2, that theNeumann Green's function plays an important role. In fact, to ?nd the MFPT ina given domain, one must ?nd the Neumann Green's function. For the unit diskand unit square we can calculate the Neumann Green's function, Gm(x;x0),and its regular part, Rm(x;x0), analytically. For more general domains, wepresent and implement a boundary element method to numerically calculateGm(x;x0) and Rm(x;x0). We use these results to ?nd the mean ?rst passagetime.3.1 The Unit DiskWe have found an expression for the mean ?rst passage time, (2.53), for a generaltwo-dimensional domain with one hole on the boundary. In addition, we havefound three di?ernt formulas for the mean ?rst passage time with N holes on theboundary, namely (2.99), (2.101) and (2.103). We would now like to investigatethe behaviour of the MFPT in speci?c geometries. We begin with the unit disk.In the next section we will consider the unit square.3.1.1 The Neumann Green's Function in a Unit DiskTo investigate the mean free path time in the unit disk, we must ?nd the Neu-mann Green's function Gm for the unit disk. For x0 2  , so that x0 is in the41Chapter 3. Numerical Realizations: The Neumann Green's Functioninterior of  , the modi?ed Green's function satis?es?Gm = 1j j ? ?(x ? x0); x 2  ; (3.1)@nGm = 0; x 2 @ ; (3.2)Z Gm(x;x0)dx = 0; (3.3)where the integral condition ensures that we obtain a unique solution to theproblem.The solution to the problem above for the unit disk is well known, and takesthe form [8]Gm(x;x0) = ? 12? log jx ? x0j? 12? log????xjx0j ? x0jx0j????+ 14?(jxj2 +jx0j2)+C(x)(3.4)where the singularity, x0, is located in the interior of the disk and C(x) is tobe determined from the integral condition. Writing the expression above in theformGm(x;x0) = ? 12? log jx ? x0j + Rm(x;x0);we see that the regular part Rm is given byRm(x;x0) = ? 12? log????x0 jx0j ? x0jx0j???? + 14?(jxj2 + jx0j2) + C(x): (3.5)For the theory in Chapter 2, the singularity is located on the boundary of thecircle, so that jx0j = 1. However, we ?rst ?nd the constant C before we takethe limit as x0 ! @ .To determine C, we multiply (3.4) for Gm(x;x00) by (3.1) and integrate over . We make use of the integral property (3.3) to obtainGm(x0;x00) = ?Z Gm(x;x00)?Gm(x;x0)dx: (3.6)We now use integration by parts and the Neumann boundary condition (3.2) to?ndGm(x0;x00) =Z rGm(x;x00) ? rGm(x;x0)dx: (3.7)We can see from this equation that Gm(x0;x00) = Gm(x00;x0). This is just aproof of the symmetry property of the Green's function, which follows from thefact that the Laplacian is a self-adjoint operator. It does allow us to deduce that42Chapter 3. Numerical Realizations: The Neumann Green's FunctionC(x0) = C(x00) = C. Before we ?nd C, we show another property. Followingthe same argument as outlined above, along with the symmetry property of theGreen's function, we ?nd that1j jZ G(x;x00)dx = Gm(x0;x00) ?Z rGm(x;x00) ? rGm(x;x0)dx: (3.8)Imposing the symmetry property, we see that the R  Gm(x;x0)dx is indepen-dent of x0. We will need this result later on when we consider a numericalapproach to ?nding the Neumann Green's function.To determine C, we integrate (3.4) over  . Since we know that the integralover G is independent of x0, we pick x0 = 0. We know that R  log jxj dx =??=2, R  jxj2 dx = ?=2 and R  Gm(x;0)dx = 0. Using this, we ?nd thatC = ? 38?.Using this, the expression for the Neumann Green's function in a unit diskwith a singularity on the boundary is obtained by letting jx0j ! 1. We obtainGm(x;x0) = ? 1? log jx ? x0j + 14? jxj2 ? 18?: (3.9)Note that this expression is consistent with the remark regarding (2.31) for asurface Green's function, in that a singularity on the boundary is twice as largeas one in the interior.Now we are ready to use expression (2.53) to ?nd the mean free path timefor a Brownian particle trapped in a unit circle with one hole on the boundary.3.1.2 The Mean First Passage Time in a Unit Disk withOne Hole on the BoundaryWe set our initial position to be the origin, x = (0;0) and we place the singu-larity at x0 = (1;0). Note that Rm(x0;x0) = 18?. We substitute for Rm(x0;x0)and Gm(x;x0) with the starting point x at the centre of the disk and the sin-gularity x0 at (1;0) into (2.53) to obtainv(x) = j j?D?? log("d) + 14?+ O(?): (3.10)where d is the logarithmic capacitance. For a length j@ j = 2", d = 12 andj j = ? we ?nd that:43Chapter 3. Numerical Realizations: The Neumann Green's FunctionProposition 3.1: (One Hole in the Unit Disk) In the limit that " ! 0,we ?nd the two-term asymptotic behaviour of the MFPT from the centre of thedisk to beE[? j x(0) = (0;0)] = v(x) = 1D?log 1" + log 2 + 14?+ O(?): (3.11)Notice that with our initial position, x, at the origin, and an x0 locatedanywhere on the boundary, we will always arrive at the same result, (3.11).This is a consequence of the symmetry of the circle, and that all points onthe boundary are located an equal distance r = 1 from the origin. This iscomparable to the result derived by Singer et al., equation (1.2) of [23]Next we ?nd the MFPT with the initial position at the antipodal point tothe absorbing arc, x = (?1;0):Proposition 3.2: (One Hole in the Unit Disk) In the limit that " ! 0,we ?nd the two-term asymptotic behaviour of the MFPT from the antipodal pointto the absorbing arc to beE[? j x(0) = (?1;0)] = v(x) = 1D?log 1" + 2 log 2?+ O(?): (3.12)This result is also comparable to that derived by Singer et al. in equation(1.4) of [23]. Lastly, we investigate the behaviour of moving the initial positionx from the antipodal point at x = (?1;0) towards the singular point, at (1,0),in ?xed increments of 140 in the x direction. We set D = 1, d = 12 and " = 10?5.We expect that the mean free path time decreases as the initial position movestowards the hole. The con?guration and results are shown in Figure 3.1 below.44Chapter 3. Numerical Realizations: The Neumann Green's Function(a) The initial position moves from (?1;0) (blue) towards the hole at (1;0)(red)(b) Plot of the MFPT, v(x), as a function of the x coordinate of the initialpositionFigure 3.1: MFPT for the unit disk as the initial position moves from (?1;0)towards the absorbing arc at (1;0) with D = 1, d = 12 and " = 10?5. 45Chapter 3. Numerical Realizations: The Neumann Green's FunctionIt is important to note that the initial position does not actually reach theabsorbing arc, at (1;0), since our expression for the MFPT is not valid in theinner region. The MFPT does indeed decrease as the initial position movestowards the hole.3.1.3 The Mean First Passage Time in a Unit Disk withTwo Holes on the BoundaryWe will now make use of formula (2.103) to ?nd the MFPT in a unit circle withtwo holes on the boundary. As mentioned, in this case the placement of theholes results in the matrix ? in the system (2.85) being cyclic.To illustrate the behaviour of the MFPT, we hold one absorbing arc ?xedat x0 = (1;0) while we move another hole in ?xed angular increments of ? = ?32from a position an increment away from (1;0), at position 1, to the antipodalpoint at (?1;0), at position 32. We keep the initial position ?xed at the origin.We set D = 1, d = 12 and " = 10?5. The simpli?ed expression for the MFPTfrom the origin isE[? j x(0) = (0;0)] = 12D?log 1" + 38 + 12 log 2 ? 12 log(1 ? cos ?)?+ O(?):(3.13)The con?guration and the behaviour of the MFPT as a function of positions 1to 32 are shown in Figure 3.2.46Chapter 3. Numerical Realizations: The Neumann Green's Function(a) Con?guration of the absorbing arcs (red) and initial position (blue)(b) Plot of the MFPT, v(x), versus the con?guration numberFigure 3.2: MFPT from the centre of the unit disk with two holes on theboundary with D = 1, d = 12 and " = 10?5.47Chapter 3. Numerical Realizations: The Neumann Green's FunctionWe see that the MFPT decreases as the two holes on the boundary movefurther apart. The MFPT reaches a minimum when the two holes are antipodalto each other, which is what we would expect.3.1.4 Equally spaced Points on a Unit DiskWe investigate one more case in the unit disk. By placing N absorbing arcs ofequal length at equally spaced points along the circumference of a unit circle,we can ?nd a closed form expression for the eigenvalue ?0 and the mean freepath time v(x). For a pattern of N identical holes located symmetrically on thecircumference of the unit disk, we havexj = e2?ij=N j = 1;:::;N; (3.14)with N > 1 and i = p?1. We start by simplifyingp(x1;x2;:::;xN) =NXk=10@Rm(xk;xk) +NXj=1;j6=kGm(xj;xk)1A: (3.15)Recall that Gm(xj;xk) = ? 1? log jxj ? xkj + 18? as stated in (3.9). We ?rststate the following lemma [8]:Lemma 3.3: Let N > 0 and n be integers, and i = p?1. Then, for y > 0,we haveNYj=1?x ? ye2?i(j?n)=N?= xN ? yN: (3.16)Proof: The interested reader is directed to the proof in [8].Now, we consider each term of Gm separately and along with the previousLemma we formulate the following Lemma:Lemma 3.4: Let N > 1 be an integer and let xj for j = 1;2;:::;N satisfy(3.14). Then we haveNXj=1?jxjj2 + jxkj2?= 2N; (3.17)NXj=1;j6=klog jxj ? xkj = log N: (3.18)48Chapter 3. Numerical Realizations: The Neumann Green's FunctionProof: The ?rst result, (3.17), is immediate. To prove the second result (3.18),we start withNXj=1;j6=klog jxj ? xkj =NXj=1;j6=klog???e2?ij=N ? e2?ik=N???;= log??????NYj=1;j6=k?1 ? e2?i(j?k)=N???????: (3.19)Now using Lemma 3:3 we ?ndlog??????NYj=1;j6=k?x ? ye2?i(j?k)=N??????? = log????xN ? yNx ? y????= log????xN?1?1 +?yx?+ ? ? ? +?yx?N?1????? (3.20)By using (3.20) with x = y = 1, and substituting this into (3.19), we obtain(3.18).?Using the Lemma above along with the Green's function (3.9) we can rewritep asp =NXk=124NXj=1? 14??jxjj2 + jxkj2?? 38??? 1?NXj=1;j6=klog jxj ? xkj35;=NXk=1?2N4? ?3N8? ?1? log N?: (3.21)Proposition 3.5: (Unit Disk) Let N > 1 be an integer, and let xj satisfy(3.14). Then p = p(N), given by (2.93) reduces top(N) = 1??N28 ? N log N?: (3.22)We can use the simpli?cation above in (2.94) for ?0 to give us:Proposition 3.6: (Unit Disk) For " ! 0, the two-term asymptotic be-49Chapter 3. Numerical Realizations: The Neumann Green's Functionhaviour of the principal eigenvalue, ?0, in a unit disk with N symmetricallylocated holes on the boundary is?0(") = ??j j?N ? ??N28 ? N log N??+ O(?3): (3.23)From (2.101) we determine the mean ?rst passage time by using (3.22) forp. This leads to the following result:Proposition 3.7: (Unit Disk) For " ! 0, the two-term asymptotic be-haviour of the mean ?rst passage time, v(x), in a unit disk with N symmetricallylocated holes on the boundary isv(x) = j j?ND???1 + N8 ? log N ? ?NXk=1Gm(x;xk)!+ O(?): (3.24)We now ?nd an explicit expression for the mean escape time from the centreof the disk with N absorbing arcs each with length 2", for which d = 12:Proposition 3.8: (Unit Disk) For " ! 0, the two-term asymptotic be-haviour of the mean ?rst passage time, v(x), in a unit disk with N symmetricallylocated holes on the boundary, from the centre of the disk isE[? j x(0) = (0;0)] = v(x) = 1ND?log 1" + log 2 + N4 ? log N?: (3.25)Notice that for N = 1 we recover the result obtained for one hole in the unitdisk given in Proposition 2.3.Lastly, we expect that as N gets larger, the MFPT tends to get smaller andsmaller. In fact as N ! 1, the MFPT given by (3.25), tends to 14D whenlog 1" >> N. This is depicted in Figure 3.3 below for D = 1, d = 12, " = 10?5and N = 100.50Chapter 3. Numerical Realizations: The Neumann Green's FunctionFigure 3.3: MFPT, v(x), from the centre of the unit disk with N symmetricallylocated holes on the boundary with D = 1, d = 12 and " = 10?5.3.2 The Unit SquareWe have found an expression for the mean ?rst passage time, given by (2.53), fora general two-dimensional domain with one hole on the boundary. In addition,we have found three formulas for the mean ?rst passage time with N holeson the boundary, namely (2.99), (2.101) and (2.103). We would now like toinvestigate the behaviour of the MFPT in the unit square.3.2.1 The Neumann Green's Function in a Unit SquareTo investigate the mean ?rst passage time in the unit square we must ?rst derivethe Neumann Green's function for the unit square. We start with a derivationof the Neumann Green's function for a rectangle.We will calculate the Neumann Green's function in a rectangle, where   =[0;L] ? [0;d] with j j = Ld. For the moment we place the singular point x0within the domain  . We will eventually take the limit as x0 ! @  so that wecan use the result to ?nd the mean escape time. We understand x = (x;y) so51Chapter 3. Numerical Realizations: The Neumann Green's Functionthat Gm(x;x0) = Gm(x;y;x0;y0), where x and y are the Cartesian coordinatesof x, and x0 and y0 are the Cartesian coordinates of x0. We need to solve(3.1)-(3.3) for the Neumann Green's function. That is, we need to solve?Gm(x;x0) = 1j j ? ?(x ? x0); x 2  ;@xGm(0;y;x0) = @xGm(L;y;x0) = 0; (3.26)@yGm(x;0;x0) = @yGm(x;d;x0) = 0:We impose the following integral condition to ensure uniqueness of the Neu-mann Green's function, Z Gm(x;x0)dx = 0: (3.27)3.2.2 Solution MethodWe expand Gm in terms of a complete set of eigenfunctions, and we use aseparation of variables technique to solve the resulting eigenvalue problem. Inother words, we use a Fourier series representation.We ?nd thatGm(x;x0) =1Xn=01Xm=0cm;n cos?n?xL?cos?m?yd?: (3.28)To enforce the integral condition, we must set c0;0 = 0. That is, we set thecoe?cient of the n = 0, m = 0 mode to zero. This ensures that the integralcondition is satis?ed.We then de?ne?0;n = cos?n?xL?;?m;0 = cos?m?yd?;?m;n = cos?n?xL?cos?m?xd?: (3.29)We decompose the n = 0, m = 0 modes as followsGm(x;x0) =1Xn=1c0;n?0;n +1Xm=1cm;0?m;0 +1Xn=11Xm=1cm;n?m;n: (3.30)52Chapter 3. Numerical Realizations: The Neumann Green's FunctionWe ?nd thatc0;n = 2L2j jcos?n?x0L ?n2?2 ;cm;0 = 2d2j jcos ?m?y0L ?m2?2 ;cm;n = 4j j cos?n?x0L?n2?2L2 +m2?2d2: (3.31)We will make use of the following identities1Xk=1cos (k?x)k2 + ?2 =?2?cosh(??(1 ? x))sinh(??) ?12?2 ; 0 ? x ? 2; (3.32)1Xk=1cos (k?x)k2 = ?2?16 ?x2 +x24?; 0 ? x ? 2: (3.33)To simplify the ?rst term in (3.30) we make use of the identity (3.33). We?nd that the ?rst term reduces toH(x;x0)= L212j j?4 ? 6?x + x0L?+ 3?x + x0L?2? 6?jx ? x0jL?+ 3?x ? x0L?2!(3.34)To simplify the third term in (3.31) we make use of identity (3.32). Thesecond term in (3.31) cancels with a term in this simpli?ed expression. Makingall the necessary simpli?cations we ?nd that (3.31) reduces toGm(x;x0) = H(x;x0) +12?1Xm=10@cosh?mL?d?1 ??jx+x0jL???+ cosh?mL?d?1 ??jx?x0jL???msinh?mL?d ?1A ??cos?m? (y + y0)d?+ cos?m? (y ? y0)d??: (3.35)We now make use of the following identitycosh(a ? b) + cosh(a ? c)sinha =11 ? e?2a?e?b + e?c + eb?2a + ec?2a?; (3.36)53Chapter 3. Numerical Realizations: The Neumann Green's Functionwith a = mL?d , b = mL?d jx+x0jL , c = mL?d jx?x0jL and q = e?2L?d .We expand the second term in (3.35) into exponentials. We de?ne the fol-lowing complex constants:r+;+ = ? jx + x0j + i(y + y0);r+;? = ? jx + x0j + i(y ? y0);r?;+ = ? jx ? x0j + i(y + y0);r?;? = ? jx ? x0j + i(y ? y0);?+;+ = jx + x0j + i(y + y0) ? 2L;?+;? = jx + x0j + i(y ? y0) ? 2L;??;+ = jx ? x0j + i(y + y0) ? 2L;??;? = jx ? x0j + i(y ? y0) ? 2L: (3.37)Then we introduce the new variablesz+;+ = er+;+?=d; z+;? = er+;??=d;z?;+ = er?;+?=d; z?;? = er?;??=d;?+;+ = e?+;+?=d; ?+;? = e?+;??=d;??;+ = e??;+?=d; ??;? = e??;??=d: (3.38)Using (3.37) and (3.38) in (3.35) and recalling that 11?qm = P1n=0 (qn)m, we?nd thatGm(x;x0) = H(x;x0) + 14?1Xm=11Xn=0(qn)m ?zm+;+ + ?m+;+?+ 14?1Xm=11Xn=0(qn)m ?zm+;? + ?m+;? + zm?;+ + ?m?;+ + zm?;? + ?m?;??+ 14?1Xm=11Xn=0(qn)m ??m+;+ + ?m+;+ + ?m+;? + ?m+;??+ 14?1Xm=11Xn=0(qn)m ??m?;+ + ?m?;+ + ?m?;? + ?m?;??: (3.39)Furthermore, we recall that Re?P1m=1 wmm currency1 = ? log j1 ? wj. Provided thateach of the z's and ?'s do not equal one, the double sum is absolutely convergent,and we can change the order of summation in (3.39). We perform the summation54Chapter 3. Numerical Realizations: The Neumann Green's Functionover m to obtainGm(x;x0) = H(x;x0)? 12?1Xn=0log j1 ? qnz+;+j j1 ? qnz+;?j j1 ? qnz?;+j j1 ? qnz?;?j? 12?1Xn=0log j1 ? qn?+;+j j1 ? qn?+;?j j1 ? qn??;+j j1 ? qn??;?j :(3.40)At this point we must extract the regular and singular part of the Green'sfunction as x ! x0. We observe that r?;? = 0 as x ! x0 causing thelog j1 ? z?;?j term to diverge when n = 0. We simplify as followslog j1 ? z?;?j = log jr?;?j + log????1 ? z?;?r?;?????: (3.41)Then we write the Green's function asGm(x;x0) = ? 12? log jx ? x0j + Rm(x;x0); (3.42)withRm(x;x0) = H(x;x0) ? 12? log????1 ? z?;?r?;????? ? 12?1Xn=1log j1 ? qnz?;?j? ? 12?1Xn=0log j1 ? qnz+;+j j1 ? qnz+;?j j1 ? qnz?;+j? 12?1Xn=0log j1 ? qn?+;+j j1 ? qn?+;?j j1 ? qn??;+j j1 ? qn??;?j : (3.43)We have derived the Neumann Green's function for a rectangle. We nowconsider the speci?c case with L = d = 1, the unit square. Now we determineexpressions for a singular point on the boundary, x0 2 @ . Using the symmetryof the square, we see that we need only ?nd the Green's function on one sideand on one corner. The method of extracting the regular and singular partsis identical to that used above. We state the results here. For a unit square,L = d = 1, with a singular point along the bottom side, that is y0 = 0 and55Chapter 3. Numerical Realizations: The Neumann Green's Functionx0 2 (0;1) we ?ndGm(x;x0) = ? 12? log jx ? x0j ? 12? logqjx ? x0j2 + jy + y0j2 + Rm(x;x0);(3.44)withRm(x;x0) = H(x;x0) ? 12? log????1 ? z?;?r?;????? ? 12? log????1 ? z?;+r?;+????? 12?1Xn=1log j1 ? qnz?;?j j1 ? qnz?;+j ? 12?1Xn=0log j1 ? qnz+;?j j1 ? qnz+;+j? 12?1Xn=0log j1 ? qn?+;+j j1 ? qn?+;?j j1 ? qn??;+j j1 ? qn??;?j : (3.45)As x ! x0 we ?ndGm = ? 1? log jx ? x0j + Rm(x0;x0); (3.46)withRm(x0;x0) = H(x0;x0) ? 1? log ? ? 1?1Xn=1log j1 ? qnj? 12?1Xn=0log ??1 ? qnz0+;+????1 ? qnz0+;???? 12?1Xn=0log ??1 ? qn?0+;+????1 ? qn?0+;?????1 ? qn?0?;+????1 ? qn?0?;???;(3.47)where z and ? are de?ned as in (3.38) withr0+;+ = ?2x0; r0+;? = ?2x0;r0?;+ = 0; r0?;? = 0;?0+;+ = 2x0 ? 2; ?0+;? = 2x0 ? 2;?0?;+ = ?2; ?0?;? = ?2: (3.48)Notice that the singularity is now twice what it was with the singular pointwithin the unit square.For a unit square, L = d = 1, with a singular point at the left end corner,56Chapter 3. Numerical Realizations: The Neumann Green's Functionthat is x0 = 0 and y0 = 0 we ?ndGm(x;0) = ? 12? log jx ? x0j ? 12? logqjx + x0j2 + jy + y0j2? 12? logqjx + x0j2 + jy ? y0j2 ? 12? logqjx ? x0j2 + jy + y0j2 + R(x;x0);= ? 2? log jxj + Rm(x;0); (3.49)withRm(x;0) = H(x;0) ? 12? log????1 ? z+;+r+;+????????1 ? z+;?r+;?????????1 ? z?;+r?;+????????1 ? z?;?r?;?????? 12?1Xn=1log j1 ? qnz+;+j j1 ? qnz+;?j j1 ? qnz?;+j j1 ? qnz?;?j? 12?1Xn=0log j1 ? qn?+;+j j1 ? qn?+;?j j1 ? qn??;+j j1 ? qn??;?j : (3.50)As x ! 0 we ?ndGm(0;0) ! ? 2? log jxj + Rm(0;0); (3.51)withRm(0;0) = H(0;0) ? 2? log ? ? 4?1Xn=1log j1 ? qnj ; (3.52)with H(0;0) = 13. Notice that the log jx ? x0j term in (3.51) is multiplied by2?. This means that singularity at a corner is four times larger than a singularpoint in the interior.Now that we have expressions for the Neumann Green's function in a unitsquare, we can proceed to the calculation of the MFPT using (2.17).3.2.3 The Mean First Passage Time in a Unit Squarewith One Hole on the BoundaryWe can use (2.53) for the MFPT for one hole along with the expressions for theGreen's function for a unit square given above. Since all of the expressions forthe Green's function involve an in?nite series, it is more instructive to plot theresults until a certain tolerance is reached. That is, until the terms in the seriesare negligible in comparison to some speci?ed tolerance. We use a tolerance57Chapter 3. Numerical Realizations: The Neumann Green's Function10?12. D = 1 and d = 12 in all cases investigated.Firstly, we consider one hole in the unit square located at the position (1;0:5).The initial position starts at (0;0:5) and moves closer to the hole in ?xed incre-ments of 180. The initial position does not reach the hole. We investigate thebehaviour of the MFPT in this case with D = 1, d = 12 and " = 10?5. Thecon?guration and the MFPT are shown in Figure 3.4.58Chapter 3. Numerical Realizations: The Neumann Green's Function(a) Plot of the initial position (blue), as it moves from the antipodal point at(0,0.5) towards the singular point (red) at (1,0.5)(b) Plot of the MFPT, v(x), versus the x-coordinate of the initial positionFigure 3.4: MFPT for the unit square as the initial position moves from (0;0:5)towards the exit window at (1;0:5) with D = 1, d = 12 and " = 10?5.59Chapter 3. Numerical Realizations: The Neumann Green's FunctionWe see that the MFPT decreases as the initial position moves towards theexit window.Now we consider the case where we hold the initial position ?xed at thecentre of the unit square, at x = (0:5;0:5). We place the exit window a smallincrement, 180, away from the corner at (1;0). We then move the exit windowin increments of 180 from (1;0:0125) to (1;0:9875) and calculate the MFPT withD = 1, d = 12 and " = 10?5. The con?guration is shown below in Figure3.5 along with the MFPT plotted against the y-coordinate of the moving exitwindow.60Chapter 3. Numerical Realizations: The Neumann Green's Function(a) Plot of the initial position (blue), and the singular point (red) as it movesalong the side(b) Plot of the MFPT, v(x), versus the y-coordinate of the exit windowFigure 3.5: MFPT from the centre of the unit square as the exit window movesfrom (1;0:0125) to (1;0:9875) with D = 1, d = 12 and " = 10?5.61Chapter 3. Numerical Realizations: The Neumann Green's FunctionWe observe symmetric behaviour in the MFPT, with a minimum where thedistance from the exit window to the initial position is shortest. Our results areconsistent with what we would expect.3.2.4 The Mean First Passage Time in a Unit Squarewith Two Holes on the BoundaryNow we consider two holes in a unit square. We ?x the initial position at thecentre of the unit square at x = (0:5;0:5). One exit window is ?xed at (0;0:5)while the other moves from (0;0:48) towards the corner at (0;0). Then it movesfrom (0;0) to (1;0) and ?nally from (1;0) to (1;0:5). We plot the behaviour ofthe MFPT as a function of the moving coordinate of the moving exit window oneach side, excluding the corner point. We then plot the total behaviour, as theexit window moves from (0;0:48) to (1;0:5), including the corner points, withthe points numbered from 1 to 100. We set D = 1, d = 12 and " = 10?5. Theresults are plotted shown in Figures 3.6 and 3.7 below.62Chapter 3. Numerical Realizations: The Neumann Green's Function(a) Plot of the MFPT as a function of the y-coordinate of the moving exitwindow along side 1, from (0,0.48) to (0,0.02)(b) Plot of the MFPT as a function of the x-coordinate of the moving exitwindow along side 2, from (0.02,0) to (0.98,0)Figure 3.6: MFPT from the centre of a unit square with two exit windows, one?xed at (0,0.5) and the other moving from (0,0.48) to (1,0.5), sides 1 and 2 withD = 1, d = 12 and " = 10?5.63Chapter 3. Numerical Realizations: The Neumann Green's Function(a) Plot of the MFPT as a function of the y-coordinate of the moving exitwindow along side 3, from (1,0.02) to (1,0.5)(b) Plot of the MFPT as a function of the moving exit window from (0,048) to(1,0.5), including the corner points. The position of the moving exit windowis labelled from 1-100Figure 3.7: MFPT from the centre of a unit square with two exit windows, one?xed at (0,0.5) and the other moving from (0,0.48) to (1,0.5), side 3 and allsides with D = 1, d = 12 and " = 10?5.64Chapter 3. Numerical Realizations: The Neumann Green's FunctionNotice that the local minima in Figures 3.6 (a) and (b) are not where thedistance of the moving exit window from the initial position is a minimum. Thisis consequence of the fact that there are two holes, and there is an interactionbetween the two holes, as well as with the initial position. We see that theminimum is shifted to the right slightly. In Figure 3.7 (a), we do ?nd a minimumfor the MFPT in the expected position. This is a result of the symmetry of theplacement of the two holes. In fact, if we continue to move the second exitwindow from (1;0:02) to (1;0:98), the plot of the MFPT is symmetric in thiscase.Also notice the discontinuities in Figure 3.7 (b) at the corner points. Thisis expected, since the boundary is non-smooth at these points. We remark thatour formula for the MFPT is not valid when the arc is at the corner of thesquare. The local geometry at the corner of a square is di?erent. Therefore, weneed a di?erent form for the outer solution.If we look at the expression for the Neumann Green's function at the corner,given by (3.49), (3.50), (3.51), (3.52), we see that the singularity is twice as largeat the corner than a singularity on the smooth portion of the boundary. Thisagrees with the result for a rectangle obtained by Singer et al. in [24]. Here,Singer et al. state that the leading order term for the MFPT is twice as large asthe case with the singular on the smooth portion of the boundary. We proceedto derive an expression for the MFPT from the centre of the unit square with asingle arc at the corner, x0 = (0;0).We alter the method outlined in Chapter 2, Section 2.2. The outer solutionfor the principal eigenfunction, ucurrency1, is written in terms of the Helmholtz Green'sfunction given by (2.31) and (2.32). We alter (2.32) so thatG(x;x0;?currency1) = ? 2? log jx ? x0j + R(x;x0;?currency1): (3.53)Thus, ucurrency1 = ??A?2 G(x;x0;?currency1). Matching to the inner solution, we ?nd thatR(x0;x0;?currency1) = ? 2?? . We expand G(x;x0;?currency1) as in (2.37). We ?nd that Gmsatis?es (2.43) with Gm(x;x0) = ? 2? log jx ? x0j + Rm(x;x0). We ?nd theprincipal eigenvalue and eigenfunction to be?0(?) = ??2j j ? ?2?24j j Rm(x0;x0) + O(?3)?0(x;?) = j j?12 ? ?? j j?122 Gm(x;x0) + O(?3): (3.54)65Chapter 3. Numerical Realizations: The Neumann Green's FunctionWe calculate the MFPT to bev(x) = 2j j?D?? log("d) + ?2 (Rm(x0;x0) ? Gm(x;x0))?+ O(?); (3.55)with the average MFPT given by? = 2j j?D?? log("d) + ?2 Rm(x0;x0)?+ O(?): (3.56)Now, for the corner point x0 = (0;0), Rm(0;0) is given by (3.52) withH(0;0) = 13 and Gm(x;0) is given by (3.49) and (3.50). Finally, the constantd in (3.55) is determined from the far-?eld behaviour of the inner problem. Inparticular, d depends on the manner in which the absorbing arc of length 2"is placed at the corner. If the arc is placed along one side so that v = 0 on0 < x < 2" with y = 0, then d = 1. If the arc is placed along two sides suchthat v = 0 on y = 0 with 0 < x < " and on x = 0 with 0 < y < ", then d = 14.Taking an arc of length " placed along one side such that v = 0 on y = 0with 0 < x < ", we ?nd that d = 12. Approximating Rm(0;0) by the ?rst termin (3.56) we ?nd that? ? 2j jD??log 1" + log 2? + ?6 + 2e?2??: (3.57)This agrees with equation (2.8) of [24], which was found by solving certainintegral equations asymptotically.3.3 Arbitrary Domains - A NumericalApproachAs we have seen, to investigate the MFPT in a particular domain we must ?ndthe Neumann Green's function and its regular part. We have succeeded in doingthis for a unit circle and square. However, other geometries are more challengingas an analytical solution for the Green's function is unattainable in most cases.In this section we introduce a numerical procedure to calculate the NeumannGreen's function and its regular part in an arbitrary domain.66Chapter 3. Numerical Realizations: The Neumann Green's Function3.3.1 The Boundary Element MethodWe introduce the boundary element method, the method we will use to cal-culate the Neumann Green's function in an arbitrary domain. We alter thestatement of the problem for the Neumann Green's function slightly, by puttingthe singularity directly on the boundary as?Gm(x;x0) = 1j j; x 2  ; (3.58)@nGm(x;x0) = ?(s ? s0); x 2 @ ; (3.59)Z Gm(x;x0)dx = 0: (3.60)where s is an arclength parameter and the normal derivative is taken withrespect to the unit outward normal. As mentioned, the condition that theintegral of the Green's function over the domain is zero ensures that the solutionto the problem is unique.We decompose Gm(x;x0) asGm(x;x0) = ? 1? log jx ? x0j + R(x;x0) + 14j j?jxj2 + jx0j2?(3.61)We will de?ne w asw(x;x0) = 14j j?jxj2 + jx0j2?: (3.62)The Neumann Green's function is given by Gm(x;x0) = ? 1? log jx ? x0j +Rm(x;x0). Notice that the regular part of the Neumann Green's function in(3.61) is given byRm(x;x0) = R(x;x0) + 14j j?jxj2 + jx0j2?: (3.63)We solve for R(x;x0) numerically.Substituting (3.61) into (3.58)-(3.60) we ?nd that R(x;x0) must satisfy?R(x;x0) = 0; x 2  ;@nR = ?@nw; x 2 @ ;Z Rdx = ?Z ?1? log jx ? x0j ? w?dx: (3.64)67Chapter 3. Numerical Realizations: The Neumann Green's FunctionNote that ?w = 1j j.We now apply Green's identity given by(u1;?u2) ? (u2;?u1) =Z@ u1@nu2dS ?Z@ u2@nu1dS: (3.65)We let u1 = R(x;x0) and u2 = g(x0;?), where g(x;?) is the free-spaceGreen's function satisfying?g(x;?) = ??(x ? ?); x 2  ;which is given byg(x;?) = ? 12? log jx ? ?j : (3.66)Substituting u1 and u2 into (3.65) we ?nd that the integral equation for R atx 2   isR(x;?) = ?Z@ R(x;?)@ng(?;?)dS(?) ?Z@ g(x;?)@nw(?;?)dS(?): (3.67)Now we discretize the boundary @  into n arcs @ 1;@ 2;:::;@ n. We canthen use the midpoint rule to approximate the integrals in (3.67). For example,we can approximate R(x;?) by R(x;?i) where ?i is the midpoint of each arc@ i. We let Rj = R(x;?j), using this we ?ndRj = ??ni=1 (aijRi + bij) ; (3.68)where aij = R@ i@ng(?;?j)dS(?), bij = R@ ig(x;?)@nw(?;?j)dS(?) and Ri =R(x;?i).Now we compute aij and bij by the midpoint rule to ?ndaij = li@ng(?i;?j); bij = g(x;?i)li@nw(?i;?j); (3.69)where li is the length of the arc @ i. Using this we ?nd that we can representR byR(x;?j) +nXi=1liR(x;?i)@ng(?i;?j) = ?nXi=1lig(x;?i)@nw(?i;?j): (3.70)Firstly, we calculate @nw(?i;?j) in (3.70). We denote the x and y component68Chapter 3. Numerical Realizations: The Neumann Green's Functionof ?i as ?ix and ?iy respectively. We obtain thatrw(?i;?j) = 12j j(?ix;?iy);@nw(?i;?j) = 12j j(?ix;?iy):^: (3.71)Hence @nw is independent of its second argument ?j. Now we calculate @ng(?i;?j),rg(?i;?j) ? ^ = ? 12? ?i ? ?j??i ? ?j??2? ^: (3.72)When calculating @ng(?i;?j), the case i = j requires special care becauseof the singularity of the Green's function. We must calculate the ?rst integralin (3.67) di?erently. We introduce a local coordinate system. We let 1?i bethe radius of curvature of @ i at ?i where ?i is the curvature. Since we haveassumed each arc @ i to be very small, we may assume that @ i is parametrizedfor t << 1 as?(t) = 1?i(cos t;sin t); ? li2?i? t ? li2?i; (3.73)with ?i =?1?i ;0?. We compute @ng(?;?) = rg:^ where ^ = (cost;sin t) is theunit outward normal in the local coordinate system. We ?ndrg = ? 12? ? ? ?ij? ? ?ij2 ;? ? ?i =? 1?i cos t ?1?i ;1?i sint?;^:(? ? ?i) = 1?i(1 ? cost);j? ? ?ij2 = 2?2i(1 ? cos t): (3.74)Therefore@ng(?;?i) = ? 12?1?i (1 ? cost)2 1?2i(1 ? cos t) = ??i4?: (3.75)Thusaii = ?li?i4? : (3.76)We do not encounter any problems when i = j for bij.69Chapter 3. Numerical Realizations: The Neumann Green's FunctionWe can represent (3.70) as a matrix system AR = ?B as follows:R =?R(x;?1) R(x;?2) ? ? ? ? ? ? R(x;?n)?T; (3.77)A =0BBBBBBBBB@?1 ? l14?r1?l2@ng(?2;?1) l3@ng(?3;?1) ? ? ? ln@ng(?n;?1)l1@ng(?1;?2)?1 ? l24?r2?l3@ng(?3;?2) ? ? ? ln@ng(?n;?2)l1@ng(?1;?3) l2@ng(?2;?3)?1 ? l34?r3?? ? ? ln@ng(?n;?3)... ... ... ... ...l1@ng(?1;?n) l2@ng(?2;?n)) ? ? ? ? ? ??1 ? ln4?rn?1CCCCCCCCCA;B =0BBBBBBB@l1g(x;?1)@nw(?1;?1) + l2g(x;?2)@nw(?2;?1) + ? ? ? + lng(x;?n)@nw(?n;?1)l1g(x;?1)@nw(?1;?2) + l2g(x;?2)@nw(?2;?2) + ? ? ? + lng(x;?n)@nw(?n;?2)l1g(x;?1)@nw(?1;?3) + l2g(x;?2)@nw(?2;?3) + ? ? ? + lng(x;?n)@nw(?n;?3)...l1g(x;?1)@nw(?1;?n) + l2g(x;?2)@nw(?2;?n) + ? ? ? + lng(x;?n)@nw(?n;?n)1CCCCCCCA;We solve for R using Gaussian elimination.Our goal is to calculate the MFPT. To this end, we must determine R(x0;x0).Using this formulation, this means that we must ?nd R(?j;?j). Thus, we pickx = ?j for some j = 1;2;:::;N in (3.70). This introduces another case wherethe possibility of i = j arises. We have to recalculate R@  g(?;?)@nw(?;?)dS(?)in (3.67). We use the same local coordinate system as outlined above. Thus,we have g(?;?) = ? 12? log jx ? ?j = ? 14? log?2?2i (1 ? cos t)?and @nw(?;?) =12j j?i . Once again,1?i is the radius of curvature in the local coordinate systemin segment i. The integral becomes? 18? j j ?2iZ li?i=2?li?i=2log? 2?2i (1 ? cost)?dt: (3.78)Notice that when t = 0 we have a singularity, so we can use gaussian quadratureto get around this by using an even number of sample points. However, forcertain domains we remove the singularity by rewriting the integral, as will bediscussed in later sections.70Chapter 3. Numerical Realizations: The Neumann Green's FunctionNotice that we have not imposed the integral condition. We must still deter-mine the constant C so that (3.64) is satis?ed. To do this, we must integrate ournumerical result for R over the domain  , which proves to be a computationallyintensive process. We now illustrate the use of the method for a circle and anellipse.3.3.2 The Unit DiskWe make use of the boundary element method described in the previous sectionto determine the Neumann Green's function for a unit disk. This allows us tocompare our numerical method to the analytical result given by (3.9).According to the formulation above, we letGm(x;x0) = ? 1? log jx ? x0j + 14j j?jxj2 + jx0j2 j?+ R(x;x0) + C (3.79)Thus, the unknown is R(x;x0) and C is to be determined from the integralcondition once R has been found. If we compare the expression above (3.79)with (3.9), we see that we already have the expression for the unit circle. Oursolution for R in this case should be zero for all x. In this section we con?rm thatthe boundary element formulation does in fact give us this result. We illustratethe results in terms of N, the number of segments used for the calculation. Weaim to attain zero to six decimal places. Note that since we are essentiallylooking for a constant, each entry of the vector R, (3.77), will be the same.We begin with the calculation of R(x;x0) for x 2  . We notice, by consid-ering Table 3.1 below, that as x moves towards the boundary of the unit disk,we require more segments to attain the desired accuracy.N 40 80 100x = (0;0) ?0:3534 ? 10?17 ?0:6184 ? 10?17 ?0:5624 ? 10?17x = (0:5;0) 0:7237 ? 10?14 0:1084 ? 10?16 0:2602 ? 10?16x = (0:75;0) 0:8003 ? 10?7 0:4024 ? 10?12 0:1016 ? 10?14x = (0:99;0) 0:0041 0:0015 0:9928 ? 10?3Table 3.1: R(x;x0) for the unit disk with N mesh points on the boundaryTable 3.1 shows the value of R(x;x0) for each x shown above, and for all N71Chapter 3. Numerical Realizations: The Neumann Green's Functionx0's on the boundary. The x0's on the boundary are the ?i's in the boundaryelement formulation. There is only one value shown for the vector R since eachentry of the vector is the same, since we are looking for a constant.It is clear that as x moves towards the boundary more mesh points are re-quired to attain the desired result of zero to six decimal places, whereas with x atthe centre of the disk we attain convergence immediately. This is a consequenceof the fact that as x tends towards the boundary, the self e?ect starts to play animportant role. For N = 800, with x = (0:99;0) we obtain R = 0:1282 ? 10?6.To calculate R(x0;x0) we need to recalculate the second integral in equation(3.67) as outlined in the previous section. We use integral (3.78) instead. Wecalculate the integral (3.78) using a four point Gaussian quadrature rule. Wedo not rewrite the integral to remove the singularity in this case. As long aswe do not use an odd number of sample points in the Gaussian quadratureprocedure, we are able to avoid singular contributions to the integral. We alsoanticipate that we will require many mesh points to attain the desired accuracy.Additionally, as a result of the formulation of the boundary element method, weare able to check the symmetry of R in this case. In other words, for x = ?1 werequire R(?1;?j) = R(?j;?1) for each j. The symmetry of the problem ensuresthat R(xj;xj) are identical for each j. We tabulate the results in Table 3.2belowN 50 100 200 400 800R(x0;x0) 0:0025 0:0012 0:6193 ? 10?3 0:3097 ? 10?3 0:1548 ? 10?3Table 3.2: R(x0;x0) for the unit disk with N mesh points on the boundaryWe see that as we double N, the value of R divides by 2. Thus, we canconclude that R(x0;x0) is tending to zero. We are left to calculate the integralover Gm, equation (3.79). It is simple to calculate the integral over the domainof each term in equation (3.79) with the exception of the integral over R. Weuse the composite trapezoidal rule to numerically compute the integral overR and log jx ? x0j. As shown in (3.8), the integral of Gm over the domainis independent of x0. We calculate the desired integrals for several x0's. Weexpect each integral to return the same value. This is another way to check thatthe numerical method is working correctly. We describe the integration method72Chapter 3. Numerical Realizations: The Neumann Green's Functionover R. We need to calculateZ 2?0Z 10R(x;x0)rdrd? (3.80)We ?rst split the integral over ? using the composite trapezoidal rule with nintervals. There will be n remaining integrals over r. We evaluate the remainingn integrals using the midpoint rule. We run our boundary element method forthe values of x that are required to compute the integral. The integration overlog jx ? x0j is similar. Notice that we compute the integrals over r using themidpoint rule to avoid any singularities.We ?nd that the integral over log jx ? x0j is zero for all x0 and we ?nd theintegral over R to be zero for all x0.We are now ready to calculate the constant C. Since the integral overR and log jx ? x0j are zero, we are left to calculate the integral over w =14j j?jxj2 + jx0j2?. For any x0 on the boundary, jx0j = 1. We ?nd C = ? 38?,which agrees with what we ?nd analytically.We have con?rmed that the boundary element method does in fact give usthe desired results for a unit circle.3.3.3 The EllipseWe are now ready to use the boundary element method to determine the regularpart of the Neumann Green's function for an ellipse. We state the propertiesthat will be required for the boundary element method. For a curve de?ned inpolar coordinates, we have r = r(?);x = r cos(?);y = r sin(?). For an ellipser = aba2 sin2(?) + b2 cos2(?); (3.81)r0 = ? absin(?)cos(?)(a2 ? b2)a2 sin2(?) + b2 cos2(?)(3=2) ; (3.82)r00 = ?ab(a2 ? b2) ??cos2(?) ? sin2(?)(a2 sin2(?) + b2 cos2(?))3=2 ?3(a2 ? b2)sin2(?)cos2(?)(a2 sin2(?) + b2 cos2(?))5=2?;(3.83)where a is the semi-major axis and b is the semi-minor axis.73Chapter 3. Numerical Realizations: The Neumann Green's FunctionThe outward normal vector, the arclength and curvature are given by^ = 1r02 + r2 (r cos(?) + r0 sin(?);r sin(?) ? r0 cos(?)); (3.84)s =Z ?2?1pr02 + r2d?; (3.85)? = r2 + 2r02 ? rr00(r02 + r2)3=2 : (3.86)To determine the length of a segment li we will need to calculate the integral(3.85) for the arclength. We compute this numerically using the trapezoidalrule.The algorithm outlined for the boundary element method holds for a generaldomain. We consider an ellipse with a = 2 and b = 1. Once again, we beginwith the calculation for R(x;x0) for x 2  . We use more mesh points as xmoves towards the boundary of the ellipse. The results are tabulated below inTable 3.3.N 80 160 320 640x = (0;0) 0:1290 0:1291 0:1291 0:1291x = (1;0) 0:1291 0:1291 0:1291 0:1291x = (1:75;0) 0:1291 0:1291 0:1291 0:1291x = (1:99;0) 0:1331 0:1305 0:1294 0:1291x = (0;0:5) 0:1290 0:1291 0:1291 0:1291x = (0;0:99) 0:1298 0:1292 0:1291 0:1291x = (1;0:5) 0:1291 0:1291 0:1291 0:1291Table 3.3: R(x;x0) for the ellipse with N mesh points on the boundaryTable 3.3 shows the value of R(x;x0) for each x within the ellipse, and forall x0 on the boundary. We ?nd that each entry of the vector R is identical,regardless of where we place x or x0. This implies that R is a constant, R =0:1291. The behaviour of the Green's function for the ellipse is thus containedwithin log jx ? x0j and w(x;x0). We have found thatGm(x;x0) = ? 1? log jx ? x0j + 14j j?jxj2 + jx0j2?+ 0:1291 (3.87)Notice that we have not imposed R  G(x;x0)dx = 0, and thus the expression for74Chapter 3. Numerical Realizations: The Neumann Green's Functionthe Neumann Green's function, (3.87), for the ellipse is unique up to a constant,C. Below, in Figure 3.8 we show the behaviour of the Neumann Green's functionfor the ellipse, Gm ? C.Figure 3.8: Plot of the Neumann Green's function for the ellipse, Gm?C, versus?0In Figure 3.8 we have x at the origin and x0 at points along the boundary inincrements of ?100. We observe the symmetric behaviour of the Green's functionfor the ellipse, which is consistent with our expectations.We proceed to the calculation of R(x0;x0). We compute the integral (3.78)as followsZ li?i=2?li?i=2log? 2?2i (1 ? cost)?dt = li?i log? 2?2i?+Z li?i=2?li?i=2log(1 ? cost)dt:(3.88)We rewrite the last integral in (3.88) as=Z li?i=2?li?i=2log(1 ? cost) ? log(t2=2)dt +Z li?i=2?li?i=2log(t2=2)dt; (3.89)=Z li?i=2?li?i=2log(2(1 ? cos t)=t2)dt + 2li?i log li?i ? 2li?i ? 3li?i log 2:75Chapter 3. Numerical Realizations: The Neumann Green's FunctionWe compute the remaining integral using Gaussian quadrature, with an evennumber of sample points. Initially we use weight one with sample points ? 1p3,two point Gaussian quadrature. Notice that since we have observed that R isa constant, R(x0;x0) should be the same for each x0. We encounter extremelyslow convergence of the numerical method in this case. We require many meshpoints in order for the method to converge to an answer that we expect, thatis, all the entries of the vector R(x0;x0) should be identical and each vectorR should be the same regardless of x0. In Table 3.4 below, we tabulate themaximum di?erence between the entries in any vector R(x0;x0).N 150 300 600 1200 2400min 0:1136 0:1202 0:1240 0:1263 0:1275max 0:1429 0:1371 0:1336 0:1316 0:1305jmax ? minj 0:0293 0:0169 0:0096 0:0054 0:0030R(x0;x0) 0:1209 0:1244 0:1264 0:1276 0:1282Table 3.4: R(x0;x0) for the ellipse with N mesh points on the boundaryThe last row in Table 3.4 shows the average value of R(x0;x0). We ?ndthat the vector R satis?es the expected properties when the mesh points areclustered very closely together, at ? = ?=2;3?=2. At ? = 0;? the points arefar apart along the boundary due to the curvature of the ellipse. It is here thatwe obtain the largest discrepancy between values of R. We have been unableto get the method to converge for N = 2400. We propose altering the mannerin which we have computed the integrals for the boundary element method byusing Gaussian quadrature instead of the midpoint rule to compute the integralsin (3.67). In addition, we did attempt to use four point Gaussian quadratureto compute the last integral in (3.90). However, there was no improvement inthe convergence. In fact, we obtained exactly the same results as we did withtwo point Gaussian quadrature. The reason for this is that the integrand is notsu?ciently di?erentiable in order for four point Gaussian quadrature to makea di?erence. We expect that R(x0;x0) = 0:1291. We see, by considering theaverage values of R, that R is tending towards this value.We are left to calculate the integral over Gm, (3.87), to determine theconstant C. We use the same procedure as we used for the unit circle. We?nd the same value for the integral over R regardless of the x0 we choose.76Chapter 3. Numerical Realizations: The Neumann Green's FunctionThis agrees with the fact that integral over Gm should be independent of x0.We ?nd that R  R(x;x0)rdrd? = 0:8108. We know that since R is a con-stant, that R  R(x;x0)dx = 2? ? 0:1291 = 0:8112. Our numerical answeris correct up to three decimal places. It is simple to calculate the integralR 14j j?jxj2 + jx0j2?dx = 14j j?7:85398 + 2? jx0j2?. However, the integralover ? 1? log jx ? x0j is not simple to calculate as a consequence of the singular-ity when x = x0. We are only able to estimate this integral up to one decimalplace accurately. As a result, we ?nd that C = ? 12? (0:8108 + 0:1) = ?0:1450We can now ?nd the MFPT for the ellipse. We know the Neumann Green'sfunction for the ellipse, (3.87). In addition, we can assume that R(x0;x0) =0:1291. To calculate the MFPT for the ellipse with one hole on the boundary weuse (2.53). It is worth noting that the constant, C, does not play a role directlyin this equation since the constant term from Rm(x0;x0) and Gm(x;x0) canceleach other. However, it is assumed that we are able to ?nd the constant C suchthat R  Gm(x;x0)dx = 0 since it is an implementation of the normalizationcondition. It is only in the case of normalized eigenfunctions that (2.53) holds.We place the initial position at the origin. We place one hole on the bound-ary, initially at an angle ?100, which moves from this position at equal angularincrements of ?50 until it reaches ? ?100. We set D = 1, d = 12 and " = 10?5. TheMFPT from the centre of an ellipse with semi-major axis 2 and semi-minor axis1 is depicted in Figure 3.9.77Chapter 3. Numerical Realizations: The Neumann Green's FunctionFigure 3.9: Plot of the MFPT, v(x), from the centre of the ellipse versus ?0with D = 1, d = 12 and " = 10?5We observe the expected symmetric behaviour in the MFPT. The MFPT isa minimum at ? = ?=2, where the distance from the centre is a minimum. Wesee a maximum in the MFPT at ? = 0;?, where the distance from the centre isa maximum.Now, we move the initial position from the origin to x = (1;0). The MFPTin this case is depicted in Figure 3.10.78Chapter 3. Numerical Realizations: The Neumann Green's FunctionFigure 3.10: Plot of the MFPT, v(x), from x = (1;0) for the ellipse versus ?0with D = 1, d = 12 and " = 10?5Since the initial position is shifted in this case, we ?nd that the MFPT hasa minimum at ? = ?19?=100 and a maximum at ? = ?, where the hole is thefurthest from the initial position.3.4 Perturbed Circular Domains - AnalyticalApproachWe now derive expressions for the regular part for the Neumann Green's functionin a perturbed circular domain. This work was completed by T. Kolokolnikov,[7]. We present an overview of his work.3.4.1 DerivationWe want to solve system (2.43) for the Neumann Green's function, with theexception of the integral condition. That is, we solve for the Green's function79Chapter 3. Numerical Realizations: The Neumann Green's Functionup to an arbitrary constant. We have the general expressionGm(x;x0) = ? 1? log jx ? x0j + R(x;x0); x0 2 @ : (3.90)We want to determine Rm(x;x0) and Rm(x;x0) in the limit as x ! x0 for do-mains that are near the unit disk. We parametrize the boundary of the domain, , in polar coordinates. We propose the following theorem:Theorem 3.9: Suppose that the domain   is de?ned in polar coordinatesbyr = r(?) = 1 + "?(?); " ? 1: (3.91)Suppose that ?(?) is given by?(?) =1Xn=1(an cosn? + bn sin n?) : (3.92)Let x0 = (r0 cos ?0;r0 sin ?0) be a point on the boundary where r0 = 1 + "?(?0).We de?ne?(?) = Rm(x;x0) and ?(?0) = Rm(x0;x0): (3.93)We have the following expression for ?0(?0)?0(?0) = "?Xn?1?n2 + n ? 2?(bn cos n?0 ? an sin n?0) (3.94)Proof: We ?nd that Rm(x;x0) satis?es?Rm(x;x0) = 1j j; x 2  ; (3.95)rRm(x;x0) ? ^ = 1? (x ? x0) ? ^n(x)jx ? x0j2 ; x 2 @ ; (3.96)where the unit outward normal ^ is de?ned by (3.84). We compute the expres-sion on the right-hand side of (3.96) for " << 1 as1?(x ? x0) ? ^(x)jx ? x0j2 =12??1 + "?? cos(? ? ?0) ? ?0 ? ?0 sin(? ? ?0)1 ? cos(? ? ?0)??+ O("2):(3.97)80Chapter 3. Numerical Realizations: The Neumann Green's FunctionThe expression in square brackets is bounded for ? ! ?0. Thus the expressionis uniformly valid for all ? 2 [0;2?). We de?ne f(?) asf(?) = ? cos(? ? ?0) ? ?0 ? ?0 sin(? ? ?0)1 ? cos(? ? ?0) : (3.98)Then we express f(?) in terms of a Fourier series asf(?) =1Xn=1(An cosn(? ? ?0) + Bn sin n(? ? ?0)); (3.99)where An and Bn are given byAn = 1?Z 2?0f(?)cos m(? ? ?0)d?; Bn = 1?Z 2?0f(?)sin m(? ? ?0)d?:(3.100)We denote I1 and I2 by I1 = ?An and I2 = ?Bn.Firstly, we consider the case where ? = cos n? = Re?ein??. We need tocalculateI1 = ReZ 2?0cos(? ? ?0)ein? ? ein?0 ? inein? sin(? ? ?0)1 ? cos(? ? ?0) cos m(? ? ?0)d?:(3.101)Let z = ei? and z0 = ei?0. The equation above becomesI =Z zn? zz0 +zz0?12?? zn0 ? nznzz0 ?zz0?122 ? zz0 + zz0 ?1??zz0?m+? zz0??m! dziz :(3.102)Let w = zz0 . The integral becomesI =Z 1?n2 wn+1 + 1+n2 wn?1 ? 1w2 ? 2w + 1?wm + w?m?dw: (3.103)The integration is performed over the boundary of the unit disk.Note that we can write 1w2?2w+1 as1(1 ? w)2 =ddw11 ? w =ddw1Xn=0wn: (3.104)Notice that jwj < 1 since w is within the unit circle. After some simpli?cation81Chapter 3. Numerical Realizations: The Neumann Green's Functionwe ?nd1?n2 wn+1 + 1+n2 wn?1 ? 1w2 ? 2w + 1= ??1 + 2w + 3w2 + ? ? ? + (n ? 1)wn?2 + (n ? 1)2 wn?1?: (3.105)We use the residue theorem to calculate the integral. We ?nd thatI = zn0? 2?m; 1 ? m < n;?(n ? 1); m = n;0; m > n:(3.106)We can now state the following results for I1 and I2. For ?0 = cosn?0:I1 = cos n?0? 2?m; 1 ? m < n;?(n ? 1); m = n;0; m > n;(3.107)andI2 = ? sin n?0? 2?m; 1 ? m < n;?(n ? 1); m = n;0; m > n:(3.108)For ?0 = sin(n?0):I1 = sin n?0? 2?m; 1 ? m < n;?(n ? 1); m = n;0; m > n;(3.109)andI2 = cos n?0? 2?m; 1 ? m < n;?(n ? 1); m = n;0; m > n:(3.110)We have now found An = 1?I1 and Bn = 1?I2. We have the following Fourierexpansions for f(?). For ?0 = cos n?0:f(?) =n?1Xm=12m[cos n?0 cosm(? ? ?0) ? sin n?0 sin m(? ? ?0)]+ (n ? 1)(cos n?0 cosn(? ? ?0) ? sin n?0 sin n(? ? ?0)) : (3.111)82Chapter 3. Numerical Realizations: The Neumann Green's FunctionFor ?0 = sin n?0:f(?) =n?1Xm=12m[cos n?0 sin m(? ? ?0) + sin n?0 cos m(? ? ?0)]+ (n ? 1)(cos n?0 sin n(? ? ?0) + sinn?0 cosn(? ? ?0)) : (3.112)Recall that ? = P1n=1 (an cos n? + bn sinn?) from (3.92). We can ?nd f(?)for ? given by (3.92) by superposition. In this way we ?ndf(?) =1Xn=1?an(n ? 1)?cn +n?1Xm=1an2m?cm!+1Xn=1?bn(n ? 1)?sn +n?1Xm=1bn2m?sm!; (3.113)where?cm = cosn?0 cos m(? ? ?0) ? sin n?0 sin m(? ? ?0); (3.114)?sm = cosn?0 sin m(? ? ?0) + sinn?0 cosm(? ? ?0): (3.115)We can rearrange f(?) and rewrite it asf(?) =1Xn=1[(n ? 1) (an cosn?0 + bn sin n?0)cos n(? ? ?0)]+1Xn=1n?1Xm=1[2m(an cosn?0 + bn sin n?0)cos m(? ? ?0)]+1Xn=1[(n ? 1) (bn cosn?0 ? an sin n?0)sin n(? ? ?0)]+1Xn=1n?1Xm=1[2m(bn cos n?0 ? an sin n?0)sin m(? ? ?0)]:(3.116)Then we interchange the order of summation as follows:1Xn=1n?1Xm=1?mn =1Xm=11Xn>m?mn =1Xn=11Xm>n?nm: (3.117)83Chapter 3. Numerical Realizations: The Neumann Green's FunctionUsing this we ?nd thatf(?) =1Xn=1An cos n(? ? ?0) + Bn sin n(? ? ?0); (3.118)An = (n ? 1)(an cos n?0 + bn sin n?0) + 2n1Xm>n(am cos m?0 + bm sin m?0);Bn = (n ? 1)(bn cos n?0 ? an sin n?0) + 2n1Xm>n(bm cos m?0 ? am sin m?0):(3.119)To proceed we de?ne S(x;x0) byRm(x;x0) = S(x;x0) + jxj24j j: (3.120)Substituting this into the system for R, (3.95) and (3.96), we ?nd that?S(x;x0) = 0; x 2  ; (3.121)@nS(x;x0) = @nRm(x;x0) ? 14j j@n jxj2 ; x 2 @ : (3.122)We calculate @n jxj2 = 2rp1+r02=r2 . Using this along with (3.96), (3.97) and(3.98) we ?nd that@nS(x;x0) = "2? (f(?) ? ?(?)) + O(")2; x 2 @ : (3.123)Note that j j ? ?. Now we introduce S0(x;x0) byS(x;x0) = "2?S0(x;x0): (3.124)To leading order we can write @nS0 = @rS0 jr=1 +O("). We obtain the followingleading order problem:?S0(x;x0) = 0; 0 < r < 1; (3.125)@rS0(x;x0) jr=1 = f(?) ? ?(?); r = 1: (3.126)84Chapter 3. Numerical Realizations: The Neumann Green's FunctionWe can solve this system using separation of variables. We ?nd thatS0 = D0 +1Xn=1rn (Dn cos n(? ? ?0) + En sin n(? ? ?0)) : (3.127)To solve for the coe?cients Dn and En we must use the boundary condition(3.126). To this end, we must rewrite ?, given by (3.92), in terms of cos n(???0)and sin n(? ? ?0). Notice that f(?) is already in this form, (3.99). We calculatefor ? that? =1Xn=1?~n cosn(? ? ?0) + ~n sin n(? ? ?0)?; (3.128)~n = an cos n?0 + bn sin n?0; (3.129)~bn = bn cos n?0 ? an sin n?0: (3.130)By imposing the boundary condition (3.126) we ?nd that@rS0 jr=1=1Xn=1n(Dn cos n(? ? ?0) + En sin n(? ? ?0)) : (3.131)Equating coe?cients with the right-hand side of (3.126) we obtainnDn = An ? ~n; nEn = Bn ? ~n: (3.132)To summarize, we have, for x 2 @ , thatRm(x;x0) = S(x;x0) + jxj24? ="2?S0(x;x0) +14? +"?2? + O("2): (3.133)Using de?nition (3.93), we calculate the derivative of ?0(?0) as follows?0(?0) = dd?0Rm(x0(?0);x0(?0)) = 2 dd?Rm(x(?);x0(?0)) j?=?0 : (3.134)Substituting (3.133) into (3.134) we ?nd?0(?0) = "?? dd?S0(x(?);x0(?0)) j?=?0 +?0(?0)?: (3.135)We di?erentiate S0(x;x0) given by (3.127) and ? given by (3.92) with respect85Chapter 3. Numerical Realizations: The Neumann Green's Functionto ?. We substitute these into the equation above to ?nd that?0(?0) = "?1Xn=1?nEn + n~n?: (3.136)Then, using equations (3.130) and (3.132) we ?nd?0(?0) = "?1Xn=1[Bn + (n ? 1)bn cos n?0 ? (n ? 1)an sinn?0]: (3.137)Next, we use the expression for BnBn = (n ? 1)(bn cosn?0 ? an sin n?0) + 2nXm>n(bm cosm?0 ? am sinm?0) ;(3.138)so that?0(?) = "?1Xn=1?2(n ? 1)?n + 2n1Xm>n?m!; (3.139)?m = bm cosm?0 ? am sin m?0: (3.140)We calculate ? = P1n=1 P1m>n n?m. We can change the order of summationsuch that ? = P1m=2 ?m Pm?1n=1 n. We ?nd that ? = P1n=1 n(n ? 1)?n. Thusour ?nal result for ?0(?0) is?0(?0) = "?1Xn=1(n + 2)(n ? 1)(bn cosn?0 ? an sin n?0) : (3.141)?We have found R0m(x0;x0).3.4.2 ExampleWe present an example here that makes use of (3.141). We determine whetherthere is a relationship between the curvature of @  and where the regular partof the Green's function, Rm, has its extrema. In addition, we compare thesepredictions with Rm from our boundary element method.Consider the shape with ? de?ned by?(?) = cos 2? + acos3? (3.142)86Chapter 3. Numerical Realizations: The Neumann Green's FunctionThis implies that a1 = 0;a2 = 1;a3 = a;bn = 0. Recall r(?) = 1 + "?(?). Thecurvature in polar coordinates is given by (3.86). Substituting the expressionfor r into the expression for the curvature we ?nd?(?) = 1 + "(3 cos2? + 8acos3?) + O(")2 (3.143)?0(?) = ?6"sin2? ? 24"asin3? (3.144)Using equation (3.141) we ?nd that?0(?) = ?4"??sin2? + 5a2 sin 3??(3.145)Using (3.143) and (3.144), we can ?nd conditions on the parameter a thatinvolves the minima and maxima of ? and ?. Note that at ? = n? , the ?rstderivative of ? and ? are zero. Thus, we have extrema at ? = n?. We calculatethe second derivatives when ? = ??00(?) = ?6"(2 ? 12a) (3.146)?00(?) = 4"??2 ? 15a2?(3.147)We see that at ? = ? , there is a maximum for ? when a < 16 and there is amaximum for ? when a < 415. Thus for a 2 (16; 415), ? has a maximum where ?has a minimum. We have the reverse result if ? is negative. Thus, the principleeigenvalue, ?0, given by (2.51), (2.89), (2.94) for one hole, N absorbing arcs andN identical absorbing arcs are not in general minimized at the maximum of thedomain curvature.Now at ? = 0, the second derivatives are?00(0) = ?12"(1 + 36a) (3.148)?00(0) = ?4"??2 + 18a2?(3.149)We see that ?00(0) and ?00(0) are both negative for a > 0. Thus, both ? and ?have maximums at ? = 0. In fact, it is where both attain their global maximums.We can relate these ?ndings to our boundary element method by plottingRm(x0;x0) versus x0(?0). The shape of the domain is shown in Figure 3.11 (a),while the curvature, ? ? 1, and ? ? C, where C is a constant of integration, areplotted in Figure 3.11 (b) for a = 0:2 and " = 0:1.87Chapter 3. Numerical Realizations: The Neumann Green's Function(a) Plot of the perturbed unit disk with boundary r = 1+"(cos2? + acos3?).(b) Plot of k(?) ? 1 (red) from (3.143) and ?(?) ? C (blue) from (3.145).Figure 3.11: Plot of the perturbed unit disk, the curvature and ? with " = 0:1and a = 0:2.We see from Figure 3.11 (a) that at ? = ? the curvature has a local minimum88Chapter 3. Numerical Realizations: The Neumann Green's Functionand ? has a local maximum. At ? = 0 and ? = 2?, we see that ? and ? haveglobal maxima here. Now we plot ? = Rm(x0;x0) from the numerical methodfor N = 600 and N = 2400 mesh points in Figure 3.12 with a = 0:2 and " = 0:1.Figure 3.12: Plot of Rm(x0;x0) for N = 600 (red) and N = 2400 (blue) with" = 0:1 and a = 0:2.We see that ? attains a local maximum at ? = ?, while at ?0 = 0, ? attainsa global maximum. Thus, our numerical results using the boundary elementmethod agree with the analytical predictions qualitatively. However, the scalingof our numerical results does not match the analytical results. We see fromconsidering the curve for N = 600 and N = 2400 that the numerical methodis converging extremely slowly since the curves are almost identical. To attainaccurate results that are quantitatively comparable to the analytical result werequire many more mesh points. We are not able to run the method for moremesh points due to computational limitations.3.5 DiscussionIn this chapter, we applied the results from Chapter 2 to calculate the MFPTin a few domains. We calculated the MFPT in the unit disk. In particular, we89Chapter 3. Numerical Realizations: The Neumann Green's Functionconsidered one hole on the boundary of the unit disk with the initial position atthe centre of the disk and the initial position at the antipodal point to the hole.Our results are comparable to that derived by Singer et al. in [23]. In addition,we investigate the behaviour of the MFPT as the initial position moved towardsthe hole at the antipodal point. These are depicted in Figure 3.1. With twoholes on the boundary, and the initial position at the centre, we found thatthe MFPT reached a minimum when the holes were furthest apart. This isdepicted in Figure 3.2. We also derived results for N equally spaced points onthe boundary of the unit disk. We found that v(x) ! 14D when N ! 1, whenlog 1" >> NWe derived the Neumann Green's function in the unit square in order tocalculate the MFPT. We considered the behaviour of the MFPT as the ini-tial position moves towards the hole, and the case where the initial position is?xed, but the hole moves. These are depicted in Figures 3.4 and 3.5. We alsoconsidered two holes on the boundary of the unit square. The behaviour inthis case was not as predictable as the other cases. We held one hole ?xed at(0;0:5), while the other hole moved from (0;0:48) to (1;0:5) along the boundary.We found that the minimum of the MFPT on the ?rst two sides were slightlyshifted as a result of the interaction between the two holes. In other words, theminimum was not where the distance between the initial position and movingexit window was a minimum. However, on the third side, this was the case as aresult of the symmetry of the con?guration. These results are shown in Figures3.6 and 3.7. It is worth noting that our expression for the MFPT is not validat the corners of the unit square. We derived an expression for the MFPT thatis valid at the corner x0 = (0;0) given by 3.55.In the last two sections of this Chapter, we considered an analytical and a nu-merical approach to solving for the Green's function in more arbitrary domains.The boundary element method was introduced as a numerical technique to solvefor the Green's function. We showed that the boundary element method didindeed give the correct results in the case of a unit disk. We then proceeded touse the boundary element method to ?nd the Neumann Green's function in anellipse with semi-major axis, a = 2 and semi-minor axis, b = 1. More precisely,we used our numerical technique to ?nd R(x;x0). We found that the numericalmethod converged for N = 640 mesh points. We found that R(x;x0) = 0:1291for all x within the   and for all x0 on the boundary. Thus, R(x;x0) is aconstant. We expect then that R(x0;x0) = 0:1291.To calculate R(x0;x0) we must rewrite integral (3.78) as we did in (3.88)90Chapter 3. Numerical Realizations: The Neumann Green's Functionand (3.90). The calculation of R(x0;x0) proved to be a challenge, and requiresmany mesh points for convergence. By considering Table 3.4 we see that thedi?erence between entries in the vector R decreases by approximately 0:57 asthe number of mesh points doubles. The convergence is thus extremely slow.We propose that the integral (3.78) be calculated a more accurate way in thiscase. We propose using gaussian quadrature in stead of the midpoint rule tocompute this integral. This is left for further work.We plotted the behaviour of the Neumann Green's function for an ellipsewith x at the origin and x0 at various points along the boundary in Figure3.8. We found that the Green's function satis?ed the symmetry propertiesthat we expected. We used the numerical results for the Green's function to?nd the MFPT for the ellipse with x at the origin and x0 moving along theboundary. The results are depicted in Figure 3.9. We then moved the initialposition to x = (1;0) with x0 moving along the boundary as before. The resultsare depicted in Figure 3.10. We see that the MFPT is a minimum when theexit window is closest to the initial position, and a maximum when the exit isfurthest from the initial position.To ?nd the Neumann Green's function in a perturbed circular domain weconsidered an analytical approach by T. Kolokolnikov [7]. Here, we derived anexpression for the regular part of the Green's function, Rm(x0;x0). One mustintegrate this expression and impose the integral condition to obtain a uniquesolution. We considered an example of a particular domain with ? given by(3.142), which is depicted in Figure 3.11 (a). We found that the regular partof the Green's function, ?, (3.145), has a local maximum while the curvature,?, (3.143) has a local minimum at ? = ? for a particular parameter regime,a 2 (16; 415). We plotted the the curvature and ? in Figure 3.11 (b) for a = 0:2and " = 0:1.We then calculated Rm(x0;x0) using our numerical boundary element methodand plotted this against ?0 in Figure 3.12 for N = 600 and N = 2400 meshpoints. By comparing these curves we see that the numerical method is con-verging extremely slowly since the two curves are similar. We see that our nu-merical results qualitatively agree with the analytical result in that R(x0;x0)has a maximum at ? = ?. However, the scaling is incorrect in comparison to theanalytical result. To obtain a more accurate numerical result we must calculatethe integral (3.78) more accurately.91Chapter 4The Narrow EscapeProblem in ThreeDimensions - The SphereWe now proceed to the narrow escape problem for the mean ?rst passage timefor a Brownian particle trapped within a three-dimensional domain. The math-ematical statement of the problem remains the same, however, the method ofapproach and results are di?erent from those in two dimensions. We will deriveresults for a sphere as the con?nement domain. It is an open problem to extendthis framework to arbitrary domains in three dimensions. The extension is notas straightforward as in two dimensions, as the local geometry and curvatureof the three-dimensional domain plays a signi?cant role in the ?nal solution. Itis, thus, necessary to solve for the Green's function for each three-dimensionaldomain before one can proceed to ?nd the mean ?rst passage time.4.1 Derivation of the Neumann Green'sfunction for a SphereWe begin by deriving the Neumann Green's function, Gs(x;x0), for a sphere,with a singularity on the boundary, which satis?es?Gs = 1j j; x 2  ; (4.1)@rGs = 1R2 ?(cos ' ? cos'0)?(? ? ?0); x 2 @ ; (4.2)Z Gsdx = 0: (4.3)92Chapter 4. The Narrow Escape Problem in Three Dimensions - The Spherewhere   = fx j jxj ? 1g and j j = 4?R33 . The angles are de?ned with 0 < ' < ?,the latitude, and 0 ? ? < 2?, the longitude. We have the property that ifg('0) = 0 and g(') is a monotone function, then ?(g(')) = ?('?'0)jg0('0)j . Thus?(cos ' ? cos '0) = ?(' ? '0)sin'0: (4.4)The singularity of the Green's function in three dimensions is a simple pole.Since the singular point is on the boundary, the singularity is twice as large.Thus, we expectGs(x;x0) = 12? jx ? x0j+ R(x;x0); (4.5)where R(x;x0) has a milder singularity than a simple pole as x ! 1. In thethree-dimensional case we show that R(x;x0) has a logarithmic singularity asx ! x0.Now we de?ne GpGp = 16j j?jxj2 + jx0j2?: (4.6)Note that jxj2 = r2 and jx0j2 = R2. To solve (4.1) and (4.2) we let Gs =Gp + ?s. Substituting this into (4.1) and (4.2) we ?nd? ?s = 0; x 2  ; (4.7)@r ?s = 1R2 ?(cos ' ? cos'0)?(? ? ?0) ? 14?R2 ; x 2 @ :(4.8)We will look for a solution written in terms of Legendre Polynomials. To thisend, we must rewrite the Neumann boundary condition in terms of Legendrepolynomials. This leads to the following Lemma:Lemma 4.1: We claim that?(cos ' ? cos '0)?(? ? ?0) = 14?1Xm=0(2m + 1)Pm(cos ?); (4.9)93Chapter 4. The Narrow Escape Problem in Three Dimensions - The Spherewherecos ? = x ? x0 = cos 'cos '0 + sin 'sin '0 cos(? ? ?0);x = (cos ?sin ';sin ?sin ';cos');x0 = (cos ?0 sin '0;sin ?0 sin '0;cos '0):Here x and x0 are written in spherical coordinates in the unit sphere, R = 1.Proof: We recall the completeness formula1Xm=0mXn=?mY currency1mn('0;?0)Ymn(';?) = ?(? ? ?0)?(cos ' ? cos'0); (4.10)where the Ymn are the spherical harmonics. The addition theorem for Legendrepolynomials states that(2m + 1)4? Pm(cos ?) =mXn=?mY currency1mn('0;?0)Ymn(';?): (4.11)Summing (4.11) from m = 0 to 1, we obtain the desired result (4.9).?Using this Lemma, the boundary condition (4.2) becomes@r ?s = 14?R21Xm=1(2m + 1)Pm(cos ?);onr = 1: (4.12)We look for a solution for ?s in the form?Gs =1Xm=1am? rR?mPm(cos ?); (4.13)which satis?es Laplace's equation by construction. Applying the boundary con-dition (4.12) we ?nd thatam = 14?R (2m + 1)m : (4.14)94Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereWe can write the solution for ?s as?Gs = 14?R1Xm=1(2m + 1)m? rR?mPm(cos ?) (4.15)= 12?R1Xm=1? rR?mPm(cos ?) + 14?R1Xm=11m? rR?mPm(cos?):(4.16)We deal with these two terms separately.We start with the ?rst term. Recall the generating function for Legendrepolynomials1p1 ? 2xt + t2 =1Xm=0Pm(x)tn: (4.17)Thus,12?R1Xm=1? rR?mPm(cos ?) = 12? 1r2 + R2 ? 2rRcos ? ? 12?R: (4.18)Now we consider the second term in (4.16). Let ? = rR andI = P1m=1 1m?mPm(cos ?). Using this we ?ndI0(?) = 1?1Xm=1?mPm(cos ?) = 1?"1p1 ? 2? cos ? + ?2 ? 1#(4.19)Also note, that I(0) = 0. We integrate the equation aboveI =Z ?0?1s1p1 ? 2scos ? + s2 ?1s!ds= log?21 ? ? cos ? + p1 + ?2 ? 2? cos ?!(4.20)Then with Gs = Gp + ?s, and substituting for ? we ?nd thatGs = 16j j?jxj2 + jx0j2?+ 12? 1jx ? x0j? 12?R+ 14?R log? 2RR ? r cos? + jx ? x0j?+ C: (4.21)95Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereWe are ready to use the integral condition (4.3), R2?0 R?0 RR0 Gs sin ?r2d?d?dr =0, to ?nd the constant C.We know thatlog? 2RR ? r cos? + jx ? x0j?=1Xm=11m? rR?mPm(cos ?): (4.22)We also have thatZ ?0Pm(cos ?)sin ?d? = 2sin m?m?(m + 1) =? 0 if m 6= 0;2 if m = 0: (4.23)Thus, the integral over ? of (4.22) is zero.Next we note that12? jx ? x0j =12?R1Xm=0Pm(cos?)? rR?m: (4.24)We can perform the integration over ? by using property (4.23). We ?nd that12?R1Xm=0Z 2?0Z ?0Z R0Pm(cos ?)? rR?msin?r2d?d?dr = 2R23 : (4.25)Using (4.22) and (4.25) we ?nd that C = ? 15?R for the integral condition, (4.3)to be satis?ed. Our ?nal expression for Gs isGs(x;x0) = 12? jx ? x0j+ 16j j?jxj2 + R2?+ 14?R log? 2RR ? jxj cos ? + jx ? x0j?? 710?R(4.26)We will need the expression for the Green's function, (4.26) in the limit thatx ! x0. We lety = x ? x0" ; currency1 = R ? r" : (4.27)We then calculate using the law of cosines thatR ? jxj cos ? = 12R?jx ? x0j2 ? (jxj2 ? R2)?? 12R ?O(")2 ? ((R ? "currency1)2 ? R2)? ? "currency1 + O(")2: (4.28)96Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereTherefore, by substituting (4.27) and (4.28) into (4.26), we obtain in the limitas x ! x0 thatGs ? 12?"jyj + 14?R log? 2R"(jyj + currency1)?? 920?R: (4.29)Now we calculate jx ? x0j in the limit as x ! x0 in spherical coordinates.We use a Taylor expansion about x0. We ?nd thatjx ? x0j =q(r ? R)2 + R2(' ? '0)2 + R2 sin2 '(? ? ?0)2 as x ! x0(4.30)We letcurrency1 = R ? r" ; s1 = Rsin '(? ? ?0)" ; s2 = R(' ? '0)" : (4.31)Using this, we ?nd thatjx ? x0j = "?currency12 + s21 + s22?1=2 as x ! x0 (4.32)This implies thatjyj = ?currency12 + s21 + s22?1=2 : (4.33)Combining (4.29) with (4.33) we obtain the far ?eld behaviour of the Neu-mann Green's function.4.2 Three-Dimensional Sphere with NAbsorbing Patches on the BoundaryWe want to ?nd the MFPT in a unit sphere,   = fx j jxj ? 1g. @ r is there?ecting part of the boundary, @ , while @ a is the absorbing part consistingof N non-overlapping circular patches de?ned by@ "j = f(';?) j (' ? 'j)2 + sin2 'j(? ? ?j)2 ? "2a2j = r2"g (4.34)Thus, the area of the circular patch centred at (1;'j;?j) is ??@ "j?? = ?"2a2j forj = 1;2;::N. Without loss of generality, there is no absorbing patch centred ata pole of the coordinate system. In Cartesian coordinates, the location of the97Chapter 4. The Narrow Escape Problem in Three Dimensions - The Spherecentre of the jth patch is given byxj = cos?j sin 'j; yj = sin ?j sin 'j; zj = cos 'j; (4.35)where jxjj = 1.The Laplacian in spherical coordinates is?u = 1r2 ?r2ur?r + 1r2 sin2 '@??u + 1r2 sin '@' (sin'@'u);= urr + 2rur + 1r2 sin2 'u?? + cot 'r2 u' + 1r2 u'': (4.36)We write the mean escape time in terms of eigenfunctions as in (2.6), as??0 + ?0?0 = 0; x 2  ;?0(x) = 0; x 2 @ a =N[j=1@ "j;@r?0 = 0; x 2 @ r: (4.37)4.2.1 Calculation of ?0 and ?0We need to solve the system (4.37) to ?nd the principal eigenvalue, ?0, and theprincipal eigenfunction, ?0.The eigenvalue is expanded as?0(") = "?1 + "2 log?"2??2 + "2?3 + ? ? ? : (4.38)In the outer region, away from the absorbing arcs, we expand the principaleigenfunction as?0(x;") = u0 + "u1 + "2 log?"2?u2 + "2u3 + ? ? ? : (4.39)Here, u0 = 1j j1, the solution to the leading order problem.Substituting these expansions into (4.37), we ?nd to order O(") that?u1 = ??1u0; x 2  ;@ru1 = 0; x 2 @ r;Z u1dx = 0: (4.40)98Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereThe integral condition enforces the normalization condition. At the next order,O ?"2 log ?"2??, we ?nd?u2 = ??2u0; x 2  ;@ru2 = 0; x 2 @ r;Z u2dx = 0: (4.41)At order O("2) we ?nd?u3 = ??1u1 ? ?3u0; x 2  ;@ru3 = 0; x 2 @ r;Z ?2u0u3 + u2?dx = 0: (4.42)In the absorbing region we introduce the following change of variables, (4.31),near the jth absorbing patchcurrency1 = 1 ? r" ; s1 = sin 'j(? ? ?j)" ; s2 = ' ? 'j" : (4.43)We letv(currency1;s1;s2;") = u?1 ? "currency1; "s1sin ' + ?j;"s2 + 'j?: (4.44)To proceed, we must transform the Laplacian, (4.36), into an expressioninvolving the inner variables using the chain rule. We ?nd that?w = "?2 (wcurrency1currency1 + ws1s1 + ws2s2) + "?1 (?2wcurrency1 + 2currency1 (ws1s1 + ws2s2))+ "?1 cot 'j (ws2 ? s22ws1s1) + O(1): (4.45)In the inner region we expand the eigenfunction asv = v0 + "log?"2?v1 + "v2 + ? ? ? : (4.46)We project all the patches onto a circular disk of radius a2j in the plane. At99Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereO("0) we ?ndLv0 = v0currency1currency1 + v0s1s1 + v0s2s2 = 0; currency1 ? 0; ?1 < s1;s2 < 1;@currency1v0 = 0; currency1 = 0; s21 + s22 ? a2j;v0 = 0; currency1 = 0; s21 + s22 ? a2j: (4.47)At O ?"log ?"2?? we ?ndLv1 = 0; currency1 ? 0; ?1 < s1;s2 < 1;@currency1v1 = 0; currency1 = 0; s21 + s22 ? a2j;v1 = 0; currency1 = 0; s21 + s22 ? a2j: (4.48)At order O(") we haveLv2 = 2 (v0currency1 ? currency1(v0s1s1 + v0s2s2)) ? cot 'j (v0s2 ? 2s2v0s1s1)@currency1v2 = 0; currency1 = 0; s21 + s22 ? a2j;v2 = 0; currency1 = 0; s21 + s22 ? a2j: (4.49)To proceed, we must match inner and outer solutions. We compare (4.39)to (4.46) in the limit as x ! xj or jyj = ?currency12 + s21 + s22?1=2 ! 1. We see thatv0 = u0 = j j?12 in this limit. We write the solution asv0 = j j?12 (1 ? vc); (4.50)where vc satis?esLvc = 0; currency1 ? 0; ?1 < s1;s2 < 1;@currency1vc = 0; currency1 = 0; s21 + s22 ? a2j;vc = 1; currency1 = 0; s21 + s22 ? a2j: (4.51)This is the well-known electri?ed disk problem in electrostatics, [6], which hasa solution given in terms of the Bessel Function J0(z) byvc = 2?Z 10sin ?? e??currency1=ajJ0???aj?d?; (4.52)100Chapter 4. The Narrow Escape Problem in Three Dimensions - The Spherewhere ? = ?s21 + s22?1=2. This solution has the asymptotic behaviourvc ? cjjyj + O?1jyj3!as jyj = ?currency12 + s21 + s22?1=2 ! 1: (4.53)Here, cj is the capacitance of a disk of radius aj in an in?nite plane given bycj = 2aj? : (4.54)The result above is obtained using Laplace's method and is shown in AppendixB.Thus v0 ? 1j j1?1 ? cjjyj?. Writing y in inner variables, as in (4.27) and(4.33), we ?nd from matching in the limit as x ! xj that1j j12+ "u1 + "2 log?"2?u2 + ? ? ? ? 1j j12? "cjj j1=2 jx ? xjj+ ? ? ? :This gives the missing singular condition on u1u1 ? ? cjj j12 jx ? xjjas x ! xj j = 1;::;N: (4.55)We must solve for u1 from (4.40) and (4.55). Recall, from (4.1) and (4.2), thata singularity on the boundary can be written as@ru1 = ?(' ? 'j)?(? ? ?j)sin 'jon r = 1: (4.56)This gives rise to a term in the solution of the form u1 ? 12?jx?xjj as x ! xj.Thus, we can write the problem for u1 as?u1 = ??1u0; x 2  ; (4.57)@ru1 = ? 2?j j12NXj=1cj ?(' ? 'j)?(? ? ?j)sin 'jon r = 1; (4.58)Z u1dx = 0: (4.59)Before we ?nd the solution for u1 we impose a solvability condition that willallow us to ?nd the eigenvalue ?1.101Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereWe integrate (4.57) over   and apply the Divergence Theorem to obtainZ ?u1dx =Z@ @u1@r jr=1 dS;??1u0 j j = ?NXj=12?cjj j12:Thus,?1 = 2?j j1Xj=1cj; cj = 2aj? : (4.60)We can write the solution for u1 asu1 = ? 2?j j121Xi=1ciGs(x;xi); (4.61)where Gs satis?es (4.1)-(4.3) with the Neumann boundary condition writtenusing (4.4). That is, Gs satis?es?Gs = 1j j; x 2  ; (4.62)@rGs = ?(' ? 'i)?(? ? ?i)sin 'ion r = 1; (4.63)Z Gsdx = 0: (4.64)We proceed by matching. In the limit as x ! xj, u1 becomesu1 ? ? 2?j j12[cjGs(x;xj)] ? 2?j j12NXi=1;i6=jciGs(xj;xi); (4.65)where the near-?eld behaviour of Gs(x;xj) as x ! xj is given by (4.29) withy given by (4.33).The outer expansion in terms of inner variables for x ! xj is1j j12? cjj j12 jyj+ 2?j j1=20@" cj4? log ?"2? + "cj4? log (jyj + currency1) + "9cj20? ? "NXi=1;i6=jciGs(xj;xi)1A+ "2 log?"2?u2 + "2u3 + ? ? ? : (4.66)102Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereThe inner expansion is1j j12? cjj j12 jyj+ "log?"2?v1 + "v2 + ? ? ? : (4.67)Comparing these expressions we ?nd thatv1 ? cj2j j12as jyj ! 1: (4.68)To proceed, we letv1 = cj2j j12(1 ? vc); (4.69)where vc satis?es (4.51). Using the behaviour of vc as jyj ! 1 given by (4.53),we ?nd thatv1 ? cj2j j12?1 ? cjjyj + O?1jyj3!!as jyj ! 1: (4.70)Substituting this into (4.67) and comparing to (4.66) we ?nd from matchingthatu2 ? ? c2j2j j12 jx ? xjjas x ! xj; j = 1;:::;N: (4.71)Thus, we must solve for u2 from (4.41) and (4.71). Proceeding as we didwhen ?nding u1, we rewrite the boundary condition of (4.41) as@ru2 = ? ?j j12NXj=1c2j ?(' ? 'j)?(? ? ?j)sin 'jon r = 1: (4.72)Before solving for u2 we use a solvability condition to ?nd ?2.We integrate (4.41) over   and use the Divergence Theorem. We ?nd that?2 = ?j jNXj=1c2j: (4.73)At this stage, as " ! 0, we have the following two-term asymptotic behaviour? ? 2?j j"NXj=1cj + ?j j"2 log?"2? NXj=1c2j + O("2); (4.74)where cj = 2aj? .103Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereThe solution for u2 is written in terms of the Green's function Gs asu2 = ? ?j j12NXi=1c2i Gs(x;xi); (4.75)where Gs satis?es (4.62)-(4.64). It is key to note that there is no need to expandu2 as x ! xj. This is a consequence of the fact that any unmatched terms inu2 are dominated by the O(") terms in the expansion (4.66).The O(") terms in u as x ! xj written in the variable y are? 2?j j12Aj + cj2j j12log(jyj + currency1); (4.76)whereAj = ? 9cj20? +NXi=1;i6=jciGs(xj;xi): (4.77)From (4.67) we see thatv2 ? ? 2?j j12Aj + cj2j j12log(jyj + currency1) as jyj ! 1: (4.78)Thus v2 must satisfy (4.49) and (4.78). We can rewrite the right-hand side ofthe ?rst equation of (4.49) by using v0s1s1 + v0s2s2 + v0currency1currency1 = 0 as2(currency1v0currency1)currency1 ? cot 'j (v0s2 ? 2s2v0s1s1): (4.79)We use an asymptotic decomposition to solve for v2 to obtainv2 = ? 2?j j12Aj(1 ? vc) + cj2j j12v2p; (4.80)wherev2p ? log (jyj + currency1) as jyj ! 1: (4.81)Here, vc is homogeneous solution, satisfying (4.51) with the behaviour as jyj !1, given by (4.53). That is vc ? cjjyj + O?1jyj3?as jyj ! 1. Furthermore,v2p is the particular solution satisfying (4.49) with the far-?eld behaviour givenby (4.81). By using the far-?eld behaviour of vc, (4.53), combined with a few104Chapter 4. The Narrow Escape Problem in Three Dimensions - The Spheresimple estimates, we ?nd that the far-?eld behaviour of v2 isv2 ? ? 2?j j12Aj?1 ? cjjyj?+ cj2j j12log(jyj + currency1); as jyj ! 1: (4.82)From matching, we must then conclude thatu3 ? 2?j j12Ajcjjx ? xjj as jyj ! 1: (4.83)This matches the unmatched term in the inner expansion v2 ? 2?j j1Ajcjjyj .The problem for u3 becomes?u3 = ??1u1 ? ?3u0; x 2  ; (4.84)@ru3 = 4?2j j12NXj=1Ajcj ?(' ? 'j)?(? ? ?j)sin 'jon r = 1; (4.85)Z u3dx = 0: (4.86)We can ?nd ?3 using a solvability condition. We integrate (4.84) over  , recall-ing that R  u1dx = 0. We use the Divergence Theorem to obtain?3 = ?4?2j jNXj=1Ajcj: (4.87)where Aj is given by (4.77).We have the following proposition:Proposition 4.2: (N Holes) In the limit as " ! 0, ?0 and ?0 have thefollowing three-term asymptotic behaviour?0(") ? 2?"j jNXj=1cj + "2 log?"2? ?j jNXj=1c2j ? 4?2"2j jNXj=1Ajcj; (4.88)?0(x;") = 1j j12? 2?j j12"NXi=1ciGs(x;xi) ? "2 log?"2? ?j j12NXi=1c2i Gs(x;xi)+ O("2): (4.89)105Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereWe can simplify the notation by introducing the Green's matrix, ?, andcapacitance vector C by? =0BBBBBBB@R Gs(1;2) Gs(1;3) ? ? ? Gs(1;N)Gs(2;1) R Gs(2;3) ? ? ? Gs(2;N)Gs(3;1) Gs(3;2) R ? ? ? Gs(3;N)... ... ... ... ...Gs(N;1) Gs(N;2) ? ? ? ? ? ? R1CCCCCCCA;C =?c1 c2 c3 ? ? ? cN?T;e =?1 1 1 ? ? ? 1?T:Note that the matrix ? is symmetric since Gs(j;i) = Gs(i;j), where Gs(j;i)denotes Gs(xj;xi). Also, we de?ne R byR = ? 920?: (4.90)Now, we can write ? as?0(") ? 2?"j j?CT e + "2 log?"2?CT C ? 2?"CT ?C?: (4.91)We can simplify this further by simplifying the expression for the Green's func-tion Gs given by (4.26). We know that jxij = jxjj = 1. Furthermore, by thelaw of cosines, 1 ? cos ?ij = jxi?xjj22 . Using this, the Green's function reducestoGs(xj;xi) = R+ 12?? 1jxj ? xij + log 2 ?12 log jxj ? xij ?12 log (2 + jxj ? xij)?:We de?ne Hij byHij =? 1jxj ? xij ?12 log jxj ? xij ?12 log (2 + jxj ? xij)?; (4.93)for i 6= j. Thus, ? can be written as?ij = R + log 22? + 12?Hij; (4.94)106Chapter 4. The Narrow Escape Problem in Three Dimensions - The Spherefor i 6= j. For i = j we see that ?ii = R.We de?ne p(x1;x2;:::;xN) byp(x1;x2;:::;xN) = CT ?C =NXi=1NXj=1cicj?ij;= RNXi=1c2i +NXi=1NXj=1;j6=icicj?ij: (4.95)We now write ? as?0(") ? 2?"j j?CT e + "2 log?"2?CT C ? 2?"p(x1;x2;:::;xN)?: (4.96)Now we consider the special case when each of the absorbing patches areof the same radius, that is aj = a. In this case, the capacitance, cj, becomesc = 2a? . We have the following proposition:Proposition 4.3: (N Identical Holes) In the limit as " ! 0;?0 and ?0have the three-term asymptotic behaviour?0(") ? 2?"Ncj j?1 + "2 log?"2?c?+ 2?"Ncj j0@"c0@9N10 ? (N ? 1)log 2 ? 1NNXi=1NXj=1;j6=iHij1A1A;(4.97)?0(x;") = 1j j12? 2?"cj j12NXi=1Gs(x;xi) ? "2 log?"2? ?c2j j12NXi=1Gs(x;xi)+ O("2); (4.98)where Hij is given by (4.93).For the case of equally sized absorbing patches, ?0("), is maximized withrespect to fx1;x2;:::;xNg when we ?nd fx1;x2;:::;xNg with jxjj = 1, for all jthat satis?es the following discrete variational problemHN = minfx1;x2;:::;xNgH(x1;x2;:::;xN) = minfx1;x2;:::;xNgNXi=1NXj=1;j6=iHij; (4.99)where Hij is given by (4.93). By using HN in (4.97) for ?0, we obtain a max-107Chapter 4. The Narrow Escape Problem in Three Dimensions - The Sphereimum for ?0 up to O("2). The ?rst term in H(x1;x2;:::;xN) is the usualCoulomb singularity in three-dimensions while the other two terms represent acontribution from surface di?usion on the boundary of the sphere.The discrete variational problem (4.99) is a modi?cation of the classicaldiscrete variational problems for ?nding optimal con?gurations of N points onthe surface of a sphere that minimize the Coulomb energyminfx1;x2;:::;xNgNXi=1NXj=1;j6=i1jxi ? xjj; jxij = 1; (4.100)or the logarithmic energy? minfx1;x2;:::;xNgNXi=1NXj=1;j6=ilog jxi ? xjj : (4.101)The problem of minimizing Hij appears to be a generalization of these twoclassic discrete variational problems. An overview of these variational problemsare given in [19] and references therein.4.2.2 Calculation of the Mean First Passage TimeWe calculate the MFPT in the case where we have N holes on the boundary ofa unit sphere, of radius aj. We use expansions (4.88) and (4.89) to compute theMFPT. Recall, that the MFPT is given byv(x) ? (1;?0)D?0?0; (4.102)where x is the initial position within the sphere. This assumes, that (?0;?0) = 1,which is satis?ed by (4.89) since the integrals of Gs over   are zero. Recall thatGs is the Green's function in a sphere of radius one.The MFPT for N holes of di?ering radius in terms of ?0 is given byv(x) = 1D?0(")?1 ? 2?"NXi=1ciGs(x;xi) ? "2 log?"2??NXi=1c2i Gs(x;xi) + O("2)!;(4.103)where ?0 is given by (4.88).By substituting for ?0 we have the following proposition for the MFPT withN holes of di?ering radius ai on the boundary of the unit sphere:108Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereProposition 4.4: (N Holes) For " ! 0, the three-term asymptotic be-haviour of the mean ?rst passage time, v(x), isv(x) = j j2?"D PNi=1 ci?1 + "2 log?2"? PNi=1 c2iPNi=1 ci? 2?"NXi=1ciG(x;xi)!+ j j2?"D PNi=1 ci?2?"PNi=1 cip(x1;x2;:::;xN) + O("2 log ")!; (4.104)where ci = 2ai? , p(x1;x2;:::;xN) is given by (4.95), Gs is given by (4.92) andj j = 4?3 . The average MFPT, ?, is? = j j2?"D PNi=1 ci?1 + "2 log?2"? PNi=1 c2iPNi=1 ci!+ j j2?"D PNi=1 ci?2?"PNi=1 cip(x1;x2;:::;xN) + O("2 log ")!:(4.105)In the case of N identical holes of radius a, so that ci = c = 2a? , we ?nd:Corollary 4.5: (N Identical Holes) For " ! 0, the three-term asymptoticbehaviour of the mean ?rst passage time isv(x) = j j2?"DNc?1 + "2 log?2"?c ? 2?"cNXi=1Gs(x;xi)!+ j j2?"DNc0@"c0@?9N10 + (N ? 1) log 2 + 1NNXi=1NXj=1;j6=iHij1A1A+ O ?"2 log "?; (4.106)where Hij is given by (4.93). The average MFPT, ?, is? = j j2?"DNc?1 + "2 log?2"?c?+ j j2?"DNc0@"c0@?9N10 + (N ? 1)log 2 + 1NNXi=1NXj=1;j6=iHij1A1A+ O ?"2 log "?: (4.107)109Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereNotice that if we minimize Hij as in (4.99), we minimize the MFPT. We can?nd an expression for one hole by setting N = 1 in (4.106) to obtain the follow-ing corollary:Corollary 4.6: (One Hole) For " ! 0, the three-term asymptotic be-haviour of the mean ?rst passage time, v(x), isv(x) = j j2?"Dc?1 + "2 log 2"c ? 2?"cGs(x;x0) ? 9"c10 + O("2 log ")?: (4.108)The average MFPT, ?, is? = j j2?"Dc?1 + "2 log 2" ? 9"c10 + O("2 log ")?: (4.109)In the case of one hole with a = 1, the MFPT from the centre of the sphere is:Corollary 4.7: (One Hole) For " ! 0, the three-term asymptotic be-haviour of the mean ?rst passage time from the centre of a unit sphere, v(x),with one absorbing patch of radius a = 1 on the boundary isE[? j x = (0;0)] = j j4"D?1 + "? log 2" ? 3"2? + O("2 log ")?: (4.110)Furthermore, for the case of one absorbing circular window of radius ", theaverage MFPT is ? ? j j4"D ?1 + "? log 1" + O(")?. This is comparable to the resultderived by Singer et al in [25] for the MFPT for one absorbing circular windowof radius ". They found the MFPT, and the average MFPT, to beE[? j x = (0;0)] = j j4"D?1 + "? log 1" + O(")?(4.111)The original result in equation (3.52) of [25] omits the ? term in (4.111) due toan omission of an extra factor of ? on the left-hand side of the equation above(3.52) in [25]. Our result, (4.110), agrees asymptotically with that of (4.111)and determines explicitly the O(") term to v(x).Our results, (4.104) and (4.106), generalizes the result in [25] to N circularabsorbing windows of di?erent radii on the unit sphere.110Chapter 4. The Narrow Escape Problem in Three Dimensions - The Sphere4.3 Three-Dimensional Sphere with NAbsorbing Patches on the Boundary - AnAlternate DerivationHere, we present an alternate derivation of the MFPT. Instead of writing theMFPT in terms of a complete set of eigenfunctions, we solve for the MFPTdirectly. That is, we want to solve?v(x) = ? 1D; x 2  ;v(x) = 0; x 2 @ a =N[j=1@ "j;@nv(x) = 0; x 2  r; (4.112)where @ "j is given by (4.34). The Laplacian is given (4.36). As before, thereis no absorbing patch centred at a pole of the coordinate system. In Carte-sian coordinates, the location of the centre of the jth is given (4.35). We willnow proceed to solve (4.112) directly using the method of matched asymptoticexpansions.4.3.1 Calculation of the MFPT DirectlyWe want to solve (4.112). In the outer region, away from absorbing patches, weexpand v(x) asv(x;") = v0" + v1 + "log?"2?+ "v3 + ? ? ? : (4.113)We obtain the following problems: At O("?1)?v0 = 0; x 2  ;@nv0 = 0; x 2  r; (4.114)at O("0)?v1 = ? 1D; x 2  ;@nv1 = 0; x 2  r; (4.115)111Chapter 4. The Narrow Escape Problem in Three Dimensions - The Sphereat O ?"log ?"2???v2 = 0; x 2  ;@nv2 = 0; x 2  r; (4.116)and at O(")?v3 = 0; x 2  ;@nv3 = 0; x 2  r: (4.117)In the inner region we introduce the change of variables given by (4.31) and(4.43). In the inner region we expandw(currency1;s1;s2;") = w0" + log?"2?w1 + w2 + ? ? ? : (4.118)We must transform the Laplacian, given by (4.36) into an expression in terms oflocal variables by using the chain rule. As before, we obtain (4.45). We obtainto O("?1)Lw0 = w0currency1currency1 + w0s1s1 + w0s2s2 = 0; currency1 ? 0;?1 < s1;s2 < 1;@currency1w0 = 0; currency1 = 0; s21 + s22 ? a2j;w0 = 0; currency1 = 0; s21 + s22 ? a2j: (4.119)At O ?log ?"2?? we haveLw1 = 0; currency1 ? 0; ?1 < s1;s2 < 1;@currency1w1 = 0; currency1 = 0 s21 + s22 ? a2j;w1 = 0; currency1 = 0 s21 + s22 ? a2j: (4.120)At O("0) we haveLw2 = 2 (currency1w0currency1currency1 + w0currency1) ? cot 'j (w0s2 ? 2s2w0s1s1);@currency1w2 = 0; currency1 = 0; s21 + s22 ? a2j;w2 = 0; currency1 = 0 s21 + s22 ? a2j: (4.121)We proceed by matching, which implies that w0 ! v0 as jyj ! ?currency12 + s21 + s22?1=2112Chapter 4. The Narrow Escape Problem in Three Dimensions - The Sphereor as x ! xj. We write the solution for w0 asw0 = v0(1 ? wc); (4.122)where wc satis?es (4.51). The solution to this is given by (4.52). The solutionhas the asymptotic behaviour (4.53) where cj is the capacitance of a disk ofradius aj in an in?nite plane given by (4.54). Writing y in outer variables asy = x?xj" and matching w0 to the outer solution we ?nd thatv1 ? ? v0cjjx ? xjjas x ! xj: (4.123)Thus, we must solve for v1 from (4.115), and (4.123).Recall from (4.1) and (4.2), that a singularity on the boundary can be writtenas@ru = ?(' ? 'j)?(? ? ?j)sin 'jon r = 1: (4.124)This gives rise to a term in the solution of the form u ? 12?jx?xjj as x ! xj.Thus, we can write the problem for v1 as?v1 = ? 1D x 2  ; (4.125)@rv1 = ?2?v0NXj=1cj ?(' ? 'j)?(? ? ?j)sin 'jon r = 1: (4.126)Before we ?nd v1, we apply a solvability condition that will allow us to ?ndv0. We apply the Divergence Theorem by integrating (4.125) and using theboundary condition (4.126). We ?nd thatv0 = j j2?D PNj=1 cj: (4.127)Now, the solution to v1 is written in terms of the surface Green's functionv1 = ?2?v0NXi=1ciGs(x;xi) + ?; (4.128)where Gs satis?es (4.62)-(4.64). Using this, notices that R  v1dx = ?j j so that113Chapter 4. The Narrow Escape Problem in Three Dimensions - The Sphere? = ?. We proceed by matching. In the limit as x ! xjv1 ? ?2?v0cjGs(x;xj) ? 2?v0NXi=1;i6=jciGs(xj;xi) + ?; (4.129)where the far-?eld behaviour of Gs is given by (4.29) along with (4.33). We let? = ?0 log?"2?+ ?1: (4.130)We write the outer expansion in terms of inner variables and match to the innersolution. That is, in the limit as jyj ! 1 we havev0" ?v0cj"jyj +v0cj2 log?"2?+ Bj + v0cj2 log (currency1 + jyj)+ ?0?"2?+ ?1 + "?"2?v2 + "v3+ ? ? ? ? v0"?1 ? cjjyj?+?"2?w1 + w2 + ? ? ? ; (4.131)whereBj = v00@9cj10 ? 2?NXi=1;i6=jciG(xj;xi)1A: (4.132)If we consider the expression for Aj, (4.77), we see that Bj = ?2?v0Aj.By matching we see thatw1 ? v0cj2 + ?0 as jyj ! 1: (4.133)The solution isw1 =?v0cj2 + ?0?(1 ? wc); (4.134)where wc satis?es (4.51). Using the far-?eld behaviour of wc given by (4.53) we?nd thatw1 ??v0cj2 + ?0??1 ? cjjyj?as jyj ! 1: (4.135)Substituting this into the matching condition (4.131) we ?nd thatv2 ? ??v0cj2 + ?0? cjjx ? xjj as x ! xj: (4.136)Thus, we must solve for v2 from (4.116) and (4.136). Proceeding as we did when114Chapter 4. The Narrow Escape Problem in Three Dimensions - The Sphere?nding v1, we rewrite the boundary condition of (4.116) as@rv1 = ?2?NXj=1?v0cj2 + ?0?cj ?(' ? 'j)?(? ? ?j)sin 'jon r = 1: (4.137)We use a solvability condition to ?nd ?0. Applying the Divergence Theoremto (4.116) and using the boundary condition (4.137) we ?nd that?0 = ?v0PNj=1 c2j2PNj=1 cj : (4.138)We can write the solution to v2 asv2 ? ?2?NXi=1?v0ci2 + ?0?ciGs(x;xi) + ?2: (4.139)At this stage, we have the following outer expansion for the MFPTv ? v0" + ?0 log?"2?+ ?1 ? 2?v0NXi=1ciGs(x;xi) + "log?"2?v2 + "v3; (4.140)where v0 is given by (4.127), ?0 is given by (4.138) and v2 is given by (4.139).The inner solution isw ? v0" (1 ? wc) + log?"2??v0cj2 + ?0?(1 ? wc) + w2 + ? ? ? ; (4.141)where wc satis?es (4.51).Now, if we consider the matching condition, (4.131), with (4.140) and (4.141)substituted, we ?nd that v2 contributes an unmatched term of O ?"2 log2 ?"2??.This term can be ignored since the O(1) unmatched terms dominate. Thus, wemust have thatw2 ? Bj + ?1 + v0cj2 log (currency1 + jyj) as jyj ! 1; (4.142)so that all the necessary terms are matched. Recall that w2 satis?es (4.121).We solve this problem for w2 by superposition. We letw2 ? (Bj + ?1) (1 ? wc) + v0cj2 w2p;: (4.143)115Chapter 4. The Narrow Escape Problem in Three Dimensions - The SphereHere wc is the homogeneous solution satisfying (4.51) with the behaviour asjyj ! 1, given by (4.53). That is wc ? cjjyj + O?1jyj3?as jyj ! 1. Further-more, w2p is the particular solution satisfying (4.121) with the far-?eld behaviourgiven by (4.81). By using the far-?eld behaviour of wc, (4.53), combined with afew simple estimates, we ?nd that the far-?eld behaviour of w2 isw2 ? (Bj + ?1)?1 ? cjjyj?+ v0cj2 log(jyj + currency1): (4.144)From matching we must then conclude thatv3 ? ?cj (Bj + ?1)jx ? xjjas x ! xj: (4.145)We must solve for v3 from (4.117) and (4.145). We rewrite the boundary con-dition of (4.117) as@rv3 = ?2?NXj=1cj (Bj + ?1) ?(' ? 'j)?(? ? ?j)sin 'jon r = 1: (4.146)We apply the Divergence Theorem to (4.117) and use the boundary condition(4.146). We ?nd that?1 = ?PNj=1 cjBjPNj=1 cj; (4.147)where Bj is given by (4.132). We let R = ? 920? as in (4.90). Substituting forv0 from (4.127) we ?nd that?1 = 2?v0PNj=1 cj0@NXj=10@c2jR + cjNXi=1;i6=jciGs(xj;xi)1A1A: (4.148)Substituting for ?1 and v2 into (4.140), we ?nd that the outer solution forthe MFPT is116Chapter 4. The Narrow Escape Problem in Three Dimensions - The Spherev(x) = j j2?"D PNj=1 cj0@1 ? "2 log ?"2??PNj=1 c2j?PNj=1 cj? 2?"NXi=1ciG(x;xi)1A+ j j2?"D PNj=1 cj0@ 2?"PNj=1 cj0@NXj=10@c2jR + cjNXi=1;i6=jciGs(xj;xi)1A1A1A+ O("2 log "): (4.149)Notice that?PNj=1?c2jR + cj PNi=1;i6=j ciGs(xj;xi)??is the same expression asp(x1;x2;:::;xN) given by (4.95). Thus, we have recovered the result that wasderived in the previous section for the MFPT with N holes on the boundary,Proposition 4.4, (4.104). Thus, by simplifying (4.149) in the same manner as wedid (4.104), we recover Corollary 4.5, 4.6 and 4.7 for N identical holes, one holeand one hole of radius a = 1. Thus, the two methods for deriving the MFPTagree exactly.Lastly, we compare the asymptotic results with the numerical results solving(4.112) computed using COMSOL [15] by R. Straube [28]. We use the expressionfor the average MFPT for N identical holes, (4.107). We set a = 1, so that wehave N identical circular windows of radius " equidistantly placed on the surfaceof the sphere. We let ?2 be the two-term asymptotic result obtained be omittingthe O(") terms in (4:107), ?3 be the full three-term asymptotic result (4.107)and ?n the full numerical result. We tabulate the results for N = 1;2 and 4 inTable 4.1 below.N = 1 N = 2 N = 4" ?2 ?3 ?n ?2 ?3 ?n ?2 ?3 ?n0:02 53:89 53:29 52:81 26:95 26:40 26:12 13:47 13:10 12:990:05 22:17 21:57 21:35 11:09 10:54 10:43 5:54 5:17 5:120:10 11:47 10:87 10:78 5:74 5:19 5:14 2:87 2:50 2:470:20 6:00 5:40 5:36 3:00 2:45 2:44 1:50 1:13 1:130:50 2:56 1:95 1:96 1:28 0:73 0:70 0:64 0:27 0:30Table 4.1: Comparison of ?2, ?3 and ?n for various values of ", N = 1;2 and 4.We see good agreement between the numerical results and the asymptoticresults, especially when using the full three-term asymptotic result, ?3.117Chapter 4. The Narrow Escape Problem in Three Dimensions - The Sphere4.4 DiscussionIn this Chapter we aimed to ?nd the MFPT in a unit sphere. To this end, webegan with a derivation of the Neumann Green's function for a sphere of radiusR. As in two dimensions, we transformed the problem for the MFPT into aneigenvalue problem. We solved for the principle eigenvalue, ?0, and the principleeigenfunction, ?0 for N holes of di?ering radius, using the method of matchedasymptotic expansions. These results are given by (4.88) and (4.89). We thensimpli?ed the expressions for ?0 and ?0 for the case of N identical holes, given by(4.97) and (4.98). We then found the MFPT in a unit sphere with N absorbingpatches of di?ering radii on the boundary given by (4.104). We found a simplerexpression in the case of N identical holes, given by (4.106). In the case ofone hole in a unit sphere we were able to compare our result, (4.110), to thatderived by Singer et al. in [25], (4.111). Our result agrees asymptotically with(4.111) and determines explicitly the O(") term. Our result, (4.104) generalizesthe result obtained in [25]. Furthermore, the result contains more terms in theasymptotic expansion for the MFPT, which depend on the initial position x,the location of the patches, xj, j = 1;::;N, and the radii of the holes.Next, we derived the MFPT directly using the method of matched asymp-totic expansions. We found that the result for the MFPT derived directly,(4.149), agrees exactly with (4.104). In addition, we compared our asymptoticresult for the average MFPT for identical holes, (4.107), with numerical resultssolving (4.112). We considered N equidistantly placed identical holes of radius" for N = 1;2 and 4. The numerical results were computed using COMSOL [15]by R. Straube [28]. We found good agreement with the numerical results whenusing the three-term expression for the average MFPT, (4.107). The results aretabulated in Table 4.1.It is worth noting that in the case of N circular absorbing windows of com-mon radius ", that the MFPT is minimized in the limit as " ! 0 at the con-?guration fx1;x2;:::;xNg that minimizes the discrete sum PNi=1 PNj=1;j6=i Hij.This also maximizes the eigenvalue ?0, thus minimizing the MFPT.Recall that our expressions for the MFPT are valid in the outer region andare not valid in the vicinity of the absorbing windows. That is, our results forthe MFPT hold when jx ? xjj >> O("), for each absorbing window j.118Chapter 5ConclusionIn this thesis we focused on the mean ?rst passage time, which is the mean timea Brownian particle takes to escape from a domain   2 Rd, d = 2;3, whoseboundary is re?ecting, @ r, except for a small absorbing window, @ a. We let" = j@ aj. This is known as the narrow escape problem, a singular perturbationproblem. As we have discussed, the mean ?rst passage time has applicationsin many ?elds including electrostatics, biology, and ?nance to name a few. Wediscussed speci?c applications in biology and relevant results found in [25], [5],[23], [24]. We compared our results with the work of these authors.As it was shown in Chapter 1, the MFPT is de?ned in terms of a conditionalexpectation, (1.10), dependent on the initial position x, which satis?es a mixedboundary value problem for the Poisson equation, (1.26)-(1.28). The focus ofthis thesis was to ?nd the MFPT in various two-dimensional domains and in aunit sphere.In Chapter 2, we transformed the problem for the MFPT, (2.3)-(2.5), intoan eigenvalue problem by writing the MFPT in terms of a complete set of eigen-functions. We found that the MFPT to leading order is inversely proportionalto the principal eigenvalue in the limit that " ! 0. This agrees with the lit-erature. Recall that the absorbing arc @ a is a perturbation that is small inextent but large in magnitude in a localized region. Thus, our methodology wasto ?nd the principal eigenvalue and eigenfunction, ?0 and ?0, using the methodof matched asymptotic expansions, a singular perturbation technique.We found the two-term asymptotic behaviour of the principal eigenvalue andeigenfunction, in the limit that " ! 0. We then proceeded to ?nd the MFPT.We determined the two-term asymptotic behaviour of the MFPT in an arbitrarytwo-dimensional domain with one hole on the boundary in the limit as " ! 0,given by (2.53). The leading order term agrees with that derived by Singer et.al in [25]. Our result provides more information since it is dependent on theinitial position, x, the location of the hole, x0, the length of the hole through dand the shape of the domain.119Chapter 5. ConclusionWe extended our results to include N absorbing holes on the boundary. Wederived results for the principal eigenvalue and eigenfunction with N holes onthe boundary and N identical holes on the boundary. In the case of N identicalholes, we derived another result that exploited the special property of a cyclicmatrix in the matrix system (2.85). We proceeded to derive results for theMFPT in these cases, given by (2.99), (2.101) and (2.103). We found in thecase of N identical holes that the leading order term was proportional to 1N ,which agrees with a result found by Holcman and Schuss in [5].With these results in hand, we proceeded to investigate the MFPT in speci?cgeometries, namely the unit disk and unit square in Chapter 3. We comparedour results for a unit disk with one hole on the boundary with those derived bySinger et al. in [23]. We found good agreement with their results. Furthermore,we investigated the behaviour of the MFPT as the initial position moves towardsthe hole on the boundary. We found that the MFPT decreases in this case. Wealso investigated two holes on the boundary, and found that the MFPT reaches aminimum when the holes are furthest apart. These results are shown in Figures3.1 and 3.2. With N symmetrically located holes on the boundary of the unitdisk, we found a special result. In the limit as N ! 1, v(x) ! 14D when x isplaced at the centre of the disk, with log 1" >> N.We derived the Green's function for a unit square. We used this result to?nd the MFPT. In the case where the initial position moves towards the hole,the MFPT decreases. When the initial position is ?xed and the hole movesalong one side of the unit square, we ?nd that the MFPT is a minimum wherethe distance between the hole and the initial position is a minimum. Theseresults are depicted in Figures 3.4 and 3.5. We then considered two holes onthe boundary, one ?xed, and one moved. The results are illustrated in Figure3.6 and 3.7. We found that the MFPT is not a minimum where the distancebetween the initial position and the moving arc is a minimum along sides 1and 2. This is a result of the interaction with the second arc which is ?xedat the midpoint of side 1. However, along side 3, the MFPT is a minimumwhen the two holes are furthest apart, and the distance from the moving arcand the initial position is a minimum. This is a consequence of the symmetryof the con?guration with the moving exit window along side 3. Notice thatour expression for the MFPT is not uniformly valid at the corners of the unitsquare, since the boundary is non-smooth at these points. By considering theexpression for the Green's function at the corner, we see that the singularityat the corner is twice as large as a singularity along the smooth portion of the120Chapter 5. Conclusionboundary. We derived an expression for the MFPT that is valid at the cornerx0 = (0;0) given by (3.55).As we have seen from the general expressions for v(x), the Neumann Green'sfunction plays an important role. To investigate di?erent geometries, we must?nd the Neumann Green's function. We proposed a numerical technique, theboundary element method, to ?nd the Neumann Green's function. We showedthat the method gave accurate results in the case of the unit circle. We thenattempted to use the method to ?nd the Green's function in an ellipse. Wefound that R(x;x0) = 0:1291 for all x 2   and x0 2 @ . Thus, we concludedthat R was a constant. This implies that R(x0;x0) = 0:1291. We plotted thebehaviour of the Green's function for x at the origin and x0 moving along theboundary in Figure 3.8. We attempted to ?nd R(x0;x0) numerically. However,the numerical method did not converge in this case. We require many moremesh points to attain the desired accuracy. At the present time, we do nothave the computational ability to reach the desired accuracy of four decimalplaces. We propose using an integration technique that converges quicker, orRichardson's extrapolation to attain a more accurate answer.In the case of perturbed circular domains, we used an analytical approach,derived by T. Kolokolnikov [7], to ?nd the Green's function. The result, givenby (3.141), is unique up to a constant. To ?nd the constant, we must imposethe integral condition R  G(x;x0)dx = 0. This analytic approach allows fora comparison with the numerical approach. We considered an example of aparticular domain with ? given by (3.142), which is depicted in Figure 3.11 (a).We found that the regular part of the Green's function, ?, (3.145), has a localmaximum where the curvature, ?, (3.143), has a local minimum at ? = ? for aparticular parameter regime, a 2 (16; 415). We plotted the the curvature and ?in Figure 3.11 (b) for a = 0:2 and " = 0:1. We then used our numerical methodto plot Rm(x0;x0) versus ?0 in Figure 3.12, for the shape given by (3.142),for N = 600 and N = 2400 mesh points along the boundary. By comparingthese curves we see that the numerical method is converging extremely slowlysince the two curves are similar. The numerical results qualitatively agree withthe analytical result in that Rm(x0;x0) has a maximum at ? = ?. However,the scaling is incorrect in comparison to the analytical result. To obtain amore accurate numerical result, or at least one that converges faster, we mustcalculate the integral (3.78) more accurately instead of using the midpoint rule.In three-dimensions, the statement of the narrow escape problem is no dif-ferent from the statement in two dimensions. Thus, the MFPT must satisfy121Chapter 5. Conclusion(1.26)-(1.28). As in two-dimensions, we transform the problem into an eigen-value problem, for the principal eigenvalue and eigenfunction. To proceed by themethod of matched asymptotic expansions, we solved for the Green's function?rst. This is di?erent to the approach in two dimensions. In three-dimensions,one must ?nd the Green's function ?rst, before ?nding the MFPT. We foundthe Neumann Green's function for a sphere of radius R. We found the principaleigenvalue and eigenfunction in a unit sphere with N holes on the boundarygiven by (4.88) and (4.89). We derived simpli?ed results in the case of N iden-tical holes. We found the MFPT for N absorbing patches of di?ering radiion the boundary, (4.104). We simpli?ed this result for the case of N identicalpatches. The result is given by (4.106). For one hole on the boundary of radiusa the MFPT is given by (4.108). In the case of one hole with radius a = 1 ina unit sphere we were able to compare our result, (4.110), to that derived bySinger et al. in [25], (4.111). Our result agrees asymptotically with (4.111) anddetermines explicitly the O(") term. Our result, (4.104), generalizes the resultobtained in [25]. Furthermore, the result contains more terms in the asymptoticexpansion for the MFPT, which depend on the initial position x, the locationof the patches, xj, j = 1;::;N, and the radii of the holes.Next, we derived the MFPT directly using the method of matched asymp-totic expansions. We found that the result for the MFPT derived directly,(4.149), agrees exactly with (4.104). We also compared our asymptotic resultfor the average MFPT for N identical holes of radius " on the boundary, (4.107),with numerical results solving (4.112). These numerical results were computedfor N = 1;2 and 4 identical holes placed equidistantly on the surface of thesphere using COMSOL [15] by R. Straube [28]. We found good agreement withthe numerical results when comparing them to the three-term asymptotic result(4.107). The results are tabulated in Table 4.1.It is worth noting that in the case of N circular absorbing windows of com-mon radius ", that the MFPT is minimized in the limit as " ! 0 at the con-?guration fx1;x2;:::;xNg that minimizes the discrete sum PNi=1 PNj=1;j6=i Hij.This also maximizes the eigenvalue ?0, thus minimizing the MFPT.We note that our expressions for the MFPT are valid in the outer regionand are not valid in the vicinity of the absorbing windows. That is, our resultsfor the MFPT hold when jx ? xjj >> O("), for each absorbing window j.We would like to extend the framework introduced in Chapter 4 to includean arbitrary three-dimensional domain. It is an open problem to ?nd the MFPTand the Green's function in an arbitrary three-dimensional domain.122Chapter 5. ConclusionWe would like to improve the convergence of the boundary element methodintroduced in Chapter 3 when ?nding R(x0;x0). As mentioned, we could cal-culate the integrals more accurately using gaussian quadrature. Alternatively,we could ?nd the error term for the technique and attempt Richardson extrapo-lation to ?nd more accurate answers for R(x0;x0). Lastly, we could change theelements for the method by using line segments instead of circular arcs. Withimproved convergence, we could use the technique to ?nd the Neumann Green'sfunction for domains with smooth boundaries more accurately.Finally, the related problem for the Dwell time would be an extension of thework presented in this thesis.123Bibliography[1] P.C. Bresslo? and B.A. Earnshaw. Di?usion-trapping model of receptortra?cking in dendrites. Physical Review E, 75:041915{1 { 041915{7, 2007.[2] P.C. Bresslo?, B.A. Earnshaw, and M.J. Ward. Di?usion of protein recep-tors on a cylindrical dendritic membrane with partially absorbing traps.SIAM Journal on Applied Mathematics, 68(5):1223{1246, 2008.[3] R.L. Burden and J. Douglas Faires. Numerical Analysis. Brooks/ColePublishing Company, 7th edition, 2001.[4] I.V. Grigoriev, Y.A. Makhnovskii, A.M. Berezhkovskii, and V.Y. Zitser-man. Kinetics of escape through a small hole. Journal of Chemical Physics,116(22):9574{9577, 2002.[5] D. Holcman and Z. Schuss. Escape through a small opening: Recep-tor tra?cking in a synaptic membrane. Journal Of Statistical Physics,117(516):975{1014, 2004.[6] J.D. Jackson. Classical Electrodynamics. Wiley, New York, 2nd edition,1945.[7] T. Kolokolnikov. Regular part of the modi?ed green's function on theboundary of a nearly circular domain. Unpublished Notes.[8] T. Kolokolnikov, M.S. Titcombe, and M.J. Ward. Optimizing the funda-mental Neumann eigenvalue for the Laplacian in a domain with small traps.European Journal of Applied Mathematics, 16(2):161{200, 2005.[9] T. Kolokolnikov and M.J. Ward. Reduced wave Green's functions and theire?ect on the dynamics of a spike for the Gierer-Meinhardt model. EuropeanJournal of Applied Mathematics, 14(5):513{546, 2003.124Bibliography[10] T. Kolokolnikov and M.J. Ward. Bifurcation of spike equilibria in thenear-shadow Gierer-Meinhardt mode. Discrete and Continuous DynamicalSystems { Series B, 4(4):1033{1064, 2004.[11] E. Kreyszig. Advanced Engineering Mathematics. John Wiley & Sons, NewYork, 6th edition, 1988.[12] R.E. Larson, R.P. Hostetler, and B.H. Edwards. Calculus of a Single Vari-able. Houghton Mi?in Company, 6th edition, 1998.[13] J.J. Linderman and D.A. Lau?enburger. Analysis of intracellular recep-tor/ligand sorting. calculation of mean surface and bulk di?usion timeswithin a sphere. Journal of BioPhysics, 50(2):295{305, 1986.[14] R.C. McCann, R.D. Hazlett, and D.K. Babu. Highly accurate approxi-mations of Green's and Neumann functions on rectangular domains. InMathematical, Physical and Engineering Sciences, volume 457, pages 767{772. The Royal Society, April 2001.[15] COMSOL Multiphysics (v3.4). COMSOL AB, Stockholm.[16] I.M. Nemenman and A.S. Silbergleit. Explicit Green's function of a bound-ary value problem for a sphere and trapped ?ux analysis in gravity probeB experiment. Journal of Applied Physics, 86(1):614{624, 1999.[17] J.W.S. Baron Rayleigh. The Theory of Sound, volume 2. Dover, New York,2nd edition, 1945.[18] S. Redner. A Guide to First Passage Processes. Cambridge UniversityPress, 2001.[19] E.B. Sa? and A.B.J. Kuijlaars. Distributing many points on a sphere.Mathematical Intelligencer, 19(1):5{11, 1997.[20] Z. Schuss. Theory and Applications of Stochastic Di?erential Equations.Wiley Series in Probability and Statistics. John Wiley & Sons, Inc., NewYork, 1980.[21] Z. Schuss, A. Singer, and D. Holcman. The narrow escape problem fordi?usion in cellular microdomains. PNAS, 104(41):16098{16103, 2007.125Bibliography[22] A.S. Silbergleit, I. Mandel, and I.M. Nemenman. Potential and ?eldsingularity at a surface point charge. Journal of Mathematical Physics,44(10):4460{4466, 2003.[23] A. Singer, Z. Schuss, and D. Holcman. Narrow escape, part II: The circulardisk. Journal of Statistical Physics, 122(3):465{489, 2006.[24] A. Singer, Z. Schuss, and D. Holcman. Narrow escape, part III: Non-smoothdomains and Riemann surfaces. Journal of Statistical Physics, 122(3):491{509, 2006.[25] A. Singer, Z. Schuss, D. Holcman, and R.S. Eisenberg. Narrow escape, partI. Journal of Statistical Physics, 122(3):437{463, 2006.[26] I. Stakgold. Green's Functions and Boundary Value Problems. John Wiley& Sons, Inc., 2nd edition, 1998.[27] J. Stewart. Multivariable Calculus. Brooks/Cole Publishing Company, 4thedition, 1999.[28] R. Straube. Personal Communication.[29] R. Straube, M.J. Ward, and M. Falcke. Reaction rate of small di?us-ing molecules on a cylindrical membrane. Journal of Statistical Physics,129(2):377{405, 2007.[30] A. Ta?ia and D. Holcman. Dwell time of a Brownian molecule in a mi-crodomain with traps and a small hole on the boundary. Journal of Chem-ical Physics, 126(23):234107{1 { 234107{11, 2007.[31] M.S. Titcombe and M.J. Ward. An asymptotic study of oxygen transportfrom multiple capillaries to skeletal muscle tissue. SIAM Journal of AppliedMathematics, 60(5):1767{1788, 2000.[32] S. Torquato, I.C. Kim, and D. Cule. E?ective conductivity, dielectricconstant, and di?usion coe?cient of digitized composite media via ?rst-passage-time equations. Journal of Applied Physics, 85(3):1560{1571, 1999.[33] M.J. Ward, W.D. Henshaw, and J.B. Keller. Summing logarithmic ex-pansions for singularly perturbed eigenvalue problems. SIAM Journal onApplied Mathematics, 53(3):799{828, 1993.126Bibliography[34] M.J. Ward and J.B. Keller. Strong localized perturbations of eigenvalueproblems. SIAM Journal on Applied Mathematics, 53(3):770{798, 1993.127Appendix ADerivation of dWe want to solve for vc(y), where vc satis?es?yvc = 0; y = @ 0;vc = 0; y 2 @ 0;vc ! log jyj as jyj ! 1: (A.1)The problem for vc(y), (2.28), has a unique solution and the behaviour at 1 is[8]vc(y) ? log jyj ? log d + p:yjyj2 : (A.2)where d, the logarithmic capacitance and p = (p1;p2), the dipole vector, aredetermined from the length of the hole, @ a. In our case, where the length ofthe hole is j@ aj = 2", d = 12. This is obtained by a special solution in terms ofelliptic cylindrical coordinates. We present a derivation here.In elliptic cylindrical polar coordinates y = (x;y) transforms asx = Rcosh ? cos currency1; y = Rsinh? sincurrency1; (A.3)with 0 ? currency1 < 2?. Furthermore we ?nd for a function f(x;y) that fxx + fyy =f?? + fcurrency1currency1 = 0.Under these coordinates we have thatx2R2 cosh2 ? +y2R2 sinh2 ? = 1: (A.4)Thus, curves of constant ? form ellipses. We chooseRcosh ?0 = a; Rsin ?0 = b: (A.5)128Appendix A. Derivation of dTherefore,R = (a2 ? b2)1=2; ?0 = 12 log?1 + b=a1 ? b=a?for a > b: (A.6)We have mapped the ellipse onto an in?nite sheet in the currency1 ? ? plane. To solvethe problem (A.1) we must solvef?? + fcurrency1currency1 = 0; ? > ?0; 0 ? currency1 < 2?;f = 0 on ? = ?0;f ? log jrj as ? ! 1; (A.7)where r = px2 + y2 and f is 2? periodic in currency1. The solution is f = C(? ? ?0),where C is a constant.As ? ! 1 the coordinates (A.3) becomex = Re?2 coscurrency1; y =Re?2 sincurrency1: (A.8)In this limit r = Re?2 and therefore ? = log r + log 2R. We ?nd the behaviour off as ? ! 1 by substituting for ? in this limit and ?0 and R from (A.6). We?nd thatf = C?log r + log 2R ? 12 log?1 + b=a1 ? b=a??;= C?log r ? log?a + b2??: (A.9)Thus, to obtain f ? log r at 1, C = 1. By setting d = a+b2 we obtain thebehaviour (A.2) at 1. In our case b = 0, thus for an arc of length 2a, d = a2 . Ingeneral, for an arc of length l, d = l4. The solution for vc is given by C(? ? ?0),where one must map back to Cartesian coordinates. The far-?eld behavior ofvc is given by (A.2).129Appendix BFar-Field Behaviour of ?cWe consider?ccurrency1currency1 + ?cs1s1 + ?cs2s2 = 0;currency1 ? 0; ?1 < s1;s2 < 1; (B.1)@currency1?c = 0; currency1 = 0; s21 + s22 ? a2j; (B.2)?c = 1; currency1 = 0; s21 + s22 ? a2j: (B.3)The exact solution from [6] is?c = 2?Z 10sin ?? e??currency1=ajJ0???aj?d?; (B.4)where ? = (s21 +s22)1=2. Now, we want to derive the asymptotics of this solutionas jyj = ?currency12 + s21 + s22?1=2 ! 1.We introduce r = ??aj . Then?c = 2aj??Z 10sin(ajr=?)ajr=? e??rJ0(r)dr; (B.5)where ? = currency1=?. We would like to determine the asymptotic behavoiur as? ! 1, but with ? ?xed. As x ! 0 we know that sin xx ? 1 ? x26 + ? ? ? . Thus,?c ? 2aj???I0 ? a2j6?2 I1 + ? ? ?!; (B.6)whereI0 =Z 10e??rJ0(r)dr; (B.7)I1 =Z 10r2e??rJ0(r)dr: (B.8)130Appendix B. Far-Field Behaviour of ?cIt follows that I0 is the Laplace transform of J0(r). HenceI0 = 11 + ?2 : (B.9)Also I1 = d2d?2 I0(?).We ?nd that?c ? 2aj?"I0? ?a2j6?3 I1 + ? ? ?#; (B.10)whereI0? =1pcurrency12 + ?2 ; (B.11)I1?3 = ?1(currency12 + ?2)3=2+ 3currency12(currency12 + ?2)5=2: (B.12)De?ning jyj = ?currency12 + s21 + s22?1=2 we ?nd that?c ? 2aj??1jyj +a2j6?1jyj3 ?3currency12jyj5!!as jyj ! 1: (B.13)This is the asymptotic behaviour at 1 uniform in currency1;s1 and s2.131Appendix CMatlab CodeHere we give the code for the boundary element method for an ellipse. We?rst give the code for ?nding R(x;x0) and then the code for ?nding R(x0;x0).Lastly we give the code for ?nding the integral of R(x;x0) over  .clear all;N=640; dth = 2*pi/N; th = 0:dth:2*pi;s = zeros(N,1); %arclength vectorm = 5000; %number of steps in the trapezoidal rule forthe arclengtha = 2; b = 1; sum2 = 0; area = pi*a*b; for j = 1:Nlowerth = th(j);upperth = th(j+1);ht = (upperth-lowerth)/m;sum2 = 0;for k = 0:mthk = lowerth + k*ht;if k == 0sum2 = ht/2*sqrt((a^6*b^2*(sin(thk))^2+a^2*b^6*(cos(thk))^2)/((a^2*(sin(thk))^2 +b^2*(cos(thk))^2)^3)) + sum2;elseif k == msum2 = ht/2*sqrt((a^6*b^2*(sin(thk))^2+a^2*b^6*(cos(thk))^2)/((a^2*(sin(thk))^2 +b^2*(cos(thk))^2)^3)) + sum2;elsesum2 = ht*sqrt((a^6*b^2*(sin(thk))^2+a^2*b^6*(cos(thk))^2)/((a^2*(sin(thk))^2 +b^2*(cos(thk))^2)^3)) + sum2;end;end;132Appendix C. Matlab Codes(j) = sum2;end;% s is the vector of lengths of each segment% Setting up the midpointsfor j = 1:Nnth(j) = (0.5*(th(j)+th(j+1))); %midpointsend; r = (a*b)./(a^2.*sin(nth).^2 + b^2.*cos(nth).^2).^(1/2); xm=r.*cos(nth); ym = r.*sin(nth);plot(xm,ym,'*r'); %plot of meshpointsx = 1; y = 0.5;c = -a*b*(a^2-b^2); d = a^2.*sin(nth).^2 + b^2.*cos(nth).^2;rdr = c*sin(nth).*cos(nth).*(a^2.*sin(nth).^2+b^2.*cos(nth).^2).^(-3/2); %first derivativerddr = c.*((cos(nth).^2-sin(nth).^2).*(d.^(-3/2))-3*(a^2-b^2).*(sin(nth).^2).*(cos(nth).^2).*d.^(-5/2)); %second derivativek = (r.^2+2*rdr.^2-r.*rddr)./((r.^2+ rdr.^2).^(3/2)); %curvaturenx = (1./sqrt(rdr.^2+r.^2)).*(xm+rdr.*sin(nth)); ny =(1./sqrt(rdr.^2+r.^2)).*(ym-rdr.*cos(nth));f = inline('sqrt((x1-x2)^2+(y1-y2)^2)');for j = 1:Nbb(j) = 0;for i = 1:Nif i == jA(i,j) = 1-(s(i)*k(i))/(4*pi);bb(j) = s(i)/(4*pi*area)*log(f(x,xm(i),y,ym(i)))*(xm(i)*nx(i)+ym(i)*ny(i))+bb(j);elseA(i,j) = -s(i)/(2*pi*f(xm(i),xm(j),ym(i),ym(j))^2)*(nx(i)*(xm(i)-xm(j))+ny(i)*(ym(i)-ym(j)));bb(j) = s(i)/(4*pi*area)*log(f(x,xm(i),y,ym(i)))*(xm(i)*nx(i)+ym(i)*ny(i))+bb(j);end;end;end; R = A'\bb';133Appendix C. Matlab CodeNow we give the code for the calculation of R(x0;x0).clear all; clfN=2400; dth = 2*pi/N; th = 0:dth:2*pi; s = zeros(N,1); m = 5000;a = 2; b = 1; area = pi*a*b; sum2 = 0; delta = 0; for j = 1:Nlowerth = th(j);upperth = th(j+1);ht = (upperth-lowerth)/m;sum2 = 0;for k = 0:mthk = lowerth + k*ht;if k == 0sum2 = ht/2*sqrt((a^6*b^2*(sin(thk))^2+a^2*b^6*(cos(thk))^2)/((a^2*(sin(thk))^2 +b^2*(cos(thk))^2)^3)) + sum2;elseif k == msum2 = ht/2*sqrt((a^6*b^2*(sin(thk))^2+a^2*b^6*(cos(thk))^2)/((a^2*(sin(thk))^2 +b^2*(cos(thk))^2)^3)) + sum2;elsesum2 = ht*sqrt((a^6*b^2*(sin(thk))^2+a^2*b^6*(cos(thk))^2)/((a^2*(sin(thk))^2 +b^2*(cos(thk))^2)^3)) + sum2;end;end;s(j) = sum2;end;% s is the vector of lengths of each segmentfor j = 1:Nnth(j) = (0.5*(th(j)+th(j+1))); %midpointsend; r = (a*b)./(a^2.*sin(nth).^2 + b^2.*cos(nth).^2).^(1/2); xm=r.*cos(nth); ym = r.*sin(nth); plot(xm,ym,'.'); c =-a*b*(a^2-b^2); d = a^2.*sin(nth).^2 + b^2.*cos(nth).^2; rdr =c*sin(nth).*cos(nth).*(a^2.*sin(nth).^2+b^2.*cos(nth).^2).^(-3/2);rddr = c.*((cos(nth).^2-sin(nth).^2).*d.^(-3/2)-3*(a^2-b^2).*(sin(nth).^2).*(cos(nth).^2).*d.^(-5/2));134Appendix C. Matlab Codek = (r.^2+2*rdr.^2-r.*rddr)./((r.^2+ rdr.^2).^(3/2)); %curvaturenx = (1./sqrt(rdr.^2+r.^2)).*(xm+rdr.*sin(nth)); ny =(1./sqrt(rdr.^2+r.^2)).*(ym-rdr.*cos(nth));f = inline('sqrt((x1-x2)^2+(y1-y2)^2)'); Ns = N; for ks = 1:Nsx = xm(ks);y = ym(ks);for j = 1:Nbb(j) = 0;for i = 1:Nif i == jA(j,i) = 1-(s(i)*k(i))/(4*pi);if (x == xm(i))& (y == ym(i))c = 0;m = (k(i)*s(i)-delta)/2;bb(j) = 1/(8*pi*area*k(i)^2)*(s(i)*k(i)*log(2/k(i)^2)+2*(s(i)*k(i))*log(s(i)*k(i)/2) - s(i)*k(i)*log(2)-2*s(i)*k(i)+m*(log(2*(1-cos(c+m/sqrt(3)))/((c+m/sqrt(3))^2))+log(2*(1-cos(c-m/sqrt(3)))/((c-m/sqrt(3))^2)))) + bb(j);else%THIS INTEGRAL IS CALCULATED IN THE GLOBAL COORD%SYSTEM WITH THE MIDPOINT RULEbb(j) = s(i)/(4*pi*area)*0.5*log((x-xm(i))^2+(y-ym(i))^2)*(xm(i)*nx(i)+ym(i)*ny(i))+bb(j);end;else%EXPRESSION USING GLOBAL COORDINATE SYSTEM AND MIDPOINT%RULEA(j,i) = -s(i)/(2*pi*((xm(i)-xm(j))^2+(ym(i)-ym(j))^2))*(nx(i)*(xm(i)-xm(j))+ny(i)*(ym(i)-ym(j)));if (x == xm(i))& (y == ym(i))c=0;m = (k(i)*s(i)-delta)/2;bb(j) = 1/(8*pi*area*k(i)^2)*(s(i)*k(i)*log(2/k(i)^2135Appendix C. Matlab Code)+2*(s(i)*k(i))*log(s(i)*k(i)/2) - s(i)*k(i)*log(2)-2*s(i)*k(i)+m*(log(2*(1-cos(c+m/sqrt(3)))/((c+m/sqrt(3))^2))+log(2*(1-cos(c-m/sqrt(3)))/((c-m/sqrt(3))^2)))) + bb(j);else%THIS INTEGRAL IS CALCULATED IN THE GLOBAL COORD%SYSTEM WITH THE MIDPOINT RULEbb(j) = s(i)/(4*pi*area)*0.5*log((x-xm(i))^2+(y-ym(i))^2)*(xm(i)*nx(i)+ym(i)*ny(i))+bb(j);end;end;end;end; R = A\bb'; finalR(1,ks) = R(1,1);%finalR(:,ks) = R;end; %ksmin = 1; max = 0; for k = 1:Nif finalR(1,k) < minmin = finalR(1,k);end;if finalR(1,k) > maxmax = finalR(1,k);end;end;Finally, we give the code for calculating R  R(x0;x0)dx using the midpointrule.clear all;%changed to use the midpoint ruleN=80; %x_0's, the integral should be independent of thism=5000; %integration steps in theta for finding the arclengthdth = 2*pi/N; th = 0:dth:2*pi; s = zeros(N,1); a = 2; b = 1; sum2= 0; area = pi*a*b;for j = 1:Nlowerth = th(j);upperth = th(j+1);ht = (upperth-lowerth)/m;136Appendix C. Matlab Codesum2 = 0;for k = 0:mthk = lowerth + k*ht;if k == 0sum2 = ht/2*sqrt((a^6*b^2*(sin(thk))^2+a^2*b^6*(cos(thk))^2)/((a^2*(sin(thk))^2 +b^2*(cos(thk))^2)^3)) + sum2;elseif k == msum2 = ht/2*sqrt((a^6*b^2*(sin(thk))^2+a^2*b^6*(cos(thk))^2)/((a^2*(sin(thk))^2 +b^2*(cos(thk))^2)^3)) + sum2;elsesum2 = ht*sqrt((a^6*b^2*(sin(thk))^2+a^2*b^6*(cos(thk))^2)/((a^2*(sin(thk))^2 +b^2*(cos(thk))^2)^3)) + sum2;end;end;s(j) = sum2;end;% s is the vector of lengths of each segmentfor j = 1:Nnth(j) = (0.5*(th(j)+th(j+1))); %midpoints/singular pointsend; r = (a*b)./(a^2.*sin(nth).^2 + b^2.*cos(nth).^2).^(1/2); xm=r.*cos(nth); ym = r.*sin(nth); plot(xm,ym,'*'); hold on; c =-a*b*(a^2-b^2); d = a^2.*sin(nth).^2 + b^2.*cos(nth).^2; rdr =c*sin(nth).*cos(nth).*(a^2.*sin(nth).^2+b^2.*cos(nth).^2).^(-3/2);rddr = c.*((cos(nth).^2-sin(nth).^2).*d.^(-3/2)-3*(a^2-b^2).*(sin(nth).^2).*(cos(nth).^2).*d.^(-5/2));kr = (r.^2+2*rdr.^2-r.*rddr)./((r.^2+ rdr.^2).^(3/2)); %curvaturenx = (1./sqrt(rdr.^2+r.^2)).*(xm+rdr.*sin(nth)); ny =(1./sqrt(rdr.^2+r.^2)).*(ym-rdr.*cos(nth));%this sets up the integration variables over xn = 100; ht = (2*pi)/n; k = 1:1:n; tk = k.*ht; uppert =(a*b)./((a^2.*(sin(tk)).^2 + b^2.*(cos(tk)).^2).^(1/2));hr =(uppert); f = inline('sqrt((x1-x2)^2+(y1-y2)^2)');%we have separated everything out into integrals137Appendix C. Matlab Codefor k = 1:nthk = tk(k);h = hr(k); % this is the length of the interval, b-ax = 0.5*h*cos(thk);y = 0.5*h*sin(thk);plot(x,y,'.r');hold on;for j = 1:Nbb(j) = 0;for i = 1:Nif i == jA(i,j) = 1-(s(i)*kr(i))/(4*pi);bb(j) = s(i)/(4*pi*area)*log(f(x,xm(i),y,ym(i)))*(xm(i)*nx(i)+ym(i)*ny(i))+bb(j);elseA(i,j) = -s(i)/(2*pi*f(xm(i),xm(j),ym(i),ym(j))^2)*(nx(i)*(xm(i)-xm(j))+ny(i)*(ym(i)-ym(j)));bb(j) = s(i)/(4*pi*area)*log(f(x,xm(i),y,ym(i)))*(xm(i)*nx(i)+ym(i)*ny(i))+bb(j);end;end;end;R(:,k) = A'\bb';end;for ks = 1:N %this is the loop over the x0,y0's%numerical integrationsum3 = 0; sumt = 0; sumw = 0;for k = 1:n %(t index)h = hr(k);thk = tk(k);sumt = sumt + 0.5*ht*0.5*h^2*log((0.5*h*cos(thk)-xm(ks))^2+(0.5*h*sin(thk)-ym(ks))^2);sum2 = 0.5*h^2*R(ks,k); %r index, midpoint rulesum3 = sum3 + ht*sum2;138Appendix C. Matlab Codesumw = sumw + 0.5*ht*h^2*(h^2/4+xm(ks)^2+ym(ks)^2);end;sumfinal(ks) = sum3;sum2final(ks) = sumt;sum3final(ks) = sumw;end; %for ks139

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
China 10 21
United States 5 1
United Kingdom 4 0
France 3 0
Canada 2 0
Poland 2 0
Palau 2 0
City Views Downloads
Unknown 11 1
Beijing 9 4
Ashburn 2 0
Redmond 2 0
Kitchener 1 0
Vancouver 1 0
Shenzhen 1 17
Oxford 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}
Download Stats

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.24.1-0066535/manifest

Comment

Related Items