Cholera-Modeling/SIR.m at master · waldmannly/Cholera-Modeling · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
%% SIR MODEL from MAA
%%
clc

% type RKStage5.m

format compact
format long

% estimation with stage 5 Runge Kutta method

% s= y(1) = susectible
% i= y(2) = infection
% r= y(2) = immune

initialPop = 10900000; %10.9 million
initialLevelofInfection = 8400;
initialImmune = 0 ;
b= 1/3; %infecting contact period (one contact every 3 days)
k= 1/7; % average period of infectiousness is about 7 days
sknot =  (initialPop)/ (initialPop + initialImmune);
iknot = initialLevelofInfection/(initialPop+initialImmune) ;
rknot = initialImmune/ (initialPop + initialImmune);
time = 150;


timex = [0:1:time];

out = RKStage5( @(Y) [-b*Y(1)*Y(2), b*Y(1)*Y(2)-k*Y(2), k*Y(2)], [sknot, iknot, rknot], 1, time, 1);

totalPopArr = ones(time+1 , 1);
percentofpeoplethatdie = .60;
totalPopArr = totalPopArr - percentofpeoplethatdie*out(:,3) ;

hold on
a1 = plot(timex,out(:,1),'m');
M1 = 'susceptible fraction of pop';

a2 = plot(timex, out(:,2), 'r');
M2 = 'infected fraction of pop   ';

a3 = plot(timex, percentofpeoplethatdie*out(:,3), 'b');
M3 = 'dead fraction of pop       ';

a4 = plot(timex, totalPopArr, 'g');
M4 = 'total fraction of pop left ';

a5 = plot(timex, (1-percentofpeoplethatdie)*out(:,3), 'y');
M5 = 'immune/recovered fraction  ';

xlabel('Time in Days'),ylabel('Fraction of People'), title('SIR MODEL with scaled variables from Haiti ( 60% deathrate )')
legend([a1; a2; a3; a4; a5], [M1; M2; M3; M4; M5]);
snapnow
hold off
clf


initialPop = 10900000; %10.9 million
initialLevelofInfection = 8400;
initialImmune = 0 ;
b= 1/3; %infecting contact period (one contact every 3 days)
k= 1/7; % average period of infectiousness is about 7 days
sknot =  (initialPop)/ (initialPop + initialImmune);
iknot = initialLevelofInfection/(initialPop+initialImmune) ;
rknot = initialImmune/ (initialPop + initialImmune);
time = 150;


timex = [0:1:time];

out = RKStage5( @(Y) [-b*Y(1)*Y(2), b*Y(1)*Y(2)-k*Y(2), k*Y(2)], [sknot, iknot, rknot], 1, time, 1);

totalPopArr = ones(time+1 , 1);
percentofpeoplethatdie = .30;
totalPopArr = totalPopArr - percentofpeoplethatdie*out(:,3) ;

hold on
a1 = plot(timex,out(:,1),'m');
M1 = 'susceptible fraction of pop';

a2 = plot(timex, out(:,2), 'r');
M2 = 'infected fraction of pop   ';

a3 = plot(timex, percentofpeoplethatdie*out(:,3), 'b');
M3 = 'dead fraction of pop       ';

a4 = plot(timex, totalPopArr, 'g');
M4 = 'total fraction of pop left ';

a5 = plot(timex, (1-percentofpeoplethatdie)*out(:,3), 'y');
M5 = 'immune/recovered fraction  ';

xlabel('Time in Days'),ylabel('Fraction of People'), title('SIR MODEL with scaled variables from Haiti ( 30% deathrate )')
legend([a1; a2; a3; a4; a5], [M1; M2; M3; M4; M5]);
snapnow