Warning: Failure at t=1.664683e+01. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (5.684342e-14) at time t.
2 次查看(过去 30 天)
显示 更早的评论
I recieve the following error:
Warning: Failure at t=1.664683e+01. Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (5.684342e-14) at time t.
> In ode15s (line 655)
I'm adjusting code for Method of Lines solving a PDE system, and the only thing that I adjusted from my previous computation is the parameter omega (located in the Initial_Parameters(t) function). I've used this same code and ran it successfully with a different omega function, but once, I use this newly "modified" omega the code doesn't produce values for ALL 29 time steps like it's mentioned in the error message above.
Here's my code: (Note: I haven't included every single function involved with the code since it would be too cumbersome for the page)
h = (x(end) - x(1))/Nx;
t0 = 9; tf = 37; tspan = [t0:1:tf]; Nt = length(tspan); % <- Time
[S0,E0,I_S0,I_A0,R0,V0] = Initial_Pop();
IC = [S0, E0, I_S0, I_A0, R0, V0].'; % <- Storing Initial Conditions u(t=0,x)
% Solver:
% options = odeset('RelTol',1e-5,'Stats','on');
[t U] = ode15s(@MOL,tspan,IC,options,Nx,h);
% Method of Lines Function:
function dUdt = MOL(t,U,Nx,h)
% Initializing Derivatives: Preallocation
dUdt=zeros(Nx,1);
dSdt = zeros(Nx,1); d2Sdx2 = zeros(Nx,1); % Susceptible
dEdt = zeros(Nx,1); d2Edx2 = zeros(Nx,1); % Exposed
dISdt = zeros(Nx,1); d2ISdx2 = zeros(Nx,1); % Infected Symptomatic
dIAdt = zeros(Nx,1); d2IAdx2 = zeros(Nx,1); % Infected Asymptomatic
dRdt = zeros(Nx,1); d2Rdx2 = zeros(Nx,1); % Recovered
dVdt = zeros(Nx,1); d2Vdx2 = zeros(Nx,1); % Environmental Infected
% Initial Population:
S(1:Nx) = U(1:Nx);
E(1:Nx) = U(Nx+1:2*Nx);
I_S(1:Nx) = U(2*Nx+1:3*Nx);
I_A(1:Nx) = U(3*Nx+1:4*Nx);
R(1:Nx) = U(4*Nx+1:5*Nx);
V(1:Nx) = U(5*Nx+1:6*Nx);
% Calling Model's Parameters:
[D1,D2,B_e,B_s,B_a,B_v,theta,lambda,gamma,Omega,rho_e,rho_s,rho_a] = Initial_Parameters(t);
% Boundary Conditions: (u(t,x))_x = 0 <- No Flux BC
S(3) = 1./3.*(4*S(2) - S(1)); % dS/dx = 0 at x = x0
S(Nx+1) = 1./3.*(4*S(Nx) - S(Nx-1)); % dS/dx = 0 at x = xf
E(3) = 1./3.*(4*E(2) - E(1)); % dE/dx = 0 at x = x0
E(Nx+1) = 1./3.*(4*E(Nx) - E(Nx-1)); % dE/dx = 0 at x = xf
I_S(3) = 1./3.*(4*I_S(2) - I_S(1)); % dI_S/dx = 0 at x = x0
I_S(Nx+1) = 1./3.*(4*I_S(Nx) - I_S(Nx-1)); % dI_S/dx = 0 at x = xf
I_A(3) = 1./3.*(4*I_A(2) - I_A(1)); % dI_A/dx = 0 at x = x0
I_A(Nx+1) = 1./3.*(4*I_A(Nx) - I_A(Nx-1)); % dI_A/dx = 0 at x = xf
R(3) = 1./3.*(4*R(2) - R(1)); % dR/dx = 0 at x = x0
R(Nx+1) = 1./3.*(4*R(Nx) - R(Nx-1)); % dR/dx = 0 at x = xf
V(3) = 1./3.*(4*V(2) - V(1)); % dV/dx = 0 at x = x0
V(Nx+1) = 1./3.*(4*V(Nx) - V(Nx-1)); % dV/dx = 0 at x = xf
% f(t,x,U):
F1 = @(S,E,I_S,I_A,V) -(B_e*E + B_s*I_S + B_a*I_A + B_v*V).*S;
F2 = @(S,E,I_S,I_A,V) (B_e*E + B_s*I_S + B_a*I_A + B_v*V).*S - lambda*E;
F3 = @(E,I_S) theta*lambda*E - lambda*I_S;
F4 = @(E,I_A) (1 - theta)*lambda*E - lambda*I_A;
F5 = @(I_S,I_A) gamma*(I_S + I_A);
F6 = @(E,I_S,I_A,V) rho_e*E + rho_s*I_S + rho_a*I_A - Omega*V;
for i = 2:Nx
% Eqn 1:
d2Sdx2(i) = 1./(h.^2).*(S(i+1) - 2.*S(i) + S(i-1)); % Centered Finite Difference
dSdt(i) = D1.*d2Sdx2(i) + F1(S(i), E(i), I_S(i), I_A(i), V(i));
% Eqn 2:
d2Edx2(i) = 1./(h.^2).*(E(i+1) - 2.*E(i) + E(i-1));
dEdt(i) = D1.*d2Edx2(i) + F2(S(i), E(i), I_S(i), I_A(i), V(i));
% Eqn 3:
d2ISdx2(i) = 1./(h.^2).*(I_S(i+1) - 2.*I_S(i) + I_S(i-1));
dISdt(i) = D1.*d2ISdx2(i) + F3(E(i), I_S(i));
% Eqn 4:
d2IAdx2(i) = 1./(h.^2).*(I_A(i+1) - 2.*I_A(i) + I_A(i-1));
dIAdt(i) = D1.*d2IAdx2(i) + F4(E(i), I_A(i));
% Eqn 5:
d2Rdx2(i) = 1./(h.^2).*(R(i+1) - 2.*R(i) + R(i-1));
dRdt(i) = D1.*d2Rdx2(i) + F5(I_S(i), I_A(i));
% Eqn 6:
d2Vdx2(i) = 1./(h.^2).*(V(i+1) - 2.*V(i) + V(i-1));
dVdt(i) = D2.*d2Vdx2(i) + F6(E(i), I_S(i), I_A(i), V(i));
end
dUdt = [dSdt; dEdt; dISdt; dIAdt; dRdt; dVdt];
end
% Calculating Initial Population Function:
function [S0,E0,I_S0,I_A0,R0,V0] = Initial_Pop()
% 1. Calculating Initial Population Distribution:
% Initial Susceptible Population Density:
S0 = Initial_Cond_1D(Lat,Pop,1,"Initial Susceptible $$S(0,x)$$ Population Distribution");
n = length(S0);
% Initial Exposed Population Density:
% Figure 2: E0 = 80.7
E0 = zeros(1,n);
% Initial Symptomatic Infected Population Density:
I_S0 = Initial_Cond_1D(Lat,Case_Data,2,"Initial Infected Symptomatic $$I_{S}(0,x)$$ Population Distribution"); % Discrete uniform distribution whose range is in the close interval [1, 10]
% Initial Asymptomatic Infected Population Density:
I_A0 = zeros(1,n);
% Initial Recovered Population Density:
R0 = zeros(1,n);
% Initial Environment Population Density:
V0 = zeros(1,n);
end
% Establishing Parameter Values
function [D1,D2,B_e,B_s,B_a,B_v,theta,lambda,gamma,omega,rho_e,rho_s,rho_a] = Initial_Parameters(t)
% Parameters:
D1 = 10^(-4); % Diffusion Coefficient
D2 = 10^(-8); % Diffusion Coefficient
B_e = 1.15*10^(-10);
B_s = 8.92*10^(-8);
B_a = 2.67*10^(-8);
B_v = 8.36*10^(-18);
theta = 0.6;
lambda = 0.21;
gamma = 0.07;
omega = Omega_Surface_Trans(7*t) % New OMEGA
% omega = Omega_Temp(7*t); % time's units are in WEEKS so (7*t) is in DAYS
% <- Original Omega (RUNS SUCCESSFULLY)
rho_e = 1223;
rho_s = 9.13*10^(8);
rho_a = 6.54*10^(5);
end
% Omega Function:
function [omega] = Omega_Temp(t)
% Omega is rate at which the virus remains infectious
% in the environment. Thus, we've added a temporal component
% since Omega would varry depending on the temperature.
T = Final_Temperature(t);
Num = log(10);
a = 1;
b = 1;
Denom = exp(a*T + b);
omega = Num/Denom;
end
% Real Omega function for Steel:
function [omega] = Omega_Surface_Trans(x)
T = Final_Temperature(x);
% Determining a and b values:
if (T <= 68)
a = -0.23624;
b = 5.5799;
elseif (68 < T && T <= 86)
a = -1.331;
b = 5.839;
else
a = -4.93;
b = 5.63;
end
% a = -0.23624; b = 5.5799;
Num = log(10); % ln(10)
Denom = exp(a*T + b);
omega = Num/Denom;
end
7 个评论
Steven Lord
2022-10-20
Can you add the option 'OutputFcn' with value @odeplot to the options structure you pass into ode15i and show us the plot created, like the one below? I'm curious if one or more of the lines starts increasing or decreasing very rapidly as the solver starts approaching t = 1.664683e+01 which could account for why the solver starts trying smaller and smaller steps as it approaches the cliff or chasm.
options = odeset('OutputFcn', @odeplot);
sol = ode45(@(t, y) y, [0 5], 1, options);
回答(1 个)
Abhijeet
2023-10-23
Hi Thomas,
I understand that you are facing the following warning message while solving a PDE system.
To resolve the warning, you need to set the absolute and relative tolerances to a higher value. In case of sharp discontinuities in the reference input such as when tracking a reference signal, it is possible that the encountered warning is generated at intermediate time steps.
MATLAB tries to reduce the time step to a small amount to counter the abrupt change due to the discontinuity in the reference signal. For sharp discontinuities, it might not be possible to avoid this warning. However, for non-discontinuous inputs, we can set relative and absolute tolerance to a higher number than the default setting.
I suggest you to use the following commands to set the tolerances to a higher value:
options = odeset('Mass',M,'RelTol',1e-4,'AbsTol',1e-6);
[t,y] = ode15s(@objectivefunction,tspan,y0,options);
Please check the following reference to know more about relative and absolute tolerance:
I hope this resolves the issue you were facing.
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!