ode15i issue for long simulation times
조회 수: 91 (최근 30일)
이전 댓글 표시
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);
%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
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
답변 (2개)
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
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]
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]
[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 !
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
h = 5e-8
tnew = t + h
h = tnew-t
Maybe there are possibilities for a kind of normalization of the simulation times.
Best,
Torben
댓글 수: 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!