Add quarantine to the model

SRI_matlab/main.m

@@ -4,6 +4,7 @@ N = 4*365;
 
 S = 410000000;
 I = 55595;
+Q = 0;
 R = 0;
 
 coeffs = options;
@@ -13,10 +14,11 @@ coeffs.ImmunisationLossRate = 2.7e-4;
 coeffs.DeathRate = 1.3e-2;
 coeffs.NaturalDeathRate = 2e-4;
 coeffs.NatalityRate = 5.3e-4;
+coeffs.QuarantineRate = 4.4e-4;
 VaccinationVector = build_vacc_rate(0, 4.4e-4, 0, N);
 % VaccinationVector = build_vacc_rate(0, 0, 50, N);
 
-m = model(S, I, R, 12, N);
+m = model(S, I, Q, R, 12, N);
 
 t = 1:N;
 
@@ -30,9 +32,9 @@ end
 % plot(t, Ilist);
 % plot(t, Rlist);
 % plot(t, Dlist);
-area([m.Slist.', m.Wlist.', m.Ilist.', m.Rlist.']);
+area([m.Slist.', m.Wlist.', m.Ilist.', m.Qlist.', m.Rlist.']);
 
-legend('Susceptible', 'Waiting', 'Infected', 'Immunized');
+legend('Susceptible', 'Waiting', 'Infected', 'In quarantine', 'Immunized');
 title('Evolution of virus in population');
 xlabel('Time');
 ylabel('Population');

SRI_matlab/model.m

@@ -2,43 +2,43 @@ classdef model < handle
     properties
         vector
         Slist
+        Wlist
         Ilist
+        Qlist
         Rlist
-        Wlist
     end
     
     methods
-        function obj=model(S, I, R, Waiting, N)
-            if nargin > 4
+        function obj=model(S, I, Q, R, Waiting, N)
+            if nargin > 5
                 repeat = N;
             else
                 repeat = 100;
             end
             
-            obj.vector = zeros(1, 3 + Waiting);
+            obj.vector = zeros(1, 4 + Waiting);
             obj.vector(1) = S;
             obj.vector(2) = I;
-            obj.vector(3) = R;
+            obj.vector(3) = Q;
+            obj.vector(4) = R;
             for i = 1:Waiting
-                obj.vector(3+i) = 0;
+                obj.vector(4+i) = 0;
             end
             obj.Slist = zeros(1, repeat);
+            obj.Wlist = zeros(1, repeat);
             obj.Ilist = zeros(1, repeat);
+            obj.Qlist = zeros(1, repeat);
             obj.Rlist = zeros(1, repeat);
-            obj.Wlist = zeros(1, repeat);
         end
         
         function step(obj, i, coeffs, VaccinationVec)
             obj.Slist(i) = obj.vector(1);
             obj.Ilist(i) = obj.vector(2);
-            obj.Rlist(i) = obj.vector(3);
-            obj.Wlist(i) = sum(obj.vector(4:length(obj.vector)));
+            obj.Qlist(i) = obj.vector(3);
+            obj.Rlist(i) = obj.vector(4);
+            obj.Wlist(i) = sum(obj.vector(5:length(obj.vector)));
             obj.vector = step(obj.vector, coeffs, VaccinationVec(i));
         end
-        
-        function goto(obj, other)
-            
-        end
     end
 end
     
\ No newline at end of file

SRI_matlab/options.m

@@ -6,6 +6,7 @@ classdef options
         DeathRate
         NaturalDeathRate
         NatalityRate
+        QuarantineRate
     end
 end
 

SRI_matlab/step.m

 function res = step(X, options, VaccinationRate)
 S = X(1);
 I = X(2);
-R = X(3);
+Q = X(3);
+R = X(4);
 
 N = sum(X);
 
-new_waiting = zeros(1, length(X) - 3);
-for i=5:length(X)
-    new_waiting(i - 3) = X(i - 1) * (1 - options.NaturalDeathRate);
+new_waiting = zeros(1, length(X) - 4);
+for i=6:length(X)
+    new_waiting(i - 4) = X(i - 1) * (1 - options.NaturalDeathRate);
 end
 
-if length(X) > 3
+if length(X) > 4
     % infection = S * (1 - 25 * options.InfectionRate / N) ^ I;
-    infecter = sum(X(4:length(X)) + I);
+    infecter = sum(X(5:length(X)) + I);
     infection = options.InfectionRate * S * infecter/N;
     new_waiting(1) = infection;
     dS = -infection;
@@ -23,8 +24,9 @@ else
 end
 
 dS = dS - VaccinationRate * S + options.ImmunisationLossRate * R + options.NatalityRate * N - options.NaturalDeathRate * S;
-dI = dI - options.ImmunisationRate * I - options.DeathRate * I;
-dR = options.ImmunisationRate * I + VaccinationRate * S - options.ImmunisationLossRate * R - options.NaturalDeathRate * R;
+dI = dI - options.ImmunisationRate * I - options.DeathRate * I - options.QuarantineRate * I;
+dQ = -options.ImmunisationRate * Q - options.DeathRate * Q + options.QuarantineRate * I;
+dR = options.ImmunisationRate * (I + Q) + VaccinationRate * S - options.ImmunisationLossRate * R - options.NaturalDeathRate * R;
 
-res = [max(S + dS, 0), max(I + dI, 0), max(R + dR, 0), new_waiting];
+res = [max(S + dS, 0), max(I + dI, 0), max(Q + dQ, 0), max(R + dR, 0), new_waiting];
 end
\ No newline at end of file