ODE45 stop after all events have fired

12 views (last 30 days)
Eoin
Eoin on 20 Nov 2012
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?

Answers (1)

Teja Muppirala
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
Eoin
Eoin on 20 Nov 2012
Yes, this would work, but isn't as elegant/efficient as I was hoping for..
I finally settled for running the particles in their own ode45 calls, but this only works because the particles are actually non-interacting.
Jan
Jan on 20 Nov 2012
@Eoin: When it works, the correct result is a good compensation for a lack of elegance.

Sign in to comment.

Categories

Find more on Programming in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!