Finding the closest point right after a value in ode45
2 ビュー (過去 30 日間)
古いコメントを表示
Hi,
I have the following code which I'm trying to solve using ode45 function.
v = y(1);
gamma = y(2);
h = y(3);
g = 9.2;
rm = 1200;
if h < 10
v_dot = D - g;
gamma_dot = 0;
h_dot = v;
elseif h > 10
v_dot = D - g*sin(gamma);
gamma_dot = v*(g - v/Rm)*cos(gamma);
h_dot = v*sin(gamma);
end
dydt(1) = v_dot;
dydt(2) = gamma_dot;
dydt(3) = h_dot;
so when h is less than 10, certain actions are done. And when h is more than 10, I want v_dot and gamma_dot take other values. So far all is nice and dandy but what I need to add to this is gamma_dot has to take the 0.1 value exactly when h is 10 or the immediate value that it takes after 10 (closest value to 10). Problem is, h may never get exactly 10 so how can I tell MATLAB that gamma_dot should be 0.1 when h is 10 or the very next immediate value that it takes after 10? Simply adding:
elseif h == 10; gamm5a_dot = 0.1;
does not help. Can somebody shed some light on this please? Thanks a bunch.
0 件のコメント
回答 (1 件)
Jan
2017 年 2 月 28 日
編集済み: Jan
2017 年 2 月 28 日
Matlab's integrators cannot handle discontinuities reliably, see http://www.mathworks.com/matlabcentral/answers/59582#answer_72047. If the function to be integrated contains discontinuities, you run ODE45 outside its specifications and you cannot expect an accurate result, or any result at all.
To switch a parameter at y(3)=10 use an event function (see https://www.mathworks.com/help/matlab/ref/odeset.html#input_argument_namevalue_Events), which stops the integeration, when this point is reached. Then change the parameter and restart the integration using the final value of the former one as start value of this one. Pseudo-code:
tIni = 0;
y0 = ...
tFin = 10;
parameter = ...
Control = odeset('Events', @myEventFcn);
while tIni < tFin
[t, y] = ode45(@(t,y), myFcn(t,y,parameter), [tIni, tFin], y0, Control);
tIni = t(end);
paramter = ... % Set new parameters depending on y(end, :)
y0 = y(end, :); % New initial value is old final value
end
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!