필터 지우기
필터 지우기

ode15i issue for long simulation times

조회 수: 91 (최근 30일)
Aline Luxa
Aline Luxa 2024년 10월 8일 15:34
댓글: Aline Luxa 2024년 10월 9일 11:19
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
Alan Stevens
Alan Stevens 2024년 10월 9일 9:28
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 2024년 10월 9일 11:19
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.

댓글을 달려면 로그인하십시오.

답변 (2개)

Shivam Gothi
Shivam Gothi 2024년 10월 9일 9:29
편집: Shivam Gothi 2024년 10월 9일 9:30
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
Aline Luxa
Aline Luxa 2024년 10월 9일 11:13
@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

댓글을 달려면 로그인하십시오.


Torben Warnecke
Torben Warnecke 2024년 10월 9일 10:11
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

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

제품


릴리스

R2023b

Community Treasure Hunt

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

Start Hunting!

Translated by