# Trying to model red and grey squirrels

• Mathematica
• LETS_GO

#### LETS_GO

TL;DR Summary
Hello, we are trying to change the max reproductive rate of grey squirrels. with a value starting at 1.2 and want it to go down to 0.3 with a step of 0.05.
Code:
ClearAll["Global*"]
(*R = reds, G = greys*)
(*S = susceptible, I = infected, R = recovered*)

tseries = {t, 0, 3};
vars = {HG[t], HR[t], SG[t], IG[t], RG[t], SR[t], IR[t], aG[t], qG[t]};

b = 0.4;    (*natural mortality rate, both species*)
\[Beta] = 0.7;    (*rate of virus transmission*)

aR = 1;       (*Reds max. reproductive rate*)
\[Alpha] = 26;       (*Reds mortaility rate due to virus*)
cR = 0.61;(*Reds competative effect on greys*)
KR = 60;   (*Reds carrying capacity, 5x5 km*)

\[Gamma] = 13;      (*Greys recovery rate from virus*)
cG = 1.65;(*Greys competative effect on reds*)
KG = 80;   (*Greys carrying capacity, 5x5 km*)

qR = (aR - b)/KR;  (*Reds crowding susceptibility*)
(*Greys crowding susceptibility*)

(*initial conditions*)
SGinit = 10;
IGinit = 2;
RGinit = 0;
SRinit = 60;
IRinit = 0;
HGinit = SGinit + IGinit + RGinit;
HRinit = SRinit + IRinit;

eqns =
(*total populations*)
{HG[t] == SG[t] + IG[t] + RG[t],
HR[t] == SR[t] + IR[t],

aG[t] == aG[t - 1] - 0.05,(*aG = Greys max. reproductive rate -
possible birth control??*)
qG[t] == (aG[t] - b)/KG,  (*Greys crowding susceptibility*)
(*SIR growth rates*)

SG'[t] == ((aG[t] - (qG[t]*(HG[t] + (cR*HR[t]))))*HG[t]) - (b*
SG[t]) - (\[Beta]*SG[t]*(IG[t] + IR[t])),
IG'[t] == (\[Beta] *SG[t]*(IG[t] + IR[t])) - (b*IG[t]) - (\[Gamma]*
IG[t]),
RG'[t] == (\[Gamma]*IG[t]) - (b*RG[t]),
SR'[t] == ((aR - (qR*(HR[t] + (cG*HG[t]))))*HR[t]) - (b*
SR[t]) - ((\[Beta]*SR[t])*(IR[t] + IG[t])),
IR'[t] == (\[Beta]*SR[t]*(IG[t] + IR[t])) - ((\[Alpha] + b)*IR[t]),

(*call initial conditions*)
HG[0] == HGinit, HR[0] == HRinit,
aG[0] == 1.2, qG[0] == (1.2 - b)/KG,
SG[0] == SGinit, IG[0] == IGinit, RG[0] == RGinit,
SR[0] == SRinit,
IR[0] ==
IRinit                                                            \
};

sol = NDSolve[eqns, vars, tseries];

{Plot[Evaluate[{SG[t], IG[t], RG[t], SR[t], IR[t]} /. sol], tseries,
ImageSize -> Large, PlotLegends -> {"SG", "IG", "RG", "SR", "IR"}],
Plot[Evaluate[{HG[t], HR[t]} /. sol], tseries, ImageSize -> Large,
PlotLegends -> {"HG", "HR"}]}
<Moderator's note: Please use CODE tags when posting code.>

Last edited by a moderator:

how do we change aG with time

You need to make it an actual function of time
Code:
aG[t] == a0 + a1 * t
with appropriate constants a0 and a1.

how do we change aG with time
It looks like it already is changing with each time step on line 37::
aG[t] == aG[t - 1] - 0.05,(*aG = Greys max. reproductive rate -
possible birth control??*)
That looks like the change that you want. I'm not familiar with NDsolve and don't see anything wrong with your code.

It looks like it already is changing with each time step on line 37::
aG[t] == aG[t - 1] - 0.05,(*aG = Greys max. reproductive rate -
possible birth control??*)
That looks like the change that you want. I'm not familiar with NDsolve and don't see anything wrong with your code.
That doesn't work, because NDSolve does not iterate the solution like that. It needs to be a proper function of time. (And note that t is a real variable, not an index, so t - 1 means "one unit of time earlier".)

FactChecker
You need to make it an actual function of time
Code:
aG[t] == a0 + a1 * t`
with appropriate constants a0 and a1.

So I want my starting value for aG to be 1.2 and I want it to go down to at least 0.35 with a decrease of 0.5 for each time step

So I want my starting value for aG to be 1.2 and I want it to go down to at least 0.35 with a decrease of 0.5 for each time step
But there is no time step! The problem is not fundamentally discrete, but continuous. You need to figure out what aG is as a function of time.

hmmm I have tried the last couple of days, and im really not sure how to do this