ode45 event trigger not responding
6 views (last 30 days)
Show older comments
Nader Mohamed
on 31 Oct 2021
Commented: Walter Roberson
on 31 Oct 2021
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
0 Comments
Accepted Answer
Walter Roberson
on 31 Oct 2021
Edited: Walter Roberson
on 31 Oct 2021
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.
3 Comments
Walter Roberson
on 31 Oct 2021
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
on 31 Oct 2021
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.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!