Deterministic SEIR ODE model running slow

2 views (last 30 days)
Gbeminiyi Oyedele on 30 Aug 2022
Commented: Gbeminiyi Oyedele on 30 Aug 2022
Good Afternoon,
I am running an ODE extended SEIR risk based model and it is running exceptionally slow. I am runnig 100 time steps and I expect it to be done in less then 5seconds but it is taking over 2 hours without running.
This is the ODE code:
%ODE SI model code
%Run ODE using ODE45
opts = odeset('RelTol',1e-2);
[t, pop] = ode45(@diff_SEIR_adhe, [0:maxtime], [ICs.S ICs.E ICs.I ICs.A ICs.H ICs.R ICs.Sa ICs.Ea ICs.Ia ICs.Aa ICs.Ha ICs.Ra],opts, para);
%Convert output to structure
Classes = struct('S',pop(:,1),'E',pop(:,2),'I',pop(:,3),'A',pop(:,4),'H',pop(:,5),'R',pop(:,6), ...
'Sa',pop(:,7),'Ea',pop(:,8),'Ia',pop(:,9),'Aa',pop(:,10),'Ha',pop(:,11),'Ra',pop(:,12),'t',t);
%Diff equations
S=pop(1);
E=pop(2);
I=pop(3);
A=pop(4);
H=pop(5);
R=pop(6);
Sa=pop(7);
Ea=pop(8);
Ia=pop(9);
Aa=pop(10);
Ha=pop(11);
Ra=pop(12);
dS = - ((para.beta0(1,1))*S*I + (para.beta0(1,2))*S*Ia)/para.Np-para.gamma*S+para.gamma_a*Sa;
dE = ((para.beta0(1,1))*S*I + (para.beta0(1,2))*S*Ia)/para.Np- para.sigma*E-para.gamma*E+para.gamma_a*Ea;
dI = para.omega*para.sigma*E - para.mu_0*I -para.eta*I - para.gamma*I +para.gamma_a*Ia ;
dA = (1-para.omega)*para.sigma*E - para.mu_1*A - para.gamma*A +para.gamma_a*Aa;
dH = para.eta*I -para.mu_2*H - para.phi*H- para.gamma*H + para.gamma_a*Ha;
dR = para.mu_1*A + para.mu_0*I + para.mu_2*H - para.gamma*R + para.gamma_a*Ra;
dSa = - ((para.beta0(2,1))*Sa*I + ((para.beta0(2,2))*Sa*Ia))/para.Np +para.gamma*S-para.gamma_a*Sa;
dEa = ((para.beta0(2,1))*Sa*I + ((para.beta0(2,2))*Sa*Ia))/para.Np +para.sigma*Ea +para.gamma*E-para.gamma_a*Ea;
dIa = para.omega*para.sigma*Ea - para.mu_0*Ia -para.eta*Ia + para.gamma*I -para.gamma_a*Ia ;
dAa = (1-para.omega)*para.sigma*Ea - para.mu_1*Aa + para.gamma*A -para.gamma_a*Aa;
dHa = para.eta*Ia -para.mu_2*Ha - para.phi*Ha+ para.gamma*H - para.gamma_a*Ha;
dRa = para.mu_1*Aa + para.mu_0*Ia + para.mu_2*Ha + para.gamma*R - para.gamma_a*Ra;
dPop = [dS; dE; dI; dA; dH; dR; dSa; dEa; dIa; dAa; dHa; dRa];
end
end
Then I am runing it with the following initial conditions and parameter values:
%Define model parameters as a structure (N.B. stick to days here for ease
%later on)
A = [10 0.1; 0.1 1];
J = 13000000+53000000;
para = struct('beta0',A,'N',13000000,'Na',53000000,'sigma',0.5,'phi',0.0018, ...
'mu_0',1/20,'mu_1',1/30,'mu_2',0.0714,'omega',0.7,'gamma',0.3,'gamma_a',0.25, ...
'eta',0.2,'l',0,'Np',J);
%Define initial conditions as a structure
ICs = struct('S',para.N-1,'E',1,'I',1,'A',0,'H',0,'R',0, ...
'Sa',para.Na-0.001,'Ea',0.001,'Ia',0,'Aa',0,'Ha',0,'Ra',0);
%Define time to run model for
maxtime = 100;
%Run model
Is there something I am doing wrongly?

Torsten on 30 Aug 2022
Edited: Torsten on 30 Aug 2022
%Define model parameters as a structure (N.B. stick to days here for ease
%later on)
A = [10 0.1; 0.1 1];
J = 13000000+53000000;
para = struct('beta0',A,'N',13000000,'Na',53000000,'sigma',0.5,'phi',0.0018, ...
'mu_0',1/20,'mu_1',1/30,'mu_2',0.0714,'omega',0.7,'gamma',0.3,'gamma_a',0.25, ...
'eta',0.2,'l',0,'Np',J);
%Define initial conditions as a structure
ICs = struct('S',para.N-1,'E',1,'I',1,'A',0,'H',0,'R',0, ...
'Sa',para.Na-0.001,'Ea',0.001,'Ia',0,'Aa',0,'Ha',0,'Ra',0);
%Define time to run model for
maxtime = 100;
%Run model
plot(Classes.t,Classes.S)
%ODE SI model code
%Run ODE using ODE45
%opts = odeset('RelTol',1e-2);
tspan = linspace(0,maxtime,100);
[t, pop] = ode15s(@(t,y)diff_SEIR_adhe(t,y,para), tspan, [ICs.S ICs.E ICs.I ICs.A ICs.H ICs.R ICs.Sa ICs.Ea ICs.Ia ICs.Aa ICs.Ha ICs.Ra]);%,opts);
%Convert output to structure
Classes = struct('S',pop(:,1),'E',pop(:,2),'I',pop(:,3),'A',pop(:,4),'H',pop(:,5),'R',pop(:,6), ...
'Sa',pop(:,7),'Ea',pop(:,8),'Ia',pop(:,9),'Aa',pop(:,10),'Ha',pop(:,11),'Ra',pop(:,12),'t',t);
end
%Diff equations
S=pop(1);
E=pop(2);
I=pop(3);
A=pop(4);
H=pop(5);
R=pop(6);
Sa=pop(7);
Ea=pop(8);
Ia=pop(9);
Aa=pop(10);
Ha=pop(11);
Ra=pop(12);
dS = - ((para.beta0(1,1))*S*I + (para.beta0(1,2))*S*Ia)/para.Np-para.gamma*S+para.gamma_a*Sa;
dE = ((para.beta0(1,1))*S*I + (para.beta0(1,2))*S*Ia)/para.Np- para.sigma*E-para.gamma*E+para.gamma_a*Ea;
dI = para.omega*para.sigma*E - para.mu_0*I -para.eta*I - para.gamma*I +para.gamma_a*Ia ;
dA = (1-para.omega)*para.sigma*E - para.mu_1*A - para.gamma*A +para.gamma_a*Aa;
dH = para.eta*I -para.mu_2*H - para.phi*H- para.gamma*H + para.gamma_a*Ha;
dR = para.mu_1*A + para.mu_0*I + para.mu_2*H - para.gamma*R + para.gamma_a*Ra;
dSa = - ((para.beta0(2,1))*Sa*I + ((para.beta0(2,2))*Sa*Ia))/para.Np +para.gamma*S-para.gamma_a*Sa;
dEa = ((para.beta0(2,1))*Sa*I + ((para.beta0(2,2))*Sa*Ia))/para.Np +para.sigma*Ea +para.gamma*E-para.gamma_a*Ea;
dIa = para.omega*para.sigma*Ea - para.mu_0*Ia -para.eta*Ia + para.gamma*I -para.gamma_a*Ia ;
dAa = (1-para.omega)*para.sigma*Ea - para.mu_1*Aa + para.gamma*A -para.gamma_a*Aa;
dHa = para.eta*Ia -para.mu_2*Ha - para.phi*Ha+ para.gamma*H - para.gamma_a*Ha;
dRa = para.mu_1*Aa + para.mu_0*Ia + para.mu_2*Ha + para.gamma*R - para.gamma_a*Ra;
dPop = [dS; dE; dI; dA; dH; dR; dSa; dEa; dIa; dAa; dHa; dRa];
end
Gbeminiyi Oyedele on 30 Aug 2022
Yes, thank you. Runs better now.

Steven Lord on 30 Aug 2022
The first thing I'd probably try is to use a stiff ODE solver instead of the nonstiff ODE solver ode45.
Gbeminiyi Oyedele on 30 Aug 2022
Thank you for your response. But It still take longer time like ODE45.

Categories

Find more on Ordinary Differential Equations in Help Center and File Exchange

R2021b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by