Notebook of fun Chapter 2 things

Getting slightly tabular output

Try running this S - i - R code and see how much easier the output is to read!  Thank Chris Gole for thinking of this.  
Notes: 
(a) N[t] gives a decimal value (it stands for "numerical") instead of a fraction.
(b) The "       " keep spaces between the columns of output.

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
Print["t",  "      S   ","       i   ", "      R  "];
S =45400;
i = 2100;
R = 2500;
tinitial=0;
tfinal = 1;
t = tinitial;
numberofsteps = 10;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat; 
deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
Print[N[t],"   ", S,"    ",i,"    ", R];
 k++]

Problem 1 (Stabilization or lack thereof)

For our usual measles epidemic, construct a table of approximations for S(13), I(13), R(13) in which numberofsteps takes the values 10, 100, 1000, 10000, 100000. At each stage, indicate how many digits of S(13), I(13), R(13) have stabilized.  
Here is some convenient code for you to use... by the way, you can preserve the output from a calculation by copying it and then pasting it into a different cell.

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
S =45400;
i = 2100;
R = 2500;
tinitial = fill me in;
tfinal = fill me in;
t = tinitial;
numberofsteps = fill me in;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat; 
deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
 k++]
Print["steps","     t",  "      S   ","       i   ", "      R  "];
Print[numberofsteps,"      ",t,"   ", S,"    ",i,"    ", R];

Problem 2 (Survivors)

Let 200 days pass, so t = 200. Determine the values of I(200) and
S(200) to one decimal place accuracy by constucting tables of approximations in
which you let numberofsteps increase by powers of 10 (as in Problem 1).  Here's the code again...

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
S =45400;
i = 2100;
R = 2500;
tinitial = fill me in;
tfinal = fill me in;
t = tinitial;
numberofsteps = fill me in;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat; 
deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
 k++]
Print["steps","     t",  "      S   ","       i   ", "      R  "];
Print[numberofsteps,"      ",t,"   ", S,"    ",i,"    ", R];

How do you interpret the values of S, i and R when numberofsteps = 10?

Are there any “survivors”, that is, susceptibles who never become ill before the illness disappears?  If so, how many are there (measured in a whole number of persons)?  State clearly the evidence on which you base your assertions.  This last part is crucial!

Was it necessary to go out 200 days to answer this question?  That is, can you argue that the infection disappears in even fewer days?  Or does it take even more than 200 days?

Problem 3 (Survivors compared to susceptibles)

Do you suspect that having a larger susceptible population will increase or decrease the number of survivors?   What about a smaller susceptible population?  
Explore how the eventual number of survivors depends on the number of susceptibles at the outset. That is, 
keep I(0) = 2100 and R(0) = 2500, but vary the number of susceptibles; let S(0) = 4540 and determine the number of survivors, then let S(0) = 20000, and then let S(0) = 80000, and then decide how many and which other values of S(0) to test.

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
S =fill me in;
i = 2100;
R = 2500;
tinitial = fill me in;
tfinal = fill me in;
t = tinitial;
numberofsteps = fill me in;
deltat = (tfinal- tinitial)/numberofsteps;
For[k=1,k≤numberofsteps,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat; 
deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
 k++]
Print["t",  "      S   ","       i   ", "      R  "];
Print[t,"   ", S,"    ",i,"    ", R];

Roughly speaking, does the number of survivors depend directly or inversely on the initial number of susceptibles? How does this compare with your earlier conjecture about this question?

Problem 4 (Peak of infection rate)

Below is a version of the S-i-R code which will run from tinitial to tfinal as long as k ≤ numberofsteps and iprime is positive.  Run the code.

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
S =45400;
i = 2100;
R = 2500;
tinitial = 0;
tfinal = 10;
t = tinitial;
numberofsteps = 1000;
deltat = (tfinal- tinitial)/numberofsteps;
For[iprime=1;k=1,k≤numberofsteps && iprime>0,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat; 
deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
 k++]
Print["t",  "      S   ","       i   ", "      R  "];
Print[t,"   ", S,"    ",i,"    ", R];

Does this tell you when i is at its highest point?  Why or why not?
What stopped the code from running?  
(Hint: there are three possibilities...  (1) k > numberofsteps (2) iprime ≤ 0 (3) t=10. Which is it?)
If this code didn't tell you when i is at its peak, then modify it so it does tell you when i is at its peak.
What is the value of i at its peak?  What is the value of S when i is at its peak?

Problem 5 (Epidemic Duration)

Here's another copy of the code.  Modify this code to make it show you when the epidemic is over, that is, when the number of infected people is less than 1.  What value of t do you get as the duration of the epidemic?

Clear[t];
Clear[tinitial];
Clear[tfinal];
Clear[numberofsteps];
Clear[S];
Clear[i];
Clear[R];
Clear[deltat];
S =45400;
i = 2100;
R = 2500;
tinitial = 0;
tfinal = 10;
t = tinitial;
numberofsteps = 1000;
deltat = (tfinal- tinitial)/numberofsteps;
For[iprime=1;k=1,k≤numberofsteps && iprime>0,
Sprime = -.00001*S*i; 
iprime =.00001*S*i - i/14;
 Rprime = i/14; 
deltaS = Sprime*deltat; 
deltai = iprime*deltat; 
deltaR = Rprime*deltat; 
t=t+deltat; 
S = S+deltaS; 
i = i + deltai; 
R = R+deltaR; 
 k++]
Print["t",  "      S   ","       i   ", "      R  "];
Print[t,"   ", S,"    ",i,"    ", R];