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

118 ビュー (過去 30 日間)
Praveen Kumar
Praveen Kumar 2019 年 11 月 15 日
コメント済み: Ioannis Matthaiou 2021 年 4 月 15 日
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 件のコメント
Praveen Kumar
Praveen Kumar 2019 年 11 月 15 日
I will try this. Thanks.
Ioannis Matthaiou
Ioannis Matthaiou 2021 年 4 月 15 日
Have you found out the source of the problem? I think by having such a switching function may cause instability in the solutions..

サインインしてコメントする。

回答 (0 件)

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by