Deterministic SEIR ODE model running slow

Question

Gbeminiyi Oyedele 2022-8-30

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1790880-deterministic-seir-ode-model-running-slow

评论： Gbeminiyi Oyedele 2022-8-30

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2022-8-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1790880-deterministic-seir-ode-model-running-slow#answer_1038135

编辑：Torsten 2022-8-30

在 MATLAB Online 中打开

%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 个评论
显示 1更早的评论隐藏 1更早的评论

Torsten 2022-8-30

编辑：Torsten 2022-8-30

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 2022-8-30

Yes, thank you. Runs better now.

请先登录，再进行评论。

Answer 2

Steven Lord 2022-8-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1790880-deterministic-seir-ode-model-running-slow#answer_1037965

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

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Gbeminiyi Oyedele 2022-8-30

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

请先登录，再进行评论。