Math 4452 (Peh Ng)

Mathematica Examples for the Predator-Prey Model :
      dH/dt=(α-βP)H
               dP/dt=(-μ+ϕH)P
where H(t) and P(t) represent the populations of 
the host-prey and parasite-predator, respectively.

Mathematical Analysis

First, input the rate equations:

Clear[α, β, μ, ϕ, H, P]
f[H_, P_] = (α- β*P)*H;  (* f is dH/dt *)
g[H_, P_] = (-μ + ϕ*H)*P;      (*  g is dP/dt  *)
Factor[f[H, P]]
Factor[g[H,P]]

Second, we can use mathematica to find the equilibrium values:

solution = Solve[{f[H,P] == 0, g[H,P] == 0}, {H,P}]

Now, we can  input the  various parametres in the above rate equations:

α = .4;   (* Intrinsic Growth Rate of Prey *)
β = .01;  (* Coefficient of Interaction of Prey *)
μ = .5;   (* Intrinsic Growth Rate of Predator *)
ϕ = .02;  (* Coefficient of Interaction of Predator *)

equilibrium = {H,P} /. solution

In class, we did not solve for H and P explicitly as a function of time, but we found an equation , that relates H and P, i.e. F(H,P)=C where C is a constant.  From Calc3, we know that plotting vector fields given by the vector <dH/dt, dP/dt> will enable us to see contours formed by the curve F(H,P)=C. So, we can use PlotVectorField to do this. (The values of maxpredator and maxprey are for the range.)

Needs["Graphics`PlotField`"]
Clear[H,P]
maxprey = 50;
maxpredator = 80;
directionField = PlotVectorField[{f[H,P], g[H,P]},
				{H,0,maxprey}, {P,0,maxpredator},
				ScaleFunction -> (1&),
				ScaleFactor -> 3,
				PlotPoints -> 15,
				Axes -> True,
				AxesLabel -> {"PREY", "PREDATORS"}]

Trajectories in time

initialprey = 50;
initialpredator = 10;
α = .4;   (* Intrinsic Growth Rate of Prey *)
β = .01;  (* Coefficient of Interaction of Prey *)
μ = .5;   (* Intrinsic Growth Rate of Predator *)
ϕ = .02;  (* Coefficient of Interaction of Predator *)
maxtime = 20;
maxprey = 100;
maxpredator = 150;

Trajectory, via numerical solutions in time,  on the (H,P)-plane

This is how you numerically solve a system of differential equations wrt time, t; and parametrically plot the graph of Predator vs Prey implicitly as a function of time.

solution = NDSolve[{H'[t] == H[t]*(α - β*P[t]),
					P'[t] == P[t]*(-μ + ϕ*H[t]),
					H[0] == initialprey,
					P[0] == initialpredator},
					{H,P}, {t,0,maxtime},
					MaxSteps -> 2000];
graphTraj = ParametricPlot[Evaluate[{H[t],P[t]} /. solution],
				{t,0,maxtime},
				PlotRange -> {{0,maxprey},{0,maxpredator}},
				AxesLabel -> {PREY,PREDATOR},
				PlotStyle -> {RGBColor[1,0,0]}]

Trajectories or level curves, via the algebraically obtained function of H and P, F(H,P)=Constant,  on the (H,P)-plane

ContourPlot[(P^(α))*(H^(μ))/(Exp[ϕ*H+β*P]), {H,0.02,100}, {P,0.02,150}, Axes->True,AxesLabel->{"H-host-prey","P-parasite-predator"}]

Trajectory showing H (host-prey)  as function of t

graphHost = ParametricPlot[Evaluate[{t,H[t]} /. solution],
				{t,0,maxtime},
				PlotRange -> {{0,maxtime},{0,maxprey}},
				AxesLabel -> {TIME,PREY},
				PlotStyle -> {RGBColor[0,1,0]}]

Trajectory showing P (parasite-predator) as function of t

graphParasite = ParametricPlot[Evaluate[{t,P[t]} /. solution],
				{t,0,maxtime},
				PlotRange -> {{0,maxtime},{0,maxpredator}},
				AxesLabel -> {TIME,PREDATOR},
				PlotStyle -> {RGBColor[0,0,1]}]

Graphs shown together

Show[GraphicsArray[{graphTraj,graphHost,graphParasite}]]

Several trajectories together in the (H,P)-plane

initialprey = 50;
initialpredator = 40;
solution = NDSolve[{H'[t] == H[t]*(α - β*P[t]),
					P'[t] == P[t]*(-μ + ϕ*H[t]),
					H[0] == initialprey,
					P[0] == initialpredator},
					{H,P}, {t,0,maxtime},
					MaxSteps -> 2000];
graphtraj1 = ParametricPlot[Evaluate[{H[t],P[t]} /. solution],
				{t,0,maxtime},
				PlotRange -> {{0,maxprey},{0,maxpredator}},
				AxesLabel -> {PREY,PREDATOR},
				PlotStyle -> {RGBColor[0,1,0]},
				DisplayFunction -> Identity]
initialprey = 30;
initialpredator = 40;
solution = NDSolve[{H'[t] == H[t]*(α-β*P[t]),
					P'[t] == P[t]*(-μ+ϕ*H[t]),
					H[0] == initialprey,
					P[0] == initialpredator},
					{H,P}, {t,0,maxtime},
					MaxSteps -> 2000];
graphtraj2 = ParametricPlot[Evaluate[{H[t],P[t]} /. solution],
				{t,0,maxtime},
				PlotRange -> {{0,maxprey},{0,maxpredator}},
				AxesLabel -> {PREY,PREDATOR},
				PlotStyle -> {RGBColor[1,0,1]},
				DisplayFunction -> Identity]
Show[graphTraj, graphtraj1, graphtraj2, DisplayFunction -> $DisplayFunction]