http://www.radford.edu/~thompson/webddes/tutorial.pdf (Problem 1 on page 32).
A model for predator-prey populations is given by:
where and are the prey and predators. These are ordinary differential equations that are straightforward to solve.
If there is a resource limitation on the prey and assuming the birth rate of predators responds to changes in the magnitude of the population y1 of prey and the population y2 of predators only after a time delay , we can arrive at a new set of delay differential equations:
these equations cannot be solved with ordinary ODE solvers, because the value of the solution at any time depends on the value of the solution some time ago. So, even at , we have to know the history to know the initial derivative is. We use the dde23 function in Matlab to solve this problem.
clear all; close all; clc
We first solve the ordinary model so we can see the impact of adding the delay.
a =0.25; b= -0.01; c= -1.00; d=0.01; y10 = 80; % initial population of prey y20 = 30; % initial population of predator tspan = [0 200]; % time span to integrate over [t,Y] = ode45(@myode,tspan,[y10 y20]); y1 = Y(:,1); y2 = Y(:,2);
You can see a stable limit cycle (purely oscillatory behavior) in the population of prey and predators.
plot(y1,y2) xlabel('y_1') ylabel('y_2')
function dYdt = myode(t,Y) % ordinary model y1 = Y(1); y2 = Y(2); dy1dt = a*y1 + b*y1*y2; dy2dt = c*y2 + d*y1*y2; dYdt = [dy1dt; dy2dt]; end
Now we solve the model with the delays. We solve this example with a delay of .
lags = ; % this is where we specify the vector of taus % we need a history function that is a column vector giving the value of y1 % and y2 in the past. Here we make these constant values. function y = history(t) y = [80; 30]; end % we define the function for the delay. the Y variable is the same as you % should be used to from an ordinary differential equation. Z is the values % of all the delayed variables. function dYdt = ddefun(t,Y,Z) m = 200; % additional variable y1 = Y(1); y2 = Y(2); % Z(:,1) = [y1(t - tau_1); y2(t - tau_2)] y1_tau1 = Z(1,1); y2_tau1 = Z(2,1); dy1dt = a*y1*(1-y1/m) + b*y1*y2; dy2dt = c*y2 + d*y1_tau1*y2_tau1; dYdt = [dy1dt; dy2dt]; end
now we solve with dde23
sol = dde23(@ddefun,lags,@history,tspan) y1 = sol.y(1,:); % note the y-solution is a row-wise matrix. y2 = sol.y(2,:);
sol = solver: 'dde23' history: @dde_example1/history discont: [0 1 2 3] x: [1x127 double] y: [2x127 double] stats: [1x1 struct] yp: [2x127 double]
hold all plot(y1,y2)
you can see the delay and resource limit significantly affect the solution. Here, it damps the oscillatory behavior, and seems to be leading to steady state population values.
figure plot(sol.x,y1,sol.x,y2) xlabel('Time') legend 'y1' 'y2'
end % categories: ODEs % tags: delay differential equation