Finding the closest point right after a value in ode45
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
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 comentarios
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 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 comentarios
Ver también
Categorías
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!