NDSolve[.] - solution of The S-I-R Problem

by
K. D.  Stroyan
University of Iowa

copyright 2007 by Keith Stroyan - All rights reserved.

Notebook Overview

This NoteBook uses NDSolve[.] to solve the S-I-R differential equations.  This is a fast reliable way to find numerical solutions to differential equations with an initial condition.  Remember one of the most important basic ideas in this course is that if you know:

1)	where a quantity starts, such as  s[0] = 1,

and

2)	how the quantity changes, such as  (ds/dt)= s ,

then you can determine the value of the quantity  s[t]  at later times.  It is not always easy to find  s[t], but the computer can help.

In the case of the S-I-R model, we use starting values for the susceptible and infectious fractions and the model equations,

s[0] = 1/3
i[0] = 1/30

(ds/dt) = - a s i		or	s-time-rate = s'[t] = -a s i

(di/dt) = a s i - b i	or	i-time-rate = i'[t] = a s i - b i

References to the Text

The S-I-R differential equations are derived and explained in Section 2.3 of the main text.  The variables are:

t = time in days from  t = 0  at the beginning of observation

s = the fraction of the population that is susceptible to infection

i = the fraction of the population that is infectious

The parameters  c, b, a, and  n are expalined in Chapter  2.  They are determined by a particular disease and population.

NDSolve[.] Explained for a simple differential equation

Throughout the course we will encounter the extremely important general problem: Given a starting value and a rule for how a quantity changes, find the quantity at later times.  The following specific instance has the simplest possible change rule in the dependent variable:

1)	s[0] = 1
and
2)	(ds/dt)= s

This initial value problem (or differential equation) has the analytical solution

s[t] = e^t = Exp[t]

The next program solves the problem numerically for the starting information (1) and the chage information (2).

Clear[s,t,soln];
sRate[s_]:=s;
startUchange={s[0]==1,s'[t]==sRate[s[t]]};
finalTime=1;
soln=NDSolve[startUchange,{s[t]},{t,0,finalTime}];
Plot[s[t]/.soln,{t,0,finalTime},PlotStyle->{Thickness[0.01`],Red}],Table[{t,s[t],N[(E^t)]}/.soln,{t,0,finalTime,0.1`}]

The computer solution "soln" to the above program is a list of "Interpolating Functions."  This means that Mathematica solves the differential equation at certain times and can interpolate to other times.  The program shows how to extract graphs and values from this funny form of the answer using "/.":

Mathematica's Description of NDSolve[.]

Enter the following command to see Mathematica's explaination of NDSolve[.]

??NDSolve

Step-by-step Computation

Define the Differential Equations as Rate Functions

Enter the next cell to define equations in Mathematica's kernel:

Clear[s,i,t,soln,sRate,iRate];
c=6.8;
b=(1/11.0);
a=b c;
sRate[s_,i_]:=-a s i;
iRate[s_,i_]:=(a s-b) i;
Print["The differential equations for rubella:"];
Print["  s'[t] = ",sRate[s[t],i[t]]];
Print["  i'[t] = ",iRate[s[t],i[t]]];

Use NDSolve

To use NDSolve[.], we need a list of equations with the start and change laws, as well as a time interval for the solution.  The start and change is called "startUchange" in the next cell.  Enter it now and notice the peculiar output.

startUchange={s[0]==(1/3.),i[0]==(1/30.),s'[t]==sRate[s[t],i[t]],i'[t]==iRate[s[t],i[t]]};
finalTime=60;
soln=NDSolve[startUchange,{s[t],i[t]},{t,0,finalTime}];
s[t_]=s[t]/.soln[[1,1]]
i[t_]=i[t]/.soln[[1,2]]

Plot and evaluate the solution

Interpolating functions are explained in the NoteBook Overview.  Enter the next cell to see the solutions.  The explanation of "/." is in the NoteBook Overview section above.

Plot[s[t],{t,0,finalTime},AxesOrigin->{0,0},PlotStyle->{Blue}]
Plot[i[t],{t,0,finalTime},AxesOrigin->{0,0},PlotStyle->{Red}]
Table[{t,s[t],i[t]}/.soln,{t,0,finalTime,10}]

A Combined Computation

Clear[s,i,t,soln,sRate,iRate];
c=5.2;
b=(1/16.);
a=b c;
Print["Spread of smallpox in a large fixed population."];
Print["Contact number, c = ",c,", Infectious period 1/b =",(1/b)];
Print["Daily contact rate, a =",Round[a,.005]];

sRate[s_,i_]:=-a s i;
iRate[s_,i_]:=(a s-b) i;

Manipulate[
soln=NDSolve[{s[0]==1-R0,i[0]==.05,s'[t]==sRate[s[t],i[t]],i'[t]==iRate[s[t],i[t]]},{s[t],i[t]},{t,0,finalTime}];
Plot[{s[t]/.soln[[1,1]],i[t]/.soln[[1,2]]},{t,0,finalTime},AxesLabel->{"t","susceptible & infectious"},AxesOrigin->{0,0},PlotRange->{0,1},PlotStyle->{Blue,Red},PlotLabel->StringForm["Initial fraction immune =`1`, final susceptible s[`2`] = `3`",R0,finalTime,Round[s[t]/.soln[[1,1]]/.t->finalTime,.0005]]]
,{R0,.05,.95},{finalTime,60,120}]

Practice Exercise

Now you try.  Copy, paste and edit the program above, changing the disease parameters to  c = 1.4,  b = 1/3, finalTime = 120.

Exercises

Problem 2.2

Copy and re-run the "Combined Computation" program above with parameters b and c as follows.

(a)

A rubella epidemic, where  1/b = 11, c = 6.8 and a = b c.

(b)

Measles, with a very high contact number c=15 and b=1/8.

(c)

Influenza, with a low contact number of c=1.4 and b=1/3.

Make the final time large enough that you feel the epidemic has run its course.

The most important part of this problem is to analyse your results.

What are the main differences between these three diseases?  
	
	Which one is highly contagious, which is medium and which least contagious?  
	
	Which one allows a larger portion of the population to escape infection?  
	
	How many escape in each case?

Write your summary analysis of the comparison between these diseases in a few brief paragraphs.  (You may use a text cell in your copy of this program.)