ode45 event trigger not responding

I'm trying to get the integration of ode45 to stop when the first column (which represents position) reaches a certain value. But even though I set an event trigger and the integration reaches that value, it continues to integrate until the final time and can't understand why. Is there some error in my code of the event function that I dont understand?

options = odeset('RelTol',1e-9,'AbsTol',1e-9,'Events',@event);
[tt,yy] = ode45(@ODEsystem,[0,3],[0,0,10e-3,0],options);
plot(tt,yy(:,1))

function [value,isterminal,direction] = event(~,yy)
  x_max = 200e-3;
  value = (yy(:,1)==x_max); % 
  isterminal = 1;  % Halt integration 
  direction = -1;  
end


function dydt = ODEsystem(t,y)
dydt(1) = y(2);  %  xa_dot = va;

dydt(2) = (1/ma)*((P4*A4) - (P5*A5) - (F_0+(k*y(1))));  %  va_dot = (1/ma)*((P4*A4) - (P5*A5) - (F_0+(k*xa)));

dydt(3) =  y(2)*A5;  %  Vt_dot = Qout;
 
dydt(4) =  y(2)*A4;  %  Vacc_dot = -Qin; 
end

Answer 1

[Walter Roberson]

  x_max = 200e-3;
  value = (yy(:,1)==x_max); % 

Remember that == between double precision values, is a bit-for-bit comparison. If the yy value was off by just 1 bit from 200e-3 then your test would not detect the value.

You should be using

   value = yy(:,1) - x_max;

or possibly the negative of that.

However, remember also that your input boundary conditions in yy are always a column vector (unless parallel computation was requested, in which case the boundary conditions go down columns.) So yy(:,1) is testing every boundary condition, not just the first boundary condition. It is unlikely that is what you want.

You also need to be creating a column vector of outputs in ODEsystem such as by first doing

   dydt = zeros(4,1);

Also, ma is not defined at the time it is needed in your function.

[Nader Mohamed]

@Walter Roberson Yes I'm sorry, I posted just a snippet of code because the full script was long and I just wanted to understand the event function. Below is the full code where i create the column vector in ODEsystem. I still don't understand your take on yy(:,1) testing every boundary condition and how can I solve this issue.

When I run the code with your change it still integrates until the final time and doesn't stop at xmax

close all
clear
clc

options = odeset('RelTol',1e-9,'AbsTol',1e-9,'Events',@event);
[tt,yy] = ode45(@ODEsystem,[0,3],[0,0,10e-3,0],options);
plot(tt,yy(:,1))

function [value,isterminal,direction] = event(~,yy)
  x_max = 200e-3;
  value = yy(:,1) - x_max; % The value that we want to be zero
  isterminal = 1;  % Halt integration 
  direction = -1;   % The zero can be approached from either direction
end


function z = command(t,t1,t2) %FUNCTION FOR COMMAND Z
   z=zeros(size(t));
   z(t>t1&t<t2)=(t(t>t1&t<t2)-t1)*t1/(t2-t1);
   z(t>=t2)=t1;
end




function dydt = ODEsystem(t,y)

t0 = 0;
t1 = 1;
t2 = 1.5;
tf = 3;

d0 = 5e-3;

gamma = 1.2;

L_23 = 2; %m

L_67 = 15;

rho = 890;

f67 = 0.035;
f23 = 0.032;

ka = 1.12;
kc = 2;
kd = 12;
kt = 1.12;

D_23 = 18e-3;
D_67 = 18e-3;

A4 = ((D_23/2)^2)*pi; 
A5 = ((D_67/2)^2)*pi; 

V_accumulator = 10e-3; %m^3
P_accumulator = 2.5e6; %Pa

% t=0:.01:3;
z=command(t,t1,t2);
alpha = 2*acos(1 - 2*abs(z));
A_d = (d0^2)/4*(alpha - sin(alpha));

dydt = zeros(4,1);

PT = 0.1e6;
P0 = 21e6;

ma = 2; %kg
x_max = 200e-3; %m
F_0 = 1; %N
k = 120; %N/m

V0 = (P_accumulator*V_accumulator)/P0;

PA = P_accumulator * (V0/(V_accumulator - y(4)))^gamma;

P4 = PA - ((rho*(y(2)*A4)^2)/2) *( (ka/A4^2)+ ((kc.*y(2))/A4^2) + (kd./((A_d+sqrt(eps)).^2)) + f23*(L_23/(D_23*A4^2))  );
P5 = PT - ((rho*(y(2).*A5)^2)/2) *( (kt./(A5.^2))+(kd./(A_d+sqrt(eps)).^2) + f67*(L_67/(D_67.*A5.^2)) );


% control valve instructions
if z == 0
    P4 = PA;
    P5 = P0; %P0
elseif z > 0
    if A_d < 1e-9
        P4 = PA;
        P5 = P0;
    else
        P4 = PA - ((rho*(y(2)*A4)^2)/2) *( (ka/A4^2)+ ((kc.*y(2))/A4^2) + (kd./((A_d+sqrt(eps)).^2)) + f23*(L_23/(D_23*A4^2))  );
        P5 = PT - ((rho*(y(2).*A5)^2)/2) *( (kt./(A5.^2))+(kd./(A_d+sqrt(eps)).^2) + f67*(L_67/(D_67.*A5.^2)) );
    end
else
    if A_d < 1e-9
        P4 = PA;
        P5 = P0;
    else
        P4 = PA - ((rho*(y(2)*A4)^2)/2) *( (ka/A4^2)+ ((kc.*y(2))/A4^2) + (kd./((A_d+sqrt(eps)).^2)) + f23*(L_23/(D_23*A4^2))  );
        P5 = PT - ((rho*(y(2).*A5)^2)/2) *( (kt./(A5.^2))+(kd./(A_d+sqrt(eps)).^2) + f67*(L_67/(D_67.*A5.^2)) );
    end
end


if y(1) < 0
    y(1) = 0;
end
if y(1) <= 0 && dydt(1) < 0
    y(1) = 0;
    dydt(1) = 0;
end
if y(1) > x_max
    y(1) = x_max;
end
if y(1) >= x_max && dydt(1) > 0
    y(1) = x_max;
    dydt(1) = 0;
end


dydt(1) = y(2);  %  xa_dot = va;

dydt(2) = (1/ma)*((P4*A4) - (P5*A5) - (F_0+(k*y(1))));  %  va_dot = (1/ma)*((P4*A4) - (P5*A5) - (F_0+(k*xa)));

dydt(3) =  y(2)*A5;  %  Vt_dot = Qout;
 
dydt(4) =  y(2)*A4;  %  Vacc_dot = -Qin; 

 if (y(1) == x_max && dydt(2) > 0) || (y(1)  == 0 && dydt(2) < 0)
    dydt(1) = 0;
    dydt(2) = 0;
 end
 if y(4) >= V_accumulator
    dydt(4) = 0;
end



end

[Walter Roberson]

Your system has if statements that appear to be implementing discontinuities. The ode*() routines rely on mathematics that is not compatible with discontinuities in the first two derivatives.

[Walter Roberson]

  value = yy(:,1) - x_max; % The value that we want to be zero

NO.

  value = yy(1,1) - x_max; % The value that we want to be zero

I'm trying to get the integration of ode45 to stop when the first column (which represents position)

The y that is passed into your ODEsystem function will not have multiple columns. It will have exactly one column, and that column will have as many rows as you have boundary conditions.

Your initial conditions are [0,0,10e-3,0] which is a vector of length 4, so the y passed into your function will be 4 x 1.

It does not matter whether you give [0,0,10e-3,0] or [0;0;10e-3;0] as your initial conditions: ode45() will always use (:) to reshape the initial conditions, and all your function will see is an input column vector.

You should be indexing yy(1,1) or just yy(1) instead of yy(:,1) because yy(:,1) will be the entire vector.