Project A: Corn, cats and rats 3
Project A: Corn, cats and rats
On an isolated island, some begins to grow. From the beginning there are 75 plants but soon
they increase in number according to dP/dt=α1P−α2P ^2.
The unit of time is years, P (t) the number of plants at time t and the constants α1 = 15 and
α2 = 1.4 * 10−5. If it remains isolated, the number of plants will reach a maximum value,
which then? To be solved without ODE solvers the program!
But on the day when the plants reach 80% of the maximum value, four rats appear. Find
with Runge-Kutta method and step length 1 day (ie 1/365) which day it is. Check if
Runge-Kutta's method with the step length 1 week combined with appropriate interpolation
a
ives at the same day. Investigate what hour Runge-Kutta's method combined with
appropriate interpolation a
ives with the different ones step lengths 1 day and 1 week.
The rats thrive on the island and eat up a lot of the plants. If there are too many of them,
there will be a shortage of plants.
The interaction between number of rats and plants can be described by
dP/dt = α1*P−α2*P^2−α3*P*R
dR/dt=−β1*R^1.4+β2*P^0.6*R^0.8
where α3 = 0.014, β1 = 1.9, β2 = 0.083 If the island now remains isolated, the number of
plants and rats will reach constant equili
ium value, which ones? (Then use your ODE
eader program to check your calculation).
But, exactly 0.5 years after the rats a
ived, 2 cats landed. The cats eat the rats but not
the plants, so that the growth equation of the plants is not affected but the equation of the
ats (and of course the cats) is:
dR / dt = −β1*R^1.4 + β2*P^0.6R^0.8 - β3*RK
dK / dt = −γ1*K + γ2*R*K^0.5
where β3 = 1.5, γ1 = 1.9, γ2 = XXXXXXXXXXCalculate the number of the plants, rats and cats in
the different scenarios (ie when the plants came, when the rats came, when the cats came,
and when it was four years later. Plot the development for all the different varieties. How
much do the final values for each system differ from each other? (Then use your ODE solve
software to control your calculation of final values).
On the first day of the fifth year, the people come there and decide to fight the rats with
poison. 65% of the rats die, but also 17% of the cats. How many plants, rats and cats do you
have after 4.5 Years, ie in the middle of summer the fifth year? And at the end of the year?
Was the spray (of poison) worth it?
I would like written solutions for this: I want a solution for the equili
ium
points using precise and efficient methods (the method used here is not
efficient), an e
or estimation for every value in the tables below. Thank you in
advance.
Table 1
When plants
each 80% of
the maximum
value (step
When plants
each 80% of
the maximum
value (step
length one day)
When plants
each 80% of
the maximum
value (step
length one
week)
When plants
each 80% of
the maximum
value (step
length one
day)
The maximum
value of plants
unit day day hour hour plants
amount XXXXXXXXXX1071429
Table 2
Time When plants and
ats reach constant
equili
ium values
When 2 cats a
ived
(six months after the
ats)
4.5 years (six months afte
humans a
ived)
unit plants rats plants rats plants rats cats
amoun
t
XXXXXXXXXX XXXXXXXXXX
Here are the codes with the answers:
main code:
clear; clc; close all;
global alpha1 alpha2 alpha3 beta1 beta2 beta3 gamma1 gamma2
global t80DSta
alpha1 = 15; alpha2 = 1.4E-5; alpha3 = 0;
eta1 = 0; beta2 = 0; beta3 = 0;
gamma1 = 0; gamma2 = 0;
tol = 1E-1;
hD = 1/365; % One day time step
hW = 1/52.1429; % Week time step
P0 = 75; % Initial plants
R0 = 0; %Initial rats
C0 = 0; % Initial cats
y0 = [P0, R0, C0]';
Tf = 1;
tspan = [0, Tf];
[tD, yD] = RK4(@myode, tspan, y0, hD);
PStab = alpha1/alpha2;
idxStab = find(abs(yD(:,1) - PStab) <= tol,1);
tD(idxStab)
% Analytical solution for plant only system
%fSol = @(t) 1125*exp(15*t)/((exp(15*t)-1)* XXXXXXXXXX);
fprintf('\n\nWhen solved analytically, the plants grow stable to %10.3f plants in %4.3f years(%5.1f days), in absense of rats.\n\n',PStab,
tD(idxStab),tD(idxStab)*365);
% Find time at which 80% of this stable value is reached
idx80DStab = find((yD(:,1)-0.8*PStab) >= 0,1);
t80DStab = tD(idx80DStab);
[tW, yW] = RK4(@myode, tspan, y0, hW);
idx80WStab = find((yW(:,1)-0.8*PStab) >= 0,1);
t80DWStab = (tW(idx80WStab XXXXXXXXXX*PStab -
yW(idx80WStab-1,1))/(yW(idx80WStab,1)-yW(idx80WStab-1,1)))*(tW(idx80WStab)-tW(idx80WStab-1)));
fprintf('\n\nDay interpolation from week time step is %6.2f and from day time step is %6.4f. Clearly, they dont
match\n\n',t80DWStab*365,t80DStab*365);
% RATS come in
tspan = [t80DStab, 2];
P01 = yD(idx80DStab,1);
R01 = 4;
C01 = 0;
y01 = [P01, R01, C01]';
alpha3 = 0.014; beta1 = 1.9; beta2 = 0.083;
[tD1, yD1] = RK4(@myode, tspan, y01, hD);
Y = cat(1,yD(1:idx80DStab-1,:),yD1);
T = cat(1,tD(1:idx80DStab-1,:),tD1);
fprintf('\n\nWhen 4 rats are introduced, the plants stabilize to %8.2f and rats stabilize to %5.2f.\n\n',Y(end,1),Y(end,2));
% And the CATS come
idxp5yr = find(tD1-tD1(1)>=0.5,1);
tspan = [tD1(idxp5yr) 4];
P02 = yD1(idxp5yr,1);
R02 = yD1(idxp5yr,2);
C02 = 2;
y02 = [P02 R02 C02]';
eta3 = 1.5; gamma1 = 1.9; gamma2 = 0.0231;
[tD2, yD2] = RK4(@myode, tspan, y02, hD);
Y = cat(1,cat(1,yD(1:idx80DStab-1,:),yD1(1:idxp5yr-1,:)),yD2);
T = cat(1,cat(1,tD(1:idx80DStab-1,:),tD1(1:idxp5yr-1)),tD2);
% And the POISON is applied by evil humans
tspan = [T(end)+hD 4.5];
P03 = Y(end,1);
R03 = Y(end,2)*0.35;
C03 = Y(end,3)*0.83;
y03 = [P03 R03 C03]';
eta3 = 1.5; gamma1 = 1.9; gamma2 = 0.0231;
[tD3, yD3] = RK4(@myode, tspan, y03, hD);
Y = cat(1,Y,yD3);
T = cat(1,T,tD3);
str = 'Population dynamics plot';
figure('Numbertitle','off','Name',str);
plot(T,Y(:,1),'Color',[1,0,0],'LineWidth',2,'DisplayName','Plants');
hold on;
plot(T,Y(:,2),'Color',[0,1,0],'LineWidth',2,'DisplayName','Rats');
plot(T,Y(:,3),'Color',[0,0,1],'LineWidth',2,'DisplayName','Cats');
grid on;
legend;
xlabel('Time(years)');
ylabel('Species population');
title(str);
fprintf('\n\nClearly, the rats are again rising in population, which is why, the poison application was really not effective.\n\n');
differential equations :
function dy = myode(t,y)
%My differential equation system
global alpha1 alpha2 alpha3 beta1 beta2 beta3 gamma1 gamma2
P = y(1); % Plants
R = y(2); % Rats
C = y(3); % Cats
dP = alpha1*P - alpha2*power(P,2) - alpha3*P*R;
dR = -beta1*power(R,1.4) + beta2*power(P,0.6)*power(R,0.8) - beta3*R*C;
dK = -gamma1*C + gamma2*R*power(C,0.5);
dy = [dP, dR, dK]';
RK4:
function [t, y] = RK4(f, tspan, y0, h)
% RK4 with input step size
% Inputs:
% f: function handle of f(t, y)
% tspan: the time period for simulation (should be a 1x2 a
ay ...
% contain start time and end time)
% y0: the initial conditions for the differential equation
% Outputs:
% t: co
esponding time sequence as a T x 1 vecto
% y: the solution of the differential equation as a T x n matrix, ...
% where T is the number of time steps and n is the dimension of y
hk = h;
% Initialize e
o
e0 = 1E-4;
tk = tspan(1);
y = y0;
tf = tspan(2);
% Compute solutions by RK4 adpative step size algorithm
while tk(end) < tf
hk = min(hk, tf-tk(end));
k1 = f(tk(end), y(:,end));
k2 = f(tk(end)+0.5*hk,y(:,end)+0.5*hk*k1);
k3 = f(tk+0.5*hk,y(:,end)+0.5*hk*k2);
k4 = f(tk(end)+hk,y(:,end)+hk*k3);
ykp1 = y(:,end) + (1/6)*hk*k1 + (1/3)*hk*k2 + (1/3)*hk*k3 + (1/6)*hk*k4;
tk = cat(1,tk,tk(end)+hk);
y = cat(2,y,ykp1);
end
t = tk;
y = y';