The Lotka-Volterra predator-prey equation is defined below.

In[187]:=

nt = n '[t]  rn n[t] - a p[t] n[t] ;

pt = p '[t]  a rp n[t] p[t] - b p[t] ;

system = {nt, pt}

Out[189]=

{n^′[t] rn n[t] - a n[t] p[t], p^′[t]  -b p[t] + a rp n[t] p[t]}

where
    n[t] is the density of the prey species,
    p[t] is the density of predator species,
    rn is the per capita growth rate of the prey,
    rp is  the per capita growth rate of the predator,
    a is the per capita attack rate of the predator,
    b is the per capita death rate of the predator without the prey.

The  predator-prey system is parameterized using the /. replacement operator.

In[190]:=

predsystem = system /. {rn ->0.5, rp ->0.2, a->0.1 , b->0.8}

Out[190]=

{n^′[t] 0.5 n[t] - 0.1 n[t] p[t], p^′[t]  -0.8 p[t] + 0.02 n[t] p[t]}

Initial conditions need to be added to the predator-prey system. The initial conditions are the densities of the prey (n[0]) and the predators (p[0]).

In[191]:=

fullsystem = Join[predsystem, {n[0] 15, p[0]  50}]

Out[191]=

{n^′[t] 0.5 n[t] - 0.1 n[t] p[t], p^′[t]  -0.8 p[t] + 0.02 n[t] p[t], n[0] 15, p[0] 50}

The NDSolve function is the heavy lifter in Mathematica for numerically solving differential equations. For NDSolve to work the variables in the differential equation need to defined ({n,p}) and the range of values for which you want to numerically solve the equations needs to be given ({t,0,100}).

In[192]:=

predprey = NDSolve[fullsystem, {n, p}, {t, 0, 50}]

Out[192]=

{{nInterpolatingFunction[{{0., 50.}}, <>], pInterpolatingFunction[{{0., 50.}}, <>]}}

NDSolve returns an InterpolatingFunction which is the numerical solution to the two inerlinked differential equations. The next three steps isolates out the functions for the predators and prey.

In[193]:=

curves = Flatten[predprey]

Out[193]=

{nInterpolatingFunction[{{0., 50.}}, <>], pInterpolatingFunction[{{0., 50.}}, <>]}

In[194]:=

prey = n[t] /. curves[[1]]

Out[194]=

InterpolatingFunction[{{0., 50.}}, <>][t]

In[195]:=

predator = p[t] /. curves[[2]] 

Out[195]=

InterpolatingFunction[{{0., 50.}}, <>][t]

Below is a plot of the number of prey (in red) and predators (in blue) versus time.

In[196]:=

Plot[{prey, predator}, {t, 0, 50}, PlotRangeAll, PlotStyle {Red, Blue}]

[Graphics:HTMLFiles/index_18.gif]

Out[196]=

⁃Graphics⁃


The ParametricPlot function allows you to plot the limit cycles for the density of predator and prey. Notice how close the density of the predators and the prey get to zero.

In[197]:=

ParametricPlot[{prey, predator}, {t, 0, 50}, PlotRangeAll]

[Graphics:HTMLFiles/index_21.gif]

Out[197]=

⁃Graphics⁃

If you are interested in the numeric values at each time step you can use the Table function. This table is generated at a step of t=1.

In[198]:=

Table[prey, {t, 0, 50, 1}]

Out[198]=

{15., 0.658334, 0.203365, 0.157416, 0.184572, 0.260947, 0.401371, 0.641264, 1.04224, 1.70713 ... 29.1825, 0.88456, 0.223554, 0.157959, 0.177746, 0.246682, 0.376276, 0.598911, 0.971755, 1.59045}

In[200]:=

Table[predator, {t, 0, 50, 1}]

Out[200]=

{50., 24.2394, 10.9675, 4.94486, 2.2293, 1.00607, 0.455002, 0.206546, 0.0943506, 0.043552, 0 ... 7, 52.7555, 27.1179, 12.2892, 5.54167, 2.49819, 1.12719, 0.509581, 0.231169, 0.105482, 0.048601}


Created by Mathematica  (May 25, 2005) Valid XHTML 1.1!