The SIR Model for Spread of Disease

Part 1. Background: Hong Kong Flu

There are no questions to be answered in Part 1.

Part 2. The Differential Equation Model

Answer the questions in Part 2 in the text cell below.

1.

(Step 2 is to be done on paper.)

3.

4.

5.

6.

7.

Part 3. Euler's Method for Systems

In this part, we first set up command lines for generating numerical and graphical solutions.  After the last set of command lines, you will find text lines in which you can write your answers to questions.

Enter sample values of b and k.

b  =  1/2;
k  =  1/3;

Enter the next block of commands to define the SIR system of differential equations.  Mathematica's notation for (ds/dt) is D[s[t], t], and similarly for the other derivatives.

DE1  =  D[s[t], t]  ==  -b*s[t]*i[t];
DE2  =  D[i[t], t]  ==  b*s[t]*i[t] - k*i[t];
DE3  =  D[r[t], t]  ==  k*i[t];
SIRsys  =  {DE1, DE2, DE3};

The next command block draws the solution curves for the SIR model. Each scene specifies one function of t.  The three separate curves are overlaid with the Show command.  The inputs start with a list of the differential equations.  Next comes the list[s(t),i(t),r(t)] identifying the dependent variables,then the range for t;  and a list of initial conditions,and finally any optional parameters.

initial  =  {s[0]  ==  1, 
          i[0]  ==  1.27*10^(-6), 
          r[0]  ==  0};
solution  =  
      Flatten[NDSolve[
          Join[SIRsys, initial], {s, i, r}, {t, 0, 150}]]; 
splot  =  
      Plot[Evaluate[s[t] /. solution[[1]]], {t, 0, 150}, 
        PlotStyle -> RGBColor[0, 0, 1], 
        DisplayFunction -> Identity]; 
iplot  =  
      Plot[Evaluate[i[t] /. solution[[2]]], {t, 0, 150}, 
        PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction -> Identity]; 
rplot  =  
        Plot[Evaluate[r[t] /. solution[[3]]], {t, 0, 150}, 
          PlotStyle -> RGBColor[0, 1, 0], 
          DisplayFunction -> Identity];
Show[splot, iplot, rplot, AxesLabel -> {"t", "s,i,r"}, 
        DisplayFunction -> $DisplayFunction];

Now we repeat our calculations of solution functions, but using Euler's Method instead of whatever Mathematica did before.  Because Mathematica does not have a built-in function for Euler's Method, we have written one below.  The same number of steps, n, is used in all three solutions.

slope[w_] :=  {-b*w[[1]]*w[[2]], 
          b*w[[1]]*w[[2]] - k*w[[2]], k*w[[2]]};
n  =  20;
Delta  =  150/n;
Clear[t];
t[k_] :=  k*Delta;
tvalues  =  Table[t[k], {k, 0, n}];
Clear[s, i, r];
sir[0] =  {1, 1.27*10^(-6), 0};
sir[k_] :=  sir[k]  =  
          sir[k - 1] + N[slope[sir[k - 1]]]*Delta;
Table[sir[k], {k, 0, n}];
data  =  
        Table[Thread[{tvalues, %[[All, i]]}], {i, 3}];
splotEuler  =  
        ListPlot[data[[1]], 
          PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.015]}, 
          DisplayFunction -> Identity];
iplotEuler  =  
        ListPlot[data[[2]], 
          PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.015]}, 
          DisplayFunction -> Identity];
rplotEuler  =  
        ListPlot[data[[3]], 
          PlotStyle -> {RGBColor[0, 1, 0], PointSize[0.015]}, 
          DisplayFunction -> Identity];
Show[splotEuler, iplotEuler, rplotEuler, PlotRange -> {0, 1}, 
        AxesLabel -> {"t", "s,i,r"}, 
        DisplayFunction -> $DisplayFunction];
Clear[s, i, r];

To get a sense of the accuracy (or lack of accuracy) of the Euler solutions, plot both sets of solutions together:

Show[splot, iplot, rplot, splotEuler, iplotEuler, rplotEuler, 
        PlotRange -> {0, 1}, 
        DisplayFunction -> $DisplayFunction];

Answer Part 3 questions in the text cell below.

2.

3.

4.

Part 4. Relating Model Parameters to Data

Here too you will find a text cell for your responses immediately after the computational and graphic commands.  Use the following block of commands for your experimentation.

Clear[s, i, r];
b  =  1/2;
k  =  1/3;
DE1  =  D[s[t], t]  ==  -b*s[t]*i[t];
DE2  =  D[i[t], t]  ==  b*s[t]*i[t] - k*i[t];
DE3  =  D[r[t], t]  ==  k*i[t];
SIRsys  =  {DE1, DE2, DE3};

initial  =  {s[0]  ==  1, 
          i[0]  ==  1.27*10^(-6), 
          r[0]  ==  0};
solution  =  
        Flatten[NDSolve[
            Join[SIRsys, initial], {s, i, r}, {t, 0, 150}]];
iplot  =  
        Plot[Evaluate[i[t] /. solution[[2]]], {t, 0, 150}, 
          PlotStyle -> RGBColor[1, 0, 0], 
          AxesLabel -> {"t", "i[t]"}];

Answer Part 4 questions in the text cell below.

1.

2.

3.

4.

5.

6.

7.

Part 5. The Contact Number

Put your answers here.  If you need Mathematica's help with any computations, use the input cell provided below or feel free to create your own.

1.

2.

3.

4.

5.

Part 6. Herd Immunity

Put your answers here.  If you need Mathematica's help with any computations, use the input cell provided below or feel free to create your own.

1.

2.

3.

4.

5.

6.

Part 7. Summary

1.

2.

3.

4.

5.

6.

7.

8.

The SIR Model for Spread of Disease

Part 1. Background: Hong Kong Flu

There are no questions to be answered in Part 1.

Part 2. The Differential Equation Model

Answer the questions in Part 2 in the text cell below.

1.

(Step 2 is to be done on paper.)

3.

4.

5.

6.

7.

Part 3. Euler's Method for Systems

In this part, we first set up command lines for generating numerical and graphical solutions.  After the last set of command lines, you will find text lines in which you can write your answers to questions.

Enter sample values of b and k.

b  =  1/2;
k  =  1/3;

Enter the next block of commands to define the SIR system of differential equations.  Mathematica's notation for (ds/dt) is D[s[t], t], and similarly for the other derivatives.

DE1  =  D[s[t], t]  ==  -b*s[t]*i[t];
DE2  =  D[i[t], t]  ==  b*s[t]*i[t] - k*i[t];
DE3  =  D[r[t], t]  ==  k*i[t];
SIRsys  =  {DE1, DE2, DE3};

The next command block draws the solution curves for the SIR model. Each scene specifies one function of t.  The three separate curves are overlaid with the Show command.  The inputs start with a list of the differential equations.  Next comes the list[s(t),i(t),r(t)] identifying the dependent variables,then the range for t;  and a list of initial conditions,and finally any optional parameters.

initial  =  {s[0]  ==  1, 
          i[0]  ==  1.27*10^(-6), 
          r[0]  ==  0};
solution  =  
      Flatten[NDSolve[
          Join[SIRsys, initial], {s, i, r}, {t, 0, 150}]]; 
splot  =  
      Plot[Evaluate[s[t] /. solution[[1]]], {t, 0, 150}, 
        PlotStyle -> RGBColor[0, 0, 1], 
        DisplayFunction -> Identity]; 
iplot  =  
      Plot[Evaluate[i[t] /. solution[[2]]], {t, 0, 150}, 
        PlotStyle -> RGBColor[1, 0, 0], 
        DisplayFunction -> Identity]; 
rplot  =  
        Plot[Evaluate[r[t] /. solution[[3]]], {t, 0, 150}, 
          PlotStyle -> RGBColor[0, 1, 0], 
          DisplayFunction -> Identity];
Show[splot, iplot, rplot, AxesLabel -> {"t", "s,i,r"}, 
        DisplayFunction -> $DisplayFunction];

Now we repeat our calculations of solution functions, but using Euler's Method instead of whatever Mathematica did before.  Because Mathematica does not have a built-in function for Euler's Method, we have written one below.  The same number of steps, n, is used in all three solutions.

slope[w_] :=  {-b*w[[1]]*w[[2]], 
          b*w[[1]]*w[[2]] - k*w[[2]], k*w[[2]]};
n  =  20;
Delta  =  150/n;
Clear[t];
t[k_] :=  k*Delta;
tvalues  =  Table[t[k], {k, 0, n}];
Clear[s, i, r];
sir[0] =  {1, 1.27*10^(-6), 0};
sir[k_] :=  sir[k]  =  
          sir[k - 1] + N[slope[sir[k - 1]]]*Delta;
Table[sir[k], {k, 0, n}];
data  =  
        Table[Thread[{tvalues, %[[All, i]]}], {i, 3}];
splotEuler  =  
        ListPlot[data[[1]], 
          PlotStyle -> {RGBColor[0, 0, 1], PointSize[0.015]}, 
          DisplayFunction -> Identity];
iplotEuler  =  
        ListPlot[data[[2]], 
          PlotStyle -> {RGBColor[1, 0, 0], PointSize[0.015]}, 
          DisplayFunction -> Identity];
rplotEuler  =  
        ListPlot[data[[3]], 
          PlotStyle -> {RGBColor[0, 1, 0], PointSize[0.015]}, 
          DisplayFunction -> Identity];
Show[splotEuler, iplotEuler, rplotEuler, PlotRange -> {0, 1}, 
        AxesLabel -> {"t", "s,i,r"}, 
        DisplayFunction -> $DisplayFunction];
Clear[s, i, r];

To get a sense of the accuracy (or lack of accuracy) of the Euler solutions, plot both sets of solutions together:

Show[splot, iplot, rplot, splotEuler, iplotEuler, rplotEuler, 
        PlotRange -> {0, 1}, 
        DisplayFunction -> $DisplayFunction];

Answer Part 3 questions in the text cell below.

2.

3.

4.

Part 4. Relating Model Parameters to Data

Here too you will find a text cell for your responses immediately after the computational and graphic commands.  Use the following block of commands for your experimentation.

Clear[s, i, r];
b  =  1/2;
k  =  1/3;
DE1  =  D[s[t], t]  ==  -b*s[t]*i[t];
DE2  =  D[i[t], t]  ==  b*s[t]*i[t] - k*i[t];
DE3  =  D[r[t], t]  ==  k*i[t];
SIRsys  =  {DE1, DE2, DE3};

initial  =  {s[0]  ==  1, 
          i[0]  ==  1.27*10^(-6), 
          r[0]  ==  0};
solution  =  
        Flatten[NDSolve[
            Join[SIRsys, initial], {s, i, r}, {t, 0, 150}]];
iplot  =  
        Plot[Evaluate[i[t] /. solution[[2]]], {t, 0, 150}, 
          PlotStyle -> RGBColor[1, 0, 0], 
          AxesLabel -> {"t", "i[t]"}];

Answer Part 4 questions in the text cell below.

1.

2.

3.

4.

5.

6.

7.

Part 5. The Contact Number

Put your answers here.  If you need Mathematica's help with any computations, use the input cell provided below or feel free to create your own.

1.

2.

3.

4.

5.

Part 6. Herd Immunity

Put your answers here.  If you need Mathematica's help with any computations, use the input cell provided below or feel free to create your own.

1.

2.

3.

4.

5.

6.

Part 7. Summary

1.

2.

3.

4.

5.

6.

7.

8.