Popular in Course
Popular in Mathematics (M)
This 66 page Class Notes was uploaded by Miss Hillary Grady on Thursday October 29, 2015. The Class Notes belongs to MATH 490 at College of William and Mary taught by Staff in Fall. Since its upload, it has received 7 views. For similar materials see /class/231145/math-490-college-of-william-and-mary in Mathematics (M) at College of William and Mary.
Reviews for Seminar
Report this Material
What is Karma?
Karma is the currency of StudySoup.
Date Created: 10/29/15
Chapter 5 Reaction diffusion systems and pattern formation 51 Reaction diffusion systems from biology ln ecological problems different species interact with each other and in chemical reactions different substances react and produce new substances System of differential equations are used to model these events Similar to scalar reaction diffusion equations reaction diffusion systems can be derived to model the spatial temporal phenomena Suppose that Ut z and Vt z are population density functions of two species or concentration of two chemicals then the reaction diffusion system can be written as 8U 82U E DUW fFU7V7 51 8V 82V 39 E w t W where DU and DV are the diffusion constants of U and V respectively and FU V and CU V are the growth and interacting functions For simplicity we only consider one dimensional spatial domains The domain can either be the whole space foo 00 when considering the spread or invasion of population in a new territory or a bounded interval say 0 L In the latter case appropriate boundary conditions need to be added to the system Typically we consider the Dirichlet boundary condition hostile exterior Ut 0 Ut L 0 Vt0 Vt L 0 52 or Neumann boundary condition no ux Umt 0 Umt L 0 Vmt 0 Vmt L 0 53 The following are a list of models which has been studied by biologists and chemists 83 84 CHAPTER 5 REACTION DIFFUSION SYSTEMS Two species spatial ecological models Predator prey Ut DUUM AU 7 BU V V DVVm 7 CV kBU V 54 Here U is a prey and V is a predator AU is the growth rate of U CV is the death rate of V k gt O and BU V is the interaction between U and V Typical forms of these functions are AU aU or aU 17 2 CV CV dU JV 55 B U V dUV 7 0r 1 5U Other ecological models are Competition Ut DUUm AU 7 BU V Vt DVVM CV 7 kBU V 56 and Mutualism Ut DUUm 7 AU BU V V DVVm 7 CV 7 kBU V 57 Chemical reaction models See more introduction in Murray s book Thomas mechanism Ut DUUm a 7 U 7 pRU V Vt DVVM 01 7 V 7 pRU V 58 where RUV Thomas model is based on a speci c reaction involving the e fU gUZ substrates oxygen and uric acid which reacts in the presence of the enzyme uricase U and V are the concentration of the uric acid and the oxygen respectively and they are supplied at constant rates a and 0 Both U and V degrde linearly proportional to their concentrations and both are used in the reaction at a rate of pRU V RU V is an empirical formula For xed V RU V is nearly linear in U for small U and it approached to zero as U 7 oo Schnakenberg reaction U DUUm a 7 bU CUZV Vt vam 7 d 7 CUZV 59 This equation is from an arti cial tri molecular chemical reaction proposed by Schnakenberg X 7 A B 7 Y 2X Y 71 3X 510 k 71 Gierer Meinhardt model 52 TURING INSTABILITY AND PATTERN FORMATION 85 U2 U DUUma7bUc7 Vt DVVmdU27eV 511 Chlorine dioxide iodine and malonic acid reaction CIMA reaction 4UV UV UtDUUmm Slta7U7Wgt7 VtDvabltU71U2gt 512 Gray Scott model U DUUm 7 UZV f17 U Vt DVVm UZV 7 f kV 513 All these models can be rewritten to a standard form 814 82u 7d l 825 3m2 91quot after a nondimensionalization procedure 52 Turing instability and pattern formation We consider the nondimensionalized system ut umfuv tgt0 m 6 01 U dvm guv tgt 0 m 6 01 7 515 71020 um 1 M7570 Mt 1 07 u0z ax U0z We suppose that 140120 is a constant equilibrium solution Le 100107110 07 and 90407110 0 516 In fact the constant solution 140120 is also an equilibrium solution of a system of ordinary differ ential equations g uv tgt 0 517 In 1952 Alan Turing published a paper The chemical basis of morphogenesis which is now regarded as the foundation of basic chemical theory or reaction diffusion theory of morphogenesis Turing suggested that under certain conditions chemicals can react and diffuse in such a way 86 CHAPTER 5 REACTION DIFFUSION SYSTEMS as to produce non constant equilibrium solutions which represent spatial patterns of chemical or morphogen concentration It is noticeable that Alan Turing was one of the greatest scientists in twentieth century He was responsible for much of the fundamental work that underlies the formal theory of computation general purpose computer and arti cial intelligence During World War II Turing was a key gure in the successful effort to break the Axis Enigma code used in all German U boats which is one of the most important contributions scientists made to the victory of the war Turing s idea is a simple but profound one He said that if in the absence of the diffusion considering ODE 517 u and v tend to a linearly stable uniform steady state then under certain conditions the uniform steady state can become unstable and spatial inhomogeneous patterns can evolve through bifurcations In another word an constant equilibrium can be asymptotically stable with respect to 517 but it is unstable with respect to 515 Therefore this constant equilibrium solution becomes unstable because of the diffusion which is called di usion driven instability Diffusion is usually considered a smoothing and stablising process that is why this was such a novel concept More information of Turing instability can be found in Murray s book Chapter 14 So now we look for the conditions for the Turing instability described above If uo U0 is stable with respect to 515 then the eigenvalues of Jacobian J f f 518 gu 91 at uo 120 must all have negative real part which is equivalent to TraceJ fu 91 lt 0 Determinant fugv 7 fvgu gt 0 519 Now suppose that uo U0 is unstable with respect to 517 Then u0v0 must also be unstable with respect to the linearized equation 1 m fuu07 U0 gt fuU07U0117 75gt 07 90 E 071 11 dwm MAW U0 gt AMUO vow tgt 07 90 E 071 520 1t70 t71 t70 17571 07 gt07 90 C90 10 90 90 or in matrix notation Ilt DkIIm JII 521 Where 514535 and D 1 2 522 From the separation of variables for single equation we expect that j and 1 are both separable and we try the following form 11t 5 Wig W gtcosk7rz Wktcosk7rm 523 52 TURING INSTABILITY AND PATTERN FORMATION 87 We shall try this form of solution and the general solution for 11 would be 00 t 5 Z Wkt cosk7m 524 k0 By substitution of this form into 521 we obtain that Wkt satis es a system of ordinary differ ential equations de T T If 520 is an unstable system then at least one of Wkt would go to in nity as t 7 00 26 one of the eigenvalues of matrix Mk J 7 52sz has positive real part The eigenvalues p1 and p2 of J 7 2sz satisfy J 7 k27r2DWk 525 22 7 Af5gu7 d1k27erM 5 26 fu9v71011972ik2W2dfugvdk47r4 07 39 TraceMk fu 1 972 7 d 1k27r27 5 27 DeterminantMk L291 7 fugu2 7 k27r2dfu 971 dk47r439 39 Since fu gv lt 0 then TraceMk lt 0 is always true So Mk cannot be a source or spiral source If Mk is unstable then it can only be a saddle Therefore we must have DeterminantMk lt 0 Remember that we only need one of Mk is unstable So let s take a look of the values of hk27r2 DeterminantMk for different k or 15271392 indeed Here we de ne the function hw 7 595 7 1295 7 mm gm dwz 528 When w 0 710 fugv 7 fugu2 gt 0 because we require fugv 7 fvgu gt 0 and when 152 is large 7 00 we also have lim hw 00 So we must hope that for some 0 lt 152 lt 00 217200 hk27r2 lt 0 Since h is a quadratic equation we must require min hw lt 0 It is easy to calculate 71 dfu l 91 2d that min hw is achieved at we This implicitly implies that w dfu g gt 0 529 otherwise hw gt 0 for all w 2 O The minimum value at we is d 2 hw0 A2 595 7 fugu 7 hwy 530 which must be negative so dfu l 912 7 4dltfu9v 7 fugu gt 0 531 Summarizing all conditions we got so far if 140120 is stable with respect to 515 but unstable with respect to 517 then 1 and J must satisfy in l 91 lt 07 fug39u 7 fvgu gt 07 2 532 dfu 91 gt 07 91 7 4610072911 7 fugu gt 88 CHAPTER 5 REACTION DIFFUSION SYSTEMS The conditions in 532 give the range of parameter at such that uo 120 is possibly unstable but only 532 can only assure that certain values of hw are negative but it cannot guarantee there is a value k27r2 such that hk27r2 lt 0 To achieve that we directly calculate the conditions for hk27r2 lt O From 528 we have W279 129 7 fuguA2 7 k27r2dfu gm dk47r4 lt 0 533 and if fu gt 0 and gv lt 0 then ADA 7 kz n39zgv d gt ham 7 k27r2 E dkA 534 where D fugv 7 fugu dkA is a family of curves with dkA de ned on kZn39Zfu 00 see example below For a pair of parameters A d if d gt dkA for some natural number k then uo 120 is unstable But ifd lt dkA for any h then it is stable Summarizing the above calculation we conclude Theorem 51 Suppose that uouo is a constant equilibrium solution of the reaction di usion system ut U Afuu t gt 0 m 6 01 1 allm Aguu t gt 0 m 6 01 535 we 7 um 1 7 mt 0 7 mt 1 7 0 M0790 495 110795 M907 and the Jacobian of fg at u0U0 is Ju0u0 lt 5 71 If J at and A satisfy U 39U fugt079ult07 fu9ult07 fugu fu9ugt07 536 din gt gt 0 din gt 7 4dltfugi 7 firm gt 0 39 and d gt dkA for some k where A DA 7 l 27r2 dk 91 D fugu 7 qum 53937 k27r2qu 7 k27r27 then u0U0 is an unstable equilibrium solution with respect to 535 but a stable equilibrium solution with respect to u Afu U 1 Aguu We use the example of Schnackenberg to illustrate the theorem Consider the following and a nondimensionalized system of reaction diffusion equations describing the chemical reaction is ut ummAa7uu2u tgt0 m 6 01 at alum Ab 7 7m 25gt 0 w 6 01 uwt0 umt 1 umt0 umt 1 0 u0x A0m u0z B0m 538 52 TURING INSTABILITY AND PATTERN FORMATION 89 where a b gt 0 There is a unique equilibrium constant solution 140120 such that b b 539 U0 1 7 U0 14 and the Jacobian at 140120 is b 7 a b a2 Ju0v0 7 5 251 540 7 71 12 b L Then from simple computation the conditions 536 become 0ltb7alt ba3 ba2gt0 lt gt lt gt 541 11 7 a gt b a3 db 7 17 b 13 gt 4db a4 For simplicity and transparency of the result we select a pair of a b 1 2 then the rst row of 541 is satis ed The formula for dk in this case is 27 7 kW dkw k27r2 7 3k27r239 542 0 200 400 x 600 800 Figure 51 Parameter space d for Turing instability The rst three dk are drawn in Fig 51 Each dk is a concave up curve with a unique minimum point From Theorem 51 if 1 falls below all of dk then 140120 is stable for that parameter choice and if d is above any of dk then 140120 is unstable When 140120 is unstable the solution of linearized equation 520 has a form 0 Ilt m 2clkemktW1k CgkeMktWZk COSk7rm 543 k 0 90 CHAPTER 5 REACTION DIFFUSION SYSTEMS where pm and ngk are the eigenvalues of AJikzsz and Wlk ng are corresponding eigenvectors Since it is unstable then one or more eigenvalues ulk or 2 are positive For example if d is a point such that d gt d2 but 1 lt dk for all other k then there is one eigenvalue 12 gt 0 or complex with positive real part and all other eigenvalues have negative real part and as t a 00 Ilt z z 0125 12tW12t cos27rm 544 So the solution grows exponentially as t a 00 and it grows in a temporal spatial pattern W1kt cos27rm which is called an unstable mode of this equilibrium solution For an unstable equilibrium there may be more than one unstable modes if the parameter pair d is below two curves dk The linear analysis of Turing instability is completed now For nonlinear systems if the initial value is near the constant equilibrium solution then in a short time the orbit would be close to that of linearized equation When the equilibrium is unstable a spatial pattern is formed following its unstable mode and the pattern grows exponentially fast for a short period But after that the nonlinear structure takes over and usually nonlinear functions force the solutions staying bounded like in the logistic model thus the solution approaches a bounded spatial pattern as t a 00 Although the mechanism of the nonlinear evolution is still not clear mathematically people believe that the spatial pattern formed in the earlier stage near equilibrium solution is kept which seems to be supported by numerical simulations Many problems are still topics of current research interests More biological discussion can be found in Murray s book page 388 414 53 Turing patterns in 2D and animal coat patterns 54 Chemotactic motion of microorganisms 55 Spatial epidemic model Chapter 5 Exercises 1 Consider the diffusive predator prey model 2 8 UDU8 UAU7BUV sgt0y 0L 59 3 DviCVDUX sgt0y 0L 545 U 8V 8V 39 ho EL 8740570 EL 07 WOW Uoy7 WOW VOW Use the substitution Chapter 4 Nonlinear reaction diffusion equations In predicting most population Malthus model can only be used in a short time Due to the limited resource the growth rate of population will slow down and the population may saturate to a maximum level Thus it is natural to use a density dependent growth rate per capita and the simplest and mostly used is the logistic growth rate Together with the diffusion of the species we consider the following diffusive logistic equation 8P P E iDAPaP lt17Ngt 41 where a gt 0 is the maximum growth rate per capita and N gt 0 is the carrying capacity Diffusive logistic equation is the simplest nonlinear reaction diffusion equation but the methods of analyzing it can be used in many other problems For the simplicity we will mainly consider the spatial dimension n 1 41 Nondimensionalization Let s start with diffusive logistic equation with Dirichlet boundary condition on 0 1 ntDnman 17 tgt0 m 0L nt 0 nt L o 42 740795 1 m 6 0L The method of separation of variables fails here to nd an explicit expression of the solution to the equation and that is not surprised due to the nonlinear nature of the equation It is possible to write a formula for the equilibrium solution of the equation in an elliptic integral but that is mainly because of the special form of the quadratic growth rate Interested readers should nd that calculation in Ske51 or Ban94 page 363 We will mainly use qualitative and numerical methods to analyzing the equation although sometimes this will still require analytic tools to achieve In both qualitative and numerical methods the dependence of solutions on the parameters plays an important role and there are always more dif culties when there are more parameters Hence 45 46 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS usually we need to reduce the number of parameters by converting the equation into nondimension ized version Here we reduce the number of parameters in the equation by a dimension analysis and a nondimensionalization substitution In the following table we list the dimension of all variables and parameters To reduce the number of parameters we introduce dimensionless variables 1 U s LZ y L N 43 Then from the chain rule we have 8a8u amp8v 8y Wu 825781 33 at L2 337 amziav 8342 a and the equation becomes 31 82v ga yzv17v sgt0y 01 44 We m1 o 39 UOvy f1y7 where A mwi mm In the following sections we will use 44 in the analysis but since we normally use variables tmu instead of 8 3412 Thus we rewrite 44 in the old variables 8 82 uu17u tgt0 356 01 at 8x2 4 6 ummumna W 140790 1 90 6 071 The way of reducing parameters is not unique The new equation 46 is sometimes called a nonlinear eigenvalue problem with as a bifurcation parameter Another nondimensionalization will result in an equation with the length of interval as the parameter see the exercises 42 NUMERICAL METHOD FOR 1 D PROBLEM 47 42 Numerical method for 1d problem Investigation of many hard problems starts from a numerical computation and a successful numer ical study will lead to the discovery of the mathematical theory behind the phenomena displayed by the numerical simulations Numerical computation of reaction diffusion equations itself can well be the subject of a whole book but here we will just introduce the most natural one which may not be the most ef cient or accurate but an easy one to learn The simplest numerical algorithm for an ordinary differential equation is Euler s method For an equation du E 10057107 14050 U07 47 the linear approximation gives the simple iteration formula tn1 tn l A757 un1 un flttn7 UnAt7 where At is the step size For reaction diffusion equation we can design a similar algorithm but now we have to take the spatial variable into the consideration Again we use a set of nite points to represent the continuous variable m Let s take the example of diffusive logistic equation for our numerical computation by considering equation 46 We divide the interval 01 to n equal subintervals mimi1 i 1 2 nn 1 with mi mlquot1 1 Am m1 0 and Am 1n which we call grid size We also denote the step size in the time direction by At The initial data now is at M u m2 is mm us Hz 49 where utjmi tj tj1 1 At mi mlquot1 1 Am jZ 1 1 i g n 1 410 Because of the Dirichlet boundary condition and iii should always be zero thus we only need to calculate um with 2 g i g ii Let W 7u1lT7 411 where T is the transpose of the matrix since we would like to have column vector instead of row vector Then the idea is very similar to Euler method as we try to get vector le from 12739 The differential equation can be approximated by difference equations from the following rela tions 1 I I I I 1 I aimw N tbs17 u 8ulttizigt N u 7 u 3t At 7 8m Am 7 412 82utj uwtjmi1 7 umtj U21 7 21 71271 4 13 8m2 Am Am2 39 Therefore the differential equation it Dum gu is now ui 7 M Ui1 7 21 M71 71 7 13 gu 414 At Am2 48 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS De ne DAt r W c 17 2 415 Then 414 can be rewritten to u21 rug 1 17 27M ru 1 1 guAt 416 Notice that for i 2 or i n 416 still hold since u 0 from the Dirichlet boundary condition If we use matrix notation we obtain assuming 71 10 14 Cr0000000 u 90 u rcr000000 u 9W u 0rcr00000 u 904 a 00rcr0000 u 904 1461 0 0 0 r c r 0 0 0 u At 9w 417 U1 OOOOTCTOO 71 00000rcr0 u 90 1 000000767 90 311 000000076 7601 971 If we denote the matrix in 417 by A then 417 can be written as VM A V AGVjAt 418 where G is the vector valued function given in 417 In general is an n 7 2 vector and A is an n 7 1 x n 7 2 matrix The algorithm can easily be implemented in Maple either based on matrix computation like 418 or the iteration computation in 416 Now we adapt the algorithm to Neumann boundary value problem with umt z 0 at the endpoints In this case the boundary values are not xed thus we should consider an n 1 vector W to include the boundary points and an n 1 x n 1 matrix B For each j u can be determined by 416 except at i 1 and i n We have to extend the solution utz in a reasonable way beyond i 1 and i 71 Since homogeneous Neumann boundary condition is re ecting then we can de ne u and u Thus 71 TU 17 27M 1 T712717 guAt 27u 17 27M 1 guAt 419 And same for ln matrix notation the matrix B is of form in the case n 10 c 2r 0 O O O O O O O O r c r O O O O O O O O O r c r O O O O O O O O O r c r O O O O O O O O O r c r O O O O O O O O O r c r O O O O 420 O O O O O r c r O O O O O O O O O r c r O O O O O O O O O r c r O O O O O O O O O r c r O O O O O O O O 0 2r 0 43 DIFFUSIVE LOGISTIC EQUATION ON AN INTERVAL 49 In numerical analysis the method we describe in this section is called nite difference method Here we use a forward difference for at in 412 and a centered difference for um in 413 A backward difference is I I I 8147573 901 N 9 971 8t N At 39 We should notice that if we choose At and Am improperly the result of the numerical approximation may be ridiculous similar to Euler s method in ordinary differential equations For example when D 1 At Am 1 then 412 becomes 421 14 14 7 2u ug l 422 Suppose that we start with an initial value ul 0 0 0 0 0 1 0 0 0 0 0 and consider the Dirichlet boundary value problem Then the iteration sequence is u2 107070707177 27170707070 us 070707 17747677417170070 u4 00176157 201576100 423 However as you remember the solution of diffusion equation with a positive initial value should keep positivel Such an iteration will not provide an approximation to the true solution thus it is called unstable It can be show that the above forward nite difference method is stable if DAt 1 M g 5 424 r The stability condition 424 means that the time step size At has to be small relative to Am2 For example when A 01 and D 1 then At must be less than 0005 To reach T 1 200 iterations are needed On the other hand more advanced numerical methods are available for the approximations of these types of partial differential equations and some do not require such stability criterion But the forward nite difference method is the simplest and most straightforward 43 Diffusive logistic equation on an interval In this section we use Maple program to do some numerical experiments to motivate some of our qualitative analysis later We rst consider the Dirichlet problem a u U17U tgt0m 01 8t 8m2 ut 0 ut 1 0 425 140790 90 90 E 071 We use the algorithm described in the last section and grid size Am 01 In the program we assume 1 and fm sin7rm The simulation shows a sequence of decreasing arches and the height of the arch decreases from u 1 when t 0 to u z 005 when t 03 see Figure 41 50 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS 1 1 09 09 08 08 07 07 06 06 05 05 04 04 03 03 02 02 01 01 0 0102 03 04 05 06 07 08 09 1 0 0102 03 04 05 06 07 08 09 1 Figure 41 a u0m b u03z when A 1 If we increase the number of time iterations we will see that decreasing trend will continue so that the solution appears to approach u 0 as t a 00 We could connect this with our early analysis of diffusive Malthus equation 276 In the Malthus model the solution will also tend to u 0 as t a 00 since the critical number of survival or extinction is a0 71392 when L 1 and A 1 lt 71392 On the other hand the logistic growth is slower than the linear growth in Malthus model thus the population in diffusive logistic equation will also become extinct for all A lt 71392 This intuitive idea can be showed in a more mathematically rigorous way in later sections and this argument is also very valuable in our later analysis So our natural question is when A gt 712 what will happen Will the population have an exponential growth or still die out Let s turn to our program again but change the parameter to A 20 which is larger than 712 and try two different initial data u0z sin7rz and u0 m 018in7rm 1 1 09 09 08 08 07 07 06 06 05 05 04 04 03 03 02 02 01 01 0 0102 03 04 05 06 07 08 09 1 0 0102 03 04 05 06 07 08 09 1 Figure 42 a u03z when u0z sin7rm b u03z when u0z 013271011 Both are when A 20 From the simulation we can see that in the rst case the solution decreases to the pro le in 42 a and in the second case the solution increases to the pro le in 42 So zero or in nity We are more confused now Let s push the computational ability of our Pentium IV computer to limititry the step number m1000 Now let s see what we will have after a long 05 time unit 43 DIFFUSIVE LOGISTIC EQUATION ON AN INTERVAL 51 1 1 09 09 08 08 07 07 06 06 05 05 04 04 03 03 02 02 01 01 0 0102 03 04 05 06 07 08 09 1 0 0102 03 04 05 06 07 08 09 1 Figure 43 a u05z when u0z sin7rm b u05z when u0z 013271011 Both are when A 20 Two almost identical functions To convince yourself you may try change the initial function to another function and you will that the solution will again tend to the function in Figure 43 a or b Thus our numerical experiment for A 20 can reach to a conclusion the solution eventually will approach a limit function with maximum value at about u 055 This limit function does not change when if changes thus it has to be an equilibrium solution of the equation and it satis es 8214 Au1iu 0 m 6 01 426 u0 u1 0 1 1 09 09 08 08 07 07 06 06 05 05 04 04 03 03 02 02 01 01 0 010203040506070809 1 0 0102 03 04 05 06 07 08 09 1 Figure 44 a u05 z when u0 z sin7rz and A 100 b u05 z when u0 z sin7rz and A 400 If we change the value of A again but to values greater than 712 then we will nd similar phenomena only the equilibrium solutions have different pro les In fact we can observe that the equilibrium solution for larger A is larger and when A is really large the equilibrium solution almost equals to 1 for all z except two narrow intervals near z 0 and z 1 Such a solution is usually called a boundary layer solution We can now summarize our numerical exploitation 1 When A lt 71392 all solutions of 425 tends to the equilibrium solution u O 52 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS 2 When A gt 71392 all solutions of 425 tends to an equilibrium solution u We recall that in the original model A aLZD where a is the maximum growth rate per capita D is the diffusion constant and L is the size of the habitat Then we can see that the concept of critical patch size is still relevant here and L0 Da7r is the critical patch size when a and D are xed But in the case of L gt L0 the population will not grow exponentially as in the diffusive Malthus model Instead the population will reach an upper limit just similar to the case of ordinary differential equation P P1 7 P Moreover when L is very large corresponding to A very large the equilibrium solution almost equals to 1 in most of interior of the habitat and not surprisingly P 1 is also the carrying capacity of P P1 7 P Thus even though some population is lost due to the outer ux at z 0 and z 1 when the habitat is large enough to reproduce new lives the population will still saturate to a level close to the carrying capacity As a comparison we now turn to the Neumann boundary value problem 2 8 xu1iu tgt0 356 01 8t amz umt0 wt 1 o 427 140790 1 z 6 071 We leave the numerical experiments to the reader We will nd that no matter what value A is all solutions will tend to the equilibrium solution u 1 that is exactly the carrying capacity A curious question is that in the cases of Dirichlet problem with A gt 71392 and Neumann problem there are at least two equilibrium solutions for the equation u 0 and u u for Dirichlet and u 0 1 for Neumann Why do all solutions tend to only one of these two solutions not the other one Or maybe there are some solutions which will approach to the other equilibrium solution and we are just cheated by the computer We shall answer this question in the next section 44 Stability of equilibrium solutions In diffusive Malthus equation there is a critical patch size L0 Da7r such that all solutions are decaying to zero when L lt L0 and all solutions are exponentially growing when L gt L0 This phenomenon can also be explained in term of the stability of the equilibrium solution u 0 We can say that u 0 is stable when L lt L0 and it is unstable when L gt L0 Now we de ne a stable equilibrium solution Suppose that vz is an equilibrium solution of 2 8 D8 gu tgt0 z OL at 8x2 ut0 ut L 0 428 u0x m E 0 L which means vz is a solution of 821 Dy gv 0 m E 0 L 110 UL 0 429 44 STABILITY OF EQUILIBRIUM SOLUTIONS 53 vz is said to be stable if for any 6 gt 0 there exists 6 gt 0 such that whenever 7 g 6 for all m E 0L 430 we have lutz 7 g 8 for all z E 0L 431 where ut z is the solution of 428 with u0 z If an equilibrium solution vz is stable and moreover tlim lutm 7 O uniformly for z E 0L 432 400 then we say vz is asymptotically stable An equilibrium solution vz is unstable if it is not stable Roughly speaking if all solutions of 428 with initial values near vz stay close to 12m then vz is stable if all solutions of 428 with initial values near vz will approach vz as t 7 00 then vz is asymptotically stable If vz is asymptotically stable then it is also stable The stability for other boundary conditions are de ned similarly To illustrate the notion of stability we consider several examples of the linear diffusion equations where solution can be explicitly found 1 Dirichlet problem of diffusion equation at Dam ut 0 ut L O u0 z fm and 00 the solution is Zonemp7Dn27r2tL2 sinmrzL The only equilibrium solution is u O n1 and it is asymptotically stable to Neumann problem of diffusion equation at Dam umt 0 umt L O u0 z fm CO and the solution is c Z onexp7Dn27r2tL2 cosmrzL where c L 1 fOL fzdm Any n1 constant function u c is an equilibrium solution and all ofthem are stable from the formula but none of these equilibrium solution is asymptotically stable since for any a 0 there is always other constant solutions nearby 03 Dirichlet problem of diffusive Malthus equation at Dum au ut0 ut L O 00 u0 z fm and the solution is Zonexpat 7 Dn27r2tL2sinmrzL The equilibrium n1 solution a 0 is asymptotically stable when L lt L0 and it is unstable when L gt L0 The stability of equilibrium solutions for linear diffusion equation can easily be determined so the question is how about nonlinear equation Let s recall the situation of ordinary differential equations 1 scalar equation Suppose that u uo is an equilibrium solution of u 9a then it is stable if g u0 lt 0 and it is unstable if g u0 gt O 2 planar systems Suppose that 1412 140120 is an equilibrium solution of u guv 1 hu u then it is stable if both eigenvalues of Jacobian Ju0 120 have negative real part and it is unstable if both eigenvalues have positive real part 54 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS In both cases the stability is determined by the linearization of the right hand side of the equa tion and that would also be true for the case of reaction diffusion equation The stability of an equilibrium solution can be determined by the following criteria Theorem 41 Suppose that 11m is an equilibrium solution of 428 Then 1 The equation D ltzgtg ltvltzgtgt ltzgti mu 433 MWMMQ has an sequence of eigenualues al a2 such that M is a real number Marl lt m and lim ai 700 1700 2 The eigenfunction changes its sign exactly 7 1 times in 0L39 3 If all eigenualues M are negatiue then 11m is a stable equilibrium solution if there is a positiue eigenualue then 11m is an unstable equilibrium solution You can certainly nd that the theorem resembles its easier counterparts for scalar equations and planar systems The fact that 433 has only real eigenvalues is because 433 is also a scalar equation and later we will learn that the linearization of solutions of reaction diffusion systems can also have complex eigenvalues as planar systems The proof of the theorem is beyond the scope of this notes thus it will not be given here lnterested readers can nd the proof of the rst two parts of the theorem in Strauss Str92 Chapter 11 and the proof of the last part in 7 In the following we will apply Theorem 41 and one of its con eqnence the L rst eigenvalue M of 433 does not change sign thus it can be assumed as a positive function b1 LULLU pundiug t0 the Now let s apply Theorem 41 to diffusive logistic equation with Dirichlet boundary condition u 0 is an equilibrium solution of 425 The linearized eigenvalue problem at u 1 is jgm Ag0 m M x 6 071 434 0 M1 07 where 9u u17 Since 90 17 then satiSEes In m 6 01 435 0 M1 0 F rom Section 21 a 7 A must be 7n7r2 n 2 1 So an A 7 n27r2 In particular a1 A 7 71392 Thus u 0 is asymptotically stable if A lt 712 and it is unstable if A gt 712 The stability of a constant equilibrium solution is easy to determine by Theorem 41 but for the non constant ones it is usually not so easy However we are able to show that the positive equilibrium solution which we observe in the last section must be asymptotically stable via the following calculation Note that we have not shown that that equilibrium solution actually exists which we will in the following sections so here we prove that it is asymptotically stable under the assumption of its existence 44 STABILITY OF EQUILIBRIUM SOLUTIONS 55 Proposition 42 Suppose that is a positiue equilibrium solution of 425 for gt 71392 Then is asymptotically stable Proof We prove it using a trick of calculus integral by parts Let al 11 be the rst eigenvalue eigenfunction pair Ial is the largest eigenvalue of lt0gt gtlt1gt 0 43936 where gu u17 From Theorem 417 we could assume that gt1z gt 0 for all z E 017 and we only need to show that m is negative hence all other eigenvalues are also negative We multiply 425 u gu 0 by 11m7 and multiply 436 by um7 subtract and integrate on 07 17 then we obtain it96 A9 U gt Mn 90 E 071 1 A u z gt190 1 1 437 A we 7 ultzgtg ltultzgtgti 1ltzgtdz 7M1 ultzgt 1ltzgtdz We handle the rst integral by integrating by parts as follows A1 u ltzgt 1ltzgtdz 1 1ltzgtdwltzgt 1 gt1zu xl i wowmm 438 0 e 1 WWW here the terms on end points vanishes because gt10 gt11 O and the second part can be manipulated similarly 1 1 7 lltzgtultzgtdz ultzgtd altzgt 0 0 1 7uz izl u ltzgt gt ltzgtdz 439 0 1 ux gtmdm 0 Thus 1 0 dimmer ltzgtultzgtidzo 440 which makes the identity 437 becoming A llgltultzgtgt e ultzgtg ltultzgtgti 1ltzgtdz 7M1 1ultzgt 1ltzgtdm 441 56 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS Since gu u17 u then gu 7 ug u u2 and the integral on the left hand side of 440 now is f01ux2 gt1xdx which is positive since we choose gt1x to be positive On the other hand if gt1x is positive then fol ux gt1xdx is also positive since is a positive solution Therefore al must be negative from the expression in 441 45 Motivation from algebraic equation In many aspects the reaction diffusion equation ut umgu is similar to the ordinary differential equation u gu and the behavior of the equilibrium solutions of the two equations is one of them In this section we consider an equation of equilibrium solutions of u gu that is an algebraic equation gu 0 But the methods we use here can be later used in more complicated reaction diffusion equation as well Consider 35355 x 17 5 0 442 where 6 is a parameter Clearly x 0 is always a solution to the equation no matter what 6 is Other solutions of the equation satisfy x5 x 1 7 6 O which is not solvable algebraically Here we apply two advanced mathematical methods perturbation and bifurcation to this problem First we discuss the solutions of xx5 x 17 8 0 when 6 is near 60 1 When 6 1 x 0 is the only solution of x xx5 x 1 7 6 O and it is a double zero of x in the sense that f 0 0 If you draw the graph of x for 6 701 0 and 01 you will nd that x 0 has another zero near x 0 when 6 344 0 but close to 0 Thus a bifurcation occurs at 6 0 when 6 lt 1 x 0 has two zeros x 0 and a negative one when 6 1 x 0 is the only one and when 6 gt 1 x 0 also has two zeros x 0 and a positive one This is called a transcritical bifurcation and it can be depicted as in Figure 45 a The question is what is the value of the other zero of x when 6 is near 1 We know it is close to x 0 but in applications we may want to know the value more precisely Figure 45 a Transcritical bifurcation b pitchfork bifurcation 45 MOTIVATION FROM ALGEBRAIC EQUATION 57 Now we use the perturbation method suppose that r 6 7 1 is small we assume that the zero point mg has a form mg 0 rml 7 sz rsmg r4m4 r5m5 1 higher order terms 443 and substitute this form into 5 z 1 7 6 0 then we obtain rml 7 sz 1 5 rm 1 7 sz rsmg r4m4 r5m5 17 1 r 0 444 thus 11 2340 571 and we get an approximate value of mg mg 1757 1755R5 446 with the remainder R is in an order smaller than 1 7 65 The equation 446 gives a very accurate estimate of mg but from a rigorous mathematical point of view the reason that R is indeed a small quantity and hence whether the expansion of mg is convergent is questionable But nevertheless we gain useful information from this simple calculation In fact this argument can be made rigorous in analytic mathematics via one of most important theorem in mathematical analysis the implicit function theorem To introduce implicit function theorem we think the function 5 z 1 7 6 is a function depending on two variables 6 and m F6m 5 z 1 7 8 Then 6 z 10 is a zero point of F The level set F6m 0 near 10 is a curve if the gradient VF 8F86 8F8z is not zero at 10 and the points on the curve are solutions of F z 0 If we use the linearization of the level curve we obtain 3F10 3F10 0F6m F10T671Tz 447 Thus a solution gm of F z near 10 has a rough form 7 8F8z10 m 4 imp 71 448 You can verify that 448 and 446 are consistent A simple version of implicit function theorem is as follows and 80 0 W0 7 80 1 F4 449 zm07 36 6780 where lim 0 Moreover 6m6 is a smooth curve in 8 m space 5750 58 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS The implicit function theorem can also be used in calculating the solution of 442 when 6 is near 3 When 6 3 z 1 is a solution Using the perturbation method or implicit function theorem we can calculate the solution as of the equation for g z 3 25 7E5732H39 450 mg1 673 The calculation of 450 is left as exercise The perturbation and bifurcation method can be used to solve the equation near a known solution and both methods only work in a local neighborhood of the known solution Neither gives information of global bifurcation diagram of F z 0 In fact there is an easy way to do that instead of express the solution x in term of parameter 8 we can solve 6 in term of z easily the solutions of 4124 are either z O or 6 5 z 1 451 451 is especially useful when you draw the bifurcation diagram of the problem although it does not provide the information of z in term of 8 So from our observation of an algebraic equation we learn the following lessons 1 Bifurcation points are the points where fm z O 2 Perturbation method can be used to nd approximate solutions near a known solution 3 When you cannot solve the solution x in term of parameter 8 you may solve 6 in term of m All three ideas will be useful when we return to the discussion of equilibrium solutions of reaction diffusion equation 46 Perturbation and local bifurcation In this section we use the ideas in the last section to study the equilibrium solutions of diffusive logistic equation so we consider the equation Wm mu 7 u 07 m 6 071 Ax52gt u0 u1 0 u 0 is always an equilibrium solution for any gt 0 so we would like to know equilibrium solutions other than u 0 First on the bifurcation points of the equation we have Proposition 44 If 0 0 0 is a bifurcation point of the equation 452 then a 0 must be an eigenvalue of M Ag0 m M m E 071 453 gtO M1 07 where gu u17 46 PERTURBATION AND LOCAL BIFURCATION 59 This is similar to F6m 0 where the bifurcation points must satisfy Fm O The rigorous proof of Proposition 44 requires some advanced knowledge which we will not try here The basic idea is the following when 0 is not an eigenvalue of 453 when A A0 then the solutions of 452 near A0 0 must on a curve A uA which must be the zero solutions A 0 By using Proposition 44 we can now determine the possible bifurcation points since g 0 1 u 7 A 7712712 thus to make u O A must be 712712 Recall that when A lt 712 u 0 is a stable equilibrium solution and when A gt 712 u 0 is unstable And from our numerical experiments we also nd that when A gt 712 452 has a stable positive equilibrium solution Now we can use perturbation method to determine the analytic property of this stable solution Since A u 7120 is a bifurcation point then we look for a solution Au which is close to 7r2 0 We introduce a small parameter 8 and we look for a solution with form A7r2a16a262a363 u06u162u263u3 454 where a is a real number and u is a function We substitute the form 454 into equation 452 then the equation becomes 0 u Au17u Eu52rLZ4 71392 l 016 l 0262 l 6u1 62142 17 6741 7 62142 455 LL1 7r2u1 2 LL2 7r2u2 alul 7 wzu and the boundary condition becomes 710 57110 52520 m 0 57111 52521 711 456 Thus ul and U2 satisfy 711 52141 0 m 6 01 7110 7111 0 457 1247 7r2u2 alul 7 71214 0 m 6 01 u20 u21 0 458 The solution of 457 must be u1z ksin7rz but the constant k can not be determined from 457 To determine k we have to look further for U2 The solvability of U2 is not a necessity Here let s play the multiply multiply subtract integrate trick in the proof of Proposition 42 we multiply 458 by ul multiply 457 by U2 subtract the two equations and integrate over 01 then we get 1 515357 wzu mdx o 459 0 Notice that the equality 459 contains only function 741 so that is an equation which you can use to solve the constant kl By substituting ul ksin7rm to 459 we obtain 1 1 52153 sin37rxdz 51k sin27rmdx 460 0 0 thus k 3a187r Notice that we can select a1 arbitrarily because of scaling so we can assume that al 1 Then U2 can also be solved but with another undetermined constant 3 i 3 3 9 u2z 7 cos7rz k1 s1n7rz 128712 1672 7 and k will be determined in the equation of ug The formula of 742 above is calculated by Maple Put together all results in this section we have the following local bifurcation picture of 452 cos27rz z cos7rz 60 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS Theorem 45 Near A7 u 7139270 the solution set of 452 consists of two parts a line of constant solutions A7 0 and a curue of non constant solutions A67 u6 which crosses the line of constant solutions at 7120 ue is positiue when A gt 7139 and it is negatiue when A lt 71392 Moreouer g 712 5 l 0282 higher order term k1 sin7rz mm 9 z cos7rz 7 W 3 cos27rz 16W2 icsinn39z 62 ficoswrz W Sr 32 47 Global Bifurcation We have found that it is hard to express the equilibrium solutions in term of parameter A In this section7 we try the opposite express the parameter in term of the solution u To do that7 we use a dynamical system approach The differential equation u Agu 0 can be viewed as a system of rst order ordinary differential equations u u 1 7Agu 462 and there is a rich geometric theory about the phase plane of the system 462 A solution of the boundary value problem 452 is the right half of an orbit of 462 such that uO u1 0 Now we show that the parameter A can be solved in term of u057 the maximum point of To derive a formula of A in term of um we consider a more general equation D 0 0 L u ltzgtgltugt z elt gt 463 uO uL O7 and we can always convert 463 to nondimensionalized version u x Agu 07 m E 017 4 64 mmmno 39 where A LZD The equation in 463 can be rewritten as u u 1 7D Igu 465 A positive solution of 463 is equivalent to a solution orbit of 465 such that l uO uL O 2 uL2 maxuz7 uLQ 0 3 gt 07 for z E 07L7 uz gt 0 for z E 0L27 and uz lt 0 for z E L27L 4 7 GLOBAL BIFURCATION 61 AmNNNNN Figure 46 Phase portrait of diffusive logistic equation From the phase portrait of 465 where gu u1 7 747 the possible values of uL2 are from u 0 to u 17 so we consider a solution of 465 with 0 lt uL2 lt 1 We multiply 463 by 1z7 and integrate it over z L2 for some z E 0L27 then we obtain L2 o Du zu x 9uzu zldz D 2 gt6 dz ALZgumdum uL2 4 ltL2gti2 7 64642 gltugtdu w 7 66 GltultL2gtgt 7 Gltultmgtgtv 466 where Cu fgtdt Since uL2 is a constant7 we denote it by s The equality 466 can be written as W Gltsgt e Guml 467 So after an integration7 we obtain a new equation 4677 which is a rst order ordinary differential equation gdae 7 C14 468 Here we take dudm to be positive since z E 0L2 We integrate 468 from u 0 to u s then we obtain LZ 5 Ma aB 46 469 2 0 2 0 CSFGW S L 2D du 470 0 LL 62 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS Using the relation A LZD we obtain 5 1 2 A 2 du 471 lt0 Gs7Gu gt for 464 Therefore if u is a solution of 464 with uL2 s then the parameter A can be expressed in term of s as in 471 In another word the global bifurcation diagram can be drawn as a graph As s where As is de ned in 471 and s 6 01 The integral in 470 is often called time mapping as it maps a solution to L which can be thought as the half of the time traveled from u 0 to the next time u O The integral in 471 is not easy to evaluate For example when gu u17u the expression of As is S 2 xE A 2 d 472 S lt0 333 7 232 7 3772 1 2u3 u This integral is improper since when u s the integrand has zero denominator but this improper integral is convergent since the integral near u s is approximately for 6 gt 0 small 9 du du 975 V333 7 232 7 3772 1 2773 7 3 du J 4 73 i 975 s7u3s3u723272su72u2 39 N 1 9 du 7 2 6s76s2 55 xSiu 6s76s2 thus the whole integral is convergent On the other hand from the phase plane the integral is proportional to the time that the solution travels from Ou 0 to s 0 which is a nite time Maple is able to calculate this integral for any 3 E 0 1 thus we can obtain numerically a graph of the function As xS 7 Even with the numerical help of Maple a little more calculation of the function As would help us to better understand this bifurcation diagram and generate a better numerical bifurcation diagram Here we use some integral calculus to show two properties of the bifurcation diagram The key here is the calculation of the timemapping S 1 Ts du 474 0 Gs 7 Cu since As 2Ts2 and As 4TsT s Thus critical points of Ts and As are equivalent and if Ts is monotone so is As To calculate Ts we use a change of variables u 3w to convert the integral in 474 to s 1 J i 1 S J m 4 might 4 Cs7Gsw 43975 The advantage of the new integral is that the bounds of the integral are now independent of 3 variable that would be helpful when we consider the derivative of Ts We differentiate Ts and 4 7 GLOBAL BIFURCATION 63 simplifying the expression we obtain 1 H s 7H 311 T Ad 476 S 2168 7 Gems2 W7 l where Hg12 C12 7 1291 Now if we let 91 1 7 122 then Hg12 16122412 7 3 and H12 2122 7 1 gt 0 if 1 6 01 In particular Hg1 is a strictly increasing function when 1 6 01 and Hgs 7 Hgsw gt 0 for all 11 6 01 and s 6 01 Therefore the integral in 476 is positive since the integrand is positive and Ts gt 0 for all s E 0 1 Figure 47 a Graph of 91 1 7 122 b Graph of Hg12 U24 U 7 3 The domain ofthe function Ts and s is s E 0 1 from our observation of the phase portrait in Figure 47 So if we can determine the behavior of the time mapping when 3 7 0Jr and s 7 1 we will have a complete qualitative picture for Ts and s We have found previously that 71392 is the bifurcation point but this conclusion can also be obtained through the calculation of time mapping In fact when 3 7 0 Cs z 10 G 0s 05G 032 0552 thus 1 2 7139 l T d 477 351 8 0 17 102 w and thus limJr s 71392 So once again we nd that 71392 is the bifurcation point The integral s70 Ts approaches to 00 as s 7 1 from the calculation in 473 Summarizing the above we nd that limJr s 7r2 lim 3 00 and s gt 0 s 6 01 478 570 971 Now we use Maple to visualize the bifurcation diagram s 8 Here is the Maple code gt withplots gt b1 fx b x de ne the nonlinearity gt F1intf x FunapplyF1 x nd the antiderivative of gt fsolvef x solve the zeros of fu which determines the domain of the timemapping gt Initial0 End1 domain of the timemapping gt numplots100 the number of plotting points 64 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS V lambdaarray1 numplots uarray1 numplots gt lambda1 evalf Piquot2 u1 O the bifurcation point gt for n from 2 to numplotsl do gt sInitialEndInitialn1numplots gt evalfInt1sqrtFsFx x0 s gt lambdan 2 oquot2 un s gt end do gtplot seq lambda n un n1 numplots1 0 100 tickmarks 10 11 Now we obtain the bifurcation diagram in A7 5 coordinate system Figure 48 Bifurcation diagram of u Au17 u O7 u0 u1 O From Figure 48 and previous discussions7 we shall conclude our investigation of the diffusive logistic equation with Dirichlet boundary condition Theorem 46 Consider 8 82 u UAu17u tgt0m 01 at 8x2 4 7g 74250 ut1 o lt gt u0x gt 07 m 6 01 1 When 0 lt A lt 71392 472 has only one non negatiue equilibrium solution u 0 which is asymptotically stable For any gt 0 tlim utm 0 400 2 When A gt 71392 472 has exactly two non negatiue equilibrium solution u 0 which is 48 ALLEE EFFECT 65 unstable and a ux which is positive and asymptotically stable For any gt 0 tlim ut m 4 00 Biologically diffusive logistic equation provides a more reasonable alternative to diffusive Malthus equation where population has an exponential growth when the size of the habitat is larger than the critical patch size In logistic model the population will approach a positive equilibrium solu tion if the size of the habitat is larger than the critical patch size which play a similar role as the carrying capacity as the ODE case Mathematically we have sampled a few powerful tools in studying nonlinear evolution equation perturbation bifurcation time mapping and numerical methods An important omission here is the maximum piinciple due to the time constrain and it is the key to complete the proof of Theorem 46 in proving the convergence of all solutions to the equilibrium solutions 48 Allee effect In this section we introduce some variants of the diffusive logistic equation In logistic equation dP P E 7 kP lt17 7 P gP 480 the growth rate per capita gP k1 7 PN is a strictly decreasing function with respect to the population density P which considers the crowding effect For some species a small or sparse population may not be favorable since for example mating may be dif cult This is called Allee e ect Mathematically gP will not have maximum value at P 0 If the growth rate per capita gP is negative when P is small we call such a growth pattern has a strong Allee e ect A typical example is dP P E412 lt17Ngt lt 71gtP9P 481 where 0 lt M lt N M is the sparsity constant and N is the carrying capacity The graph of an example of growth rate fP and growth rate per capita gP is in Figure 410 If the growth rate per capita gP is smaller than the maximum but still positive for small P we call such a growth pattern has a weak Allee e ect A typical example is 412 17 lt 71gtP9P 482 where M lt 0 lt N The graph ofa example of growth rate fP and growth rate per capita gP is in Figure 49 The dynamics of strong and weak Allee effects are different the weak Allee effect has the qualitatively same dynamics as the logistic equation and all positive solutions tend to P N as t a 00 but in the strong Allee effect case there is a positive unstable equilibrium point P M which is a threshold value and there are two non negative stable equilibrium points P 0 and P N In the latter case if 0 lt P0 lt M then limtH00 Pt O and if M lt P0 lt N then limtH00 Pt N 66 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS In this section we consider the corresponding spatial model Suppose that we make the same nondimensionalization as in Section 31 then the nondimensionalized equation is 814 8214 a wkuliuuia tgt0 m6 01 ut 0 ut1 0 740790 U07 90 E 071 483 where foo lt a lt 1 In particular we have strong Allee effect if 1 gt a gt O and we have weak Allee effect if a lt O Figure 410 Srong Allee effect growth rate left growth rate per capita right The mathematical tools developed in the last few sections can be used here for a complete study of the equation which will be done over a series of exercises But we want to point out two new phenomena which does not happen in the case of diffusive logistic equation Existence of equilibrium solutions First we consider the case when a is close to 1 This is the case of strong Allee effect with the sparsity constant close to the carrying capacity We can use the phase portrait of the rst order dynamical system to study the existence of equilibrium solutions The equilibrium solutions are equivalent to an orbit starting from and ending on v O The system in this case is 1 U 1 7u17 7 a 484 4 8 ALLEE EFFECT 67 But from the phase portrait Figure 411 a there is no orbit starting from v 0 will hit 1 0 again So there is no positive equilibrium solutions Figure 411 Phase portraits u u17 7 02 0 left u u17 7 08 0 right So the question is when the equation has an positive equilibrium solution We recall the time mapping formula 5 1 2 2 du 485 0 Gs 7 Cu where s u12 To make this formula meaningful we must have Cs7Gu 2 0 for all u E 0 s In this case it means 5 13 7 Cu gtdt Z 0 where gt u17 7 a 486 U Let s check the graph of By integration we obtain i14a13a27212a1a Cu7 414 3 u 214 7U 414 3 u 2 487 From the graph of Cu when a 08 Figure 412 a one can see that there is no any 8 satisfying 486 thus there is no any valid 8 for u12 On the other hand when a is closer to 0 than 1 for example a 02 see Figure 412 b then there are 8 satisfying 486 and from the phase portrait Figure 411 a there are orbits starting from and also ending on v O The valid values for s u12 is from where Cs 0 to s 1 For 8 gt 1 95 lt 0 there is no such solution with u12 3 Therefore summarizing the observation above we can conclude Proposition 47 Suppose that is a positive solution of iimm gu o m 6 01 710 711 0 43988 and u12 3 Then 3 satis es 93 gt o Sgtdt 2 o for u e 0 s 489 68 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS When 0 lt a lt 1 this is only possible when the positive hump of the function has larger area than the negative hump 239e 0 lt a lt 05 Hence when 05 lt a lt 1 there is no positive equilibrium solution no matter what A is When 0 lt a lt 05 there exists positive equilibrium solutions for some A which we discuss below Figure 412 a graph of Cu when a 08 b graph of Cu when a 02 095 09 085 08 075 07 065 06 055 05 045 04 035 80 160 120 140 160 180 260 220 Figure 413 bifurcation diagram of u Au17 7 02 0 140 141 0 Bistable equilibrium solutions When 0 lt a lt 05 the global bifurcation diagram Figure 413 shows that there exists Al gt 0 such that there is no positive equilibrium solutions when A lt Al there is one when A Ar and there are two positive equilibrium solutions when A gt At Remember that 0 is always an equilibrium solution Then we have three non negative equilibrium solutions when A gt Al This bifurcation diagram reveals several biologically interesting phenomena First the existence of Ar implies the existence of a critical patch size Ll xAj but this Ll is not obtained from linearization like diffusive logistic equation From the bifurcation diagram we estimate Ar 70 which is much larger than the one in diffusive logistic equation Ar 71392 z 10 That shows the negative impact of allee effect to the population The critical patch size is larger so the species needs a larger living environment to survive 48 ALLEE EFFECT 69 Secondly even when gt Al the species still may not survive In this case there are three equilibrium solutions 0 u1z with larger u12 and uz with smaller You can check that u 0 is always an asymptotically stable equilibrium solution so when the initial population is too small in some sense then the population will still drop to zero But if the initial population is large then it will tend to an positive equilibrium solution which is in fact the larger one among the two positive equilibrium solutions Therefore there are two stable equilibrium solutions and either one can be the asymptotic limit of the population evolution This is similar to the ODE case but more complicated To understand the threshold phenomenon here we need to introduce the concept of the basin of attraction Figure 414 phase portrait of m z17 z1 m y 7y Suppose that is an asymptotically stable equilibrium solution Then the basin of attraction of is the set ofinitial functions u0z such that ut z Mm where ut z is the solution of reaction diffusion equation with u0z u0m From numerical experiment we can easily nd initial values belonging to the basins of attraction of 0 or From a deep result in the dynamical system if u0z lt u2m then u0z belongs to the basin of attraction of u O and if no gt u2m then it belongs to the basin of attraction of But there are initial values which may not larger or smaller than u2m and there is no complete results regarding these solutions But it is known that the two basins of attraction are separated by a surface in the space of all continuous functions If the initial value is one of the functions on this surface then the limit of the solution is neither 0 nor u1m but u2ml This surface is called stable manifold of the equilibrium solution Although in the experiment it is vary hard to catch a solution on the stable manifold but the stable manifold serves as the role of the threshold here This situation is similar to the following ODE system see Figure 414 m M17 x1 z y 7y 490 On this phase plane the line z 0 is the stable manifold of the saddle point 00 and it separates the basins of attraction of stable equilibrium points 71 0 and 10 Biologically the population for strong Allee effect has a conditional persistence for gt A1 in 70 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS which the survival of the population depends on initial population distributions In the logistic case the population has a unconditional persistence for A gt 712 The surface77 mentioned above is called threshold manifold This threshold manifold separates the set of all non negative initial distributions into two disconnected subsets which we can call above threshold77 A and below thresold77 B sets For any initial distribution in A the asymptotic state is the stable steady state U1 On the upper branch of the bifurcation diagram for any initial distribution in B the asymptotic state is the stable steady state 0 To conclude this section we mention that the bifurcation diagram for the weak Allee effect case u1 7 7 a where a E 71 0 is a combination of the logistic and strong Allee effect cases There are two critical values A1 7r2f 0 and Al lt A1 which divide the parameter space into three parts extinction regime A E 0A conditional persistence regime A E AA1 and unconditional persistence regime A 6 A1 00 In the following numerical bifurcation diagram A1 7r2f 0 10712 x 9870 and Al x 3353 Figure 415 Bifurcation diagram of a Au17 01 O u0 u1 O The bifurcation diagram in Figure 415 allows for the possibility of hysteresis as the diffusion coef cient A or the habitat size is varies Suppose we start with a large size habitat and then slowly decrease the size Then initially the population will stabilize at the unique steady state solution u1A However when the habitat is too small when A lt Al the population collapses quickly to zero To salvage the population we may attempt to restore the habitat by slightly increasing A so that A gt Al But if the population has dropped below the threshold u2A at that moment then the population cannot be saved since it is now still in the basin of attraction of the stable steady state a O 49 Pattern formation Reaction diffusion equation is used to model natural phenomena of individuals moving and repro ducing in the habitat Models are formulated and computed and then we compare the simulated 49 PATTERN FORMATION 71 data with experimental data If a particular population distribution appears in the experimental data often then it may be of important value since this distribution could occur in future again The model based on the phenomenon is more reliable if the same density distribution can be gen erated by the model Such particular density distribution which persists and can be observed in nature or experiment is called a pattem Reaction diffusion equation is one of the simplest model which considers the spatial distribution of a substance or species Thus it is also thought as a basic mechanism of generating patterns by physicists chemists biologist and engineers Suppose there is a reaction diffusion model What is a pattern In natural science any observ able phenomena could be called a pattern which could be temporal spatial or spatial temporal Mathematically the de nition of pattern is not clear In this section we use a naive way to de ne it De nition 48 A pattern is a non constant stable equilibrium solution of a reaction diffusion equation For a reaction diffusion equation on a bounded domain this de nition seems quite reasonable because from our studies so far the solution of the equation always has a limit as t a 00 unless it goes to in nity In the case of Dirichlet problems the limit is either a constant u O which we do not consider as a pattern or an equilibrium solution which is a pattern In the case of Neumann problems the limit is always a constant for all numerical experiments we did so far In this section we con rm these observations mathematically by establishing the following results 1 The limit of a bounded solution of a reaction diffusion equation is always an equilibrium solution which is either a constant or a pattern 2 There is no pattern for the Neumann boundary value problem of reaction diffusion equation The rst result is similar to that of ordinary differential equation 1 For the ODE any solution is monotonic and if it is also bounded then it has to approach an equilibrium point For the reaction diffusion equation the solutions also have a certain monotonicity which can be best explained using the energy method The equation u can be viewed as describing a physical process with a potential energy 7Fu f ftdt Suppose that ut is a solution of the equation then the rate of change of energy along the solution orbit is d ElFutl F Utu t futu t flut S 0 491 Thus the potential energy decreases along an orbit and each solution will eventually fall to a local minimum of the energy function 7Fu which is equivalent to a zero point of 7fu 7F u with ifu 7F u gt 0 We can nd a similar formulation for the reaction diffusion equation Suppose that ut z is a 7 2 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS solution of 2 8 fltugt tgt0melto1gt at 8x2 74050 ut71 07 492 u0x u0z gt 07 m 6 01 The energy function is Eutz 01 lt uitx 7 Fun gt dag 493 where f ftdt Then the rate of change of energy along the solution orbit is 1 Eutx A Qumtmumt7 m 7 Fut7 xut7 dm 1 494 umt7mumtw 7 fltultmgtgtutltmgtgtdm 0 lntegral by parts again provides a big help 1 1 umt7 Mum t7 zdz umt7 mdut t7 m 0 0 1 uwwm w i Utt7mdtw 0 495 1 07 uttmumtmdm 0 1 7 uttzumtmdm 0 Here the boundary terms are zero since utt7 0 utt7 1 0 for the Dirichlet boundary condition From 494 and 4957 we obtain 1 Eut7 7 uttmumtm ut7 1 496 7 ut2txdz 0 0 Here in the second step we use the fact that u is a solution of ut um ful Therefore the energy is decreasing along an orbit utz For this reason7 sometimes is called a Lyapunov function If utz is also bounded7 then the energy Eut7 is also bounded7 and it must have a limit7 which in fact must be a critical point of the energy function The critical point here means a function ut7 x such that dEutmdt O From 4967 this critical point must satisfy fol ut2tmdz 07 thus it must be an equilibrium solution Thus we have proved though not completely rigorously Proposition 49 Suppose that utm is a bounded solution of 492 Then there exists an equilibrium solution 11m such that tlim utm 4 00 49 PATTERN FORMATION 73 This property of reaction diffusion equation tells us why the equilibrium solutions are so impor tant In fact the limit equilibrium solution is usually a local minimum point of the energy function which is equivalent to stable equilibrium solutions So stable equilibrium solutions patterns or sta ble equilibrium solutions are even more important However in the next result we show that there is no pattern in Neumann boundary value problems Theorem 410 Suppose that is a non constant solution of DUH 90 07 6 ELI7 43997 u0 uL 0 Then u is unstable Proof Since is a non constant solution on 0 L then must be one of the following cases a is increasing on 0 L b is decreasing on 0L or c is neither increasing nor decreasing lf is increasing on 01 then u z 2 0 for z E 0 L and 710 0 and uL 0 498 because of Neumann boundary condition Also u z 2 u 0 for any x gt 0 so u 0 2 0 is increasing at z 0 and similarly u L 0 Suppose that gt1z is the eigenfunction of the largest eigenvalue 1 then maz g ltultzgtgt 1 m z e M 410 altLgt 0 From Theorem 41 gt1z gt 0 for z E 0 L On the other hand if we differentiate 497 with respect to m we obtain 499 DoWm g ltultzgtgtu 07 x 6 lt07 4100 u0 uL 0 We apply the multiply multiply subtract integrate trick to 11 and u on 0 L then we obtain L L lu u gt1dxm gt1zu zdz 4101 0 0 The left hand side does not completely disappear In fact by using the boundary conditions we have L A 147 u gt1dz iu L gt1L u 0 gt10 2 0 4102 since u L 0 gt1L gt 0 u 0 2 0 and gt10 gt 0 On the other side L gt1xumdx gt 0 4103 0 because both gt1z and u z are positive Thus by comparison 1 gt 0 and is unstable 74 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS If is decreasing the proof above can also be carried out in an exact same way If is neither increasing nor decreasing then there exists a subinterval 12 of 0L such that u z1 u m2 0 u z1 2 0 u z2 0 and u z 2 0 in 12 Then we can apply the proof above to the interval 1 2 Hence is unstable in any case I 410 Traveling waves Let s consider the Neumann boundary problem 814 8214 EWU17U tgt0m 01 wcmwcna u0x u0x m 6 01 4104 From the results in last section the asymptotic behavior of the solution is pretty clear a the limit of any solution is an equilibrium solution b since any non constant solution and u 0 are all unstable equilibrium solutions then most likely the limit is u 1 which is stable However we can observe some interesting dynamics if we choose the initial function u0z to be a function which is almost 0 on the whole interval 01 but is positive near z 1 The particular function we choose here is uo 25 40m2 We could assume that this is a population concentrating near z 0 and completely absent on 05 1 see Figure 416 a Once the evolution according to the diffusive logistic equation begins the peak of the population quickly drop to near u 1 but the remaining part of the pro le hardly changes except the tail extends from x 04 to z 045 this is only approximation since we only takes 10 grid points Since then the shape of the solution does not change muchithe maximum is at z 0 near z 0 the solution is nearly at at u 1 the solution transits from u 1 to u 0 in an interval of length 03 then it is almost 0 on the subinterval near z 1 But the transition part which we call a transition layer moves as time elapse initially t 0001 the center of the transition layer is at z 024 Figure 416 b then it moves to z 043 at t 0005 Figure 416 c and to z 068 at t 001 Figure 416 Here we de ne the center of the transition layer to be the value of x such that 05 But after that the transition layer breaks off Figure 417 e and the solution converges to 1 quickly Figure 417 Biologically the evolution of the solution from t 0 to about t 002 simulates an invasion of a species from its original habitat to the whole region The pro le of the solution during the the invasion is almost same and it looks like a wave propagating from the left to the right Now we derive more analytical information about this wave propagating phenomenon We notice that during the invasion the transition layer moves at an almost constant velocity Suppose this velocity is c then at the point 251 1 and 252 m2 if zgtg zltl c then ut1 1 ut2 m2 So utz should be in a form of wm 7 ct We notice that during the invasion the impact of boundary condition is not important so we consider just the equation and on the whole space 814 8214 a W u17 u tgt 0 m E 700 00 4105 x 4 10 TRAVELING WAVES 7 5 2 2 1 1 0 0102 03 04 05 06 07 08 09 1 0 0102 03 04 05 06 07 08 09 1 2 2 1 1 0 0102 03 04 05 06 07 08 09 1 0 0102 03 04 05 06 07 08 09 1 Figure 416 Solution pro le a t 0 b t 0001 c t 0005 d t 0010 and we looks for a special solution with form ut z wz7ct We call such a solution a traveling wave solution Traveling wave solutions are similar to equilibrium solutions in fact equilibrium solutions can be regarded as traveling wave solution with velocity 0 If your observation point moves with the traveling solution at the same speed then in your eyes the traveling wave does not move just like an equilibrium solution Let y z 7 ct Then 2 wy 7c and w y 4106 Thus equation 4105 becomes an ordinary differential equation 7010 w w17 w 4107 but with constrains wy lt 0 lim wy 1 lim wy 0 4108 37700 117100 The equation can be converted to a planar system 10 U 1 7w17 w 7 31 4109 and the constrains now are U lt 07 hm wy7vy 1707 07039 439110 y7gt7oo 76 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS 0 0102 03 04 05 06 07 08 091 0 0102 03 04 05 06 07 08 091 Figure 417 Solution pro le 13 t 0015 f t 0020 By phase plane analysis we nd that 10 and 0 0 are the only equilibrium points of system 4109 thus the traveling wave solution is a connecting orbit between two equilibrium points The question is is there such an orbit and what is the velocity 0 We can achieve that by a little more phase plane analysis We linearize the system at both equilibrium points then the Jacobians are J00 BA EC J10 2 1 4111 7C Thus 10 is always a saddle point and the type of 0 0 depends on the values of and c The eigenvalues of J0 0 are 7 2 Ci4 4112 I lf 02 7 4 lt 0 then 0 0 is a spiral sink And we do not have an orbit satisfying 4110 since all solution approaching 00 are in an oscillatory fashion then vy cannot be always negative lf 02 7 4 2 0 then 0 0 is a sink and the unstable solution coming out of the saddle point 10 hits 00 directly see Fig 418 This solution satis es 4110 so uz 7 ct is a desired traveling wave solution For the existence of the traveling wave solutions we have the following result from phase plane analysis above Theorem 411 4105 has a family of traveling wave solutions uctz c 2 4 and for each 25 E 700 00 lim utm 1 lim utm 0 4113 17700 17700 In biological applications an important question is the speed of the population wave From Theorem 411 we nd that the minimum speed is c 2xX But in the application which speed does the population prefer Let s go back to the example in the beginning From the coordinate of the center of the wave we can estimate that the wave speed is about 0 z 52 and the theoretical minimum wave speed from Theorem 411 is c 14000 z 63 since 1000 Considering the calculation errors and the impact of the boundary conditions we can guess that this population is invading at the minimum speed c Indeed this fact was proved in a pioneer paper by Kolmogolov Petrovki and Pisconov in 1936 411 STANDING WAVE 77 W v v qwliwicv w Figure 418 Phase portrait of traveling wave We could compare this to the diffusive Malthus equation 814 8214 a y 1 Au tgt 0 m 6 70000 4114 In Section 34 we have obtain that the front of the population wave approximately advances at the speed R t v 4 which is same as the minimum wave speed of the logistic case lntuitively the wave speed of logistic equation should always be smaller than the Malthus case the growth rate is smaller here we show that asymptotically the two speeds are same But in the Malthus case population has an exponential growth in the metro area and the population tends to the carrying capacity in the center part for the logistic equation The traveling wave solution here can be used to model the biological invasion of a foreign species into a new territory 411 Standing Wave 412 Summary In this chapter we studied the behavior of the solution of reaction diffusion equation 8U 82u EWfu tgt0 m6lt0717 w 0 ut10r 1417570 M221 07 4115 u0x 6 07139 Mathematically our investigation can be summarized as follows H to CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS For a solution utz if it remains bounded then it approaches an equilibrium solution of the equation and most likely this equilibrium solution is a stable one For Dirichlet problem the stable equilibrium solution is either a positive one or zero for Neumann problem the stable equilibrium solution is a positive constant Biologically we can observe several phenomena for different growth rates different diffusion constants different size of the domains or even different initial values H to DJ to 03 F Extinction the population approaches an equilibrium zero Persistence the population approaches a positive equilibrium solution Invasion the population invades the new territory at a more or less constant velocity Chapter 4 Exercises Use the substitution sat y4x U 4116 to deduce a new equation from utDumault1igt tgt0 m 0L ut0 ut L 0 4117 M0790 1 m E 07L Express the new parameters and initial values in terms of old ones Find a substitution different from 43 and 4116 to deduce a nondimensionalized equation from 4117 Consider diffusive population model with Allee effect utDumault1igt 71 tgt0 x 0L ut0 ut L 0 1407 90 1 4118 m E OL Use a table to list the dimensions of all quantities and parameters and convert the equation to a dimensionless equation Modify the numerical algorithm in Section 42 for non homogeneous boundary condition umt0 cl umt 1 62 4119 U a T 00 to H O SUMMARY 79 Modify the numerical algorithm in Section 42 for Robin boundary condition u t 0 but 0 umt 1 7but 1 b gt 0 4120 Write a Maple program to implement your algorithm Write a Maple program to numerically compute the solutions of 8 82 u UAu1iu tgt0 x6 01 at 8x2 4 121 4mmmna ltgt MQMfL memA Use positive initial data x gt 0 Con rm that for all positive initial data and all A gt 0 which you try the solution will always tend to u 1 Write a Maple program to numerically compute the solutions of 8 u Au17uu7005 tgt0 x 6 01 at 8x2 w 0 ut1 0 4122 MQMfL memA Use the following parameters a A 20 any x gt O b A 80 x O3sin7rx x 06 sin7rx and x sinp x Describe the qualitative behavior of the solutions For the Dirichlet problem of diffusive Malthus equation 814 8214 Au tgt07m 0717 8t 8 2 ummfwmww Mm Max x 6 01 What is the stability of u 0 when A 712 Also show that sin7rx is also an equilibrium solution when A 712 What is the stability of this equilibrium solution Use the technique in the proof of Proposition 42 to show that when A g 712 4121 has no positive equilibrium solution Hint suppose there is one say ux and use the multiply multiply subtract integrate trick to and sin7rx remember that v sin7rx satis es U 7r 110 Consider the equation x6x2x76x0 4124 For 6 near 6 3 the equation has a solution with a form x5 1rx1r2x2higher order terms where r 6 7 3 a Use perturbation method to calculate x1 and x2 b Use the formula in implicit function theorem to calculate x1 80 11 H to H 03 H F H U H a H 1 CHAPTER 4 NONLINEAR REACTION DIFFUSION EQUATIONS Consider u Au17 u 0 m 6 01 u0 u1 0 4125 A u 97r2 0 is also a bifurcation point Use the perturbation method to nd the approx imate form of solutions near the bifurcation point A 971392 1 ale 1 0282 1 6u1x 62u2z 1 higher order terms 4126 Show that u1z k sin37rm for some constant k and determine k by solving the equation of U2 and assuming a1 1 Consider again 4125 but at A 471392 The perturbation method which works at A 71392 and A 971392 does not work exactly same way at A 4712 By using the same perturbation method in 4126 with 971392 replaced by 471392 show that an 0 and ul 0 Use Maple to further expand the series show that a2 gt 0 and assuming that a2 1 calculate the rst nonzero 71 Use perturbation method to calculate the approximate solutions of u l Au17 u 0 m 6 01 uO u1 0 4127 at a bifurcation point Au 7r2 0 Prove the following Fredholm alternative theorem by using the multiply multiply subtract integrate77 technique suppose that y Wzy f7 240 07 241 07 4128 has a solution ym then 01 fz sin7rmdz must equal to zero Consider 81 821 w EWAu1iui 1bu tgt0 z 01 4129 ut0 ut 1 0 140796 f907 96 6 071 where a b gt 0 This equation models a species which is harvested or predated and the rate of harvesting has a saturation limit as u a 00 which is also called Holling s type II functional response Write a Maple program to numerically compute the solutions of 4129 Use positive initial data fz gt O Show that for any A gt O u 1 is an asymptotically stable equilibrium solution of 427 and u 0 is an unstable equilibrium solution of 427 1 Show that if u c is an unstable stable equilibrium of 1 fu then c is an unstable asymptotically stable equilibrium of 814 8214 a i W f 10007 wmmwcna u0x m 6 01 tgt0 m6 01 4130 412 18 H to SUMMARY 81 a Use Maple program to draw the bifurcation diagrams of the equilibrium solutions of 483 when the parameter is the following values 1 a 04 2 a a 7 4 a 702 5 a 71 6 a 710 4131 We have known that a 05 is a bifurcation point where the bifurcation diagram has a qualitative change from no positive solutions to exist positive solutions From the various bifurcation diagrams7 guess other bifurcation points b When a lt 0 weak Allee effect7 the positive solutions may bifurcate from the zero solutions for some a Determine for which 017 this is possible7 and where is the bifurcation point a for each possible a For which a the bifurcation is supercritical7 and subcritical c Explain the biological signi cance for each different bifurcation diagram When a type ll functional response is added to the diffusive logistic equation7 the new equation is 82u au 8u 7 u17u71bu7 54 tgt0x 01 ut0 ut1 O7 u0x gt 07 m 6 01 4132 a Determine for which values of 017 b that the growth function is of strong Allee effect7 or weak Allee effect b Choose two sample parameter sets 017 b which are of strong Allee effect7 and weak Allee effect respectively7 and use Maple to draw the corresponding bifurcation diagrams of equilibrium solutions Consider 2 8u 3 u EW ffuv tgt07 6 ioovoob where fu u when 0 lt u lt12 and fu17 u when 12 lt u lt1 Solve the exact traveling wave solution ut7 z 12z 7 ct 4133 Consider 8 82 8 9 7 Zbu7u37 tgt07 m6 70000 This equation has a traveling wave solution with speed c 0 Find this solution by solving the corresponding ordinary differential equation 4134 Consider 2 8u 3 u 2 a 7 W u 17 u7 tgt 07 m 6 70000 4135 a rite t e or 1nary 1 erent1a equation sat1s e y t e trave1ng wave equation 0 W h d d ff l h d b h l f 4114 b Find the range of speed c for which 4114 has a traveling wave with that speed Chapter 2 Diffusion equations 21 Separation of variables intervals Diffusion equation is a linear partial differential equation7 since the functions related to u in the equations ut and Au are both linear Recall for the linear ordinary differential equations 2 kP and DE kP 21 13 it is well known that 5 and Ac0s kth Bsin kDz are the solutions of the respective equations Indeed in the second case7 the exponential function skimD is a complex solution Thus we can make a lucky guess for the solution of the diffusion equation 8P 82F a Tm 2392 that the solution is an exponential function Pt z swim 23 for some constants 11 which could be complex numbers By substitution7 we nd that the function in 23 is a solution of 22 if a D1927 and depending on the values of b7 we obtain three families of solutions Ptz eDbZtebm b gt 0 24 Ptz eDV e bw b gt 0 25 Ptz 513172517 b gt 0 26 The rst two families of solutions are not reasonable solutions7 since when x a 00 or foo Pt7 z A 007 which implies an unlimited growth at z 00 or foo The third family of solutions are complex7 but their real and imaginary parts are both solutions of 22 too Pt x 5 13th cosbm and Pt x 5 13th sinbx 27 15 16 CHAPTER 2 DIFFUSION EQUATIONS These solutions are little more reasonable as they are bounded as x a 00 but still they do not satisfy natural boundary conditions P Pm a 0 as x a 00 Nevertheless solutions with above forms are solutions of diffusion equations and we notice that they are in form of Ptz UtVz This motivates us to consider the solutions which is separable on the two variables 25 and z utz UtVz 28 With the form of u in 28 the diffusion equation 22 becomes U t VW U WDU Wam Dk Ut 7 We 7 where k is a constant independent of both if and z Thus U satis es an equation U t DkUt 25gt 0 29 and Ut 015W 210 and V satis es V z 211 The solutions of 211 are Cge m with any complex number k Thus we obtain solution of 22 CeDkte m 212 just like the three families of solutions we obtain by making lucky guess To obtain more speci c solutions we must take the boundary and initial conditions into considerations We start with the one dimensional diffusion equation on an interval with length L 8P 7 82F 3t 7 8x2 7 Here the individuals of the the species inhabits a one dimensional patch 0 L and the diffusion coef cient of the population is D gt 0 To get a de nite solution of the problem we need to have an initial population distribution 35 e 0 L 213 P mean no and a boundary condition We assume that the exterior is hostile so homogeneous Dirichlet bound ary condition is satis ed Pt0 Pt L 0 215 We look for solution with form Pt z UtVz Then we obtain 29 and 211 and 215 implies that V0 VL 0 216 211 and 216 become a boundary value problem for an ordinary differential equation which is a special case of Sturm Lz39oum39lle problems We shall discuss this kind of problems in a separate section next 22 EIGENVALUES AND EIGENFUNCTIONS 17 22 Eigenvalues and eigenfunctions of boundary value problems We consider the Dirichlet boundary value problem 2 kg m e M ylt0gt yltLgt 0 217 We can see that 0 is always a solution of the equation So we shall look for solutions other than the trivial zero solution From elementary ODE knowledge we know that the general solution of y kg is one of the following three forms clam mm if k gt 0 218 01 022 if k 0 219 01 cos7kt 02 sin7kt if k lt O 220 One important characteristic of a boundary value problem is that it may not have a solution while the initial value problem always has a solution and the solution is unique if all terms in the equation are differentiable For 217 there is no solution when k gt 0 Indeed if there is a solution it must be yt cle t cge t for some constants 01 and 02 and we can solve the constants by 0 y0 cl 1 2 0 yL cle L Cge L But this system of equations is unsolvable if one gets 01 702 from the rst equation and plugs it into the second one one has e L LCz O and one must have 62 0 and 01 0 Thus 217 has only zero solution when k gt 0 It is also easy to check that 217 has only zero solution when k 0 So the only meaningful case is when k lt O 75 In this case 217 may still have only the zero solution but for some number k lt 0 it can have a solution yt which is not the zero solution Note that from the linear principle if yt is a solution of 217 so is 0 yt for any constant c and if y1t and y2t are solutions of 217 so is y1t y2t So for k which 217 has non zero solutions we call these numbers eigenvalues and the solutions of 217 eigenfunctions In fact historically the solutions of 217 were rst called eigenvalues and eigenfunctions Now let s search for the eigenvalues and eigenfunctions of 217 the eigenfunction must be yt 01 coshfist 62 sin7kt for some constants 01 and 02 by the boundary conditions we ave 0 y0 01 0 yL 01 cos7kL 6281H7kL so 01 0 and if 02 is not zero sinx7kL 0 thus xikL mm the eigenvalue k must be m27r2 k 77 221 2 t and the eigenfunctions associated with k imL are csin Therefore 217 has a mz 392 7 L2 sequence of eigenvalues km 7 and the set of eigenfunctions are the sine functions which 18 CHAPTER 2 DIFFUSION EQUATIONS can t in77 the interval 0 L When the interval is 0 7r so L 7139 the expression of eigen pairs is much simpler km 7m2 fmt sinmt 222 Back to 213 and 215 we nd a sequence of solutions 2 2 umtm exp lt7DmL t sin m is a positive integer 223 Because of the linear principle the general solution of 213 is 00 mz 392 7mm 142532 gcmemp 71372 s1n 25gt 0 m e 0 L 224 Example 21 y kg 0 220 we yltwgt 0 225 This problem is same as previous example except the boundary conditions So we still consider the cases of k lt O k 0 and k gt 0 When k lt 0 yt cle m Cge t from the boundary conditions we get y0 y 0 01 02 77 k01 x7 kcg 0 and y7r 1equot 1W 02 0 So 01 02 satis es the equations i ik 1x7k 01 7 0 226 e7x7 k7r ex7 k7r CZ 7 0 39 i i i 1 7 7k 1 7k i The equation has non zero solution only if det lt sixj EH O The determinant equals to 17 pe 7r 7 1 pe pquot if we denote p x7k thus det 0 is equivalent to 52pquot The 7 P 1 functions f1 p 52p and f2p 2 have no common points when p gt 0 in fact they only 1729 HP intersect at p 0 thus 52 quot is not possible for p gt O and there is no eigenvalue such that k lt 0 When k O yt 01 0225 from the boundary conditions we get y0 y 0 01 02 0 and y7r 01 027r 0 since the determinant of the matrix lt1 W is not zero then only solution for 01 02 is 00 So 0 is not an eigenvalue either For k gt O yt 01 cosh t 62 sinEt The boundary conditions y0 y 0 01 1 02 O y7r 01 coshEr 02 sinE7r O The matrix is MiEr sinkgwgt and the determinant of the matrix is 0 if sinE7r 7 xEcosE7r O or tan7rp p for p The graphs of f1p tan7rp and f2p p intersect at a sequence of points pm m 1 23 and 05 lt p1 lt 15 15 lt p2 lt 25 m 7 05 lt pm lt m 05 km p3 are the eigenvalues of the problem 225 and the eigenfunctions are ymt McosQEt 7 sinEt 23 FOURIER SERIES 19 23 Fourier series To completely solve 213 215 we need to determine the constants cm in 224 by setting t 0 in 224 00 u w 0 sin m 26 OL 227 olt gt m L lt gt lt gt The expression in the right hand side of 224 is the Fourier series of u0z on 0 L Fourier series is widely used in mathematical physics and signal process All sine and cosine functions are considered to be regular smooth waves and Fourier series convert any noise or signal into a sum of more regular waves and that is a big help for the signal processing Here we collect a few facts of Fourier series Theorem 22 Let and be piecewise continuous functions on the interval 0L Then can be expanded in either a pure sine series 102 bmsin 228 2 L where bm sin dz m 1 2 or a pure cosine series 0 102 22am cos 229 m n w 2 L where am ZJ fz cos lt L gtdxm012 The Fourier series is similar to the Taylor series n 1 2 f 090 230 m0 39 71 but instead of expanding the function in polynomials Fourier series expands a function in sine or cosine functions But Taylor series converges to the original function only in a small neighborhood of z 0 though some time it is large and the neighborhood depends on the function f Fourier series converges to the original functions for all functions and on the whole interval 0 L But it may not converge to fz at z 0 and z L Theorem 22 shows all sound wave can be decomposed into several regular waves plus a small noise The approximation of Fourier series is sine or cosine Fourier polynomials N N fzzmz1bmsinltmmgt fx amCOSlt m n w a 231 20 CHAPTER 2 DIFFUSION EQUATIONS Now we can complete the solution formula of 213 215 00 mz 392 7mm utz Z cmezp lt7D L2 25 s1nlt L gt quotF1 L 232 2 7mm gt where cm LO u0z s1n lt L gtdmm 71 25gt O z E 0 L Example 23 utummy tgt0 07r ut 0 ut 7r 0 233 u0z sin m From the solution formula 232 we get 2 2 W utz Z cme m t sinmz and cm sinzsinmzdm 234 7quot 0 m1 7r 7r However sinzsinmzdz 0 for any m 2 2 and sin2 mdz 7r2 Therefore the solution is 0 simply ut z 5 sin m 235 Example 24 For the homogeneous Neumann boundary value problem ut Du tgt Om E 0L umt0 umt L 0 236 M07 95 14095 the solution is 2 0 2 ut z 630 Z cmezp lt7DmL t cos m1 77mm 2 L where cm Zj u0z cos lt L 237 gtdz mZO tgt0 x 0L The derivation of this formula is left as exercise For these examples the asymptotic behavior of the solution can be read from the formulas For homogeneous Dirichlet problem 213 215 the asymptotic limit is i i i 7r2t i 7m lim utz O and the asymptotic pro le is utz N clamp 7 2 sin 238 taco L L On the other hand for homogeneous Neumann problem 236 the asymptotic limit is 2t ut z g and the asymptotic pro le is ut z N 630 clamp cos 239 24 RELAXATION TO EQUILIBRIUM SOLUTIONS 21 L The constant co above is 2L u0zdm thus the limit is 0 L A mam 240 the average value of the initial value u0m The limits of both cases are also solutions of diffusion equation and in fact they are solutions which are independent of time 25 Such solutions are equilibrium solutions For a general reaction diffusion equation 8 8 dAufmu tgt0 m E 9 Mom 140907 m e n 241 Bu0 tgt0m 8 2 where Bu is the boundary condition the equilibrium solutions satisfy dAu 7 0 m E Q 242 Bu 0 m E 89 A solution u1z of 242 is a solution of 241 in the sense that the function U1t z u1z is a solution of 241 with u0 z It is not just a coincidence that for two examples we have considered so far the limit is always an equilibrium solution We will show in next section that the limit of a diffusion equation is usually an equilibrium solution Equilibrium solutions of diffusion equation are easier to nd and that will help us understand the limit of solution when the solution itself is not so easy to nd The equilibrium equation for 213 215 is u 0 u0 uL 0 243 The general solution for u 0 is ax b and from the boundary conditions we obtain 0 is the only equilibrium solution The equilibrium equation for 236 is u 0 uO uL 0 244 For this problem b for any I E R is an equilibrium solution However the limit in 236 is the average value of the initial value u0z only Why does the solution of diffusion equation selects this particular equilibrium solution but not the others We will answer this question in next section 24 Relaxation to equilibrium solutions We have obtained analytic solutions of 1 d diffusion equation in previous sections The solutions are in a form of in nite series However in general even such formulas cannot be obtained when the spatial dimension is higher and the domain is arbitrary Later we will formulas of solutions of diffusion equations for some special cases Here we use some integral formulas to show some qualitative properties of solutions of general diffusion equation 22 CHAPTER 2 DIFFUSION EQUATIONS As a calculus warmup we consider the diffusion equation with no ux boundary condition jag 7DAu tgt0z 2 74071 was x e o 245 Vutm nx 0 t gt 0 m E 89 From the biological modelling point of view the species only randomly moves in Q the growth rate is zero and the boundary is sealed Thus the total population should not change since no one moves in or out and no one is born or dead lndeed that is the basic character of the equation and it also explains why the limit of the solution of 236 is the average value of the initial value function Proposition 25 Suppose that utm is the solution of 245 then utzdz u0mdm 246 9 9 Proof Integrating the left hand side of the equation in 245 over 9 we obtain 9 tmdm ltQutmdmgt 247 and integrating the right hand side we have D Autzdm D was x nmds 0 248 0 39 from the divergent theorem and the no ux boundary condition Therefore from 247 and 248 we nd that d E utxdmgt o 249 which implies 246 III Next we show that the solution ut z always tends to a limit constant c and because of 246 that constant must be 1 uo m dx 250 m 9 lt gt lt gt A rigorous proof of this fact needs many more advanced mathematical knowledge so we only sketch the ideas here We multiply the equation in 245 by the function ut m 1314 U5 7 DuAu 251 Similar to the proof of Proposition 25 we integrate the left hand side of 251 ea 7 a 12 7d 12 71d 2 Quatdmi at lt2ugtdmidt Qlt2ugtdmi2dt lutml dx 252 24 RELAXATION TO EQUILIBRIUM SOLUTIONS 23 and the right hand side D uAudx D uVu nds 7 D Vu Vudx 253 D 39 9 Here we use Green s identity a consequence of divergence theorem uAudx uVu nds 7 Vu Vudx 254 D 39 D From the boundary condition fanuVu nds 0 thus from 252 and 253 we nd 1 d ut m2dm 7D qut ml2dz 255 26h Q 9 Let Et ut m2dz Then 255 shows that E t lt 0 for all t gt 0 and obviously Et gt 0 Q for all t gt 0 Thus limtH00 Et 2 0 exists and as t a 00 E t a 0 so i 2 7 qutml dm 7 0 256 256 implies limtH00 qut 0 for any x E 9 thus the limit of utz when t a 00 must be a function with zero gradient thus the limit function must be a constant 0 Since elm tlim utzdm u0zdm then c Notice that we 9 4 00 Q Q implicitly assume utz has a limit when t a 00 While such a limit indeed exists the proof of its existence is beyond the scope of this notes that is why this is only a sketch of proof Nevertheless we put our ndings into another proposition Proposition 26 Suppose that utm is the solution of 245 then utm tends to a constant xdx function W uo 9 Similar result can also be shown for Dirichlet problem in which the limit function is u O From this property of diffusion equation we can get a conclusion that the diffusion equation will make the initial value atter and bring it eventually to an equilibrium state without spatial pattern This can also be indicated from the differential equation itself Let the average of the initial value be c For the points where ut z gt c the function is likely to be convex down with Au lt 0 thus the function value will decrease over these points and for the points where ut z lt c the opposite occurs This is called smothering effect of diffusion In the following we use a Maple program to demonstrate the smothering effect and the asymptotic tendency to the equilibrium We consider ut 4um t gt 0 m 6 04 ut0 ut 4 0 257 u0z 24 CHAPTER 2 DIFFUSION EQUATIONS 100 x 6 13 From 232 we get the 0 x E 01U34l where x is given by the formula x formula of the solution of 257 7 00 m27r2 m7rx utx 7 gcmexp lt7 4 t sin 4 258 1 where cm 5 xsin dxm 21 25gt 0 x 6 04 0 The coef cients cm in the Fourier series can easily be calculated but we will let Maple to do that Here is the Maple program diffusionmws and some of the results gt restart gt withplots n n gt fx gt100Heavisidexl100Heavisidex3 gt plotfx x04 gt cn gtO5intfxsinnPix4x04 Cn cos14 n Pi cos34 n Pi 2000 n Pi gtuapprox1xt gtsumcnsinnPix4exp4nquot2t n1 30 30 uapprox1 xt gt Z cn sin14 n Pi x exp 4 n2 t n1 gt tvals seq0 0720i i0 20 gt toplotsequapprox1xt ttvals gt plottoplotx0 4 gt densityplotuapprox1x t x0 4 t0 O 2 colorstyleHUE gt animateuapprox1xt x0 4 t0 O 1 frames50 This is an animation generated by approximating the solution by the partial sum with n 30 of the series solution The approximation is accurate since the remainder is an exponential decaying function with smaller exponent Note that the initial condition x can be written as the difference of two Heaviside functions which is de ned as 1 x Z 0 0 m lt 0 259 It is an interesting experiment to observe the changes of the spatial patterns of ut x when 25 gets forward Here are a few snap shots of the animation 24 RELAXATION TO EQUILIBRIUM SOLUTIONS 25 Figure 21 a fm b partial sum of Fourier series of fz 100 100 80 80 so so 40 40 2o 20 0 1 g 3 4 0 1 g 3 4 Figure 22 a u0002m b u0010z Because of the truncation of the series solution utz above indeed is the solution of 257 with initial data 2 cm sin dm cm cos 7 cos 260 4 m H Figure 21 is a comparison of f and From the last section we know that when t a 00 ut z a 0 exponentially fast However the road to distinction is worth a more careful lookingiit can be dismantled into the following stages 1 Forming a plateau From 25 0 Figure 21b tot 0002 Figure 22 a we can observe that the initial data has been smoothed by the diffusion equation and all the small wiggles in u0z are gone at t 0002 On the other hand although u0z is not a constant for z 6 13 u0002z is almost a constant on an interval slightly smaller than 13 So in the initial stage of the evolution of 1 under diffusion a plateau or at core of height 100 forms with a sharp but smooth drop off to 0 near z 1 and z 3 CHAPTER 2 DIFFUSION EQUATIONS Figure 2 x 03 q Figure 24 a u0100m b u0600z Erosion of the plateau From 25 0002 Figure 22 a tot 001 Figure 22 b the at core keeps shrinking though the top of the core is still at the level near 100 At t 001 the erosion nally dissolves the whole at top and the maximum point of ut at z 2 looks more like a regular local maximum point The graph of u001 z is close to a bell shape In that time interval the in ection points of ut are still near z 1 and z 3 but at t 001 the interface between u 0 and u 100 is no longer as sharp as when t 0002 Bell getting round In the next stage smoothing effect brings the concave part of the graph down and the convex part of the graph up And at the same time the convex part is shrinking as population loss via boundary gets bigger At t 004 the graph becomes almost concave for all z and the bell shape becomes an arch see Figure 23 Collapse of the arch In the nal stage the arch is close to the leading term of the series clezp77r2t4 sin7rz4 and an exponential collapse is obvious After 25 1 the graph can hardly be seen see Figure 24 25 ONE DIMENSIONAL CHEMICAL MIXING PROBLEM 27 25 One dimensional chemical mixing problem Here we consider a simpli ed version of the chemical mixing problem which is discussed in Chapter 1 A tube with length L contains salt water The cross section yz direction of the tube is so small so we can assume that the concentration of the salt water is same for any point on a cross section Let the concentration of the salt water in the tube be ctm kgms t gt 0 and 0 lt z lt L We assume that the diffusion of the salt is onedimensional in the tube and the velocity of uid is ignorable thus 005 z satis es the diffusion equation 36 82 E 261 At the left hand side of the tube a salt water solution with constant concentration C0 kgm3 enters the tube at a rate of R0 77138 and on the right hand side of the tube the mixed solution is removed at the rate of R0 mSs We also assume that the area of the cross section of the tube is A Let A1 and A2 be the cross sections at z 0 and z L correspondingly From the Fick s law the total in ux at z 0 is J inds 713mg 0 1ds iDAcmt0 00120 262 A1 A1 and the total out ux at z L is J nds 713mg L 1ds iDAcmt 0 CtLR0 263 A2 A1 thus we obtain the boundary conditions 8675 0 CoRo 8675 L 7 7675 LR0 7 d 7 264 8m DA an 8m DA To simplify the notations we de ne 7 CoRO 7 R0 B7 DAandE7DA 265 Thus the initial boundary value problem is 2 D tgt0 ze 0L 13 80t 0 8025 L 2 66 7B E t L 0 39 8 8 clt gt ctmcoz tgt0 m 0L Before attempting to solve the equation we could nd the equilibrium solutions which satisfy equation 82 138 0 m e 0L 13 860 7 B 80L 8m 8m 267 EcL 0 28 CHAPTER 2 DIFFUSION EQUATIONS From a simple calculation we nd the unique equilibrium solution B CORO 01m7EBL7m 700 which is a linear positive function with negative slope on 0L and 01L Co To solve the equation we de ne bt z Ct z 7 01z Then bt z satis es L 7 x 268 2 D8 I tgt0 z 0L w www bt 00m 7 ciz tgt 0 m E 0 L 269 The reason we consider 269 instead of 266 is that the boundary conditions in 269 is homoge neous thus we can use the technique of separation of variables We assume that bt z UtVz Then Ut V W DVz Dk 270 Ut 5D and Vm satis es V 7 kV 0 m e 0L V 0 0 V L EVL 0 271 This is a half Neumann and half Robin boundary condition We use the method in Section 22 to solve 271 If k gt 0 then Vm 1equotEm cge m From the boundary conditions we get 77 7E detlt Ef V L E e L gt 0 272 and k must satisfy 27H 7 E 5 7 xE E But the function on the left is greater than 1 when k gt O and the one on the right is always less than 1 when k gt 0 So 273 has no solution thus there are no positive eigenvalues Similarly k 0 is not eigenvalue either 273 If k lt 0 Vm 01 cosx7km 2 sinx7km From the boundary conditions we get 2 0 and 01cosx7kL 7 x7ksinx7kL 0 Thus k satis es cot7 kL 77 17 274 From the graph of the functions f1 cotLm and f2z m there are exactly one intersection point mm of f1 and f2 in m717rLm7rL for any positive integer m gt 0 Then this km 7z is the m th eigenvalue of 271 and the corresponding eigenfunction is cos Wm Thus the solution of 269 is m717r M 00 bt z Z bme Dmgnt cosmmz where mm 6 L L satis es cotsz mm 275 m1 26 CRITICAL PATCH SIZE AND BlFURCATIONS 29 Hence the solution is the sum of the equilibrium solution and an exponential decaying Fourier series and the solution converges to the equilibrium solution exponentially fast It shows that in equilibrium state the exiting solution also has the concentration CO same as the incoming solution thus the amounts of salt entering and exiting the tube are balanced but the concentration of solution in the tube is higher than CO and the concentration at z 0 is the highest The example in this section can be thought as a model of a continuously polluted river The left end point z 0 is the source of the pollution where the pollutant enters the river at a constant rate and the right end point z L is where the river merges to a bigger river or the ocean Our mathematical result shows that the concentration of pollutant in a continuously polluted river is even higher than that in the source As a consequence we can see that at the con uence of two polluted rivers the pollution can be most serious since the concentration there is much higher than the ones of either branches We should be cautious that we ignore the effect of drifting here which could make a big difference for rivers 26 Critical patch size and bifurcations The simplest population growth model is the Malthus model dPdt aP for a positive constant a gt 0 In this section we consider Dirichlet boundary value problem of diffusive Malthus equation 2 8 13 Da PaP tgt0 m 0L at 8x2 ut 0 ut L 0 276 u0x u0x m E 0 L We have seen that population are at a loss from the hostile boundary but on the other hand the population will increase through linear reproduction So the loss and gain are in a competition Let s see what factor will determine the outcome of this battle Equation 276 can be solved in a similar way as 213 215 In fact the form of the solution only changes slightly adding a term 5 to re ect the exponential growth The solution of 276 is co 2 2 170595 Z cmezp at 7 DmL t sin 7 m1 L where cm u0zsin dmm 21 25gt 0 m e 0 L 0 277 However the fate of the solution may be different since at is added to the exponent now the exponent may be a positive one then ut x will have a part which is exponentially increasing In fact when the index m in the sum increases the exponent gets more negative thus whether ut z is exponentially increasing is determined by the rst exponent coef cient a 7 D39n39ZLZ We have D712 D i whenagt 7 orLgt 471 tlimou m foo 2 278 D7r D i whenalt70rLltHZ7n tlimoutz70 30 CHAPTER 2 DIFFUSION EQUATIONS If we assume that the growth rate a and D are intrinsic parameters determined by the nature of species and environment then the size of the habitat L will play the deciding role here The number L0 Dmr is called the critical patch size as the population cannot survive if the habitat is too small the outgoing ux wins over the growth but the population will not only survive but thrive to an exponential growth if the habitat is large enough growth outpaces the emigration We say that a bifurcation occurs at L L0 as the equilibrium solution u 0 changes stability type from stable for L lt L0 to unstable for L gt L0 We notice that at L L0 the equation has other equilibrium solutions sin7rmL the eigenfunction associated with rst eigenvalue k1 77r2L2 We will have more discussion on bifurcation problems when the growth rate is nonlinear From 278 we can also use a or D as bifurcation parameter instead of L sometime that is more convenient since we do not need to change the domain but we get equivalent results The concept of critical patch size is related to habitat fragmentation in ecological studies Habitat fragmentation1 is the breaking up of a continuous habitat ecosystem or land use type into smaller fragments which is considered to be one of several spatial processes in land transfor mation It is commonly used in relation to the fragmentation of forests Habitat fragmentation is mainly caused by human activities such as logging conversion of forests into agricultural areas and suburbanization but can also be caused by natural processes such as re From our mathematical result in this section if the the fragmentation of habitat limits the spatial movement of the plants and animals then the species can become extinct 27 Separation of variables rectangles Most of this chapter we have exclusively considered one dimensional spatial domain for the sim plicity of the mathematical analysis However any interesting application happens in 2 D or 3 D spaces In this section we discuss the separation of variables in a two dimensional rectangle We rst consider a linear diffusion reaction equation 8U 82 LL 82 LL a ltW8 y2gtu7 tgt07 79ER07 1X0757 Vun0 my 8R Mommy uo7y7 any 6 R 279 where D a b gt O and 8B is the boundary of rectangle R From the method of separation of variables we assume that ut my UtVmy 280 Then similar to onedimensional case we obtain U t Dk Ut 281 1From TEMS Terrestrial Ecosystem Monitoring Sites httpwwwfaoorggtostems 27 SEPARATION OF VARIABLES RECTANGLES 31 and 82 82 V V kV 25 0 R 0 0 1 WWW gt we mm gt 2827 VVn07 my 68R For 2827 we need a further separation of variables Waxy WZy7 283 Then we obtain that W77 Z77 z y 2847 k7 WW Zy so both and Z yZy must be constants On the other hand7 the Neumann boundary condition implies that W 0 W a 0 and Z 0 Z b 0 Therefore W and Z satisfy W x k1Wm7 m E 07a7 WO Wa 07 285 Z y k2Zy7 y E 075 ZO Zb 07 286 and k k1 kg 287 Both 285 and 286 are familiar eigenvalue problems in one dimensional space which have been studied in Chapter 2 The eigenvalues and eigenfunctions of 285 are 774271392 n n39m km 777 27 W748 cos 7 n 071727 7 288 and the eigenvalues and eigenfunctions of 286 are m27r2 2 m n39y kzmi 7Zmycoslt m071727 289 Therefore the eigenvalues and eigenfunctions of 282 are 71271392 m27r2 kmm 7 2 7 b2 7 Vnmz7y cos cos 7 290 where n7m 07 17 27 Thus we nd that the solution of 279 is ut7z7y i cmmemkmw w cos cos 7 291 n0m0 where 0mm can be determined by the initial conditions The spatial patterns of the eigenfunctions for rectangle are much more complex than those of onedimensional Here we consider a square R 077139 x 077r7 then the rst ve distinctive eigenvalues are see the contour graph of eigenfunctions below k00 07 k01 k10 17 km 27 k20 02 47 k2 k1 57 292 32 CHAPTER 2 DIFFUSION EQUATIONS Figure 25 From left to right Figure 26 From left to right a V0 cos2y b V10 cos2z c V12 cosz cos2y d V11 cos2z cosy We notice that many eigenvalues km have multiplicity more than 1 for example km kmm 7m2 712 In such case the eigenspace is of higher dimensional so more possible spatial patterns are generated For example km k0gt1 71 and the eigenspace is spancosmcosy thus k1 cosz k2 cosy is a spatial pattern for any k1 and k2 In particular the patterns generated by cosz i cosy are symmetric with respect to one of the diagonal lines Figure 27 From left to right a cosz cos y b cosy 7 cos m The solution 291 is always an exponential growth if A gt 0 Indeed from the equation the no ux boundary condition prevents the emigration through the boundary and gt 0 implies a positive growth rate thus the total population in R has an exponential growth We can consider corresponding Dirichlet boundary problem 8 82 82 8 7 Wye QM tgt07ltz7ygteRlto7agtxlto7bgt7 Mammy 07 any 68R u0my u0zy my 6 R 293 27 SEPARATION OF VARIABLES RECTANGLES 33 Similar to the method above the solution is 00 71713913 m7r w z y Z bn ne mmWt sin sin 294 a n1m1 where bmm can be determined by the initial conditions and kmm is de ned in 290 Similar to last section the persistenceextinction of the population can be determined h Lgtii 1 t 7 wen D712 a2 b2tir 10u 700 295 h Alt 11 t0 wenm a Z tir 10u 7 Chapter 2 Exercises 1 Find the and f t of Neumann boundary problem 24 247 95 6 07M 20 yL 0 296 2 Find the eigenvalues and eigenfunctions of periodic boundary problem 2 kg m e 0 gm m 210 W 297 3 Find the eigenvalues and eigenfunctions of Robin boundary problem I gt 0 y 1 95 6 07 MW by07 2L byL 298 4 Find the eigenvalues and eigenfunctions of the problem y 1 231 kg 0 m E 07r y0 y7r 0 299 0 ln Example 21 use Maple to nd the numerical values of the rst three eigenvalues and plot the graphs of the corresponding eigenfunctions a Show that the solution of ut Du tgt 0 m E 0L umt 0 umt L 0 2100 M0790 uo7 z 6 07L7 t cos L where cm u0mcos dm m 2 0 tgt 0 z E 0 L 0 co co mzn392 utm 3 Z cmezp lt7D L2 m1 2101 00 to H O H H H to CHAPTER 2 DIFFUSION EQUATIONS a Find the solution of ut 472m tgt 0 m 6 01 um 0 Mt 1 0 2102 1407 90 Hm 1 m e 0 05 71 m e 051 b Modify the Maple program in Section 24 to simulate the solution of 2102 where fz is given by the formula fz Find the solutions of equilibrium equation u 0 m 6 01 u0 310 5 u1 7 511 7 2103 Find the solutions of equilibrium equation 0 31 0 m 6 01 00 5 71 4 2104 Suppose that ut z is the solution of 2 7DAu tgt0 69 u07 22 mm m e 0 2105 Vutm nz a utm 0 tgt 0 m E 89 where a gt O Show that the mean squared value 25 ut z2dm is strictly decreasing 9 unless ut z E 0 Suppose that ut z is the solution of the convective diffusion equation D gt 0 and V gt 0 814 8214 814 a 7 DW 7 Va m E 07r u0x u0x m 6 07139 2106 814 814 t0 a r 0 1 Show that the total population ut mdz is a constant if u0z sinz but not a constant if u0z cosz Suppose that in the chemical mixing problem in Section 25 the salt solution drifts to the right with velocity V Then 005 z satis es the convectivediffusion equation 36 82 36 D 7 V 2107 at 8x2 8x7 with DV gt 0 Suppose that ctm satis es the same initial and boundary conditions in 266 Find the equilibrium solutions in this situation and explain the effect of the drifting to the related pollution problem 27 SEPARATION OF VARIABLES RECTANGLES 13 H F H U H a 35 Suppose that in the chemical mixing problem in Section 25 the salt is consumed in a chemical reaction at a rate k Then 005 z satis es a linear diffusion equation 31 D 8t 8213 amz 7 kc 2108 with Dk gt 0 Suppose that ctm satis es the same initial and boundary conditions in 266 Find the equilibrium solutions in this situation and discuss the possible spatial patterns of equilibrium solution Determine the critical patch for a diffusive Malthus equation with Robin boundary conditions 8P NP 5 DWaP tgt0 m 0L 7 x 7 7 2109 umt0 7 but0 umt L 7 butL u0x u0x m E OL Consider the doubly periodic boundary value problem 814 8214 8214 a D 4147 tgt07 79ER07 1X0757 140711 Why 712079 ma y 2110 ux 0 ux b uyz0 uyxb u0my u0xy xy E R Find the eigenvalues and eigenfunctions of the problem and the series representation of the solution Consider a convectivediffusion equation 8214 814 8t 4E8utgt0m ut0 0 ut L 0 u0x u0x m E OL 07L7 2111 Find the solution of the equation with the following steps a Use separation of variables method to show that if ut z UtVz is a solution then for some constant k U and V satisfy U t kUt V z 4Vz 8Vz kVz V0 VL 0 b Find the eigenvalues and eigenfunctions of V z 4Vm 8Vz kVx V0 VL 0 Hint treat the cases of k gt 4 k 4 and k lt 4 separately