The S-I-R Model: A Comaprison of
Euler's Recursion with NDSolve[.]

by
K. D.  Stroyan
University of Iowa

copyright 2007 by Keith Stroyan - All rights reserved.

Notebook Overview

The purpose of this notebook is to compare the recursive solution with Mathematica's built-in solver NDSolve[.].  The recursive method is explained for the second S-I-R model in the program SecondSIR (also see the NoteBook EulerApprox).  The NDSolve[.] method is explained in more detail in SIRsolver.ma

References to the Text

In Section 2.3 of the main text we derive the continuous model of an epidemic as a system of differential equations.

Step-by-step Computation

The next  4 subsections show you how to compare solutions done by the recursive method of Chapter 2  and with the built-in Mathematica solver for differential equations.   You must enter the computations in order from frist to last.

The Differential Equations

Enter the next cell to define the differential equations in the Mathematica kernel:

Clear[s,i,t,soln];
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]]];

Euler's Method for the Continuous Variable Model

Enter the next cell to solve the differential equations by "Euler's Method" or the recursive method you used by hand and in the SecondSIR program:

Print["Spread of rubella in fractions of a population."];
Print["Contact number, c = ",c,", Infectious period 1/b =",(1/b)];
Print["Daily contact rate, a =",a];
t0=0.`;
s0=(1/(3.`));
i0=(1/(30.`));
initialValues={t0,s0,i0};
dt=N[(1/2)];
finalTime=N[90];
steps=Floor[(finalTime/dt)];
Print["Increments of dt = ",dt," for ",steps," recursive steps."];
G[{t_,s_,i_}]:={t+dt,s+sRate[s,i] dt,i+iRate[s,i] dt};

solnList=NestList[G,initialValues,steps];
suscs=Flatten[Take[Transpose[solnList],{2}],1];
times=Flatten[Take[Transpose[solnList],{1}],1];
infes=Flatten[Take[Transpose[solnList],-1],1];
sPairs=Transpose[{times,suscs}];
iPairs=Transpose[{times,infes}];

eulerPlot=ListPlot[sPairs,Joined->True,AxesLabel->{"t","s"},PlotStyle->{Red}]

ListPlot[iPairs,Joined->True,AxesLabel->{"t","i"},PlotStyle->{Red}]

Print["The last susceptible value:"];Spair=Flatten[Take[sPairs,-1]];
eulerSdata=Spair[[2]];
Print[eulerSdata];

The NDSolve[.] Command

Mathematica has a built-in initial value solver that works as follows.

Enter this computation

startUchange={s[0]==(1/(3.`)),i[0]==(1/(30.`)),s'[t]==sRate[s[t],i[t]],i'[t]==iRate[s[t],i[t]]};
finalTime=90;

soln=NDSolve[startUchange,{s[t],i[t]},{t,0,finalTime}];
{s[t_],i[t_]}={s[t]/.soln,i[t]/.soln};

accuratePlot=Plot[s[t],{t,0,finalTime},AxesLabel->{"t","s"},PlotStyle->{Blue}]

Plot[i[t],{t,0,finalTime},AxesLabel->{"t","i"},PlotStyle->{Blue}]

Print["Spread of rubella in fractions of a population."];
Print["Contact number, c = ",c,", Infectious period 1/b =",(1/b)];

Print["Daily contact rate, a =",a];
Print["NDSolve[.] value of s[finalTime] =",s[finalTime]];

The Comparison

The next cell compares the results of the two computations of  s[t]:

Show[{eulerPlot,accuratePlot},PlotLabel->"Solutions Together"]
Show[GraphicsRow[{eulerPlot,accuratePlot}]]
Print["t = ",finalTime];
Print[" E-Approx. = ",eulerSdata," for dt = ",dt];
Print[" Accurate Approx. = ",s[finalTime]];

Exercises

Exercise 2.3.3

Run the recursive Euler's method program with step sizes of  dt = 1, 1/2, 1/4, 1/8  over the time interval from   t = 0 to  t = 120 and compare each of these solutions with the NDSolve[.] solution over the time interval from  t = 0 to  t = 120.  (Make copies of both programs and the comparison program  above and edit them to compute for the different  dt  values.)