ode15i issue for long simulation times

9 Ansichten (letzte 30 Tage)
Aline Luxa
Aline Luxa am 8 Okt. 2024
Kommentiert: Aline Luxa am 9 Okt. 2024
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
  4 Kommentare
Alan Stevens
Alan Stevens am 9 Okt. 2024
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
Aline Luxa am 9 Okt. 2024
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.

Melden Sie sich an, um zu kommentieren.

Antworten (2)

Shivam Gothi
Shivam Gothi am 9 Okt. 2024
Bearbeitet: Shivam Gothi am 9 Okt. 2024
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 !
  1 Kommentar
Aline Luxa
Aline Luxa am 9 Okt. 2024
@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

Melden Sie sich an, um zu kommentieren.


Torben Warnecke
Torben Warnecke am 9 Okt. 2024
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

Produkte


Version

R2023b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by