Deterministic SEIR ODE model running slow

2 views (last 30 days)
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
function [Classes] = ODE_SEIR_adhe(para,ICs,maxtime)
%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
function dPop = diff_SEIR_adhe(t,pop,para)
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);
%The Non-adherent population
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;
%The Adherent population
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
[Classes] = ODE_SEIR_adhe(para,ICs,maxtime);
Is there something I am doing wrongly?

Accepted Answer

Torsten
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
[Classes] = ODE_SEIR_adhe(para,ICs,maxtime);
plot(Classes.t,Classes.S)
%ODE SI model code
function [Classes] = ODE_SEIR_adhe(para,ICs,maxtime)
%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
function dPop = diff_SEIR_adhe(t,pop,para)
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);
%The Non-adherent population
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;
%The Adherent population
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
  3 Comments

Sign in to comment.

More Answers (1)

Steven Lord
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.
  1 Comment
Gbeminiyi Oyedele
Gbeminiyi Oyedele on 30 Aug 2022
Thank you for your response. But It still take longer time like ODE45.

Sign in to comment.

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by