ME 406Predator-Prey ExamplePart 2 - Limit CycleIn[327]:= sysidMathematica 4.1.2, DynPac 10.66, 3ê8ê2002In[328]:=intreset;In[329]:= plotreset;ü IntroductionThis is a simple model of an ecological system with three components: a plant, a small mammal called the Murat which eats the plant, and a carnivorous predator called the Vekton which eats the mammal. We are using the model to test our belief that we can control this ecosystem by controlling the plant population. Our goal is to maintain both species at a healthy numbers. We denote the population of the plant eater by M and the population of the carnivore by V. The basic equations governing this system were set-up in class. They are dMÅÅÅÅÅÅÅÅÅÅÅÅdt= rM J1 -MÅÅÅÅÅÅÅÅAN -bMVÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅM + H,dVÅÅÅÅÅÅÅÅÅÅÅdt=bMVÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅM + H- cV .The parameters were explained in class, so we just review them very briefly here. The parameter A is the maximum sustainable population of the plant-eaters in the absence of the carnivores, and the model incorporates this via a logistic law. The parameter A is the one we control, by controlling the plant population, and this is our only management tool for this system. The birth rate of the plant-eaters is r. The death rate is in effect incorporated in the logistic law. The second, negative, term in the M-equation models the effects of the predators. The loss is proportional to the product of the populations, because this product is proportional to the encounter rate. The denominator is a saturation effect, accounting for the fact that if the population of plant eaters is large, predators won't be as hungry and there will be fewer kills per encounter. In the predator equation, the first term models the birth rate, and accounts for the fact that it will be higher if food is more plentiful. If food is very plentiful (M much larger than H), this birth term saturates with a rate b*V, so b is the natural birth rate with ample food. It is an artificiality of the model that the saturation parameter has the same value H in both equations. ü Defining the System for MathematicaWe define the system for DynPac.In[330]:= setstate@8M, V<D;predpt2.nb 1In[331]:= setparm@8A<D;In[332]:= slopevec = 8r * M * H1 - M ê AL - Hb*M * VLêHM + HL, Hb * M * VLêHM + HL - c * V<; With so many parameters, the exploration of the system could take a very long time. To keep things simple, we choose values for all of the parameters except A (our management variable), and then study the system as we vary A. This is why we have included only A in the parameter vector. The units for the population quantities V, M, H, and A are in millions of individuals. The units of the time constants r, beta, b and c are in inverse years. We will assume that A can be varied from 0 up to a maximum of 15. We specify now all of the parameters that will be fixed throughout this study.In[333]:= r = 12 H** yr-1**L;In[334]:= b=20 H** yr-1**L;In[335]:= b = 4H** yr-1**L;In[336]:= c = 8 ê 5 H** yr-1**L;In[337]:= H = 4 H** millions **L;In our evaluations below of the effects of changing A, we adopt the following criterion for a healthy ecosystem: if the population of each species remains above 1 million, we count this as a successful preservation strategy. ü Summary of Part 1In Part 1, we did a detailed of the equilibrium states as we varied the parameter A. Here is a summary of the results. (a) For 0 < A < 8/3, the only stable equilibrium state is the all-Murat state at the maximum sustainable population A.(b) For 8/3 < A < 4ÅÅÅ3I-5 +è!!!!!!74 M= 4.8031, the only stable equilibrium state is the state coeq in which members of both species are present. This equilibrium is a stable node.(c) For 4ÅÅÅ3I-5 +è!!!!!!74 M< A < 28/3, the only stable equilibrium state is coeq, and the state is a stable spiral.(d) For 28/3 < A, there are no stable equilibrium states.Because the system is bounded, these results sharply pose the question of what happens to the system for A > 28/3. The results of Part 1 for A below but near 28/3 suggest the answer. The stable equilibrium, which is a stable spiral, becomes highly oscillatory. If we go just beyond 28/3, we would expect the spiral to be unstable, and we would expect to find a limit cycle around it. In other words, we would look for a Hopf bifurcation. Let us verify all of this with our computer work.ü The Equilibrium StatesWe first find the equilibrium states in terms of the parameter A. In[338]:= eqstates = findpolyeqOut[338]= 980, 0<, 8A, 0<, 98ÅÅÅÅ3,4 H-8 + 3ALÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ3A==predpt2.nb 2We see that there are three equilibrium states. The first is {0,0}, in which no animals of either species are present. The second is {A,0}, in which the Vektons are absent, and the Murats are living free of fear at their maximum sustainable population A. In the third state, both species are present. As we saw in Part 1, it is the third state which is the stable equilibrium for A under 28/3. We name this state coeq.In[339]:= coeq = eqstates@@3DDOut[339]= 98ÅÅÅÅ3,4 H-8 + 3ALÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅÅ3A=Let's check the stability of this state just below and just above 28/3.In[340]:= parmval = 27 ê 3;In[341]:= classify2D@coeqDAbbreviations used in classify2D.L = linear, NL = nonlinear, R2 = repeated root.Z1 = one zero root, Z2 = two zero roots.This message printed once.strictly stable - spiralIn[342]:=parmval = 29 ê 3;In[343]:= classify2D@coeqDunstable - spiralü First Look At The Limit CycleWe set A to a value above the bifurcation value, and look at some typical orbits.In[344]:= parmval = 810<;In[345]:= nsteps = 1000;In[346]:= h = 0.02;In[347]:= t0 = 0.0;In[348]:= initset = 881, 1<, 84, 3<, 86, 6<, 810, 10<<;In[349]:= plrange = 88-1, 12<, 8-1, 12<<;In[350]:= asprat = 1.0;predpt2.nb 3In[351]:= graph1 = portrait@initset, t0, h, nsteps, 1, 2D;246810 12M24681012VSystem 8A<=8 10.00<The limit cycle is a little fuzzy because of the approach orbits. Clearly it is an attractor for the first quadrant. The orbit leaving the saddle point, which in Part 1 connected with the equilibrium coeq, now merges into the limit cycle. Let's look at the pure limit cycle. In[352]:= nsteps = 2000;In[353]:= sol2 = limcyc@85, 3<, t0, h, nstepsD;predpt2.nb 4In[354]:= graph2 = phaser@sol2D;246810 12M24681012VSystem 8A<=8 10.00<We get the