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()
Online computer courses, code, programming tutorial and sidebar information for monitoring Canadian S&P/TSX index. Build friendship and networking. Welcome to visit my blogs often!!! I also have two other sites: YouTube Channel and Google site.
Adsense
Popular Posts
- PHPWind-- A PHP forum script applcaition in China
- How to blend adsense inside your post?
- Formatting my post
- Install PHPMailer 5.2.4 and use smtp gmail
- Datatable export excel wraptext and newline
- phpexcel toggle expand and hide column in EXCEL and summary
- Linux, automatically backup MySQL database daily
- SPSS job interview questions
- Notepad++ - Add C++ compiler
- Install PHP ibm_db2 extension in Linux (redHat)
Friday, June 25, 2010
Solve SIR model -C++/MatLab/Python
Subscribe to:
Post Comments (Atom)
good day. Might I explain in detail the code that you generated in matlab for the SIR model or function that meets this term "pop" and "DPOP"
ReplyDeleteHello There,
ReplyDeleteGrazie! Grazie! Grazie! Your blog is indeed quite interesting around Jiansen Lu's Computing Blog with you on lot of points!
A subset of a natural language contains these sentences:
1) Triangle ABC.
2) Line Segment AB.
3) Angle A.
In this set, names of triangles are formed by putting three letters together,
all drawn from the alphabet {A,B,C,D,E}. Line segment names are defined by
putting two letters together from the previous alphabet. And angle names are
(suprise!) given by any letter in that alphabet.
Great effort, I wish I saw it earlier. Would have saved my day :)
Thank you,
Uday
Python Training, SVR Technologies offers Online & Corporate Classes with Free Live Demo, We cover practical sessions and all the modules in Python.
ReplyDeletepython training course