ODE45 stop after all events have fired
12 views (last 30 days)
Show older comments
I have a code in which I'm integrating the motion of several particles. I want the integration to terminate when all of the particles have had at least one crossing point in their velocity. So I tried an event code that looks like
function [val, isterm, dir] = event(t,y)
global stopped;
val = y(length(y)/2+1:end);
stopped(val==0)=1;
val = [val; all(stopped)];
isterm = 0*val; isterm(end)=1;
dir = 0*val;
But the line
stopped(val==0)=1
doesn't indicate when an event happened, because events is looking for a zero crossing. So, how can I tell (inside the event function) that the event is going to fire? Or is there an easier way to do this?
0 Comments
Answers (1)
Teja Muppirala
on 20 Nov 2012
One way of doing this would be to call the ODE solver repeatedly in a loop, until you find that all particles have set off events. Each iteration will start off at the previous time and state. This has 3 particles moving with different velocities, and integrates until they all cross zero.
function my_ode_fun
global elist
options = odeset('Events',@event_function);
elist = zeros(1,3);
tstart = 0;
tend = 100;
ystart = [1;2;3];
t_all = [];
y_all = [];
figure;
hold on;
while ~all(elist)
[tout,yout,te,ye,enum] = ode45(@(t,x) [-1;-.25;-.5], [tstart tend],ystart,options);
elist(enum(end)) = 1;
ystart = ye(end,:);
tstart = te(end,:);
t_all = [t_all; tout];
y_all = [y_all; yout];
plot(tstart*[1 1], [-5 5],'k');
end
plot(t_all, y_all);
grid on;
function [value,isterminal,direction] = event_function(t,x)
global elist
value = [x(1) x(2) x(3)]; % when value = 0, an event is triggered
isterminal = ~elist;
direction = [0 0 0];
2 Comments
Jan
on 20 Nov 2012
@Eoin: When it works, the correct result is a good compensation for a lack of elegance.
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!