% Time Delay Differential Equations
% Consider the following model for the growth of an animal population.
% \frac{\textrm{dN}\left(t\right)}{\textrm{dt}}=\textrm{rN}\left(t\right)\left(1-\frac{N\left(t\right)}{K}\right)
% Without the parenthetical term at the end of this equation, a model such as this will predict eponential growth of the population. The growth constant, r, determines the rate of this growth. The parenthetical term introduces the idea of resource limitations, indicating that the ultimate size of the population will be limited by the available resources in the habitat of these animals. A sample solution would be as follows

trange=[0 20];
init=1;
[t,y]=ode45(@rkfile, trange, init);
yexact=exp(t.^2);
plot(t,y, 'LineWidth', 2)
xlabel('t')
ylabel('N')
set(gca,'FontSize',20)

% You can see in the figure that N(t) stops increasing as it approaches a population of 100 and it takes on the order of 10 time units to get reasonably close to that value.
% Now, suppose that we introduce a delay into the system. Specifically, suppose that the growth rate at a certain time depends on what the population was in the past. In other words, what if the growth rate today depends on how close we were to a population of 100 30 days ago. In this case, our equation becomes:
% \frac{\textrm{dN}\left(t\right)}{\textrm{dt}}=\textrm{rN}\left(t\right)\left(1-\frac{N\left(t-\tau \right)}{K}\right)
% Here, \tau represents the delay. We now must account for N at times less than t=0 in order to solve the equation.
% To solve this equation in Matlab, we can use the dde23 function. Let's assume r=1, K=100, \tau =30, and N=50 for times less than 0.

tau=1;
trange2=[0 20];
sol = dde23(@ddex1de,tau,@ddex1hist,trange2);
plot(sol.x,sol.y)
xlabel('t');
ylabel('N');
set(gca,'FontSize',20)

%

function s = ddex1hist(t)
% Constant history function (t<0)
s = 50;
end
% ------------------------------------------------
function dydt = ddex1de(t,y,Z)
r=1;
K=100;
Nlag=Z(:,1);
N=y(1);
dydt = r*N*(1-Nlag/K);
end

% Notice that the lag introduces oscillations. In this case, the oscillation is damped. If we change r to a value of 2, then the oscillations will persist.

function dydt=rkfile(t,y)
  r=1;
  K=100;
  N=y(1);
  dydt=r*N*(1-N/K);
end

% © James Blanchard, University of Wisconsin-Madison, blanchard@engr.wisc.edu