Deterministic SEIR ODE model running slow

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?

Answer 1

[Torsten]

%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

[Gbeminiyi Oyedele]

Thank you this helps.

Though I might have to check if it's not my device. I try to run it up to 100 iteration because I need to get some values for the infectious class, but it take too much time. That is, beyong 70 iteration it takes very very long to run. I don't know why this is like that

[Torsten]

That is, beyong 70 iteration it takes very very long to run. I don't know why this is like that

As Steven Lord already suggested, switch to ode15s instead of ode45.

The results will come up within a second.

This indicates that your system is stiff.

I thought you already tried it without success - that's why I didn't suggest it.

I incorporated the necessary changes in the code above.

[Gbeminiyi Oyedele]

Yes, thank you. Runs better now.

Answer 2

[Steven Lord]

The first thing I'd probably try is to use a stiff ODE solver instead of the nonstiff ODE solver ode45.

[Gbeminiyi Oyedele]

Thank you for your response. But It still take longer time like ODE45.