0

Usually any ode like ode23, ode 45 will do the integration from initial time to final time [t0 tf]. Is there a way, the integration can be stopped on some other parameter not dependent on time? For example, I have a linear damper.

Initial Pressure p1 = some value
Initial Pressure p2 = some value (not = p1)
time = [t0 tf]
some other constants
options = odeset
y0 = [initial conditions for some parameters containing p1 and p2]
[t,y] = ode45(@func,[t0 tf],y0,options,other constants to carry)

and in func code:
equations for integration for p1 and p2 and some other variables

How would it be possible to not run the ode from t0 to tf but stop it when p1 = p2? Or some way I can pre decide the limits for p1 and p2 so that the ode does not exceed them? Please help. Thanks

Chris Toomer
  • 23
  • 4
  • Is it a computational effort problem, i.e. does it take too long otherwise? You can simply remove the parts of the "regular" solution from the point where p1=p2 – rst Jan 28 '16 at 07:18
  • I guess it will take a lot of effort to find and remove the values from the regular solution. Thats why I think I systematic way out of this would be better. But thanks anyway. – Chris Toomer Jan 28 '16 at 11:54

1 Answers1

1

You can use the event option for your ode.

Use

opts=odeset('Events',@events);
[t,y]=ode45(@func,[t0 tf],y0,opts);

And

function [value,isterminal,direction] = events(t,y)
% check p1 = p2
% and stop
value = y(x1)-y(x2)
isterminal = 1; % stop the integration
direction = 0; % all events

Ofc you need to set p1 and p2 first, but I don't know where in your y there are located.

Here is a working example

function mytest()
clear T
clear Y
y0 = [1 3];
options = odeset('RelTol', 1e-4, 'AbsTol', [1e-4 1e-4 1e-5], 'Events', @events, 'OutputFcn', @odeplot);
[Tn,Yn,T,Y,I] = ode45(@func, [0 12], [0 1 1], options);
odeplot([],[],'done')
hold on
figure
options = odeset('RelTol', 1e-4, 'AbsTol', [1e-4 1e-4 1e-5], 'OutputFcn', @odeplot);
[Tn,Yn,T,Y,I] = ode45(@func, [0 12], [0 1 1], options);
odeplot([],[],'done')
end


function dy = func(t,y)
dy=zeros(3,1);
dy(1) = y(2)*y(3);
dy(2) = -y(1)*y(3);
dy(3) = -0.51*y(1)*y(2);
end

function [value,isterminal,direction] = events(t,y)
value = y(1)-y(2);
isterminal = 1; % stop the integration
direction = 0; % all events
end
rst
  • 2,510
  • 4
  • 21
  • 47
  • So this events code should be integrated with the func code that I have to run the computation for ode or should it be a separate code and if separate where do I set the p1 and p2 over here because it should know those values first from the func code in order to see if they are equal or not right? Thanks – Chris Toomer Jan 28 '16 at 11:51
  • A seperate code, the p1 and p2 sould be somewhere in the y vector. – rst Jan 28 '16 at 11:52
  • No. It did not work. Just the old result coming out. – Chris Toomer Jan 28 '16 at 12:05
  • function [value,isterminal,direction] = events(tt,xx,Ap,Aorifice,Cd,rho,speed,S) % check p1 = p2 % and stop p1 = xx(3,1); p2 = xx(4,1); if(p1==p2) value = 0; else value = 1; end isterminal = 1; % stop the integration direction = 0; % all events – Chris Toomer Jan 28 '16 at 12:06