Deterministic SEIR ODE model running slow
16 次查看(过去 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 个评论
采纳的回答
Torsten
2022-8-30
编辑:Torsten
2022-8-30
%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 个评论
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.
更多回答(1 个)
Steven Lord
2022-8-30
The first thing I'd probably try is to use a stiff ODE solver instead of the nonstiff ODE solver ode45.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!