计算科学 - 如何停止负数和错误消息：“ t=3.562559e+03 失败。无法满足积分容差 - 吾爱随笔录

我在 Matlab 中使用 8 个 ODE 来模拟无症状感染在媒介传播疾病流行病学中的影响。在某些设置下搜索参数空间会为人口产生负数，并在命令 consol 中出现以下警告：

“警告：在 t=3.562559e+03 时失败。如果不将步长减小到低于时间 t 允许的最小值 (7.275958e-12)，则无法满足积分容差。> 在 ode45 中 308 在 Mosquito_Framework_2_human_plots 中 79”

是否有一些设置我可以使用求解器来解决问题。当前的求解器是“RelTol”下的 ode45，1e-6 选项。

以下代码产生负面人口和警告消息：

function [t,Sh,Ah,Ih,Rh,Se,Sv,Iv] = Mosquito_Framework_2_human_plots(betaI,betaV,c,...
theta,d,muh,p,muv,omega,epsilon,Sh0,Ah0,Ih0,Rh0,Sv0,Ev0,Iv0,MaxTime)

% Sets up default parameters if necessary.
if nargin ~= 18
    disp('defaults')
    %displays message default to say simulation is using defaults if not
    %enouth argumants is given
   betaI = (0.3+1)/2*(0.1+0.75)/2;
   betaV = (0.3+1)/2*(0.5+1)/2;
   c=1;
   theta=8;
   d=0.25;
   muh=1/(60*365.0);
   p=0.9;
   muv=1/6;
   omega=1/11;
   epsilon=1/((8+12)/2);
   Sh0=(1e6)-1;
   Ah0=1;
   Ih0=0;
   Rh0=0;
   Sv0=9.5e6;
   Ev0=0;
   Iv0=0;
   MaxTime=10*(365);
end

gammaI=1/theta;
gammaA=1/(d*theta);
betaA=c*betaI;
Nh0=Sh0+Ah0+Ih0+Rh0;
Nv0=Sv0+Ev0+Iv0; 
Se0=Nv0*(muv/omega);

 % Checks all the parameters are valid
 if Sh0<=0 
     error('Initial level of susceptibles (%g) is less than or equal to  zero',Sh0);
 end

if Ah0+Ih0<=0 
     error('Initial level of infecteds (%g) is less than or equal to zero',Ah0+Ih0);
end

if betaA+betaI<=0  
    error('Transmission rate betaA+betaI (%g) is less than or equal to zero',betaA+beta);
end

if betaV<=0  
    error('Transmission rate betaV (%g) is less than or equal to zero',betaA+beta);
end

if gammaA+gammaI<=0 
     error('Recovery rate gammaA+gammaI (%g) is less than or equal to zero',gammaA+gammaI);
end

if muh<=0 
    error('Death and birth rate gamma (%g) is less than or equal to zero',muh);
end

if MaxTime<=0 
    error('Maximum run time (%g) is less than or equal to zero',MaxTime);
end

if Sh0+Ah0+Ih0>Nh0
warning('Initial level of susceptibles+infecteds (%g+%g=%g) is greater than human population size (%g)'...
,Sh0,Ah0,Ih0,Sh0+Ah0+Ih0,Nh0);
end

if Sv0+Ev0+Iv0>Nv0
warning('Initial level of susceptibles+infecteds (%g+%g=%g) is greater than mosquito population size (%g)'...
,Sv0,Ev0,Iv0,Sv0+Ev0+Iv0,Nv0);
end

Sh=Sh0; Ah=Ah0; Ih=Ih0; Rh=Rh0; Se=Se0; Sv=Sv0; Ev=Ev0; Iv=Iv0;

% The main iteration 
options = odeset('RelTol', 1e-6);
[t, pop]=ode45(@Diff_Framework_1,[0 MaxTime],[Sh Ah Ih Rh Se Sv Ev Iv],options,...
    [betaA betaI betaV gammaA gammaI muh p muv omega epsilon]);

Sh=pop(:,1); Ah=pop(:,2); Ih=pop(:,3); Rh=pop(:,4); Se=pop(:,5); Sv=pop(:,6); Ev=pop(:,7); Iv=pop(:,8);

%plot the graphs with scaled colours
plot(t,Ah,'Color',[250/255,228/255,32/255],'LineWidth',1.4);
hold on
%hold on makes matlab plot to the same plot.
plot(t,Ih,'-r','LineWidth',1)
set(gca, 'FontSize', 14)
xlabel('Time in Days','FontSize',18);
ylabel('Human Population','FontSize',18);
% ylim([0,1e5])
name = strcat('Vector_dep_M2@_c=',num2str(c),'_d=',num2str(d),'_p=',num2str(p));
Human_pop_plot_name = strcat('Human_pop_',name,'.fig');
savefig(Human_pop_plot_name);
hold off
%hold off stops matlab ploting to the same plot.
%Saves graph
%creates string based around parameters used called dataset_name 

Both_pops_csv_names = strcat('Both_pops_',name,'.csv');

Both_pops_csv = dataset({t,'Time_in_Days'},{Sh+Ah+Ih+Rh,'Total_Human_Population'},...
{Sh,'Susceptibles_Humans'},{Ah,'Asymptomatic_Humans'},{Ih,'Symptomatic_Humans'},...
{Rh,'Recovered_Humans'},{Se+Sv+Ev+Iv,'Total_Mosquito_Population'},...
{Sv+Ev+Iv,'Total_Adult_Mosquito_Population'},{Se,'Pre_Adult_Mosquito_Population'},...
    {Sv,'Susceptible_Mosquito_Population'},     {Ev,'Incubating_Mosquito_Population'},...
    {Iv,'Infectious_Mosquito_Population'});
%
export(Both_pops_csv,'File',Both_pops_csv_names,'Delimiter',',')
%creats csv file of outputs



% Calculates the differential rates used in the integration.
function dPop=Diff_Framework_1(t,pop, parameter)

betaA=parameter(1); betaI=parameter(2); betaV=parameter(3); gammaA=parameter(4);... 
    gammaI=parameter(5); muh=parameter(6); p=parameter(7); muv=parameter(8); ...
    omega=parameter(9); epsilon=parameter(10);

Sh=pop(1); Ah=pop(2); Ih=pop(3); Rh=pop(4); Nh=Sh+Ah+Ih+Rh; ...
    Se=pop(5); Sv=pop(6); Ev=pop(7); Iv=pop(8); Nv=Sv+Ev+Iv; 

dPop=zeros(8,1);

dPop(1)= muh*Nh - Sh/Nh*(betaV*Iv) - muh*Sh;
dPop(2)= (1-p)*Sh/Nh*(betaV*Iv)-(gammaA+muh)*Ah;
dPop(3)= p*Sh/Nh*(betaV*Iv) - (gammaI+muh)*Ih;
dPop(4)= gammaA*Ah+gammaI*Ih - muh*Rh;
dPop(5)= muv*Nv - omega*Se;
dPop(6)= omega*Se - (Sv*betaA*Ah/Nh + Sv*betaI*Ih/Nh) - muv*Sv;
dPop(7)= (Sv*betaA*Ah/Nh + Sv*betaI*Ih/Nh) - epsilon*Ev - muv*Ev;
dPop(8)= epsilon*Ev - muv*Iv;