MATLAB Answers

Finding the closest point right after a value in ode45

4 views (last 30 days)
Mehdi
Mehdi on 28 Feb 2017
Edited: Jan on 28 Feb 2017
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.

Answers (1)

Jan
Jan on 28 Feb 2017
Edited: Jan on 28 Feb 2017
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

Community Treasure Hunt

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

Start Hunting!