T H E STABILITY A N D D Y N A M I C S OF S P I K E - T Y P E SOLUTIONS T O T H E G I E R E R - M E I N H A R D T M O D E L by DAVID IRON • B.A.Sc. The University of Toronto, 1988. .. -.M.Sc. The University of British Columbia, 1997 f . A THESIS S U B M I T T E D IN P A R T I A L F U L F I L L M E N T OF ".: T H E R E Q U I R E M E N T S F O R T H E D E G R E E O F D O C T O R O F P H I L O S O P H Y in T H E F A C U L T Y OF G R A D U A T E STUDIES Department of Mathematics Institute of Applied Mathematics We accept this thesis as conforming to the required standard T H E U N I V E R S I T Y OF BRITISH C O L U M B I A ; Apr i l 2001 . © David Iron, 2001 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for refer-ence and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Mathematics The University of British Columbia Vancouver, Canada Abstract A well-known system of partial differential equations, known as the Gierer-Meinhardt system, has been used to model cellular differentiation and morphogenesis. The system is of reaction-diffusion type and involves the determination of an activator and an inhibitor concentration field. It is believed that long-lived isolated spike solutions for the activator model the localized concentration profile that is responsible for cellular differentiation. In a biological context, the Gierer-Meinhardt system has been used to model such events as head determination in the hydra and heart formation in axolotl. This thesis involves a careful numerical and asymptotic analysis of the Gierer-Meinhardt system in one dimension and a limited analysis of this system in a multi-dimensional setting. We begin by studying a reduced model, referred to as the shadow system, which results from simplifying the Gierer-Meinhardt model in the limit of inhibitor diffusivity tending to infinity. This reduced model is studied in both one and in several spatial dimensions. In §2 we study the stability and dynamics of interior spike profiles for this reduced model. We find that any n-spike profile, with n > 1, is unstable on a fast time scale. Profiles with a single interior spike are also unstable but on an exponentially slow time scale. In this case the spike tends towards the closest point on the boundary. In §3 we examine the behaviour of a spike profile in which the spike is confined to the boundary. This scenario is studied in the case of a two and a three dimensional domain. It is found that the spike moves in the direction of increasing boundary curvature and increasing boundary mean curvature in two and three dimensions, respectively. Stable spike equilbria correspond to local maxima of these curvatures. We then study the case of a spike confined to a flat portion of the boundary in two dimensions. In this case it is found that the spike moves on an exponentially slow time scale. The remainder of this thesis examines the full Gierer-Meinhardt system in a one-dimensional spatial domain. In §4 we study the stability properties of n-spike equilibrium solutions to the full system. A necessary and sufficient condition is found for the linear stability of an n-spike solution. In §5 we study the dynamics of spike profiles. We derive a system of ordinary differ-ential equations which govern the motions of the spikes in one spatial dimension. Numerical computations of this asymptotic system are compared with numerical computations of the full system. In §6 we study the effects of precursor gradients. The mathematical result of these spatial inhomogeneities in the chemical reaction is that some of the coefficients in the equations are no longer constant in space. We study the effects of spatially varying activator and inhibitor decay rates as well as a spatially inhomogeneous activator diffusivity. It is found that these spatial inhomogeneities can affect both the dynamics and equilibrium position of the spikes. ii Table of Contents Abstract ii Table of Contents iii List of Tables v List of Figures vi Chapter 1. Introduction 1 1.1 The Turing Instability 2 1.2 Examples of Reaction-Diffusion Models 3 1.2.1 Long Range Inhibition 4 1.2.2 Depletion of Substrate 4 1.2.3 Microwave Heating of Ceramic Fibers 5 1.2.4 Combustion on a Slowly Diffusing FuelField - The Grey-Scot Model 5 1.3 Scaling of Gierer-Meinhardt Equations for Spike Solutions 6 1.3.1 The Simplified Gierer-Meinhardt Equations - The Shadow System 7 1.4 A Brief History of the Analysis of the Gierer-Meinhardt Equations 8 Chapter 2. The Shadow System 13 2.1 A Spike in a One-Dimensional Domain 14 2.1.1 The Nonlocal Eigenvalue Problem 16 2.1.2 An Exponentially Small Eigenvalue ; 21 2.1.3 The Slow Motion of the Spike 27 2.2 An n-Spike Solution 31 2.3 A Spike in a Multi-Dimensional Domain 34 2.3.1 The Nonlocal Eigenvalue Problem 36 2.3.2 An Exponentially Small Eigenvalue 40 2.3.3 The Slow Motion of the Spike : 46 Chapter 3. Spike Motion on the Boundary 50 3.1 Spike Motion on the Boundary in Two Dimensions 52 3.1.1 A Few Explicit Examples 56 3.2 Spike Motion on the Boundary in Three Dimensions 59 3.3 Qualitative Properties of the Associated Eigenvalue Problem 65 3.4 A Spike on a Flat Boundary in Two Dimensions 68 3.4.1 The Translation Eigenvalue 70 3.4.2 The Slow Spike Motion 72 Chapter 4. Stability of n-spike equilibrium solutions to (1.19) ' 76 4.1 An Asymptotic Analysis of the Equilibrium Solution 78 4.2 Analysis of the Large Eigenvalues 82 4.2.1 Analysis for s - 0 : 82 iii 4.2.2 Analysis for s > 0 88 4.3 Analysis of the Small Eigenvalues 91 4.3.1 Deriving the Matrix Eigenvalue Problem 91 4.3.2 Analyzing the Matrix Eigenvalue Problem 95 4.4 The Dynamics of a One-Spike Solution 103 4.4.1 The Differential Equation for the Spike Location 104 4.4.2 The Stability of a One-Spike Solution for D -+ oo 108 Chapter 5. Spike Dynamics 111 5.1 The Dynamics of Quasi-Equilibrium Solutions 112 5.2 Symmetric and Asymmetric Equilibria 118 5.3 Stability of the Profile: The Large Eigenvalues 122 5.4 Comparison of Asymptotic and Numerical Results 125 5.4.1 Two Spikes n = 2 127 5.4.2 Three Spikes n = 3 135 5.4.3 Conclusions 136 Chapter 6. Spike Pinning for the Gierer-Meinhardt Model 140 6.1 One-Spike Dynamics: The Perturbed Shadow Problem 142 6.1.1 Exponentially Small Eigenvalue 143 6.1.2 The Metastable Spike Motion 145 6.1.3 An Example of the Theory 147 6.2 One-Spike Dynamics: The Perturbed Gierer-Meinhardt Model 148 6.2.1 A Spatially Varying Inhibitor Decay Rate When D = 0(1) 148 6.2.2 Case 1: p(x) > 0 depends on x when D large 151 6.2.3 Case 2: u(x) > 0 depends on x when D is small 153 6.2.4 A variable activator decay rate 154 Chapter 7. Conclusions 157 7.1 Overview of Results 157 7.2 Possible Extensions 159 Appendix A. Proof of Theorem 2.1 160 Appendix B. The Laplacian in the Boundary Layer Coordinate System 168 B. l An Asymptotic Estimation of an Inner Product 170 Appendix C. Calculation of E and V 171 C. l Calculation of Matrix Eigenvalues of B 173 C.2 Calculation of Bg and Vg 174 C.3 Calculation of Matrix Eigenvalues of Bg 175 Bibliography 176 iv List of Tables 2.1 5 and Ao for the case (p, q, r, s) = (2,1,2,0) 21 2.2 Comparison of asymptotic and fully numerical estimates of Ai for various values of e with the parameter set (p, q,r,s) — (2,1,2,0) 25 2.3 Height of spike 1 centered at XQ = —0.5 and of spike 2 centered at x\ = 0.5 35 2.4 Logarithmic Interpolation of AQ and 0o(-5) 36 2.5 5 and An in R2 and M3 for the case of (p, <?, r, s) = (2,1,2,0) 41 3.1 Numerical results for the principal eigenvalue ni of the local problem on a flat boundary (3.41) 66 4.1 A comparison of the asymptotic and numerical results for xo{t) corresponding to the parameter values shown in the caption of Fig. 4.2 107 5.1 The numerical and asymptotic results for x\ and x2 versus r for Example la 129 5.2 The numerical and asymptotic results for x\ and x2 versus r for Example lb 130 5.3 The numerical and asymptotic results for x\ and x2 versus r for Example 2a 133 5.4 The numerical and asymptotic results for x\ and x2 versus r for Example 2b 135 5.5 The numerical and asymptotic results for x\ and x2 versus r for Example 4 136 5.6 The numerical and asymptotic results for xy, x2 and £ 3 versus r for Example 5 138 5.7 The numerical and asymptotic results for x\, x2 and £ 3 versus r for Example 6 139 List of Figures 2.1 Numerical solution for uc(y) when p = 2,3,4 15 2.2 Ao and A2 versus 5 for the parameter set (p, q, r, s) = (2,1,2,0) 20 2.3 Ao and A2 versus 5 for the parameter set (p, q, r, s) = (3,2,2,0) 22 2.4 Ai versus e for the parameter set (p, q, r, s) = (2,1,2,0). 26 2.5 xo versus t for e = .05, xo(0) = —0.4 and the parameter set (p, q, r, s) = (2,1,2,0). .. 28 2.6 Xo versus t for e = .06, xo(0) = —0.4 and the parameter set (p, q, r, s) = (2,1,2,0) .. 29 2.7 xo versus t for e = .07, xo(0) = —0.4 and the parameter set (p, q, r, s) = (2,1,2,0). .. 30 2.8 Numerical solution for uc(p) when N — 1,2,3 and p = 2 37 2.9 Ao(#) and XN+I(S) versus 8 in E 2 for the parameter set (p,q,r,s) = (2,1,2,0) 42 2.10 Xo{S) and \N+I(5) versus S in R3 for the parameter set (p,q,r,s) = (2,1,2,0) 43 3.1 A plot of u versus x at different times showing a spike merging with the boundary . 51 3.2 For Example 1 we plot the motion of the center of the spike on the boundary 58 3.3 For Example 1 we plot the solution #0 versus r 59 3.4 For Example 2 we plot the motion of the center of the spike on the boundary 60 3.5 For Example 2 we plot the solution 60 versus r 61 3.6 Plot of a two-dimensional domain with a flat boundary 69 3.7 Plot of part of a domain boundary, <9fi, upon which the center of the spike is at an unstable steady state. Ki > 0, KR > 0 for this domain 75 3.8 Plot of part of a domain boundary, 3f2, upon which the center of the spike is at a stable steady state. K~L < 0, KR < 0 for this domain 75 3.9 Plot of part of a domain boundary, dti, upon which the center of the spike moves towards the right. KL < 0 and KR > 0 for this domain 75 4.1 Plot of the activator and inhibitor concentration for a three-spike asymptotic symmetric equilibrium solution 81 4.2 Plot of the trajectory x0{t) of the center of the spike for a one-spike solution 106 4.3 Plot of the initial condition for a one-spike solution 107 4.4 Plot of a one-spike solution at two different time 109 5.1 Plot of the activator and inhibitor concentration 119 5.2 Plot of |a|i defined in (5.35) versus D for solutions with three or fewer spikes 121 5.3 Xj and em versus r for Example 1 • 128 5.4 Numerical results for a versus r for Example 1(a) and (b) 131 5.5 Xj and em versus r for Example 2 132 5.6 Numerical results for a versus x for Examples 2(a) and (b) 132 5.7 Numerical results for a versus x for Example 3 134 5.8 Xj and em versus r for Example 4 134 5.9 Xj and em versus r for Example 5 137 5.10 Xj and em versus r for Example 6 137 vi Chapter 1 Introduction Developmental biology is the study of the mechanisms by which a single cell develops into a complex organism with many different cell types. This process is extremely complex and involves both mechanical and biochemical processes. We will be focusing on the biochemical processes related to organ formation, or organogenesis, here. In this process, the initial cell divides into a large number of identical cells. Then a variety of mechanical and biochemical events occur. One such event is the formation of locally elevated levels of a substance, which we refer to as an activator, which cause the cells in a neighborhood of the elevated levels to differentiate from the surrounding cells. Thus a particular organ, such as a heart, is formed. In this thesis, we will analyze a mathematical model, known as the Gierer-Meinhardt model, of a proposed mechanism that may be responsible for the localization of the activator. This model uses the combination of a catalytic reaction with diffusion to produce an equilibrium in which the concentrations form a spatial pattern. The possible existence of such systems was first proposed by Turing in 1952 (see [46]). An outline of this chapter is as follows. In §1.1 we present an overview of the initial work done by Turing. We will then present some explicit models of processes that are known to yield equilibrium solutions with spatial patterns. Specifically, in §1.2.1, §1.2.2, §1.2.3 and §1.2.4 we present models using long range inhibition, depletion of substrates, a model of microwave heating and a model of combustion occurring on a slowly diffusing fuel source, respectively. A full discussion of the first two of these models, as well as some striking numerical simulations, may be found in [32] and [33]. The remainder of this thesis will then focus on the long range 1 inhibition model, which is commonly referred to as the Gierer-Meinhardt equations (see [12]). In §1.3 we begin the mathematical analysis of the Gierer-Meinhart equations, by finding a suitable scaling as well as deriving a simplified model commonly referred to as the shadow system [35]. In §1.4 we overview some of the history of the mathematical analysis of the Gierer-Meinhartd equations. Finally, we discuss the goals of the thesis and give an outline of the remaining chapters. 1.1 The Turing Instability Turing examined systems of the form, AT = DAAA + F{A,H), in ft, (1.1a) HT = DHAH+ G{A,H), in ft, (1.1b) with Neumann boundary conditions on 9ft. Here A represents the activator concentration, H represents the inhibitor concentration, DA is the activator diffusivity, DH is the inhibitor diffusivity and F and G are nonlinear reaction terms. For simplicity we will consider ft to be the one dimensional domain [0,1], To determine if stable spatial patterns are possible, we examine the stability of the spatially homogeneous solution. We let AE, HE be such a solution (so that F(AE,HE) = 0 and G(AE,HE) = 0), and we consider a sinusoidal perturbation from this state. We let, A = A E + a cos{uix)ext, (1.2a) H = HE + hcos{uJix)ext. (1.2b) Here WJ = ni in order to satisfy the boundary conditions. We substitute (1.2) into (1.1) and expand F and G in a Taylor series about AE and HE- This results in the eigenvalue problem, Aa = -u;fDAa + K\a + K2h, (1.3a) Xh = -ujfDHh + K3a + KAh, (1.3b) 2 where, dF dF K I = QX(AE,HE), K2 = —(AE,HE), (1.4a) dG dC K3 = —(AE,HB), K4 = —(AE,HB). (1.4b) We must assume that the spatially homogeneous solution (AE,HE) is always stable in the absence of diffusion. This assumption results in several restrictions on the values of Ki. If this previous assumption is true, and the following condition is met, (KiDff + K4DA)2 - \DADH(KlK4 - K2KZ) > 0, (1.5) then the eigenvalue problem (1.3) will have eigenvalues with Re(A^ ) > 0 when, 4((^! - ufDA)(K4 - u?DH -K2Kz) < [Kx +K4- u*(DA + DH)f . (1.6) Thus, for a certain range of values, sinusoidal perturbations may grow. The end result is that only certain modes become unstable. Through numerical simulations, such as those in [12] and [16], these instabilities have been shown to yield spike-type solutions when DA is small. The non-linear effects will then cause the spikes to stabilize. For an in-depth analysis of the parameter regime for spike growth see [28]. Numerical methods have now made it possible to simulate such systems thus verifying the existence of spike-type solutions. .1.2 Examples of Reaction-Diffusion Models We now present some mathematical models of reaction-diffusion systems that are known to produce solutions with spatial patterns. Each system has different properties making them suitable for modeling different phenomena. For example, the depletion of substrate system seems to produce periodic patterns and thus may be of use in the modeling of structures with a periodic nature such as the spinal cord. The long range inhibition model tends to produce isolated structures and thus may be more appropriate for the modeling of isolated structures such as the heart. The microwave heating equation models the formation of hot-spots that can occur during sintering. The Grey-Scott model of combustion on a slowly diffusing fuel field produces pulses of locally elevated temperature that can split into traveling pulses which may again divide. We now present some details of the models. 1.2.1 Long Range Inhibition In this two component reaction-diffusion system, the activator is a slowly diffusing substance which promotes its own formation with autocatalysis. The inhibitor is a rapidly diffusing substance which uses the activator as a catalyst in its formation and itself as an inhibitor of both its own formation and that of the activator. The formation of an isolated peak happens as follows. If the system is at a spatially homogeneous equilibrium, small perturbations in the concentration of the activator will grow due to the autocatalysis. These elevated levels of activator will promote a localized increase in the concentration of inhibitor. The elevated levels of inhibitor then diffuse rapidly preventing the formation of spikes in the activator concentration elsewhere. We now present a mathematical model of this system in a one-dimensional domain and in dimensionless form. This model, known as the Gierer-Meinhardt model, is AP At = e2Axx-A + — , -1<X<1, t>0, (1.7a) rHt = DHXX - pH + , - K x < l , t>0, (1.7b) Ax{±l,t) = Ha{±l,t) = 0, (1.7c) where the exponents (p, q, r, s) are assumed to satisfy the relations, p > l , q>0, r>0, s>0, 0 < ^—^ < . (1.8) q s + 1 Here D and e2 are the diffusivities of the inhibitor and activator, respectively. This model is a scaled and slightly simplified version of that presented in [12]. Simulations of this system show that it supports extremely robust isolated spike solutions. This is the model that is studied in this thesis. 1.2.2 Depletion of Substrate In this model, the activator is autocatalytic but instead of an inhibitor, a rapidly diffusing substrate is required and used by the catalytic reaction. Thus, a localized slight increase of the concentration of the activator will grow due to autocatalysis. The autocatalytic reaction will use up the substrate required for the growth of the activator and thus regulate the size of the 4 spike. The details of the mathematical model of this reactions are as follows. At = DaAxx - \xA + cA2S, (1.9a) St = DSSXX -vS + co- cA2S. (1.9b) This model is commonly referred to as the Brusselator (see [40]). As previously mentioned, this model has a tendency to produce periodic patterns. 1.2.3 Microwave Heating of Ceramic Fibers We now give a partial differential equation modeling the microwave heating of thin ceramic cylinders in single-mode highly resonant cavity. An analysis of the equations in the limit of small Biot number results in the following non-local reaction diffusion equation (see [4]), Ut = DUxx-2(U + B[(U + l)*-l])+P- ^ } , 0 < x < l , (1.10a) Ux{0) =US{1) =0, (1.10b) where U is the scaled temperature of the ceramic cylinder, D is the thermal diffusivity, x and 3 are dimensionless physical parameters, P is dimensionless power parameter and f(U) is the effective electrical conductivity. In [4] it shown that (1.10) supports highly localized hot-spot solutions. The forms of (1.10) and an asymptotic reduction of (1.17) known as the shadow system, given below in (1.18), are very similar. Thus, the analysis for the shadow system given in §2 should be useful for analyzing (1.10). 1.2.4 Combustion on a Slowly Diffusing Fuel Field - The Grey-Scot Mode l This model is used for the simulation of brush fires or small forest fires. One of the most interesting features of this model is that it has spike solutions that may split into two spikes which then travel in opposite directions. This splitting process may then continue. In one spatial dimension, the model in non-dimensional form is (see [9]), Ut = Uxx - UV2 + 5a(l - U), -Kx<l (1.11a) Vt = 52VXX + UV2 - 8^bV - K x < l , (1.11b) UX(±1) = VX(±1)=0. ' (1.11c) 5 Here U is a scaled temperature field and V is a fuel concentration field. The model is studied in [9] in the limit of small 5, for various ranges of the positive parameters a, b and fi. 1.3 Scaling of Gierer-Meinhardt Equations for Spike Solutions Before we begin the scaling of (1.7), we make the simplification of setting r = 0, as finite values of T can lead to complicated oscillatory behaviour which is beyond the scope of this thesis (see [34]). From numerical computations, it is evident that for sufficiently small values of r, this simplification will cause no difficulties. Thus, for simplicity we r = 0 in (1.7). The amplitude of a spike solution to (1.7) will tend to infinity as e —>• 0. Therefore, we introduce new variables for which the spike solution has an 0(1) amplitude as e -> 0. To this end, we introduce a and h by A = e~"aa, H = e~Vhh, (1.12) where va and are to be found. To balance the terms in (1.7a), we require, -va = -vap + qvh. (1-13) We will construct a solution in which the spike has its support in an 0(e) region near some point in Ct. Therefore, to obtain an additional equation relating ua and Vh we consider an average balancing of (1.7b). Specifically, we integrate (1.7b) over the domain to get, /l rl AT Hdx + J j^dx = 0. (1.14) Since A will be localized to an O(e) region about the spike center XQ, we scale x in the last term by y = <Tl{x — xo). Balancing the terms in this equation we get -uh = -var + uhs + 1. (1.15) The solution of (1.13) and (1.15) yields, 9 = JP- M 1 6 y V a (l-p){l + s)+rq, U h {l-p)(l + s)+rq-6 In terms of these new variables (1.7) becomes 2 aP j o>t = ^ axx - a+—, - 1 < X < 1 , £ > 0 , (1.17a) T 0 = Dhxx-fj,h + e~1^-, -Kx<l, t>0, (1.17b) h3 ax{±l,t) = hx(±l,t) =0. (1.17c) 1.3.1 The Simplified Gierer-Meinhardt Equations - The Shadow System In §1.2.1 it was stated that the inhibitor must be a rapidly diffusing substance. Thus, it seems logical to study this system in the limit D —> oo. In this hmit, h is a constant in space which may be determined from a solvability condition. In this limit, (1.17) reduces to the following nonlocal ordinary differential equation, at = eaxx - a + — , -1 < x < 1, t>0, (1.18a) h = { ^ £ / d x y + i ' ^ ax{±l,t) = 0. (1.18c) All of the analysis above can also be carried out in M.N. We omit the details as the analysis is almost identical to that given above except for some minor changes in the scaling. As the analysis of the full system in RN is beyond the scope of this thesis, we only present the N-dimensional version of the shadow system. It is given by, bfl" o aP at = e2Aa-a + — in Q, t>0, ' (1.19a) dna = 0 on dQ, (1.19c) where ft is a bounded subset of RN and dn refers to the normal derivative. 7 1.4 A Brief History of the Analysis of the Gierer-Meinhardt Equations There are very few results for (1.17). Those that we have found will be discussed in the main body of the thesis at the appropriate times. Here we will restrict ourselves to results found previously for (1.19). The first results are for equilibrium solutions to (1.19). If we set at = 0 and a = h^u in (1.19), then (1.19) reduces to, e2Au -u + up = 0, in ft, (1.20a) dnu = 0, on dti. (1.20b) A spike solution for (1.20) in the limit e —>• 0 is a solution of the form, u(x) ~ uc[e~l\x — XQ\] where Xrj is a point in ft or on <9ft to be determined, and uc(p) is the unique positive radially symmetric solution to, uc —> ae p as p—> oo, (1.21b) for some a > 0. An n-spike solution is defined as, n U ~ /\^uc[e~1\x ~ xi\] > (1.22) where uc is defined in (1.21) and the spikes are centered at points Xi, i = 1,... ,n, to be determined. In [35] the existence of a single spike solution centered on the boundary at the global maximum point of the boundaries mean curvature is proven using variational methods, in particular the Mountain Pass Lemma. In axially symmetric domains, multiple peak solutions, with the peaks all resting on the boundary at local critical points of the boundaries mean curvature are constructed in [36]. In [13] multiple peak boundary spike solutions are found in a general domain with a smooth boundary. In [14] solutions with multiple interior peaks are constructed where the locations of the peaks coincide, in the limit e ->• 0, with centers of spheres that solve a ball packing problem in the domain. In [15] solutions with spikes resting in the interior as well as on the boundary are found. 8 There axe some key results in [48], which we will examine in more detail, as many of the ideas in this paper are important to an understanding of the analysis in this thesis. The result which concerns us states that all of the solutions found for equation (1.20) are dynamically unstable when the nonlocal effects from (1.18) are ignored. We briefly explain how this result is obtained. If we consider U in (1.20) to be all of RN, then uc ( ^ " f 0 ' ) from (1.21) will satisfy (1.20) for any choice of XQ. If we consider now a finite domain, uc will still satisfy (1.20a) and will fail to satisfy the boundary conditions by exponentially small amounts for any choice of XQ in the interior of 51. Next, we linearize about this approximate solution and obtain the following eigenvalue problem, Le<j) = e2A$ + (-l+pup-1)<t> = X(j>, (1.23a) ^ = 0 on. 3(1. (1.23b) We note that the function fa = for i = .1 , . . . ,N satisfies Lefa = 0 and fails to satisfy the boundary conditions by only exponentially small amounts. This suggests that in the limit e -» 0, (1.23a) has N exponentially small eigenvalues. We note that the eigenfunctions associated with these eigenvalues each have exactly one nodal line. Thus, these are not the principal eigenfunctions. Therefore, (1.23a) must have a proper positive eigenvalue. Since, by a rescaling of the spatial variables, we may eliminate the e from (1.23a), this eigenvalue must be 0(1). Thus, the linearization about a spike solution has a positive 0(1) eigenvalue and therefore, these spike solutions are not stable. These large and small eigenvalues will reappear many times in the course of this thesis. One of the main tasks in this thesis is to show what conditions are necessary to move the positive 0(1) eigenvalue into the left half of the complex plane. We now present the main goals of this thesis and an overview of the remainder of the chapters. Turing did his initial work in 1952 using the techniques of linear analysis. The only questions that this form of analysis has been able to answer are what are some of the necessary conditions for spikes to form from spatially homogeneous initial data. The conditions arriving from this analysis are far from sufficient. Once the initial Turing instability is triggered it is assumed that the nonlinear effects will stabilize the spike profile. This linear analysis can say nothing about the fully formed spike profile, which is far from the linear regime. Until recently, the only 9 way to verify that the fully formed spike profiles are stable is to use numerical analysis. One of the main goals of this thesis is to use modern analytical techniques to find explicit criteria, both necessary and sufficient, determining the stability of a particular spike profile. There are many other questions that this thesis will address. We will consider the behavior of spikes after they are fully formed. We will also address such questions as how do the spikes move and how do they interact dynamically. We will also examine the differences in the behaviour of the full system (1.17) with the shadow system (1.18) and determine under what circumstances it is appropriate to use the reduced model (1.18). The remainder of the thesis proceeds as follows. In §2 we consider equilibrium and quasi-equilibrium solutions to (1.18) and (1.19). In §2.1 we construct a canonical spike solution to (1.18). A linearization about this solution leads to a non-local eigenvalue problem. The principal eigenvalue of this spectrum is shown to be exponentially small. A projection method is then used to derive a solvability condition resulting from the exponentially small eigenvalue. This leads to an ordinary differential equation governing the exponentially slow motion of one spike. These results are then favorably compared to numerical simulations of (1.18). We now present the main result from this chapter. For e -» 0, metastable spike solution for (1.18), is represented by a(x,t) = CLE{X; xo(t)), where a# is defined below in (2.1) and xo(t) satisfies, ^ ^Ll \p-2(l-xo)/e _ p-2(l+x0)/e dt p (1.24) Here a and $ are positive constants defined below in (2.5) and (2.20), respectively. In §2.2 we demonstrate than an n-spike solution, with n > 1, to (1.18) is always unstable with an 0(1) positive eigenvalue. In §2.3 we repeat the analysis of §2.1 in an iV-dimensional setting. The differential equation governing the motion of the spike shows that the spike will drift exponentially slowly towards the closest point on the boundary. In §3 we examine the dynamics of the spike once it has merged with the boundary. This is carried out in two and three dimensions in §3.1 and §3.2 respectively. In both cases it is found that the spike moves in the direction of increasing curvature (mean curvature in the case of three dimensions) until reaching a stable equilibrium at a local maximum of the (mean) curvature of 10 the boundary. We now provide the main results of this chapter. For e -> 0, the motion of a spike confined to the smooth boundary of a two-dimensional domain follows the trajectory * o ( * ) ~ | ^ 3 « W , (1.25) where 7 = q/(p — 1) and b > 0 defined below in (3.20c). Here s is the arclength along the boundary and «(s) > 0 is the curvature of the boundary. This differential equation has stable equilibrium when we are at a local maximum of K(S). Correspondingly in three dimensions, for e -> 0 the motion of a spike confined to a smooth boundary is described by £'(t)~llK?VH{ZJ. (1.26) Here £ = (£1, £2)7 7 = l/(p ~ l ) i b > 0 is defined in (3.37b) and H is the mean curvature of <9ft, with H > 0 for a sphere. Stable equilibria are at local maxima of H. In §4 we examine the stability of n-spike solutions to (1.17). We find explicit conditions, both necessary and sufficient, to determine the stability of an n-spike solution. This is accomplished by linearizing about an n-spike solution and finding criteria that guarantee that both the large and small eigenvalues discussed previously lie in the left half of the complex plane. The main result of this chapter is a stability result which we now provide. An n-spike equilibrium solution to (1.17) is stable when £><£>* where, as e -> 0, -1 (1.27) An n-spike solution is unstable when D > D*. In §4.4 these results are used to construct an ordinary differential equation governing the motion of a single spike profile. These results are compared with full numerical simulations of (1.17). These results predict D* — 00. A more refined analysis shows that a one-spike solution is stable only when D < 0(ecle). Thus the shadow system (1.18) and the the full system (1.17) give similar stability conclusions for a one spike solution only when D > 0(ec/e). In §5 we examine the dynamic properties to n-spike profiles under (1.17). A solvability con-ditions resulting from the small eigenvalues results in a differential algebraic system governing 11 D* ~ r1 2 ' qr LP + the motion of the spikes. However, the validity of this system will depend on the location of the large eigenvalues in the complex plane. Thus, a criterion is developed to ensure the validity of the derived system at any given time. We now provide some of the details of this result. For e <§C 1, the quasi-equilibrium solution for a and h is given by n a(x,t)~ac = 'Y^h'jUc[e~1(x-Xj)]J (1.28a) i=i h{x,t)^hc = brYJhJ~SG{x;xj), • (1.28b) 3=1 where hi = /ij(r) and Xi = Xi(r) satisfy the differential-algebraic system for i = 1,.., n n hi = brY,h]T~sG(xi;xj), (1.29a) 3=1 ( \ dxi 2qbT , _x dr p — 1 (1.29b) h?-'{Gs)i + YfhT~sGx(xi;xj) \ 3% J Here uc satisfies (5.5), bm is defined in (5.11a), r = e2t, and (Gx)% = [Gx{x%+\Zi) + Gx{xi-',Xi)] /2. Here G is the Green's functions satisfying, DGXX - fiG = -5{x - xk), - 1 < X < 1 , (1.30a) Gx{±l;xk) =0. (1.30b) In §6 we examine the effects of precursor gradients on the equilibrium position of a single spike. A precursor gradient is a pre-existing spatially inhomogeneity which is imposed on the model. It could come into existence from the inherent polarity of the developing cells or due to localized sources of a precursor which diffuse through the group of cells. The mathematical consequence of a precursor gradient is that the coefficients in the model will now depend on space. This can have a variety of effects depending on the nature of the inhomogeneity. Several different scenarios will be explored. 12 Chapter 2 The Shadow System We begin by considering the shadow system (1.18) which is a simplified model of the Gierer-Meinhardt system. In [35] it was shown that for sufficiently large values of D the solutions to the shadow system and that of the full system will coincide. However, we will show in §4.4.2 that D will have to be asymptotically exponentially large as e -> 0 for this to be the case and thus outside the range of any physical model. The shadow system still warrants investigation as techniques used in its investigation will be useful in the investigation of the full system. In addition, the shadow system provides an interesting example of a non-local differential equation and a novel application of the projection method. The outline of this chapter is as follows. In §2.1 we consider the one-spike quasi-equilibrium solution to the one-dimensional problem (1.18). We examine the stability and dynamics of this solution by analyzing the spectrum of the linear operator resulting from a linearization of (1.18) about this non-constant solution. This eigenvalue problem is a non-local Sturm-Liouville problem of the type considered in [11]. A combination of analytical and numerical techniques will be used to demonstrate that the principle eigenvalue of this operator is exponentially small. The non-local term in the Sturm-Liouville operator is essential for this conclusion. The exponentially small eigenvalue will be estimated asymptotically. A differential equation characterizing the motion of the center of the spike will be derived in the limit e -»• 0 by using a limiting solvability condition, which requires that the solution to the quasi-steady linearized problem has no component in the eigenspace associated with the exponentially small eigenvalue. This procedure is known as the projection method and has been used in other contexts (see 13 [48], [49]). The resulting ODE for the motion of the center of the spike shows that the spike drifts exponentially slowly towards the point on the boundary closest to the initial location of the spike. This metastable behavior is verified by calculating full numerical solutions to (1.17) in §2.2. Solutions with n-spikes, where n > 1, will be shown to be unstable. In §2.3 we give a similar analysis of metastable spike-layer motion for the multi-dimensional problem (1.19). 2.1 A Spike in a One-Dimensional Domain We first construct a one-spike quasi-equilibrium solution a# for (1.18) in the form a = aE{x;x0) = h?uc[e~l{x - x0)], j = q/(p-l). (2.1) Here xo, with. |xo| < 1, is the center of the spike. The function uc(y), called the canonical spike solution, is the unique positive solution satisfying, u'c - uc + u\ = 0, 0 < y < oo, (2.2a) u'c(0)=0; uc(y) ~ ae~y, as y->-oo, (2.2b) with a > 0. The existence and uniqueness of the solution to (2.2a) is proven in [45]. In terms of this solution, h = hE, where p - i / 1 Cl \ 0>+l)(p-l)-9r h'~ • <2'3' Since uc is localized near xn, we estimate as e —> 0, that ( ^ ^ ^ , J . j T W * . (2.4) To determine numerical values for certain asymptotic quantities below we must compute uc(y), 3, and a numerically. The constant a is obtained by integrating (2.2) log(2±i) f{*py> log(a)- V ; ' ' p - i Jo-- i 2 _ _2_„P+1 V dt). (2.5) To compute uc numerically, we use the asymptotic boundary condition uc + uc = 0 at y = yi, where yr, is a large positive constant. To compute solutions for various values of p, we use a 14 continuation procedure starting from the special analytical solution, uc(y) = ^sech2 ( | ) , ' (2.6) which holds when p — 2. The boundary value problem solver COLNEW (see [3]) is then used to solve the resulting boundary value problem. In Fig. 2.1, we plot the numerically computed uc(y) when p = 2,3,4. X 0 — 1 I i 1 1 1 1 1 1 1 4 A -1 2 -1 - \ \ \ .P = 3 -0 8 - • -0 6 • \ \ \ -0 4 - '••\ \ -0 2 0 / P = 2 -^ T^t—r.i 1 1 1 1 | 1 0 2 4 6 8 10 12 14 16 18 20 y Figure 2.1: Numerical solution for uc(y) when p = 2,3,4. We note that, for any XQ with \XQ\ < 1, the solution a£(x;xn) -will satisfy the steady-state problem corresponding to (1.18a), (1.18b), but will fail to satisfy the boundary conditions in (1.18c) by only exponentially small terms as e -> 0. Thus, we expect that the spectrum of the eigenvalue problem associated with the linearization about as will contain an exponentially small eigenvalue. 15 2.1.1 The Nonlocal Eigenvalue Problem Let XQ G (—1,1) be fixed and linearize (1.18) around aE, hE. We obtain the eigenvalue problem for the linearization by introducing (j) and n by a(x,t) = aE{x;x0) + e A V (2-7a) h(t) = hE + extn, . (2.7b) where 0 <?; 1 and n <C 1. Substituting (2.7) into (1.18) we obtain the following non-local eigenvalue problem of Sturm-Liouville type on [—1,1]: —ip ri Le<f> = e2<f>xx + (-1 +pup-1)<l> ~ J T / ure-1<l>dx = X(/>, (2.8a) 2/3(s + 1) y_x ^«(±1) = 0. (2.8b) The non-local integral term in (2.8) will drastically change the nature of the eigenvalue problem. In (2.8), uc = uc [e_1(x — XQ)]. Therefore, we will only seek eigenfunctions that are localized near x = XQ. These eigenfunctions are of the form 4>(y) = 4>(xQ + ey), y = eTl{x-xQ). (2.9) Therefore, we can replace the finite interval by an infinite interval in the integral in (2.8) and impose a decay condition for 4>(y) as y —• ±oo. This gives us the non-local eigenvalue problem for the infinite domain — oo < y < oo: 4> -> 0 as y -> ±oo. (2.10b) Now we examine the spectrum of the operator in (2.10). To demonstrate that this operator has no eigenvalues with positive real part, we apply a theorem from [52]. 16 Theorem 2.1 (Wei [52]) Consider the eigenvalue problem for 70 > 0 <fr - * + j n ? e ^ - 7 0 ( p - l K ^ JZo[Uc{y)]r d y )=™' - O C < y < C C , (2.11a) $->0 as | y | -»oo , (2.11b). corresponding to eigenpairs for which A 7^ 0. iTere uc satisfies (2.2). Let Ao 7^ 0 6e the eigenvalue of (2.11) with the largest real part. Then, i/70 < 1, we conclude that Re(\0)>0. (2.12) Alternatively, if 70 > 1 and if either of the following two conditions hold (i)r = 2,l<p<5, or (ii) r = p +1, p > 1, (2.13a) then Re{X0)<0. (2:13b) Thus for all the parameter sets satisfying (2.13a) we have that the operator Le has no eigenvalues with positive real part. The non-constant coefficients of both the operators Le and Le are both localized to a small region about xo, so we will only consider eigenfunctions which are also localized to an 0(e) region about xo- Such eigenfunctions of the operator L€ satisfy (2.8a), but fail to satisfy (2.8b) by exponentially small terms. Thus, we expect the eigenvalues of Le to differ from those of Le by exponentially small terms. The operator L€ has one zero eigenvalue with eigenfunction u'c. We will carefully examine below the effect of a finite domain on this eigenvalue for the operator Lt. The important point is that Theorem 2.1 shows that the operator Le has no 0(1) eigenvalues with positive real part. The proof of this theorem is given in Appendix A. It is instructive to see how the presence of non-local term in (2.8) effects the operators spectrum. We may use a numerical solution to see the precise effect the non-local term has on the spectrum 17 of L€. To treat the non-local eigenvalue problem numerically, we split the operator Le into two parts, Acj> = e2cj>xx + {-l+p^-l)fa Bt=2(3(s + 1) f ^ * d x • (2-14) We define a new operator Ls by Ls4> = A(f> - 8B(j). where 8, with 0 < 8 < 1, is a continuation parameter. When S = 0 we have a simple Sturm-Liouville problem. At 5 = 1 we have our full non-local eigenvalue problem (2.8). We define Ls, A and B in a similar fashion, but on the extended domain — oo < y < oo with the appropriate decay boundary conditions at ±oo. The zero eigenvalue of the operator L€ corresponds to the eigenfunction u'c, which decays exponentially as |y| —> oo. To see this, we differentiate (2.2a) with respect to y, to show that Au'c — 0. This is translation invariance. In addition, due to the symmetry of uc(y), we also have Bu'c — 0. For the finite domain problem (2.8), the function u'c [e - 1(x — XQ)] fails to satisfy the equation and boundary conditions in (2.8) by only exponentially small terms as e —> 0. Therefore, as estimated carefully below, the presence of the finite domain will perturb the zero eigenvalue and the corresponding eigenfunction of the extended problem (2.10) by only an exponentially small amount. Thus, Le has an exponentially small eigenvalue. The function uc(y) has a unique maximum at y — 0 and thus the eigenfunction u'c{y) has exactly one zero at y = 0. This implies that uc(y) corresponds to the second eigenfunction of A. The principal eigenvalue of A is simple, positive and independent of e. The principal eigenvalue of A is exponentially close to the principal eigenvalue of A. Hence, in the absence of the non-local term, the operator L€ has an 0(1) positive eigenvalue and no metastable spike motion can occur. Since Lg has a positive eigenvalue when 5 = 0, we must consider what happens to this eigenvalue as 5 ranges from 0 to 1. If this eigenvalue remains positive then, since the eigenvalues of Ls and Ls will differ by only exponentially small amounts as e —• 0, we can conclude that the one-spike quasi-equilibrium solution is unstable. Alternatively, if this eigenvalue crosses through zero at some finite value of 5 < 1, then the principal eigenvalue of Ls when 8 = 1 (which corresponds to our eigenvalue problem (2.8)) will be exponentially small. Hence, if this occurs, the one-18 spike solution is anticipated to be metastable. Applying theorem 2.1 to the operator Le, we conclude that for any parameter set satisfying (2.13a), the eigenvalue should cross through zero The calculation of the eigenvalues of Lg will require some numerical analysis. Thus, we will work with specific parameter sets. We first consider the set (p,q,r,s) = (2,1,2,0), which is commonly used in simulations and satisfies (2.13a). For this parameter set, we begin by reviewing some exact results for the spectrum of the local eigenvalue problem A4> - 4>yy + (-1 + pupc~l)4> = X<j> - oo < y < oo, (2.15a) 0 as y->-±oo. (2.15b) This problem has three isolated eigenvalues, and a continuous spectrum in the left half plane. When p = 2, these three isolated eigenpairs are (see [26]), A0 = 5/4, = sech3(y/2), (2.16) Ai = 0, fa = tanh(y/2)sech2(y/2), (2.17) A2 = -3/4, 02 = 5sech3(y/2)-4sech(y/2). (2.18) Since these eigenfunctions, written in terms of y = e~l(x — xo), will fail to satisfy the boundary conditions in (2.8) by only exponentially small terms as e ^ 0, we expect that the eigenvalues of A will be only slightly perturbed from those of A. As we have previously noted, the zero eigenvalue of (2.15) will persist for L$ as 8 ranges from zero to one. Hence, there is an eigenvalue of (2.8) that is exponentially small as e —» 0. To numerically compute the eigenvalue branches Xo(S) and A2(o~) of Lg for which \o{8) -> 5/4 and \2(8) -» —3/4, as 8 —• 0, we use the initial guesses, provided above for 8 = 0 and a continuation procedure to compute these eigenvalues as 8 increases. The computations are done using COLNEW. The analysis of [11] showed that these eigenvalue branches are smooth functions of 8, and they cannot terminate suddenly at some value of 8. Hence, 8 is a natural homotopy parameter. In Fig. 2.2 we plot the numerically computed XQ(8) and X2(8) versus 8. As can be seen from this graph, A0 « 0 for 8 = 1/2, agreeing with Theorem 2.1. As 8 19 increases past 1/2, Ao becomes negative and then complex. At this point, COLNEW is no longer able to track the eigenvalue. As 8 increases from 0 to 1, Ao decreases and A2 increases. At a value of 8 « 0.65 the two eigenvalues collide and split into complex conjugates eigenvalues with negative real parts. To track the eigenvalues beyond 8 « 0.65 one must employ a different numerical technique. We accomplish this by discretizing the finite domain problem (2.8), which has eigenvalues exponentially close to those of Ls- This is done using a centered difference approximation applied to the second derivative and Simpson's rule applied to the integral. Thus, the operator Ls is approximated by a discrete linear operator Ls- The eigenvalues of the continuous problem may then be approximated by the eigenvalues of this matrix. Numerical calculations of the eigenvalue Ao of Ls are shown in Table 2.1. As seen in Table 2.1, the real part of Ao remains negative as 8 tends to one. Similar computations, with similar conclusions, can be performed for other values of p, q, r and s. In particular, Ao and A2 are shown in Fig. 2.3 for the parameter set (p, q, r, s) = (3 ,2 ,2 ,0 ) . 1.4 1 r 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 8 Figure 2.2: AQ and A2 versus 8 for the parameter set (p, q, r, s) = (2 ,1 ,2 ,0 ) . 20 6 Ao 0.0 1.2518 0.1 1.0073 0.2 0.76149 0.3 0.51345 0.4 0.26158 0.5 0.0052548 0.6 -0.28247 0.7 -.59237 + 0.15315* 0.8 -.71522 + 0.23035* 0.9 -.84093 + 0.23008* 1.0 -.98551 +0.14507* Table 2.1: 8 and Ao for the case (p, q,r, s) = (2,1,2,0). 2.1.2 A n Exponentially Small Eigenvalue In the previous section, we explained qualitatively why the principal eigenvalue of Le is expo-nentially small. The non-local term in (2.8) was found to be essential to this conclusion. In this section we calculate the exponentially small eigenvalue precisely. We denote the eigenpair corresponding to the exponentially small eigenvalue by Ai, 4>i. To predict the dynamics of the quasi-equilibrium solution, we must obtain a very accurate estimate of Ai. We expect that 01 ~ C\uc (e~l(x — XQ)) in the outer region away from 0(e) boundary layers near x = ±1. The behavior of </>i in these regions will be analyzed using a boundary layer analysis. The eigenfunction <f>\ has the boundary layer form fa{x) = Ci (u'c [e'^x-xo)] +<j>i [e_1(x + l)] + 0r [e^l-a:)]) . (2.19) Here <pi (77) and (pr (77) are boundary layer correction terms and Ci is a normalization constant 21 - I 1 1 1 r - i r 2.5 ^1.5 0.5 h oh--0.5 -A2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 5 ^ Figure 2.3: An and A2 versus 6 for the parameter set (p, q, r, s) = (3,2,2,0). given by (*) -1/2 where /3 /oo , •oo (2.20) In the boundary layer region near x Thus, as e -» 0, <&(") satisfies —1, [e x(x — XQ)] is exponentially small as e -> 0. fa — fa = 0, 0 < 77 < 00. ^(0 ) ~ - a e - e _ 1 ( 1 + I ° ) . (2.21a) (2.21b) Similarly, the boundary layer equation for <pr{il) I S 4>r — (f>T = 0, 0 < TJ < OO . 4>r(0) ~ ae -e-^l-xo) (2.22a) (2.22b) 22 Here a is defined in (2.5). Solving the boundary layer equations we get <f>r(T]) = -ae-e~1(l-x^e-r>. (2.23a) (2.23b) To estimate Ai we first derive Lagrange's identity for (u,L e v), where (it,u) = §_xuvdx. Using integration by parts we derive (v, L€u) — e2 (uxv - vxu) | l x + (u, L*v), (2.24) where L\v = e2vxx + (-1 + pvP-l)v - f v?cv dx. (2.25) We now apply this identity to the functions u'c[e~l(x — XQ)] and fa(x) to get Ai(u'c)<t>x) = - € 0 i < | L i + (fa,L*u'c) • (2.26) We examine each of the terms in (2.26). We begin with (u'c, 4>\). The dominant contribution to this integral arises from the region near x = XQ where fa ~ C\u'c[e~l{x — XQ)}. Therefore, the inner product can be estimated as {u'c,fa) = Ci{u'c,u'c) -deft. (2.27) Next, to estimate —efau'^\}_1, we use our asymptotic estimates of uc and fa. Since uc(z) ~ ae - ' 2 ' as z ->• ±oo we have that u'c [ e _ 1 ( ± l — XQ)] ~ ae~€~l^x°\ In addition, using the previous boundary layer results for fa we get the following estimate for fa(±l): fa{±l) ~ ^ d o c - ' " ^ 1 ^ ' . (2.28) Using these results, we get ~ 2eCia 2 (e-*~Hi+*o) + e - 2 6 - i ( i - x 0 ) ) _ ( 2 2 9 ) 23 The only term left to examine is (<f>i, L*u'c). Since u'c is a solution to the local operator, we have rqe lurc 1 f 1 2/3(s+ l ) ( p + l) (2.30) where -pe ^ l i i o ) (2.31) In a similar way, the term (<pi,L*u'c) is approximated by ( ^ , L : O ~ -2/3(a + l ) ( p + l) r g e - 1 a P + r + 1 C i 20(s + l ) ( p + l)r r g e - V + 1 C a (2.32) Since p > 1 and r > 0, upon comparing the terms in (2.32) and (2.29), it is clear that the second term on the right side of (2.26) is asymptotically negligible compared to the first term. Finally, substituting (2.27) and (2.29) into (2.26), we get the following asymptotic estimate for the exponentially small eigenvalue Ai as e —» 0: In (2.33), a and /? are defined in (2.5) and (2.20), respectively. The estimate (2.33) holds for p, q, r and s satisfying (1.8). Since Ai > 0, the quasi-equilibrium profile is unstable. However, since A i is exponentially small, the instability is extremely weak and the quasi-equilibrium profile can persist for an exponentially long time interval. To verify the estimate for A i , we also numerically estimate Ai by solving (2.8) using COLNEW. In Fig. 2.4, we compare the numerically computed values of Xi (the dots) with the asymptotic es-timate (2.33) (dashed curve) for various values of e for the parameter set (p, q, r, s) = (2,1,2,0). The asymptotic and numerical results are also shown in Table 2.2. Similar favorable compar-isons can be made for other parameter sets. (2.33) 24 e Ai(full numerics) Ai (asymptotic) .220000 .12005x10" -l .135223x10" -l .210000 .80010x10" -2 .877088x10" -2 .200000 .50829x10" -2 .544799x10" -2 .190000 .30582x10' -2 .321855x10" -2 .180000 .17289x10" -2 .179344x10" -2 .170000 .90945x10" -3 .932899x10" -3 .160000 .43964x10" -3 .447198x10' -3 .150000 .19224x10" -3 .194352x10" -3 .140000 .74490x10" -4 .74985x10-4 .130000 .24894x10" -4 .249878x10" -4 .120000 .69198x10" -5 .69333x10-5 .110000 .15224x10" -5 .152376x10" -5 .100000 .24725x10" -6 .247338x10" -6 .90000 xlO"1 .26800x10" -7 .268036x10" -7 .80000 xlO"1 .16665x10" -8 .166655x10" -8 .70000 x lO - 1 .46857x10" 10 .468562x10" 10 .60000 xlO"1 .40156x10-12 .400589x10- 12 Table 2.2: Comparison of asymptotic and fully numerical estimates of Ai for various values of e with the parameter set (p, q, r, s) = (2,1,2,0). 25 0 .009 0.001 Figure 2.4: Ai versus e for the parameter set (p, q, r, s) = (2,1,2,0). We end this section with a few remarks. Firstly, we recall that Ai and fa ~ C\uc (e-1(x — xo)) are an eigenpair of Ls when 8 = 0. To within negligible exponentially small terms this eigenpair remains an eigenpair of Ls as 8 ranges from 0 to 1. To see this, we note that the only difference between the calculations of the eigenvalue for the local problem and for the non-local problem, is that the term {L*tu'c, fa) in (2.26) would be replaced by (A fa, fa) = 0, since A is self-adjoint. In the final calculation of Ai the term (L*uc, fa) was ignored since it is asymptotically exponentially smaller than the other terms in (2.26). Secondly, we note that Ai, fa is exponentially close to an eigenpair X\, <j>\ of the adjoint operator, i * . For the same reasoning as above, 4>* would have the same interior behavior near x = XQ as fa and the same boundary layer correction terms near x = ±1. Repeating the calculation to find X\, we would arrive at the same estimate as in (2.33). 26 2.1.3 The Slow Motion of the Spike The quasi-equilibrium solution fails to satisfy the steady-state problem corresponding to (1.18) by only exponentially small terms for any value of XQ in \XQ\ < 1. Moreover, the linearization about this solution admits a principal eigenvalue that is exponentially small. Therefore, we expect that the one-spike quasi-equilibrium profile evolves on an exponentially slow time-scale. We will now derive an equation of motion for the center of the spike corresponding to the quasi-equihbrium profile. To do so we first linearize (1.18) about a(x, t) = / i 7 u c [e~l(x — xo(t))], where the spike location XQ — xo{i) is to be determined. For a fixed XQ we have shown that the linearization around this solution has an exponentially small principal eigenvalue as e -» 0. By eliminating the projection of the solution on the eigenfunction corresponding to this eigenvalue, we will derive an equation of motion for xo(t). We begin by linearizing around a moving spike solution by writing, a(x, t) = aE(x; xo{t)) + w(x, t), (2.34) where aE is defined in (2.1) and xo(t) is the trajectory of the spike. Since (2.8) does not have an 0(1) positive eigenvalue, we may assume that w « a £ and wt <S dtaE. Substituting (2.34) into (1.18), we get Lew = dtaE, - 1 < x < 1, .t > 0 (2.35a) wx(±l,t) = -dxaE(±l;x0). (2.35b) Next, we expand w in terms of the eigenfunctions </V of Lt as oo w = Y^Ei{t)<i>i. (2.36) The solvability condition for w is that w is orthogonal to the eigenspace of L* associated with the exponentially small eigenvalue. Let (f>* be the z t n eigenfunction of L*. Then, since (fa, (p*) = dij, we integrate by parts to show that Ei(t) = (w, $) = 1 [(Lew, ft) - e2wx<p*\U] , (2.37) 27 -0.395 -0.4 -0.405 XO -0.41 -0.415 -0.42 -0.425 0 28+07 4e+07 6e+07^ 88+07 1e+08 1.2e+08 1.4e+08 Figure 2.5: xo versus t for e = .05, xn(0) = — 0.4 and the parameter set (p, q,r, s) = (2,1,2,0). where L*4>1 = A|0*. Using (2.35), we have Ei(t) = ± [(dtasA*) + ^d^Elli] • (2.38) As discussed previously, when e <C 1, the nonlocal term in the eigenvalue problem (2.8) is insignificant in the asymptotic estimation of the eigenspace associated with the exponentially small eigenvalue of Le. Therefore, we can replace <j>\ and \ \ by <f>\ and Ai in (2.38), where 4>\ and Ai are given in (2.19) and (2.33), respectively. Since Ai -)• 0 exponentially as e —> 0, we must impose the limiting solvability condition that Ei = 0. This projection step yields the following implicit differential equation for xo(t): (dtaE, <t>\) = -e24>idxaE\-i • (2.39) The dominant contribution to the left side of (2.39) arises from the region near XQ. For e -> 0, 28 -0.56 ' 1 1 1 1 1 1 1 1 1 0 500000 1e+06 1.56+06 2e+0^ 2.5e+06 , 3e+06 3.5e+06 4e+06 4.5e+06 Figure 2.6: xo versus t for e = .06, xo(0) = —0.4 and the parameter set (p, q, r, s) = (2,1,2,0) we calculate {o\aE, <h) ~ -CihEx0pe-1, (2.40) where XQ = dxo/dt. Finally, we can evaluate the right side of (2.39) using our estimates for 4>\ (±1) in (2.28) and for uc(z) as z —> oo. This yields our main result of this section. Proposition 2.1 (Metastability) For e -> 0, metastable spike solution for (1.17), is repre-sented by a(x,t) — aE{x;xo(t)), where aE is defined in (2.1) and xo(t) satisfies, e -2 ( i - i 0 )A _ e-2(i+x0)Al _ (2.41) Here a and Q are defined in (2.5) and (2.20), respectively. x0{t) For a given initial condition x0{0) € (-1,1), this ODE shows that the spike drifts towards the endpoint that is closest to the initial location xn(0). As a consistency check on our solvability 29 -0.45 -0.46 --0.47 ' 1 1 1 1 1 1 1 L 0 20000 40000 60000 80000 100000 120000 140000 160000 180000 200000 Figure 2.7: XQ versus t for e = .07, XQ(0) = —0.4 and the parameter set (p, q,r, s) = (2,1,2,0). condition E\ = 0, we note from (2.20), (2.33), (2.40) and 2.41), that if the solvability condition were not imposed, then Eifa = 0(e~1^2), which would violate our linearization assumption. To verify the asymptotic result (2.41) we computed numerical solutions to (1.17) for various values of e for the parameter set (p, q, r, s) = (2,1,2,0). The computations were done by using a variable coefficient variable time step backward-differentiation (BDF) scheme to integrate in time (see [44] where a similar scheme is used). The boundary value problem solver, COLNEW (see [3]), was then used to solve the resulting boundary value problem at each time step. At each time step the solution is calculated using a third and fourth order BDF scheme to estimate the error and determine the next maximum allowable time step. Comparisons of these results with a numerical integration of the asymptotic ODE (2.41) may be found in Fig. 2.5-2.7. In the figures the solid curves are the numerical solutions and the dots are obtained from the asymptotic ODE. In computing numerical solutions to (1.18) the initial condition as(x;XQ(0)) 30 was used for certain values of XQ(0) as can be seen from these figures. 2.2 An n-Spike Solution We will now examine the properties of an n-spike quasi-equilibrium solution. The analysis will proceed in the same manner as for the case of the one-spike quasi-equilibrium solution. The stability of an n-spike quasi-equilibrium solution will be examined by linearizing about this solution and studying the resulting spectrum. We begin by defining an n-spike quasi-equilibrium solution by n - l -anMx) = hn,E^2uc [e 1ix~ xi)] » (2.42a) i-0 {e~1hLia"'EdxY+l' (2-42b) hN,E = ( e where 7 = q/(p — 1). Substituting (2.42a) into (2.42b), we can determine hn^ as e -> 0 as, /n/3\ («+i)(p-i)-«r h ^ = \ - j ) » (2-43) where ]3 was defined in (2.4). In (2.42a), the spike locations x\ for i = 0, ..,n — 1 satisfy —1 < XQ < x\,.., < xn-\ < 1. They correspond to local maxima of On,E-We now linearize (1.18) about aUtE and hnjE by introducing <f> and 77 defined by a{x,t)'= OnMx) + • (2.44a) h{x, t) = hntE + extn(x). (2.44b) Here <f> an,E and 77 hn,E- Substituting (2.44) into (1.18) we get the following eigenvalue problem, p—i p fy ^ 7 1 E E e 4>xx ~ <t> + P7j—<t> ~ Ir^iV = (2.45a) nn,E hlE 31 Since each spike of the quasi-equilibrium solution is localized to within an 0(e) region near x = Xi for some i , we look for an eigenfunction <f>(x) of the form n - l <j>(x) = YJ~4>iV-\x-Xi)) • (2-46) i=0 Therefore, we need to introduce local coordinates near each spike. In particular, the ith set of inner variables are defined as 4>i(y%) = <Pixi + ey»), Vi = e~l(x - xi). (2.47) Substituting (2.43) and (2.45b) into (2.45a) and switching to the localized coordinate system (2.47), we get the following system of eigenvalue problems, n - l fayiVi ~ & + P 53 Ucl€ 1(X~ Xj)] & 3=0 rqixul ^ '°° 2n0n(s + 1) " * poo 53/ u p 1 = A <&> ll/t|'< oo, (2.48a) $ 0 as yi ->• ± o o . (2.48b) Now we note that if each were independent of i (i. e. <j>i(yi) = $(yi)) for i = 0, . . ,n — 1, then E z^To1 X!ccuc-1'(yi)4>i(yi) dV% = n /!c^«c _ 1(y)*(y) T n e factor of n would cancel in (2.48) and we would be left with the same eigenvalue problem as (2.10). Thus, for the parameter set we have used previously, we would conclude that an n-spike solution has no 0(1) positive eigenvalue. However, we now show that this conclusion is erroneous. To see this we note that we can construct a global eigenfunction by taking (f>i(yi) = 6j$(y) for some constant 6$. The non-local term in (2.48) then becomes poo poo / n - l \ 53 / ur-\yi)Uvi)dyi = / ur-l{yMy)dy 5> . (2.49) i=0J-°° J-°° \i=0 J Then, if we impose the constraint that n - l 5> = 0, (2.50) i=0 32 the non-local term vanishes. Hence, with this constraint, <5(y) satisfies the local eigenvalue problem $" -§+pup-l§ = A 0 $ . (2.51) This problem has exactly one positive eigenvalue A n . Whenp = 2, we found that A n = 5/4 with corresponding eigenfunction $o(y) = sech2(y/2). Hence, under the constraint (2.50), An is also a positive eigenvalue of (2.48). This then leads to an instability. In summary, when there is more than one spike we may always construct an eigenfunction of the form <j)(x) = Yli^o ^ t 6 " 1 ^ — xi)] where J2i=o — 0. This eigenfunction has a positive eigenvalue. Therefore, it is impossible to find a stable multiple spike solution for the shadow problem. We now illustrate this instability result numerically for a two-spike solution for the parameter set (p,q,r,s) = (2,1,2,0), /j = 1, r - 0.01, D = 40, and e = 0.05. We took the quasi-equilibrium solution as our initial condition. The first spike (Spike 1) is centered at rrn = —0.5 while the second spike (Spike 2) is centered at x\ = 0.5. In Table 2.3 we tabulate the numerically computed amplitudes of the two spikes as a function of time. We now use this data to estimate the positive eigenvalue. We remark that the data in Table 2.3 is taken after the simulation has been run approximately t = 20 units to eliminate any transients and to ensure that the positive eigenvalue is dominant. After this time the solution at the spike locations x = XQ and x = x\ will be approximately given by, a{xi,t) « 02,E{XX) + eXot4>a(xi), i = 0,1. (2.52) This relation will only govern the linear instability of a.2,E- For the parameter set we have used Q2,E(%i) — 6-25. Then, we can re-write (2.52) as, Ao< +log [^o(a:i)]«log(|a(a:i,t)-6.25|), * = 0,1. (2.53) To estimate Ao from the data in Table 2.3, we take x\ = 0.5 and evaluate (2:53) at two different values of time, labeled by t\ and t2. Using the numerically computed values for a(0.5,t) at 33 t = t\ and t = £2 gives us two equations for the two unknowns 0n(O.5) and Ao- In this way, Ao can be estimated. In Table 2.4 we give the numerical results for Ao and </>o(0.5) using various values of t\ and £2- For this parameter set, we would expect that the principal eigenvalue is 1.25. The interpolated values, obtained by our numerical procedure, are all close to 1.25 as expected. 2.3 A Spike in a Multi-Dimensional Domain We now construct a quasi-equilibrium solution aE for (1.19). This is done in a similar manner as in the one-dimensional case, except that here the quasi-equilibrium solution will be radially symmetric about the center of the spike. Thus, we look for a steady-state solution to (1.19) in all of RN of the form a = a £ ( x ; x 0 ) =/i7u c(e_ 1|x - x 0 | ) , l = q/(j>-l), (2.54) where x<j is arbitrary. The function uc(p), called the canonical spike solution, is a non-negative radially symmetric function, which decays exponentially as p —> 00. It satisfies N — 1 v." + u'c -uc + upc = 0, p>0, (2.55a) u'c(0) = 0; uc(p) ~ ap(l-N)l2e-p, as p - > 0 0 , (2.55b) where a > 0 is some constant. In dimension N > 2, we require that p < pc, where pc is the critical Sobolev exponent for dimension N. In terms of this solution, h = hE, where P - I /' f~N C \ ( a + l ) ( P - l ) - ? ' -h*={mJa**) • (2-56) Here is the volume of ft. Since uc is localized near xo, as e ->• 0 we get, hE where fix is the surface area of the unit iV-dimensional sphere. Recall that in the one-dimensional case and with p = 2 we have the exact solution uc(p) = |sech2(|), and hence a — 6. To find numerical solutions for uc(p) and for a in other dimensions, 34 time Spike 1 Height Spike 2 Height 19.5 6.2738663390032 6.2635545772640 19.8 6.2761841264723 6.2612374542097 20.1 6.2795439171902 6.2578790593718 20.4 6.2844142872978 6.2530116219365 20.7 6.2914746492378 6.2459574213615 21.0 6.3017102226534 6.2357347923334 21.3 6.3165498360813 6.2209223724487 21.6 6.3380658088900 6.1994635220925 21.9 6.3692634678024 6.1683858288480 22.2 6.4144988374999 6.1234022942033 22.5 6.4800753144382 6.0583539884737 22.8 6.5750761790975 5.9644587136514 23.1 6.7124619645111 5.8293778760449 23.4 6.9103141671528 5.6362910167926 23.7 7.1926098041313 5.3636802602846 24.0 7.5876099073842 4.9877041819062 24.3 8.1196649411629 4.4906888214834 24.6 8.7900467006609 3.8780493118962 24.9 9.5540932480348 3.1943794842134 25.2 10.3232394145038 2.5150430348060 25.5 11.006159488840 1.9098035904776 Table 2.3: Height of spike 1 centered at XQ = —0.5 and of spike 2 centered at x\ = 0.5. 35 h *2 a(.5,<i) a(.5,*2) Ao to(-5) 22.8 23.4 5.9644587136514 5.6362910167926 1.275223721 —e -30.32846949 23.1 23.7 5.8293778760449 5.3636802602846 1.242238171 —e" -29.56172215 22.5 23.7 6.0583539884737 5.3636802602846 1.276189822 —e -30.36637629 22.2 23.4 6.1234022942033 5.6362910167926 1.315422050 —e" -31.26911038 Table 2.4: Logarithmic Interpolation of Ao and <f>o(.5). we will treat N as a real parameter, and use N (and p for p ^ 2) as continuation parameters. We can use the far field asymptotic behavior (2.55b) to obtain the boundary condition u'c = ^ 2 ^ uc, which we impose at some large value p = PL in our numerical computations of (2.55). The computations are done using COLNEW. In Fig. 2.8 we plot the numerically computed solutions uc(p) for N = 1,2,3 when p = 2. For the finite domain problem we restrict xo to be strictly contained in ft so that dist(xo, 9ft) >> 0(e). Then, under this restriction we note that, aE will satisfy the steady-state problem for (1.19a), but will fail to satisfy the no flux boundary condition (1.19c) by only exponentially small terms for any value of xo in the interior of ft. Thus, we expect that the spectrum of the eigenvalue problem associated with the linearization about aE contains exponentially small eigenvalues. 2.3.1 The Nonlocal Eigenvalue Problem Let xo 6 ft be fixed, and linearize (1.19) about aE, hE. We obtain the eigenvalue problem for this linearization by introducing <p and n defined by a(x, t) = a £ (x; x 0 ) + eA t (£(*), (2.58a) h(t) = hE + extn, (2.58b) 36 4.5 3.5 3 h 2.5 0.5 Figure 2.8: Numerical solution for uc(p) when N = 1,2,3 and p = 2. where (j> <S ag and rj <C h. Substituting (2.58) into (1.19) we obtain, after a lengthy calculation, the following non-local eigenvalue problem; Le(j) = e'A<f> + (-1 +pv%-1)<j> -rqe u£ f <~l Jn /3NnN(s +1) 4>n = 0 on dQ,. 0dx = A0, in ft, (2.59a) (2.59b) Here uc = uc [e *|x — xo|], and 3N is defined by 3N= / ur-xpN-ldp. Jo (2.60) Since uc is localized near xo, we will only seek eigenfunctions that are localized near x = xo-These eigenfunctions are of the form 0(y) = 0(xo + ey), y = e x (x - x 0 ) (2.61) 37 Therefore, we can replace Cl by W1 in (2.59a) and impose a decay condition for 4> as |y| -» co. This gives us the non-local eigenvalue problem for the infinite domain L4 = Ay~4>+{-l+pupc'l)~4>- R a ^ \ u / <~14>dy = \fa in RN, (2.62a) /5]Vl2iv(S + l j JRN 4> -> 0 as |y| -> co. (2.62b) In this problem uc = uc(\y\). If, in addition, we consider an eigenfunction that is radially symmetric (i. e. 0 = fap), where p = |y|), then (2.62) reduces to l€j> = Ap4>+{-l+pv?-l)j>- rqUl rUc-14>pN-1dp=~X4>, p>0, (2.63a) 4> 0 as p -¥ oo, (2.63b) where = </>" + (N — l)p~l(j>. We now analyze the spectra of these eigenvalue problems. We first note that, for each i = 1,.., iV, the function 0j = dyiuc(\y\) satisfies (2.62) with A = 0. Here yi is the coordinate of y. This follows from the combined effects of translation invariance and the vanishing of the integral in (2.62) by symmetry considerations. Thus, (2.62) has a zero eigenvalue of multiplicity N with corresponding eigenfunctions fa = dyiuc(\y\) for i = 1,.., N. Each of these eigenfunctions has one nodal line. These eigenpairs will be perturbed by only exponentially small terms as a result of the finite domain. Hence, there are N eigenvalues of (2.59) that are exponentially small, and they are estimated below. The goal is to determine whether these are the principal eigenvalues of (2.59). We claim that these are not the principal eigenvalues for (2.59) when the non-local term in (2.59) is absent. To see this, suppose that the non-local term in (2.59), (2.62) and (2.63) is absent. The corresponding eigenvalue problems are then local and self-adjoint, and several key properties follow. In particular, since fa = dyiuc(\y\) is an eigenfunction of (2.62) with a zero eigenvalue and has one nodal line, the local eigenvalue problem must have a simple, positive eigenvalue, which is independent of e. The corresponding positive, radially symmetric, eigenfunction satisfies the local version of (2.63). The effect of the finite domain in (2.59) is to 38 perturb this eigenvalue by only exponentially small terms. Thus, when the non-local term is absent no metastable behavior can occur. The effect of the non-local term will be to ensure that the exponentially small eigenvalues are the principal eigenvalues for the non-local eigenvalue problem (2.59). Thus, the quasi-equilibrium solution will be metastable if we can show that the principal eigenvalue of (2.63) has a negative real part. We may apply a multi-dimensional version of Theorem 2.1 to show that (2.63) has no eigenvalues with positive real part. Theorem 2.2 (Wei[52]; Multi-Dimensional) Consider the eigenvalue problem for 70 > 0, A4> - <f> + pvPc V - 7o(p - 1 K -^-7 T w = X & 2- 6 4) (f>'-*0 as \y\ 00. (2.65) corresponding to eigenpairs for which A ^ 0. Here uc satisfies (2.55). Let Ao ^ 0 be the eigenvalue of (2.64) with ^ e largest real part. Then i/70 < 1, i?e(A0)>0. (2.66) Alternatively, if jo > 1 and if either of the following two conditions hold, r = 2, l<p<5, or r = p+l,p>l. (2.67). Then, Re{X0)<0. (2.68) As in the case for one-dimension, for any parameter set satisfying (2.67), the operator Lt will have no 0(1) eigenvalues with positive real part. Again it is instructive to consider how the non-local terms in (2.63) effect the spectrum. To do so, we compute the eigenvalues and eigenfunctions of the radially symmetric problem (2.63), 39 where a continuation parameter 8, with 0 < 8 < 1, multiplies the nonlocal term L5j> = Apj> + (-1 + pup-l)4> - 5 (p 0 as p oo . (2.69b) We compute the eigenvalues of this problem as a function of 8, and in particular track the first eigenvalue Xo{8). We will show that the positive principal eigenvalue Ao(0), which occurs when the non-local term is absent, will cross through zero into the left half plane as 8 increases. Thus, we must show that the first eigenvalue Xo(S) has a negative real part when 8=1. For the parameter set (p,q,r,s) = (2,1,2,0), in Fig. 2.9 and Fig. 2.10 we plot the first two eigenvalues XQ(S) and A/v +i(5) of (2.69) as a function of 8 for N = 2 and N = 3, respectively. Here AJV+I is the first eigenvalue in the sequence for (2.62) following the zero A ,^ i = 1,..,N. These computations were done using COLNEW. These plots clearly indicate that Xo(8) crosses through 0 before 8 = 1. At some value of 8, XQ and X^+i collide and become complex. To track the eigenvalues past the point where they become complex, we use the same technique as in the one-dimensional case. The differential operator is approximated by a matrix and the eigenvalues of the matrix are then approximations of the eigenvalues of the differential operator. Using this numerical procedure, we give numerical values for the real and imaginary part of An(<5) in Table 2.5. This table shows that the real part of Ao is negative when 8 = 1. Similar computations, with similar conclusions, can be performed for other values of p, q, r and s. 2.3.2 A n Exponentially Small Eigenvalue We will now use a boundary layer analysis to construct a composite approximation to the eigenfunctions corresponding to the exponentially small eigenvalues of (2.59). The correspond-ing eigenfunctions are well approximated by dXiuc, for i = 1,.., N in the interior of the domain and each of these eigenfunctions has a boundary layer correction term near 8Q. in order to satisfy the no-flux boundary condition on dCl. In order to resolve the boundary layer we must define a local coordinate system. Let r) represent the distance from a point in Cl to dQ, where fj < 0 corresponds to the interior of Cl. Let £ correspond to the other N — 1 orthogonal coordinates. 40 6 A 0 in R2 A 0 in R3 0.00000 1.6388 2.3703 0.05000 1.4814 2.1588 0.10000 1.3231 1.9456 0.15000 1.1638 1.7304 0.20000 1.0030 1.5125 0.25000 0.84032 1.2910 0.30000 0.67516 1.0646 0.35000 0.50641 0.83098 0.40000 0.33218 0.58554 0.45000 0.14857 0.31741 0.50000 -.055026 -.019898 0.55000 -.37526 -.33843 + 0.29744? 0.60000 -.48239 + 0.24569? -.44368 + 0.45028? 0.65000 -.56115 + 0.33165? -.54978 + 0.54508? 0.70000 -.64059 + 0.38475? -.65696 + 0.60964? 0.75000 -.72097 + 0.41770* -.76550 + 0.65310? 0.80000 -.80268 + 0.43510? -.87584 + 0.67970? 0.85000 -.88640 + 0.43886? -.98857 + 0.69170? 0.90000 -.97333 + 0.42959? -1.1045 + 0.69037? 0.95000 -.10657 + 0.40726? -1.2249 + 0.67652? 0.10000 -.11678 + 0.37248? -1.3513 + 0.65089? Table 2.5: 5 and A 0 in R2 and R 3 for the case of (p, q, r, s) = (2,1,2,0). 41 1.5 A, 1 0.5 -0.5 0.6 Figure 2.9: An(o") and XN+I(5) versus 5' in IR2 for the parameter set (p,q,r,s) — (2,1,2,0). To localize the region near 8Q., we let rj = e~lfj. The eigenfunction on the finite domain can then be approximated by, fa = Ci (dXiuc (e x |x - x 0 | ) + 4>i} , (2.70) where Ci is a normalization constant and fa is a boundary layer correction term. Using the fact that uc is exponentially small near dCl, we get the following boundary layer problem 8mfa - fa = 0, 7] < 0, 8vfa = -edij(dXiuc)\v=0, on r? = 0, v v ' a function of ( fa -> 0 as rj —>• —oo. (2.71) (2.72) (2.73) We require that fa 0 as rj —oo to match to the outer solution. Then, the solution for fa is fa = m(()ev, where &(£) = -df,(dXiuc)\v=Q. (2-74) 42 2.5 1.5 0.5 h 0 h 0.3 0.4 0.5 0.6 Figure 2.10: An(£) and X^+i(S) versus 5 in R 3 for the parameter set (p,q,r,s) = (2,1,2,0). Thus, the composite asymptotic solution for the eigenfunction is <t>i = Ci dXiuc + egi(0e^ , i = l...,N. (2.75) Below we need an estimate for fa on dft. To do so we need to calculate <?j. Let xoi represent the i t n coordinate of X Q . S O , setting r = |x — xn|, we apply the chain rule, which gives 9i [uc(r/e)r • nj , (2.76) where n is the outward unit normal to Cl. Since uc(p) ~ ap^1 N)l2e~p as p -» oo we get that, g i --ae^-^ixi - x0i)r-^N^2e-^v n, on dQ. (2.77) Combining (2.75) with (2.77), we get an asymptotic approximation for fa on dQ, fa~-Ciae(N-W2a(xi-x0i)r-(l+N»2e-^(l + r.v), on 0f i . (2.78) 43 In order to complete our asymptotic estimate of the exponentially small eigenvalues, we apply Green's identity to fa and dXiuc to get thefollowing relationship: \{dXiUc, fa) = - e 2 / fadn(dXiuc)dS + [L* [dXiuc], fa). (2.79) Jan > an Here L* is the adjoint of Le, -iV„,r-l L*v = e2Av - v + up-1v - fc s f upv d x . (2.80) We will now estimate each term in (2.79). Since dXiuc is an exact solution to the local problem, we have that, L*e (dXiuc) = - ^ ' Y C u f « ? 3 X i u c d x . (2.81) Next, since uc is radially symmetric and localized to a small region in the interior of ft, it is clear that j^updXiucdK ~ fnupdXjucdx., as e —> 0 = 1...N. Thus, we may write the expression above as, L- ^ ~-N7£X» L (2-82» An application of the Divergence Theorem results in, L<^)~-N7~*1lln / (^)dS- (2-83) NpNnN{s +1) jdn yp + l y On the boundary of ft, u c [e _ 1 |x - x 0 | ] ~ a e ( ^ - i ) / 2 r ( i - ^ ) / 2 e - € - 1 | x - x 0 | Therefore, the integral in (2.83) will be exponentially small. We then estimate the integral in (2.83) to get the following bound: where = rqap+l , g 5 ) Nf3N£lN(S + l)(p + l) ' { ' Here po = dist(xo,dft). Therefore, with fa ~ CidXiuc, we have | ( I * (dXiuc),fa)\ < Fde^m^'1^'2 x p ( i - A r ) ( P + i ) / 2 e - e - i ( p + i ) P 0 f u r - i d x i U c d x _ { 2 M ) Jn 44 Estimating the integral term in (2.86) by a similar procedure, (2.86) reduces to the following, \{L*edXincAi)\ < aTmlFNCie'N^e(^-l)(?+r+l)/2p-(iV-l)(P+r+l)/2e-6-lpo(p+r+i) _ ^ Therefore, we conclude that \(L*dXiuc, <f>i)\ = 0 ^be-^Po(P+r+i)y i ( 2 . 8 8 ) for some b. We will show that this term is exponentially smaller than the first term on the right side of (2.79), and therefore, we can ignore it. Now we estimate the left hand side of (2.79). Since 0j and dXiuc are exponentially small outside of a neighborhood of x = xo, this inner product is dominated by the contribution from the region near x = xrj. Using a Laplace-type approximation, we can approximate the inner product to get {dXiucAi) ~^fQ[<(r/e)}2 ^^)2d*~^W^fRN [<{p)?PN-ldPd6, (2.89) where 6 represents the N — 1 angular co-ordinates. Since the integrand is independent of 9, rco {dXiuc,fa) ~CieN-2nNpN/N, where pN = / [u'c{pyf pN~l dp. (2.90) Jo Here ft JV is the surface area of the iV-dimensional unit sphere. Then we determine Cj by using the normalization relation J n (f>2 dx. = 1 to obtain, Ci = (JL-]l/\{2-N)/2_ ( 2 9 1 ) Finally, we get our asymptotic estimate of A* by substituting (2.90) and (2.78) into (2.79), and using the estimate dn(dXiuc) ~ ae^N~5^2r~^N+1^2e~r^, on 9ft. In this way, we get Xi ~ I ^ ~ ^ ) 2 r - ( 1 + A r ) e - 2 r / £ ( r • n)(l + r • n) dS, (2.92) where r = (x - xrj)/V and r = |x - xo|, with x E dft. As a consistency check we use (2.88) to observe, by comparing the asymptotic orders of the two terms on the right side of (2.79), that the second term is asymptotically negligible compared to the first term, since the exponents satisfy p + m + 1 > 1. 45 The surface integral in (2.92) can be evaluated asymptotically by using a multi-dimensional Laplace technique. Assume that there exists a unique point x m € 9ft where rm = dist(x0,9ft) is minimized. If we parameterize the boundary near £ m (where x(£ m ) = x m ) such that each Q corresponds to arclength along one of the principal directions through £ m , then for any smooth F{r), we have (see [48]), / rl-NF{r)e-2Tl< dS = [ — ) F(rm)H(rm)e-2r^, (2.93) Jan \rmJ where H(rm) = (1 - rm/Rx)-ll2{l - rm/R2)-1/2 • • • (1 - rm/RN^)-^2 . (2.94) Here Rj > 0, for j = 1,.., N — 1 are the principal radii of curvature of 9ft at x m . This result assumes that the non-degeneracy condition Rj > rm, j = 1,... ,N — 1 holds. In this way, we obtain the following explicit asymptotic estimate for the exponentially small eigenvalue, X t ^ ( ^ f - 1 ) , 2 f ^ ) 2 m r m ) e - ^ t „ 9 5 ) 6N£lN VmJ V rm J where ej is the standard unit basis vector in the direction and r m = (xm — xo)/rm. 2.3.3 The Slow M o t i o n of the Spike Now we derive an ODE characterizing the metastable spike dynamics. We first linearize (1.19) about a moving spike by writing, a(x,i) = a^(x;x 0(i)) + w{x,t), (2.96) where aE(x;xo) is defined in (2.54). Since, (2.59) does not have an 0(1) positive eigenvalue, we may assume that w < aE, wt <C dtaE uniformly in time. Substituting (2.96) into (1.19), we obtain Lew = dtaE, in ft, (2.97a) dnw = —dnaE, on 9ft. (2.97b) 46 Next, we expand w in terms of the eigenfunctions fa of Le as oo w = Y/Ei(t)fa. (2.98) i=0 We assume that the eigenfunctions form a complete set. However, this is not required for the construction of the solvability condition as the key requirement is that w is orthogonal to the eigenspace of L* associated with the exponentially small eigenvalues. Let <j>* be the eigenfunction of L*. Then, since (fa,<f>j) = we integrate by parts to show that Ei(t) = (w,fa) A* {Ltw,fa)-e2 f wnfadS Jan (2.99) where L*$ = \%<j)*. Using (2.97), we have 1 Ei(t) = A? (8taE,fa) + e2 [ 8naEfadS Jon (2.100) As seen in (2.79)-(2.88), the nonlocal term in the eigenvalue problem Le(j> = \<f) is insignificant when e <§; 1 in the asymptotic estimation of the eigenspace associated with the exponentially small eigenvalues of Le. Therefore, for i = l,..,N and e -> 0, we can replace fa{ and A* by fa and Aj in (2.100), where fa and \ are given in (2.75) and (2.95), respectively. Since Aj —> 0 exponentially as e -> 0, for i = 1, ..,N, we must impose the limiting solvability conditions that E{ = 0, for i = l , . . , iV . This projection step yields the following implicit differential equation for xrj(i): (dtaE, fa) f dn Jan aEfa dS. (2.101) The dominant contribution to the left side of (2.101) arises from the region near x 0 , and we calculate Ch? [dtaE, fa) j^x0ittNJ3NeN-2 . (2.102) Finally, we can evaluate the right side of (2.101) using our estimates for fa on 80, in (2.78) and for uc(p) as p —y co. This yields the main result of this section. 47 Proposition 2.2 (Metastability) For e ->• 0, a metastable spike solution for (1.19), is rep-resented by a(x,t) = a£ ; (x ;xo ( t ) ) , where aE is defined in (2.54) and xn( i ) satisfies, XL0~4^-[ f r 1 - i V 2 r / 6 ( l + f -n)f -ndS. (2.103) Here f = (x — xo) r - 1 , r = |x — xo|, x £ 9ft, and n j's the unit outward normal to 9ft. In addition, a and PN are defined in (2.55b) and.(2.90), respectively. There are a few corollaries that follow from this result. Corollary 2.1 (Equilibrium) For e -> 0, an equilibrium solution for (1.19), is represented by a(x, i ) = a # ( x ; x o e ) , where x n e is a root of I(XQ), where 7(x 0 )= f rrl-Ne-2r/'{l+r -nJf-ndS (2.104) Jan It was shown in [48] that for a strictly convex domain xoe is unique and is centered at an 0 ( e ) distance from the center of the uniquely determined largest inscribed sphere for ft. This equilibrium solution is unstable. Assuming that there is a unique point x T O G 9ft closest to the initial center xo (0) of the spike, we can evaluate the surface integral in (2.103) using Laplace's method to get the following explicit result: Corollary 2.2 (Explicit Motion) Let x m be the point on 9ft closest to x n ( 0 ) . Then, for t > 0, and e —> 0, the spike moves in the direction ofxm and the distance rm(t) = | x m — xn(i)|, satisfies the first order nonlinear differential equation (JV+l)/2 £ r m ^ j H(rm)e-2r™/£, (2.105) where (2.106) ftiv/3jv Here (3N is defined in (2.90) and H(rm) is determined in the terms of the principal radii of curvature o/9ft at xm given in (2.94). 48 This result is valid up until the spike approaches to within an 0(e) distance of x m . If the initial condition for (2.105) is rm(0) = rn, then the time T needed for rm(T) = 0, is readily found for e -> 0 to be, e ( l - i V ) / 2 r ( " - l ) / 2 (2.107) Once the spike reaches the boundary, it moves in the direction of increasing mean curvature until it reaches an equilibrium point where the mean curvature of the boundary has a local maximum (see [20]). The existence of such equilibrium solutions, where the spike is located at these special points on the boundary, is demonstrated rigorously in [13], [35]. 49 Chapter 3 Spike Motion on the Boundary When the spike approaches to within an 0(e) distance from the boundary, the analysis leading to Proposition 2.1 is no longer valid, and the spike presumably begins to merge with the boundary. This process is difficult to study in the multi-dimensional case without a full-scale numerical analysis. In Fig. 3.1 we show this merging process in the simpler case of one dimension by solving (1.18) numerically with p = 2 and e = .07 using the method described in §2.1.3. Results concerning the stability of an equilibrium boundary spike in one dimension for (1.18) are given in [54]. The merging process of a spike with the boundary should probably be similar in higher dimensions. For the equilibrium problem, the existence of boundary spike solutions to the multi-dimensional problem (1.19) has been proved in [13] and [35]. In particular, the result of [13] proved that there exists a solution to (1.19) where the spike is centered at a local maximum of the mean curvature of the boundary of a three-dimensional domain. The goal of this chapter is to analyze the motion of a spike solution for the non-local shadow problem (1.19) when the spike is confined to the smooth boundary of a two or three-dimensional domain. Since all spikes that are initially located in the interior will tend towards the boundary of the domain (cf. Proposition 2.1), it is natural to study the motion of a spike on the boundary as it tries to reach an equilibrium where its associated energy can be minimized. We assume that the merging process of an interior spike with the boundary has taken place and thus our initial condition is a spike located at an arbitrary point on the boundary. From using a formal, asymptotic analysis combined with imposing appropriate solvability conditions on the linearized problem, we derive differential equations characterizing the motion of a boundary spike. This 50 0.6 Figure 3.1: A plot of u versus x at different times showing a spike merging with the boundary in one dimension. Times from the simulation are t = 389,1354,1375,1383,1386 from right to left. motion generically occurs on a slow time scale of 0(e3). From this differential equation for the spike motion, we show that the spike drifts towards a local maximum of the curvature in two dimensions and a local maximum of the mean curvature in three dimensions. Again the non-local term in (1.19) is essential for ensuring the existence of this slow motion. In the derivation we assume the boundary is sufficiently smooth so that the derivatives of the curvatures exists. The differential equations mentioned above predict no motion when a spike is on a segment of the boundary having constant curvature. To illustrate the spike motion in this case, we analyze the motion of a spike on a flat boundary segment of a two-dimensional domain (see Fig. 3.6 below for the geometrical configuration). For this case, we show that the motion is metastable and depends critically on the local behavior of the boundary near the corner points at the two ends of the flat segment. A similar analysis for the constrained Allen-Cahn equation of material 51 science has been studied in [43] (see also the references therein). The remainder of this chapter proceeds as follows. In §3.1 and §3.2 we analyze the motion of boundary spikes for (1.19) in two and three dimensions, respectively. In §3.3, we examine the stability properties of the equilibrium boundary spike solutions found in §3.2 and §3.2 by using the results of [52] and [53]. Finally, in §3.4 we analyze the metastable behavior of a boundary spike in two dimensions that lies on a flat segment of the boundary. In this chapter we will be using r as a radial variable for polar and spherical co-ordinates. In order to avoid confusion with the parameter r in the equation (1.19), for this chapter only we write the equations as, The only difference between (1.19) and (3.1) is the substitution of m for r. 3.1 Spike Motion on the Boundary in Two Dimensions We now derive an asymptotic differential equation for the motion of a spike confined to the boundary of a two-dimensional domain. The boundary of the domain is assumed to be suffi-ciently smooth so that the curvature and its derivatives are differentiable functions. To derive this differential equation we first transform (3.1) to a localized boundary layer coordinate system centered near the spike. The solution is then expanded in powers of e. A nontrivial solvability condition for the 0(e2) equation in this expansion is obtained when we choose the time scale of the spike motion to be 0(e 3). The differential equation for the motion of the spike is obtained from this solvability condition. We now give the details of the analysis. We first introduce a boundary layer coordinate system, where n > 0 denotes the distance from x 6 0 to 80 and where s is arclength along 80. In (3.1b) (3.1a) 8na 0 on 80. (3.1c) 52 terms of these coordinates, (3.1a) transforms to at = e2 {a™ ~ r^a»+ ih^ds {\~^as) )~a+aP/hq - (3-2a) av = 0, on n = 0. (3.2b) Here K = K(S) is the curvature of the boundary and h = h(t) is given in (3.1b). Next, we introduce u(r], S, t) by a = / i 7 u , 7 = g/(P~l)- (3-3) Substituting (3.3) into (3.2) we obtain u , = 0, on ?| = 0. (3.4b) Suppose that the spike is initially located on the boundary of the domain. Then, its subsequent location is given by s = so(t) and 77 = 0, where so{t) is to be determined. Since the spike has a support of order 0(e) near so(t), we introduce local variables v, s and 77 by 5 = e _ 1 [s - s0(t)} , 77 = e - 1 7 7 , v(fj,s) = u(efj,s0 + es,t). (3.5) Then, (3.4) transforms to 'yvht _i ' en 1 / 1 \ — e ^ Q V J = Vffi - -Vfj + zd-s zVs - v + vp, (3.6a) Al 1 — 6K77 1 — € « 7 7 \ 1 — eKT] J Vjj = 0, on 77 = 0, (3.6b) where s 0 = dso/dt and K = K ( S O + es). Since the boundary was assumed to have a well-defined tangent plane at each point, the domain of definition for (3.6) is n = {(77,s) | -oo<s<oo ,77>0}. (3.6c) From differentiating the integral in (3.1b) with respect to t, it follows that ht = 0 ( e - 1 ^) . (3.7) 53 Thus, the two terms on the left hand side of (3.6a) are of the same order in e. However, as we show below, only one of these terms will contribute a nonzero term to a solvability condition. The solution to (3.6a), (3.6b) is expanded as v — VQ + ev\ + e2V2 H , (3.8) and the curvature is expanded as K(SQ + es) = ACO + CSK'Q + 0(e 2). (3.9) Here we have defined KQ = K(SQ) and K 0 = K (SQ). Since a nontrivial solvability condition arises at order 0(e2), we must choose a slow time scale r by r = e3t. (3.10) Substituting (3.8)-(3.10) into (3.6), and collecting powers of e, we obtain the following sequence of problems that are to be solved on the half-plane Q: vofjfj + voss + Q{v0) = 0, (3.11a) Cvi = Vyffi + VlSs + Q' (vo)Vl = KQVOjj ~ 2fJK0V0ss , (3.11b) CV2 = ^ T j + V2ss + Q' (V())V2 = -SQVOS + + Fe + F0 . (3.11c) Here so = dso/dr. The boundary conditions for (3.11a)-(3.11c) are vofj = vifj - v2rj = 0, on TJ = 0. (3.lid) In (3.11) we have defined Q{VQ) and its derivatives by Q>o) - -vo + vp0 , Q'(v0) = -1 , Q"(^0) = Pip - 1 K ~ 2 • (3-12) The terms Fe and F0 in (3.lie) are defined by Q"ivo) Fe = K0vifj + T?KOUOTJ - ZrJKQViss ^—vl - Snffivoss > (3.13a) F0 = sKQvofl - TJKQVQS - 2TJSK'0VOSS • (3.13b) 54 The leading-order problem (3.11a) has a unique positive radially symmetric solution uc(p) that satisfies (see [35]) u'c + -u'c + Q{uc) .= 0, 0 < p < o o , (3.14a) P u'c(0) = 0; uc{p) ~ acp~l/2e~p, as p -> co. (3.14b) Here a c is some positive constant and p = (rp + s2) . Thus, we take v0(v,s) = uc [(fj2 + S 2 ) 1 / 2 ] , (3.15) which also satisfies the boundary condition in (3.lid). To obtain our solvability condition, we notice that the tangential derivative vos satisfies CVQS = 0 and the boundary condition (3.lid). Hence, upon defining the inner product (/, g) = f^fgdx, we must have that the right-hand sides of (3.11b) and (3.11c) are orthogonal to VQS with respect to this inner product. An important observation is that vos is an odd function of s. This solvability condition for (3.11b) yields, «o (vofj, vos) - 2K0 {fjvoss, v0s) = 0. (3.16) Since U O T J and VQSS are even functions of s, the integrands associated with the inner products in (3.16) are odd, and hence the left hand side of (3.16) vanishes identically. Then, we can solve (3.11b) for ui and obtain that v\ is even in s. Next, upon applying the solvability condition to (3.11c), we obtain so (vos, von) = (F0, VQS) + (Fe, vos) + -j- (vo, vos) • (3.17) The significance of the decomposition in (3.11c) is that Fe is even in s, whereas F0 is odd in s. Hence, the last two terms on the right-hand side of (3.17) are zero. Next, substituting (3.13b) into (3.17), we get so (vos, vos) = «o (5yo»j, vos) ~ K'Q (rjvos, vos) - 2K'q {fjsvossiVos) • (3.18) 55 Finally, the inner products in (3.18) are evaluated exactly using polar coordinates to get (5U(tf},U0s) = (W)S,l>0s) = 3 / P2 «'c(P) («o», «os) = 2 Jo P[uc(p)\ DP' 2(T?5V0IS,VOI) = . - - y P2 U'C(P) dp. dp, Substituting (3.19) into (3.18) we obtain the following main result: (3.19a) (3.19b) (3.19c) Proposition 3.1 (Boundary Motion in 2 Dimensions) For e 0, the motion of a spike confined to the smooth boundary of a two-dimensional domain is described by a ~ K>uc (e" 1 [(s - s0(t))2 + n2]1/2) + 0(e), SoW ~ « («o) , w/iere 7 = q/(p — 1) and 6 > 0 is defined by I0°°P2U(P) (3.20a) (3.20b) 6 = dp /O°°PK(P)] DP (3.20c) Here uc(p) is the positive solution to (3.14) and « is the curvature of the boundary, with K > 0 for a circle. From (3.20b), we observe that the spike will move on the boundary in the direction of increasing curvature until a local maximum of the curvature is reached. The stable steady-states of (3.20b) are at the local maxima of the curvature. A similar differential equation has been derived in [1] for small bubble solutions of the constrained Allen-Cahn equation. 3.1.1 A Few Explicit Examples We now illustrate (3.20b) in a convex domain. Let the origin be contained in O. and let (xi,x2) be a point on dCl. Let £ denote the perpendicular distance from the origin to the tangent line to dQ, that passes through (x\,x2). Let 0 denote the angle between this perpendicular line and 56 the positive x\ axis. Then, when 9 ranges over 0 < 9 < 2 T T , we sweep out a closed domain Q, whose boundary is given parametrically by (see [17]) Xl(9) = C(0) cos(0) - C'(0) sin(0), x2{9) = ({9) sin(0) + ('{9) cos{9), (3.21) Here ((9) is 2TT periodic. Next, we transform the ODE (3.20b), written in terms of arclength, to one involving 9. Let s — /(#) be the mapping between 9 and the arclength s. Then, f'(0) and the curvature of the boundary K(9) are given by - l /(0) = C(*) + C (<?), «(*)= C(*) + C (*) (3.22) Hence, (3.20b) transforms to 90(T) ~ -46 C(0o) + C">o) (3.23) 3 7 r [CW + Oo)] 4 Here #0(7") is the value of 0 at the center of the spike, and r = e3t is the slow time variable. Using the boundary value problem solver COLSYS [3] we can solve (3.14) numerically to de-termine the constant b in (3.20b). In the examples below we took p = 2. For this value, we compute that f «C(P) 1 2 dp = 4.23, Jo uc{p) dp = 2.47. (3.24) Hence, when p = 2, we get b = 1.71. In the examples below, solutions to (3.23) were computed using the Sandia ODE solver [42]. Example 1: Let ({9) = 3 + 1.2sin3(0), and take the initial condition for (3.23) as 0O(O) = -1.5, for which x\ = 0.339 and x2 = —1.790 at r = 0. The curvature has three local maxima for this domain. In Fig. 3.2 we plot the domain bounded by ((9) and show snapshots of the motion of the center of the spike (labeled by the starred points) towards the nearby local maximum of the curvature. In Fig. 3.3 we plot the numerical solution to (3.23) with 0n(O) = —1.5. For this initial value and with e = 0.1, this figure shows that it takes a time t = r /e 3 « 77500.0 to reach the steady-state value. 57 -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 Xi Figure 3.2: For Example 1 we plot the motion of the center of the spike on the boundary at different times as it tends to its steady-state limit. The initial point is labeled and the times corresponding to the other points (in counterclockwise order) are r = 34.49, T = 58.71, r = 74.01, r = 76.98 and r = 77.50. 58 -0.50 -0.75 h 0O -1.0 --1.25 h -1.5 0.0 20.0 40.0 60.0 80.0 r Figure 3.3: For Example 1 we plot the solution 6>0 versus r to (3.23) showing the behavior towards a local maximum of the curvature. Example 2: Let £(0) = 3 + cos(50)/lO, and take 0O(O) = 0.6, for which x\ = 2.430 and x2 = 1.567 at r = 0. For this case, the ODE (3.23) becomes Hence there are five local maxima of the curvature. In Fig. 3.4 we plot the domain bounded by ((0) and show snapshots of the spike motion towards the nearby local maximum of the curvature at Oo = 0. In Fig. 3.5 we plot the numerical solution to (3.25) with 0O(O) = 0.6. The apparent nonsmoothness of the graph of 9o versus r near the equilibrium point results from the constant. A similar explanation hold for the apparent nonsmoothness in Fig. 3.3. 3.2 Spike Motion on the Boundary in Three Dimensions We now derive an asymptotic differential equation for the motion of a spike confined to the boundary of a three-dimensional domain. The boundary of the domain is assumed to be suffi-ciently smooth so that the principal curvatures and their partial derivatives are differentiable (3.25) fact that the linearization of (3.25) near 0o = 0 has the form 0' « —cf?, where c > 0 is a large 59 -4.0 -3.0 -2.0 -1.0 0.0 1.0 2.0 3.0 4.0 Xi Figure 3.4: For Example 2 we plot the motion of the center of the spike on the boundary at different times as it tends to its steady-state limit. The initial point is labeled and the times corresponding to the other points (in clockwise order) are r = 23.40, r = 31.41, r = 35.38, and r = 35.5. 60 0.8 0.6 #0 0.4 h \ 0.2 0.0 -0.1 0.0 10.0 20.0 30.0 40.0 50.0 T Figure 3.5: For Example 2 we plot the solution #0 versus r to (3.25) showing the behavior towards a local maximum of the curvature. The initial condition was #o(0) = 0.6. functions. In order to evaluate the Laplacian on the boundary we will use boundary layer co-ordinates. Lines of curvature form a local orthonormal basis for a coordinate system restricted to the boundary. We may then extend this system locally using the normal to the boundary as our third coordinate. The formulation of the Laplacian operator using these coordinates is not as simple as for the two-dimensional case. However, since the spike is localized on the boundary, we do not need an exact expression for the Laplacian in terms of these coordinates as only the first few terms in the local expansion will suffice for the analysis. We introduce a boundary layer coordinate system (si,S2,r]), where n > 0 is the distance from x € ft to dQ and where si and s2 correspond to coordinates through the two principal directions at the center of the spike. The boundary spike is assumed to be located at s\ = £i(i) 61 and 52 = with 77 = 0. Using Appendix B, we obtain that (3.1a) transforms to 2 2 ( K i K2 \ e2 Q (\,— T\K2 A T = E A ^ - E + T ^ J a " + ( i - ^ ) ( i - , K 2 ) * • > 1 n—a 1 ""Si 1 - 7/Kl 71 T71 r9 S 2 ( f - ^ a S 2 ) - a + 0?lb? , (3.26a) av = 0, on n = 0. (3.26b) Here K\ = K i ( s i , s 2 ) and K 2 ( S I , S 2 ) are the two principal curvatures at each point on the boundary. As in (3.3) we set a = h^u, where u = U(T],SI,S2,£). Then, we introduce local coordinates v, s"i, S2 and 77 by «i, «2 , t) = u(efj, £1 + esi, 6 + es2, t). (3.27) The estimate (3.7) and the time scale (3.10) still apply in the three-dimensional case. Next, we expand v as in (3.8). The principal curvatures are also expanded in the Taylor series «i(fi + esi,£2 + th) = «i + e s i K n 4- es2«i2 + 0(e2), (3.28a) «2(£i + esi,£ 2 + "2) = K 2 + esiK2i + C S 2 K 2 2 + 0(e2). (3.28b) Here on the right-hand side of (3.28a) and (3.28b) we have defined K\ = K I ( £ I , £ 2 ) > «2 = «2(£i ,6), and Kij = dSiKi{si,S2) evaluated at (si,s 2) = (£1 ,6) . Substituting (3.7), (3.8), (3.10), (3.27) and (3.28) into (3.26), we obtain the following sequence of problems upon collecting powers of e: V07777 + vo i ia i + VQS2S2 + Q(v0) = 0, (3.29a) Cvi = vifjfj + v U l h + vu2s2 + Q'(V0)VI = («i + K2)V0TJ - 2r ]«; iVos 1 i 1 - 2fJK2v0s2s2, (3.29b) Cv2 = v2fjjj + V2S1S1 + V2s2s2 + Q' (VO)V2 = FEE + FEO + FOE . (3.29c) The boundary conditions for (3.29) are voij = ni} = v2j) = 0, on 77 = 0. (3.29d) 62 In (3.11), Q{VQ) is defined in (3.12). The terms Fee, Feo and Foe in (3.29c) are defined by 1_», . o , . o oN evojh Fee = - ( u o) v l + ( K l + «2)«lJj + « ( « 1 + « 2 ) u 0 T J H - 2f]Kivu1h - ^K2Vis2s2 ~ 3«i7?2uo5iSi - ^l^^hh > (3.30a) •Feo = ( « 1 2 + « 2 2 ) s 2 U 0 7 7 ~ S ^ K i a ^ V O s i S ! ~ ^V^hVOhh - rK«22 - K 1 2 ) V 0 S 2 - i\VQs2 , (3.30b) •Foe = ( « 2 1 + Kll)SlVOfj - ^TJKxiSiVOs^! ~ 2fJK2lSlVo3252 - - K2l)v0s1 - &VQ5i • (3.30c) Here = d£j/dr for j = 1,2. The problems in (3.29) are to be solved in the half-space ft defined by ft = {(TJ, SI, s2) | — co < 5 i < oo, — co < S2 < oo , fj > 0} . (3.31) When p is less than the critical Sobolev exponent p c = 5, there is a unique positive radially symmetric solution uc(p) to (3.30a) that satisfies (see [35]) u'c + -u'c + Q(uc) = 0, 0 < p < oo, (3.32a) u'c(0) — 0; uc(p) ~ acp~le~p, as p oo, (3.32b) for some ac > 0 where p = (r]2 + s 2 + s2) • Therefore, our leading-order spike solution is given by vo{f], 3i,h) = uc [{fj2 + 3? + 4)l/2] • (3-33) To obtain our solvability condition, we first define the inner product (/, g) by (/, g) = fQ dx where ft is the half-space defined in (3.31). Then, we note that CVQSI — 0 and £VQS2 = 0, where £ is the operator defined in (3.29b), and that VQS-1 and VQ§2 satisfy the boundary condition in (3.29d). Therefore, the solvability condition is that the right hand sides of (3.29b) and (3.29c) must be orthogonal to both vo^ and VQS2 with respect to this inner product. In imposing the solvability condition on (3.29b), we note that the right-hand side of (3.29b) is even in both s~i and S2, whereas VQS1 is odd in Si and even in S2, while uos2 is even in si and 63 odd in S2 - Hence, the solvability condition for (3.29b) is automatically satisfied. Then, from (3.29b) together with (3.29d) we can calculate a function vi, which is even in both 3\ and sV Next, the solvability condition for (3.29c) yields the two equations {Fe^Vos,) + {Feo,V0-Sl) + t>05i) =0, {Fee,v0s2) + {Feo,v0s2) + (Foe,v0s2) =0 . (3.34) The significance of the decomposition in (3.29c) is that Fee is even in both Si and s2, Fe0 is even in si and odd in 52, and Foe is odd in 3\ and even in s2. Therefore, the inner products involving Fee vanish, and also (i^o, uosj = {F0e,vos2) = 0. Using these results, and substituting (3.30b)' and (3.30c) into (3.34) we obtain £l (voaj.VOSi) = ( « 2 1 + « l l ) (hV0rj,V0Sl) - 2 K U {fjSiVOsih^Osi) - 2 K 2 I (fjsivohh^oh) ~ ( K n - K 2 l ) (fjVOsi,VQh) -£ 2 ( vos 2 , uos 2 ) = ( « 1 2 + K 2 2 ) (s2vofi,vos2) - 2 « ! 2 (vhvoSih, V0s2) - 2 K 2 2 (vhVQS2S2,V0h) - ( K 2 2 ~ K 1 2 ) (VV0s2,V052) . We now evaluate the inner products in (3.35). We first integrate by parts to get 2 (fjsjVosjSj, VOSJ ) = - {fj, [VQSJ ]2) , j = 1,2 • Next, we use spherical coordinates to obtain ^ poo . -2 {fjvos^vosj) = 7 / p 3 [«c(p)J dp, J = l , 2 , [VOij.VoSj) = Y JQ p ( ^ V o i j S i . U O J a ) = (^ 5 l u 05 2 Sa,W0«i) = g y p 3 t*'c(p) " c(p)J dp, J = 1,2, u'c(p)l dp, j = l , 2 , dp. Then, substituting (3.36) into (3.35) we obtain, ki = ^ ( « 2 1 + «l l ) , 6 = ^ ( « 2 2 + K 1 2 ) , where b > 0 is defined by 6 = / o ° ° P 3 « ' C ( P ) dp (3.35a) (3.35b) (3.36a) (3.36b) (3.36c) (3.36d) (3.36e) (3.37a) Jo™ P2 Kip)} d P (3.37b) 64 Finally, upon introducing the mean curvature H(£\,£2) denned by H — (K\ + K2)/2, we obtain the main result. Proposition 3.2 (Boundary Mot ion in 3 Dimensions) For e -> 0, the motion of a spike confined to the smooth boundary of a three-dimensional domain is described by a ~ K*ue (e"1 [(8l - Ut))2 + (S2 - 6 W) 2 + " 2 ] 1 / 2 ) + 0(e), (3.38a) £ ' ( £ ) - ^ e 3 V t f ( £ ) . (3.38b) Here £ = ( £ L , £ 2 ) > 7 — l/(p — 1)> b > 0 is defined in (3.37b), uc(p) is the positive solution to (3.32), and H is the mean curvature of dQ,, with H > 0 for a sphere. Notice that the stable equilibrium points of (3.38b) are at local maxima of the mean curvature. Using the boundary value problem solver COLSYS [3] we can solve (3.32) numerically to de-termine the constant b in (3.38b). In particular, when p = 2 we compute that r°° r 1 2 r°° r 1 2 / p 3 u'c{p) dp =17.36, / p2 u'c(p) dp = 10.42. (3.39) Jo 1 J Jo L J Hence, when p = 2, we get b = 1.67. 3.3 Qualitative Properties of the Associated Eigenvalue Prob-lem In this section we qualitatively explain why the non-local term in (3.1) is essential for ensuring the existence of slow boundary spike motion. We first consider the local problem corresponding to (3.1) in which we delete (3.1b) and fix h > 0. Let us suppose for the moment that the boundary is flat and is given by the coordinate line XN = 0. Hence we take Q, = {x = (x\,... ,XN) | XN > 0}. Setting a = hq^p~^u, we then get ut = e2Au - u + up , x /y>0; dXNu = 0, on xN = 0. (3.40) 65 Let u = uc (e - 1 |x |) , where uc is the canonical spike solution defined in (3.14) and (3.32) for N = 2 and N = 3, respectively. We linearize (3.40) around uc by writing u = uc (e _ 1 |x|) + (e _ 1 |x|) e^. This leads to the local eigenvalue problem for ip(y) and p on a flat boundary A > + (-l+pup-l)i> = nil;, yN>0, (3.41a) dyN1> = 0, on y,v = 0. (3.41b) Here y = e _ 1 x, uc = uc(jyj) and A ' denotes the Laplacian in the y variable. The eigenvalues Hj of (3.41) satisfy Mi > 0, M2 = • • • = UN = 0, MJV+I < 0. (3.42) The positivity of / i i was shown in [27] and [48]. The eigenfunctions V>i and ipN+i are radially symmetric (i. e. tp\ = ipi{\y\)). A numerical procedure to calculate / i i was given in [48] and the results are shown in Table 1. The translation eigenfunctions corresponding to the zero eigenvalues are given by ipj = 9J,^_1uc(|y|) for j = 2,... ,N. p Mi(iV = 2) M i ( iV = 3) 2 1.65 2.36 3 5.41 15.29 4 13.23 144.18 Table 3.1: Numerical results for the principal eigenvalue pi of the local problem on a flat boundary (3.41). The local eigenvalue problem on a curved boundary has the form A€i>€ + (-1 + puP'1) ipe = pfrlf, ??>0, (3.43a) dfjipe = 0, on fj = 0. (3.43b) Here Ae is the operator that results from converting the Laplacian into local boundary coor-dinates as explained in §2 and §3. Clearly, At A ' as e —>• 0. Hence, we would expect that the eigenvalues of (3.41) and (3.43) are close as e ->• 0. In fact, it was proved in [53] that, for 66 e < l , the eigenvalues pj of (3.43) satisfy ^• = ^ • + 0(1), ; = 1,... ,iV + l , . . . . (3.44) Hence, n\ > 0 for e sufficiently small. This shows that a boundary spike solution for the local problem corresponding to (3.1) will not drift slowly along the boundary of the domain. The eigenvalues • • • >M/V corresponding to the near translation modes were calculated for e <C 1 in [53]. Next, we linearize the non-local problem (3.1) around a spike solution centered on a flat bound-ary. In place of (3.41), we obtain the non-local eigenvalue problem for A and </> given by A > + (-l+pup-l)<f> + I{4>) = \<f>, yN>0, (3.45a) dVN<t> = 0, on yjv = 0, (3.45b) where I(<p) is defined by Here Cl is the half-space y^ > 0, 0,^ is the surface area of the unit iV-dimensional sphere, and uc = « c(|y|)-As is similar to (3.44), the eigenvalues of the non-local problem defined on a curved boundary should be asymptotically close to within o(l) terms to the eigenvalues of (3.45). Hence, to ensure the existence of slow boundary spike motion for (3.1) we need only show that all of the eigenvalues of (3.45) satisfy Re(A) < 0. To accomplish this, we will apply Theorem 2.2. Before we may apply this theorem, we note that the boundary conditions, (3.45b) and (2.65) don't match. However, the eigenfunctions we are concerned with are radially symmetric and thus will satisfy (3.45b) automatically. Thus, for any parameter set satisfying (2.67), we may conclude that the real part of all of the eigenvalues satisfy Re(A) < 0. For any parameter set not satisfying (2.67), we must apply numerical techniques, as in §2.1.1, to determine if the spectrum of (3.45) contains an eigenvalue with positive real part. 67 Finally, we note that the problem (3.45) preserves the translation eigenvalues associated with the local problem (3.41), since by symmetry the non-local term I(<j>) satisfies I ( ^ . ^ c ) = 0 for j = 2, . . . , N. As is similar to (3.44), these translation eigenvalues are perturbed by o(l) as e —>• 0 when the non-local eigenvalue problem is defined on a curved boundary. The resulting small eigenvalues are responsible for the slow boundary spike motion derived in §3.2 and §3.3. 3.4 A Spike on a Flat Boundary in Two Dimensions In §3.2 we showed that the motion of a spike centered on the boundary of a two-dimensional domain is in the direction of increasing curvature. This leads us to the problem of determining the motion of a spike when the curvature is constant. In particular, we will analyze the motion of a spike on a flat boundary where the curvature vanishes. Our analysis below shows that this motion is metastable. To obtain this result, in §3.5.1 we show that the principal eigenvalue associated with the linearization of a spike solution centered on a flat boundary is exponentially small. This establishes the metastability. Then, in §3.5.2 we derive an asymptotic differential equation for the metastable spike motion on the flat boundary by imposing a limiting solvability condition on the solution to the linearized problem. This condition ensures that this linearized solution is orthogonal to the eigenfunction associated with the exponentially small eigenvalue. For the analysis we let x = (a;, y) and we suppose that the spike is located on the straight-line boundary segment joining the points ( X L , 0 ) and {XR,0) as shown in Fig. 3.6. The flat portion of dQ, is taken to be the straight-line segment between (XL, 0) and (XR, 0). The spike is centered at xo = (£,0) where xi < £ < XR. We decompose dQ, as dCl = dQc U dCls where d£ls refers to the straight-line segment of the boundary and dVLc denotes the remaining curved part of dQ,. The distance between the spike and dClc is assumed to be a minimum at either of the two corners (xi,0) or (XR, 0). The local behavior of the boundary near the corner points is critical to our analysis. Near the 68 (xL,0) (£,0) dSls (xR,0) Figure 3.6: Plot of a two-dimensional domain ft with a flat boundary segment. The spike is centered at x = £ on the flat segment. The dotted line indicates an approximate equipotential for u. corner points, dftc is assumed to have the local behavior near (xL,0); y = ipL(x), ip'L{x) KL(xL - x)aL , as x -> s£ , (3.46a) near(xH ,0); y = ipR{x), tp'R{x) ~ KR(x - xR)aR, asx->xR~, (3.46b) where az, > 0 and aR > 0. When Of, = a;? = 1, KL and . S T R are proportional to the curvature of dftc at the left and right corners, respectively. The spike solution to (3.1) is given asymptotically by a{x,t) ~ ae = hlttP-Vuc (e-^x - x0(*)|) , (3.47a) ~ e=\j4h~\Jo [ U c { ] J ' ( 3 ' 4 7 B ) where xn(t) = (£(*), 0) is to be determined. Here uc(p) is the radially symmetric solution defined in (3.14), and |ft| denotes the area of ft. We first linearize (3.1) around ae by writing a = ae + v, where v < ae. We get that u satisfies, £ e u = e2Av + (—1 + p«£ _ 1 )t; — —§T~~7^ \ dx = + dtae, in ft, (3.48a) •Kfl{s + 1) JQ dnv = -dnae on dft. (3.48b) 69 Here 0 is defined by /•CO 0= W{P)]mPdP (3.49) Jo Now we let v = eA t0 to get the eigenvalue problem, C€(j) — \<j>, in ft ; 9 n 0 = 0 on <9ft. (3.50) 3.4.1 The Translation Eigenvalue Suppose for the moment that ft is the half-space y > 0 so that 9ft3 is the entire x-axis. Then, the function 0 = dxuc satisfies £ e 0 = 0 and the normal derivative boundary condition dn4> = 0. Hence, for this case, 0 is an eigenfunction of C with a zero eigenvalue. This corresponds to translation invariance in the x direction. For our geometry, 0 is localized near (£, 0) on the flat segment dCls and 0 decays exponentially away from this point. The interaction of the exponentially small far-field behavior of 0 with the corner regions, where <9fts and <9ftc meet, perturbs the zero eigenvalue by exponentially small terms. This shift in the zero eigenvalue is calculated below. The non-local term in the operator C€ is asymptotically negligible in the calculation of this shift. However, as shown in §3.4, the non-local term is essential for ensuring that the translation eigenvalue is the principal eigenvalue of the linearization. Now we calculate Ai and 4>\. Since 0 fails to satisfy the boundary condition on 9ftc, the principal eigenfunction <f>\ has the form Here C\ is a normalization constant and <f>i is a boundary layer correction term localized near 9ftc. Let 77 < 0 be the distance between x 6 ft and 9ft and let 77 = e~xr) be the localized coordinate. Then, from (3.50), we get that <PL satisfies 0 i ~ C\ {dxuc -I- 0 L ) . (3.51) 7]<0, (3.52a) djj0L = -e9jj[dsuc]|jj=o , on 77 = 0; <pi —> 0, as 77 —>• — 0 0 . (3.52b) The solution to (3.52) is <f>L = -e (drj[dxuc]\ij=o) enl' (3.53) 70 Since uc is localized near x = £ € 9fts, we can calculate fa and fa on 9ftc by using the far-field behavior of uc given in (3.14b). In this way, we get an estimate for fa on 9ftc from (3.51) and (3.53) fa ~ -Cxace-ll2r-zl2{x - £)e~r/e (1 + r • fi) , on 9ft c. (3.54) Here ac is defined in (3.14b), r = \x — xo|, r = (x — xo)/r, and n is the unit outward normal to 9ftc. Applying Green's identity to fa and dxuc, and using the facts that dnfa = 0 on 9ft and dn[dxuc\ — 0 on 9ft s, we obtain Ai {dxuc, fa) = -e 2 / fadn[dxuc] dS + (£* [dxuc] ,fa). (3.55) Here (/, g) =. J"fi /<?<ix and C* is the adjoint operator defined by = e 2 A0 - $ + t ig-V ~ T," 2^ 1 / <^ dx- (3-56) /j7r(s +1) y n We now estimate each term in (3.55). Using polar coordinates, we calculate {fa,dxuc) ~ —^—, where 7 = y p [MP)J DP- (3-57) Next, we use (3.54) and the far-field form of uc in (3.14b) to get fadn [dxuc] ~ - C l Q r ' £ 2 2 e _ 2 r A ? • n (1 + r • n) , on 9ft c. (3.58) Substituting (3.58) into the boundary integral in (3.55), we observe that the dominant contri-bution to this integral arises from the corner regions of 9ftc, where r is the smallest. Near the corner regions we use the local behavior (3.46) to calculate { KL(XL— x)aL , as x -» x7 L (3.59) KR(x - xR)aR , as x ->• x\. Substituting (3.58) and (3.59) into the boundary integral in (3.55), and using Laplace's method, we get B = -e 2 f fadn(dxuc] dS ~ dal H K l { * l ~ ^ e- 2(?-*)A dx, •/anc 7 -oo ? - X L t ^ f * ^ ! t - M / ! f c (3.60) 71 The integrals in (3.60) are evaluated explicitly by using / e \ O H - i (5) r(<, + i), (3.61) where T(z) is the Gamma function. In this way. (3.60) becomes B ~ Cia\ (3.62) Finally, in Appendix B . l we give asymptotic estimates to show that (£* [dxuc], 0i) = o(B), as e —>• 0. Hence, we can neglect the last term on the right side of (3.55). Substituting (3.57) and (3.62) into (3.55), we get the following key asymptotic formula for the principal eigenvalue of (3.50): Proposition 3,3 (Eigenvalue) Assume that the distance between £ and dDc is a minimum at either of the two corners (XL,0) or (XR,0). Then, for e -> 0, the principal eigenvalue \\ of (3.50) has the asymptotic estimate Here ac is given in (3.14b), 7 is defined in (3.57), and Ki, KR, OIL and OLR are defined in (3.46) in terms of the local behavior of dQ.c near the corners. 3.4.2 The Slow Spike Mot ion We now derive a differential equation for £(£) for the time-dependent problem. We assume that the spike is initially on <90s. Then, since the spike motion is metastable we have vt < dtae in (3.48a). Multiplying (3.48a) by fa, and using dnfa = 0 on dQ, we obtain upon integration by parts that -2{xR-t.)/e (3.63) (3.64) 72 From (3.48) we have £ev ~ dtae in ft, and dnv = —dnae on 80,, where dnae = 0 on 8QS. Thus, (3.64) reduces to (dtucfa) ~ -e2 /" ^ ^ ^ 5 + V / ( P _ 1 ) • (3-65) Similar estimates to those given in Appendix B . l , which we omit, shows that the last term on the right-hand side of (3.65) is negligible as compared to the boundary integral term. Hence, (dtue,<(>i) e2 / 4>xdnucdS, (3.66) where uc = uc [e—1 jai — xo|] and xo{t) = (£(£), 0). . The remaining part of the analysis is very similar to the derivation of the eigenvalue estimate for Ai , and hence we omit many of the details. For e <C 1, we calculate using (3.51) that (dtuM--??^1, (3.67) where 7 was defined in (3.57). Next, using the far-field behavior of uc given in (3.14b), the estimate for <j>\ on 8£lc given in (3.54), and Laplace's method, the boundary integral in (3.66) can be evaluated asymptotically as in (3.58)—(3.62). The following main result is obtained from this calculation: Proposition 3.4 (Spike Motion) Assume that the distance between £ and 8Q,C is a minimum at either of the two corners ( Z L , 0 ) or (XR,0). Then, for e —• 0, the x-coordinate of the center of the spike along the flat segment dQ,s, denoted by £(£), satisfies the asymptotic differential equation, ^ ( ( ) . M ( ^ ( i ) a « + 1 r ( Q S + 1 ) e - 2 ( I R - f l A 7T7 lxR — ? Here ac is given in (3.14b), 7 is defined in (3.57), and Ki, KR, CLL and OLR are defined in (3.46) in terms of the local behavior of dQ,c near the corners. A similar differential equation for the motion of a straight-line interface in a constant width neck region of a dumbbell-shaped domain has been derived in [25] for the Allen-Cahn equation. Using 73 the boundary value problem solver COLSYS [3] we can solve (3.14) numerically to determine the constants ac and 7 in (3.63) and (3.68). In this way, when p = 2 we compute that ac = 10.80 and 7 = 2.47. The result (3.68) shows that the motion of the spike along the straight-line boundary segment between (XL,0) and (XR,0) is determined by the shape of the boundary at (XL,0) and (XR,0) and by the distance between the spike and the corner regions. The spike will move according to (3.68) until a stable steady state is reached or until the spike touches (XL,0) or (xR,0). Once the spike reaches the curved part of the boundary dQ.c, it will subsequently evolve according to (3.20b). From (3.68), the steady-state spike-layer location £ e on 8QS satisfies & - S L , 4 g . / « _ KLT{aL + 1) re\OiL-aR 2 { x r + X l ) / , xR-^e ~ KRY{aR + l) \2) • [6-™> Since the left hand side of (3.69) increases from 0 to co as £ e ranges from XL to XR, a unique steady-state solution to (3.69) exists on XL < £ < XR whenever KL and KR have the same sign. This solution is stable when KL < 0 and KR < 0, and is unstable when K1 > 0 and KR > 0. In particular, this implies that if 0 is convex near (XL,0) and (£#,0), then there is no stable equilibrium spike location on dCls. A simple calculation using (3.69) shows that the equilibrium spike-layer location £ e, when it exists, has the expansion <. xL + xR e te j + 4 g (§) + • • • . (3.70) Thus, the equilibrium location, £ e is located at an O(e) distance from the midpoint of the straight-line boundary segment. The following dynamical behavior can be deduced from (3.68) and (3.70) when the initial condition is £(0) = £o- When KL > 0 and KR > 0, f (i) will move monotonically towards XL if £0 < £ e, or monotonically towards XR if £0 > £e (see Fig. 3.7). When KL < 0 and KR < 0, £(i) will approach the stable steady-state at £ e (see Fig. 3.8). If KL < 0 and KR > 0, then £(i) will move towards xR (see Fig. 3.9). Similarly, £(t) will move towards XL if KL > 0 and KR < 0. When the spike touches ( X L , 0 ) or (XR, 0), its subsequent evolution is determined by (3.20b). 74 Figure 3.7: Plot of part of a domain boundary, dCl, upon which the center of the spike is at an unstable steady state. KL > 0, KR > 0 for this domain. Figure 3.8: Plot of part of a domain boundary, dQ,, upon which the center of the spike is at a stable steady state. KL < 0, KR < 0 for this domain. Figure 3.9: Plot of part of a domain boundary, d£t, upon which the center of the spike moves towards the right. KL < 0 and KR > 0 for this domain. 75 Chapter 4 Stability of n-spike equilibrium solutions to In this chapter, we study the stability of n-spike equilibrium solutions to (1.17). As will be demonstrated below, the analysis of (1.17) is considerably more involved than the previous anal-ysis of (1.18). The motivation for carrying out this difficult procedure are certain discrepancies between the behaviour of (1.18) found in §2 and numerical simulations of the full system (1.17) which (1.18) is supposed to mimic. In §2.2 we found that n-spike solutions to (1.18), where n > 1 and the spikes are all strictly within the interior of the domain, are unstable with an 0(1) positive eigenvalue. Solutions with one interior spike are also unstable but with an expo-nentially small principle eigenvalue as e —>• 0. However, numerical computations, such at those in [16], [32] and [33] suggest that equilibrium solutions with n > 1 stable interior spike solutions to (1.19) may be possible. We conjecture that for sufficiently small values of D, solutions to (1.19) with n > 1 interior spikes may be stable. The goal of this chapter is to investigate this conjecture analytically in the simple case of a one-dimensional spatial domain for equilibrium solutions with spikes of equal height. Equilibrium solutions with asymmetric equilibrium spike solutions are studied in [50] and the case of spikes in a two dimensional domain are studied in [30]. It is important to mention that our stability analysis is very different from the classical Turing-type stability analysis that is based on linearizing a reaction-diffusion system around a spatially homogeneous steady-state equilibrium solution. Our analysis is based on the study of the linearization of (.1.17) around an n-spike equilibrium solution, which has a very high degree of spatial inhomogeneity. A similar analysis for the Fitzhugh-Nagumo model has been carried out in [38]. Some stability results for the case of one spike with T ^ 0 is given in [37]. 76 We now give an outline of the chapter and summarize some of the key results obtained. In §4.1 we use the method of matched asymptotic expansions to construct equilibrium solutions to (1.17) in the limit e —> 0 that have n > 1 spikes of equal amplitude in the activator concentration. In §4.2 and §4.3 we study the spectrum of the eigenvalue problem associated with linearizing (1.17) around the equilibrium solution constructed in §4.1. In §4.2 we study the large eigenvalues of order A = 0(1) in the spectrum, while in §4.3 we study the small eigenvalues of order A = 0(e2). The n-spike solution is stable when both sets of eigenvalues lie in the left half-plane. For n > 2 and e -¥ 0, in §4.2 we obtain an explicit critical value Dn such that the large 0(1) eigenvalues are in the left half-plane only when D < Dn. When this condition on D is satisfied, we say that the equilibrium solution is stable with respect to the 0(1) eigenvalues. In §4.3, for n > 2 and e -> 0, we show that the small eigenvalues are always real and that they are negative only when D < D*. An explicit formula for D* is derived and it is found that D* < Dn. Thus, for n > 2 and e —> 0, an n-spike symmetric equilibrium spike pattern is stable when D < D*n and is unstable otherwise. The results for Dn and £>* are given below in Propositions 4.7 and 4.11, respectively. The main stability results, summarized in propositions 4.5, 4.7, 4.8, 4.10, and 4.11, are obtained from a careful but formal asymptotic analysis. It would be of interest to establish these results rigorously. Finally, in §4.4 we study the stability and dynamics of a solution to (1.17) with exactly one spike. For a certain range of exponents (p, q, r, s), we show that a one-spike equilibrium solution to (1.17) will be stable when D < Di(e), where D\(e) is exponentially large as e ->• 0. It is unstable when D > D\(e). An asymptotic formula for D\(e) is given in Proposition 4.13 of §4.4.2. This result is consistent with the result of [18] for the shadow problem (1.18) for which D = oo, where it was shown that a one-spike equilibrium solution is unstable, but with an asymptotically exponentially small positive eigenvalue. In §4.4.1, we study the dynamics of a one-spike solution to (1.17) for finite D by deriving an asymptotic differential equation for the trajectory of the center of the spike using the method of matched asymptotic expansions. The asymptotic differential equation is given in Proposition 4.12 of §4.4.1 and is favorably compared in §4.4.1 with results from full numerical computations. 77 4.1 An Asymptotic Analysis of the Equilibrium Solution For e ->• 0, we construct an n-spike equilibrium solution to (1.17) with equal amplitude using the method of matched asymptotic expansions. A solution with three spikes is shown in Fig. 4.1. The locations Xj, for j = 0,... , n — 1, of the spikes for an n-spike solution, which follows from symmetry considerations under the Neumann conditions (1.17c), satisfy XJ = -1 + ^-^-, i = 0 , l , . . . , n - l . (.4.1) n At these points the equilibrium solution satisfies a'(xj) = 0 and h(xj) = H, where H is independent of j . For an n-peak equilibrium solution to (1.17), the activator concentration is localized in the inner regions defined near each Xj, and is exponentially small in the outer regions defined away from the spike locations. In the inner region near the spike we introduce new variables by yj = e~l(x-Xj), h(yj) = h(xj + ey), d(yj) = a(xj + ey), (4.2a) and we expand h{yj) = ~h0{yj) + ehl(yj) + --- , d(yj) = d0(yj) + 0(e). (4.2b) Substituting (4.2) into the equilibrium problem for (1.17), and collecting powers of e, we get d0 — do + ag/^o = 01 — °° < Vj < ° ° , (4.3a) />o=0, (4.3b) Dhl = -dl/hl. (4.3c) The conditions at yj = 0 are that o'0(0) = 0, ho(0) = H, and hi(0) = 0. The conditions needed to match to the outer solution are that ho is bounded as \yj\ —>• oo and ao -> 0 as \yj\ -> oo. Thus, the solution to (4.3b) is ho = H. Next, we introduce uc by ao - H^Uc, where .7 = q/(p — 1). (4.4) 78 Then, (4.3a) and (4.3c) become uc — u c + u£ = 0, —co < yj < oo, (4.5a) uc ->• 0 as \yj\ -> co; u'c(0) = 0, (4.5b) DK[ = -urcH^~s. (4.5c) From phase-plane considerations, there is a unique positive solution to (4.5). In particular, when p = 2 we have tiefo) = j^ sech2 (I) . (4.6) Upon integrating (4.5c) from yj = —oo to yj = oo we obtain lim /ii — lim h\ = ——H1 T~ sbT , where br = \uc(y)]r dy . (4.7) This equation yields a jump condition for the outer solution. In the outer region, defined away from 0(e) regions near each Xj, a is exponentially small and h is expanded as h(x) = h0(x) + o(e). (4.8) Here ho satisfies Dh'Q' — pho = 0 on the interval [—1,1] with suitable jump conditions imposed across the Xj. Upon matching to the inner solution constructed above, we obtain that ho is continuous across each Xj and that the jump in h'0 is given by the right-hand side of (4.7). Therefore, ho satisfies n - l Dh'0' - ph0 = -IT^-'br ^2^(x-xk), - 1 < x < 1, (4.9a) fc=0 h'0(±l)=0, (4.9b) where 8(y) is the Dirac delta function. To solve (4.9) we introduce the Green's function G(x; Xk) satisfying DGXX - pG = -6(x - xk), - K x < l , (4.10a) • G x (± l ;a*) =0. (4.10b) 79 A simple calculation gives, Afecosh [0(1 + x)] I cosh [0(1 + xk)], -1 < x < xk , G{x;xk) Here (4.11a) [ Afecosh [0(1 - x)] I cosh [0(1 - xk)}, xk < x < 1. Ak = --L= (tanh [0(1 - xk)} + tanh [0(1 + xk)})-1 , 0 = (/i/D) 1 / 2 . (4.11b) In terms of G(x;xk), the solution to (4.9) is n - l /i0(x) = f T r - V j ] G(x; xfc). (4.12) k=0 Finally, to determine H we set ho(xj) = H and use the fact that Y^l=o G(XJ\ xk) is independent of j when the locations satisfy (4.1). This can be shown directly either by using (4.11) and summing certain geometric series, or by using the matrix analysis given following Proposition 4 below. In either way, we get 1 n — 1 H7T-(s+i) = — ^ w h e r e ag = Y^ G{xj-xk). (4.13) k=0 This leads to the following equilibrium result: Proposition 4.1 For e —> 0, an n-spike equilibrium solution to (1.17), which we label by ae(x) and he(x), is given asymptotically by n - l ae(x) ~H~<J2uc b'Hx ~ *k)] , (4.14a) fc=0 M S ) ~ - £ < 7 ( X ; S A ) . (4.14b) Here uc(y) is the positive solution to (4-5), H and ag are defined in (4-13), G is given in (4-H), and xk satisfies (4-1)-The three-spike equilibrium solution plotted in Fig. 4.1 is obtained from (4.14). To determine the stability properties of the equilibrium solution we introduce the perturbation a(x, t) = ae(x) + ext(p(x), h(x, t) = he(x) + extn(x), (4.15) 80 Figure 4.1: Plot of the activator and inhibitor concentration for a three-spike asymptotic sym-metric equilibrium solution with e = .02, D = .10, \x = 1, and (p, q, r, s) = (2,1,2,0). The solid curve is the activator concentration and the dotted curve is the inhibitor concentration. where 77 <§C 1 and <f> <C 1. Substituting (4.15) into (1.17) and linearizing, we obtain the eigenvalue problem p—i p e 2 0 M - 0 + ^ - 0 - = A0 , - l < x < l , (4.16a) T — 1 T Drixx- M^ = -e _ 1T%r0 + € _ 1 s -^ T T7, - l < x < l , (4.16b) 0 X (±1) =T?x(±1) = 0. (4.16c) In §3 we analyze the spectrum of (4.16) corresponding to those eigenvalues that are bounded away from zero as e -> 0. These eigenvalues are referred to as the large eigenvalues. In §4 we analyze the spectrum of (4.16) corresponding to those eigenvalues that approach zero as e -> 0. These eigenvalues, referred to as the small eigenvalues, are shown to be 0(e2) as e ->• 0. The goal is to determine the range of D as a function of n for which the large and the small eigenvalues both have negative real parts. Qualitatively, the small eigenvalues arise from the near translation invariance property of the system. When D = co, then he and 77 are constants in (4.16a). In this special case, the resulting 81 eigenvalue problem has n exponentially small eigenvalues. These exponentially small eigenvalues arise as a consequence of the near translation invariance property and an exponentially weak interaction between adjacent spikes (mediated by their tail behavior) and between the spikes and the boundary. The corresponding eigenfunction is, to within exponentially small terms, a linear combination of the first spatial derivative of uc \e~l(x — Xj)\. However, when D is finite so that 77 is a slowly varying function of x near each spike, then these exponentially small eigenvalues are dominated by an algebraically small spike interaction mediated by the function f](x). The leading term in the eigenfunction is still a linear combination of the first spatial derivative of uc, but the expansion of the eigenfunction proceeds in powers of e. When D = co and 77 = 0, the operator in (4.16) has exactly one positive eigenvalue in the vicinity of each spike, and this eigenfunction is of one sign. Hence, when D = 0 0 and 77 = 0, an n-spike solution is unstable on an 0(1) time scale. However, when D is decreased from infinity, the 0(1) positive eigenvalue near each spike can be pushed into the left-half plane owing to the dependence of 77 on D. This is the origin of the large 0(1) eigenvalues. 4.2 Analysis of the Large Eigenvalues In this section we analyze the eigenvalues of (4.16) that do not approach zero as e 0. In §4.2.1 we consider the case where s = 0 and in §4.2.2 we extend the analysis to treat 5 > 0. For ease of notation, the subscripts such as T J ^ shall indicate derivatives with respect to x, whereas the primes will generally refer to differentiation with respect to the stretched variable y. 4.2.1 Analysis for s = 0 To study the eigenvalue problem (4.16) it is convenient to introduce scaled variables defined by ae = H1u, he = Hv, tf> = ff7<£, 77 = Hf7, (4.17) where 7 = q/(p — 1). From (4.14a), we conclude that u ~ IJLO U C {e~lix ~ xk)] • Substituting (4.17) into (4.16) with s = 0, using (4.13) for H17"1, and dropping the overbar notation, we 82 get e <Pxx — <P + p-i qW VQ *~^HV = >«I>, -1 < X < 1, &(±1) =jfc(±l) = 0 . r - l ebrag •Kx <1, (4.18a) (4.18b) (4.18c) Using the symmetry of the equilibrium solution and the localization of the coefficients in (4.18), we look for an eigenfunction for (4.18) in the form n - l (4.19) 4>(x) ~ 53 CFE$ [€ 1 ( x ~ x k ) ] , jfc=0 for some ck, where $(y) 0 as \y\ —> co. The right-hand side of (4.18b) behaves like a sum of delta functions as e —> 0. Thus, for e —> 0, we calculate that ?y satisfies Dvxx-m = -7—/ [My)} ${y)dy22Ck6(x~Xk">> -i<x<i, ( 4 . 2 0 a ) 7fe(±l) = 0 . (4.20b) The solution to (4.20) is written in terms of the Green's function G(x;xk) defined in (4.11) as v(x) = — [uc{y)}T-l§{y)dy G(x; **)<*. (4.21) a S ° r • / - ° ° fc=0 Then, we substitute (4.19) and (4.21) into (4.18a), and use the fact that v = 1 + O(e) when \x — Xj\ = 0(e). The resulting eigenvalue problem, when written in terms of the stretched variable y = e~1(x — Xj), becomes for j = 0,... , n — 1, P roo n—1 Cj (V - $ +pup-l$) - ^ ~ \ \uc(y)X'X $(y) dy V Gfo; xfc)cfc = Cj\$, -oo < y < oo, v ' ag°r J-co k = 0 (4.22) with $ —>• 0 as \y\ -» oo. This eigenvalue problem is the same for each j when co,... , c„_i are the components of the eigenvector for the matrix problem Qc = ac, y Cn-l J (4.23) 83 Here Q is the n x n symmetric matrix whose entries are the coefficients G(XJ; Xk). The eigenval-ues of Q are real. Then, using (4.7) for br, we get that (4.22) becomes the nonlocal eigenvalue problem «• - , - ^ ( ^ f ' ^ f iV)=™, - o o < , < o o , (4.24a) a9 V J - o o M y ) J dy J $ -> 0 as \y\ -> oo. (4.24b) The goal is to determine conditions on £>, a and n for which the eigenvalue Ao # 0 of (4.24) with the largest real part satisfies Re(Ao) > 0 for any eigenvalue a of the matrix problem (4.23). The outline of the rest of the analysis is as follows. First, we obtain explicit formulae for the eigenvalues ctj and the eigenvectors Cj of Q. These eigenpairs depend on the values of D, a and n. The next step is to use a key result of [52], which we restate below, that proves that the principal eigenvalue of (4.24), in the restricted subset for which A ^ 0, has a positive real part when a < ac and a negative real part when a > ac. Here ac > 0 is some specific threshold value. Hence, we conclude that there is no eigenvalue of (4.24) with a positive real part when the minimum eigenvalue cti of the matrix problem (4.23) satisfies c*i > ac. We show explicitly the range of parameter values D, /j and n for which this relation holds. We now carry out the details of this analysis. We first calculate the eigenvalues of the full symmetric matrix Q. This is readily done since Q~l is a symmetric tridiagonal matrix. To see this, in Appendix C we solve (4.18b) on each subinterval [XJ-\,XJ] and impose the following jump conditions across each x = Xj that are associated with (4.186): rr r°° fo],. = 0, [£>%] • = -u,j , ^ = ~r [uc(y)}r-l$(y) dy. (4.25) O,g0r 7—00 Here [a]j = a(xj+) — a(xj-). This procedure then leads to a linear system for n(xj),j = 0,... , n — 1 of the form Bri = {pDyl/2uj, (4.26a) 84 where the n x n tridiagonal matrix B and the n-vectors t] and u; are defined by / J £ n n r\ n \ B d f 0 / e f 0 / e 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 e / 0 / e / 0 / d V = ^ "(x n _i) j (4.26b) Here d, e and / are defined by d = coth(20/n) + tanh(0/n), e = 2 coth(20/n), / = -csch(20/n), (4.26c) where 9 = {n/D)1/2. Note that d = e + f. Thus, r/ is given by rj = B~xu (pD)~ 1 / 2 . Another way to determine rj is to evaluate (4.21) at x = Xj, for j = 0,... , n — 1. The equivalence of these two representations of r) yields ZT 1 g = fpD (4.27) In Appendix C . l we show the explicit calculation that yields the following result for the eigen-values Kj and the eigenvectors q3 of B: Proposition 4.2 The eigenvalues Kj, ordered as 0 < «i < . . . < nn, and the normalized eigenvectors qd of B are (4.28a) «i = e + 2/ ; Kj = e + 2/ cos ( ) , j — 2, . . . , n , n 1 g x = - ^ ( 1 , . . . ,1) ; qij = \j-cos <j ~ 1) n ( J - l / 2 ) , j = 2,. . . , n . (4.28b) Here q* denotes transpose and q^ = (51 j,... , g n j - ) . Therefore, from (4.27) and since e > 0 and / < 0, the smallest eigenvalue of Q is proportional to K~1 and the corresponding eigenvector is proportional to qn. Relabeling this eigenpair we obtain: 85 Proposition 4.3 The smallest eigenvalue ct\ ofQ and the corresponding (unnormalized) eigen-vector qi are - e — 2/cos (7r/n) . fir(n-l)\ fw(n - l)(l -l/2)\ , A n n l . Here e and f are defined in (4.26c) and q\ = (q\ti,. • • ,qi,n)-We now apply Theorem 2.1 to find a criterion to ensure that the 0(1) eigenvalue has negative real part. Comparing (4.24) with (2.11), the theorem above yields the following key result on the spectrum associated with (4.24): Proposition 4.4 Let Ao 7^ 0 be the eigenvalue of (4-24) with the largest real part and assume condition (2.13a) holds. Then, i?e(Ao) > 0 when [v — l)a„ ,, a\ < ac where ac = — — . (4.30) qr Also Re(Xo) < 0 when a\ > ac. Here a.\ is the minimum eigenvalue of Q given in (4-29a) and ag is the constant row sum of Q defined in (4-13): To get an explicit stability criterion we must calculate ag. Since q\ = (1,... , 1) is an eigenvector of B with eigenvalue n\ = e + 2/ we can multiply both sides of (4.27) by q1 to get Gqx = ag (1, . • . . , ! ) * = - i ^ q i = — i = r f (!>•••> 1)* • (4-31) Hence, a" ~ yfilD ( e + 2 / ) 1 2^D coth (26/n) - csch (26jn) ' ( 4 ' 3 2 ) Substituting (4.29a) and (4.32) into (4.30) we obtain that Re(A0) = 0 when e + 2/ p-1 (4.33) e — 2/ cos (n/ri) qr Using the definition (4.26c) for e and / , we calculate e/(2f) = - cosh (26/n). Substituting this expression into (4.33) and solving for the critical value of 9 we get the following main result: 86 Proposition 4.5 Let Ao 7^ 0 be the eigenvalue of (4-24) with the largest real part and assume condition (2.13a) holds. Then, Re(\o) < 0 when D<Dn = £ , n = l , 2 , . . . , "n 0„ = £ln a + Va2 - 1 o = l + 1 + » ( - ) ] - 2 T - 1 \ n / J \p — 1 Alternatively, when D > Dn then i?e(Ao) > 0. (4.34a) (4.34b) This result gives the stability criterion for the large eigenvalues of (4.16) when s = 0. For example, from this result we can conclude that a three-spike equilibrium solution is stable with respect to the large 0(1) eigenvalues only when when D < Dz- To stabilize one additional spike we need to decrease D below D4. We now examine (4.34) for the G M parameter set (p, q, r, s) = (2,1,2,0) for which 2 + C O S l n J + V l 2 + C ° H n J J - 1 We then calculate the following sequence of critical values of Dn: D1= n/6f = oo, 0 i = O , D2 = a/9l = 0.5766/i , 92 = In (2 + v ^ ) , D3 = n/el = 0 .1810A* , D 4 = fx/ej = 0.0915/x, In the limit n 2> 1, we get 3, 5 V21 * 3 = 2 l n 2 + — Dn ~ 4 /xn" 2 (ln [3 + v^]) * + O (rT 4) . (4.35) (4.36a) (4.36b) (4.36c) 04 = 2 In 2^ + ^ + yi7~2\/2~j . (4.36d) (4.37) For the analysis leading to (4.34) to be valid we require that D/e2 » 1 in order to ensure that h is slowly varying in the inner regions near each spike. Hence it follows that we require n < 1/e, which limits the range of validity of (4.37). 87 For the other common parameter set (p, g, r, 5) = (4,2,2,0) we get the critical values (4.38a) (4.38b) (4.38c) 01 = CO , 6\ = 0 , = nlQ\ = 0.2349M , 02 D3 = 0.0778M , 03 = ^ 4 = 0.0401M , 04 3, M l V1T7 2 12 2 3 ^ , /39~ 2 (4.38d) The results in (4.36a) and (4.38a) suggest that the principle eigenvalue of (4.24) for a one-spike equilibrium solution will always have a negative real part for any value of D. This conclusion is true when D is independent of e, but needs to be modified if we allow D to depend on e. More specifically, we show in §5 that a one-spike equilibrium solution is stable only when D < D\(e), where D\(e) is exponentially large as e ->• 0 and satisfies D\(e) — O (e2e2/£) for e <C 1. 4.2.2 Analysis for s > 0 For s > 0 we again introduce the new variables (4.17) into (4.16) and use # T " - ( S + 1 ) = l/(brag) from (4.13), with the result e 2 0xx-0 + ^ - 0 - ^ r r 7 = : A 0 , -1< x < 1, (4.39a) ^ - M " - ^ = - ^ ^ - ^ 0 , - K x < l , (4.39b) C l / j * Uig U C1/7* U/g U 0 x ( ± l ) = T f e ( ± l ) = O. (4.39c) Here « ~ ]Cfc=d u c [ e - 1 ^ _ x f c ) ] - Substitute the form for 4> given in (4.19) into (4.39b) and use the facts that u is localized and that v = 1 + O(e) near each xk. Then, in place of (4.20), (4.39b) and (4.39c) become fDrix n - l S M + — Y ] S(x - xk) a„ — 9 k=0 r f°° _ i ' \—i 77 = - / [uc{y)]r §(y)dy^Tck5(x-xk), \x\<\, a9°r J-00 £ = i k=0 (4.40a) r]x(±l)=0. (4.40b) 88 Thus the term proportional to s in (4.40a) acts as a psuedo-potential and hence will modify the jump condition for rjx across each Xj. Since u is localized near each Xj, and n(x) is slowly varying with respect to e near each Xj, we need only calculate rj(xj) and substitute into (4.39a) to obtain the eigenvalue problem. To calculate n(xj) we proceed as follows. We introduce 77 and u; as defined in (4.25) and (4.26). We then solve (4.40a) analytically on each subinterval in terms of hyperbolic functions and then patch the subinterval solutions together using the appropriate jump conditions fo]- = 0 , = -UJ + —r,(Xj), (4.41) ag where Uj was defined in (4.25). This calculation, which we omit, shows that the solution for 77 can be written in the form BsV = (uD)-V2 u, (4.42) where the matrix Bs is given by "B, = B+—?=I. (4.43) Here I is the n x n identity matrix and B is the matrix defined in (4.26b) and (4.26c). Therefore, using (4.25) and (4.26b), we obtain V = /oo [uc{y)}r-l${y)dyB:lc, (4.44) •00 CLgbry/nD where c is defined in (4.23). In place of (4.22) we get, for j = 0,... , n — 1, that P poo c . ( $ " _ $ + j n i P - i $ ) _ _ ^ l _ j^JUc(y)Y-^(y)dy ( 5 7 ^ = ^ , - o o < y < c o , (4.45) with $ -> 0 as \y\ -> 00. Here (Bjlc). denotes the j * * 1 component of the vector B~lc. Now let c be an eigenvector of the matrix eigenvalue problem Bsq = Kq. , (4.46) 89 Then, using (4.7) for br, (4.45) becomes 7-°°oo My)rx <%) V qrupc agK*/jlD y [uc{y)]r dy $ -> 0 as |y| oo . A<&, -oo < y < oo, (4.47a) (4.47b) Let Ao 7^ 0 be the eigenvalue of (4.47) with the largest real part. Then, from comparing (2.11) and (4.47), we conclude from the theorem of [52] stated above that Re(Ao) < 0 only when > (P -I ) (4.48) K,agy/(j,D qr To obtain a condition in terms of the minimum eigenvalue of Q, we use (4.27) to get that Qq = aq, where Kagy/pD =" s + ag/a. Substituting this relation between K and a into (4.48), we obtain the following result in terms of the smallest eigenvalue a\ of Q: Proposition 4.6 Let Ao ^ 0 be the eigenvalue of (4-24) with the largest real part and assume condition (2.13a) holds. Then, i2e(Ao) > 0 when - l ai < qr P~^T (4.49) Also Re(Xo) < 0 when the inequality in (4-49) is reversed. The right-hand side of (4.49) is always positive by the assumption (1.8) on the exponents. Setting ct\/ag = [qr/(p — 1) — s] _ 1 , and using (4.29) and (4.32) for ct\ and ag, respectively, we get the following main result for the stability of the equilibrium solution with regards to the large 0(1) eigenvalues: Proposition 4.7 Let XQ ^ 0 be the eigenvalue of (4-4V with the largest real part and assume condition (2.13a) holds. Then, Re(Xo) < 0 when D<Dn = £ , n = l,2,. 6n = ^ln a+ Va2 - 1 a = 1 + 1 + cos qr P~^T (s + l) (4.50a) (4.50b) Alternatively, when D > Dn then Re(Xo) > 0. 90 From (1.8) we get that a > 1 since qr/(p — 1) > (s +1). In addition, Dn decreases as s increases, and so for each fixed n it follows that D must be made smaller as s increases in order to stabilize an n-spike equilibrium solution. 4.3 Analysis of the Small Eigenvalues The results in §4.2 establish conditions for which the equilibrium solution is stable on an 0(1) time scale. Now, we examine the more difficult problem of determining conditions guaranteeing that the small eigenvalues with A = 0(e2) lie in the left half-plane. The first step, done in §4.3.1, is to reduce (4.16) to the study of a matrix eigenvalue problem. In §4.3.2 we analyze this matrix eigenvalue problem to determine the small eigenvalues and their signs explicitly. 4.3.1 Deriving the Matrix Eigenvalue Problem We begin by writing (4.16) in the form p L£(t) - = Xcj), - K x < l , (4.51a) a T - l ar Dr}xx-un = -e-lr-^4> + e-1Sj-£Ir), -1 < x < 1, (4.51b) <^(±1) =77 x (± l ) =0, (4.51c) where p - i L€4> = e2(j>xx-(t) + ^^-(j>. (4.51d) Here ae and he are given by ra—1 n—l ae-'^Htuk; fte~-^G(I;Ii)i #7r-(i+.) = _ L . (4.52) We have defined Uk{y) = uc [e_1(o; — x^], where uc(y) satisfies (4.5). The equilibrium positions for XJ are such that {hex)j = 0, ; = 0 , . . , n - l . (4.53) Here and below we have defined (Qj = (C(xj+) + Cixj-))/^ an-d. = Cixj+) — C^j-), where C(XJ±) are the one-sided Hmits of £(x) as x -» Xj±. 91 If the inhibitor diffusivity was infinite and there only one spike, then by translation invariance we would obtain L € a e x = 0. Here we expect that L € a e x is still small. To show this, we differentiate the equilibrium problem for (1.17a) with respect to x to get Thus, for x near Xj we get Leaex - ~^TThex • (4.54) , eqH«uP LeUj g + 1 h e x • (4.55) This fact suggests that we expand 0 = 00 + e^ i H , v(x) = eVo(x) H , (4.56a) where n—1 n—1 00 = 5Z cJu'j [e~l(x ~ xj)] > 0i = 5Z C J ^ ' [ e _ 1 ( x ~ xi)] ' (4.56b) i=o j=o and the Cj are arbitrary coefficients. We substitute (4.56a) into (4.51a) and use (4.55) and A = 0{e2). For x near XJ, we get that <t>\j{y) satisfies quPH" CjLefaj ~ — [ c j h e x ( x j + ey) - H^rjoixj + ey)] . (4.57) Before solving this equation for <p\j we need to determine an important continuity property of the right-hand side of (4.57). Substituting (4.56a) into (4.51b), we get that 770 satisfies ar-i ar Dvoxx-Mo = -e-2r^—{(t>0 + e<t>1) + e-1s—^rjo, - K x < l . (4.58) Since 0o is a linear combination of u'j, it follows that the term multiplied by 0o on the right-hand side in (4.58) behaves like a dipole. Hence, for e < 1, this term is a linear combination of (x - Xj) for j = 0, ..,n - 1, where S(x) is the delta function. Thus, 170 will be discontinuous across x = Xj. However, if we define the function f(x) by f{x) = Wr}0{x)-CjheX{x), (4.59) 92 then / is continuous across x = Xj. To see this, we differentiate (1.17b) for he with respect to x and subtract appropriate multiples of the resulting equation and (4.58) to find that the dipole term cancels exactly. Thus, / is continuous across x = Xj, and we have {f)j = f{x3). However, (hex)j — 0 from (4.53). Hence, f(xj) = . f f 7 ( T 7 o)j . Therefore, for e <C 1, we get from (4.57) that 4>ij satisfies cjLecj>lj^qu]H''-l{r}Q)j. (4.60) Since L€Uj — {p — l)uP- + 0(e), (4.60) is easily solved to get c^ l j (y ) = ^ T n j ( y ) i f 7 - 1 ( r ? 0 ) 1 + O ( e ) . (4.61) This condition shows that 4>\j is continuous across x = Xj and has the form of a spike. This implies that the term in (4.58) proportional to <f>\ behaves like a linear combination of S(x — Xj) when e « l and, most importantly, is of the same order in e as the dipole term proportional to 0o- This shows the fact that we need to determine the approximate eigenfunction for <f> to both the 0(1) and 0(e) terms in order to calculate an eigenvalue of order 0(e 2). Next, let e -> 0 and use (4.56b) and (4.13) for ffTMs+i) t o calculate for x near x$ that "e 2r~h~(t>0 ~ _ ~ ci5 ( x ~ XJ) ' (4.62a) tig Qg T-l -e ' X r ^ - < h j ~ T ~ / < - dy 5(x - X j ) . (4.62b) rle ag0r J-oo , r - l j ; iguT J—oo Substituting (4.62) into (4.58), and using the formula (4.61) for faj, we get n - l S a9 Tri—ry n—1 n—1 T 7 0 - H 7 V - - ^ - - ^ — E c/(x - )^ - („_•>)„ Y,(*>)j8(x - Xj). a9 j = 0 \P j = Q (4.63) This problem is equivalent to Dvoxx ~ Wo = 0, - 1 < x < 1; r? 0 x(±l) = 0, (4.64a) P"l, = -(^) f l 7- ,< ' ^ = ^ (3-5?il)w- <««>> 93 For convenience we introduce 770 defined by % = H^fjo • (4.65) Next, we estimate the small eigenvalue. Substitute (4.56) and (4.65) into (4.51a) and multiply both sides of (4.51a) by u'j. Integrating the resulting equation across the domain, we get n—1 n—1 / p~ \ n—1 ^2 CiLeu^ + e (u'j, CiL^ii^ - eqH1'1 ( u'j, J ~ A ^ u'j) • (4-66) i=0 t=0 \ he / i = Q Here we have defined (/, g) = f^x f(x)g(x) dx. To within negligible exponentially small terms,, the dominant contribution in the sum comes from i = j since u'j is exponentially localized near x = Xj. Thus, (4.66) becomes Cj (u'j,L€u'j^ + ecj (u^Lefaj^) - eqHl+q ^u'j, ~ Xcj (u'j'u'j) • (4-67) Since Lf is self-adjoint, we integrate by parts on the second term on the left-hand side of (4.67) and use (4.55) for L t u y The integrands are localized near x = Xj. Thus, writing the resulting integrals in terms of the stretched variable y = e~1(x — Xj), we get -^heXdy-e2qH1^ / - ^ T ? 0 dy -oo he J —oo he /°° UP(j)-\ • f 0 0 / \ 2 -J—LheXdy~t\cj K) dy. (4.68) -00 he J— oo v ' In this expression 770 = 770(XJ + ey), he = he(xj + ey), and hex = hex(xj + ey). We now estimate each of the terms in (4.68). Since [4>\j]j = 0, (hex)j = 0, and u'j is odd, it follows that COO *,P Uj<t>lj -00 he^ , 3 \ /oo u (£>-[ • - J^ -hexdy = o(l) as e->0. (4.69) - Hence, the third integral on the left-hand side of (4.68) will be o(e3) and can be neglected. Next, we combine the first two terms on the left-hand side of (4.68) to get /°° UPU f°° u p u - f°° u'u^ ^ hex dy - e2qH l+i J ^ ^ 7 7 0 dy = -e2qH' J ^ + ey) dy. (4.70) 94 Here /(a;), defined in (4.59), is given in terms of 770 by f(x) — Hfjo(x) — Cjhex(x). The function / is continuous across x = Xj but its derivative is not. For e < 1, we calculate -e'qHi ^f(X3 + ey)dy~e*qCjhex*{Xj) f ° yu'tf dy - e'q^j ["\u)u]dy. J—00 h,e J—00 J—00 . (4.71)-Upon integrating by parts in (4.71), and using h e x x ( x j ) = pH/D, we get ^Ilwi'ir-^*~,^M-DL[uMrliv- (472) Substituting (4.69) and (4.72) into (4.68), we obtain a formula for A. We summarize the result (redefining 770 for convenience) as follows: Proposition 4.8 The eigenvalues of order 0(e2) for (4-16) satisfy /co r - 2 n r°° / r 11 \ M [«; (V)] dy~^j JucMY+Uy^-^) , j = 0 , . . , n - l . (4.73) Here {r]x}j is obtained from the solution to the boundary value problem Dr)xx-nrj = 0, -Kx<l; r)x{±l) = 0, (4.74a) [Dn]3 = ; [Dr,x]j = U(n)3, 5 S s - ^ . (4.74b) 4.3.2 Analyzing the Matr ix Eigenvalue Problem The next step in the derivation is to calculate {r)x)j from the solution to (4.74). The solution to (4.74) can be decomposed as I / n - l n - l \ v(x) = — {j2ckg{x;xk) + ^2mkG{x;xk)\, (4.75) 9 \k=0 k=0 ) for some coefficients mk, for k = 0,... , n - 1. Here G satisfies (4.10), and g(x; xk) is the dipole Green's function satisfying Dgxx - P9 = S {x - xk), - 1 < x < 1, (4.76a) gx(±l;xk) = 0. (4.76b) 95 Satisfying the jump conditions in (4.74b) we get the following matrix problem for the coefficients —g + l)m = -—VqC. ag ) ag Here Q is the Green's function matrix defined in (4.23) with entries G(xy,Xk), and \ (g(x0;xo))o •••• g(xQ;xn-i) y g(xn-i;x0) ••• (g{xn-i;xn-i))n^i J m = / m 0 ^ c = (4.77) . (4.78) The problem (4.77) determines m in terms of c. Then, using (4.75), we can calculate (r)x)j, for j = 0,... , n — 1, from the matrix problem (Vx) = — [Qgc + Vm) , (4.79) where Qg is the Green's dipole matrix defined by G9 = and V=. ^ (Gx{x0;x0))o y Gx(xn-i;xQ) gx(x0;x0) ••• gx{xo;xn-i) \ (4.80) Gx(xQ\xn-\) (Gx{xn-i;xn-\))n-i J ( (Vx)o ^ \ {Vx)n-1 ) (4.81) Next, we define a by A «V (f-ooMy)}p+l dy Substituting (4.79) and (4.82) into (4.73), we get a matrix eigenvalue problem for a and c (4.82) Qgc + Vm= (<j + ^ ) c . (4.83) 96 Here m is determined in terms of c by (4.77). The next step in the analysis is to reduce (4.77) and (4.83) to an equivalent generalized eigen-value problem. This analysis involves matrices associated with G and g. To avoid confusion we have indicated with a subscript g those matrices associated with the dipole Green's function g. In the analysis below we must find the eigenvalues of Qg explicitly. This is done as in §3 by showing that Qq 1 is a symmetric tridiagonal matrix. More specifically, in Appendix C.2 we show that 9 Q9 = ^ V ,. (4-84) where Bg is a tridiagonal matrix with exactly the same form as in (4.26b), except that here the definitions of d, e and / in (4.26c) are to replaced with d = coth(20/n) + coth(0/n), e = 2 coth(20/n), / = -csch(20/n), (4.85) where d — e — f. In Appendix C.3 we calculate the eigenvalues and eigenvectors of Bg analyti-cally. The result is summarized as follows: Proposition 4.9 The eigenvalues ordered as 0 < £i < . . . < £ n , of Bg and the associated normalized eigenvectors Vj of Bg are £j = 2coth (jpj -2csch cos • j = l , - - - , n , (4.86a) v*=^(1, -1'lj • • •'(-1)n+1); vi>j=vfsin (/~1/2))}' j = h • • •'n"1' (4.86b) Here vf denotes transpose and — (u i j , . . . ,vnj). Other key relations that we need are derived in Appendices C and C.2, where we show that 97 Here the matrix C is defined by C = V From (4.87a) we obtain the result that 1 1 0 •• • 0 0 0 -1 0 1. •• • 0 0 0 0 -1 0 ' • 0 0 0 0 0 0 ' • 0 1 0 0 0 0 •• • -1 0 1 0 0 0 •• • 0 -1 -1 VB = -(VgBgf (4.87b) (4.87c) We begin by solving (4.77) for m. The matrix in (4.77) is invertible if s (ai/ag) + 1 < 0, where a\ is the minimum eigenvalue of Q. From (4.49) and the definition of s in (4.74b), we see that this condition is satisfied when the large 0(1) eigenvalues are in the left half-plane. We will assume that D < Dn so that this condition holds. Let KJ be the normalized eigenpairs of B as given in Proposition 4.2 for j = 1,... , n. Then, B = QKQ1, (4.88)" where Q is the orthogonal matrix whose columns are the normalized qj and K is the diagonal matrix of the eigenvalues of B. Since Q = B~llyf\xD from (4.27) and Q}Q = I, we get -g + i Q / c _ 1 +1 •agSnD Using 9 = (u/D)1/2, we can solve for m in (4.77) in the form ~s6 (4.89) m = a„ -Q £ - 1 + / Q*VgC. "9 \affA* We then substitute (4.90) and (4.84) into (4.83). This yields, ( 9 )VQ (^ K - 1 + i f * Q'7>,c = ( ^ - + (4.90) B g C S D \agaj ' " \aga (Da jj,ag (4.91) 98 In (4.91) we use (4.87c) and (4.88) to replace VQ with * r \ * — l VQ = VBB~lQ = - (VgBgy QIC We then introduce u and the diagonal matrix V defined by u = Bg-lc, V = sD^IC'1 (sjIC'1 + (4.92) (4.93) Here we have defined 7 by 7 = — = 2 , . 29\ , (29 coth — — csch — = 2 tanh n (4.94) n J \ n Equation (4.94) is obtained from using the expression for ag in (4.32). Using Proposition 4.2 for the eigenvalues KJ of JC we then calculate V as V = di 0 0 '•• \ where Kj + 57 V 0 0 • • • dn ) Substituting (4.92) and (4.93) into (4.91), we obtain the eigenvalue problem j = l,...,n. (4.95) BgU = Lj{I + Ti)U. (4.96a) Here we have defined ui and 1Z by 'Da 1 9 7 TZ = (VgBgf QVQ'VgBg (4.96b) The assumption that the solution is stable with respect to the large 0(1) eigenvalues is equiv-alent to the condition that KJ + 37 < 0 for j = 1,.., n. Under this condition, and since s < —1, we conclude that V > 0. Hence, / + TZ is a positive-definite and symmetric matrix. This means that the eigenvalues Uj, and consequently Xj, for j — 1,... , n are real. The generalized eigenvalue problem (4.96) is equivalent to the combined problem (4.77) and (4.83). The next step is to determine the spectrum of (4.96) analytically. To do so, we show that TZ has the same eigenvectors as Bg. Hence, we claim that 1Z can be written in terms of some positive 99 diagonal matrix E as K = QgZQg . (4.97) Here Qg is the eigenvector matrix associated with Bg. The column of Qg is the eigenvector Vj given in Proposition 4.9. Using the formula for VgBg in (4.87a), we can write 1Z in (4.96b) as U = W> csch2 (^ ) C'QVQ'C. (4.98) This is equivalent to K = ^ c s c h 2 QgMVM^l, (4.99a) where the matrix M. is defined by M^Qi^Q. (4.99b) 5 9 Comparing (4.99a) and (4.97), we then define E by E = -4oCsch2 (—^ . M D M * . (4.100) We now show that E defined in (4.100) is a diagonal matrix. To show this, we first calculate the matrix M in (4.99b) using the explicit formulae for the eigenvectors of B and Bg given in Propositions 4.2 and 4.9. Let m^j be the i,j entry of the matrix M. Then, we calculate rriij explicitly using (4.99b) and the definition of C in (4.87b), to get n m. 1=1 Here we have defined qoj = qi j and qn+i,j = qn,j, where qij- and vi j are defined for / = 1,... , n and j = 1,... , n in (4.28b) and (4.86b), respectively. A tedious, but straightforward, calculation shows that m^j = 0 for i ^ j — 1. However, the entry W j - i j is non-zero. We calculate, for j = 2, . . . , n, that n - l rrij-ij = vij-i {qi,j - gaj) + u„ j_ i (qn-ij - qn,j) + vU-i (ll-ij ~ Ql+ij) • (4.102) 1=2 100 Using (4.28b) and (4.86b), and standard trigonometric identities, we can reduce (4.102) to (4.103) 4 . f i v ( j . 2 f i r { j - l / 2 ) \ . n m ^ = n S m ( ^ ) g S m ( n ) ' J = 2 ' ' • • ' n . Therefore, we get the key result that rrij-ij = 2 sin n j = 2,. . . n ; ?n j j = 0 otherwise. (4.104) Therefore, it is clear that the matrix product MVM1 in (4.100) is a diagonal matrix. This shows that Bg and 1Z have the same eigenvectors. Then, by using (4.95) for the diagonal entries of V, we calculate S in (4.100) explicitly as s = zi 0 o '•• 0 0 0 0 (4.105a) where 1 / 20\ = ID 2 " c s c h 2 (TT J ^ m j ' J ' + 1 ^ 2 d j + 1 ' i = • • • - n - 1 ; «h = o. (4.105b) Finally, we use (4.104) for m,jj+i, (4.95) for dj+i and the result that KJ+I = £j for j = 1,... , n — 1, as obtained by comparing (4.28a) and (4.86a). In this way, we find that Zj = Zj(s), where Zj = r-csch — sm — 6 + \ n / \ n j = 1,... , n - 1; 0. (4.106) Here 7 was defined in (4.94). When D < Dn, so that the solution is stable with respect to the large 0(1) eigenvalues, then + 37 < 0 for j = 1,.., n — 1. Since we have shown the crucial result that Bg and TZ have the same eigenvectors, it is easy to calculate the spectrum of (4.96). The eigenvalues CJJ of (4.96) are <"j = £i/(l + *j), i = l , . . . , n , (4.107)" where £j and ZJ are given in (4.86a) and (4.106), respectively. Then, substituting (4.106) into the expression for a in (4.96b), we get that Uj = o-j(s), where 0 fij - 1 •3 ' j = 1,... , n. (4.108) 101 Finally, we substitute (4.108) into (4.82) to obtain explicit formulae for the small eigenvalues A = 0(e2). The main result is summarized as follows: Proposition 4.10 For e « 1, consider the eigenvalues of (4-16) of order A = 0(e2). The corresponding eigenfunction has the form (4-56) where Cj = Vj, with Vj defined in (4-86b). The explicit formula for the small eigenvalues is X e V ( f-oo iuc(y)]P+1 dy\ n - cos (nj/n) - Z j (cosh (26/n) - 1)' j ~ D{p + 1) y [u'c(y)f dy J [ cosh (26/n) - cos (irj/n) for j = 1,... , n. Here Zj = Zj(s) is defined in (4-106). (4.109) The final step in the analysis is to determine the sign of aj with respect to the parameter D. The condition aj < 0 for j = 1,.., n holds when ^ - l - ^ > 0 , ; = l,.. ,n. (4.110) Defining Wj = we calculate from (4.86a) and (4.94), and from some standard trigonometric identities, that = 1 W (g) csch* (£) . (4.111) Since zn = 0 and wn > 1, (4.110) holds when j = n. Substituting (4.106) and (4.111) into (4.110), we see that aj < 0 when sin2 (g) csch2 ( f + l) > sin2 ( ^ ) csch2 , for j = 1, ..,n - 1. (4.112) Using (4.111) and some standard identities, (4.112) reduces to ( 1 + .- + - d , » (£))(!-«• (£))<o. (4.113) The second bracketed term on the left-hand side of (4.113) is always positive for any j = 1,.., n. The first bracketed term is negative when D is very small since s < — 1 and However, this term will switch sign when D crosses through the critical value where csch2 (^j = -(1 + 5). (4.114) 102 Hence, n — 1 of the small eigenvalues switch sign at the same value of D. Let D*n be the value of D satisfying (4.114). By solving (4.114) we obtain the following main result for the stability of the solution with respect to the small eigenvalues: Proposition 4.11 For e « l , consider the eigenvalues of (4-16) of order A = 0(e 2). These eigenvalues are negative only when D < D*n, where (4.115) There are n — 1 small positive eigenvalues when D > D*n. When D = D*, then A = 0 is an eigenvalue of algebraic multiplicity n — 1. It is a simple exercise to show that, in general, these critical values are smaller than the critical values Dn given in Proposition 4.7 for the stability of the solution with respect to the large 0(1) eigenvalues. Thus, our final conclusion is that an n-spike equilibrium solution will be stable only when D < D*. For the parameter sets (p, q, r, s) = (2,1,2,0) we get 0 = 1, and for (p,q,r,s) = (4,2,2,0) we get 0 = 3. From (4.115) we then calculate the critical values n = 2 -> D2 = 0.3218/u 0 = 1; D2 = 0.1441/j 0 = 3, (4.116a) n = 3 -»• D3 = 0.1430/i 0 = 1; £>3 = 0.0641/j /? = 3, (4.116b) n = 4 -> Di = 0.0805/i 0 = 1; D4 = 0.0361/i 0 = 3. (4.116c) The numerical computations of [19] of the time-dependent problem (1.17) with (p,q,r,s) = (2,1,2,0), starting with initial conditions close to an asymptotic equilibrium solution, suggested that D2 « 0.33 and D% w 0.14. The detailed analysis presented above gives the theoretical basis for these numerical predictions. 4.4 The Dynamics of a One-Spike Solution In this section we analyze the dynamics of a one-spike solution to (1.17). For finite inhibitor diffusivity D, in §4.4.1 we derive a differential equation determining the location xo(t) of the maximum of the activator concentration for a one-spike solution to (1.17). By linearizing this 103 D* - A* \nln(J0 + J0TT)] 2 ' 0 qr J^l (1 + *) differential equation around the stable equilibrium location XQ = 0, we show that the decay rate of infinitesimal perturbations coincides with the small eigenvalue result (4.109) when n = 1. Alternatively, when D = oo, we know from [18] that the equilibrium solution XQ = 0 for a one-spike solution is unstable. When D = oo, the spike drifts exponentially slowly towards the closest endpoint of the domain (cf. [18]). To reconcile the finite D result with the infinite D analysis of [18], we show in §4.4.2 that the equilibrium location XQ = 0 for a one-spike solution is stable when D < -Di(e), where D\ is exponentially large as e —¥ 0. 4.4.1 The Differential Equation for the Spike Location In the inner region near the spike we introduce the new variables y = e _ 1 [x - X0(T)] , h{y) = h{x0 + ey), a(y) = a(x0 + ey), r = e2t, (4.117a) and we expand h(y) = ho(y) + eh^y) + • • • , a(y) = aQ(y) + ea^y) + • • • . (4.117b) Substituting (4.117) into (1.17), we find from the leading terms that ao and ho satisfy (4.3a) and (4.3b), respectively. Hence, oo(y) = Wuc(y), ho{y) = H, 7 = « / ( p - l ) , (4-118) where uc(y) satisfies (4.5). Here H = H(T) is a function to be determined. Setting a\ = H^ui, we get to next order that " n - l QUchi II .. . . ux — u\+pvfc u\ = —— x0uc, —oo<y<oo, (4.119a) H D~hl = -WT-suTc, (4.119b) with ui ->• 0 exponentially as \y\ —> co. The right-hand side of (4.119a) must be orthogonal to the homogeneous solution u'c of (4.119a). From this solvability condition we obtain x'o = —~ 2 / ucu'Mdy. (4.120) H LooKiy)] dyJ-oo 104 If we integrate (4.120) by parts twice, and use the facts that h'[ and uc are even functions, we get -Q (TooM^r1 dy\ Xn — 2H(p+l) \ [n'M? dy lim hx + lim h, y—t+oo y—¥—oo (4.121) In the outer region away from the spike, a is exponentially small and, similar to the analysis in §2, we expand h = ho + • • •, where ho satisfies Dh0 - pho = -H1T $br6(x - x0), - K x < l , h'0(±l) = 0. Here br is defined in (4.7). To match with the inner solution we require that h0{x0) = H, lim hx+ Um hx = hQx(xo+) + h0x(xo-) 3/->+oo y-t—oo (4.122a) (4.122b) (4.123) The solution to (4.122) is h0(x) =Hir-sbTG{x;xo), (4.124) where the Green's function G(x; xo) satisfies (4.10). Substituting (4.124) into (4.123), and using (4.13) for IP*-', we get H lim hi + lim hY = y->+oo y-+-co G(xo;xo) [GX{XO+;XQ) + Gx(x0-;xo)] , H = 1 l/7r-(a+l) (4.125a) (4.125b) _bRG(xo;xo)_ The solution G{X;XQ) was given in (4.11). Using this solution we can calculate the right-hand side of (4.125). Then, substituting the result into (4.121) and (4.124), we obtain Proposition 4.12 For e « 1 , the dynamics of a one-spike solution to (1.17) is characterized by a(x,t) ~ Wuc (e 1[x-x0(t)]) , h{x, t) ~ EG [x; x0{t)] /G [xQ(t); x0(t)\, (4.126a) (4.126b) 105 where H = H(t) is given in (4.125b). The spike location xo(t) satisfies the differential equation where C is defined by tanh( A / g ( l + xo) tanh [ - x0) C = 1* ( f-oo iuc(y)]P+l dy 2(p + l ) V D V S^K{y)?dy (4.126c) (4.126d) The equilibrium solution xo = 0 for this differential equation is stable for any D. The decay rate of infinitesimal perturbations around xo = 0 is F (0) ~ — eV {£»M)rX *y\ ^ 2 ( jp, D(p + V\ IZoK(y)fdy J~ VV This result agrees precisely with the small eigenvalue result (4.109) when n = j = 1. (4.127) log 1 0 (i + 1) Figure 4.2: Plot of the trajectory xo{t) of the center of the spike for a one-spike solution with e = .03, /j = 1.0, D = 1.0 and (p,q,r,s) = (2,1,2,0). The solid curve is the full numerical result and the dotted curve is the asymptotic result. To verify (4.126c) for the parameter set (p,q,r,s) = (2,1,2,0) we compared the asymptotic result (4.126c) for Xo(t) with the corresponding full numerical result computed from (1.17). The problem (1.17) was solved numerically using the routine D03PCF from the NAG library [39]. The initial condition was taken to be of the form (4.126a) and (4.126b) with x0{0) = 0.6 i06 0.4 h -1.0 -0.5 0.0 0.5 1.0 Figure 4.3: Plot of the initial condition for a one-spike solution corresponding to the parameter values shown in the caption of Fig. 4.2. The solid curve is the activator concentration and the dotted curve is the inhibitor concentration. t log 1 0 (l + <) xo{t) (num.) xQ(t) (asy.) 12.0 1.114 0.5937 0.5942 96.0 1.987 0.5524 0.5552 204.0 2.312 0.5039 0.509.1 486.0 2.688 0.3974 0.4073 864.0 2.937 0.2905 0.3035 1314.0 3.119 0.2008 0.2148 1884.0 3.275 0.1262 0.1392 2274.0 3.357 0.0919 0.1035 Table 4.1: A comparison of the asymptotic and numerical results for xo(t) corresponding to the parameter values shown in the caption of Fig. 4.2. 107 and e = .03, a = 1.0, and D = 1.0. An interpolation scheme was then used to locate the position of the maximum of a on the computational grid. In Fig. 4.2 and in Table 1 we compare this numerical result for xo with the corresponding asymptotic result obtained^ by solving the differential equation (4.126c) with the initial condition XQ(0) = 0.6. In solving the differential equation, the integrals in (4.126c) were evaluated using Romberg integration on a large but finite interval using the form for uc given in (4.6). We find a close agreement between the asymptotic and numerical results for xo{t). In Fig. 4.3 we plot the initial condition used and then in Fig. 4.4(a) and Fig. 4.4(b) we plot the numerical solution to (1.17) at two different times showing the slow convergence to a one-spike equilibrium solution. 4.4.2 The Stability of a One-Spike Solution for D —)• oo When D — oo it was shown in §2.1.3 that a one-spike solution is metastable and that the center xo (t) of the spike satisfies the asymptotic differential equation ^ £ = G { x o ) „ / -2(i-x0)A _ c - 2 ( i + xo ) M ( 4 1 2 8 ) (LooMy)] dv) provided that XQ is not within 0(e) of the endpoints, i . e. (4.128) is valid when 1 — XQ 2> O(e) and 1 + XQ 3> O(e). Here a is defined by the limiting behavior uc(y) ~ ae - ' 2 ' ' as \y\ —>• oo. This result shows that xo = 0 is unstable and that there is a metastable drift of the spike towards the closest endpoint of the domain. When D is asymptotically large we can superimpose the result (4.128) with (4.126c) to obtain D ^ = G(xo)+F(x0). (4.129) at Here F(XQ) and G(XQ) are defined in (4.126c) and (4.128), respectively. This superposition is valid since the metastable interaction between the tails of the spike and the boundaries x = ±1 results in an additive term to the solvability condition that we impose on (4.119a). The stability property of the equilibrium solution for this differential equation is then given in the following result: 108 0.5 i -0.4 -x (a) t = 960 (b) t = 2400 Figure 4.4: Plot of a one-spike solution at two different times corresponding to the parameter values shown in the caption of Fig. 4.2. The solid curve is the activator concentration and the dotted curve is the inhibitor concentration. 109 Proposition 4.13 For e < 1, o one-spike equilibrium solution to (1-17) is stable when D < Di(e) and is unstable when D > Di(e), where nnf2 2/e /-oo D^~WTTwLlu'iv)r *• (4'130) Here a is defined by the limiting behavior uc ~ ae~^ as \y\ -¥ co, where uc(y) satisfies (4-5). For the special case with p = 1 and (p,q,r,s) = (2,1,2,0), where uc(y) = |sech2 (y/2) and a = 6, we can calculate analytically that Di(e) ~ e2e2£/125.0. 110 Chapter 5 Spike Dynamics The goal of this chapter is to study the dynamics of multi-spike solutions to (1.17). In §4 we found criteria determining the stability of a multi-spike equilibrium solution to (1.17). These criteria were derived by ensuring that the spectrum of the operator associated with a lineariza-tion about an equilibrium multi-spike solution contains no eigenvalues with positive real part. In §4 we examine two different types of eigenvalues. The stability of the equilibrium solution on an 0(1) time scale was determined by the sign of the real part of the large eigenvalues, and the stability on an 0(e - 2 ) time scale was determined by the sign of the 0(e2) eigenvalues. The 0(e2) eigenvalues were real. The stability of the equilibrium solution with respect to both sets of eigenvalues gave different ranges of D. Values of D that satisfy both ranges yield stable equilibrium spike solutions. In this chapter we linearize (1.17) around a quasi-equilibrium solution consisting of a sequence of spikes of different heights. As with the one-spike case treated in §4, the motion of such a multi-spike quasi-equilibrium solution is on a slow 0(e - 2 ) time scale. The quasi-equilibrium solution at a fixed time is stable on an 0(1) time scale when the large eigenvalues associated with the linearization are in the left-half plane. The outline of this chapter is as follows. In §5.1 we linearize (1.17) about a quasi-equilibrium n-spike solution where the height and the spike centers are unknown functions of time. By imposing a solvability condition, we obtain a differential-algebraic system of ordinary differential equations governing the motion of the spikes. In §5.2 we discuss the equilibrium solutions to this system which yield both the symmetric solution discussed in §4 and new asymmetric solutions 111 found in [50]. In §5.3 the results of §4.2 axe used to determine when a given spike profile is stable with respect to the large eigenvalues. In §5.4 numerical simulations of the full system are compared to numerical simulations of the system found in §5.1, and these results are compared with the stability results of §5.3 . 5.1 The Dynamics of Quasi-Equilibrium Solutions We derive a system of ordinary differential equations describing the dynamics of the spike locations for an n-spike quasi-equilibrium solution to (1.17). The spike locations X j are assumed to satisfy —1 < xz- < Xj+i < 1 for i = 1,..,n — 1. In §4.4.1 a one-spike solution was analyzed in detail and it was found that the spike evolves on a long time-scale t = 0(e~2). Hence, we expect that X j = Xj(r), where T = e2t. In the inner region near each #j, we introduce the new variables Hi = e_1(x - X j ) , hiiyur) = h(xi + eyi,e~2r), ai{yi, r) = a(x* + ey;, e _ 2r), r = e2t. (5.1) We then expand hi(yi,r) = hi0(yi,T) + ehn(yi,T) -I , (Hiyur) = ai0(yi) + ean(yi,T) -\ . (5.2) Substituting (5.2) into (1.17), we get to leading order that a'io ~ ai0 + aP0/hqi0 = 0, -co < y , < co, (5.3a) . h-0 = 0. (5.3b) Here the primes indicate derivatives with respect to y j . In order to match to the outer solution below, we need that hio is bounded and that aiQ -> 0 as | y j | -¥ co. Thus, from (5.3b), we get that hiQ = Hi(r), (5.4a) for some unknown Hi. The solution to (5.3a) is ai0(yi) =H]uc(yi), where 7 = 9/(p-l)> (5-4b) 112 and uc(y) satisfies uc — uc + upc = 0, —oo < y < co, uc 0 as |y|->oo; x/c(0) = 0, u c (0)>0. In particular, when p = 2, we have « c(y) = 2 s e c h 2 (y/ 2) • (5.5a) (5.5b) (5.5c) The 0(e) problems, obtained from substituting (5.2) into (1.17), are p—l p " , Paio Qaio u ' • a-il — ail T 7a—Oil = T^+T^l ~~ ai0 xi > rti0 «io "iO Here i j = dxi/dr. Substituting (5.4) into (5.6), we get L(an) = a'n - an +pupc 1aii = qH] luvchix - H]ucii, Dh'i = -HT-'ul. (5.6a) (5.6b) (5.7a) (5.7b) Since L(u'c) = 0, and u'c -> 0 exponentially as |y| -> oo, the right-hand side of (5.7a) must satisfy the solvability condition that it be orthogonal to u'c. In this way, we get |2 u'c(y) dy ~ — / t i ?«c /m dy* • o o L J -"t y—oo (5.8) Integrating (5.8) by parts twice, and using the facts that fin and uc are even functions, we obtain xi / uc(y) J—oo dy 2Hi(p + T ) ( / _ > W r l dy lim hix + lim . (5.9) Now in the outer region, defined away from 0(e) regions near each X{, a is exponentially small and we expand h(x) = h0(x) + o(e). (5.10) 113 Since e < 1, and a is localized, the term e~lam/hs in the outer region behaves like a linear combination of delta functions. Substituting (5.4) and (5.10) into (1.17b), and letting e -+ 0, we find that ho satisfies n poo Dh0xx - ftho = -K ]T H]T~s5(x - X i ) , br= [uc(y)]r dy, (5.11a) i=i J - ° ° h0x{±l)=0. (5.11b) The solution to (5.11) is n ho{x) = brYJHT~SG{x;xi), (5.12) i=l where G(x; Xi) satisfies DGXX - fiG = -S(x -Xi), -1 < x < 1, (5.13a) Gx{±l;xi) = 0. (5.13b) To match with the inner solutions near each Xi, we require for i = 1, ..,n that h0{xi) = Hi, (5.14a) lim h'a + lim h'n = hox(xi+) + hox(xi-) • (5.14b) yi->oo 2/j->—oo From (5.14a) and (5.12), we obtain a nonlinear algebraic system for Hi, i = 1,.., n. The final step in the derivation is to calculate the integral / defined by f°° M y ) F + 1 dy f=J~Z { J 2 (5.15) I-ooiu'M] dv To do so, we first multiply (5.5a) by uc. Upon integrating the resulting equation over the domain, we obtain /oo poo poo ucuc dy- u2cdy+ up+l dy = 0. (5.16) -oo J —oo J —CO Upon integrating the first term in this equation by parts, we get P I M y ) ] 2 dy -l = e-f, where e = 1 , . 5.17 . H o K O / ) ] 2 ^ 114 To obtain an additional equation, we multiply (5.5a) by u'c and integrate over the domain to fix the constant of integration. We then integrate the resulting expression again to get 2/ R e -solving (5.17) and (5.18), we obtain / = p + 1 2(p + l) p-1 (5.18) (5.19) The dynamics of the n-spike pattern is obtained by substituting (5.14) and (5.19) into (5.9) and (5.12). The main result is summarized as follows (relabeling Hi by hi in the notation below): Proposition 5.1 For e <C 1, the quasi-equilibrium solution for a and h is given by n a(x,t) ~ ac = Y^hjUc [e~l\x - Xj)] , i.=i . n h(x,t) ~ hc = br yy^KY~sG(x]Xj), where hi = hi(r) and x\ = Xi(r) satisfy the differential-algebraic system for i = 1, ..,n hi = br ^ 2 h]T sG(xi\Xj), 3=1 ( dxi 2qbr x hT~'(Gx)i +YJhJ-sGx(xi-Xj) (5.20a) (5.20b) (5.21a) (5.21b) dr p — 1 V Here uc satisfies (5.5), bm is defined in (5.11a), r = e2t, and {Gx)i = [Gx(xi+; Xi) + Gx(xi-\Xi)] /2. J The system (5.21) can also be written in matrix form as h = brgwr-s, dr p — 1 where dxjdr = (xx,... ,xn)'. Here we have defined 9 = G(xx;xi) ••• G(xi)xn) \y G(xn;x{) ••• G(xn;xn) J U = hi 0 0 '•• 0 0 0 0 ••• h. (5.22) (5.23a) 115 and (Gx(xi;xx))o ••• Gx(xi;xn) y Gx(xn;xi) ••• (Gx{xn;xn))n J \K J hTr-s _ / h y r s \ (5.23b) Alternatively, we can write (5.22) in terms of certain tridiagonal matrices. In Appendix C, we show that fjID (5.24) where B is the tridiagonal matrix B = ci di 0 di c2 0 "•• o ••. 0 0 0 0 0 0 0 0 0 0 0 0 C n _ l d n - l 0 dn-l Cn J (5.25a) with matrix entries defined by ci = coth [0(x2 — xi)} + tanh [0(1 + x\)] ; c^ = coth [0(xn — xn-i)] + tanh [0(1 — xn)} , (5.25b) Cj = coth [6(xj+i — Xj)] + coth [9(XJ - XJ-I)] , j = 2,.., n — 1, (5.25c) dj = — csch [0(XJ+I — Xj)] , j = l,..,n — 1. (5.25d) Here 0 = y/JTJD. Next, we calculate the matrix product VB using the procedure as outlined in Appendix A. We 116 find that VB is a tridiagonal matrix of the form VB = —Vb, where Vb ei fx 0 -fx e2 0 '•• o ••. 0 0 0 0 0 0 0 0 0 0 0 e n-i / n - l 0 -/n - i en J The matrix entries are defined by (5.26a) ei = tanh [0(1 + x\)] - coth [0(x2 - xi)] ; en = coth [0(:rn - xn_i)] - tanh [0(1 - xn)] , (5.26b) = coth [9(XJ - a;j_i)] - coth [0(XJ+\ - Xj)} , j = 2,.., n - 1, (5.26c) /j = csch [9(xj+i - XJ)] , j = 1,.., n - 1. (5.26d) Substituting (5.24) into (5.22), we obtain the following result equivalent to (5.21): Corollary 5.1 The differential-algebraic system (5.21) is equivalent to the matrix system (5.27) d x q l^n-'nh, Bh = -±=hr-. dr p- 1 V D" ' ' " yfpD Here K, B and Vb are given in (5.23a), (5.25), and (5.26), respectively. The advantage of this formulation over (5.21) is that (5.27) is expressed only in terms of tridiagonal matrices. Starting with certain initial data a;(0), in §5.4 we give numerical examples showing the evolution of the quasi-equilibrium solution (5.20) and (5.27) towards a stable equilibrium solution. These asymptotic results are also compared with full numerical results computed from (1.17). 117 5.2 Symmetric and Asymmetric Equilibria From (5.27), the equilibrium values of h and x satisfy Bh = -^=h^-s, Vbh = 0. (5.28) In this section we review some results obtained in §4 and in [50] for the existence and stability of symmetric and asymmetric spike patterns for (1.17) respectively. As shown in §4, for the symmetric spike patterns where hj = H for j = l , . . ,n , the Xj are located at the symmetry points Xj = -1 + ^-^-, j = l,..,fc. (5.29) In this case, it was shown in §4 that 2 tanh (0/n), with 0 = y/JI/D, is an eigenvalue of B with associated eigenvector v = (1,1,.., 1)*. In addition, V\,v = 0. Hence, from (5.28), the common spike value hj = H is ^ - ( • + D = H ^ t a n h ( r - V (5.30) The symmetric n-spike solution is obtained by substituting (5.29) and hj = H into (5.20). In [50], asymmetric n-spike equilibrium patterns were constructed asymptotically. The type of asymmetric patterns that were constructed consisted of n\ > 0 small spikes of type A and ri2 = n — ni > 0 large spikes of type B arranged in any particular order from left to right across the interval [—1,1] as A B A A B . . . . B , m - A's , n 2 - B's. (5.31) A plot of such a solution with five spikes in an ABBAB pattern is shown in Fig. 5.1. The following result for asymmetric equilibrium spike patterns was obtained in [50]. Proposition 5.2 (see [50]) Let e —> 0 and D < Dm, where Dm is some critical value. Then, there exists an asymmetric n-spike equilibrium solution to (1-17) of the form (5.20), where hj 118 0.10 0.075 h a, h 0.05 0.025 Figure 5.1: Plot of the activator and inhibitor concentration for a five-spike asymptotic asym-metric equilibrium solution with e — .02, D = .04, a — 1, and (p, q, r, s) = (2,1,2,0). The solid curve is the activator concentration and the dotted curve is the inhibitor concentration. satisfies h j r H , + l ) = (l^R^ t a n h { l . e ) > Q = y ^ / p . (5.32) Here for each j, lj = I or lj = I, where I and I are determined in terms of n\, n2 and y/nJD by the nonlinear system nil + n2l = l, b l^/a/D = b l^/a/D (5.33a) where 6(z) = tanh r z r = (5.33b) cosh z ' ' -yr — (s + 1) ' The value lj = I must occur n\ > 0 times, while lj = I must occur n2 = n — n\ > 0 times. The small and large spikes can be arranged in any sequence. Finally, the spike locations Xj are found from x\ = h - 1, Xk = l- lk, Xj+i = Xj + + lj , j = 1,.., k - 2. (5.34) 119 Detailed numerical computations for the critical value Dm and further more refined results were obtained in [50]. For our purposes, the key point concerns the relationship between the symmetric and the asymmetric spike patterns. Define an Z^-type norm of the equilibrium solution for a by n |a|i = ^ / i j , (5.35) j=i where |a|i is a function of D. Label the symmetric branch with n spikes by sn. Then, the following result was shown numerically in [50]: Proposition 5.3 (see [50]) An n-spike asymmetric solution branch with n\ small spikes of type A provides the connection as D is varied between the symmetric branches sn and sn-m. All of the asymmetric branches with n spikes bifurcate from the symmetric branch sn at the critical value D = given by Db = — r r 2 - (5-36) n 2 log ( V r + V r + l ) Here r is defined in (5.33b). The stability properties of asymmetric spike patterns was studied in [50]. In the analysis of [50], there were two classes of eigenvalues that needed to be considered. The first class are the large 0(1) eigenvalues, resulting from global spike interactions, that correspond to strong instabilities of small perturbations of the equilibrium solution on an 0(1) time-scale. This instability property, referred to here as profile instabilities, is not contained in the system (5.27). The second class of eigenvalues are the small 0(e2) eigenvalues that are associated with near translation invariance and slow dynamics near the equilibrium solutions. These small eigenvalues arise from the linearization of (5.27) about equilibrium values for h and x. Thus, these eigenvalues are contained in the system (5.27). Critical ranges of D that ensure the stability of the equilibrium solution with respect to both classes of eigenvalues were derived in [50]. 120 The main qualitative stability conclusions from the analyses of §4 for symmetric solutions and [50] for asymmetric solutions are as follows. For a certain nontrivial range of D near the bifurcation value Df,, each n-spike asymmetric solution branch that bifurcates from the symmetric branch sn is stable with respect to the large 0(1) eigenvalues. However, based on numerical evidence, each of these bifurcating asymmetric branches is always unstable with respect to the small 0(e2) eigenvalues. The symmetric branch is stable with respect to the large eigenvalues when D < Dn and is stable with respect to the small eigenvalues when D < Dj . Here < D n , and Dn is given by 1-2 Dn = Aun'2 ln (@ + y/(32 - l j , 6 = 1 + (1 + cos {%/k)) r, (5.37) for n > 2, where r is defined in (5.33b). This results predicts that D\ is infinite, but as shown in §4, D\ is exponentially large as e —> 0. In Fig. 5.2 we plot a bifurcation diagram of |a|i versus D, showing the asymmetric and symmetric branches with fewer than four spikes and their stability ranges with respect to the large 0(1) eigenvalues. 0.35 0.30 0.25 0.20 0.15 0.10 0.05 0.0 1 i i 1 101 - / / 01 ' /yooi<\".-\ Si -i 1 -0.0 0.1 0.2 0.3 D 0.4 0.5 Figure 5.2: Plot of |a|i defined in (5.35) versus D for solutions with three or fewer spikes. Here /j = 1 and (p,q,r,s) = (2,1,2,0). The symmetric branch with A; spikes is labeled by Sfe. The asymmetric patterns AB, BAB, and AAB are labeled by 01, 101, and 001, respectively. The portions of the branches that are solid (dotted) are stable (unstable) with respect to the large 0(1) eigenvalues. The implication of these results is that there are many equilibria of the quasi-equilibrium 121 dynamics (5.27) with exactly n interior spikes. One of them corresponds to the symmetric spike pattern, and the rest correspond to asymmetric patterns. However, only the symmetric branches will be stable with respect to both the large and small eigenvalues when D < Df,. Hence, it is reasonable to expect that when D < the quasi-equilibrium dynamics (5.27) starting from some x(0), will tend to a symmetric equilibrium with n spikes. Based on the numerical evidence of §5.4, this scenario will only occur if the quasi-equilibrium solution is stable with respect to the large 0(1) eigenvalues throughout the slow dynamics. For symmetric equilibria, this stability threshold is given by Dn in (5.37). For the quasi-equilibrium solution, the stability threshold depends on the values of h and a; at a given r. 5.3 Stability of the Profile: The Large Eigenvalues We now examine the stability, at a fixed value of T , of the quasi-equilibrium profile constructed in §5.1. The quasi-equilibrium profile of §5.1 varies on a slow time-scale r = e2t. We would like to determine whether this profile can undergo an instability on a fast time-scale of 0(1). Hence, since there is a time-scale separation, in the eigenvalue analysis below we can treat T as being a fixed parameter. To derive the eigenvalue problem, we substitute a(x, t) = ac + eM4>{x), h{x, t) = hc + extv(x), (5.38) into (1.17) where ac and hc are defined in (5.20), and rj 1 and cf> <C 1. This leads to p — l p e2fax-<P + P-^<t>-^r) = \<t>, - \ < x < \ , (5.39a) T — 1 T DVxx-m = -e~1r^^-4> + e~1s-£^Ti, - 1 < x < 1, (5.39b) <t>x{±l) = 7 f c ( ± l ) =0. ' (5.39c) The spectrum of (5.39) contains large eigenvalues that are 0(1) as e ->• 0. In §5.4, we show that the quasi-equilibrium profile becomes unstable on an 0(1) time-scale when the Xi(r), hi(r), for i = 1,.., n, are such that (5.39) has an eigenvalue with Re(A) > 0. 122 We look for an eigenfunction of (5.39) in the form k j=l (5.40) where 4>j(y) -» 0 exponentially as \y\ -> co. Then, the right-hand side of (5.39b) with 5 = 0 behaves like a sum of delta functions when e <^ 1. Substituting (5.20) and (5.40) into (5.39b), we get Dr)x Vx(±l) = 0, j] = - ^TOJJ5{X - Xj). 3=1 where /•oo Uj = rhj~l'~s / uT~l4>j dy. J — O O This problem is equivalent to Drjxx - /J77 = 0 , -1<X<1; nx(±l) = Q, [v)j = 0, [DVx]j = -UJ + sbrh]r-{s+l)v(xj), (5.41a) (5.41b) (5.41c) (5.42a) (5.42b) where [v]j = V(XJ+) — V(XJ-). By solving this system on each subinterval as in Appendix A, we can show that ( v(xi) ^ \ n{xn) ) U7 1 ^ (5.43) Here T-L and B are defined in (5.23a) and (5.25), respectively. Since s > 0 and ri is a positive matrix, we can solve for 77 as rbr („ , sbT r ] - ^ D \ b + V H D n ) U I Sturdy - 1 'fZoUc W (5.44) Next, we substitute (5.20) and (5.40) into (5.39a) to obtain for j = 1, ..,n that 4'j - <t>j +Pupc~1(pj - qh^vP^ixj) = \(f>j , -co < y < 00 , (5.45) 123 with 4>j(y) —> 0 as y -» oo. We can write (5.45) in matrix form as <p" - 0 + pu p - x 0 - qulW-lT) = \<p, (5.46) Substituting (5.44) into (5.46), we obtain the eigenvalue problem <p -<t> + pvP-l<l>-TqvPA joo\rdy J=W, - o o < y < o o , (5.47a) 4> —> 0, as \y\ —> oo , (5.47b) where the matrix 5 is defined by £ = -^.W~x (l3 + -^U^-^A 1 Kr-b+Vu1-"!. (5.48) Since Q = B~l/y/pD is positive definite and It is a positive diagonal matrix, we conclude that £ has real positive eigenvalues. We decompose £ as £ = S~1AeS, (5.49) for some nonsingular matrix S. Then, upon defining ip = S<f>, we obtain from (5.47) that xp" - ^ +pug"V ~ rgug J-°° c y y)=\*l>, - o o < y < o o , (5.50a) ip -¥ 0, as \y\ oo. . " (5.50b) Since A e is a diagonal matrix we obtain n uncoupled problems from (5.50). The next step is to determine the conditions for which Re(A) < 0 in (5.50). For this we use Theorem 2.1. By comparing (5.50) with (2.11), we obtain the following result on the spectrum of (5.50): Proposition 5.4 Let XQ ^ 0 be the eigenvalue of (5.50) with the largest real part and assume condition (2.13a) holds. Let a\ be the minimum eigenvalue of the matrix £ defined in (5.1^8). Then, .Re(Ao) > 0 when . a x < ^ . (5.51) qr Also, Re(\o) < 0 when a\> (p — l)/qr. 124 A more convenient criterion is obtained by writing £ in (5.48) as £ = W where £ = 4^~7 H" ( 3 + 1 )£- (5-52) The eigenvalues of £ and are identical. We can then rewrite Proposition 4.1 as the following simple criterion: Corollary 5.2 Let Ao ^ 0 be the eigenvalue of (5.50) with the largest real part and assume condition (2.13a) holds. Let em be the maximum eigenvalue of the tridiagonal matrix £ defined in (5.52). Then, Re(Xo) > 0 when e m > £ ? T ) " s : ( 5 - 5 3 ) Also, Re(Xo) < 0 when em < qr/[p — 1) — s. The results of §4 and [50] for the stability of symmetric and asymmetric equilibrium spike patterns, respectively, with respect to the large 0(1) eigenvalues were obtained from a criterion such as Corollary 4.1. In our numerical examples in §5.4 of the evolution of the quasi-equilibrium solution, we track the maximum eigenvalue of £ as a function of r and determine the behavior of the solution if the threshold in (5.53) is exceeded. 5.4 Comparison of Asymptotic and Numerical Results We now compare the asymptotic results for the spike motion with corresponding numerical results computed from (1.17). Unless otherwise stated, in the comparisons below we have taken the parameter values p, = 1, the exponent set (p,q,r,s) = (2,1,2,0), and e = 0.02. The comparisons are made for various values of the inhibitor diffusivity D, of the number of spikes n, and of the initial spike locations Xj(0), for j = 1, ..,n. The time variable given in the plots below correspond to the slow time variable r denned in (5.1) by r = e2t. With e = .02, we get t = 2500T. 125 To compute the full numerical results from (1.17), we use the NAG library code [39] with 2001 equidistant meshpoints. For given values of D, n, and Xj(0) for j = 1, ..,n, we take the quasi-equilibrium solution (5.20) to be the initial condition for a and h. To compute a(ar, 0) and h(x,0), we must determine the initial values hj(0), for j = 1, ..,n, from the nonlinear algebraic system for hj in (5.27). This is done using Newton's method. In the Newton iteration we require an initial guess for the hj. This is quite a nontrivial task, and to do so we use one of two different methods. The first method is to perform a homotopy in the spike locations starting from the symmetric equilibrium solution. The second method is a homotopy in the value of D, starting from a large value of D where hj « h^, for some easily determined /loo, and then to decrease D to the desired value. Sometimes one of these homotopy approaches failed. This is because the linearized system for the hj is not invertible at certain specific parameter values, as we explain below. Once the nonlinear system for the hj in (5.27) is solved, the initial condition for a(x, 0) and h(x, 0) is known from (5.20) and the NAG solver [39] is used to compute the solution to (1.17) at later times r. The locations of the spikes from these numerical results were obtained by a local quadratic interpolation. The asymptotic results were obtained by solving the differential-algebraic system (5.27) for XJ(T) and hj[r) using the ODE solver [42] coupled together with Newton's method to determine the hj. Another method, which we found to work equally well, is to use the differential-algebraic solver DDASSL [5] directly. The initial values for hj (0) were obtained using one of the homotopy methods described above. It is very important to calculate the maximum eigenvalue em of (5.52) as a function of T . The curve em = eTO(r) is computed at each step by using the eigenvalue solver [2]. The asymptotic theory predicts that if em > 2 at any value of r, the quasi-equilibrium profile develops an instability on an 0(1) time scale. We would like to check this asymptotic prediction with the full numerical results. In §5.4.1 and §5.4.2 we give asymptotic and numerical results for the case of n = 2 and n = 3, respectively. Other cases, and some general results, are discussed briefly in §5.4.3. 126 5.4.1 Two Spikes n = 2 When n = 2, the nonlinear system in (5.27) for hj is ci/ i i +dih2 = brhJT~s / y/JiD, (5.54a) di / i i + c 2/i 2 = brhY~" I\fpD, (5.54b) where Cj and dj are defined in (5.25a). From (5.54) we readily obtain an equation for h\jh2. The differential equations for Xj in (5.27) depend only on the ratio h\/h2. In this way,, we obtain the next result. Corollary 5.3 When n = 2, the differential-algebraic system (5.27) is equivalent to xi =-C{tanh[0(l + xi)]-coth[e{x2-xi)] +csch[9{x2-x1)]/^) , . (5.55a) X2 = - C {-csch{6{x2 - xx)\ji + coth[0(x2 - ari)] - tanh[0(l - x2)}) , (5.55b) where £ = h\jh2 satisfies f(0 = csch[6(x2 - an)] (p™-** - 1) + coth[0(x2 - xi)] (£ - Cm~s) + tanh[0(l + xi)]£ - tanh[0(l - x2)]Cm~' = 0. (5.55c) .Here £ = <7#/(p — 1), «//iere # = ^Jn/D. The behavior of the spike dynamics depends on the value of D with respect to Dj , the initial position of the spikes, and the initial value of em at T = 0. Recall from §4 and (4.116a) that a two-spike symmetric equilibrium solution with spikes located at xx = —x2 = —1/2 will be unstable if D < Db « 0.3218. In the special case of symmetric initial data where xi(0) = —x2(0), we have that £ = 1 is a root to (5.55c). In this symmetric case, £ = 1 and x i = - x 2 for all r. Combining (5.55a) and (5.55b) we obtain the next result. Corollary 5.4 Let n = 2 and assume that the initial data is symmetric in the sense that 127 #i(0) = — 2:2(0). Let y = X2 — x\. Then, (5.55) is equivalent to the single differential equation y = 2((csch(9y)-coth(9y) + ta,nh[9(l-y/2)}) , (5.56) where £ = q9/(p — 1) and 9 = y/p/D. We now gives some examples of the theory. Example 1 (Symmetric Initial Data): We consider two different cases of symmetric initial data with xi(0) — — 2:2(0). The parameter values are Example la: xx (0) = -0.4 = -2:2(0), 0 = 0.4, (5.57a) Example lb: Xl{0) = -0.2 = -x2{0), 0 = 0.1, (5.57b) For Examples l a and lb we compute numerically that em(0) = 1.944 and em(0) = 1.793, respectively. Thus, for both examples, the initial profile is stable with respect to the large 0(1) eigenvalues. In Example lb, D < Db and in Example la, D > I V Thus, in Example la the symmetric equilibrium solution is unstable with respect to the small eigenvalues. I 1 1 2.25 1 : 1 1 (a) Xj versus r (b) em versus r Figure 5.3: In Fig. 5.3(a) we plot Xj versus r for the two symmetric parameter sets in Example 1. The solid and heavy solid curves correspond to the full numerical results for Example la and Example lb, respectively, given in (5.57). The asymptotic results computed from (5.56) correspond to the dotted curves. In Fig. 5.4(b) we plot the maximum eigenvalue em of (5.52) versus T . The solid and heavy solid curves refer to Example la and Example lb. 128 r x\ (num) x2 (num) xi (asy) x2 (asy) 0.06 -0.4163 0.4163 -0.4157 0.4157 0.09 -0.4238 0.4238 -0.4225 0.4225 0.18 -0.4419 0.4419 -0.4399 0.4399 0.27 -0.4557 0.4557 -0.4535 0.4535 0.36 -0.4662 0.4662 -0.4639 0.4639 0.45 -0.4742 0.4742 -0.4721 0.4721 0.54 -0.4803 0.4803 -0.4783 0.4783 0.63 -0.4850 0.4850 -0.4832 0.4832 Table 5.1: The numerical and asymptotic results for x\ and x2 versus r for Example la. The asymptotic result for these two examples is obtained by integrating (5.56). From (5.56) we observe that, for any initial data y(0) > 0 and any value of D, the solution satisfies y ->• 1 as T —^ oo. Hence, for symmetric initial data, x\ —> —1/2 and x2 -> 1/2 as r —>• oo for both Example la and Example lb. To attempt to explain this paradox, we refer to Proposition 4.11 of §4 where it was shown that when D > D\, there is exactly one positive and one negative small eigenvalue for the linearization of (1.17) about the two-spike symmetric equilibrium solution. The positive eigenvalue becomes negative when D crosses below Db- Thus, when D > Db, the symmetric equilibrium corresponds to a saddle point in phase-space with one only stable direction towards the equilibrium. Thus, we conjecture that by taking symmetric initial data we approach the equilibrium solution along the stable manifold. In Example 4 below we show that we do not tend to a symmetric equilibrium solution when we make a slight perturbation to symmetric initial data when D > Db-For Examples l a and lb, in Fig. 5.3(a) we compare the trajectories Xj(r) for j = 1,2 computed from the asymptotic result (5.56) with the full numerical results computed from (1.17). Selected values for Xj from the asymptotic and full numerical results are shown in Table 5.1 and Table 5.2 for Examples la and lb, respectively. The agreement is found to be close. In Fig. 5.3(b) we plot em = em(r) for both examples. The threshold for an 0(1) instability is given by the 129 r x\ (num) x 2 (num) xi (asy) x 2 (asy) 0.06 -0.2762 0.2762 -0.2670 0.2670 0.12 -0.3283 0.3283 -0.3146 0.3146 0.24 -0.3946 0:3946 -0.3779 0.3779 0.36 -0.4337 0.4337 -0.4177 0.4177 0.42 -0.4472 0.4472 -0.4321 0.4321 0.48 -0.4579 0.4579 -0.4439 0.4439 0.66 -0.4786 0.4786 -0.4681 0.4681 0.78 -0.4863 0.4863 -0.4781 0.4781 Table 5.2: The numerical and asymptotic results for x\ and x 2 versus r for Example lb. dotted line in Fig. 5.3(b). Notice that em decreases as r increases, so that no 0(1) instability occurs at later times. In Fig. 5.4(a) and Fig. 5.4(b) we plot the numerical solution a versus x at different values of r for Example la and Example lb, respectively. Example 2: (Generic Initial Data) We consider two sets of generic initial data. The parameter values are taken to be Example 2a: xi(0) = -0.2, x2(0) = 0.32, D = 0.1, (5.58a) Example 2b: xi(0) = -0.8, x2(0) = 0.2, D = 0.2, (5.58b) For Examples 2a and 2b we compute numerically that em(0) = 1.4846 and em(0) = 1.4106, respectively. Thus, for both examples, the initial profile is stable with respect to the large 0(1) eigenvalues. For both examples D < Dj so that the symmetric equilibrium solution is stable. For Examples 2a and 2b, in Fig. 5.5(a) we compare the trajectories XJ(T) for j = 1,2 computed from the asymptotic result (5.55) with the full numerical results computed from (1.17). Selected values for Xj from the asymptotic and full numerical results are shown in Table 5.3 and Table 5.4 for Examples 2a and 2b, respectively. As seen in Fig. 5.5(b) the maximum eigenvalue once again decreases as r increases and so no 0(1) instability is triggered. In Fig. 5.6(a) and 130 0.10 0.05 1— 1 >! jj J 1,1 (a) a versus x at different r (Example la) (b) a versus x at different r (Example lb) Figure 5.4: The numerical results for a versus x at different values of r are plotted in Fig. 5.4(a) for Example l a and in Fig. 5.4(b) for Example lb. In Fig. 5.4(a) the solid, dotted, and heavy-solid curves correspond to r = 0.0, r = .202, and r = 1.202, respectively. In Fig. 5.4(b) the dotted, solid, light dotted, and heavy solid curves correspond to r = 0.0, r = .123, r = .321, and r = 1.002, respectively. Fig. 5.6(b) we plot the numerical solution a versus x at different values of r for Example 2a and Example 2b, respectively. Example 3: (An Initial 0(1) Instability) We take the parameter values a;1(0) = -0.2, x2(0) = 0.32, £ = 0.25. (5.59) In this case, we compute numerically that em(0) = 2.228, so that the initial profile is unstable with respect to the large positive eigenvalues. Since em(0) > 0, the asymptotic result (5.55) does not apply. In Fig. 5.7(a) we plot the numerically computed a versus x for r = 0 and for several very small values of r. In terms of t, the plots in Fig. 5.7(a) are for t < 20.0. What we observe is that the smaller of the initial two spikes starts collapsing but that its locations remains quite fixed during the collapse. The other spike grows during the collapse of its neighbor. This type of collapse behavior is qualitatively very different than internal layer collapse behavior for phase transition models (see [47]). After one of the spikes has collapsed, the other one moves very slowly towards the symmetric one-spike equilibrium solution centered at the origin. This is 131 em (a) Xj versus r (b) em versus r Figure 5.5: In Fig. 5.5(a) we plot Xj versus r for the two parameter sets in Example 2. The solid and heavy solid curves correspond to the full numerical results for Example 2a and Example 2b, respectively, given in (5.58). The asymptotic results computed from (5.55) correspond to the dotted curves. In Fig. 5.5(b) we plot the maximum eigenvalue em of (5.52) versus r. The solid and heavy solid curves refer to Example 2a and Example 2b. 0.25 (a) a versus x at different r (Example 2a) (b) a versus x at different r (Example 2b) Figure 5.6: The numerical results for a versus x at different values of T are plotted in Fig. 5.6(a) for Example 2a and in Fig. 5.6(b) for Example 2b. In Fig. 5.6(a) the dotted, solid, light dotted, and heavy solid curves correspond to r = 0.0, r = .202, r = .403, and r = 1.402, respectively. In Fig. 5.6(b) the dotted, solid, light dotted, and heavy solid curves correspond to T = 0.0, r = .404, r = .804, and r = 3.200, respectively. 132 T x\ (num) X2 (num) xi (asy) x2 (asy) 0.032 -0.2329 0.3526 -0.2284 0.3480 0.08 -0.2724 0.3906 -0.2634 0.3817 0.30 -0.3782 0.4789 -0.3622 0.4667 0.40 -0.4051 0.4954 -0.3892 0.4853 0.48 -0.4218 0.5036 -0.4063 0.4953 0.60 -0.4410 0.5104 -0.4265 0.5048 0.80 -0.4619 0.5138 -0.4500 0.5118 1.28 -0.4856 0.5100 -0.4788 0.5117 Table 5.3: The numerical and asymptotic results for xi and x2 versus r for Example 2a. shown in Fig. 5.7(b). Example 4: (An O(l) Instability at a Later Time) We consider a case with near sym-metric initial data, where we take Xl(0) = -0.4, x2(0) = .401, D = 0A. (5.60) In this case, we compute that em(0) = 1.941, so that the initial profile is stable with respect to the large eigenvalues. However as time increases, the profile will become unstable with respect to the large eigenvalues before reaching equilibrium. The results from these computations are plotted in Figures 5.8(a) and 5.8(b). Initially, the numerical results from the full system and the asymptotic system agree. However, at a later time these results start to diverge and one of the spikes from the numerical computation becomes unstable before the predicted asymptotic time. We speculate that the reason for this is that the nearly symmetric initial data is very close to the stable manifold and thus remains close to the stable manifold for some time. As the solution moves away from the stable manifold, errors in the simulation as well as differences between the asymptotic and full system begin to grow rapidly. We conjecture that it is this growth in errors that is responsible for the discrepancies between the simulation of the full system and that of the asymptotic differential equations. 133 (a) o versus x (collapse event) (b) a versus x (one-spike motion) Figure 5.7: The numerical results for a versus x at different values of r are plotted for Example 3. In Fig. 8a, where we show the collapse event, the solid, dotted, and heavy solid curves correspond to r = 0.0, r = 0.00504, and r = 0.00704, respectively. In Fig. 8b, where one spike drifts towards the origin, the solid, light dotted, and heavy solid curves correspond to r = 0.4, T = 1.6, and T = 4.8, respectively. The initial condition at r = 0 is the dotted curve. (a) Xj versus T (b) e m versus r Figure 5.8: In Fig. 5.8(a) we plot Xj versus r for Example 4. The solid curve corresponds to the full numerical results and the dotted curves corresponds to the asymptotic results. In Fig. 5.8(b), we plot the maximum eigenvalue e m of (5.52) versus r. 134 T xi (num) X2 (num) xi (asy) x2 (asy) 0.06 -0.6556 0.2161 -0.6528 0.2157 0.66 -0.6171 0.3400 -0.6136 0.3383 0.84 -0.6003 0.3678 -0.5981 0.3657 1.00 -0.5866 0.3891 -0.5854 0.3868 1.40 -0.5585 0.4299 -0.5589 0.4273 1.60 -0.5476 0.4445 -0.5485 0.4420 1.90 -0.5347 0.4611 -0.5359 0.4588 2.50 -0.5180 0.4808 -0.5192 0.4793 Table 5.4: The numerical and asymptotic results for xi and x2 versus r for Example 2b. 5.4.2 Three Spikes n = 3 Example 5 (Symmetric Initial Data): We now consider the case of three spikes with symmetric initial data (i.e. £1(0) = —£3(0) and x2(0) = 0). In this case the two outer spikes will be of the same height, however the middle spike will generally be of a different height. Thus, unlike the case of two spikes with symmetric initial data, we can not reduce this differential algebraic system to a single ordinary differential equation. Another important difference between the case of two and three spikes with symmetric initial data is in the form of the stable and unstable manifolds with respect to the small eigenvalues. In the case of two spikes both the stable and unstable manifolds have one dimension. In the case of three spikes, we have a one dimensional stable manifold and a two dimensional unstable manifold. In this case, it appears that symmetric solutions are not confined to the stable manifold. To illustrate this, in this example we take the values xi(0) = -0.5, x2(0) = 0, x 3(0) = 0.5, D = 0.01. (5.61) The results of the simulations are given in Fig. 5.9. 135 T xi (num) X2 (num) xi (asy) x2 (asy) 0.004 -0.40126 0.40226 -.40762 0.42297 0.032 -0.40894 0.41053 -.40618 0.42048 0.052 -0.41390 0.41614 -.40971 0.42662 0.08 -0.42032 0.42351 -.41423 0.43476 0.1 -0.42455 0.42844 -.41717 0.44026 0.3 -0.45404 0.46566 -.43574 0.48466 0.6 -0.47031 0.49756 -.43768 0.53190 0.988 -0.46514 0.52649 -.40640 0.59999 Table 5.5: The numerical and asymptotic results for x\ and x2 versus T for Example 4. Example 6 (Generic Initial data) We now compare simulations of (1.17) with (5.27) for the case of generic initial data. We use the parameters, xi(0) = -0.5, ar2(0)=0.1, x 3(0)=0.5, £> = 0.04. (5.62) For these parameters, we calculate em(0) = 1.6935. As seen in Fig. 5.10(b) em once again decreases as r increases, thus the profile is stable with respect to the large eigenvalues and no 0(1) instability is triggered. Asymptotic and numerical results for Xj are compared in Table 5.7. 5.4.3 Conclusions In this section we give two results related to the dynamics of n-spike profiles which are not of sufficient importance to warrant a section of thier own. Firstly, we note that the simulation in Fig. 5.4.1 stops just as em approaches the critical value. The simulation stops here not only because the system is no longer valid, but at this value the nonlinear solver for the h system fails to converge. We will demonstrate that when em is at the critical value the system is no longer solvable. Finally we discuss the role of the 0(1) eigenvalues on the dynamics of the system. 136 r r (a) Z j . versus r (b) em versus r Figure 5.9: In Fig. 5.9(a) we plot x3 versus r for Example 5. The solid curve corresponds to the full numerical results and the dotted curves corresponds to the asymptotic results. In Fig. 5.9(b), we plot the maximum eigenvalue em of (5.52) versus T . (a) Xj versus r (b) em versus r Figure 5.10: In Fig. 5.10(a) we plot Xj versus r for Example 6. The solid curve corresponds to the full numerical results and the dotted curves corresponds to the asymptotic results. In Fig. 5.10(b), we plot the maximum eigenvalue e m of (5.52) versus r. 137 T x\ (num) X2 (num) xz (num) xi (asy) x2 (asy) £3 (asy) 0.004 -0.50134 0.0 0.50134 -0.50053 0.0 0.50053 0.1 -0.52484 0.0 0.52484 -0.51239 0.0 0.51239 0.2 -0.54443 0.0 0.54443 -0.52340 0.0 0.52340 0.3 -0.56064 0.0 0.56064 -0.53328 0.0 0.53328 0.4 -0.57434 0.0 0.57434 -0.54224 0.0 0.54224 0,5 -0.58615 0.0 0.58615 -0.55041 0.0 0.55041 0.6 -0.59638 0.0 0.59638 -0.55792 0.0 0.55792 0.7 -0.60535 0.0 0.60535 -0.56485 0.0 0.56485 Table 5.6: The numerical and asymptotic results for x\, x2 and X3 versus r for Example 5. Solvability of (5.21) at critical values of em We now demonstrate why the system (5.21a) is not solvable at the critical value of em. We begin by considering the system linearized about a solution. We begin by linearizing the system about h = h + 6, (5.63) where h is a solution at the critical value of em and 6 C / i . Substituting into (5.21) we get B _ b m ( 7 m - 5 ) ^ 7 m - s _ A 0 = o ( 5 6 4 ) y/uD J From (5.52) and (5.53) we can see that this equation has a non-trivial solution at exactly the critical value of em. This implies that the Jacobian of (5.21a) is not invertible when e = em. Thus, this system does not have a unique solution at this point. Role of the O(l) Eigenvalues in the Dynamics From the numerical experiments shown above, and further computations, we speculate that when we are above the threshold of 0(1) stability from (5.53) for an n-spike profile, at least one of the spikes collapse as the remaining spikes grow. At the end of this collapses time interval, 138 T xi (num) x 2 (num) X3 (num) x i (asy) x 2 (asy) £ 3 (asy) 0.004 -0.50186 0.09228 0.40912 -0.50112 0.09354 0.40539 0.3 -0.62535 -0.00145 0.61632 -0.60449 0.00270 0.59417 0.6 -0.66074 -0.00230 0.65767 -0.65005 -0.00103 0.64704 0.9 -0.66631 -0.00143 0.66480 -0.66286 -0.00103 0.66150 1.2 -0.66692 -0.00818 0.66609 -0.66597 -0.00072 0.66520 1.5 -0.66688 -0.00046 0.66642 -0.66665 -0.00047 0.66617 1.8 -0.66680 -0.00026 0.66654 -0.66676 -0.00031 0.66645 Table 5.7: The numerical and asymptotic results for x i , x 2 and X3 versus r for Example 6. the profile again becomes stable, and the remaining spikes will then follow the dynamics of (5.21). 139 Chapter 6 Spike Pinning for the Gierer-Meinhardt Model In this chapter we examine the effects of spatially varying coefficients in (1.18) and (1.17). The spatial variations can results from pre-existing polarities in the cells and can have dramatic effects on the dynamics and equilibrium positions of the cell. We now give a detailed outline of this chapter. In §6.2 we examine the effect of a spatially variable inhibitor decay rate / i = (j,(x) > 0 and a spatially variable activator decay rate V(x) on the dynamics and equilibrium position of a one-spike solution to (1.17). In the biological context, these terms are examples 'of precursor gradients. The significance of precursor gradients from a biological viewpoint is discussed in [16]. From a mathematical viewpoint, we show that the effect of a spatially varying n is to perturb the equilibrium location for a one-spike solution away from the midpoint of the interval. The exact equilibrium location now depends on certain global properties of p,(x) over the domain. The effect of a spatially varying activator decay rate also perturbs the equilibrium location, but the perturbation depends on local properties of V(x). In §6.1 we consider a related problem with a spatially inhomogeneous coefficient. Specifically, we now let /j, be spatially uniform and we consider the shadow problem (1.18). A weak spatially inhomogeneous diffusivity for the activator equation is then introduced. This leads to the 140 perturbed shadow problem, defined by the scalar nonlocal problem e2 ap o-t = — - a + TT > - 1 < X < 1 , t>0, (6.1a) h={^-j ardxy+l , ax(±l,t)=0. (6.1b) The motivation for this form of K is mentioned below in context with a related problem studied in [44] for the Ginzburg-Landau equation. When K = 1 in (6.1a) it was shown using formal asymptotic in §2, and then proved later in [7], that (6.1) admits a one-spike solution that drifts exponentially slowly towards the closest endpoint of the domain. In addition, the equilibrium one-spike solution, which is centered at x = 0, is unstable. The metastable behavior associated with (6.1) when K = 1 suggests that exponentially small changes in K — 1 should influence the metastable dynamics greatly. Therefore, in §2 we study the dynamics of a one-spike solution for (6.1) as e -> 0 for a K(X; e) of the form K(x;e) = l + eug{x)e-*~ld . (6.1c) Here v and d > 0 are constants and g(x) is smooth. When g" (x) < 0 and 0 < d < dc, where dc is some constant, we show in §2 that (6.1) can have a stable spatially inhomogeneous equilibrium one-spike solution where the spike is centered at a zero of g (x). Results for this problem are given in Propositions 6.1 and 6.2 below. This effect whereby a localized structure is stabilized by a weak but spatially inhomogeneous coefficient in the differential operator is called pinning. The effect of weak spatially inhomogeneous terms has been examined in other contexts, includ-ing the pinning of vortices in superconductivity (cf. [6], [29]) and the pinning of an interface for the Ginzburg Landau equation posed in a thin cylinder of revolution, modeled by (cf. [44]) e 2 ut = ^ (Aux)x + 2(u - u3), -1< x < 1; u x (±l ,£) = 0. (6.2) This equation was derived in the thin channel limit in [44]. In this context, x is the direction along the axis of the cylinder and A = A(x; e) denotes the slowly varying cross-sectional area of the cylinder. For this problem, it was shown in [44] that an internal layer solution can be stabilized by an exponentially weak non-convex perturbation of a straight cylindrical domain 141 by taking A(x;e) to have the form A(x;e) = 1 + el/e~d/ig(x), where g"(x) < 0 and d > 0. A similar pinning result is obtained for spike solutions of the related problem (6.1). 6.1 One-Spike Dynamics: The Perturbed Shadow Problem In this section we construct a one-spike solution to (6.1). When K = 1 in (6.1), the resulting unperturbed problem admits a quasi-equilibrium one-spike solution of the form (see §2) a~ae(x;xQ) = hi^-^uc[e-1(x-x0)] , (6.3a) / br \ (3+l)fp~- l ) -«r poo e = {i) ' w h e r e br=J ^ c ( y ) ] dy- ( 6 - 3 b ) Here uc(y) is the unique positive solution to // uc -uc + upc = 0, -oo < y < oo, (6.4a) u'c(0) = 0; y,c{y) ~ a e~l yl, as y - > ± o o , (6.4b) for some a > 0. The function a e is called the quasi-equilibrium one-spike solution to (6.1) since ae has exactly one localized maximum and it fails to satisfy the steady-state problem corresponding to (6.1) by only exponentially small terms as e —> 0 for any value of XQ in \xQ\ < 1. For the time-dependent problem, our goal is to derive an asymptotic differential equation for the trajectory x = xo(t) of the localized maximum of ae, which represents the center of the spike. Since K — 1 is exponentially small, we begin by linearizing (6.1) in the form a(x,t) = ae[x;xo(t)] + w(x,t), (6.5) where w <^ ae. Substituting (6.5) into (6.1), we get the linearized problem for w Lew = dtae + dtw — e2 KXK~1 dxae, —1 < rc < 1, t > 0 , (6.6a) wx{±l,t) = -dxae{±l;x0), (6.6b) where the nonlocal operator L€ is defined by Le4>= e-(K(f>x)x-<f> + pUP-icf>- r g 6 f1 u'-'tdx. ' (6.6c) K x br{S + 1) J_i 142 Here br is defined in (6.3b). The coefficients in the operator L e are localized near x = XQ since uc = uc [e _ 1(x - x 0)]. 6.1.1 Exponentially Small Eigenvalue Let xo S (—1,1) be fixed and consider the eigenvalue problem Le<f> = \<f> - K x < l ; & . ( ± 1 ) = 0 . (6.7) When K = 1, it was shown in §2 that, under the conditions on p and r given in (2.13a) below, the eigenvalue of (6.7) with the largest real part is exponentially small as e —> 0 and the corresponding eigenfunction is localized near x = XQ. We will estimate the change in this exponentially small eigenvalue as a result of the exponentially small perturbation K — 1. To do so, we let K — \ and proceed as in §2 by introducing the localized eigenvalue problem for $ = $(y) given by £o$ - cr$ — co < y < oo; 3>—>•() exponentially as y -> ±co , (6.8) where y = e - 1 (x — xo). Here £Q is defined by Co^^^yy-^+pu^-r^^f^ul-^dy, (6.9) where uc = uc(y). Results on the spectrum of (6.9) were given in Theorem 2.1 in §2. The effects of the exponentially small perturbation K — 1 and of the finite domain in (6.7) perturb this eigenvalue <JQ by exponentially small terms as e -> 0. Let the perturbed eigenpair be denoted by Ao and <po- In order that (po satisfy the boundary condition on x = ±1 , we proceed as in §2 by constructing <po in the boundary layer form <p0 ~ u'c [e-\x - x0)] + (PLO [e _ 1 ( l + *)] + <Pm [e _ 1 ( l - , (6.10a) where 4>LO(V) =ae-€~^1+X0U-\ cj>w(r]) = - a e - ^ l - ^ e ^ . (6.10b) 143 To estimate Ao we first define the inner product (u,v)K = UVKCLX. Integrating by parts, we obtain for any two functions <p and ip that ( I e ^ ) K = eM -^#*)|li + (&£eV) , (6.11) where L* is the adjoint operator defined by = - (Kipx)x - iP+pvP-1^ - [9€~1U2~\ f Kufydx. (6.12) Now let ip = u'c and </> = fa. Then, since 0o:c(±l) = 0, we get from (6.11) that A0 (fa,u^ + €Kfauc\l_x = 6 ^ " ^ T 2 " ' ^oj - b r ^ € + ^ ^ 0o"c _ 1 ^ «u£u'ccfo^ .. (6.13) The two terms on the left side of (6.13) can be estimated as in §2 by using (6.10) for fa and (6.4) for u"c(y) as y —>• ±co. The exponentially small perturbation K — 1 does not affect the leading order asymptotic estimates for these quantities. Thus, for e —> 0, we get V 0 ' " ' 0 ) * ~ e / S ° ' = uc(y)j rfy. (6.14a) H i ~ 2ea2 (g-^-Hi+^o) + e - 2 e - i ( i - * 0 ) ) ( 6 1 4 b ) Next, we use (6.1c) to asymptotically estimate the first term on the right side of (6.13) as ' ( ^ T ' * 0 ) - ^ ~ d / \ £ 9 \ x 0 + ey)U[(y)u:(y)dy. (6.15) Expanding g in a Taylor series, and using the fact that uc(y) is even, we integrate by parts to get (KXU" \ e^ e-'/'Ae2^1 _ . [°° / K 3=0 w ,2^ + 1 rOO , ,2 £ ^-)T5{2j+2)(*o)/%, P3 = J ^ v2j [«c(v)] <*y • (6-16) For the exponentially small perturbation K — 1 to have a significant effect on the eigenvalue for at least some values of XQ we must balance the exponential orders of (6.14b) and (6.16). Thus, for e —> 0 we take d to satisfy 0<d<2-celoge, (6.17) 144 for some c independent of e. With this restriction on d, the last term on the right side of (6.13) is estimated for e -> 0 as rqe b (s + i) ^<~ldx) (f KuPcu'c<te) ^ O ^ i n a x l e - ^ ' 1 ^ 0 ) ^ - ^ ) 6 - ^ 1 - ^ ) } ) , (6.18) for some q. Since r > 0, we conclude that the term in (6.18) is asymptotically smaller as e —> 0 than the inner product term in (6.16). Finally, substituting (6.14) and (6.16) into (6.13), and neglecting the last term on the right side of (6.13), we obtain the following main result for Ao: Proposition 6.1 (Exponentially Small Eigenvalue) Let (2.13a) and (6.17) be satisfied. Then, for e —> 0 the eigenvalue Ao of (6.7) with the largest real part is exponentially small and it has the asymptotic estimate Ao = A0(xo) ~ ^ ( e - * - 1 ^ ) + e - ^ ( i - * o ) ) _ <^!i £ ^ . g ^ ) { x ^ . . ( 6 . 1 9 ) Here a is defined in (6.4) and 3j for j> 0 is defined in (6.16). From the calculations above we observe that the nonlocal term in Le, which is not self-adjoint, does not influence the leading order asymptotic estimate for the exponentially small eigenvalue. Hence, the adjoint operator L* also has an exponentially small eigenvalue with the same estimate as that given in (6.19). 6.1.2 The Metastable Spike Motion We now derive an ODE for xo(t) from (6.6). Assume that w(x, 0) = 0 in (6.6) so that the initial condition for (6.1) is a spike-layer solution with spike center initially located at xo(0) = XQ € (—1,1).. Since the principal eigenvalue of Le is exponentially small we can assume that the motion is quasi-steady by neglecting dtw in (6.6). To ensure that w is uniformly small over exponentially long time intervals, we must impose the condition that w is orthogonal to 4>Q. To derive an ODE for xo{t) we begin by using (6.11) with ij) = 4>Q and w = <f> to get (L€w, 0O)K = e 2k(?Wx|li + KL*cfr 0)K . (6.20) 145 Using the remark following (6.19) above, we get (w,L*cj)o) ~ A 0 (w,</>Q) = 0 by our condition of orthogonality. From this condition, and by using (6.3a) and (6.6), (6.20) reduces to the following implicit differential equation for xo(t): i (00, dtuc)K + e^ ou'dlx ~ e <fa] . (6.21) K / K The two terms on the left side of (6.21) can be estimated as in §2 by using (6.10) for </>o and (6.4) for u"c{y) as y -» ±oo. For e —>• 0, we get ((po,dtuc)K x080, (6.22a) -e«^o« ' c |L i ~ 2ea 2 (e~2^l+x^ - g ^ d - s o ) ) . ( 6 . 2 2 b ) Next, for e -» 0, we calculate as in (6.15) that ( ' \ poo o 0 0 2 j A> j ~ £»"e-"< j J (xo + ey) [u'c(y)} dy ~ e2^e~^ g ^Jy9W+1) (*<>)/% • (6.23) Finally, substituting (6.22) and (6.23) into (6.21), we obtain the following metastability result for xo(t): Proposition 6.2 (Metastability) Let (2.13a) and (6.17) be satisfied. Then, for e -> 0, the trajectory Xo(t) of the center of the spike for a one-spike solution to (6.1) satisfies the asymptotic differential equation Here a. is defined in (6.4) and 8j for j > 0 is defined in (6.16). A n important remark, as observed from (6.24), is that the behavior of xo(t) depends only on pointwise values of certain derivatives of the perturbation K — 1. This will be different from the behavior that we wil l observe in §6.2. 146 6.1.3 A n Example of the Theory We now discuss the qualitative behavior associated with (6.19) and (6.24). From (6.24) we see that h(—l) < 0 and h(l) > 0 as e —> 0 when d > 0. Thus, there exists at least one equilibrium value XQ for XQ (t). The existence of any other equilibria for XQ depends on the constants d and v and the function g'{x). In particular, when d > 0 is sufficiently small, then (6.24) has an equilibrium point near each zero of g (x). As shown in the example below, we can have other equilibrium points in the interval [—1,1] when d is near some critical value. The next result concerns the stability of the equilibria of (6.24). Let XQ satisfy h(xQ)=0. Then, by comparing (6.19) and (6.24), we find that h'(x%) = 2XO(XQ) . This shows that the decay rate for the differential equation (6.24) associated with infinitesimal perturbations about XQ is 2 A O ( X Q ) . This leads to the next result. Corollary 6.1 (Stability of Equilibrium) Let x% satisfy H(XQ) = 0. Then (6.1) has a one-spike equilibrium solution of the form given in (6.3) and this solution is stable (unstable) if \O{XQ) < 0 (Ao(X(j) > 0) • Here uc(y), Ao(xo) and h(xo) are given in (6.4) , (6.19) and (6.24). Thus a one-spike equilibrium solution centered at x\ is unstable when ( / ' ( X Q ) < 0 • Since g"(x) < 0 corresponds to a weakly convex domain, this result predicts that there is no stable spike-layer solutions in such domains. However, when g"(xo) > 0, then A O ( X Q ) c a n be negative for certain choices of v and d, resulting in a stable spike-layer solution centered at x% . The key point to construct such a stable equilibrium solution is to guarantee that (6.24) has multiple equilibria corresponding to simple zeroes of h(xo). Then, we must have exactly one stable equilibrium of (6.24) between every two consecutive unstable equilibria. To illustrate the result, let p = 2 for which uc(y) = (3/2) sech2 (y/2). Then, we calculate from (6.4) and (6.14a) that 0Q = 6/5 and a = 6. Thus, to leading order (6.19) becomes Ao = Ao(xo) ~ 60 (e-*- 1^) + c-fc-1d-*o)) _ ^ V ^ 1 dg"(xQ) + •••, (6.25) 147 and (6.24) reduces to x0 ~ h(x0) = 60e ( e - ^ ^ d - x o ) _ E-2e-i(i+*o)) _ ^ V ^ V (*o) + . . . . ( 6 2 6 ) Choose #(:r) = x 2 / 2 , which corresponds to a non-convex domain. In this case, XQ = 0 is an equilibrium solution to (6.26) for any v and d. This solution can be stable if d is small enough. Prom (6.25), we calculate A0(0) = -^e~e~ld + 120e- 2 e _ 1 . Let dc be the zero of A0(0) as a function of d. Then, dc = 2 - e log 240 + (v + 2)eloge. From Corollary 6.1 it follows that the equilibrium x% = 0 is stable (unstable) when d < dc(d > dc). When d < dc there are two other equilibrium points for (6.26) on either side of x = 0 that are unstable. This example clearly shows the effect of pinning whereby the exponentially weak non-convex perturbation of the original domain leads to a stable one-spike equilibrium solution to the shadow problem (6.1) centered at the midpoint of the domain. 6.2 One-Spike Dynamics: The Perturbed Gierer-Meinhardt Model We now analyze the effects of spatially varying terms on the full system (1.17). We first examine the effects of a spatially varying inhibitor decay rate \i = u(x). We consider two cases. If D = 0(1) we repeat the process used to derive the differential-equation governing the motion of the spike from §5.1. The resulting differential equation depends on the global properties of u(x) and clearly demonstrates the qualitative effects of the spatially varying coefficient. However, as demonstrated in §4.3, for n-spike solutions to be stable we require that D to be small. Thus we also examine the effect of a spatially varying inhibitor decay rate when D is small. In this case we must use W K B theory to solve the inhibitor equation, (more here). Finally we examine the effect of a spatially varying activator decay rate. 6.2.1 A Spatially Varying Inhibitor Decay Rate When D = O(l) In this section we analyze the dynamics of a one-spike solution to (1.17). For finite inhibitor diffusivity D and for /z = a(x) > 0, we derive a differential equation determining the location xo(t) of the maximum of the activator concentration for a one-spike solution to (1.17). 148 In the inner region near the spike we introduce the new variables y = e - 1 [x - XQ(T)] , h{y) = h(x0 + ey), a(y) = a(x0 + ey), T = e2t. (6.27a) We then expand the inner solution as My) = ho(y) + ehi(y) + --- , d(y) = d0{y) + eai{y) + • • • . (6.27b) The spike location is to satisfy a0(0) = 0. Substituting (6.27) into (1.17), and collecting terms that are 0(1) as e -> 0, we get the leading order problem for ao and JIQ: a0' — do + cft/hl = 0, -oo < y < oo , (6.28a) h'o = 0, (6.28b) with a0(0) = 0. In order to match to the outer solution, to be constructed below, we require that ho is independent of y. Thus, we set ho = H, where H. = H(T) is a function to be determined. We then write the solution to (6.28a) as do = JTuc, where 7 = q/(p-l), (6.29) where uc satisfies (6.4). Collecting the 0(e) terms in the inner region expansion, we obtain the problem for d\ and hi, ~p—i ~p a-i - ai + yq oi = ~^+ihi ~ xo aOy > -co < y < oo, (6.30a) Dhiyy = -dro/hs0 . (6.30b) Here x'0 = dxo/dr. Next, we write ai as di = Wui. (6.31) Using (6.29), (6.31), and ho = H, (6.30) becomes P L(ui) = Ui - ui +pu% ui = —hi-xQuc, - c o < y < o o , (6.32a) Dhiyy = , (6.32b) 149 where u\ is to decay exponentially as |y| —» co. Since Lu'c = 0 and u'c -> 0 exponentially as \y\ -» co, the right side of (6.32a) must satisfy the solvability condition that it is orthogonal to u'c. From this condition we obtain the differential equation Xn — /O O upcu'ch •oo idy. (6.33) H I-ooiu'c(y)}2 dyJ-co If we integrate (6.33) by parts twice, and use the facts that h"x and uc are even functions, we get q ('HoM^r1 dy 2H(p + l) \ f^Kiy)}2 dy lim h.x + lim hx y—>+oo y—>—oo (6.34) In the outer region defined at an 0(1) distance away from the center of the spike, a is expo-nentially small and we expand h as h = ho[x) + o(l) as e —> 0. Then, from (1.17b), we obtain that ho satisfies Dh0 - /j,h0 = -H1T sbr5{x - x0), -1 < x < 1, h'0(±l) = 0. Here br is defined in (6.3b). Solving for ho we get h0(x) = H"<r-sbrG(x; x0), where the Green's function G(x; XQ) satisfies DGXX — uG = —5(x — xo), — 1 < x < 1, G s ( ± l ; x 0 ) = 0 . To match with the inner solution we require that hoixo) = H, lim hi + lim h, = h0(xo+) + hQ(xo-) • y—>-+oo y—>—oo Substituting (6.36) into (6.38), we get lim hl + lim Ai1 H = H y^+oo j/->-oo G(XQ;XQ) 1 [Gx(XQ+;X0) + Gx{x0-;x0)] , l / [ 7 r-(s+l)] brG{x0;xo) (6.35a) (6.35b) (6.36) (6.37a) (6.37b) (6.38) (6.39a) (6.39b) 150 Finally, substituting (6.39) into (6.29), (6.34) and (6.36) and letting r = e2t, we obtain the main result of this section: Proposition 6.3 For e <C 1, the the dynamics of a one-spike solution to (1-17) is characterized by a{x,t) ~ Wuc (e _ 1[s - x0(t)]) , (6.40a) h(x, t) ~ HG [x; x0{t)] jG [x0{t); x0(t)], (6.40b) where H = H(t) is given in (6.39b). The spike location xo(t) satisfies the differential equation dxp ^ _e2Q (GX(XQ+;X0) + Gx(XQ-;XQ) dt \ G(x0;xo) where C > 0 is defined by J-F^i' ( M O d ) In calculating the integral in (6.40d) we used (5.19) 6.2.2 Case 1: n(x) > 0 depends on x when D large In general, when p, depends on x we must compute the Green's function satisfying (6.37) to determine the dynamics as described in (6.40c). However, to illustrate qualitatively the effect of a spatially varying p(x), we now derive an approximate differential equation for XQ in the limit D > 1 with D independent of e. In the limit D S> 1, we expand G as G{x; x0) = G0(x; x0) + D^G^x; x0) + O (D~2) . (6.41) Substituting (6.41) into (6.37) and collecting powers of D~ l, we get GoXx = 0; Gixx = pGo - 5{x - XQ) , (6.42) (6.40c) 151 with Gjxx — 0 at x = ±1 for j = 0,1. The problem for G\ does not have a solution unless Go satisfies a solvability condition. In this way, we calculate that Go = {2par1; Glx = (2aa)~1 ^ a(y) dy - j ° 0 —1 < x < xo , (6.43) XQ < x < 1 . Here / j a is the average of /j over the interval, defined by 1 rl M a 5 5 2 j x ^ dX' ^6'44^ Substituting (6.43) into (6.40c) we obtain the following result: Corollary 6.2 For e < 1 and J) » 1, with D independent of e, the differential equation for the spike location (6.40c) reduces to dx0 2elC ~dt ~ iJ where C is defined in (6.40d). (J\(y)dy-^ , (6.45) From (6.45) we observe that the pinning effect induced by /J(X) depends on global properties of the spatial inhomogeneity /J(X), in contrast to the pointwise values as obtained in §6.1 for perturbations of the shadow problem. Since /J(X) > 0, there is a unique equilibrium spike-layer location xoe for (6.45) satisfying clOe I. Ii(y) dy = ua. (6.46) This equilibrium is a stable fixed point for (6.45). Notice that if JQ \xdx < J^fidx, then the equilibrium location satisfies xoe € (—1,0). Alternatively, if there is more mass of a on the right side of x = 0, then xoe E (0,1). As an example, let u > 0 and consider the profile "<*H(i+;£)- (M7) It is easy to see that aa = 1 for any LJ > 0. Also, as u -> 0 we have /J(X) -> 1. As u' -> co we have //(x) -t 1/2 + ue~w(~x+1\ which has a boundary layer near x = - 1 . For any u > 0 there 152 is more mass of /J to the left of x = 0 than to the right of x = 0. From (6.46), xoe satisfies the algebraic equation x0e - —. r = 1 - ~ r • (6.48) sinhw smhw It is easily seen that —1 < xqe < 0 when <j > 0 and XQe -» — 1 as u> -» oo. 6.2.3 Case 2: /^(x) > Q depends on x when D is smal l We now examine the effect of setting p = p(x), in (1.17), to the dynamics of a one-spike solution when D is small. In the previous subsection we considered a spatially varying inhibitor decay rate in the limit D —> oo. The motivation for considering the limit D -> 0 is that multi-spike solutions are only stable when D is sufficiently small. Although here we only consider the one-spike case when D is small, a similar analysis could be done for a multi-spike solution. We now solve (6.40c) in the limit D 0 using the W K B method. Since D is small, we make the following ansatz as to the form of G, G{x; x 0) = exp (??j=l + Qx (x) + • • •) . (6.49) Substituting (6.49) into (6.37) and collecting powers of D gives us the following set of equations, e'i = H , (6.50a) 0o' + 20 o 0i=O. (6.50b) This leads to the following leading order solution for G, G(I.Io)=f^"1/4OTK^i>1/2*)cosh(*/-i''1/2«''')' -'<«<*•• (M1) | ^ - V 4 c o s h (_>_ / » pl/ljy) c o s h £ ^ 1 / 2 ^ x „ < X < 1 . To solve for A we use the jump condition associated with (6.37), Gx{x0+;x0) - Gx{x0-;x0) = - - ^ . (6.52) Applying the condition above, results in the following value for A: 153 where ^ = \f\M(y)]1/2dy. (6.54) Substituting (6.51), (6.53) and (6.54) into (6.40c gives us the following differential equation for 4 / z - 1 f * n W f i : n V r n s W J - i i . i 7 2 - ' \ 4-mh(-LM*1/2^ -I- titlt ainh(-L.(1*1/2^ dxo _ - |M- 1 (^o)^(xo)(cosh(j- MV2) + co8h( j - / Q ) + ^ s i n h ( ^ X o / 2 ) (6.55) d T c o s h ( ^ m J i / 2 ) + c o s h ( ^ ^ / 2 ) where / i * / 2 is defined in (6.54 and a*xlJ2 is given by, ltl\y)dy- / A* 1 / 2 (y)* / - (6-56) -1 ./xo 6.2.4 A variable activator decay rate We now consider the case where the activator decay rate varies spatially and the inhibitor decay rate is a constant. The effect of adding spatial variation to the activator decay rate is to change both the equilibrium location of the spike and its dynamical behaviour. The Gierer-Meinhardt system with a spatially varying activator decay rate is given by, at = e2axx - V(x)a + T- , -1 < a; < 1, i > 0, (6.57a) aT 0 = Dhxx -uh + — , -Kx <l, t>0 (6.57b) hs o x ( ± l ) = hx(±l) = 0. (6.57c) We now extend the analysis given is §4.4.1 to derive an ordinary differential equation for the location of the spike for a one-spike solution. As we are examining a one spike solution, with the spike centered at x = XQ, we make the following change of variables, y = (x- x0)e~l, a(y) = a(x0 + ey), h(y) = h(x0 + ey), (6.58) where XQ = XQ(T) with T = e2t. We then expand the inner solution as follows, Hv) = M y ) + ehi + ••• , a(y) = a0(y) + ea,i{y) + ... . (6.59) 154 The spike location is to satisfy a'0(0) = 0. Substituting (6.58) and (6.59) into (6.57), and collecting terms that are 0(1) as e —> 0, we get the leading order problem for ao and JIQ: a0 - V(XQ)CLQ + a^/hl = 0, -oo < y < co, (6.60a) h'Q' = 0, (6.60b) with a0(0) = 0. In order to match to the outer solution, to be constructed below, we require that ho is independent of y. Thus, we set ho = H, where i f = H(T) is a function to be determined. We then write the solution to (6.60a) as ^ = F [ 7 ( i f l ) ] 1 / ^ \ ( v / V ? W y ) . , where 7 = a / ( p - l ) , (6.61) where uc satisfies (6.4). Collecting the O(e) terms in the inner region expansion, we obtain the problem for a\ and hi, ~p—i ~p a[ - V{xo)a\ + Pa~g ox = yV'(x0)ao + ^ r ^ i _ ^Oy , -co < y < co, (6.62a) K K Dhiyy = -fio/fto • (6.62b) Here x'0 = dxo/dr. Next, we write ax as di = . (6.63) Using (6.61), (6.63), and ho = H, (6.62) becomes L(m) = yV'(xo)V(xo)^uc(VV{x7)y) + ^ x ^ ~ l ) u h i ti - x0[V(xo)f+p)Mp-l)u'c, -oo < y < oo, (6.64a) Dhiyy = -H1T-Surc , (6.64b)' where ui is to decay exponentially as |y| -> co and L(ui) = u[ - V(xo)tii +pV{x0)up-1ui. (6.65) Since, Lu'c{y/V(xo)y) = 0 and u'c-*Q exponentially as |y| -> oo, the right side of (6.64a) must satisfy the solvability condition that it is orthogonal to uc(^/V(xo)y). From this condition we 155 obtain the differential equation qV(x0) XQ — /OO u?(y/V(^y)uc(^vW)y)hi dy -oo Hf-oolu'M] dyJ-oo + T /, A 12 /" uc{y)u'c(y)ydy. (6.66a) If we repeat the process from §6.4.2, we then arrive at the following differential equation for xo, dxo e2q . (Gx(X0+;XQ) + Gx{xQ-;x0)\ e2 (P + %\ V(XQ) ,G 6 R ) dt ~ p - i v G(x 0 ;x 0 ) 7 A P - I J ^ O ) ' L ' J If the second term on the right hand side of (6.67) were absent, the spike would tend towards xo = 0. If the first term on the right hand side of (6.67) were absent, the spike would tend towards a local minimum of V. The end result of the spatially inhomogeneous term is a competition between these two terms. 156 Chapter 7 Conclusions In this chapter we give a brief review of the main results found in this thesis and an intuitive explanation of these results. The equations and details of these results will be omitted as these details are summarized in §1.4. We will also attempt to provide a framework to place these re-sults in a useful context for those who use the Gierer-Meinhardt equations for modeling. Finally, we will discuss some possible extensions to the analysis as well as some possible connections to other systems. Before we proceed with the discussion, we must clear up some notation. In this chapter, when we refer to an n-spike profile, we are referring to a solution of the form (2.1) for the shadow system and a solution of the form (4.14) for the full system. 7.1 Overview of Results We will begin by reviewing the results from §2 and §3. These chapters analyzed the reduced system (1.19) commonly referred to as the shadow system. By constructing spike-type solutions and linearizing about these solutions we find that any solution with a spike in the interior of the domain is unstable. Profiles with more than one spike are unstable with an 0(1) eigenvalue and one-spike profiles are unstable with an exponentially small eigenvalue. The presence of an exponentially small eigenvalue is used to construct an equation of motion for the spike which demonstrates that an internal spike will move towards the closest point on the boundary on an exponentially slow time scale. Once the spike begins to approach the boundary, the analysis becomes invalid. The spike is presumed to merge with the boundary. By linearizing about a spike confined to the boundary, we find that a spike confined to boundary will move in the direction of increasing boundary curvature in two dimensional domains and increasing 157 mean curvature in three dimensional domains. The analysis of the dynamics and stability of these spike profiles for (1.19), all rely on the interaction of the exponentially small tails of the spikes with the boundary. The results are thus very delicate. Inconsistencies with the shadow system and numerical calculations of the full system (1.17) suggest that this delicate structure is broken by finite values of inhibitor diffusivity. This limits the applicability of the results in these chapters to modeling of the full system. However, the methods used in this chapter may be of use in analyzing models other systems, specifically the microwave heating model presented in [4]. The results in this section are also interesting from a mathematical point of view. Specifically, the formal asymptotics in §2 have been verified vigorously in [7] using invariant manifold techniques. In §3 and §4 we examine the stability and dynamics of spike profiles under the full system (1.17) with finite values of inhibitor diffusivity. The results from (1.19) do not coincide with numerical calculations of the full system (1.17) even for large values of inhibitor diffusivity. It is for this reason that we undertake the more difficult task of analyzing the full system. Again we construct an n-spike profile for the full system and linearize about this solution. Analyzing the spectrum of the operator resulting from a linearization about an n-spike profile results in a necessary and sufficient condition for the stability of any given profile. We find that there is a decreasing sequence of numbers Di such that an n-spike solution is stable if the inhibitor diffusivity is less than Dn. The existence of an algebraically small eigenvalue is used to construct a solvability condition which results in a differential equation governing the motion of the spikes. This equation reveals that the shape of the inhibitor profile overrides the effects of the exponentially small tail in the activator profile and takes control of the dynamics. This is the main reason for discrepancies between the behavior of the numerical calculations of the full system and the shadow system. While these results are all in one spatial dimension, and most models are in two or three dimensions, we believe that similar qualitative behavior exist in higher dimension. Again these results are also interesting from a purely mathematical point of view and the formal asymptotics of §4 have been vigorously verified in [51]. In §6 we consider the effects of spatially varying coefficients. Three different scenarios are 158 considered, a spatially varying activator diffusivity for the shadow system, a spatially varying inhibitor decay rate and a spatially varying activator decay rate. In all three cases the effect of the spatially varying coefficient is to add an extra term to the solvability condition which provides the deferential equation governing the motion of the spike. This chapter may be of particular interest to those who use the Gierer-Meinhardt equation for modeling as many of the models use spatially varying coefficients to account for preexisting patterns. 7.2 Possible Extensions There are still many open ares of research for the Gierer-Meinhardt model. The most obvious being the extension of the results in §4, §5 and §6 from one to higher dimensions. The main difficulty to be overcome for this extension is the complexity and domain reliance on the Green's function for the inhibitor equation. Another important extension is to include the effects of non-zero values of r. Numerical computations suggest that for reasonable small values of r, the simplification of setting r to 0 presents no problems. However, for larger values of T , many interesting phenomena, such as oscillatory behavior, are possible. For examples of such behavior see [31]. The methodology used in this thesis may also be extended to other systems. At present the analysis in §4 is being modified to analyse the Schnakenberg model, another reaction diffusion system used in developmental biology modeling. The analysis in §2 may also be of use in analyzing the microwave heating model presented in [4]. 159 Appendix A Proof of Theorem 2.1 In this appendix, we prove Theorem 2.1. Although this has been proved in by J. Wei in [52], we include a proof here for the convenience of the readers. We present a proof which works in the general case of M.N. Let 4 N + 2 r = 2, 1< P <1 + — • or r=p+l, l<p<{^—^)+, where = if N > 3 and = +co if N = 1,2. Define w(\y\), with y = (yx.-.i/jv)*, to be the unique positive solution to " N-l , _ n • w H w —w + wp = 0, p>0, P u/(0)>0, to(p) ~ a p ( 1 - i V ) / 2 e _ p , as p - » o o . When iV = 1, then w = uc, where uc satisfies (3.6). Suppose that (</>, Ao), with Ao ^ 0, satisfies the following eigenvalue problem: - <f> + pvf~V - 70(P ~ 1) % ^ ^ t i f = , 4>£H2{R}), 70 >1- (A.l) When N = 1 this problem reduces to (3.39). Thus, the proof of (2.13) is complete once we show that Ee(A0)<0. (A.2) Let Ao = A/j -+- 4> = 4>R + We first introduce some notations and make some preparations. Set L 0 : = L O ^ - 7 O ( P - 1 ) % ^ ^ , <t>£H2(&), where 70 > 1 and Lo := A - 1 + pwp~l. Note that L is not self-adjoint if r ^ p + 1. It is well-known that Lo admits the following set of eigenvalues: Pi > 0, p2 = 0,...uN+i = 0, fj,N+2<0, (A.3) 160 where the eigenfunction corresponding to p\ is of constant sign (see Theorem 2.1 of [27]). Let duo' XQ := kernel(Lo) = span{-—,j = 1, ...,N}. Then and LQW = (p - l)wp , LQ(——-w + \-xVw) = w, (A.4) [ {LQ1W)W = f w(-^—w + \xVw) = (-L- - y) /" ™2> (A-5) JR* JR" p - 1 2 p - 1 4 y j^v / {LQ1W)WP = f wp(—^—-w + \xVw) JUN JM." p - 1 2 = / (LQ1W)-^—L0W = -^—[ W2. , (A.6) yRjv p - i p - i y R ^ r We divide our proof into three cases: Case 1: r = 2 ,1< p < 1 + ± . Since L is not self-adjoint, we introduce a new operator as follows: := - (p - 1 ) 4 ^ - (p - + (P - l ) i k = ^ % ^ w . (A.7) V ™ 2 V™ 2 ( V ^ 2 ) 2 We have the following important lemma: Lemma A . l (1) L \ is self-adjoint and the kernel of L \ (denoted by X\) = span {w, = 1,...,N}. (2) There exists a positive constant a\ > 0 such that :=r {m2+? - P Wp-v) + 2 { p~ i ) fr t v wP(f> - (p - DT^CC / •/R" V™ 2 O R " ™ 2 ) 2 . / R " > aid2L2^N}{4>,Xi), for all <f> G Hl(RN), where d/J2(RAR) denotes distance in the L2-norm. 161 Proof: By (A.7), L\ is self-adjoint. Next we compute the kernel of L\. It is easy to see that w, %yj,j = 1 , N , £ kernel(Li). On the other hand, if ^ G kernel(Li), then by (A.4) LQ4> = ci(<p)w + C2(<P)VJp = ci(4>)L0(—^—-w + ^ V w ) + c 2(</0Lo(—r). p — 1 2 p — 1 where Hence Note that f „„2 w 77 2 V 2 ' — \f ± r i (f> - cx(<f>)( * it; + Lvio) - c2((f>) 1 w € kernel{L0). (A.8) ci(0) = (p - l)ci(^) r — 5 ( p - l) C l (0) ' = C 1 W - , 1 W ) ( ? 3 T - T ) ^ ^ by (A.5) and (A.6). This implies that ci(</>) = 0. By (A.8), this proves (1). It remains to prove (2). Suppose (2) is not true, then by (1) there exists (a, (j>) such that (i) a is real and positive, (ii) <f) ±w, 4> -1 j = I,... ,N, and (iii) L\<$> = a<j>. We show that this is impossible. Prom (ii) and (iii), we have ( L o - a ) ^ = ( p - l ) A ^ t i ; . (A.9) We first claim that fRN wp4> / 0. In fact if JRN wp(f> = 0, then a > 0 is an eigenvalue of LQ. But by (A.3), a = and (f> has constant sign. This contradicts with the fact that </> ± w. Therefore a y£ HI,Q, and hence Lo — a is invertible in XQ. SO (A.9) implies ^ = ( p - l ) # ^ ? ( L o - a ) - 1 u , : Thus VR" JRJV W JR" / w2 = (p-l) [ ( ( L o ' - a J - ^ K , 162 / w2 = / ((Lo —a) lw)((Lo — a)w + aw), JRN • JRN 0 = f ((L0 - a)~1w)w. (A.10) Let hi(a) = fRN((L0 - a)~lw)w. Then, hx(0) = JRN(LQ1W)W = fRN(^w + \x • Vw)w = (p=l ~ f ) IR"™2 > 0 s i n c e 1 < P < 1 + W- Moreover h[(a) = JRN((Lo - a)~2w)w = fRN((Lo — a)~lw)2 > 0. This implies h\(a) > 0 for all a G (0,/zi). Clearly, also hx(a) < 0 for a G ( /JI,CO) (since lima_>+00 h\(a) = 0). A contradiction to (A.10)! This completes this part of the proof. We now finish the proof of (2.13) in Case 1. Since Ao 7^ 0, we can choose (p J_ kernel(Lo)- Then we obtain two equations L0<t>R - (p - 1 ) 7 0 % ^ ^ = A « < ^ ~ A'<^> ( A - n ) Lo4>i - (p - 1)70 V Wirwp = XR<P! + \KPR. (A.12) JRN W Multiplying (A.ll) by <J>R and (A.12) by (pi, and adding the resulting expressions together, we obtain (<t>R + $) = LI((J)R,4>R) + Li(<pi,<pi) +(p-l)(70-2) JT—-5 JRW «> (JlBAT , / V /R" (/R^ 2) 2 L V Multiplying (A.ll) by u; and (A.12) by u; we obtain (p-1)/" wP^R-joip-l)^^ [ vf+l = XR [ w<pR-Xi[ w<t>i, (A.13) . / V JKiv ^ JR" , / V • jRN 1 (p-1)/" f^-7o(p-l)V",tf / = XR f w<pi + Xj j wcpR. (A.14) V R " J R A T ™ . / R " JRN JR" Multiplying (A.13) by fRN w<pR and (A.14) by jRN wcpi, and adding them together, we obtain (p-1) W(f)R / WP(f)R + (p - 1) / uxpi \ wp<pi JRN JRN JRN JRN 163 JRN WZ JRN JRN (4>R + <fi) = Lx (4>R, 4>R) + Li (fa, fa) Therefore we have +(p - l ) ( 7 0 - 2)(-I-A* + 7 0 % L ^ ) ( V - ^ ) 2 + + ( p - l ) p-1 fRNw2 fRNw2 j f ^ i i f ^R)2 + ([ Wfa)2}. UR"W) JR" JR" Set Then <!>R = CRW + <pR-, <J)R- 1 Xi, fa = ciw + <j>r, (/>i ± Xi. / w4>R = cR w2, / wfa = ci w2, JRN JRN JRN JRN ^ ( R A o O f o j - X i ) = II^IIL2> dL2(RN){<t>i,Xi) = .\\<j>i\\2L2. By some simple computations we have LI{4>RAR) + LI(4>I,^I) +(70 - 1 ) A * ( 4 + cj) f w2 + (p- 1)( 7 0 - 1?(C2R + cf) f v?+1 + A*( | |4 | | 2 2 + \\<t>j-\\\.) = 0. By Lemma A . l (2) ( 7 & - 1 ) A H ( 4 + ^ ) / « / R ' VJ2 + ( p - l ) ( 7 o - l ) 2 ( 4 + c2-) /* ^ ^ ( A f l + aOdl^ l^ + l l^ l l iO^O. ./RJV Since 70 > 1, we must have XR < 0, which completes the proof of (2.13) in Case 1. Case 2: r = 2,p = 1 + ^ . 164 In this case we have / (LQ1W)W= I w(—^—w + l-xWw) = 0. (A.15) JRN JRN p - 1 I Set wo = ——-w + 1-xVw. (A.16) p — 1 2 We will follow the proof in Case 1. We just need to take care of WQ- We first have the following lemma which is similar to Lemma A.l: The proof is omitted. Lemma A.2 (1) The kernel of L\ is given by X\ = span {w,wo, = 1, —,N}. (2) There exists a positive constant a 2 > 0 such that Ll(<fi,<f>)= [ ( I V ^ + ^ - p ^ - V 2 ) JRN + 2(p - 1) / „ J R , w^ / ^)2 IRNW2 ( V ™ 2 ) 2 JRN > a 2 ^ a ( R W ) ( 0 , X i ) , V ^ 5 ' ( R N ) . Now we can finish the proof of (2.13) in Case 2. Similar to Case 1, we obtain two equations (A.ll) and (A.12). We now decompose (f>R = cRw + bRwo + <t>R, <f>R±Xi, 4>i = aw + biw0 + 4>f-, (pf LXX. Similar to Case 1, we obtain Li{(pR: <PR) + Li ((pi, 4>i) +(7o - 1)A«(4 + cf) / u-2 + (p - 1)(70 - l)2(c2, + cf) / ./R^ JRN • +\R(bR([ wif + bjif ^o2)2 + ll4ll 2 i 2 + ll^lli 2)<o. JRN JM.N By Lemma A.2 (2) (70 - l)Aa(4 + cf) / w2 + (p - l)( 7 o - 1)2(4 + cf) / yR^ yRjv +A*(4(/ w2o)2 + b}([ «;o2)2) + (Afi + a 2)(||4|| 2 2 + ||^|| 2 2)<0. V R ^ JRN. 165 If XR > 0, then necessarily we have CR = a = 0, 4 = 0, ^/ = 0. Hence 4>R = bRWo, <p[i = biWQ. This implies that bRLoWQ = (bR - bi)wQ , biL0wQ = (6R -+- 6/)u;o, which is impossible unless bR = 6/ = 0. A contradiction! This completes this part of the proof. Case 3: r = p + 1 , 1 < p < ( # ± § ) + . Let r = p + 1. L becomes We will follow the proof of Case 1. We need to define another operator. L ^ ^ L ^ - i p - l ) ^ ^ - ^ . - (A.17) We have the following lemma: Lemma A.3 (1) L% is self-adjoint and the kernel of L$ (denoted by X%) consists ofw, J ^ - , j ' = l , . . . , iV . (2) There exists a positive constant a$ > 0 such that JRN V ^ P + >a3d2LHRN)(<j>,X3), V<f>eHl(RN). Proof: The proof of (1) is similar to that of Lemma A . l . We omit the details. It remains to prove (2). Suppose (2) is not true, then by (1) there exists (A, (p) such that (i) A is real and positive, (ii) cj> ± w, <f> J_ |^-, j = 1,... , JV, and (iii) L3(p = Xcp. We show that this is impossible. From (ii) and (iii), we have fr-*>»-V&"'V (A-18) Similar to the proof of Lemma A . l , we have that J R N wv<p / 0, A / / / i , 0, and hence LQ — A is invertible in XQ-. So (A.18) implies 166 Thus • f ^ = ( p - l ) i ^ / ( d o - A ) - V K JRN JR" W V + JR» [ W P + 1 = (p - 1) /" ( ( L 0 - A ) - 1 ^ ) ^ . (A.19) Let / J 3 ( A ) =• (p - 1) V ( ( L 0 - A ) " 1 ^ ) ^ - / R „ then h3(0) = (p - 1) / R j » ( £ 0 - V ) ^ -=0. Moreover / » ' 3 ( A ) = (p- 1) fRN((L0 - X)-2wp)vf = (p - 1) J R J V ( ( L 0 - A ) " 1 ^ ) 2 > 0. This implies / i 3 ( A ) > 0 for all A G (0,^ :). Clearly, also / i 3 ( A ) < 0 for A G (p-i,oo). A contradiction to (A.19)! This completes this part of the proof. We now finish the proof of (2.13) in Case 3. Similar to case 1, we obtain two equations LocpR - (p - 1)70 W^PtR1™p = XR<pR - Xjfa, . (A.20) RN Lo<fii -ip- 1)70 {RJV WPl[wp = \R<f>i + \I4>R. (A.21) Multiplying (A.20) by CJ>R and (A.21) by 4>i and adding them together, we obtain ~XR I i<l>R + 4>i)=L3{(l>R,<f>R) + L3((i>T,<pI) + ( p - i ) ( 7 o - i ) ( ^ ^ f + ( / r ^ j ) 2 -By Lemma A.3 (2) A, / (02 + + + (p - l)(7o - D ^ ^ + U H * ^ ) 8 < o, J*S JRNWP+1 which implies XR < 0 since 70 > 1. Thus, (2.13) in Case 3 is proved. 167 Appendix B The Laplacian in the Boundary Layer Coordinate System The derivation here is similar to that in [41]. We begin with a description of the boundary. Let z = H(pi,p2) define the local height of the boundary at the point (pi,p2) on the boundary. For convenience here we will center the coordinate system around the point p\ = 0 and p2 = 0, where z = 0. In terms of these coordinates, H is given locally by H{p\,P2) = + ^K2p22 + O (pl+p22) , (B.l) where K I and K2 are the two principal curvatures at the center of the coordinate system. Our boundary layer coordinate system is (pi,P2,v)y where n > 0 is the distance from x G O to dVL. Therefore, we define the following change of coordinates, x(pi,P2,»?.) = (pi,P2,H{sx,s2)) + r)n(pi,p2), where n(pi,p2) is the unit normal to the boundary defined by, (-HPl,-HP2,1) (B.2) n(pi,P2) - V1 + Hl + Hh In order to find the Laplacian in these new coordinate, we must calculate the scale factors (B.3) dx the the VP1 = dpi 1 UP2 ~ dp2 dn (B.4) To evaluate these expressions, we substitute (B.l) and (B.3) into (B.2) and differentiate to find, un = l- nm + 0{p\ +p22 + n2), l-nK2 + 0(p\+p22 + n2), In general curvilinear coordinates, the Laplacian is given by, A<j> Vp1 Vp2 dpi 'Pi " T J vn d4> + d (vVxvn d<p uPl dpi) dp2 \ uP2 dp2) dn + d<j> dn (B.5a) (B.5b) (B.5c) (B.6) 168 Thus, in these variables, the Laplacian is given to within quadratic terms by (1 — r7«x)(l — T 7 / S 2 ) P 2 V l - " K 2 P V This coordinate change is then used in (3.1a) to obtain (3.26). 169 B . l An Asymptotic Estimation of an Inner Product Now we bound the term (£*[dxuc],(j>x) in (3.65), where £* is defined in (3.56). Since dxuc satisfies the local part of the operator in (3.56), we obtain —2 m— 1 r £*e[dxuc] = - m T ( " C 7 vPcdxucd*: (B.8) In §3.5 we assumed that the distance between x = f and the curved part of the boundary dClc is minimized at one of the two corner points. Let rm = min(xR — £, £ — x£) denote this minimum distance. Let BT denote the semi-circle whose diameter is the interval [£ — r m , £ + rm] along the x-axis. Then, by our assumption, Br must be strictly contained within $X We then decompose the integral in (B.8) as / updxucdx.= / updxucdx + I updxucdx.. (B.9) Jo. JBt Jci\BR Since the point x = £, y = 0 is the center of the semi-circle and the integrand is an odd function about the line x = £, the first integral on the right side of (B.9) is identically zero. Next, since uc decays exponentially away from the point x.= £, y = 0, the second integral on the right side of (B.9) is bounded by the maximum of the integrand on the boundary of Q\Br multiplied by the area Ar of Q\Br. In this way, using the far-field behavior (3.14b), we get | f updxucdx\<Ceqe-(p+1)rm/e. (B.10) Jn\BR Here q and C are constants independent of e. Hence, in Q, we have the estimate \C*e[dxUc]\ < Ce«- 2u™- 1e-fr + 1> r™/«, (B.ll) for some new constant C. Since <pi ~ dxuc, we can then use the same reasoning as described above to estimate (£*[dxuc],<p\), where C*[dxuc\ is given in (B. l l ) . We find | (£*[dxuc], <j>i) | < Ce"e-^+m+l>^'. (B.12) for some new constants q and C independent of e. Finally, we compare (B.12) with the asymptotic order of the boundary integral in (3.65). The boundary integral in (3.65) is clearly O (e 9 e - 2 r m / € ) . However, since p > 1 and m > 0, it follows that the inner product term in (B.12) is asymptotically exponentially smaller than the boundary integral term in (3.65). Hence, we were justified in asymptotically neglecting the inner product term in (3.65). 170 Appendix C Calculation of S and V Consider the boundary value problem Dy"-w = 0, y ' ( ± l ) = 0 , [Dyh = Ot Dy -U>j , (C.la) (C.lb) for j = 0,... , n — 1, where [vh = V(XJ+) — V(XJJ) and Xj satisfies (4.1). The solution i sn - l y(x) = ^>2G{x;xk)ujk , (C.2) fc=0 where G satisfies (4.10). Define the n-vectors y and (y) by y' = (yo,.-. , vn - i ) , <y'>*= (<y'>o,..-,(y')n-i) , (C.3) where yj = y(x_?) and (y')j = ^y'(xj +) + y'(a;j_)j /2. Then, we obtain from (C.2) that v = Go, (y') = Vu, - (CA\ where u;' = (UQ, . . . , w n _i). Here the matrices f? and 7> are defined in (4.23) and (4.81), respec-tively. To determine these matrices explicitly we solve (C.l) analytically on each subinterval and impose the continuity of y to get —1 < x < XQ , y(x) = < cosh[fl(l+x)] y° cosh[0(l+xo)] ' sinhW(x,-+i—x)l . sinh[6(x-Xj)] ^ ^ • n n yJOnhMx^-xJl +yj+^[e{x-+1-xi)\ » Xj<X<Xj+l, J = 0 , . . . , n - 2 , coshffl(l-i)] ^ 1 I ^ - l c o s h [ e r i -Xn - l ) ] ' < X < 1 ' (C.5) 1 /2 where 9 = (p/D) 1 . To determine the relationship between y and a>, which yields Q, we use (C.5) and the jump condition Dy. Bv = —<*>, —> y D9 ' Wj in (C.lb) to get 1 Q = —B~l, D9 ' (C.6) 171 where B is defined in (4.26b). Now using (C.5). we can calculate (y) in terms of y in the form where C is defined in (4.87b). Comparing (C.4) and (C.7), and using (C.6), we get the key relation T = - ^ c s c h (24) C^'1 • (C8) 172 C . l Calculation of Matrix Eigenvalues of B In this appendix we calculate the eigenvalues KJ and eigenvectors q^ of the matrix problem Bq = nq, (C.9) where the tridiagonal matrix B is defined in (4.26) and q* = {qi,... ,qn). The calculation below is similar to that given in [22]. From (4.26c) it follows that d = e + f. Therefore, we get the following recursion relation for the coefficients qi of the eigenvector q: fqi-i + (e-K)ql + fql+1 = Q, Z = 2 , . . . ,n -1 , (C.lOa) /gi + ( e - « ) g i + /92 = 0, (C.lOb) fqn + (e - n)qn + fqn-.x = 0. (C.lOc) Hence, to solve for the qi we can use the relation (C.lOa) for I = 1 and I = n and then impose the end conditions 9o=-9i, qn+i = qn- (C.lOd) The solution to (C.lOa) is qi = aCl+ + Ki, C± = jj (* - e ± [(K - ef - 4/ 2] 1 / 2) . (Cll) The end conditions (C.lOd) yield a + b = aC+ + &C_ , (C.12a) aQ + bC = + bC+l • (C.12b) From (C.12) we get ( + = (_ = 1 or = which yields = C-exp(2irij/n), for j = 1,... , n - 1. If (+ = C- = 1 we get K = e + 2/ and ql = (1,... ,1). The other eigenvalues are calculated as in [22] to get KJ = e + 2/ cos (ir(j — l)/n) for j = 2,... , n and £± = exp (±7ri'(j — 1)/n). From (C.12b) we get (1 — £+)a + (1 — C-)6 = 0. Substituting this relation into (Cll), and after rearranging the result, we obtain the unnormalized eigenvectors qij = cos 1 ] (I - 1/2) J , ; = 2,... , n. (C.13) Here qij is the component of the eigenvector q^. These eigenvectors can be normalized and the result is summarized in Proposition 4.2. 173 C.2 Calculation of Bg and Vg Consider the boundary value problem Dy" - ay = 0, (±1) = 0, 3 ' Dy = 0, (C.14a) (C.14b) for j = 0,... ,n — 1, where [u] • = U ( X J + ) - V(XJJ) and x.,- satisfies (4.1). The solution is n - l y ( x ) = Yl9(x'Xk)CJk' (c.i5) fc=0 where g satisfies (4.76). In terms of the matrices Qg and Vg, defined in (4.80) and (4.78), respectively, we have that where a;* = (UQ, ... ,w n _i) . Here y and (y) are defined by /'* = (y'o,• • • ,y'n-i) > <y>* = «y)o,• • •,(y)n-i), y (C.16) (C.17) where y'j = y (XJ) and (y)j = (y(x J +) + y(xj_)) /2. To determine £ 5 and Vg explicitly, we solve (C.14) analytically on each subinterval and impose the continuity of y to get y{x) ( y0 cosh[fl(l+x)] 9 3 inh [ f l ( i + i 0 ) ] ' Vj+\ cosh[fl(x—Xj)} Vj cosh[g( i j +i— x)] 9 sinh[0(xj+ 1— Xj)] 9 sinh[0(xj +i — Xj)\ yn-i cosh[fl(l-x)] 9 sinhlO(l-x n_i)] ' — 1 < X < Xo , Xj < x < Xj+i, j = 0,.. . , n — 2, X n _ i < X < 1 , (C.18) where 0 = {/j,/D)1/2. We then impose the jump condition [DyL = — t o obtain 0 9 D (C.19) where B g has the tridiagonal form given in (4.26b) with matrix entries defined in (4.85). Now we use (C.18) to calculate (y)j, and in this way we get (V) = ^csch ) Cy , (C.20) where C is defined in (4.87b). Substituting (C.19) into (C.20), and comparing with (C.16), we obtain the key result 1 csch (-) CBZ1. 2D nj (C.21) 174 C.3 Calculation of Matrix Eigenvalues of 13g In this appendix we calculate the eigenvalues £j and eigenvectors Vj of the matrix problem Bgv = £v , (C.22) where the tridiagonal matrix Bg has the form given in (4.26b) with the coefficients d, e and / satisfying (4.85). From (4.85) it follows that d = e — / . Therefore, we get the following recursion relation for the coefficients vi of the eigenvector v: fvl-.i + (e-£)vl'+fvl+1=Q, i = l , . . . , n , (C.23a) VQ = -vi, vn = -vn+i. (C.23b) The solution to (C.23a) is vt = aCl+ + Kl_, C± = jj(C~e± [(£ - e)2 - 4/ 2] 1 / 2) . (C.24) The end conditions (C.23b) yield a + b= -aC+ - 6 C - , (C.25a) a(l + bC = -aC+ + 1 - • (C.25b) From (C.25) we get £ + = = -1 or = C-exp (2mj/n), for = 1,... ,n - 1. Substituting into (C.24) we get that the eigenvalues are & = e + 2f cos (ivj/n) , j = l , . . . , n , (C.26) which are ordered as 0 < £i < . . . < £ra since / < 0. The corresponding unnormalized eigenvec-tors are found to be V j , n = (1,-1,1, . . . , ( - l ) n + 1 ) ; v«j = s i n ^ ( / - l / 2 ) ) . , j = l , . . . , n - l . (C.27) Here vij is the J**1 component of the eigenvector Vj. These eigenvectors can be normalized and the result is summarized in Proposition 4.9. 175 Bibliography [1] N . Alikakos, X . Chen, G. Fusco, Motion of a Drop by Surface Tension Along the Boundary, (1998), preprint. [2] E. Anderson et al. Lapack User's Guide: Third Edition, SIAM Publications (1999). [3] U. Ascher, R. Christiansen, R. Russell, Collocation Software for Boundary value ODE's, Math. Comp., 33, (1979), pp. 659-679. [4] A. Bose, G. Kriegsmann, Stability of Localized Structures in Non-Local Reaction-Diffusion Equations, Methods Appl. Anal., 5, (1998), no. 4, pp. 351-366 . [5] K. Brenan, S. Campbell, L. Petzold, Numerical Solution of Initial-Value Problems in Dif-ferential Algebraic Equations, (1989), North Holland. [6] S.J. Chapman, G. Richardson, Vortex Pinning by Inhomogeneities in Type-2 Superconduc-tors, Physica D, 108, (1997), pp. 397-407. [7] X . Chen, M . Kowalczyk, Slow Dynamics of Interior Spikes in the Shadow Gierer-Meinhardt System, Center for Nonlinear Analysis report No. 99-CNA-002, (1999), preprint. [8] M . Del Pino, P. Felmer, M . Kowalczyk, Boundary Spikes in the Gierer-Meinhardt System, Center for Nonlinear Analysis report No. 99-CNA-003, (1999), preprint. [9] A. Doelman, R. Gardner, T. Kasso, Stability analysis of singular patterns in the ID Gray-Scott model: a matched asymptotics approach, Phys. D, 122,(1998), no. 1-4, pp. 1-36. [10] A. Doelman, private communication. [11] P. Freitas, A Nonlocal Sturm-Liouville Eigenvalue Problem, Proc. Roy. Soc. Edinburgh, 12A, (1994), pp. 169-188, [12] A. Gierer, H. Meinhardt, A Theory of Biological Pattern Formation, Kybernetik, 12, (1972), pp. 30-39. [13] C. Gui, Multi-peaked Solutions to a Semilinear Neumann Problem, Duke Math. J., 84, No.3, (1996), pp. 739-769. [14] C. Gui, J. Wei, Multiple Interior Peak Solutions for some Singularly Perturbed Neumann Problems, J. Diff. Eq., 158, No. 1, (1999), pp. 1-27. [15] C. Gui, J. Wei, Multiple Mixed Boundary and Interior Peak Solutions for some Singularly Perturbed Neumann Problem, Submitted JDE [16] L. Harrison, D. Holloway, Order and Localization in Reaction-Diffusion Pattern, Physica A, 222, (1995) pp. 210-33. [17] C. Hsiung, A First Course in Differential Geometry, Wiley Interscience Series in Pure and Applied Mathematics, Wiley, New York, (1981). 176 D. Iron, M . J. Ward, A Metastable Spike Solution for a Non-Local Reaction-Diffusion Model, SIAM J. Appl. Math., 60, No. 3, (2000), pp. 778-802. D. Iron, Metastability of the Gierer-Meinhardt Equations, M . Sc. thesis in the Institute for Applied Mathematics, University of British Columbia, (1997). D. Iron, M . J. Ward, The Dynamics of Boundary Spikes for a Nonlocal Reaction-Diffusion Model, European J. of Appl. Math., 11, NO. 5, (2000), pp. 491-514. D. Iron, M. J. Ward, J. Wei, The Stability of Spike Solutions to the One-Dimensional Gierer-Meinhardt Model, to appear, Physica D, Jan. 2000. A. Iserles, A First Course in the Numerical Analysis of Differential Equations, Cambridge Texts in Appl. Math. 15, (1996), Cambridge University Press, pp. 197-199. J. Keener, Activators and Inhibitors in Pattern Formation, Stud. Appl. Math., 59, (1978), pp. 1-23. M . Kowalczyk, Multiple Spike Layers in the Shadow Gierer-Meinhardt System: Existence of Equilibria and Approximate Invariant Manifold, Duke M . Journal, 98, No. 1, (1999), pp. 59-111. M . Kowalczyk, Exponentially Slow Dynamics and Interfaces Intersecting the Boundary, J. Diff. Eq. 138, No. 1, (1997), pp. 55-85. G. L. Lamb, Elements of Soliton Theory, Wiley Interscience, (1980). C. -S. Lin, W.-M. Ni, On the Diffusion Coefficient of a Semilinear Neumann Problem in Calculus of Variations and Partial Differential Equations (Trento, 1986) pp. 160-174, Lecture Notes in Math., 1340, Springer, Berlin-New York, 1988. T. C. Lacalli, L. G. Harrison, Turing's Condition and the Analysis of Morphogenic Models, J. Theo. Biol., 76, (1979), pp. 419-436. F. H. Lin, Q. Du, Ginzburg-Landau Vortices: Dynamics, Pinning and Hysteresis, SIAM J. Math. Anal. 28, (1997), pp. 1265-1293. D. Mclnerney, M . J. Ward, The Dynamics and Pinning of a Spike for a Reaction-Diffusion Model, (2000), preprint W. Ni, Diffusion, Cross-Diffusion, and their Spike-Layer Steady-States, Notices of the AMS, Vol. 45, No. 1, (1998), pp. 9-18. H. Meinhardt, Models of Biological Pattern Formation, Academic Press, London, (1982) H. Meinhardt, The algorithmic beauty of sea shells, Springer Verlag, Berlin, (1995). W. Ni, Diffusion, Cross-Diffusion, and their Spike-Layer Steady-States, Notices of the AMS, Vol 45, No 1, (1998), pp. 9-18. [35] W. Ni, I. Takagi, On the Shape of Least-Energy Solutions to a Semilinear Neumann Prob-lem, Comm. Pure Appl. Math, Vol. XLIV, (1991), pp. 819-851. 177 [36] W. Ni, I. Takagi, Point-Condensation Generated by a Reaction-Diffusion system in Axially Symmetric Domains, Japan J. Indust. Appl. Math., Vol. 12, No. 2, (1995), pp. 327-365. [37] W. Ni, I. Takagi, E. Yanagida, preprint, submitted to Tohoku Math J., (1999). [38] Y . Nishiura, Coexistence of Infinitely Many Stable Solutions to Reaction-Diffusion Equa-tions in the Singular Limit, in Dynamics Reported: Expositions in Dynamical Systems Volume 3 (editors: C. K. R. T. Jones, U. Kirchgraber), Springer-Verlag, New York, (1995). [39] NAG Fortran library Mark 17, routine D03PCF, Numerical Algorithms Group Ltd. Oxford, United Kingdom (1995). [40] I. Prigogine, R. Lefever, Symmetry-Breaking Instabilities in Dissipative Systems, J. Chem. Phys., 48 (1968), pp. 1695-1700. [41] D. Sarocka, A. Bernoff, An Intrinsic Equation of Interfacial Motion for the Solidification of a Pure Hypercooled Melt, Physica D, 85, (1995), pp. 348-374. [42] L. Shampine, M . Gordon, Computer Solution of Ordinary Differential Equations, the Initial Value Problem, W. H. Freeman publishers, San Fransisco (1975). [43] D. Stafford, M . J. Ward, B. Wetton, The Dynamics of Drops and Attached Interfaces for the Constrained Allen-Cahn Equation, submitted, Europ. J. Appl. Math. [44] X . Sun, M . J. Ward, Metastability for a Gereralized Burgers Equation with Applications to Propagating Flame Fronts, Europ. J. Appl. Math., 10, (1999), pp. 27-53. [45] I. Takagi, Point-Condensation for a Reaction-Diffusion System, J. Diff. Eq., 61, (1986), pp. 208-249. [46] A. Turing, The Chemical Basis of Morphogenesis, Phil. Trans. Roy. Soc. B, 327, (1952), pp. 37-72. [47] M . J. Ward, Exponential Asymptotics and Convection-Diffusion-Reaction Models, book chapter in "Analyzing Multiscale Phenomena Using Singular Perturbation Methods", Pro-ceedings of Symposia in Applied Mathematics, Vol. 56, AMS Short Course (1998). pp. 151-184 [48] M . J. Ward, An Asymptotic Analysis of Localized Solutions for Some Reaction-Diffusion Models in Multi-Dimensional Domains, Stud. Appl. Math., 97, No. 2, (1996), pp. 103-126. [49] M . J. Ward, Metastable bubble solutions for the Allen-Cahn equation with mass conserva-tion, SIAM J. Appl. Math., 5, (1996), pp. 1247-1279. [50] M . J. Ward, J. Wei, Asymmetric Spike Patterns for the One-Dimensional Gierer-Meinhardt Model: Equilibria and Stability, submitted, Europ. J. Appl. Math., July 2000. [51] J. Wei, On the Interior Spike Layer Solutions to a Singularly Perturbed Neumann Problem, Tohoku Math. J., 50, (1998), pp. 159-178. [52] J. Wei, On Single Interior Spike Solutions for the Gierer-Meinhardt System: Uniqueness and Stability Estimates, Europ. J. Appl. Math., 10, No. 4, (1999), pp. 353-378. 178 [53] J. Wei, Uniqueness and Eigenvalue Estimates of Boundary Spike Solutions, (1998), preprint. [54] E. Yanagida, Stability of Stationary Solutions of the Gierer-Meinhardt System, in China-Japan Symposium on Reaction-Diffusion Equations and their Applications and Computa-tional Aspects (Shanghai 1994), World Sci. Publishing, River Edge, NJ (1997), pp. 191-198. 179
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- The stability and dynamics of spike-type solutions...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
The stability and dynamics of spike-type solutions to the Gierer-Meinhardt model Iron, David 2001
pdf
Page Metadata
Item Metadata
Title | The stability and dynamics of spike-type solutions to the Gierer-Meinhardt model |
Creator |
Iron, David |
Date Issued | 2001 |
Description | A well-known system of partial differential equations, known as the Gierer-Meinhardt system, has been used to model cellular differentiation and morphogenesis. The system is of reaction-diffusion type and involves the determination of an activator and an inhibitor concentration field. It is believed that long-lived isolated spike solutions for the activator model the localized concentration profile that is responsible for cellular differentiation. In a biological context, the Gierer-Meinhardt system has been used to model such events as head determination in the hydra and heart formation in axolotl. This thesis involves a careful numerical and asymptotic analysis of the Gierer-Meinhardt system in one dimension and a limited analysis of this system in a multi-dimensional setting. We begin by studying a reduced model, referred to as the shadow system, which results from simplifying the Gierer-Meinhardt model in the limit of inhibitor diffusivity tending to infinity. This reduced model is studied in both one and in several spatial dimensions. In §2 we study the stability and dynamics of interior spike profiles for this reduced model. We find that any n-spike profile, with n > 1, is unstable on a fast time scale. Profiles with a single interior spike are also unstable but on an exponentially slow time scale. In this case the spike tends towards the closest point on the boundary. In §3 we examine the behaviour of a spike profile in which the spike is confined to the boundary. This scenario is studied in the case of a two and a three dimensional domain. It is found that the spike moves in the direction of increasing boundary curvature and increasing boundary mean curvature in two and three dimensions, respectively. Stable spike equilbria correspond to local maxima of these curvatures. We then study the case of a spike confined to a flat portion of the boundary in two dimensions. In this case it is found that the spike moves on an exponentially slow time scale. The remainder of this thesis examines the full Gierer-Meinhardt system in a one-dimensional spatial domain. In §4 we study the stability properties of n-spike equilibrium solutions to the full system. A necessary and sufficient condition is found for the linear stability of an n-spike solution. In §5 we study the dynamics of spike profiles. We derive a system of ordinary differential equations which govern the motions of the spikes in one spatial dimension. Numerical computations of this asymptotic system are compared with numerical computations of the full system. In §6 we study the effects of precursor gradients. The mathematical result of these spatial inhomogeneities in the chemical reaction is that some of the coefficients in the equations are no longer constant in space. We study the effects of spatially varying activator and inhibitor decay rates as well as a spatially inhomogeneous activator diffusivity. It is found that these spatial inhomogeneities can affect both the dynamics and equilibrium position of the spikes. |
Extent | 7438123 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
FileFormat | application/pdf |
Language | eng |
Date Available | 2009-09-23 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079985 |
URI | http://hdl.handle.net/2429/13103 |
Degree |
Doctor of Philosophy - PhD |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
GraduationDate | 2001-05 |
Campus |
UBCV |
Scholarly Level | Graduate |
AggregatedSourceRepository | DSpace |
Download
- Media
- 831-ubc_2001-611175.pdf [ 7.09MB ]
- Metadata
- JSON: 831-1.0079985.json
- JSON-LD: 831-1.0079985-ld.json
- RDF/XML (Pretty): 831-1.0079985-rdf.xml
- RDF/JSON: 831-1.0079985-rdf.json
- Turtle: 831-1.0079985-turtle.txt
- N-Triples: 831-1.0079985-rdf-ntriples.txt
- Original Record: 831-1.0079985-source.json
- Full Text
- 831-1.0079985-fulltext.txt
- Citation
- 831-1.0079985.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:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0079985/manifest