Systems of Differential Equations

A Predator-Prey Model

This is a model for what fraction of the MIT population would know a particular joke, the joke is told a fraction of the time when a person who has heard it meets someone who hasn't. The rate of joke spreading is proportional to how funny it is.

It is a closely held secret that Chuck Vest, MIT's president, has an odd sense of humor--he will not let anyone graduate who cannot tell him the joke on graduation day.  This often is a source of embarassment and confusion to the happy families.

JokeFunniness  =  0.66;
JokeSpreadRate[Naive_ ,  Jaded_]  := 
      JokeFunniness*Naive*Jaded

Some people who have heard the joke forget it, depending on how forgetful they are:

PopulationForgetfullness  =  0.2;
ForgotJoke[Naive_ ,  Jaded_]  := 
      PopulationForgetfullness*Jaded

A certain fraction of new students added to the population each day who presumably don't know the joke, the graduation rate is proportional to the entrance rate and the fraction of those who know the joke.

FreshmanEntranceRate  = 0.24;
GraduationRate[Naive_ ,  Jaded_]   =  
        FreshmanEntranceRate*Jaded;

The fraction of students who remain naive tomorrow is

TomorrowsNaive[TodaysNaive_, 
        TodaysJaded_]  :=  TodaysNaive + (FreshmanEntranceRate/365) - JokeSpreadRate[
              TodaysNaive ,  TodaysJaded]/365

The fraction of students who will know the joke tomorrow is

TomorrowsJaded[TodaysNaive_, TodaysJaded_]  :=  
      TodaysJaded  +  
        JokeSpreadRate[TodaysNaive ,  TodaysJaded]/365  -  
        ForgotJoke[TodaysNaive ,  TodaysJaded]/365   -  
        GraduationRate[TodaysNaive ,  TodaysJaded]/365

The population evolves each day as

TomorrowsPopulation[{TodaysNaive_, TodaysJaded_}] := {TomorrowsNaive[
          TodaysNaive ,  TodaysJaded], 
        TomorrowsJaded[TodaysNaive ,  TodaysJaded]}

Calculate the population after 5000 days from now if currently 1/5 of the students are naive and 4/5 of the students are in the know.

NestList[TomorrowsPopulation, { .2,  .8}, 5000] // Short

ListPlot[NestList[TomorrowsPopulation, { .24,  .75}, 20000], 
      PlotRange -> All, PlotJoined -> True, 
      PlotStyle -> {Hue[Random[]], Thickness[0.005]}, 
      AxesLabel -> {"Naive", "Jaded"}]

graphicslist  = 
      Table[ListPlot[
          NestList[TomorrowsPopulation, {Random[], Random[]}, 20000], 
          PlotRange -> All, PlotJoined -> True, 
          PlotStyle -> {Hue[Random[]], Thickness[0.005]}, 
          DisplayFunction -> Identity], {i, 1, 10}]

Show[graphicslist, DisplayFunction -> $DisplayFunction, 
      AxesLabel -> {"Naive", "Jaded"}]

Analysis of Stability of the Predator-Prey Model

The model above can be written as the system of ODEs:
(dN/dt) = α - ρ N J
(dJ/dt) = -αJ  +  ρ N J  - ϕJ

Find the fixed points:

Ndot  =  α - ρ N J
Jdot  =  -α J + ρ N J - ϕ J

fixedpoint  =  
      Solve[{Ndot == 0,  Jdot == 0}, {N, J}] // 
        Flatten

Define the general Jacobian

Jacob  =  {
{D[Ndot, N],  
            D[Ndot, J]}, 
{D[Jdot, N],  
            D[Jdot, J]}
};
Jacob // MatrixForm

Find the Jacobian at the fixed point:

JacobFixedPoint  =  
        Simplify[Jacob /. fixedpoint];
JacobFixedPoint // MatrixForm

evals  =  Eigenvalues[JacobFixedPoint]

The Eigenvalues will always have a negative real part because the input parameters are positive real coefficients

Considering the parameters used above, the solution is attracted to the fixed point but oscillates

MitEvals  =  
      evals /. {α ->  FreshmanEntranceRate,  ρ  ->  
            JokeFunniness,  ϕ  ->  
            PopulationForgetfullness}

The Oscillation perioid is about 18 years

2 π /Im[MitEvals[[1]]]