Finding the closest point right after a value in ode45
1 view (last 30 days)
Show older comments
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 Comments
Answers (1)
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
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!