Trying to model red and grey squirrels

[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.

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.>

 

[DrClaude]

What is your question?

 

[LETS_GO]

how do we change aG with time

 

[DrClaude]

You need to make it an actual function of time

aG[t] == a0 + a1 * t

with appropriate constants a0 and a1.

 

[FactChecker]

how do we change aG with time

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.

 

[DrClaude]

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".)

 

[LETS_GO]

You need to make it an actual function of time

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

 

[DrClaude]

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.

 

[LETS_GO]

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