Changing initial condition for ODE
1 view (last 30 days)
Show older comments
I am simulating a disease model using ODE45. What I want to achieve is that: as infection grows, I want to implement some lockdown measures that will reduce infection growth and relax the measure when infection growth as reduced to a particular value.
The system of ODE I am using is below. Which is fine.
%ODE SI model code
function [Classes] = ODE_SIR(para,ICs,maxtime)
%Run ODE using ODE45
opts = odeset('RelTol',1e-5);
[t, pop] = ode45(@diff_SIR_open_model, [1:maxtime], [ICs.S ICs.I ICs.R ICs.Sa ICs.Ia ICs.Ra], opts, para);
%Convert output to structure
Classes = struct('S',pop(:,1),'I',pop(:,2),'R',pop(:,3),'Sa',pop(:,4),'Ia',pop(:,5),'Ra',pop(:,6),'t',t);
%Diff equations
function dPop = diff_SIR_open_model(t,pop,para)
S=pop(1);
I=pop(2);
R=pop(3);
Sa=pop(4);
Ia=pop(5);
Ra=pop(6);
%The Non-adherent population
dS = - (para.beta0(1,1)*S*I + para.beta0(1,2)*S*Ia)/para.N;
dI = (para.beta0(1,1)*S*I + para.beta0(1,2)*S*Ia)/para.N -para.gamma*I ;
dR = para.gamma*I;
%The Adherent population
dSa = - (para.beta0(2,1)*Sa*I + para.beta0(2,2)*Sa*Ia)/para.Na;
dIa = (para.beta0(2,1)*Sa*I + para.beta0(2,2)*Sa*Ia)/para.Na -para.gamma*Ia;
dRa = para.gamma*Ia;
dPop = [dS; dI; dR; dSa; dIa; dRa];
end
end
However, I created a function to reduce Beta such that it impact the steady growth of the infection dynamic and reduce infection.
%% Create a function that changes rate on transmission base on NPI measure
% Input :
% Lockdown : The strenght of lockdown
% School Closure : Either school is opened or closed during the period in
% question
% Beta: The initial value of Beta when there is no control
% Output:
% Beta : The new value of Beta after control is implemented
function [new_Beta] = npi(Lockdown, School_closure,Beta)
ReductionLockdown = [1.0]; %Lockdown values
ReductionSchool = [1.0]; %School closure value
if Lockdown
ReductionLockdown = 1-Lockdown;
end
if School_closure
ReductionSchool = 1-School_closure;
end
new_Beta = Beta * (ReductionLockdown * ReductionSchool)*10;
Which is also working fine.
However when I tried to output my result, I am experiencing a sharp change in my output meanwhile I expect to see a sinuosodial kind of curve (That is the values will should not get to zero and pick back up again that is wrong). This is basically because I want to take the last value of the previous output as the next Initial condition.
How am i outputing this wrongly.
function [I, Ia,t] = npichange(para,ICs,maxtime,Beta)
%% Run model
[Classes] =ODE_SIR(para,ICs,maxtime);
for i = 1:maxtime
if (Classes.I(end) >0) && (Classes.I(end) <=500)
[new_Beta] = npi(0.5, 0.5,Beta);
para.beta0 = new_Beta;
ICs = struct('S',Classes.S(end),'I',Classes.I(end), ...
'R',Classes.R(end),'Sa',Classes.Sa(end),'Ia',Classes.Ia(end), ...
'Ra',Classes.Ra(end));
[Classes] =ODE_SIR(para,ICs,maxtime);
end
if (Classes.I(end) >500) && (Classes.I(end) <=15000)
[new_Beta] = npi(0.75, 0.75,Beta);
para.beta0 = new_Beta;
ICs = struct('S',Classes.S(end),'I',Classes.I(end), ...
'R',Classes.R(end),'Sa',Classes.Sa(end),'Ia',Classes.Ia(end), ...
'Ra',Classes.Ra(end));
[Classes] =ODE_SIR(para,ICs,maxtime);
end
end
t=[Classes.t; Classes.t(end)+Classes.t(2:end)];
Classes.S=[Classes.S; Classes.S(2:end,:)]; Classes.Sa=[Classes.Sa; Classes.Sa(2:end,:)];
I=[Classes.I; Classes.I(2:end,:)]; Ia=[Classes.Ia; Classes.Ia(2:end,:)];
Classes.R=[Classes.R; Classes.R(2:end,:)]; Classes.Ra=[Classes.Ra; Classes.Ra(2:end,:)];
The values I am using to run the model
%Define model parameters as a structure (N.B. stick to days here for ease
%later on)
Beta = [0.5 0.1; 0.001 0.5];
para = struct('beta0',Beta,'N',0.2*62000000,'gamma',1/9, 'Na',0.8*62000000);
%Define initial conditions as a structure
ICs = struct('S',para.N-1,'I',1,'R',0,'Sa',para.Na,'Ia',0,'Ra',0);
%Define time to run model for
maxtime = 80;
[I, Ia,t] = npichange(para,ICs,maxtime,Beta);
%% Plots
figure(1)
clf
set(gca,'fontsize',16)
plot(t,I,'r',t,Ia,'k')
xlabel('Time (days)')
ylabel('Infection')
legend('Non-adherent Population', 'Adherent Population')
0 Comments
Answers (1)
SANKALP DEV
on 21 Sep 2023
Hi Gbeiniyi,
I understand that you are getting abrupt changes in your output.
This occurs because the values of the ‘Classes’ variable is being overwritten at each iteration of the loop. As a result, the previous values are lost, and only the values from the current iteration are retained.
To address the problem, you can store the results of each iteration separately and concatenate them at the end.
This can be done by introducing a helper function which takes the initial results and the results of each iteration, and concatenates them field by field, preserving the complete set of results.
To understand more, refer to following code:
% Helper function to concatenate the results
function results = concatenateResults(results1,results2)
fields = fieldnames(results1);
for i = 1:numel(fields)
field = fields{i};
results1.(field) = [results1.(field); results2.(field)(2:end)];
end
results = results1;
end
%You can call the above helper function after the variable %Classes’ is being updated in the following manner:
[Classes] = ODE_SIR(para, ICs, maxtime);
results = concatenateResults(results, Classes);
Hope this helps you achieve the desired output!
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!