Asymptotic Analysis of First Passage Processes With Applications to Animal Movement by Venu Kurella B.Tech., The Indian Institute of Technology, Guwahati, 2009 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in The Faculty of Graduate Studies (Mathematics) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) August 2011 ➞ Venu Kurella 2011 Abstract Understanding the dependence of animal behaviour on resource distribution is a central problem in mathematical ecology. In a habitat, the distribution of food resources and their accessibility from an animal’s location together with the search time involved in foraging, all govern the survival of a species. In this work, we investigate various scenarios that affect foraging habits of animals in a landscape. The work, unlike previous studies, analyzes the first passage quantities on complex prey-predator distributions in a given domain in order to derive simple analytical problems that can readily be solved numerically. We use standard stochastic models such as the Kolmogorov equations of first passage times and splitting probability, to model both the foraging time of a predator and the chances of survival of prey on a landscape with prey and predator patches. We obtain an asymptotic solution to these Kolmogorov equations using a hybrid asymptotic-numerical singular perturbation technique that utilizes the fact that the ratio of the size of prey patches is small in comparison to the overall landscape. Results from this hybrid approach are then verified by undertaking full numerical simulations of the governing partial differential equations of the first passage processes. By using this hybrid formulation we identify the underlying parameters that affect the search time of a predator and splitting probability of prey, which are otherwise difficult to ascertain using only numerical tools. This analytical understanding of how parameters influence the first passage processes is a key step in quantifying foraging behavior in model ecological systems. ii Table of Contents Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii Table of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii List of Tables vi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . ix 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Analysis of animal movement . . . . . . . . . . . . . . . . . . 1 1.2 Stochastic modelling of animal movement . . . . . . . . . . . 1 1.3 Need for asymptotic methods . . . . . . . . . . . . . . . . . . 2 1.4 Application of asymptotic methods . . . . . . . . . . . . . . 3 1.5 Mean first passage time in ecology . . . . . . . . . . . . . . . 4 1.6 Asymptotic analysis of MFPT in ecology . . . . . . . . . . . 4 1.7 Splitting probability in ecology . . . . . . . . . . . . . . . . . 4 1.8 Numerical methods implemented . . . . . . . . . . . . . . . . 5 1.9 Brief discussion of chapters . . . . . . . . . . . . . . . . . . . 5 2 Asymptotic Analysis of the Mean First Passage Time . . 7 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Asymptotic solution for the MFPT without drift . . . . . . . 9 2.3 Numerical verification for the MFPT without drift . . . . . . 17 2.3.1 The case of one prey patch . . . . . . . . . . . . . . . 18 2.3.2 The case of several prey patches . . . . . . . . . . . . 21 iii Table of Contents 2.4 Asymptotic solution for the MFPT with drift . . . . . . . . . 22 2.5 Numerical verification for the MFPT with drift . . . . . . . . 26 2.5.1 27 The cases of one and several prey patches . . . . . . . 2.6 Comparison of the cases with and without drift . . . . . . . 2.7 Numerical approximation (regularization scheme) of Green’s 30 function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 Distribution of the predator in the domain . . . . . . . . . . 37 3 Asymptotic Calculation of the Splitting Probability . . . . 40 2.8 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.2 Asymptotic analysis . . . . . . . . . . . . . . . . . . . . . . . 42 3.3 Numerical verification . . . . . . . . . . . . . . . . . . . . . . 46 3.4 Results for rings . . . . . . . . . . . . . . . . . . . . . . . . . 47 3.4.1 . . . . . . . . . . . . . . . . 49 3.5 Effect of distance from the prey patch . . . . . . . . . . . . . 50 3.6 The effect of closely spaced patches . . . . . . . . . . . . . . 52 3.6.1 Two closely spaced patches . . . . . . . . . . . . . . . 52 3.6.2 Two close patches and a distant predator . . . . . . . 61 Results with larger rings 4 Mean First Passage Time with Variable Diffusivity . . . . 64 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.2 Asymptotic solution for one patch . . . . . . . . . . . . . . . 65 4.3 Solution for multiple patches . . . . . . . . . . . . . . . . . . 67 4.4 Numerical verification . . . . . . . . . . . . . . . . . . . . . . 68 4.4.1 One patch case . . . . . . . . . . . . . . . . . . . . . . 69 4.4.2 Multiple patches case . . . . . . . . . . . . . . . . . . 70 5 Second Moment Analysis and Estimation of Variance . . . 72 5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 5.2 Asymptotic solution . . . . . . . . . . . . . . . . . . . . . . . 73 5.3 Regularization scheme of the function, W0pk . . . . . . . . . 77 5.4 Numerical verification . . . . . . . . . . . . . . . . . . . . . . 79 5.4.1 A single patch in a circular domain . . . . . . . . . . 80 5.4.2 Multiple patches in a circular domain . . . . . . . . . 81 iv Table of Contents 6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 Appendix Numerical Methods . . . . . . . . . . . . . . . . . . . . . . . . . . 89 v List of Tables 2.1 The logarithmic capacitance, or shape-dependent parameter, dj , for some cross-sectional shapes of Ωj = ε−1 Ωεj . . . . . . . . 2.2 Spatial averages of MPFT in a circular domain of radius 1km without drift containing three prey patches . . . . . . . . . . 2.3 48 Spatial averages of the splitting probability when the target is located inside/outside the four predator ring . . . . . . . . 3.3 29 Spatial averages of the splitting probability when the target is located inside/outside the four predator ring . . . . . . . . 3.2 21 Distribution averages of MPFT in a domain with drift with three patches . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 12 48 Spatial averages of splitting probability when the target is located inside/outside an eight-predator ring . . . . . . . . . 48 3.4 Spatial averages of splitting probability with large patches. . 49 4.1 Spatial averages of the MPFT in a domain without drift with three patches with variable diffusivity . . . . . . . . . . . . . 5.1 Spatial averages of SMPFT in a domain without drift with three patches . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 71 81 Spatial averages of VFPT in a domain without drift with three patches . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 vi List of Figures 2.1 Domain of animal movement for the MFPT with localized prey patches . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9 2.2 Inner region for MPFT . . . . . . . . . . . . . . . . . . . . . . 10 2.3 SP Model domain for numerics . . . . . . . . . . . . . . . . . 17 2.4 Full numerical result of MFPT with one patch and no drift . 18 2.5 MFPT with one patch and no drift . . . . . . . . . . . . . . . 19 2.6 Spatial average of distance to patch vs patch location . . . . 20 2.7 Full numerical result of MFPT with three patches and no drift 21 2.8 Full numerical result of MFPT with one patch and with drift 27 2.9 MFPT with one patch and with drift . . . . . . . . . . . . . . 28 2.10 Full numerical result of MFPT with three patches and with drift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 2.11 Comparision of MFPT with and without drift . . . . . . . . . 31 2.12 Comparision of MFPT with and without drift . . . . . . . . . 32 2.13 Probability distribution . . . . . . . . . . . . . . . . . . . . . 39 3.1 Domain of animal movement for SP . . . . . . . . . . . . . . 41 3.2 Inner region for SP . . . . . . . . . . . . . . . . . . . . . . . . 43 3.3 SP Model domain for numerics . . . . . . . . . . . . . . . . . 46 3.4 Full numerical results for the splitting probability for a ring of four predators . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 47 Full numerical results for the splitting probability for a ring of eight predators . . . . . . . . . . . . . . . . . . . . . . . . . 49 3.6 Effect of distance on splitting probablity . . . . . . . . . . . . 51 3.7 SP: Conformal mapping for two circles . . . . . . . . . . . . . 54 3.8 Concentric circles in w-space . . . . . . . . . . . . . . . . . . 56 vii List of Figures 3.9 SP: Numerical solution for close patches . . . . . . . . . . . . 59 3.10 Plots of constant solutions as a function of radii and the distance between the closely spaced patches . . . . . . . . . . . . 60 3.11 SP: Numerical solution for close patches and a distant predator 63 4.1 MFPT with one patch and variable diffusivity D . . . . . . . 70 5.1 SMFPT with one patch and no drift . . . . . . . . . . . . . . 80 5.2 Spatial avgerage of distance to patch . . . . . . . . . . . . . . 81 viii Acknowledgements I would like to thank my supervisors Daniel Coombs and Michael Ward for such a rewarding and enjoyable experience in the course of this project. Thanks to Dan for introducing me to applied mathematics and for all the guidance and mentoring during these two years. Thanks to Michael for such a splendid experience in asymptotic analysis and for all those rewarding conversations during the course of the project. I really enjoyed working under their supervision. I would like to thank Neil Balmforth for taking time to examine this work and for his thorough comments on the thesis. I also thank all my friends and staff at IAM who have been supportive for the past two years and for sharing some memorable moments with me. Special thanks to my friend and colleague, Mike for his patient discussions which were helpful towards the completion of the project. I also thank all my friends and office staff in the mathematics department for making my stay unforgettable. Finally, I thank my family for being a constant source of love, support and encouragement so far in my life. ix Chapter 1 Introduction 1.1 Analysis of animal movement Predator-prey population dynamics are highly affected by the distribution of resources in the shared habitat. For example, the distribution of predator dens across the landscape affects prey survival, while prey foraging sites can impact the distribution of predators. Furthermore, predators hunting for prey have physical limitations (energy and various anatomical constraints) that lead to limitations on their speed as they hunt for food. So time and distance covered before successful hunting are factors that determine the fitness and survival of a predator. Similarly, prey foraging for food should avoid predator threats in order to survive and reproduce. This limits the distance prey can travel to safely find food. The time it takes for a predator to locate a prey, and the travel distance of predators involved in hunting are directly related to the location of the prey in a habitat. Similarly, the chances of prey avoiding predator patches are also related to the location of the predator in that habitat. So studying animal movement in a given landscape, with a particular prey-predator distribution is essential for building our understanding of the feeding patterns of animals, and to comment on their chances of survival or extinction. 1.2 Stochastic modelling of animal movement Like any other biological or ecological system, the dynamics of prey-predator systems are complicated to understand and ecologists depend heavily on modelling to gain insights into the system. Much of the data used for modelling comes from recording the GPS tracks of tagged animals. Studying the 1 1.3. Need for asymptotic methods GPS tracks has shown that animal movements have inherent stochasticity that cannot be ignored. So many stochastic models have been put in place to understand, analyse and predict animal movement patterns. In one of the earliest works, Skellam showed that diffusion can be used a good approximation for random walk problems (cf. [17]). Skellam used a diffusive approximation for spatial distribution and dispersal of species in one and two dimensional habitats. Stochastic methods to analyze animal movement have been later devised and improved over years. Recent models accommodate various internal factors in a landscape like resource disparities (e.g. prey rich or predator rich regions), spatial inhomogeneities like mountains, rivers, and human disturbances such as electricity transmission line cuts. A different set of models have been developed to understand the behaviour of denning animals which move only in home range (a confined region around the home where the animal is found to move most of the time). Currently pioneering mathematical ecology research is being undertaken by Professor Mark Lewis at the University of Alberta, by applying the techniques of stochastic modelling to explain the various GPS datasets generated from tagging animals in different landscapes. This thesis is partly inspired by this work. 1.3 Need for asymptotic methods Though stochastic modelling has been in place for a few decades, it has not been sufficiently robust to predict patterns. The reason for this difficulty is partly due to our inability to identify and estimate the parameters that affect the data and partly due do a lack of computational tools. Most of the stochastic formulations, including those discussed in this thesis, are finally expressed as Kolmogorov forward or Fokker- Planck equations, which are partial differential equations (PDEs) that must be solved on two-dimensional domains representing the landscape, and under various constraints. As technology has developed, there has been a growth in computational facilities to solve the PDEs numerically but numerical methods are cumbersome when it comes to analyzing problems with fine spatial structure, and are unsuitable for assessing the effects of parameter variations. Therefore, analyti2 1.4. Application of asymptotic methods cal approaches are extremely useful, allowing us to cope with parametric variation and understand their particular roles in affecting predicted animal behaviour in a landscape. However, finding analytical solutions to the governing PDEs is a daunting task for many realistic problems in ecology. When such problems have cropped up in other fields, researchers have directed their focus towards using asymptotic methods, and such methods have proven their worth over the decades. They have also been in use to verify numerical results and thus validate the results so obtained. The application of asymptotic methods has been most fruitful in solving problems related to engineering and physical systems rather than biological and ecological problems. In this thesis we aim to bring asymptotic techniques to bear on particular ecological problems which contain tightly localized spatial structures (e.g. individual animal dens in a large landscape), in an effort to supplement the field of mathematical biology with this technique, and ultimately, to deal with current challenges in foraging theory and ecology. 1.4 Application of asymptotic methods Crucially for our purposes, Ward et al. have developed a hybrid method to tackle singular perturbation problems in two dimensional domains(cf. [21]). The crux of the method lies in summing the higher order logarithmic terms with excellent accuracy. This simple approach has unleashed the power and accuracy of asymptotic analysis on PDEs in two-dimensional domains with localised traps. The benefits of the method were immediately reaped, when it was applied to biological problems like oxygen transport from capillaries to skeletal muscle and cellular signal transduction through localised traps on the cell surface (cf. [20], [4]). However, the application of this robust technique to ecological problems is novel and is the main theme of this thesis. In this work, we mainly focus on using asymptotic approaches to solve Kolmogorov equations with localized traps, for mean first passage time and splitting probability problems. 3 1.5. Mean first passage time in ecology 1.5 Mean first passage time in ecology The mean first passage time (MFPT) is defined as the mean time taken by a dynamic particle to reach a defined target for the first time. This concept finds application in wide range of fields. From stock exchange to neuron firing it is used in every field that has some level of inherent stochasticity (cf. [15]). The concept of MFPT has been in use in biology in processes like electroporation to calculate the time of mechanical break down(cf. [6]) and switching of flagellar motion in bacteria (cf. [1]). The application of MFPT analysis to ecology and in particular to understanding animal movement across the landscape is an emerging field that is yet to be fully explored (cf. [12], [7]). 1.6 Asymptotic analysis of MFPT in ecology Exact analytical calculations of the MFPT can be obtained by solving a particular Kolmogorov equation, but this is usually tedious and often impossible. But so far asymptotic analysis has not been applied to estimate MFPT in the context of ecological systems and it is our goal to explore the field and aid in improvement. Part of this thesis also shows how we can achieve improved accuracy over recent analytical works of Benichou et al., in solving the MFPT equations on two-dimensional domains (cf. [2]). 1.7 Splitting probability in ecology Splitting probability, by definition, is the probability that a dynamic particle reaches one particular target patch before reaching another, unwanted regions. For example, if we want to calculate the chances of a single prey reaching its den safely, in a landscape with predator patches, we need to calculate a splitting probability. Similar problems of calculating splitting probability in two-dimensional domains with small patch targets has been addressed by Benichou et al. but the solution has a limited accuracy (cf. [2]). Here in this work we apply hybrid asymptotic techniques to estimate the 4 1.8. Numerical methods implemented splitting probability of a prey reaching its home den in a domain with predator patches, to a higher level of accuracy. This application is novel and yields excellent results in the context of animal movement. The application of the concept of splitting probability to animal movement is, to our knowledge, completely novel. Also it is a concept well-tailored to investigate animal movement in the presence of favourable and unfavourable localities across the landscape. 1.8 Numerical methods implemented Along with asymptotic methods, we also show that multi-patch problems can be effectively handled by numerical methods in MATLAB PDE Toolbox. All through the analysis, we used the parameters estimated numerically by McKenzie et al. using animal movement GPS data (cf. [12]). Our numerical and analytical results are used to verify each other. 1.9 Brief discussion of chapters Chapter 2 Here we first introduce the singular perturbation method and use the approach to solve Kolmogorov equation for MFPT to a particular target patch in the landscape, in cases with and without a defined drift term that continuously drives the otherwise diffusing predator towards its home den. We also introduce the numerical methods used in thesis. We reproduce the results similar to those by McKenzie et al.(cf. [12]) for one target patch and then analyse the problem for more than one target patch. Chapter 3 We calculate splitting probabilities for a prey to return to its home den while not reaching predator patches, for various types of predator patch distributions across the landscape. We provide explicit solutions for an interesting case with a ring of predators surrounding the home den. We also 5 1.9. Brief discussion of chapters use different asymptotic techniques to investigate special cases where two patches are located extremely close to each other. Chapter 4 This chapter deals with the situation where the rate of animal motion is different in different locations (i.e. variable diffusivity across space). This situation is of considerable ecological interest as animal motion can certainly be retarded or accelerated by different conditions on the ground, including situations where human activity has modified the landscape. We perform analytical and numerical calculations for such a problem with one or more than one patch, confirming that our asymptotic methods have a lot of promise for analyzing this kind of situation. Chapter 5 In the chapters so far, we have calculated only the mean of the first passage time. In this chapter we extend our analysis to calculate also the variance of the first passage time, for problems where the MFPT was estimated in chapter 2. Chapter 6 In this last chapter we have a brief description of the conclusions drawn from the work and discussion of some interesting future problems. The ultimate goal of this work is to introduce the application of the powerful technique of hybrid-asymptotic method in potential areas in the field of mathematical ecology and motivate its implementation in more ecological problems in future. 6 Chapter 2 Asymptotic Analysis of the Mean First Passage Time 2.1 Introduction In a habitat, animal behaviour is dependent on resource distribution. The distribution of food resources can afffect the growth and survival of a species. The location of prey and the time required to search and hunt the prey are key aspects that affect the population dynamics of a predator. The distribution of prey in a particular habitat is subject to change due to many natural and human causes. The search time of a predator can be interpreted in terms of the mean first passage time (MFPT), which is the mean time taken by a predator to reach a specific patch (with prey) for the first time, starting from a given location in some two-dimensional spatial landscape. More specifically, the MFPT, T (X) is the mean time taken for the predator to reach its target for the first time, starting from the location, X. Earlier investigations of the MFPT with ecological applications relied heavily on a full numerical approach to solve the underlying PDE for the MFPT (cf. [12]). Here we use singular perturbation techniques to analytically derive an algebraic system that determines the MFPT in terms of a certain Green’s function, which in general must be computed numerically. This overall hybrid asymptotic-numerical approach has the benefit of eliminating the difficulty with resolving small spatial scales in a full numerical treatment of the PDE. Our semi-analytical approach also gives considerable insight into the dependence of the MFPT on the system parameters. 7 2.1. Introduction It is well-known that the MFPT, T (X), satisfies the following elliptic PDE in a two-dimensional domain (cf. [15]) D∆ T (X) + c(X) · ∇ T (X) = −1 , (2.1) where D is the diffusion constant or diffusivity associated with the underlying Brownian motion, and c(X) is the drift velocity. We will assume throughout that the drift velocity c is a conservative vector field, and so can be written as the gradient of some scalar potential. We then make the system dimensionless by introducing the new variables u, x, and ψ(x) as u = D T, L2 L . The x = L−1 X and ∇ψ = c(X) D primed derivatives are defined with respect to X while the un-primed ones are defined with respect to the non-dimensional variable x. Here 2L is the characteristic diameter of the domain under consideration. In particular, for a circular landscape, L is the radius of the lanndscape. The non-dimensional PDE problem for the MFPT then takes the form N ∆u(x) + ∇ψ(x) · ∇u(x) = −1 , Ωεj , x ∈ Ω\ j=1 ∂u = 0, ∂n u = 0, (2.2) x ∈ ∂Ω , x ∈ ∂Ωεj . The non-dimensional MFPT, u(x) satisfies the PDE (2.2) in the landscape Ω, which is perforated by N prey patches of small area denoted by Ωε1 , . . . , ΩεN . By definition, the MFPT is zero on the boundary of each patch, and it takes reflecting boundary condition on the domain boundary. Here we also assume that the domain Ω has size of order one while the nondimensional radius of each patch is O(ε), where ε assume that Ωεj 1. As ε → 0, we → xj , so that each patch shrinks to a point as ε → 0. We also assume that the distance between any two patches is O(1) as ε → 0. A schematic plot of the domain is shown in Fig. 2.1. 8 2.2. Asymptotic solution for the MFPT without drift Figure 2.1: patches. 2.2 A schematic plot of the landscape Ω with five small prey Asymptotic solution for the MFPT without drift We first analyze (2.2) without drift, i.e. with ∇ψ(x) = 0, so that the PDE for the MFPT becomes N ∆u(x) = −1 , Ωεj , x ∈Ω\ j=1 ∂u = 0, ∂n u = 0, (2.3) x ∈ ∂Ω , x ∈ ∂Ωεj . This problem is solved asymptotically by the method of matched asymptotic expansions (cf. [10]). Related problems involving elliptic PDE’s in perforated 2-D domains have been studied in [4, 11, 13, 16, 18–21]. Following [19] and [20], the outer expansion away from the prey patches is taken to have the form U ∼ U0 (x, ν) + σ(ε)U1 (x, ν) + · · · . Here ν is defined by ν = (ν1 , . . . , νN ), in terms of the logarithmic gauge 9 2.2. Asymptotic solution for the MFPT without drift functions νj ≡ −1/ log(εdj ) for j = 1, . . . , N . As shown below, the constant dj is obtained from a certain canonical inner problem defined near the j-th patch. This constant is called the logarithmic capacitance (cf. [14]), and it depends only on the shape of the patch. In the outer expansion above, the correction term σ(ε) is assumed to satisfy σ(ε) νjk for each j = 1, . . . , N , and for any positive power k, so that the correction term induced by U1 is beyond-all-orders with respect to all of the logarithmic terms captured by U0 . Upon substituting the outer expansion into (2.3), we obtain that U0 satisfies ∆U0 = −1 , ∂U0 = 0, ∂n x ∈ Ω \ {x1 , . . . , xN } , (2.4) x ∈ ∂Ω . In the outer region, the patches shrink to the points x1 , . . . , xN . This outer problem for U0 must be supplemented by appropriate singularity conditions as x → xj , for each j = 1, . . . , N . These singularity conditions are derived below upon matching the outer expansion to certain inner or local expansions that are to be constructed near each of the N patches. For the inner problem near the j-th patch, we define an inner variable y = ε−1 (x − xj ), and we define the corresponding magnified patch Ωj in terms of y by Ωj = ε−1 Ωεj (see Fig. 2.2). Figure 2.2: In the inner region the outer variable x is replaced by an inner, or local, variable y. 10 2.2. Asymptotic solution for the MFPT without drift In the inner region near the j-th patch, we introduce the inner solution qj by qj (y) = u(xj + εy), and we pose the following inner expansion: qj ∼ νj γj (ν)q0j (y) + α(ε, ν)q1j (y) + · · · . (2.5) Here γj is an unknown constant to be determined. The gauge function α is assumed to be beyond-all-orders with respect to the logarithmic terms, and so satisfies α νjk for any positive integer k as ε → 0. We impose that q0j grows logarithmically at infinity, and from the original PDE (2.3) for the MFPT we obtain that q0j satisfies ∆y q0j = 0 , q0j = 0 , y∈ / Ωj , where Ωj = Ωεj /ε , y ∈ ∂Ωj , q0j ∼ log |y| − log dj + O(|y|−1 ) , (2.6) as |y| → ∞ . We remark that the behavior qoj ∼ log |y| as |y| → ∞ is sufficient to determine the solution for q0j uniquely. In terms of this solution, the O(1) term in the far-field behavior, is uniquely determined. The O(|y|−1 ) unspecified term is the dipole term in the far-field behavior. The constant dj is known as the logarithmic capacitance (cf. [14]) and it depends on shape of Ωj but not on its orientation. Numerical values for dj for different shapes of Ωj are given in [14], and some of these are reproduced in Table 2.1. A boundary integral method to numerically compute dj for arbitrarily-shaped domains Ω1 is described and implemented in [5]. Upon substituting the far-field behavior of q0j as |y| → ∞ into (2.5), and re-writing the result in terms of the outer variable, we obtain from the matching condition that the outer solution U0 , which satisfies the PDE (2.4), must have the following singularity structure as x → xj for each j = 1, . . . , N : U0 (x, ν) ∼ νj γj (ν) log |x − xj | + γj (ν) as x → xj . (2.7) For each j = 1, . . . , N , (2.7) specifies both the regular and singular part of a outer solution. As such, for each j = 1, . . . , N , we have one constraint for the determination of the γj , j = 1, . . . , N . Overall, in this way, these 11 2.2. Asymptotic solution for the MFPT without drift Shape of Ωj ≡ ε−1 Ωεj Logarithmic capacitance dj circle, radius a ellipse, semi-axes a, b dj = a dj = a+b 2 √ 3 3 Γ( 31 ) h dj = ≈ 0.422h 8π 2 2 33/4 Γ( 14 ) h dj = 27/2 π3/2 ≈ 0.476h 2 Γ( 1 ) h dj = 4π43/2 ≈ 0.5902h equilateral triangle, side h isosceles right triangle, short side h square, side h Table 2.1: The logarithmic capacitance, or shape-dependent parameter, dj , for some cross-sectional shapes of Ωj = ε−1 Ωεj . constraints will lead to a linear algebraic system for the unknown γj for j = 1, . . . , N . By incorporating the strength of the logarithmic singularities, as given in (2.7), we can write the problem for U0 in the unpunctured domain Ω in terms of singular Dirac functions as N ∆U0 = −1 + 2π νk γk (ν)δ(x − xk ) , ∂U0 = 0, ∂n x ∈ Ω, (2.8) k=1 x ∈ ∂Ω . Upon applying the divergence theorem to this problem, we obtain the constraint that N νk γk (ν) = |Ω| , 2π (2.9) k=1 where |Ω| denotes the area of the landscape Ω. This is the first constraint on the unknown γj for j = 1, . . . , N . Since the right hand side of (2.8) is expressed in terms of the sum of Dirac forces, the solution for U0 can be conveniently written in terms of the 12 2.2. Asymptotic solution for the MFPT without drift Neumann Green’s functions as the superposition N U0 = 2π νk γk (ν)G(x; xk ) + χ . (2.10) k=1 Here χ is an arbitrary constant to be determined. In (2.10), the Neumann Green’s function G(x; ξ) is the unique solution of ∆G(x; ξ) = − 1 + δ(x − ξ) , |Ω| ∂G = 0, ∂n x ∈ Ω, x ∈ ∂Ω , (2.11) G(x) dx = 0 , Ω G∼ 1 (log |x − ξ| + R(ξ)) + o(1) 2π as x → ξ. Here R(ξ) is called the regular part of the Green’s function. The Neumann Green’s function can be calculated analytically in closed form for some simple landscapes Ω such as circles and rectangles (cf. [11], [13]), whereas for more general domains it must be computed numerically. Since G has a zero spatial average, we readily obtain that the constant χ in (2.10) can be interpreted as χ= 1 |Ω| U0 (x) dx , (2.12) Ω Although this does not help us in determining the value of χ, this simple observation is useful to interpret χ as the spatial average of U0 . Under the assumption that the starting point for the Brownian motion is uniformly distributed in the punctured domain, then up to negligible O(ε2 ) terms representing the areas of the prey patches, the constant χ can be interpreted as an asymptotic approximation to the average MFPT. Finally, we determine the linear algebraic system for the γj for j = 1, . . . , N . To do so, we must expand the solution (2.10) as x → xj and equate this near-field behavior with the prescribed singularity behavior in 13 2.2. Asymptotic solution for the MFPT without drift (2.7). For each j = 1, . . . , N , this yields that N νk γk (ν)G(xk ; xj ) + νj γj (ν)(log |x − xj | + Rj (xj )) + χ 2π k=1,k=j ∼ νj γj (ν) log |x − xj | + γj (ν) . In this expression, the logarithmic terms in |x − xj | agree identically (as they should), and from the non-singular terms we obtain a linear algebraic system for the γj for j = 1, . . . , N . We summarize our result in the following statement. Principal result 2.1: For ε 1, the asymptotic solution for the MFPT PDE (2.3) in the outer region is given by N U ∼ 2π νj γj G(x; xj ) + χ , (2.13) j=1 where the γj for j = 1, . . . , N and the constant χ are the solution to the N + 1 dimensional linear algebraic system consisting of the N equations N γj (νj Rj (xj ) − 1) + 2π νk γk G(xj ; xk ) = −χ , j = 1, . . . , N , (2.14) k=1,k=j coupled to the constraint N νj γj = |Ω| . 2π (2.15) j=1 Here νj = −1/ log(εdj ), where dj is the logarithmic capacitance of the j-th patch, as defined by the solution of (2.6), while G is the Neumann Green’s function with regular part R satisfying (2.11). We remark that this linear system is asymptotically diagonally dominant when νj is sufficiently small, and so is solvable for some range of νj small enough. 14 2.2. Asymptotic solution for the MFPT without drift Lemma: The linear system (2.14) has a solution for sufficiently small values of νj . Proof:Let R1 (x1 ) G(x2 ; x1 ) G(xN ; x1 ) G(xN ; x2 ) .. . Rj (xj ) ··· R2 (x2 ) ··· G(x2 ; x1 ) G= . . .. .. .. . G(x1 ; xN ) G(x2 ; xN ) · · · ν1 0 η= .. . 0 ··· 0 ν2 · · · .. . . . . 0 .. . ··· νN 0 0 γ = (γ1 , γ2 , · · · , γN )T , e = (1, 1, · · · , 1)T ; Then (2.14) and (2.9) can be re-written using matrices as γ − 2πGηγ = eχ eT ηγ = |Ω| 2π (2.16) (2.17) Multiplying eT η on both sides of (2.16) we get χeT ηe = eT ηγ − 2πeT ηGηγ Now defining ν ≡ eT ηe N = ν1 +ν2 +...+νN N (2.18) and using (2.17) in the above equation we get χ= |Ω| 2π T − e ηGηγ 2πν ν (2.19) Plugin (2.19) in (2.16) to get the following system for unknown γ (I − 2πGη + 2π T |Ω|e ee Gη)γ = ν 2πν (2.20) where I is N XN identity matrix. The coefficient matrix in 2.20 is diag15 2.2. Asymptotic solution for the MFPT without drift onally dominant for sufficiently small νj . Hence, for such small νj ’s a unique solution to (2.14) is guaranteed. χ can be evaluated from the γ values so obtained using (2.19). This system incorporates all of the logarithmic gauge functions in the asymptotic solution for the MFPT, leaving an error term that is beyondall-orders in −1/ log(εdj ). This error term, which we do not calculate here, arises from the local gradient behavior of G as x → xj as well as from the dipole far-field behavior of the canonical inner solution (2.6). Our asymptotic formulation of the MFPT is a hybrid asymptotic-numerical method since it uses the asymptotic analysis as a means of reducing the original problem (2.3) with N patches to the simpler asymptotically related problem (2.4) with singularity behavior (2.7). For circular prey patch shapes, for which dj is known analytically, the only numerics required to implement the hybrid formulation involves the computation of the Neumann Green’s function and its regular part at the singular point, together with the numerical solution of a linear algebraic system. For circular and rectangular shaped landscapes this Green’s function is available analytically, and hence only the solution of a linear algebraic system is required to determine the MFPT for these simple domains. An advantage of the hybrid method over the traditional method of matched asymptotic expansions is that the hybrid formulation is able to effectively “sum” the infinite logarithmic series and thereby provide an accurate approximate solution. We now compare numerical realizations obtained from this hybrid approach with full numerical results obtained by solving the PDE (2.3) for the MFPT directly. 16 2.3. Numerical verification for the MFPT without drift 2.3 Numerical verification for the MFPT without drift In this section we confirm our analysis by comparing results obtained from the analysis with corresponding full numerical results. In our comparisons, we first use the non-dimensionalised PDE to get the non-dimensional mean first passage time. The results are then rescaled by a factor of (L2 /D), to get the dimensional MFPT (roughly) according to the parameter values as given in [12]. In [12], the diffusivity of the predator in the domain was taken to be D = 0.41km2 /hr, and this value is used throughout this subsection. Figure 2.3: Model domain for numerical experiments. Prey patch is located (0,0.8) in the circular-shaped landscape. Although our asymptotic analysis is valid for any arbitrary-shaped landscape and arbitrary patch shapes, for simplicity, we will only compare asymptotic and numerical results for a circular-shaped landscape domain with circular-shaped prey patches. In the numerical comparisons below, the circular-shaped landscape has a radius, L = 1km. In addition, the prey patches are circles of radius 0.0067km. This leads to a non-dimensional prey-patch radius of ε = 0.0067. More details on the procedures used in the numerical experiments are available in the Appendix. The results quoted below from the hybrid and full numerical approaches are reported in terms 17 2.3. Numerical verification for the MFPT without drift of the MFPT averaged over all possible starting points on the domain (i.e. the spatial average of MFPT), defined by Tavg = 1 |Ω| T (X) dX . (2.21) Ω This quantity is, indirectly, the distributional average of the MFPT, assuming that the probability density of the predator follows a uniform distribution in the domain. As such, the constant χ in Principal result 2.1 must be calculated to obtain the asymptotic result predicted by the hybrid formulation. Although, explicit Neumann Green’s function is available for circular landscape, Greens function for the case with drift (2.31), studied in the next section, is not available. So, for consistency, we use a regularization scheme discussed in ➜2.7 to obtain (with and without drift) Green’s functions and then the asymptotic results given in Principal result 2.1. 2.3.1 The case of one prey patch Figure 2.4: (0,0.5) Full numerical result of MFPT when prey patch is located at 18 2.3. Numerical verification for the MFPT without drift We first study the effect of the patch location on the MFPT of the predator. To this end, we first consider the case where the landscape contains only one patch. Fig. 2.4 show the full numerical result of MFPT problem when prey patch is located at (0,0.5). In Fig. 2.5 we plot the average MFPT as the patch moves away from the centre of the circular disk. Figure 2.5: One prey patch: Spatial averages of MFPT obtained from Principal result 2.1 (solid curves) are compared with corresponding full numerical results computed from the PDE (2.3). The horizontal axis is the distance of the prey patch from the centre of a disk of radius 1km. From Fig. 2.5 we observe that the full numerical results (discrete points) and the results from the hybrid formulation (solid curve) are in very close agreement. The error between the numerical and asymptotic solution results range between 0.01% to 0.77%, which is O(ε) = 0.66%. The value of the average MFPT increases with the distance of the patch from the centre of the disk. The reason for this increase can be found by looking at the distance of the patch from starting points, x. 19 2.3. Numerical verification for the MFPT without drift Figure 2.6: Spatial average of distances between the patch and all points on the domain is plotted against distance of the patch from the centre. As the patch moves away from the center of the disk, its average distance from all the points in the domain increases as shown in Fig. 2.6. Hence, the farther the patch is from the centre, the less reachable it is, and so the average MFPT must increase. 20 2.3. Numerical verification for the MFPT without drift 2.3.2 The case of several prey patches Figure 2.7: patches Full numerical result of MFPT when prey is located in three Next, we compare results from our hybrid formulation with corresponding full numerical results for the case of three circular patches inside a circular domain. We consider two realizations of this configuration, corresponding to different locations of the three patches. In the non-dimensional formulation each patch is a circle of radius ε = 0.0067 in a circular domain of radius unity. The specific locations of the prey patches and a very favorable comparison of results from the hybrid formulation and the full numerical results are shown in Table 2.2. Patch1 (0.5,0.3) (0.3,0.8) Patch2 (-0.2,0.6) (0.1,-0.6) Patch3 (-0.4,-0.7) (-0.5,-0.7) avg. MFPT (hr) (asy) 1.6878 1.9921 avg. MFPT (hr) (num) 1.688 1.9935 Table 2.2: Spatial averages of the MPFT from the hybrid formulation and from full numerical simulations for the case of a circular domain of radius 1km with three circular prey patches each of radius 0.0067km. 21 2.4. Asymptotic solution for the MFPT with drift The numerical result corresponding to the second configuration in the table is shown in Fig. 2.7. From Table 2.2 we observe that, for each of the two different spatial arrangements of the three patches in the domain, the asymptotic solution is in very close agreement with the full numerical results. The maximum error between the numerical and asymptotic solution result is roughly 0.001, which is O(ε) when ε = 0.0067. The value of the average MFPT increases when the patches become further from the center of the domain. 2.4 Asymptotic solution for the MFPT with drift Next, we consider the PDE (2.2) for the MFPT allowing the presence of drift. The drift term in the operator models the centralizing tendency of a predator to its den site. Since the analysis of this problem with drift is similar to the case without drift we shall only outline the key steps in the derivation. In particular, in the outer region we will have the same outer expansion as the one in the case without drift. However, in this case a new Green’s function problem will be a central feature in the analysis. In contrast, the leading order inner problem near each patch (which accounts for all of the logarithmic correction terms) will be precisely the same as for the case of no drift studied in the previous subsection. The outer expansion is U ∼ U0 (x, ν) + σ(ε)U1 (x, ν) + . . . , so that U0 satisfies ∆U0 + ∇ψ · ∇U0 = −1 , ∂U0 = 0, ∂n x ∈ Ω \ {x1 , . . . , xN } , (2.22) x ∈ ∂Ω . The appropriate singularity structures for U0 as x → xj for each j = 1, . . . , N will be found by matching to the inner solution. In the inner region near the j-th patch, we introduce the inner variables 22 2.4. Asymptotic solution for the MFPT with drift y and q(y) by y = ε−1 (x − xj ) and q(y) = u(xj + εy). As similar to the case of no drift, we expand the inner solution q(y) in the form qj ∼ νj γj (ν)q0j (y) + α(ε, ν)q1j (y) + · · · , where γj is an unknown constant to be determined and the gauge function α is beyond-all-orders with respect to the logarithmic terms. From (2.2), we then obtain that q0j satisfies the same inner problem as for the case of no drift ∆y q0j = 0 , q0j = 0 , y∈ / Ωj , where Ωj = Ωεj /ε , y ∈ ∂Ωj , (2.23) −1 q0j = log |y| − log dj + O(|y| ) as |y| → ∞ . We remark that the drift term would influence the beyond-all-orders term q1j . Upon matching the inner and outer solutions, we obtain that U0 satisfies the PDE (2.22) subject to the singularity conditions U0 (x, ν) ∼ νj γj (ν) log |x − xj | + γj (ν) as x → xj , j = 1, . . . , N . (2.24) By examining the strength of the required singularity, we can express the outer problem (2.22) for U0 in the entire domain Ω in terms of the singular Dirac delta function forces as N ∆U0 + ∇ψ · ∇U0 + 1 = 2π νk γk (ν)δ(x − xk ) , ∂U = 0, ∂n x ∈ Ω, (2.25) k=1 x ∈ ∂Ω . Upon defining P (x) by P = eψ , the PDE (2.25) can be re-written in 23 2.4. Asymptotic solution for the MFPT with drift divergence form as N ∇ · (P ∇U0 ) = −P + 2πP νk γk (ν)δ(x − xk ) , x ∈ Ω, (2.26) k=1 ∂U = 0, ∂n x ∈ ∂Ω . From the divergence theorem we obtain the constraint that N 2π P (xk )νk γk (ν) = P (x) dx . (2.27) Ω k=1 Next, we express U0 as a superposition of a smooth function U0H (x) and a certain Green’s function in the form N U0 = U0H + 2π νk γk (ν)G(x; xk ) + χ . (2.28) k=1 Here χ is an arbitrary constant and the smooth part UOH satisfies ∇ · (P ∇UOH ) = −P + Pave , ∂UOH = 0, ∂n x ∈ Ω, x ∈ ∂Ω , (2.29) UOH dx = 0 . Ω Here Pave is the average of P over Ω, defined explicitly by Pave ≡ 1 |Ω| P (x) dx . (2.30) Ω Notice that the zero average condition on U0H is the condition that ensures that U0H is uniquely determined. The singular part in this decomposition consists of the Green’s function 24 2.4. Asymptotic solution for the MFPT with drift G(x; ξ), which satisfies ∇ · (P ∇G(x; ξ)) = − P (ξ) + P (x)δ(x − ξ) , |Ω| ∂G = 0, ∂n x ∈ Ω, x ∈ ∂Ω , (2.31) G(x; ξ) dx = 0 , Ω G∼ 1 (log |x − ξ| + R(ξ)) + o(1) 2π as x → ξ . Solving for this Green’s function is challenging. To numerically solve the PDE for the Green’s function, we need a combination of regularization and numerical techniques. This approach is discussed later in the chapter in ➜2.7. The constant χ in (2.28) can again be interpreted as χ= 1 |Ω| U0 (x)dx (2.32) Ω Finally, to determine the linear algebraic system for γj , we expand the solution in (2.28) as x → xj and equate the resulting expression with the required singular behavior in (2.24). This leads to N νk γk (ν)G(xk ; xj ) + νj γj (ν)(log |x − xj | + Rj (xj )) + χ U0H (xj ) + 2π k=1,k=j = νj γj (ν) log |x − xj | + γj (ν) . From the non-singular terms in the expression above we obtain a linear algebraic system for the γj . We summarize our result as follows: Principal result 2.2: For ε 1, the asymptotic solution for the MFPT PDE (2.2) in the outer region is given by N U ∼ U0H + 2π νk γk (ν)G(x; xk ) + χ , (2.33) k=1 where the γj for j = 1, . . . , N and the constant χ are the solution to the 25 2.5. Numerical verification for the MFPT with drift N + 1 dimensional linear algebraic system consisting of the N equations N γj (νj Rj (xj ) − 1)+2π νk γk (ν)G(xj ; xk ) = −U0H (xj )−χ , j = 1, . . . , N k=1,k=j (2.34) coupled to the constraint N 2π P (x) dx . P (xk )νk γk (ν) = (2.35) Ω k=1 Here P = eψ , νj = −1/ log(εdj ), where dj is the logarithmic capacitance of the j-th patch, U0H is the smooth solution satisfying (2.29), while the Green’s function G with regular part R under the assumption of drift is the solution of (2.31). 2.5 Numerical verification for the MFPT with drift In this section we numerically verify our analysis for the case of drift by using Matlab. We first compute the non-dimensional MFPT and then rescale it to get the dimensional version. Together with the parameter values used for the case without drift, we now introduce the drift velocity. We take L ψ = −0.085 D |x|, hence the drift velocity c(X) has a magnitude 0.085 km/hr ( given in [12]) and is directed towards the centre of the circular domain, giving the predator a centralizing tendency. L and D values are same as the ones used in the case without drift. The domain and prey patch sizes are the same as those in previous section. However, instead of taking spatial averages for the MFPT, in this subsection we take the distribution average of MFPT over all possible starting points on the domain defined by Tavg ≡ S(X)T (X) dX . (2.36) Ω 26 2.5. Numerical verification for the MFPT with drift Here S(X) is the probability distribution of the location of predator in the domain at steady state. The precise expression for S(X) is derived towards the end of the chapter in ➜2.8. For the circular landscape with the above parameters we obtain U0H = 1 βr + Pave eβr (βr − 1) + E1 (−βr) + log(−βr) +33.6091 (2.37) β2 where β = 0.085 0.41 , r = |x| and E1 (z) = ∞ e−t t z dt is the exponential integral. We use this U0H for evaluating the outer solution. 2.5.1 The cases of one and several prey patches Figure 2.8: Full numerical result of MFPT for when prey patch is located at (0,0.5) Similar to the case without drift, we first study the effect of the location of the patch in the domain on the MFPT of the predator. Fig. 2.8 shows the full numerical result of MFPT when prey patch is located at (0,0.5). The results for various patch locations are shown in Fig. 2.9. 27 2.5. Numerical verification for the MFPT with drift Figure 2.9: Distribution averages of the asymptotic and full numerical results for the MFPT are plotted versus the distance of the patch from the centre of the circular disk From this figure we observe that the numerical and asymptotic results are in very close agreement. The maximum relative error between the numerical and asymptotic results is 3.3%, which is of O(ε) = 0.67%. Similar to the case without drift, here MFPT increases with the distance of the patch from the centre. 28 2.5. Numerical verification for the MFPT with drift Figure 2.10: Full numerical result of MFPT for when prey is located in three patches Next, we compare the numerical and asymptotic solution averages when the domain contains three circular patches each of radius ε. Two configurations of three patches are considered. The results are shown in Table 2.3. Patch1 (0.5,0.3) (0.3,0.8) Patch2 (-0.2,0.6) (0.1,-0.6) Patch3 (-0.4,-0.7) (-0.5,-0.7) avg. MFPT (hr) (asy) 1.5990 1.8993 avg. MFPT (hr) (num) 1.7249 2.1121 Table 2.3: Spatial averages of MPFT in a domain with three patches, obtained using both the hybrid asymptotic formulation in Principal result 2.2 and the full numerical results computed from (2.2). The numerical result corresponding to the second configuration in the table is shown in Fig. 2.10. As we can see from the results above, for two different arrangements of the three patches in the domain, the asymptotic solution is in close agreement with the numerical solution. The maximum relative error between the numerical and asymptotic distribution averages is 10%, which is a little higher than the case of no drift. But the maximum 29 2.6. Comparison of the cases with and without drift error between spatial average of same MFPT is just 4%. The higher error in distribution averages could be due to the fact that the distribution function, S has a singularity near the centre of the domain, and we use S in the calculation of the distribution averages of MFPT. The estmation of S near the singularity is prone to small errors. 2.6 Comparison of the cases with and without drift In this subsection, we compare the spatial average of MFPT without drift and the distribution average of MFPT for the case with drift. We consider a single circular prey patch, where the distance of the prey patch to the boundary of the circular landscape is allowed to vary. As shown in Fig. 2.11, for the case where c(X) has magnitude 0.085 km/hr and is directed to the origin there is very little difference between these two curves and they essentially overlap. 30 2.6. Comparison of the cases with and without drift Figure 2.11: Averages of MFPT are plotted against distance of the patch from the centre. Here the drift velocity has magnitude 0.085 km/hr. However, when the magnitude of the drift velocity is increased by a factor of 10 to 0.85km/hr, there is an interesting cross-over effect between the two curves. This is shown in Fig. 2.12. 31 2.6. Comparison of the cases with and without drift Figure 2.12: In the plot, averages of MFPT are plotted against distance of the patch from the centre. The magnitude of the drift velocity is 10 times larger than in Fig. 2.11. From this figure we observe that with larger drift and a centralizing tendency, the MFPT with drift is significantly smaller initially than the MFPT with no drift. However, as the patch moves away from the centre, it eventually increases and overtakes the no-drift MFPT result. This suggests that denning animals, which have a centralizing tendency, have a larger search time when the prey is located farther from the den. It also suggests that there exists a critical drift velocity at which random diffusion is favoured over diffusion with drift. 32 2.7. Numerical approximation (regularization scheme) of Green’s function 2.7 Numerical approximation (regularization scheme) of Green’s function The Green’s function in the presence of drift is the unique solution of ∇ · (P ∇G(x; ξ)) = − P (ξ) + P (x)δ(x − ξ) , |Ω| ∂G = 0, ∂n x ∈ Ω, x ∈ ∂Ω , (2.38) G(x; ξ) dx = 0 , Ω G∼ 1 (log |x − ξ| + R(ξ)) + o(1) 2π as x → ξ . For the special case where P ≡ 1 or ψ ≡ 0, this Green’s function reduces to the Neumann Green’s function, which is relevant to the analysis of the MFPT with no drift. Although the exact analytical solution of the ’no-drift’ Neumann Green’s function for circular and rectangular domains is available ([11], [13]), there are no similar results for the Green’s function with drift. The difficulty arises from the domain shape, location of the singularity and the arbitrary nature of the function, ψ. Therefore, we must in general compute the Green’s function for the case with drift numerically. However, determining the numerical solution for G also imposes a considerable challenge. This difficulty in solving for G numerically arises from one of the constraints which states that the spatial average of the Green’s function is zero. It is a cumbersome task to add this rank-one constraint, while solving the PDE using conventional numerical methods. Without this constraint, the solution for G is only known up to an arbitrary constant. Therefore, as a way to tackle this problem we employ a regularization method motivated by the analysis in [20]. The approach involves solving a numerically solvable ’Master’ PDE using a combination of asymptotic and numerical techniques. The required Green’s function will arise as a by-product in the process of finding the solution to the ’Master’ PDE. 33 2.7. Numerical approximation (regularization scheme) of Green’s function We now discuss this regularization approach for computing the required Green’s function. Consider an auxiliary function φ(x, ξ) that satisfies the following PDE: ∇(P ∇φ(x; ξ)) − µφ = P δ(x − ξ) , x ∈ Ω, ∂φ = 0, x ∈ ∂Ω , ∂n 1 φ∼ (log |x − ξ| + R0 (ξ)) + o(1) 2π as x → ξ . (2.39) Here µ is a small positive parameter. Lemma: For µ > 0 the PDE (2.39) has a unique solution. Proof: Let φa and φb be two different solutions of the PDE then, ∇ · (P ∇φa (x; ξ)) − µφ = P δ(x − ξ) , ∂φa = 0, ∂n x ∈ Ω, x ∈ ∂Ω , and ∇ · (P ∇φb (x; ξ)) − µφ = P δ(x − ξ) , ∂φb = 0, ∂n x ∈ Ω, x ∈ ∂Ω . Subtracting one PDE from the other and defining φc ≡ (φa − φb ), we obtain ∇ · (P ∇φc ) − µφc = 0 , x ∈ Ω, ∂φc = 0, ∂n x ∈ ∂Ω . (2.40) ⇒ P ∆φc + ∇P.∇φc − µφc = 0 ⇒ φc P ∆φc + φc ∇P.∇φc − µφ2c = 0 34 2.7. Numerical approximation (regularization scheme) of Green’s function ⇒ φ2c dx − φc P ∇φc .ˆ nds = µ ∇(φc P ).∇φc dx + φc ∇P.∇φc dx Ω Ω ∂Ω Ω φc ∇P.∇φc dx Ω Ω Ω − φ2c dx − φc P ∆φc dx = µ (∵ From Green’s second identity) Ω Ω (∵ φc .ˆ n= φ2c dx + µ Ω ∇(φc P ).∇φc dx − Ω φ2c dx − ∇(φc P ).∇φc dx = µ − φc ∇P.∇φc dx Ω ∂φc = 0) ∂n φc ∇P.∇φc dx = 0 Ω φ2c dx + P µ Ω |∇φc |2 dx = 0 Ω We know that P > 0, and if µ > 0, then φc must be identically zero. If φc = 0 then φa = φb which means that the PDE (2.39) has a unique solution. As µ → 0, the solution to the PDE (2.39) can be expressed as the following asymptotic series: φ∼ φ0 + φ1 + µφ2 + µ2 φ3 + · · · . µ (2.41) Upon substituting this expansion into the PDE (2.39) and equating the coefficients of µ at each order, we obtain a set of PDEs for each order. For the O(µ−1 ) term, we obtain that φ0 satisfies ∇(P ∇φ0 ) = 0 , ∂φ0 = 0, ∂n x ∈ Ω, (2.42) x ∈ ∂Ω , with φ0 smooth in the domain. The solution to this problem is that φ0 is an arbitrary constant. 35 2.7. Numerical approximation (regularization scheme) of Green’s function Similarly, from the O(1) term we obtain the PDE: ∇ · (P ∇φ1 ) = φ0 + P δ(x − ξ) , x ∈ Ω, ∂φ1 = 0, x ∈ ∂Ω , ∂n 1 φ1 ∼ log |x − ξ| + O(1) 2π (2.43) as x → ξ. To solve for φ1 we must have, from the divergence theorem, that Ω φ0 dx + P (ξ) = 0 ⇒ φ0 = −P (ξ) |Ω| . Similarly, from the O(µ) terms, we obtain the following PDE for φ2 : ∇(P ∇φ2 ) = φ1 , x ∈ Ω, (2.44) ∂φ2 = 0, ∂n x ∈ ∂Ω , where φ2 is pointwise bounded in Ω. From the divergence theorem, we obtain that the solvability condition for the PDE for φ2 is Ω φ1 dx =0 Now using this constraint we rewrite the PDE for φ1 as ∇ · (P ∇φ1 ) = − P (ξ) + P (x)δ(x − ξ) , |Ω| ∂φ1 = 0, ∂n x ∈ Ω, x ∈ ∂Ω , (2.45) φ1 dx = 0 , Ω φ1 ∼ 1 (log |x − ξ| + A1 ) + o(1) 2π as x → ξ. The PDE for φ1 is same as the PDE for the Green’s function to be found. Hence the O(1) term in the expansion (2.41) of the solution to the regularized PDE (2.39) is the required Green’s function. The solution φ, to (2.39), can be calculated numerically using conventional techniques such as the finite element method. Once we have φ for two small values of µ, we 36 2.8. Distribution of the predator in the domain can use Richardson extrapolation to approximate φ1 and hence the Green’s function with an accuracy of O(µ2 ). 2.8 Distribution of the predator in the domain For animals with a centralizing tendency the initial probability distribution, S(X) in a 2-D domain is given by the following Kolmogorov equation, ∂S = −∇ .(cS + ∆ (DS) ∂t ∂S =0 on the domain boundary ∂n (2.46) where D is the diffusion constant or diffusivity, and c(X) is the drift velocity. The primed operators are defined with respect to X. In the context of the problem, we find the solution of the PDE on a circular domain of radius L. The solution we find for a circular domain will be used as approximate initial distribution of predator in a circular landscape with prey patches. The steady state PDE corresponding to (2.46) will be −∇ .(cS) + ∆ (DS) = 0 (2.47) The PDE is solved in dimensional circular coordinates with X = (r, θ) such that r=|X|. If we take c=0, this PDE gives uniform distribution in the domain and S will be a constant. Let c= −c0 ˆr where c0 is a non-negative real number and ˆr is the outward unit vector in radial direction. Now we plug c= −c0 ˆr in (2.47) solve for S as follows 37 2.8. Distribution of the predator in the domain ∆ (DS) = −c0 (S∇ .ˆr + ˆr.∇ S) ∂ D ∂ ∂S ˆ ∂S ∆ S = S(ˆr + θˆ ).ˆr + ˆr.(ˆr +θ ) c0 ∂r r∂θ ∂r r∂θ ∂ˆr ˆ ∂ˆr ∂S = S(ˆr. + θ. )+ ∂r r∂θ ∂r 2 2 ∂S S ∂S D ∂ S ∂ S )= + =⇒ − ( 2 + + c0 ∂r r∂r r2 ∂θ2 r ∂r (2.48) ∂ˆ r ∂ˆ r ˆ θ.ˆ ˆ r = 0. Now if we look Here we used the identities: = 0, = θ, − ∂r ∂θ for a radially symmetric solution, the simplified PDE will be Srr + 2π c0 S Sr + (Sr + ) = 0 r D r S=0 L r=L (2.49) L S(r) rdrdθ = 2π 0 0≤r≤L 0 S(r)rdr = 1 0 The third condition comes from the constraint on probability distribution. Now solve the PDE to find S (rSr ) c0 (rS) =− r D r c0 =⇒ rSr = − rS + k1 D k1 c0 =⇒ Sr + S = D r (k1 is a constant) (2.50) Now if we use the integrating factor, c0 k1 c0 S) = e D r D r c r D0 r c0 e r =⇒ e D S(r) = k1 dr + b r L c0 =⇒ e D r (Sr + (2.51) (b is a constant) c0 Since at r=L, S = 0, we get b = 0. Hence S(r) = k1 e− D r c0 r eDr L r dr where k1 can be found from the third constraint of 2.49. The non-dimensional distribution, S corresponding to the non-dimensional MFPT can be obtained 38 2.8. Distribution of the predator in the domain by taking r∗ = |x| = r/L and β = c0 L/D and s(r∗ ) = S(r) such that ∗ r∗ −βr∗ s(r ) = k1 e 1 ∗ eβr dr∗ r (2.52) s(r∗ ) in a 2D circular domain is plotted in Fig. 2.13. Figure 2.13: s(r∗ ) is plotted on 2D domain by taking c0 = 0.085 km/hr. 39 Chapter 3 Asymptotic Calculation of the Splitting Probability 3.1 Introduction In this chapter we use the hybrid asymptotic-numerical method for summing logarithmic expansions to investigate a problem relating to the calculation of a splitting probability. Here we consider a prey animal diffusing randomly in a landspace with prey and predator patches. The landscape contains one prey patch (also referred to as the target patch) and more than one predator patch. In this context, the splitting probability, u(x) is defined as the probability that the prey reaches the target patch before hitting the predator patches. The initial location of the prey is at an arbitrary location, x inside the domain, representing the 2-D landscape. Essentially, the splitting probability is the probability of the prey reaching home safely starting from a given point inside the domain. In this chapter we analyse the system with no drift velocity and a constant diffusivity. A schematic plot of the domain is shown in Fig. 3.1. 40 3.1. Introduction Figure 3.1: In the domain Ω shown in the figure, we have four predator patches and one prey/target patch As a remark, a mathematically equivalent problem leading to the calculation of a splitting probability is to consider prey undergoing free diffusion in a domain with a collection of N localized patches, including N-1 predator patches and one prey patch. The goal is to calculate the probability that the prey, starting from an arbitrary position x in the domain, reaches the targeted prey site first, before encountering any of the other N − 1 predator sites. In this context, as we shall show below, the effect of the “sea” of other predator sites that can act as a shield, which ultimately lowers the probability of reaching the target. The probability, u(x) of reaching the target patch before the predator patches satisfies the following elliptic PDE (cf. [15]): N ∆u(x) = 0 , Ωεj , x∈Ω\ j=1 ∂u = 0, ∂n u = δ1j , (3.1) x ∈ ∂Ω x ∈ ∂Ωεj j = 1, . . . , N , where δ1j is the Kronecker delta function, with δ1j = 1 if j = 1 and δ1j = 0 otherwise. 41 3.2. Asymptotic analysis The splitting probability satisfies Laplace equation on the perforated domain. For convenience we denote the prey/target patch by Ωε1 and all the other N − 1 predator patches by Ωε2 , . . . , ΩεN . By definition, the splitting probability takes a value of unity on the prey/target patch and zero on predator patches. Here we also assume that the domain Ω has size of order one while the patches have a size of O(ε) where ε every patch Ωεj 1. Also, as ε → 0, → xj , i.e. every patch shrinks to a point strictly inside the domain. We also assume that the distance between the patches is of O(1) as ε → 0. 3.2 Asymptotic analysis We begin analysis of the PDE (3.1) by formulating an appropriate outer expansion. Following [20] and [19], the outer expansion, away from the patches, has the form U ∼ U0 (x, ν) + σ(ε)U1 (x, ν) + . . . . (3.2) Here the vector ν is defined by ν = (ν1 , . . . , νN ), where νj is given in terms of the logarithmic capacitance dj of the patch by νj = −1 log εdj . It is defined in terms of the far-field behavior of the inner problem (3.5). As in the previous chapter, the gauge function σ is beyond-all-orders with respect to the logarithmic terms. Upon substituting the expansion (3.2) into (3.1) we obtain that the outer solution U0 satisfies the PDE ∆U0 = 0 , x ∈ Ω \ {x1 , . . . , xN } , ∂U0 = 0, ∂n x ∈ ∂Ω . (3.3) Since, in the outer region, the patches shrink to the points x1 , . . . , xN as ε → 0, this problem must be supplemented by appropriate singularity conditions as x → xj for each j = 1, . . . , N . These singularity conditions are 42 3.2. Asymptotic analysis determined below only upon matching the outer expansion to inner expansions, one near each patch. For the inner problem near the j-th patch, we define an inner variable y = ε−1 (x − xj ), and we define the corresponding magnified inner domain Ωj in terms of y is given by Ωj = ε−1 Ωεj (see Fig. 3.2). Figure 3.2: In the inner region the outer variable x is replaced by an inner or local variable y. In terms of the local variable y, we define the inner solution q(y) as q(y) = u(xj + εy), and we expand it as qj ∼ δ1j + νj γj (ν)q0j (y) + α(ε, ν)q1j (y) + . . . , (3.4) where γj is an unknown constant to be determined. After substituting this expansion into (3.1), we obtain the following problem for q0j : ∆y q0j = 0 , q0j = 0 , y∈ / Ωj where Ωj = Ωεj /ε , y ∈ ∂Ωj , q0j = log |y| − log dj + O(|y|−1 ) , (3.5) as |y| → ∞ . As discussed in the previous chapter, this problem uniquely defines the logarithmic capacitance, dj . An exact analytical solution for q0j can be found for some simple patch shapes such as circular patches, elliptical patches, etc. 43 3.2. Asymptotic analysis For more complex patch shapes, q0j must be computed numerically. Next, by using the far-field behavior of q0j as |y| → ∞, we obtain the required singularity behavior for the corresponding outer solution U0 in the form U0 (x, ν) ∼ δ1j + νj γj (ν) log |x − xj | + γj (ν) as x → xj , (3.6) for each j = 1, . . . , N . The problem for the outer solution is then (3.3) subject to the singularity behaviors (3.6) for each j = 1, . . . , N . Since both the regular and singular parts of the singularity structure in (3.6) are prescribed, the problem for U0 is in general under-specified, and is solvable only if the γj for j = 1, . . . , N satisfy a certain linear algebraic system. By incorporating the strengths of the logarithmic singularities (γj ’s) at the patch centres, we can re-express the right hand side of the outer PDE (3.3) in the unperforated domain Ω as the sum of Dirac delta functions as N νk γk (ν)δ(x − xk ) , ∆U0 = 2π x ∈ Ω, (3.7) k=1 ∂U0 = 0, ∂n x ∈ ∂Ω . By applying the divergence theorem to (3.7), we obtain the constraint that N νk γk (ν) = 0 . 2π (3.8) k=1 This is the first constraint on the as yet unknown γj for j = 1, . . . , N . Since the right hand side of the reformulated PDE is a sum of delta functions, the solution for U0 can be expressed as a superposition of the Neumann Green’s function as N U0 = 2π νk γk (ν)G(x; xk ) + χ . (3.9) k=1 Here χ is an arbitrary constant to be determined and the Neumann Green’s 44 3.2. Asymptotic analysis function G(x; ξ) is the unique solution of ∆G(x; ξ) = − 1 + δ(x − ξ) , |Ω| ∂G = 0, ∂n x ∈ Ω, x ∈ ∂Ω , (3.10) G(x; ξ) dx = 0 , Ω G= 1 (log |x − ξ| + R(ξ)) + o(1) 2π as x → ξ . Here |Ω| is the area of the domain and R(ξ) is called the regular part of the Green’s function. The Neumann Green’s function is known analytically for simple domains such as circles and rectangles (cf. [11], [13]), while we need numerical approximation techniques to compute it for other domain shapes. The constant χ in (3.9) has the interpretation that χ= 1 |Ω| U0 (x)dx . (3.11) Ω Finally, we determine the linear algebraic system for the γj by accounting for the regular or non-singular parts of the singularity structures in (3.6). By expanding (3.9) as x → xj , and comparing the resulting expression with (3.6) we obtain the matching condition N νk γk (ν)G(xk ; xj ) + νj γj (ν)(log |x − xj | + Rj (xj )) + χ 2π k=1,k=j ∼ δ1j + νj γj (ν) log |x − xj | + γj (ν) . The singular terms match identically, while the non-singular terms provide the linear algebraic system for the γj for j = 1, . . . , N . In this way we obtain the following main result: Principal result 3.1: For ε → 0, the asymptotic solution to the splitting probability PDE (3.1) is given in the outer region |x − xj | O(ε) for j = 45 3.3. Numerical verification 1, . . . , N by N U0 ∼ 2π νk γk (ν)G(x; xk ) + χ , (3.12) k=1 where γj for j = 1, . . . , N and χ are to be determined from the N + 1 dimensional linear algebraic system N γj (ν)(νj Rj (xj ) − 1) + 2π νk γk (ν)G(xj ; xk ) + χ = δ1j , j = 1, . . . , N , k=1,k=j N 2π νk γk (ν) = 0 . k=1 (3.13) This linear system is solvable for νj sufficiently small. 3.3 Numerical verification Figure 3.3: Model domain for numerical experiments. Predator patches are located on a concentric ring within a circular-shaped landscape and prey patch in located outside the ring In this section we verify results of the asymptotic analysis with full numerical results computed from (3.1) using Matlab. Although the analysis is valid for any arbitrary domain and patch shapes, it is easy to verify numerically for a circular landscape and for circular patch shapes for which the loga46 3.4. Results for rings rithmic capacitance is known analytically. As in the previous chapter, the Neumann’s Greens function is calculated using the regualarization scheme detailed in ➜2.7 . In each of the experiments below we use a circular domain of radius 1km. The predator and prey patches are taken to be circles of radius 0.0067km. Therefore, we take ε = 0.0067 in the dimensionless formulation. Most of the results from analytical and numerical approaches are quoted as splitting probabilities averaged over all possible starting points on the domain (i.e. spatial averages of the probability). More details on the procedures used in numerical experiments are available in the Appendix. 3.4 Results for rings Figure 3.4: Full numerical results for the splitting probability when the target is located outside the ring of four predators For our first example, we consider the case where there are either four or eight predator patches equally spaced on a concentric ring inside the circular domain. The prey/target patch is located either outside or inside the ring. The case of eight predator patches on the ring is shown schematically in Fig. 3.3. For our first experiment we take four predator patches located at (0, 0.3), 47 3.4. Results for rings (0, −0.3), (0.3, 0), and (−0.3, 0), on a ring of radius 0.3. For two locations of the prey patches, in Table 3.1 we show a very favorable comparison between the results of the hybrid formulation in Principal result 3.1 and results from full numerical simulation of the PDE (3.1). The computational result for u(x) is shown in Fig. 3.4. Prey patch location (0,0.001) (0.7,0.7) Asymptotic 0.14716 0.12799 Numerical 0.14737 0.12965 error 0.00021 0.00166 Table 3.1: Spatial averages of the splitting probability when the target is located inside/outside a four-predator concentric ring. As shown in Table 3.2, a similar very favorable agreement between the asymptotic and full numerical results occurs for the case where the four predator patches are located on circle of radius 0.5 at the locations (0.3, 0.4), (−0.3, 0.4), (0.3, −0.4), and (−0.3, −0.4). Prey patch location (0,0.001) (0.5,0.5) Asymptotic 0.1720239 0.15896 Numerical 0.1720191 0.158992 error 0.000005 0.000032 Table 3.2: Spatial averages of the splitting probability when the target is located inside/outside the four predator ring. 1 1 , ± 2√ ), In our next experiment, we take eight predator patches located at (± 2√ 2 2 (±0.5, 0), and (0, ±0.5). The very close agreement between the asymptotic and full numerical results for this case is shown in Table 3.3. Prey patch location (0.001,0) (0.65,0.65) Asymptotic 0.08805 0.08379 Numerical 0.08804 0.08374 error 0.00001 0.00005 Table 3.3: Spatial averages of splitting probability when the target is located inside/outside an eight-predator ring. The asymptotic results from our hybrid formulation are in very close agreement with the numerical results and the error is O(ε). Intuitively, the 48 3.4. Results for rings probability should be higher when the patch is outside the ring of predators than inside. That is, it should be safer for a randomly diffusing prey, if its home is not surrounded by the predators. But, as we can see, the splitting probability is higher when the prey patch is inside the ring. This suggests that the predator patches may not be large enough to shield the prey patch. 3.4.1 Results with larger rings Figure 3.5: Full numerical results for the splitting probability when the target is located outside the ring of eight predators when the radii of each predator patch is 5ε where ε = 0.0067. To examine the shielding effect we increased the predator and home patch radii by a factor of five and performed the analysis again. We again have 1 1 eight predator patches located at (± 2√ ,± 2√ ), (±0.5,0) and (0,±0.5). The 2 2 results are shown in Table 3.4 and in Fig. 3.5. Prey patch location (0,0.001) (0.6,0.6) Asymptotic 0.07304 0.08726 Numerical 0.07423 0.08695 error 0.00119 0.00031 Table 3.4: Spatial averages of splitting probability with large patches of radii 5ε, where ε = 0.0067. 49 3.5. Effect of distance from the prey patch When we increased the patch radii by a factor of five, to 5ε where ε = 0.0067, we were able to observe the shielding phenomenon, i.e. the probability of reaching the target is higher if the prey starts inside the ring. It should be noted that, although we have used larger patches, the error is still O(ε). 3.5 Effect of distance from the prey patch In this section we calculate the splitting probability, u(x) as a function of the prey’s starting position, x. For this purpose we study a scenario in which we use a circular domain of unit radius and eight predator patches are located on a ring of radius 0.5, while the prey/home patch is located at (0,0.1). In Fig. 3.6, we plot the splitting probability as the prey’s starting location changes from (−1, 0) to (1, 0) on the horizontal axis on the domain. 50 3.5. Effect of distance from the prey patch Figure 3.6: Top figure shows the case under consideration. Bottom figure shows the splitting probability as a function of the x-coordinate of the starting point of the prey on horizontal axis in the domain. From this figure, we observe that the splitting probability is largest when the prey starts close to its home/prey patch and almost zero when it starts near the ring of predators. Also, the farther the starting point from the predator ring, the higher the probability that the prey reaches home. A similar result was shown in [3] for the three-dimensional case. 51 3.6. The effect of closely spaced patches 3.6 The effect of closely spaced patches Throughout our analysis so far, we have assumed that the patches are located at O(1) distances from each other. The analysis illustrated so far requires modification when any two patches become separated by distances of O(ε). In such a situation, the two patches can be effectively combined into a patch of size O(ε), and we can use an analysis similar to the one in the previous section with some changes in the formulation and solution of the inner problem. 3.6.1 Two closely spaced patches First we consider a scenario where two patches (a prey patch and a predator patch) are close to each other. Here Ω1 and Ω2 represent prey and predator patches respectively. The governing PDE for the case of two closely spaced patches is 2 ∆u(x) = 0 , Ωεj , x∈Ω\ j=1 ∂u = 0, ∂n u = 1, u = 0, x ∈ ∂Ω , (3.14) x ∈ ∂Ωε1 , x ∈ ∂Ωε2 . When there are two patches (one prey patch and one predator patch) that are close to each other, from the outer region, they look like one point located near the patches. For convenience, this reference point is chosen at the center of the one of the patches. Thus, the outer correction satisfies ∆U0 = 0 , x ∈ Ω \ {x1 } , ∂U0 = 0, ∂n x ∈ ∂Ω , (3.15) subject to an appropriate singularity condition as x → x1 . Here x1 is the center of the bigger patch. Similar to the analysis in the previous section, 52 3.6. The effect of closely spaced patches the inner expansion in terms of the inner variable y = ε−1 (x − x1 ) has the form q ∼ q0 (y) + α(ε, ν)q1 (y) + · · · . (3.16) The function q0 is then decomposed into two separate terms as q0 = V0 + νγ(ν)Vc . (3.17) Here γ is an unknown constant to be determined and ν = −1 log εd , where d is the effective logarithmic capacitance of the two patches, as defined below. The gauge function σ in (3.16) is assumed to be beyond-all-orders with respect to the logarithmic terms. Upon substituting (3.16) and (3.17) into (3.14), we obtain that V0 satisfies y∈ / Ω0 ∆y V 0 = 0 , where Ωj = Ωεj /ε and Ω0 = (Ω1 ∪ Ω2 ) , V0 = 0 , y ∈ ∂Ω1 , V0 = 1 , y ∈ ∂Ω2 , V0 ∼ V∞ as |y| → ∞ , (3.18) where V∞ is a constant to be determined. In contrast, Vc satisfies ∆y V c = 0 , Vc = 0 , y∈ / Ω0 , y ∈ ∂Ω0 , Vc ∼ log |y| − log d + O(|y|−1 ) (3.19) as |y| → ∞ , where the far-field behavior determines d uniquely. It should be noted that in the definition of the inner problem the patches Ωε1 and Ωε2 are magnified to Ω1 and Ω2 , respectively, and their union is used. An exact solution for q0 can be found for the case where both patches are circles. In this case we can use conformal mapping to solve for V0 and then evaluate V∞ from it. In addition, for this situation we know the far-field behaviour of Vc . This allows us to determine the far-field behavior of q0 , and then write the appropriate singularity condition for the outer solution 53 3.6. The effect of closely spaced patches U0 . Using conformal mapping to solve the inner problem with circular patches In this subsection we solve (3.18) for V0 for the case with two circular patches, and from this solution we compute the required far-field constant V∞ . A plot of the inner geometry is shown in Fig. 3.7. Figure 3.7: The inner problem with two circular patches. We now outline the geometry of the inner problem more precisely. Let a1 be the radius of patch I located at y = (0, 0), and a2 be the radius of patch II located at y = (β, 0). Here without loss of generality we assume that the centres of both circles lie on the horizontal axis in the y coordinate plane. We let β > 0 denote the distance between the centers of the patches measured in the y coordinate. We solve this problem by conformal mapping. We let (x1 , 0) and (x2 , 0) be the real symmetric points of the circles. Then by definition of symmetric 54 3.6. The effect of closely spaced patches points we have x1 x2 = a21 (x1 − β)(x2 − β) = a22 ⇒ x1 x2 − β(x1 + x2 ) + β 2 = a22 a21 − a2 (x1 + x2 ) + β 2 = a22 x1 + x2 = η where η= a21 − a22 + β 2 . β Upon solving for x1 and x2 we get x1 = η+ η 2 − 4a21 2 x2 = η− η 2 − 4a21 . 2 We now show that, provided that the patches do not touch each other, i.e. when β > (a1 + a2 ), x1 and x2 are real. Lemma: If β > (a1 + a2 ) then x1 and x2 are real. Proof: For x1 and x2 to be real, the required condition is η 2 − 4a21 > 0 =⇒ (η − 2a1 )(η + 2a1 ) > 0 for β > a1 + a2 =⇒ η > 0 Hence we need η − 2a1 > 0 a21 − a22 + β 2 − 2a1 > 0 β a21 − a22 + β 2 − 2a1 β >0 β a21 − a22 + β 2 − 2a1 β > 0 (β − a1 )2 − a22 > 0 =⇒ (β − a1 − a2 )(β − a1 + a2 ) > 0 This inequality holds when β > a1 + a2 . 55 3.6. The effect of closely spaced patches Mobius transformation Next we use the following Mobius transformation thats maps the z-plane to the w-plane: w = f (z) = z − x1 . x2 − z Under this transformation, let φ be the image of V0 , in w-space. We denote the images of circles I and II by I and II respectively. Let R1 and R2 denote the radii of I and II respectively. A schematic plot of the mapping is shown in Fig. 3.8. Figure 3.8: concentric circles in w-space We calculate the values of R1 and R2 as R1 = |w| = |f (a1 )| = R1 = 2a1 − η − η− a1 − x1 x2 − a1 η 2 − 4a21 η 2 − 4a21 − 2a1 (3.20) . Similarly we calculate R2 as R2 = |w| = |f (a2 + β)| = R2 = 2a2 + 2β − η − η− a2 + β − x1 x2 − (a2 − β) η 2 − 4a21 η 2 − 4a21 − 2a2 − 2β (3.21) . 56 3.6. The effect of closely spaced patches Solving for V0 Case 1: In this case we assume that patch I is the target/prey while patch II is the predator patch. Then V0 = 1 on I and V0 = 0 on II. We now evaluate V0 to find V∞ as z → ∞. Since I and II are concentric circles, the solution for φ can be written in terms of constants A and B as φ = A log |w|+B. Since φ = 0 on |w| = R2 and φ = 1 on |w| = R1 , we calculate φ= log |w| R2 1 log R R2 . In terms of the original z-plane we get V0 = 1 log | z−x x2 −z | − log R2 1 log R R2 = 1 /z log | 1−x x2 /z−1 | − log R2 1 log R R2 . Therefore as |z| → ∞, we can extract the far-field limiting behavior of V0 . In this way, we calculate V∞ as V∞ = log R12 , 1 log R R2 (3.22) where R1 and R2 are defined in (3.20) and (3.21). Case 2: Now we choose patch II to be the target and patch I to be predator patch. Then we have V0 = 0 on I and V0 = 1 on II. Therefore, we calculate that φ= R1 log |w| 1 log R R2 . Then upon letting |z| → ∞ we obtain after a short calculation that V∞ = log R1 . 1 log R R2 (3.23) Now that V∞ is determined we can again match the inner and outer 57 3.6. The effect of closely spaced patches solutions to get that U0 satisfies (3.15) subject to the singularity condition U0 ∼ V∞ + νγ log |x − x1 | + γ as x → x1 . (3.24) We remark that the constant d associated with ν = −1/ log(εd) and the solution Vc to (3.19) is as yet undetermined. Since we will show below that γ = 0, it follows that, for these two-patch cases, the value of d is not needed. By comparing the strengths of the logarithmic singularity, we can rewrite the problem for U0 in terms of a Dirac delta function as ∆U0 = 2πνγδ(x − x1 ) , x ∈ Ω \ {x1 } ∂U0 = 0, x ∈ ∂Ω , ∂n U0 ∼ V∞ + νγ log |x − x1 | + γ (3.25) as x → x1 . From the divergence theorem we conclude that γ = 0. Therefore, when the prey and predator patches are O(ε) close there is no singularity in the outer solution.This implies that the solution to (3.25) is simply U0 = V∞ . We now compare this very simple result of the two-close-patch asymptotic analysis with corresponding results computed from the full PDE (3.14) for the splitting probability. 58 3.6. The effect of closely spaced patches Figure 3.9: Full numerical solution for two closely spaced patches for case 2: Patch II is the predator and Patch I the prey. Parameter values as in the text. We choose a1 a2 =3 and β = 2(a1 + a2 ). Then, we calculate R1 and R2 from (3.20) and (3.21) as R1 = √ 3+ 5 2 ; R2 = 4 √ . (3+ 5)2 From (3.22) and (3.23) we compute from the asymptotic theory that the two values of V∞ are 2 3 and 1 3 in the first and second cases, respectively. The corresponding full numerical results from (3.14) (domain average of the solution) are 0.675 and 0.325, which agree very closely with the asymptotic results. The full numerical result for the splitting probability for case 2, for which patch II is the target and patch I is the predator, is shown in Fig. 3.9. Hence the asymptotic analysis is in good agreement with the numerical solution as the maximum error, 0.01, is O(ε) where ε = 0.0067. It is also interesting to note that when both circles have the same radii then, irrespective of the distance between the patches, we calculate that V∞ = 1/2 from (3.22) and (3.23). We can also plot the outer solution V∞ from (3.22) and (3.23) as a function of the radii a1 and a2 of the patches, and the distance between the patches, β. This leads to the plots as shown in Fig. 3.10. 59 3.6. The effect of closely spaced patches 0.9 0.8 v−infty 0.7 0.6 0.5 beta=2*(a1+a2) beta=3*(a1+a2) beta=4*(a1+a2) 0.4 0.3 0.2 0 2 4 6 8 10 a1/a2 0.8 0.75 a1/a2=1 a1/a2=3 a1/a2=6 v−infty 0.7 0.65 0.6 0.55 0.5 0.45 1.5 2 2.5 3 3.5 4 4.5 beta (values on axis are multiplying factors of (a1+a2)) 5 Figure 3.10: In the top figure, the constant outer solution is plotted against the ratio of radii of the patches for various distances between the patches. Here a1 and a2 are the radii of prey and predator patches respectively. In the bottom figure, we plot the constant outer solution versus the distance between the patches where the distances are represented on the horizontal axis as mutiplying factors of (a1 + a2 ). 60 3.6. The effect of closely spaced patches As we can see in the plots in the Fig. 3.10, as the radii ratio increases, the outer solution increases. Also for all the three different patch separations, the constant outer solution is exactly 1/2 when the radii ratio is unity. 3.6.2 Two close patches and a distant predator One more interesting problem occurs when we have one prey patch and one predator patch close to each other, with an additional predator patch located away from the two patches. This problem is a combination of the problems discussed in the preceding two sections and it is interesting to see what difference the presence of an extra predator in the distant field makes to the otherwise constant outer solution associated with the two closelyspaced patches. As defined in the previous section let Ωε1 and Ωε2 denote the two close patches, with one being a prey patch and the other a predator patch. These two patches are centred at x1 and x2 respectively. In addition to these two closely spaced patches, we assume that there is a distant patch Ωε3 centred at x3 . The analysis of the splitting probability leads to a PDE similar to (3.15) for the outer solution U0 with singularity condition at x1 , corresponding to the closely spaced patches, and at x3 corresponding to the distant or isolated patch. The matching condition at x1 , for the close patches, will be the same as that for the two closely-spaced patches case described in the previous section, The matching condition near the distant patch will be the same as the one described in the first section of this chapter. In this way, we obtain for the closely spaced patches that U0 ∼ V∞ + ν1 γ1 log |x − x1 | + γ1 as x → x1 . Alternatively, for the distant patch, U0 ∼ ν3 γ3 log |x − x3 | + γ3 as x → x3 . The problem for (3.15) for U0 with these singularity conditions can be 61 3.6. The effect of closely spaced patches written in terms of Dirac delta functions in Ω as ∆U0 = 2πν1 γ1 δ(x − x1 ) + 2πν3 γ3 δ(x − x3 ) , ∂U0 = 0, ∂n x ∈ Ω, (3.26) x ∈ ∂Ω . We then express U0 as a superposition of the Neumann Green’s function as U0 = 2πν1 γ1 G(x; x1 ) + 2πν3 γ3 G(x; x3 ) + χ . From the divergence theorem and from the matching conditions, we obtain the following three equations for γ1 , γ3 and χ: 2πν1 γ1 + 2πν3 γ3 = 0 γ1 (ν1 R(x1 ) − 1) + 2πν3 γ3 G(x1 ; x3 ) + χ = V∞ (3.27) γ3 (ν3 R(x3 ) − 1) + 2πν1 γ1 G(x3 ; x1 ) + χ = 0 , where R is the regular part of the Neumann Green’s function. In this system V∞ can be calculated by using conformal mapping and was given previously in (3.22) and (3.23). For the distant patch, ν3 = −1/ log(εd3 ), where d3 is known for a circular patch. The only new quantity that needs to be calculated is d1 , the logarithmic capacitance associated with two closely-spaced patches. Recall that it is obtained from the farfield behavior of the solution to (3.19). For the special case when the two closely spaced circular patches have a common radius a, this logarithmic capacitance of the two-disk cluster was calculated in Appendix B of [4] as log d1 = log (2ζ) − ξc + 2 ∞ m=1 e−mξc , m cosh(mξc ) (3.28) where half-separation between the disk centers, l = β/2 (measured in the local y coordinate) while ζ and ξc are determined in terms of the disk radius 62 3.6. The effect of closely spaced patches a and l by ζ= l 2 − a2 ; cosh(ξc ) = l , a l ξc = log + a l a 2 − 1 . (3.29) The series (3.28) converges fast and we obtained that d1 = 2.0613. In this way, all of the quantities in (3.27) can be evaluated, and the resulting system solved for γ1 , γ3 , and χ. Figure 3.11: Full numerical solution for two closely spaced patches and a distant predator. We now compare the results predicted by this asymptotic theory with full numerical results. We consider an example for which the prey patch is located at at (0.25, 0) and the two predator patches located at (0.2768, 0) and (0, 0.6) such that the distance between the close patches is 0.0268 = 4ε. Each patch is taken to a circular disk of radius ε = 0.0067. For this scenario, our full numerical results for the spatial average of the splitting probability is computed as 0.28973, while the asymptotic spatial average is χ = 0.28905. This gives an error of 0.0007, which is O(ε) where ε = 0.0067. The results of the numerical simulation is shown in Fig. 3.11. 63 Chapter 4 Mean First Passage Time with Variable Diffusivity 4.1 Introduction In any habitat, the presence of terrestrial inhomogenities like hills, mountains, lowlands, marshes, and rivers can affect animal movement in the habitat. The simplest model of mean first passage time that accounts for these inhomogenities in the landscape is to allow for a variable diffusivity in the PDE for the MFPT without drift. In this chapter, we investigate a scenario in which this diffusivity varies as a function of the distance from the centre of a circular landscape. Before specializing the results to this specific scenario, we show how the previous MFPT analysis can be modified to allow for an arbitrary variable diffusivity D in an arbitrary 2-D domain. The non-dimensional PDE for the MFPT can be written as N ∆u(x) = −1/D(x) , Ωεj , x∈Ω\ j=1 ∂u = 0, ∂n u = 0, (4.1) x ∈ ∂Ω , x ∈ ∂Ωεj , j = 1, . . . , N . The problem is the same as the one without drift analyzed in the second chapter, with the only difference that the diffusivity D here is spatially dependent. Most of the analysis performed is similar to the case with constant diffusivity. First, we will analyse the situation with one prey patch and then 64 4.2. Asymptotic solution for one patch generalize it to the case of N patches. 4.2 Asymptotic solution for one patch The MFPT on domain with one patch satisfies the following PDE x ∈ Ω \ Ωε1 , ∆u(x) = −1/D(x) , ∂u = 0, ∂n u = 0, (4.2) x ∈ ∂Ω , x ∈ ∂Ωε1 . The outer expansion has the form U ∼ U0 (x, ν) + σ(ε)U1 (x, ν) + . . . . Here σ is beyond-all-orders with respect to the logarithmic terms, and U0 satisfies ∆U0 = −1/D(x) , ∂U0 = 0, ∂n x ∈ Ω \ {x1 } , (4.3) x ∈ ∂Ω , subject to a singularity condition as x → x1 . The inner expansion has the form q ∼ νγ(ν)q0 (y) + α(ε, ν)q1 (y) + . . . , (4.4) so that q0 satisfies ∆ y q0 = 0 , q0 = 0 , y∈ / Ω1 where Ω1 = Ωε1 /ε , y ∈ ∂Ω1 , q0 = log |y| − log d + O(|y| (4.5) −1 ) as |y| → ∞ . From the far-field behavior of the inner solution, we obtain that U0 must satisfy (4.3) subject to the singularity condition U0 ∼ νγ(ν) log |x − x1 | + γ(ν) , as x → x1 , (4.6) 65 4.2. Asymptotic solution for one patch where ν ≡ −1/ log(εd1 ). The problem (4.3) with singularity strength as in (4.6) is equivalent to ∆U0 = − 1 + 2πνγ(ν)δ(x − x1 ) , D ∂U0 = 0, ∂n x ∈ Ω, (4.7) x ∈ ∂Ω . The divergence theorem then yields the constraint that 2πνγ = Ω 1 dx , D(x) (4.8) which determines the γ. Under the assumption of variable diffusivity the solution of (4.7) is represented in terms of the decomposition U0 = U0H + 2πνγG(x; x1 ) + χ , (4.9) where χ is a constant, and G(x; ξ) is the Neumann Green’s function satisfying ∆G = − 1 + δ(x − ξ) , |Ω| ∂G = 0, ∂n x ∈ Ω, x ∈ ∂Ω , (4.10) Gdx = 0 , Ω G= 1 (log |x − ξ| + R(ξ)) + o(1) 2π as x→ξ The smooth part U0H of the decomposition in (4.9) is taken to be the unique solution of ∆U0H = − 1 + D(x) ∂U0H = 0, ∂n 1 |Ω| x ∈ ∂Ω , Ω 1 dx , D(x) x ∈ Ω, (4.11) U0H dx = 0 . Ω 66 4.3. Solution for multiple patches In addition, we can interpret χ as the spatial average of U0 χ= 1 |Ω| U0 dx . (4.12) Ω Finally, we expand the solution U0 from (4.9) as x → x1 and equate the resulting behavior with the required singular behavior (4.6). In this way, we obtain that χ is given by χ = γ(1 − νR(x1 )) − U0H (x1 ) , (4.13) where γ is given in (4.8). Depending on the choice of D(x) and the shape of the domain, we can solve for U0H either analytically or numerically. In terms of the Neumann Green’s function and its regular part, and the smooth solution U0H , we can use the γ as determined in (4.8) to evaluate χ in (4.13). This determines U0 from the decompositon (4.9). 4.3 Solution for multiple patches In a similar way, we can obtain the approximate solution for the case of N patches. The inner problem near the patches is precisely the same and the corresponding outer solution admits the decomposition N U0 = U0H + 2π νk γk (ν)G(x; xk ) + χ , (4.14) k=1 where U0H satisfies (4.9) and G is the Neumann Green’s function. Here νk = −1/ log(εdk ), where dk is the logarithmic capacitance of the k-th patch. The first constraint on the γk for k = 1, . . . , N , as obtained by the divergence theorem, is that N 2π νk γk (ν) = k=1 Ω 1 dx . D(x) (4.15) 67 4.4. Numerical verification From the regular part of the singularity conditions, we obtain the additional N constraints that N νk γk (ν)G(xj ; xk )+χ = −U0H (xj ) , (4.16) γj (ν)(νj (Rj (xj ))−1)+2π k=1,k=j for each j = 1, . . . , N . The system (4.15) and (4.16) consists of N + 1 linear algebraic equations for the unknowns χ and γj for j = 1, . . . , N . 4.4 Numerical verification In this subsection, we compare results predicted by our asymptotic theory with full numerical results obtained by solving the MFPT PDE with variable diffusivity using Matlab. In our computations we took a circular domain of radius 1km. The prey patches are all circles of radius 0.0067km, so that ε = 0.0067. For simplicity, in order to be able to solve the problem for U0H analytically, we chose a variable diffusivity function that is radially symmetric of the form D= 1 . 1 + |x|2 (4.17) Thus, the diffusivity decreases with increasing distance from the center of the domain. One key advantage of assuming a circular-shaped landscape is that the Neumann Green’s function is explicitly available. In particular, when Ω is the unit disk, and writing x and ξ as points in complex notation inside the unit disk, it is well-known that the solution to (4.10) is (cf. [11]) 1 2π 1 R(ξ) = − 2π G(x; ξ) = − 1 3 ξ + (|x|2 + |ξ|2 ) − |ξ| 2 4 ξ 3 − log ξ|ξ| − + |ξ|2 − . |ξ| 4 − log |x − ξ| − log x|ξ| − , In this chapter and the next one, these formulae are used to calculate the Neumann Green’s function, G and its regular part, R. For the choice of D (4.17), the problem (4.11) for U0H admits a solution 68 4.4. Numerical verification of the form U0H (r), with r = |x|, where 1 1 U0H + U0H = − r2 . r 2 The solution to this problem with U0H = 1 0 U0H r dr = 0 is simply r4 1 r2 − − . 8 16 19 The results from the asymptotic theory and the full numerical approach are given in terms of the MFPT averaged over all possible starting points on the domain (i.e. spatial averages of MFPT), where the initial starting point has a uniform distribution: Tavg = 1 |Ω| T (X) dX . (4.18) Ω More details on the procedures used in the numerical experiments are available in the Appendix. 4.4.1 One patch case We first study the effect of patch location in the domain on the MFPT of the predator. For this purpose we consider the situation where the domain contains only one patch, and we will vary the patch distance from the centre. 69 4.4. Numerical verification Figure 4.1: Spatial averages of the asymptotic and full numerical results for the MFPT are plotted versus the distance of the patch from the centre. The variable diffusivity has the form in (4.17). The asymptotic and numerical results are shown in Fig. 4.1. From this figure we observe that the numerical and asymptotic solutions are in close agreement. The maximum error between the asymptotic and full numerical results is 0.019 which is O(ε), where ε=0.0067. As to be anticipated from the monontone decreasing nature of D = D(|x|), the value of the average MFPT increases with the distance of the patch from the centre. 4.4.2 Multiple patches case Finally, we compare the numerical and asymptotic results when the circular domain contains three patches. For three spatial configurations of the three patches the results are shown in Table 4.1. These results again confirm the asymptotic theory. 70 4.4. Numerical verification Patch1 (0.1,0) (0.5,0.3) (0.5,0) Patch2 (-0.3,-0.7) (-0.2,0.6) (-0.3,-0.7) Patch3 (0,0.95) (-0.4,-0.7) (-0.6,0.4) avg. MFPT (hr) (asy) 1.1526 1.0589 1.0289 avg. MFPT (hr) (num) 1.1596 1.0364 1.0052 Table 4.1: Spatial averages of MPFT in a domain with three patches, obtained using both the asymptotic theory and full numerical simulations. The variable diffusivity is given in (4.17). 71 Chapter 5 Second Moment Analysis and Estimation of Variance 5.1 Introduction In a previous chapter we developed a singular perturbation method to asymptotically calculate the mean first passage time (MFPT) for a predator to catch a prey. In general, MFPT estimation is important for determining the search time of a predator. However, it may also be important to estimate second-moment information, such as the variance, when the first passage time has a significant spread around the mean. A simple PDE problem for the variance of the first passage time is not available. Instead, one needs to determine the MFPT and the second moment of the first passage time (SMFPT) in order to calculate the variance of the first passage time (VFPT). It is well-known that the SMFPT of a predator in a two-dimensional domain satisfies the following PDE in dimensional form ([15]): D∆ M (X) + 2T = 0 . (5.1) Here D is the constant diffusivity, T is the mean first passage time, and M is the second moment of first passage time. The PDE is made dimensionless by taking w = D2 M L4 and x = L−1 X. The primed operators are defined with respect to X, while the un-primed ones are with respect to the non-dimensional variable x. In addition, 2L is the characteristic length of the domain under consideration. Also we take T = D T L2 , which is a non-dimensionalization in accordance with the analysis of the MFPT in the 72 5.2. Asymptotic solution second chapter. The corresponding non-dimensionalized PDE for the second moment takes the form N ∆w(x) = −2T , Ωεj , x ∈ Ω\ j=1 ∂w = 0, ∂n w = 0, (5.2) x ∈ ∂Ω , x ∈ ∂Ωεj . This PDE is similar to the MFPT PDE, and it can be solved asymptotically by using an approach similar to that used to calculate the MFPT. The results from the asymptotic analysis are then verified with full numerical simulations of the PDE for a few test cases. 5.2 Asymptotic solution We now construct the solution to (5.2) asymptotically. The outer expansion has the form w ∼ W0 (x, ν) + σ(ε)W1 (x, ν) + . . . , where σ is assumed to be beyond-all-orders with respect to the logarithmic terms. The problem for W0 is ∆W0 = −2T , ∂W0 = 0, ∂n x ∈ Ω \ {x1 , . . . , xN } , (5.3) x ∈ ∂Ω , with singularity conditions to hold as x → xj for j = 1, . . . , N . Up to negligible O(ε) terms, the inner problem near the j-th patch is identical to the one obtained for the MFPT without drift. By using the far-field behaviour of the inner solution, we match the inner and outer solutions to derive that the solution W0 to (5.3) must have the singularity behavior W0 ∼ νj γj (ν) log |x − xj | + γj (ν) as x → xj , (5.4) 73 5.2. Asymptotic solution where νj = −1/ log(εdj ) and dj is the logarithmic capacitance of the j-th patch. Next, by incorporating the correct strength of the logarithmic singularity, we can write the outer problem (5.3) in terms of Dirac delta functions as N ∆W0 = −2T + 2π νk γk (ν)δ(x − xk ) , x ∈ Ω, (5.5) k=1 ∂W0 = 0, ∂n x ∈ ∂Ω . From the MFPT analysis we know that T has the outer asymptotic expansion N T ∼ 2π νk Γk (ν)G(x; xk ) + χ , (5.6) k=1 where Γk for k = 1, . . . , N are constants and the Neumann Green’s function, G(x; ξ), satisfies ∆G(x; ξ) = − 1 + δ(x − ξ) , |Ω| ∂G = 0, ∂n x ∈ Ω, x ∈ ∂Ω , (5.7) G dx = 0 , Ω G∼ 1 (log |x − ξ| + R(ξ)) + o(1) 2π as x → ξ . Next, by substituting (5.6) into (5.5), we obtain that N ∆W0 = −4π N νk Γk (ν)G(x; xk ) − 2χ + 2π k=1 ∂W0 = 0, x ∈ ∂Ω , ∂n W0 ∼ νj γj (ν) log |x − xj | + γj (ν) νk γk (ν)δ(x − xk ) , x ∈ Ω, k=1 as x → xj ; j = 1, . . . , N . (5.8) To obtain the first constraint on the unknown γj for j = 1, . . . , N , we 74 5.2. Asymptotic solution integrate (5.8) and use the divergence theorem to get π |Ω| N νk γk (ν) = χ . (5.9) k=1 Then, the problem (5.8) can be conveniently decomposed into two terms as W0 = W0H + W0p , (5.10) where W0p is given by N W0p = −2π νk Γk (ν)W0pk , (5.11) k=1 and where W0pk satisfies ∆W0pk = G(x; xk ) , ∂W0pk = 0, ∂n x ∈ Ω, x ∈ ∂Ω , (5.12) W0pk dx = 0 . Ω In contrast, the other term, W0H , satisfies N νk γk (ν)δ(x − xk ) , ∆WOH = −2χ + 2π x ∈ Ω, k=1 ∂WOH = 0, x ∈ ∂Ω , ∂n W0H ∼ −W0p (xj ) + νj γj (ν) log |x − xj | + γj (ν) as x → xj . (5.13) The solution to (5.13) can be written as N W0H = 2π νk γk (ν)G(x; xk ) + χ0H , (5.14) k=1 where G is the Neumann Green’s function defined earlier in (5.7) and χ0H 75 5.2. Asymptotic solution is an unknown constant. This constant can be interpreted as the spatial average of W0H and W0 1 |Ω| χ0H = since Ω W0p (x) dx W0H (x) dx =⇒ χ0H = Ω W0 (x) dx , Ω = 0. In order to verify that (5.14) satisfies (5.13), we calculate from the constraint (5.9) that N ∆W0H = 2π N νk γk (ν)∆G(x; xk ) = 2π k=1 N k=1 νk γk δ(x − xk ) − = 2π νk γk (− k=1 N 2π |Ω| 1 + δ(x − xk )) , |Ω| N νk γk k=1 νk γk δ(x − xk ) − 2χ , = 2π k=1 (5.15) which agree exactly with (5.13). Finally, we expand W0 as x → xj , by using (5.10), (5.11), and (5.14), and compare the result with the required singularity behavior in (5.4). This yields that N νk γk (ν)G(xk ; xj ) + νj γj (ν)(log |x − xj | + Rj (xj )) + χ0H 2π k=1,k=j ∼ νj γj (ν) log |x − xj | + γj (ν) − W0p (xj ) . Since the logarithmic terms agree automatically, we obtain from the remaining terms that the γj for j = 1, . . . , N and the constant χ0H must satisfy N γj (ν)(νj Rj (xj ) − 1) + 2π νk γk (ν)G(xj ; xk ) = −χ0H − W0p (xj ) , k=1,k=j (5.16) 76 5.3. Regularization scheme of the function, W0pk for j = 1, . . . , N together with the constraint (5.9). In setting up this linear system we must calculate W0p (xj ), where W0p is defined in (5.11) in terms of the solution W0pk to (5.12). Although the Neumann Green’s function, representing the right-hand side of (5.12), is known explicitly for a circular or rectangular domain, it is still difficult to analytically calculate the solution to (5.12) in an explicit tractable form. As such, we employ a numerical regularization, similar to the one used in second chapter to regularize Green’s function (2.7), to compute the solution W0pk to (5.12). 5.3 Regularization scheme of the function, W0pk Similar to the scheme used in second chapter, we consider an auxiliary function φ(x, ξ) that satisfies the following PDE: ∆φ − µφ = G , ∂φ = 0, ∂n x ∈ Ω, (5.17) x ∈ ∂Ω , Here µ is a small positive parameter. We can use an approach similar to the one used in Green’s function regularization (2.7) to prove that (5.17) has unique solution for µ > 0. As µ → 0, the solution to the PDE (5.17) can be expressed as the following asymptotic series: φ∼ φ0 + φ1 + µφ2 + µ2 φ3 + · · · . µ (5.18) Upon substituting this expansion into the PDE (5.17) and equating coefficients of µ at each order, we obtain a set of PDEs for each order. For the O(µ−1 ) term, we obtain that φ0 satisfies ∆φ0 = 0 , x ∈ Ω, ∂φ0 = 0, ∂n x ∈ ∂Ω , (5.19) 77 5.3. Regularization scheme of the function, W0pk with φ0 smooth in the domain. The solution to this problem is that φ0 is an arbitrary constant. Similarly, from the O(1) term we obtain the PDE: x ∈ Ω, ∆φ1 = φ0 + G , ∂φ1 = 0, ∂n (5.20) x ∈ ∂Ω , To solve for φ1 we must have, from the divergence theorem, that Ω φ0 dx + Ω Gdx = 0 ⇒ φ0 = 0 (∵ Ω Gdx = 0) . Similarly, from the O(µ) terms, we obtain the following PDE for φ2 : ∇(P ∇φ2 ) = φ1 , x ∈ Ω, (5.21) ∂φ2 = 0, ∂n x ∈ ∂Ω , where φ2 is pointwise bounded in Ω. From the divergence theorem, we obtain that the solvability condition for the PDE for φ2 is Ω φ1 dx =0 Now using this constraint we rewrite the PDE for φ1 as ∆φ1 = G , x ∈ Ω, ∂φ1 = 0, ∂n x ∈ ∂Ω , (5.22) φ1 dx = 0 , Ω The PDE for φ1 is same as the PDE for W0pk to be found. Hence the O(1) term in the expansion (5.18) of the solution to the regularized PDE (5.17) is the required function. The solution φ, to (5.17), can be calculated numerically using conventional techniques such as the finite element method. Once we have φ for two small values of µ, we can use Richardson extrapolation to approximate φ1 and hence the Green’s function with an accuracy of O(µ2 ). 78 5.4. Numerical verification 5.4 Numerical verification In this section we numerically verify the results of the asymptotic analysis with full numerical solutions of the PDE (5.2) computed using Matlab. To find SMPFT numerically using Matlab, we solve the following coupled PDEs: one PDE is that SMFPT (5.2) and the other one is the MFPT PDE from second chapter. N ∆w(x) = −2T ; ∆T (x) = −1 , Ωεj , x ∈ Ω\ j=1 ∂T ∂w = 0; = 0, ∂n ∂n w = 0; T = 0 , (5.23) x ∈ ∂Ω , x ∈ ∂Ωεj . While evaluating and comparing the numerical and asymptotic solutions, we solve the non-dimensional version of the PDE to get the non-dimensional SMFPT. The results are, then, rescaled by a factor of (L4 /D2 ), to obtain the SMFPT with dimensions. Here D is the diffusivity of predator in the domain and is taken to be 0.41 km2 /hr. In the numerical experiments below, we use a circular domain of radius 1km. The prey patches are all circles of radius 0.0067 km, so that ε= 0.0067. The exact analytical expressions given for G and R in previous chapter are used in evaluating the asymptotic solution. The results from the asymptotic theory and the full numerical simulations are reported in terms of spatial averages of the second moment as Mavg = 1 |Ω| M (X) dX . (5.24) Ω This indirectly is the distribution average of MFPT, in the case when the probability density of the predator follows a uniform distribution in the domain. We then use the SMFPT calculated here and the MFPT calculated in the second chapter to estimate the variance of first passage time (VFPT). Variance can be calculated using V F P T = SM F P T − M F P T 2 The asymptotic and numerical solutions are compared in terms of the 79 5.4. Numerical verification spatial average of VFPT. More details on the procedures used in the numerical experiments are available in the Appendix. 5.4.1 A single patch in a circular domain We first study the effect of patch location in the domain on the SMFPT and VFPT. For this purpose we consider the situation where the domain contains only one patch, and we compare asymptotic and full numerical results as the patch moves away from the centre of the circle. The results are shown in Fig. 5.1. Figure 5.1: Spatial averages of analytical and numerical solutions of second moment are plotted against distance of the patch from the centre From this figure we observe that the full numerical and the asymptotic results are in very close agreement with each other. The relative error between the numerical and asymptotic solution average ranged from 0.02% to 0.05%, which is O(ε) (= 0.67%). The value of the average SMFPT increases with the distance of the patch from the centre. Next, we use the results for the MFPT as obtained from Principal Result 2.1 of Chapter 2 to plot the variance of first passage time. The results are shown in Fig. 5.2. 80 5.4. Numerical verification Figure 5.2: Spatial averages of analytical and numerical solutions of the variance are plotted against the distance of the patch from the centre From this figure we can see that the numerical and asymptotic solutions are still in good agreement with each other. The value of the average VFPT also increases with the distance of the patch from the centre. 5.4.2 Multiple patches in a circular domain Here we compare the numerical and asymptotic solution averages when the domain contains three patches. The results are shown in Table 5.1 for two configurations of three patches. Patch1 (0.5,0.3) (0.3,0.8) Patch2 (-0.2,0.6) (0.1,-0.6) Patch3 (-0.4,-0.7) (-0.5,-0.7) avg. SMFPT (hr2 ) (asy) 5.6884 7.6893 avg. SMFPT (hr2 ) (num) 5.7081 8.1220 Table 5.1: Spatial averages of SMPFT in a domain with three patches, obtained using both the asymptotic theory and full numerical simulations. From this table, we observe that, for two different arrangements of the three patches in the domain, the asymptotic result for the SMFPT is in good agreement with the full numerical solution. 81 5.4. Numerical verification Using the SMPT and MFPT values, asymptotic and full numerical results for the average of the variance is shown in Table 5.2 for the same two configurations of three traps. Patch1 (0.5,0.3) (0.3,0.8) Patch2 (-0.2,0.6) (0.1,-0.6) Patch3 (-0.4,-0.7) (-0.5,-0.7) avg. VFPT (hr2 ) (asy) 2.8398 3.7206 avg. VFPT (hr2 ) (num) 2.8579 4.1476 Table 5.2: Spatial averages of VFPT in a domain with three patches, obtained using both the asymptotic theory and full numerical simulations. 82 Chapter 6 Conclusions The first chapter in the thesis stresses the need for asymptotic and other analytical tools in the field of mathematical ecology where complicated stochastic models are being implemented. In the second chapter, we introduce the concept of the mean first passage time (MFPT) and the implementation of hybrid asymptotic technique to solve the Kolmogorov equation for the MFPT. Using these analytical techniques along with numerical verification, we reproduce the results similar to those obtained in [12] using only full numerical simulations. It was shown that the MFPT increases with patch distance from the centre of a circular landscape. We also show that the asymptotic approach yields results that are in close agreement with those obtained by numerical methods, even though ε is not very small. Here, we also explain reason for such increase in MFPT by arguing that as patch moves away from the centre, the average distance of the patch from all the points on the domain also increases, leading to increase in average MFPT. We build on these results to calculate the MFPT for multiple patches in the domain. We also show that when the drift velocity is very small, MFPT will be approximately the same as for the case without drift. However, for larger centralizing speeds, for every patch location, there exists a critical point at which random diffusion is more advantageous than diffusion with drift. Future work can be focussed finding an analytical estimation of this critical drift value for various patch distributions. Another open problem would be to calculate the MFPT in a two-dimensional domain without drift where the patch concentrates not at a point but along some finite O(1) length curve of small width contained within the domain. The recent asymptotic analysis of [9] is relevant to such 83 Chapter 6. Conclusions a situation. The third chapter deals with the asymptotic and numerical estimation of splitting probability. In this chapter, we show that the chances of prey reaching target are smaller when a target is surrounded by a ring of predator patches, as compared to the situation when the target is located outside the ring of patches. We also investigate the effect of the starting location of the prey on the splitting probability. The analysis shows that the farther the starting points of the prey are from the ring of predators, the higher the splitting probability. We prove that when a prey patch and a predator patch are located close to each other, it is possible to analytically estimate the splitting probability and we show analytically that it depends on the ratio of the patch sizes and distance between the patches. The outer solution in this situation is a constant and can be evaluated using only the inner solution. When one more predator patch is located far away from these two close patches, the outer solution is no longer a constant. Future work on this subject can be focussed on finding optimal positioning of given predators so that the average splitting probability of the prey is the least in a domain. This can help us estimate the extinction limit of a prey species in a landscape. In fourth chapter on variable diffusivity, we extend our analysis on the MFPT in the second chapter to allow for a spatially variable diffusivity. Here we perform the analysis for a simple test function of diffusivity, D, which varies radially in a circular landscape. For this scenario, the asymptotic results are readily implemented analytically. An implementation of the analysis using various functions of variable diffusivity, to allow for more physical relevant inhomogenities such as highlands, rivers, hills etc., can be performed in the future. A further extension of this work would be to allow for a non-isotropic diffusivity modeled by a diffusion matrix (cf. [12]), or to consider the homogenization theory limit whereby the diffusivity varies periodically on some microscale (cf. [8]). 84 Chapter 6. Conclusions The fifth chapter concerns the calculation of the second moment of the first passage time. We used the second moment evaluated in this chapter and MFPT in second chapter to estimate the variance of first passage time (VFPT) in various scenarios. Point-wise analyses of the SMFPT and MFPT in the domain can yield new insights into the best hiding spots for the prey patches in the domain. The ultimate goal of this work is to introduce the application of the powerful technique of a hybrid-asymptotic method to compute asymptotic solutions of physically relavant quantities, in first passage processes, in the field of mathematical ecology and motivate its implementation in the study of more ecological problems in the future. Despite the fact that we have largely considered simplified model problems in our choices of the landscape, variable diffusivity etc., we fully believe that our hybrid method can effectively handle more biologically relavant problems, typical of those encountered by a mathematical ecologist, and we look forward to more such implementations in the future. 85 Bibliography [1] J. Wagner C. Fall, E. Marland and J. Tyson. Computational Cell Biology. Springer, 2002. [2] C. Chevalier, O. Benichou, B. Meyer, and R. Voituriez. First passage quantities of brownian motion in a bounded domain with multiple targets: a unified approach. J. Phys A: Math. Theor., 44:025002 (24 pp), 2011. [3] A. Cheviakov and M. J. Ward. Optimizing the fundamental eigenvalue of the laplacian in a sphere with interior traps. Mathematical and Computer Modeling, 53:1394–1409, 2011. [4] D. Coombs, R. Straube, and M. J. Ward. Diffusion on a sphere with localized traps: Mean first passage time, eigenvalue asymptotics, and fekete points. SIAM J. Appl. Math., 70(1):302–332, 2009. [5] W. W. Dijkstra and M. E. Hochstenbach. Numerical approximation of the logarithmic capacity. CASA Report 08-09, Technical U. Eindhoven, 2008. [6] A.E. Sowers E. Neumann and C.A. Jordan. Electroporation and Electrofusion in Cell Biology. Plenum, 1989. [7] P. Fauchauld and R. Tveraa. Using first passage time in the analysis of area-restricted search and habitat selection. Ecology, 84:282288, 2003. [8] B. Froese and A. Oberman. Numerical averaging of non-divergence structure elliptic operators. Commum. Math. Sci., 7(4):785–804, 2009. 86 Bibliography [9] R. Gadylshin and A. Ilin. On the spectrum of the neumann problem for laplace’s equation in a domain with a narrow slit. Asymptotic Anal., 67(3-4):785–804, 2010. [10] J. Kevorkian and J. Cole. Multiple Scale and Singular Perturbation Methods. Applied Mathematical Sciences, Vol. 114, Springer-Verlag, 1953. [11] T. Kolokolnikov, M. Titcombe, and M. J. Ward. Optimizing the fundamental neumann eigenvalue for the laplacian in a domain with small traps. Europ. J. Appl. Math., 16(2):161–200, 2005. [12] H. McKenzie, M. Lewis, and E. Merrill. First passage time analysis of animal movement and insights into the functional response. Bulletin of Mathematical Biology, 71:107–129, 2009. [13] S. Pillay, M. J. Ward, A. Peirce, and T. Kolokolnikov. An asymptotic analysis of the mean first passage time for narrow escape problems: Part i: Two-dimensional domains. SIAM J. Multiscale Modeling, 8(3):803– 835, 2010. [14] T. Ransford. Potential Theory in the Complex Plane. London Mathematical Society Student Texts Vol. 28, Cambridge University Press, 1995. [15] S. Redner. A Guide to First-Passage Processes. Cambridge University Press, 2001. [16] A. Singer, Z. Schuss, and D. Holcman. Narrow escape, part ii: The circular disk. J. Stat. Phys., 122(3):465–489, 2006. [17] J. Skellam. Random dispersal in theoretical populations. Biometrika, 38:196218, 1951. [18] R. Straube, M. J. Ward, and M. Falcke. Reaction rate of small diffusing molecules on a cylindrical membrane. J. Stat. Phys., 129(2):377–405, 2007. 87 [19] M. Titcombe and M J. Ward. Summing logarithmic expansions for elliptic equations in multiply-connected domains with small holes. Canadian Appl. Math. Quart., 7(3):313–343, 1999. [20] M. Titcombe and M J. Ward. An asymptotic study of oxygen transport from multiple capillaries to skeletal muscle tissue. SIAM J. Appl. Math., 60(5):1767–1788, 2000. [21] M. J. Ward, W. D. Henshaw, and J. B. Keller. Summing logarithmic expansions for singularly perturbed eigenvalue problems. SIAM J. Appl. Math., 53(3):799–828, 1993. 88 Appendix: Numerical Methods All numerical experiments were performed in Matlab using the PDE Toolbox. The toolbox uses finite element methods to solve the two dimensional elliptic PDEs. Programming was done using a combination of GUI (Graphical User Interface) and command-line functions. Mesh implementation The domains are drawn by hand using disc shapes provided in the GUI. For numerical verification the mesh generated is refined 3–4 times to obtain at least four digits of accuracy in the results. Further refining, will lead to poor conditioning and memory problems and also prompts Matlab to give the “out-of-memory” error. Typical asymptotic results are in agreement up to three to four digits. The mesh used to find the numerical solution changes with the problem. The mesh generated is triangular; the numerical mesh has close to 115000 nodes (obtained after three refinements) while the numerical mesh after four refinements has about 529000 nodes and almost twice as many triangles. Numerical solution To find the numerical solution, the domain is taken to be a unit circle with circular patches removed from it. After specifying the PDE and boundary conditions, the mesh is generated and refined. The PDE is then solved with the help of GUI and command line functions to generate the results and 89 Asymptotic solution plots of the numerical solution. Command line functions are used to extract solution data from the figure file generated by GUI. This solution data is then used to re-plot figures, calculate averages etc. Asymptotic solution For the asymptotic solution, a circle is used as the domain, and then the mesh is generated and refined. The outer solution (in most cases, an analytical formula for the outer solution is available) is evaluated at each of the mesh points and plotted. The same mesh is used for outer solution calculation for all of the problems considered in this thesis. Evaluating regular part of a Green’s function Let G(x, x0 ) be a Green’s function with singularity at x0 . Let x be the node/grid point closest x0 on the mesh. Then regular part of the Green’s function, R(x0 ) is approximated as R(x0 ) ≈ 2πG(x , x0 ) − log |x − x0 | (1) Regularization scheme of Green’s function The implementation of the regularization scheme was a tricky task due to the presence of singularities (Dirac delta functions) within the domain. In the finite element method used in the PDE toolbox, every known function in the PDE (including delta function) is assumed to a constant over each element of the mesh. The constant value is taken to be the value of the function at the centeroid of the element. However, the delta function is zero at all points except the point of singularity where it is infinity. So the toolbox normally takes delta function to be a zero function since the singularity is almost never exactly centred at the centeroid. Taking zero or infinity on the element will not give the required results. Understanding 90 Regularization scheme of Green’s function how to overcome this problem was a long cumbersome task and it consumed considerable time to figure out an optimal scheme. Dirac delta function implementation We tried three different approaches. The first one was to approximate the delta function by a sharply peaked normal distribution function. For maximum accuracy the function must be centred on the centeroid of the element with singularity. In Matlab, since each point on the domain is defined up to more than 16 decimals (e.g 2.5000000000405041), precise centering of the peak on a grid points is not possible. The second approach was to manually edit the matrices, in the final algebraic system to be solved, in the finite element scheme. The matrix entries are changed to values that one would obtain theoretically, without the assumption that the known functions are constants over elements. This method proved to be better than the first one. However, its implementation is complicated and prone to errors. Optimal method This prompted us to look for a third method of handling the delta function, which turned out to be the best of the methods we tried. In this method we analytically take care of the singularity by taking G = 1 2π log |x − ξ| + B, where B is a regular smooth function. Starting with the PDE for G we have, P (ξ) + P (x)δ(x − ξ) |Ω| 1 =⇒ ∆G + ∇ψ · ∇G = − + δ(x − ξ) (∵ P = eψ ) P |Ω| ∇ · (P ∇G(x; ξ)) = − If we write G = 1 2π log |x − ξ| + B, then the PDE for the smoother part B takes the form ∆B + ∇ψ · ∇B = − 1 1 ∇ψ · ∇ log |x − ξ| − . 2π P |Ω| (2) This PDE can be solved directly using the PDE toolbox without having to 91 Calculation of averages worry about the delta function. The Green’s functions in the second and third chapters are obtained using the above regularization scheme. It turned out that Richardson extrapolation did not improve accuracy in the second chapter, so we chose not to use it to improve the computational speed. On performing various iterations, it was found that µ must the O(10−2 ) to get optimal results. So in the numerical experiments we use µ= 0.01 or 0.02. In second chapter the results of the regularization scheme are verified against the exact solution available for the case when the patch is at the centre and we had an accuracy up to fourth decimal. Since, we know the exact solution for the Neumann Green’s function for a circular domain, we use it directly in the fourth and fifth chapters. Note that the PDE toolbox is unable to perform the numerical experiments when a patch is exactly at the centre of the domain (0,0). So we had to slightly shift the centre of the patch to (0,0.001) to perform the experiments using the PDE toolbox. Calculation of averages The following integrals that are used to calculate spatial and distribution averages of the functions T are estimated using Reimann sum on mesh elements: Tavg = 1 |Ω| Tavg = T (X) dX (3) Ω S(X)T (X) dX (4) Ω Second moment calculations In the second moment calculations, we had to use a system of coupled elliptic PDEs: one PDE in the system is the Kolmogorov equation for the mean first passage time and other one is the PDE for second moment with the first term appearing on the right-hand side of the operator. To calculate W0pk we used 92 Second moment calculations a regularization scheme explained in fifth chapter. Richardson extrapolation was then used in µ to estimate W0pk . 93
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Asymptotic analysis of first passage processes : with...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Asymptotic analysis of first passage processes : with applications to animal movement Venu, Kurella 2011
pdf
Page Metadata
Item Metadata
Title | Asymptotic analysis of first passage processes : with applications to animal movement |
Creator |
Venu, Kurella |
Publisher | University of British Columbia |
Date Issued | 2011 |
Description | Understanding the dependence of animal behaviour on resource distribution is a central problem in mathematical ecology. In a habitat, the distribution of food resources and their accessibility from an animal's location together with the search time involved in foraging, all govern the survival of a species. In this work, we investigate various scenarios that affect foraging habits of animals in a landscape. The work, unlike previous studies, analyzes the first passage quantities on complex prey-predator distributions in a given domain in order to derive simple analytical problems that can readily be solved numerically. We use standard stochastic models such as the Kolmogorov equations of first passage times and splitting probability, to model both the foraging time of a predator and the chances of survival of prey on a landscape with prey and predator patches. We obtain an asymptotic solution to these Kolmogorov equations using a hybrid asymptotic-numerical singular perturbation technique that utilizes the fact that the ratio of the size of prey patches is small in comparison to the overall landscape. Results from this hybrid approach are then verified by undertaking full numerical simulations of the governing partial differential equations of the first passage processes. By using this hybrid formulation we identify the underlying parameters that affect the search time of a predator and splitting probability of prey, which are otherwise difficult to ascertain using only numerical tools. This analytical understanding of how parameters influence the first passage processes is a key step in quantifying foraging behavior in model ecological systems. |
Genre |
Thesis/Dissertation |
Type |
Text |
Language | eng |
Date Available | 2011-08-29 |
Provider | Vancouver : University of British Columbia Library |
Rights | Attribution-NonCommercial-NoDerivatives 4.0 International |
IsShownAt | 10.14288/1.0072070 |
URI | http://hdl.handle.net/2429/36961 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2011-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Rights URI | http://creativecommons.org/licenses/by-nc-nd/4.0/ |
AggregatedSourceRepository | DSpace |
Download
- Media
- 24-ubc_2011_fall_kurella_venu.pdf [ 1016.52kB ]
- Metadata
- JSON: 24-1.0072070.json
- JSON-LD: 24-1.0072070-ld.json
- RDF/XML (Pretty): 24-1.0072070-rdf.xml
- RDF/JSON: 24-1.0072070-rdf.json
- Turtle: 24-1.0072070-turtle.txt
- N-Triples: 24-1.0072070-rdf-ntriples.txt
- Original Record: 24-1.0072070-source.json
- Full Text
- 24-1.0072070-fulltext.txt
- Citation
- 24-1.0072070.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
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>
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-0072070/manifest