ode15i issue for long simulation times

We discovered an issue with long simulation times in ode15i, which is causing large minimal step sizes and we assume therfore an error.

We assume the error to occurs in ode15i line 198 (hmin = 16*eps*abs(t);) and with the same command in line 272. Also we see in issue in 291 tnew = t + h; and 295 h = tnew - t; with the double precision.

The large simulation time t causes large hmin and triggers the warning MATLAB:ode15i:IntegrationTolNotMet for the minimum allowed step size. Also, ode15i does not allow to set minimum step size.

Does anyone know a good prcatice to overcome this?

We tried to express hmin as hmin = 16*eps*(t~=0); and commented out line 295 for the minimal example. This works but there might be different issues

Below see the a minimal example:

%minimal example based on documentation in https://de.mathworks.com/help/matlab/ref/ode15i.html

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], ...
   'Jacobian',{[],[1 0 0; 0 1 0; 0 0 0]});
%Use decic to compute consistent initial conditions from guesses. Fix the first two components of y0 to get the same consistent initial conditions as found by ode15s in hb1dae.m, which formulates this problem as a semi-explicit DAE system.

y0 = [1; 0; 1e-3];
yp0 = [0; 0; 0];
[y0,yp0] = decic(@robertsidae,0,y0,[1 1 0],yp0,[],options);
%Solve the system of DAEs using ode15i.

tspan = [0 1e3]+1e9; 
[t,y] = ode15i(@robertsidae,tspan,y0,yp0,options);

Warning: Failure at t=1.000000e+09. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-06) at time t.

%Plot the solution components. Since the second solution component is small relative to the others, multiply it by 1e4 before plotting.

y(:,2) = 1e4*y(:,2);
semilogx(t,y)
ylabel('1e4 * y(:,2)')
title('Robertson DAE problem with a Conservation Law, solved by ODE15I')

function res = robertsidae(t,y,yp)
res = [yp(1) + 0.04*y(1) - 1e4*y(2)*y(3);
   yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2;
   y(1) + y(2) + y(3) - 1];
end

The simulation time [0 1e3]+1e9 triggers:

Warning: Failure at t=1.000000e+09. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-06) at time t.

With our modifacion we get a solution.

PS. for our actual problem the simualtion time was way lower but switches triggered this issue earlier at t=682200

[Francisco J. Triveno Vargas]

Hi,

I change the tspan and run normaly, please verify your time interval.

clear,

close all,

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], ...

'Jacobian',{[],[1 0 0; 0 1 0; 0 0 0]});

%Use decic to compute consistent initial conditions from guesses. Fix the first two components of y0 to get the same consistent initial conditions as found by ode15s in hb1dae.m, which formulates this problem as a semi-explicit DAE system.

y0 = [1; 0; 1e-3];

yp0 = [0; 0; 0];

[y0,yp0] = decic(@robertsidae,0,y0,[1 1 0],yp0,[],options);

%Solve the system of DAEs using ode15i.

% tspan = [0 1e3]+1e9;

tspan = [0:0.001:1e3];

[t,y] = ode15i(@robertsidae,tspan,y0,yp0,options);

y(:,2) = 1e4*y(:,2);

semilogx(t,y)

ylabel('1e4 * y(:,2)')

title('Robertson DAE problem with a Conservation Law, solved by ODE15I')

function res = robertsidae(t,y,yp)

res = [yp(1) + 0.04*y(1) - 1e4*y(2)*y(3);

yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2;

y(1) + y(2) + y(3) - 1];

end

[Aline Luxa]

Dear Francisco J. Triveno Vargas,

thank you for your answer!

We were referring to the absolute value of the time vector. Therefore we added 1e9 to the time interval: tspan = [0 1e3]+1e9

By using tspan = [0:0.001:1e3]+1e9 it also triggers the issue.

In that case the time is really high but for our actual problem the issue (we assume with hmin) is way earlier at t=682200

[Alan Stevens]

Removing "options" from ODE15i gives the following:

%minimal example based on documentation in https://de.mathworks.com/help/matlab/ref/ode15i.html

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], ...
   'Jacobian',{[],[1 0 0; 0 1 0; 0 0 0]});
%Use decic to compute consistent initial conditions from guesses. Fix the first two components of y0 to get the same consistent initial conditions as found by ode15s in hb1dae.m, which formulates this problem as a semi-explicit DAE system.

y0 = [1; 0; 1e-3];
yp0 = [0; 0; 0];
[y0,yp0] = decic(@robertsidae,0,y0,[1 1 0],yp0,[],options);
%Solve the system of DAEs using ode15i.

tspan = [0 1e3]+1e9; 
[t,y] = ode15i(@robertsidae,tspan,y0,yp0);
%Plot the solution components. Since the second solution component is small relative to the others, multiply it by 1e4 before plotting.

y(:,2) = 1e4*y(:,2);
semilogx(t,y)
ylabel('1e4 * y(:,2)')
title('Robertson DAE problem with a Conservation Law, solved by ODE15I')

function res = robertsidae(~,y,yp)
res = [yp(1) + 0.04*y(1) - 1e4*y(2)*y(3);
   yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2;
   y(1) + y(2) + y(3) - 1];
end

[Aline Luxa]

Dear @Alan Stevens,

thanks for the suggestions. I think by removing the options the conditions changed e.g. the tolerances and that is maybe why a larger time step could result into an accepted solution?

In our problem this didn't work out because relaxing the tolerances resulted into bad solutions. Although a solution could be provided and we didn't get the warning anymore.

Answer 1

[Shivam Gothi]

Hello @Aline Luxa,

As per my understanding, you are trying to slove the ordinary diffrential equation for a longer time span, but it throws some warnings and nothing shows up on the plot. Here is one possible solution:

when you write:

tspan = [0 1e3]+1e9

tspan = 1×2

1.0e+09 * 1.0000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

you can see that the tspan vector is not correctly defined (it adds 1e9 to "0" and "1e3") . Instead of it, change the tspan vector to:

t_start = 0;
t_end = 1e9;
tspan = [t_start t_end]

tspan = 1×2

1.0e+09 * 0 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

I made the above mentioned changes and executed the code as shown below:

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], ...
   'Jacobian',{[],[1 0 0; 0 1 0; 0 0 0]});
%Use decic to compute consistent initial conditions from guesses. Fix the first two components of y0 to get the same consistent initial conditions as found by ode15s in hb1dae.m, which formulates this problem as a semi-explicit DAE system.
y0 = [1; 0; 1e-3];
yp0 = [0; 0; 0];
[y0,yp0] = decic(@robertsidae,0,y0,[1 1 0],yp0,[],options);
%Solve the system of DAEs using ode15i.
t_start = 0;
t_end = 1e9;
tspan = [t_start t_end] 

tspan = 1×2

1.0e+09 * 0 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[t,y] = ode15i(@robertsidae,tspan,y0,yp0,options);
%Plot the solution components. Since the second solution component is small relative to the others, multiply it by 1e4 before plotting.
y(:,2) = 1e4*y(:,2);
semilogx(t,y)
ylabel('1e4 * y(:,2)')
title('Robertson DAE problem with a Conservation Law, solved by ODE15I')

function res = robertsidae(t,y,yp)
res = [yp(1) + 0.04*y(1) - 1e4*y(2)*y(3);
   yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2;
   y(1) + y(2) + y(3) - 1];
end

I hope this resolves the issue !

[Aline Luxa]

@Shivam Gothi Thank you for your proposition!

I think our display of tspan was confusing and not ideal.

We wanted to show that ode15i cannot solve at too high time steps, for our problem it is probably because of hmin (which gets too large). Apparently, it is also connected to the starting time. We have a time vector of more points in our problem, this is why we didn't start this example at t_start = 0;. This should not effect a time invariant problem.

By starting later, we wanted to show what happens during the simulation of our problem when it is advanced in simulation time.

Like for this example:

options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], ...
    'Jacobian',{[],[1 0 0; 0 1 0; 0 0 0]});
%Use decic to compute consistent initial conditions from guesses. Fix the first two components of y0 to get the same consistent initial conditions as found by ode15s in hb1dae.m, which formulates this problem as a semi-explicit DAE system.
y0 = [1; 0; 1e-3];
yp0 = [0; 0; 0];
[y0,yp0] = decic(@robertsidae,0,y0,[1 1 0],yp0,[],options);
%Solve the system of DAEs using ode15i.
t_start = 1e9;
t_end = 1e10;
tspan = [t_start t_end] 

tspan = 1×2

1.0e+10 * 0.1000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[t,y] = ode15i(@robertsidae,tspan,y0,yp0,options);

Warning: Failure at t=1.000000e+09. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-06) at time t.

%Plot the solution components. Since the second solution component is small relative to the others, multiply it by 1e4 before plotting.
y(:,2) = 1e4*y(:,2);
semilogx(t,y)
ylabel('1e4 * y(:,2)')
title('Robertson DAE problem with a Conservation Law, solved by ODE15I')

function res = robertsidae(t,y,yp)
res = [yp(1) + 0.04*y(1) - 1e4*y(2)*y(3);
    yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2;
    y(1) + y(2) + y(3) - 1];
end

Answer 2

[Torben Warnecke]

Hey,

I think the Problem really lies in large simulation times tspan.

In ode15i the absolute value of the current simulation time t is used to calculate hmin (see line 198 of ode15i: hmin = 16*eps*abs(t)) which leads to large hmin when t is becoming large.

Additionally there is the problem of double precision of t. E.g. when in line 291 (tnew = t + h) t is bigger than 1e9 and h is less then 5e-8, the difference will be cut off and tnew will be similar to t. Then, with line 295 (h = tnew - t) h will be set to 0 leading to inf and NaN in the procedure afterwards (in line 307, where division by h occurs).

t = 1e9

t = 1.0000e+09

h = 5e-8

h = 5.0000e-08

tnew = t + h

tnew = 1.0000e+09

h = tnew-t

h = 0

Maybe there are possibilities for a kind of normalization of the simulation times.

Best,

Torben