Solve Predator-Prey Equations

Open Live Script

This example shows how to solve a differential equation representing a predator/prey model using variable step size Runge-Kutta integration methods. The ode23 method uses a 2nd and 3rd order pair of formulas for medium accuracy, and the ode45 method uses a 4th and 5th order pair for higher accuracy.

Consider the pair of first-order ordinary differential equations known as the Lotka-Volterra equations, or predator-prey model:

dxdt=x-αxydydt=-y+βxy.

The variables x and y measure the sizes of the prey and predator populations, respectively. The quadratic cross term accounts for the interactions between the species. The prey population increases when no predators are present, and the predator population decreases when prey are scarce.

Code Equations

To simulate the system, create a function that returns a column vector of state derivatives, given state and time values. The two variables x and y can be represented in MATLAB® as the first two values in a vector y. Similarly, the derivatives are the first two values in a vector yp. The function must accept values for t and y and return the values produced by the equations in yp.

yp(1) = (1 - alpha*y(2))*y(1)

yp(2) = (-1 + beta*y(1))*y(2)

In this example, the equations are contained in a file called lotka.m. This file uses parameter values of α=0.01 and β=0.02.

type lotka
function yp = lotka(t,y)
%LOTKA  Lotka-Volterra predator-prey model.

%   Copyright 1984-2014 The MathWorks, Inc.

yp = diag([1 - .01*y(2), -1 + .02*y(1)])*y;

Simulate System

Create an ode object to define the problem and initial conditions. Use an initial condition of x(0)=y(0)=20 so that the populations of predators and prey are initially equal. Specify that the ode23 solver should be used. Then use the solve method to simulate the system over the time interval 0<t<15.

y0 = [20; 20];   
F = ode(ODEFcn=@lotka,InitialValue=y0,Solver="ode23")
F = 
  ode with properties:

   Problem definition
               ODEFcn: @lotka
          InitialTime: 0
         InitialValue: [2×1 double]
         EquationType: standard

   Solver properties
    AbsoluteTolerance: 1.0000e-06
    RelativeTolerance: 1.0000e-03
               Solver: ode23

  Show all properties



t0 = 0;
tfinal = 15;
S = solve(F,t0,tfinal);
t = S.Time;
y = S.Solution;

Plot Results

Plot the resulting populations against time.

plot(t,y,"-o")
title("Predator/Prey Populations Over Time")
xlabel("t")
ylabel("Population")
legend("Prey","Predators",Location="north")

Figure contains an axes object. The axes object with title Predator/Prey Populations Over Time, xlabel t, ylabel Population contains 2 objects of type line. These objects represent Prey, Predators.

Now plot the populations against each other. The resulting phase plane plot makes the cyclic relationship between the populations very clear.

plot(y(1,:),y(2,:))
title("Phase Plane Plot")
xlabel("Prey Population")
ylabel("Predator Population")

Figure contains an axes object. The axes object with title Phase Plane Plot, xlabel Prey Population, ylabel Predator Population contains an object of type line.

Compare Results of Different Solvers

Solve the system a second time using the ode45 solver, instead of ode23. The ode45 solver takes longer for each step, but it also takes larger steps. Nevertheless, the output of ode45 is smooth because by default the solver uses a continuous extension formula to produce output at four equally spaced time points in the span of each step taken. (You can adjust the number of points with the Refine name-value argument of the solve method.) Plot both solutions for comparison.

F.Solver = "ode45";
S2 = solve(F,t0,tfinal);
t2 = S2.Time;
y2 = S2.Solution;

plot(y(1,:),y(2,:),"o",y2(1,:),y2(2,:),".");
title("Phase Plane Plot")
legend("ode23","ode45")

Figure contains an axes object. The axes object with title Phase Plane Plot contains 2 objects of type line. One or more of the lines displays its values using only markers These objects represent ode23, ode45.

The results show that solving differential equations using different numerical methods can produce slightly different answers.

See Also

|

Topics