ODE45 solver results changes. How to select optimal Reltol and Abstol?

122 vues (au cours des 30 derniers jours)
Praveen Kumar
Praveen Kumar le 15 Nov 2019
Hi,
I am trying to solve ODE45 for solving an initial value problem. Part of the code is shown here to explain, if required I can attach the full code.
My doubt is, the output from the mentioned state equation seems changing depending on RelTol and Abstol. i.e X. When I was using default, some values are wrong from the required results. Some values are increasing instead of decreasing. Then I have increased tolerance to 1e-6 and 1e-8 , and found the results improving.
How can I decide exactly which tolerance I should change, Reltol or Abstol? Also, how to decide the tolerance value?
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
Tu = linspace(t0,tf,t_segment);
[Tx,X] = ode45(@(t,x) stateEq_s(t,x,u,Tu), Tu, initx, options); % State Equation
% Value of X increases where it should have decreased for tolerances such as 1e-3
% stateEq_S function
function dx = stateEq_s(t,x,u,Tu)
global C Rs Rp Vsmax
dx = 0;
load PV_Preq_data.mat;
load PV_3kW_2005_2016M.mat; % Load demand Power
kW = PV_2016M(25200:25200+32400-1);
Preq = Pdc - kW;
Preq = Preq (1:32400);
Preq(1)=0;
u = interp1(Tu,u,t); % Interploate the control at ode45 time t
Preq = interp1(Tu,Preq,t); % Interpolate the Preq at ode45 time t
for j=1:1:length(Preq)
if Preq(j) > 0
dx = -1/C*(( 1/(2*Rs)+ 1/Rp )*x(j) - sqrt((x(j)).^2 - 4*Rs*(Preq(j) - u(j))/(Vsmax^2))/(2*Rs));
else if Preq(j) < 0
dx = -1/C*(( 1/(2*Rs) - 1/Rp )*x(j) - sqrt((x(j)).^2 - 4*Rs*(Preq(j) - u(j))/(Vsmax^2))/(2*Rs));
else
dx = -1/C*(( 1/(2*Rs))*x(j) - sqrt((x(j)).^2 - 4*Rs*(Preq(j) - u(j))/(Vsmax^2))/(2*Rs));
end
end
end
Also, sometimes I get 'NaN' values depending on tolerances.
  5 commentaires
Praveen Kumar
Praveen Kumar le 15 Nov 2019
I will try this. Thanks.
Ioannis Matthaiou
Ioannis Matthaiou le 15 Avr 2021
Have you found out the source of the problem? I think by having such a switching function may cause instability in the solutions..

Connectez-vous pour commenter.

Réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by