Trying to model red and grey squirrels

  • Context: Mathematica 
  • Thread starter Thread starter LETS_GO
  • Start date Start date
  • Tags Tags
    Model population sir
Click For Summary

Discussion Overview

The discussion revolves around modeling the populations of red and grey squirrels using a system of differential equations. Participants explore how to represent the reproductive rate of grey squirrels as a function of time within the context of a mathematical model, focusing on the implementation of this change in a coding environment.

Discussion Character

  • Technical explanation
  • Mathematical reasoning
  • Debate/contested

Main Points Raised

  • One participant presents a model using differential equations to describe the dynamics of red and grey squirrel populations, including variables for susceptible, infected, and recovered individuals.
  • Another participant asks how to change the maximum reproductive rate of grey squirrels over time, indicating a need for a time-dependent function.
  • A suggestion is made to define the reproductive rate as a linear function of time, with constants to be determined.
  • Some participants express uncertainty about the use of NDSolve, with one noting that the current approach of defining aG as a function of the previous time step does not work in this context.
  • There is a proposal to set initial values for aG and define its decrease over time, but a participant points out the continuous nature of the problem, suggesting that aG should be expressed as a continuous function rather than a discrete one.
  • One participant expresses frustration and uncertainty about how to implement the desired changes to aG in the model.

Areas of Agreement / Disagreement

Participants do not reach a consensus on how to properly define the reproductive rate of grey squirrels as a function of time, with multiple competing views on the correct approach to take.

Contextual Notes

There are unresolved issues regarding the proper formulation of aG as a continuous function and the implications of using discrete versus continuous time in the model.

LETS_GO
Messages
4
Reaction score
0
TL;DR
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:
Physics news on Phys.org
What is your question?
 
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.
 
LETS_GO said:
how do we change aG with time
[CODE title="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??*)[/CODE]
That looks like the change that you want. I'm not familiar with NDsolve and don't see anything wrong with your code.
 
FactChecker said:
[CODE title="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??*)[/CODE]
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".)
 
Reply
  • Like
Likes   Reactions: FactChecker
DrClaude said:
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
 
LETS_GO said:
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