Solve SIR model -C++/MatLab/Python

1) C++
------------------------------------------------------
/* Solve simple SIR Model
Author: Jiansen Lu
Date: June 25, 2010
Description:
dS/dt = - beta*S*I
dI/dt = beta*S - gamma*I
dR/dt = gamma*I
where S: Susceptibel
I: Infectious
R: Recovered
beta: transmission rate
gamma: recover rate:
set initial condition S(0)=1e-6, I(0)=1.0-S(0)
beta = 1.5, gamma = 0.2
Algorithm: using Runge-Kutta method
dy/dt = f(t,y), y(t_0) = y_0
t_(n+1) = t_n+h
y_(n+1) = y_n + (k_1+2*k_2+2*_k3+k_4)/6
where: k_1=f(t_n,y_n), k_2=f(t_n+h/2,y_n+h*k_1/2)
k_3=f(t_n+h/2,y_n+h*k_2/2), k_4=f(t_n+h,y_n+h*k_3)
*/

#include<iostream>

#include<cmath>

using namespace std;

class SIR{
private:

double t,S,I,R,Pop[3];
double dPop[3],step;

double beta, gamma;
double tmax;
public:


SIR();
SIR(double beta0, double gamma0, double step0, double S00, \
          double I00, double tmax0);
~SIR();

void Diff(double Pop[3]);
void Runge_Kutta();

void Solve_Eq();
};
// Initialise the equations and Runge-Kutta integration
SIR::SIR(double beta0, double gamma0, double step0,double S00, \
       double I00, double tmax0)
{

beta = beta0;
gamma =gamma0;
step = step0;

S =S00;
I = I00;
R = 1 - S - I;

tmax = tmax0;
}
SIR::~SIR(){
cout <<"delete SIR"<<endl;
}

void SIR::Diff(double Pop[3])
{

// The differential equations


dPop[0] = - beta*Pop[0]*Pop[1];              // dS/dt

  dPop[1] = beta*Pop[0]*Pop[1] - gamma*Pop[1];   // dI/dt

  dPop[2] = gamma*Pop[1];                    // dR/dt

}

void SIR::Runge_Kutta(){
int i;
double dPop1[3], dPop2[3], dPop3[3], dPop4[3];

double tmpPop[3], initialPop[3];

/* Integrates the equations one step, using Runge-Kutta 4
Note: we work with arrays rather than variables to make the
coding easier */

initialPop[0]=S; initialPop[1]=I; initialPop[2]=R;

Diff(initialPop);
for(i=0;i<3;i++)
{

dPop1[i]=dPop[i];
tmpPop[i]=initialPop[i]+step*dPop1[i]/2;
}

Diff(tmpPop);
for(i=0;i<3;i++)
{


dPop2[i]=dPop[i];
tmpPop[i]=initialPop[i]+step*dPop2[i]/2;
}

Diff(tmpPop);
for(i=0;i<3;i++)
{

dPop3[i]=dPop[i];
tmpPop[i]=initialPop[i]+step*dPop3[i];
}

Diff(tmpPop);

for(i=0;i<3;i++)
{

dPop4[i]=dPop[i];
tmpPop[i]=initialPop[i]+(dPop1[i]/6 + dPop2[i]/3 + dPop3[i]/3 + dPop4[i]/6)*step;
}


S=tmpPop[0]; I=tmpPop[1]; R=tmpPop[2];


}

void SIR::Solve_Eq(){
t=0;
cout <<"t    S    I       R"<<endl;

do
{
Runge_Kutta();
t+=step;
cout<<t<<"   "<<S<<"   "<<I<<"   "<<R<<"   "<<endl;
}

while(t<tmax);


}
int main(int argc, char** argv)
{

double beta0 = 1.5;
double gamma0 =0.2;

double I00 = 1e-6;
double S00 =1-I00;

double tmax0 = 50;
/* Find a suitable time-scale for outputs */

double step0=0.01/((beta0+gamma0)*S00);

SIR mySIR(beta0, gamma0,step0,S00,  I00,  tmax0);

mySIR.Solve_Eq();
return(0);

}

2) Matlab
in SIR.m
--------------------------------------------------------
function [t,S,I,R] =SIR()
%Solve SIR equation in Matlab

beta=1.5; gamma=0.2; I0=1e-6;
S0=1-I0; tmax=50;
S=S0; I=I0; R=1-S-I;

% The main iteration
[t, pop]=ode45(@Diff_SIR,[0 tmax],[S I R],[],[beta gamma]);

S=pop(:,1); I=pop(:,2); R=pop(:,3);

% plots the graphs with scaled colours
subplot(2,1,1)
h=plot(t,S,'-g',t,R,'-k');
xlabel 'Time';
ylabel 'Susceptibles and Recovereds'

subplot(2,1,2)
h=plot(t,I,'-r');
xlabel 'Time';
ylabel 'Infectious'

% Calculates the differential rates used in the integration.
function dPop=Diff_SIR(t,pop, parameter)
beta=parameter(1); gamma=parameter(2);
S=pop(1); I=pop(2); R=pop(3);
dPop=zeros(3,1);
dPop(1)= -beta*S*I;
dPop(2)= beta*S*I - gamma*I;
dPop(3)= gamma*I;
-------------------------------------------------
3. Python
in SIR.py (chmod  +x SIR.py, ./SIR.py)
------------------------------------------------
#!/usr/bin/env python



import scipy.integrate as spi
import numpy as np
import pylab as pl

beta=1.5
gamma=0.2
TS=1.0
tmax =50
I0=1e-6
S0=1-I0
INPUT = (S0, I0, 0.0)

def diff_eqs(INP,t):
'''The main set of equations'''
Y=np.zeros((3))
V = INP
Y[0] = - beta * V[0] * V[1]
Y[1] = beta * V[0] * V[1] - gamma * V[1]
Y[2] = gamma * V[1]
return Y   # For odeint

t_start = 0.0; t_end = tmax; t_inc = TS
t_range = np.arange(t_start, t_end+t_inc, t_inc)
RES = spi.odeint(diff_eqs,INPUT,t_range)

print RES

#Ploting
pl.subplot(211)
pl.plot(RES[:,0], '-g', label='Susceptibles')
pl.plot(RES[:,2], '-k', label='Recovereds')
pl.legend(loc=0)
pl.title('SIR.py')
pl.xlabel('Time')
pl.ylabel('Susceptibles and Recovereds')
pl.subplot(212)
pl.plot(RES[:,1], '-r', label='Infectious')
pl.xlabel('Time')
pl.ylabel('Infectious')
pl.show()