-
Notifications
You must be signed in to change notification settings - Fork 0
/
choleraData.asv
118 lines (112 loc) · 4.16 KB
/
choleraData.asv
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
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
function data=choleraData(high, mid, low)
time = 1000; %time in days
timex = [0:1:time]; %time array
Narr= [high;mid;low]; %population
initS = [Narr(1)*.999;Narr(2)*.999;Narr(3)*.999]; %initial susecptible pop
initI = [Narr(1)*.001;Narr(2)*.001;Narr(3)*.001]; %initial infected pop
initR = [Narr(1)*0;Narr(2)*0;Narr(3)*0]; %initial recovered pop
initB = [0;0;0]; %initial bateria in water
N = sum(Narr) ;
Lambda = [.0000638356;.0000638356;.0000638356]; %birth rates
mu= [.0000210958;.0000210958;.0000210958]; %mortality rate
beta = 10^-8*[1 1 1;
1 1 1; %human to human transmission rate
1 1 1];
epsilon = [.0002739726;.0002739726;.0002739726]; %loss of immunity rate
lambda = 10^-8* [1;1;1]; %water to human transmission rate
gamma =[.2;.2;.2]; %recovery rate
alpha =[.0005;.0005;.0005]; %disease induce mortality
xi = [7;7;7]; %shedding rate
delta =[.33;.33;.33]; %baterial death rate
syms
flist = @(Y) [+ Lambda(1)*Narr(1)-mu(1)*Y(1)-beta(1,1)*Y(1)*Y(2)-beta(1,2)*Y(1)*Y(6)-lambda(1)*Y(1)*Y(4)+epsilon(1)*Y(3),...
beta(1,1)*Y(1)*Y(2)+beta(1,2)*Y(1)*Y(6)+lambda(1)*Y(1)*Y(4)-(gamma(1)+mu(1)+alpha(1))*Y(2),...
gamma(1)*Y(2)-mu(1)*Y(3)-epsilon(1)*Y(3),...
xi(1)*Y(2)-delta(1)*Y(4),...
Lambda(2)*Narr(2)-mu(2)*Y(5)-beta(2,2)*Y(5)*Y(6)-beta(2,1)*Y(5)*Y(2)-beta(2,3)*Y(5)*Y(10)-lambda(2)*Y(5)*Y(8)+epsilon(2)*Y(7),...
beta(2,2)*Y(5)*Y(6)+beta(2,1)*Y(5)*Y(2)+beta(2,3)*Y(5)*Y(10)+lambda(2)*Y(5)*Y(8)-(gamma(2)+mu(2)+alpha(2))*Y(6),... %6Im
gamma(2)*Y(6)-mu(2)*Y(7)-epsilon(2)*Y(7),...
xi(2)*Y(6)-delta(2)*Y(8),...
Lambda(3)*Narr(3)-mu(3)*Y(9)-beta(3,3)*Y(9)*Y(10)-beta(3,2)*Y(9)*Y(6)-lambda(3)*Y(9)*Y(12)+epsilon(3)*Y(11),...
beta(3,3)*Y(9)*Y(10)+beta(3,2)*Y(9)*Y(6)+lambda(3)*Y(9)*Y(12)-(gamma(3)+mu(3)+alpha(3))*Y(10),...
gamma(3)*Y(10)-mu(3)*Y(11)-epsilon(3)*Y(11),...
xi(3)*Y(10)-delta(3)*Y(12)];
init = [initS(1) initI(1) initR(1) initB(1) initS(2) initI(2) initR(2) initB(2) initS(3) initI(3) initR(3) initB(3)];
h=1;
steps=time;
freq=1;
out = RKStage5( flist ,init, 1, time, 1);
totS = out(:,1)+out(:,5)+out(:,9);
totI = out(:,2)+out(:,6)+out(:,10);
totR = out(:,3)+out(:,7)+out(:,11);
totB = out(:,4)+out(:,8)+out(:,12);
totals = [totS totI totR];
endTotal = totals(time,:, :, :)
totals = [totS totI totR, totB];
hold on
a1 = plot(timex,totS,'m--');
M1 = 'susceptible pop';
a2 = plot(timex,totI, 'k-.');
M2 = 'infected pop ';
a3 = plot(timex,totR , 'b');
M3 = 'recovered pop ';
xlabel('Time in Days'),ylabel('Total People'), title('SIWR GDP STRUCTURED MODEL total')
legend([a1; a2; a3], [M1; M2; M3]);
snapnow
hold off
% clf
% hold on
% a1 = plot(timex,out(:,1),'m--');
% M1 = 'susceptible pop';
% a2 = plot(timex,out(:,2), 'k-.');
% M2 = 'infected pop ';
% a3 = plot(timex,out(:,3) , 'b');
% M3 = 'recovered pop ';
% xlabel('Time in Days'),ylabel('Total People'), title('SIWR GDP STRUCTURED MODEL high')
% legend([a1; a2; a3], [M1; M2; M3]);
% snapnow
% hold off
% clf
% hold on
% a1 = plot(timex,out(:,5),'m--');
% M1 = 'susceptible pop';
% a2 = plot(timex,out(:,6), 'k-.');
% M2 = 'infected pop ';
%
% a3 = plot(timex,out(:,7) , 'b');
% M3 = 'recovered pop ';
% xlabel('Time in Days'),ylabel('Total People'), title('SIWR GDP STRUCTURED MODEL mid')
% legend([a1; a2; a3], [M1; M2; M3]);
% snapnow
% hold off
% clf
% hold on
% a1 = plot(timex,out(:,9),'m--');
% M1 = 'susceptible pop';
% a2 = plot(timex,out(:,10), 'k-.');
% M2 = 'infected pop ';
% a3 = plot(timex,out(:,11) , 'b');
% M3 = 'recovered pop ';
% xlabel('Time in Days'),ylabel('Total People'), title('SIWR GDP STRUCTURED MODEL low')
% legend([a1; a2; a3], [M1; M2; M3]);
% snapnow
% hold off
% clf
% hold on
% a1 = plot(timex,out(:,4),'m--');
% M1 = 'bacteria high ';
% a2 = plot(timex,out(:,8), 'k-.');
% M2 = 'bateria mid ';
% a3 = plot(timex,out(:,12) , 'b');
% M3 = 'bateria low ';
% xlabel('Time in Days'),ylabel('Total bateria'), title('bateria differences')
% legend([a1; a2; a3], [M1; M2; M3]);
% snapnow
% hold off
peaksinfected = [max(out(:,2)); max(out(:,6)); max(out(:,10))]
peakTotInfected = max(totI)
peaksbacteria = [max(out(:,4)); max(out(:,8)); max(out(:,12))]
peakTotBateria = max(totB)
sumtotals = sum(totals,1)
data = out;
end