Models for Population Dynamics MAPLE worksheet by J. Reinhardt, 2009 The development of populations of animals (or plants) competing for resources or preying upon each other can be modelled using using systems of coupled differential equations.This assumes that the number of animals is large and can be treated as a continuous variable. Also development in time is assumed to be continuous. Literature: J.D. Murray, Mathematical Biology, Springer Verlag

1.1 Tangent field plot restart: with(plots): with(DEtools): ps1:=diff(x(t),t)=x(t)*(1-y(t)); ps2:=diff(y(t),t)=alpha*y(t)*(x(t)-1); alpha:=1: # system parameter dfieldplot([ps1,ps2],[x,y],t=0..1,x=0..+5,y=0..+5,dirgrid=[20,20],color=blue,axes=boxed,title=`Lotka-Volterra system: Tangent field plot`); 1.2 Integration of the differential equations restart: with(plots): ps1:=diff(x(t),t)=x(t)*(1-y(t)): ps2:=diff(y(t),t)=alpha*y(t)*(x(t)-1): alpha:=1: # system parameter x0:=0.5: y0:=0.5: # initial conditions sol:=dsolve({ps1,ps2,x(0)=x0,y(0)=y0},{x(t),y(t)},type=numeric,output=listprocedure): plot1:=odeplot(sol,[t,x(t)],0..30,numpoints=500,color=blue): plot2:=odeplot(sol,[t,y(t)],0..30,numpoints=500,color=red): display({plot1,plot2},view=[0..30,0..3],axes=boxed,labels=[t,x], title= `Lotka-Volterra System: Trajectory`); The two populations oscillate in synchrony but with a phase shift: First the rabbits multiply, then the population of foxes grows, who eat most of the rabbits and then die from starvation, and so on. 1.3a Phase space plot, by integration of the equation of motion restart: with(plots): ps1:=diff(x(t),t)=x(t)*(1-y(t)): ps2:=diff(y(t),t)=alpha*y(t)*(x(t)-1): alpha:=1: # system parameter x0:=array(1..5): # initial condition, setting x(0)=y(0)=x0[i] x0[1]:=1.1: x0[2]:=2.: x0[3]:=3.: x0[4]:=4.: x0[5]:=5.: for i from 1 to 5 do sol:=dsolve({ps1,ps2,x(0)=x0[i],y(0)=x0[i]},{x(t),y(t)},type=numeric,output=listprocedure): outplot(i):=odeplot(sol,[x(t),y(t)],0..100,numpoints=10000,color=blue): od: display([seq(outplot(i),i=1..5)],view=[0..+10,0..+10],axes=boxed,labels=[x,y], title= `Lotka-Volterra System: Phase space trajectories`); 1.3b Phase space plot, using contour lines for the conserved quantity C Both methods lead to the same results. restart: with(plots): alpha:=1: f(x,y):= alpha*x+y-log(x^alpha*y): contourplot(f(x,y),x=0.01..10,y=0.01..10,axes=boxed,grid=[50,50],color=blue);

2.1 Tangent field plot restart: with(plots): with(DEtools): ps1:=diff(x(t),t)=x(t)*(1-x(t)-a*y(t)); ps2:=diff(y(t),t)=rho*y(t)*(1-y(t)-b*x(t)); Example case a>1 and b>1: a:=4/3: b:=3/2: rho:=2/3: # system parameters dfieldplot([ps1,ps2],[x,y],t=0..1,x=0..+1.5,y=0..+1.5,dirgrid=[30,30],color=blue,axes=boxed,title=`Competition model: Tangent field plot`); 2.2 Determination of the fixed points and stability analysis restart: f:= (x,y) -> x*(1-x-a*y); g:= (x,y) -> rho*y*(1-y-b*x); with(VectorCalculus): with(LinearAlgebra): Jac:=Jacobian([f(x,y),g(x,y)],[x,y]); fixpoints:=solve({f(x,y)=0,g(x,y)=0},{x,y}); There are 4 fixed points. Let us evaluate the eigenvalues of the stability matrix: for i from 1 to 4 do # loop over the four solutions unassign('x','y'); # Free the fixpoint variables print(`Fixed point number`, i); assign(fixpoints[i]); # Assign output no. i of 'solve' to the variable x,y x,y; simplify(Jac); lambda:=simplify(Eigenvalues(Jac)); od; Fixed point LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUklbXN1YkdGJDYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYnLUkjbW5HRiQ2JVEiMUYnRjJGNUYyLyUscGxhY2Vob2xkZXJHRjQvJTBmb250X3N0eWxlX25hbWVHUSVUZXh0RidGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnL0krbXNlbWFudGljc0dGJFEnYXRvbWljRictSSNtb0dGJDYuUSJ+RidGMkY1LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZPLyUpc3RyZXRjaHlHRk8vJSpzeW1tZXRyaWNHRk8vJShsYXJnZW9wR0ZPLyUubW92YWJsZWxpbWl0c0dGTy8lJ2FjY2VudEdGTy8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRmhuRjJGQEY1is always instable. Fixed point LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUklbXN1YkdGJDYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYnLUkjbW5HRiQ2JVEiMkYnRjJGNUYyLyUscGxhY2Vob2xkZXJHRjQvJTBmb250X3N0eWxlX25hbWVHUSVUZXh0RidGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnL0krbXNlbWFudGljc0dGJFEnYXRvbWljRictSSNtb0dGJDYuUSJ+RidGMkY1LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZPLyUpc3RyZXRjaHlHRk8vJSpzeW1tZXRyaWNHRk8vJShsYXJnZW9wR0ZPLyUubW92YWJsZWxpbWl0c0dGTy8lJ2FjY2VudEdGTy8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRmhuRjJGQEY1is stable if a>1, otherwise hyperbolic. Fixed point LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUklbXN1YkdGJDYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYnLUkjbW5HRiQ2JVEiM0YnRjJGNUYyLyUscGxhY2Vob2xkZXJHRjQvJTBmb250X3N0eWxlX25hbWVHUSVUZXh0RidGNS8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnL0krbXNlbWFudGljc0dGJFEnYXRvbWljRictSSNtb0dGJDYuUSJ+RidGMkY1LyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZPLyUpc3RyZXRjaHlHRk8vJSpzeW1tZXRyaWNHRk8vJShsYXJnZW9wR0ZPLyUubW92YWJsZWxpbWl0c0dGTy8lJ2FjY2VudEdGTy8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRmhuRjJGQEY1is stable if b>1, otherwise hyperbolic. Fixed point LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUklbXN1YkdGJDYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYmLUkjbW5HRiQ2JVEiNEYnRjJGNUYyLyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRjUvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy9JK21zZW1hbnRpY3NHRiRRJ2F0b21pY0YnLUYvNiNRIUYnRjJGPkY1 is more difficult to analyze. (a) For either (a<1 and b>1) or (a>1 and b<1) LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzZTLUkjbWlHRiQ2JVEkdGhlRicvJSdpdGFsaWNHUSZmYWxzZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1JI21vR0YkNi1RIn5GJ0YyLyUmZmVuY2VHRjEvJSpzZXBhcmF0b3JHRjEvJSlzdHJldGNoeUdGMS8lKnN5bW1ldHJpY0dGMS8lKGxhcmdlb3BHRjEvJS5tb3ZhYmxlbGltaXRzR0YxLyUnYWNjZW50R0YxLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGSS1GLDYlUSZmaXhlZEYnRi9GMkY1LUYsNiVRJnBvaW50RidGL0YyRjUtSSVtc3ViR0YkNiYtRiw2JVEieEYnL0YwUSV0cnVlRicvRjNRJ2l0YWxpY0YnLUYjNiQtSSNtbkdGJDYkUSI0RidGMkYyLyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvSSttc2VtYW50aWNzR0YkUSdhdG9taWNGJ0Y1LUYsNiVRJWRvZXNGJ0ZYRlotRjY2LkY4RlhGWkY5RjtGPUY/RkFGQ0ZFRkdGSi1GNjYvUSRub3RGJy8lJWJvbGRHRjFGWEZaRjlGO0Y9Rj9GQUZDRkVGR0ZKRmVvLUYsNiVRJmV4aXN0RidGWEZaRmVvLUY2Ni5RI2luRidGam9GMkY5RjtGPUY/RkFGQ0ZFRkdGSkY1RitGNS1GLDYlUSthZG1pc3NpYmxlRidGL0YyRjUtRiw2JVEncmVnaW9uRidGL0YyRjUtRiw2JVEjb2ZGJ0YvRjJGNS1GLDYlUStjb29yZGluYXRlRidGL0YyRjUtRiw2JVEmc3BhY2VGJ0YvRjJGNS1JJ21zcGFjZUdGJDYmLyUnaGVpZ2h0R1EmMC4wZXhGJy8lJndpZHRoR1EmMC4wZW1GJy8lJmRlcHRoR0ZmcS8lKmxpbmVicmVha0dRKG5ld2xpbmVGJy1GYnE2JkZkcUZncUZqcS9GXXJRJWF1dG9GJy1GLDYlUShiZWNhdXNlRidGL0YyRjUtRiw2JVEkb25lRidGL0YyRjVGaHBGNUYrRjUtRiw2JVEsY29vcmRpbmF0ZXNGJ0YvRjJGNS1GLDYlRldGL0YyLUY2Ni1RIixGJ0YyRjkvRjxGWUY9Rj9GQUZDRkVGRy9GS1EsMC4zMzMzMzMzZW1GJy1GLDYlUSJ5RidGL0YyRjUtRiw2JVEjaXNGJ0YvRjJGZW8tRiw2JVEpbmVnYXRpdmVGJ0ZYRlotRjY2LVEiOkYnRjJGOUY7Rj1GP0ZBRkNGRS9GSFEsMC4yNzc3Nzc4ZW1GJy9GS0ZhdEYy assume(a>1,b<1); signum(x*y); For other values we evaluate the stability coefficients using specific sample values for the parameters. (b) LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1YkdGJDYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYmLUkjbW5HRiQ2JFEiNEYnL0Y2USdub3JtYWxGJy8lLHBsYWNlaG9sZGVyR0Y0LyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRj4vJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy9JK21zZW1hbnRpY3NHRiRRJ2F0b21pY0YnLUkjbW9HRiQ2LVEifkYnRj4vJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRlEvJSlzdHJldGNoeUdGUS8lKnN5bW1ldHJpY0dGUS8lKGxhcmdlb3BHRlEvJS5tb3ZhYmxlbGltaXRzR0ZRLyUnYWNjZW50R0ZRLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGam5GQkY+is a hyperbolic point for a>1 and b>1. Interpretation: Depending on the initial condition one of the two species wins. a:=4/3: b:=3/2: rho:=2/3: # system parameters evalf(x); evalf(y); evalf(lambda); LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1YkdGJDYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYmLUkjbW5HRiQ2JFEiNEYnL0Y2USdub3JtYWxGJy8lLHBsYWNlaG9sZGVyR0Y0LyUwZm9udF9zdHlsZV9uYW1lR1ElVGV4dEYnRj4vJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy9JK21zZW1hbnRpY3NHRiRRJ2F0b21pY0YnLUkjbW9HRiQ2LVEifkYnRj4vJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRlEvJSlzdHJldGNoeUdGUS8lKnN5bW1ldHJpY0dGUS8lKGxhcmdlb3BHRlEvJS5tb3ZhYmxlbGltaXRzR0ZRLyUnYWNjZW50R0ZRLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGam5GQkY+is a stable node for a<1 and b<1: Coexistence of both species! a:=2/3: b:=1/2: rho:=2/3: # system parameters evalf(x); evalf(y); evalf(lambda); JSFH