Finding the closest point right after a value in ode45

1 visualización (últimos 30 días)
Mehdi el 28 de Feb. de 2017
Editada: Jan el 28 de Feb. de 2017
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);
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.

Respuestas (1)

Jan el 28 de Feb. de 2017
Editada: Jan el 28 de Feb. de 2017
Matlab's integrators cannot handle discontinuities reliably, see 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, 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


Más información sobre Ordinary Differential Equations en Help Center y File Exchange.


Community Treasure Hunt

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

Start Hunting!

Translated by