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.

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

[Torsten]

You should at least include the function "Final_Temperature" from which omega is computed.

Most probably, the solver stops at a time where the parameters a and b discontinuously change:

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

Use the event option of the ODE integrators to stop and restart the solver when T reaches 68 and 86.

And you mustn't approximate values at the third, but at the first discretization point if you want to set the boundary condition d/dx = 0:

S(3) = 1./3.*(4*S(2) - S(1)); % dS/dx = 0 at x = x0
E(3) = 1./3.*(4*E(2) - E(1)); % dE/dx = 0 at x = x0
I_S(3) = 1./3.*(4*I_S(2) - I_S(1)); % dI_S/dx = 0 at x = x0
I_A(3) = 1./3.*(4*I_A(2) - I_A(1)); % dI_A/dx = 0 at x = x0
R(3) = 1./3.*(4*R(2) - R(1)); % dR/dx = 0 at x = x0
V(3) = 1./3.*(4*V(2) - V(1)); % dV/dx = 0 at x = x0

[Thomas Rodriguez]

@Torsten

Here's the final temperature function:

function [T] = Final_Temperature(x)
f_s = 64.712 + -6.9178*sin((2*pi/365)*x + 0.99605);
f_d = -4.3372*sin((2*pi/1)*x + -5.0872);
T = f_s + f_d;
end

How would I implement the event option for odes15()? The only documentation that I could find mentioning that is here: https://www.mathworks.com/help/matlab/math/summary-of-ode-options.html

As for the boundary conditions, how could I apply the finite method at the first discretization point? Would it use the "ghost" point idea? (i.e. assigning values for U(-1) & U(0) to approximate the no flux boundary value at U(1) where U(x,t) is the discretization of the solution for an arbitrary PDE).

[Davide Masiello]

Regarding the ODE Event Location

https://ch.mathworks.com/help/matlab/math/ode-event-location.html

[Torsten]

@Thomas Rodriguez

And is there a coincidence between the time the solver quits (t = 16,64...) and T = 68 or 86 ?

Concerning the approximation of the derivatives at the boundary:

https://en.wikipedia.org/wiki/Finite_difference_coefficient

Chapter "Forward Finite Difference"

[Thomas Rodriguez]

@Torsten

I'm not sure why the solver quits at the specific time. I agree with your previous comment that it probably arises from the if statement in Omega_Surface_Trans(x). Instead of constructing the if statement as that. Would it be easier by creating 3 different temperature arrays and measuring the error of T with each of those arrays. Then if the error reaches a predetermined tolerance then choose the related a and b values from there.

[Steven Lord]

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);

[Thomas Rodriguez]

@Steven Lord

This is the result using the original code seen above. I also turned on Statistics to see if anything is percular.

1689 successful steps

685 failed attempts

83021 function evaluations

428 partial derivatives

929 LU decompositions

5520 solutions of linear systems

Answer 1

[Abhijeet]

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:

• https://www.mathworks.com/matlabcentral/answers/26743-absolute-and-relative-tolerance-definitions

I hope this resolves the issue you were facing.