ode45 event trigger not responding

2021 10 31
1 个回答
回答已采纳
更新时间：2021 10 31
8 次查看（30 天）

0 个投票

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

请先登录，再回答此问题。

 采纳的回答

Walter Roberson
Walter Roberson 2021-10-31
编辑：Walter Roberson 2021-10-31

0 个投票

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 个评论
显示 1更早的评论 隐藏 1更早的评论

Nader Mohamed
Nader Mohamed 2021-10-31
@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
Walter Roberson 2021-10-31
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
Walter Roberson 2021-10-31
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Environment and Settings 的更多信息

产品

版本

R2018b

标签

Community Treasure Hunt

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

Start Hunting!

Translated by