eventfunction for stopping solver when output becomes complex doesn't work

Question

0 个投票

I am trying to solve the following DE;

function [v] = DE(t, y, Vb0, lambda, mb, rhog0, Cw)
% vaste parameters
LoadConstants
% tussenberekeningen
% krachten
Fa = rhol0*Vb0*g;
Fzk = lambda*y*g;
Fzb = (rhog0*Vb0+mb)*g;
% hoogteafhankelijke parameters
Vb = Vb0*(1-y*L/Tl0).^(1-(g*Ml/(R*L)));
rhol = rhol0*(1-y*L/Tl0).^((g*Ml/(R*L))-1);  
rb = ((3*Vb)/(4*pi)).^(1/3);
Ab = pi*rb.^2;
% snelheid
v = sqrt(2/3)*sqrt((Fa-Fzb-Fzk)./(Ab*Cw.*rhol));
end

I am doing it using the following function:

function [t, y] = SolveDE(Vb0, lambda, mb, rhog0, Cw)
tspan = [0 5000];
y0 = 0;
options = odeset('Events', @stop);
[t, y] = ode45(@(t, y) DE(t,y, Vb0, lambda, mb, rhog0, Cw),...
    tspan, y0, options);
end

with the stop-function defined as:

function [value, isterminal, direction] = stop(t, y)
value = isreal(y) && y <= 11000;
isterminal = 1;
direction = 0;
if imag(y) ~= 0
    disp(['STOP: y is imaginair geworden bij t = ', num2str(t)]);
elseif y > 11000
    disp(['STOP: y is groter dan 11km geworden bij t = ', num2str(t)]);
end
end

the eventfunction works perfect when is it y going above 11km that stops the solver, however when y becoming complex happens first something weird happens. I get the message that the stopfunction detected y becoming complex but is still continues for a bit resulting in a complex outputvector y. The values are in the form x + 0.0000000000000i. Can i fix this?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Follow Question

Answer 1

Walter Roberson 2025-3-31

编辑：Walter Roberson 2025-3-31

在 MATLAB Online 中打开

0 个投票

Try

value = [real(y) - 11000;
         imag(y)]
isterminal = [1;
              1];
direction = [1;
             0];

Note that integration will not stop as soon as it encounters an event: instead, once it encounters an event, it searches around trying to find the exact boundary so that it can stop at the boundary.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Arne 2025-4-1

This gives the same results

请先登录，再进行评论。

Answer 2

Torsten 2025-3-31

编辑：Torsten 2025-3-31

在 MATLAB Online 中打开

0 个投票

Try

options = odeset('Events', @(t,y)stop(t, y, Vb0, lambda, mb, rhog0, Cw));
[t, y] = ode45(@(t, y) DE(t,y, Vb0, lambda, mb, rhog0, Cw),...
    tspan, y0, options);
function [value, isterminal, direction] = stop(t, y, Vb0, lambda, mb, rhog0, Cw)
  % vaste parameters
  LoadConstants
  % tussenberekeningen
  % krachten
  Fa = rhol0*Vb0*g;
  Fzk = lambda*y*g;
  Fzb = (rhog0*Vb0+mb)*g;
  % hoogteafhankelijke parameters
  Vb = Vb0*(1-y*L/Tl0).^(1-(g*Ml/(R*L)));
  rhol = rhol0*(1-y*L/Tl0).^((g*Ml/(R*L))-1);  
  rb = ((3*Vb)/(4*pi)).^(1/3);
  Ab = pi*rb.^2;
  tol = 1e-3;
  value = [1-y*L/Tl0-tol;Fa-Fzb-Fzk-tol;y-11000]
  isterminal = [1;1;1];
  direction = [0;0;0];
end

If this does not work, we must have executable code for testing.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Torsten 2025-4-2

在 MATLAB Online 中打开

Vb0 = 10000;

lambda = 0.125;

mb = (20*10^-3)*(36*pi*Vb0^2)^(1/3);

rhog0 = 0.169;

Cw = 0.285;

[t,y] = SolveDE(Vb0, lambda, mb, rhog0, Cw)

t = 49×1

1.0e+00 * 0 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0001 0.0001 0.0001

<mw-icon class=""></mw-icon>

y = 49×1

0 0.0001 0.0001 0.0002 0.0002 0.0005 0.0007 0.0010 0.0012 0.0025

<mw-icon class=""></mw-icon>

plot(t,y)

function [t, y] = SolveDE(Vb0, lambda, mb, rhog0, Cw)

tspan = [0 5000];

y0 = 0;

options = odeset('Events', @(t,y)stop(t, y, Vb0, lambda, mb, rhog0, Cw));

[t, y] = ode45(@(t, y) DE(t,y, Vb0, lambda, mb, rhog0, Cw),...

tspan, y0, options);

end

function [v] = DE(t, y, Vb0, lambda, mb, rhog0, Cw)

% vaste parameters

rhol0 = 1.225;

L = 0.0065;

Tl0 =288.15;

R = 8.314;

g = 9.81;

Ml = 0.02896;

Pl0 =101325;

% tussenberekeningen

% krachten

Fa = rhol0*Vb0*g;

Fzk = lambda*y*g;

Fzb = (rhog0*Vb0+mb)*g;

% hoogteafhankelijke parameters

Vb = Vb0*(1-y*L/Tl0).^(1-(g*Ml/(R*L)));

rhol = rhol0*(1-y*L/Tl0).^((g*Ml/(R*L))-1);

rb = ((3*Vb)/(4*pi)).^(1/3);

Ab = pi*rb.^2;

% snelheid

v = sqrt(2/3)*sqrt(abs(Fa-Fzb-Fzk)./(Ab*Cw.*rhol));

end

function [value, isterminal, direction] = stop(t, y, Vb0, lambda, mb, rhog0, Cw)

% vaste parameters

rhol0 = 1.225;

L = 0.0065;

Tl0 =288.15;

R = 8.314;

g = 9.81;

Ml = 0.02896;

Pl0 =101325;

% tussenberekeningen

% krachten

Fa = rhol0*Vb0*g;

Fzk = lambda*y*g;

Fzb = (rhog0*Vb0+mb)*g;

value = [1-y*L/Tl0;Fa-Fzb-Fzk;y-11000];

isterminal = [1;1;1];

direction = [0;0;0];

end

Arne 2025-4-2

Yes for big values of Vb0 the output of the DE never bevomes complex. However i managed to fix the problem by adding Reltol and Abstol arguments with high requirements in the odeset function.

Torsten 2025-4-2

Yes for big values of Vb0 the output of the DE never bevomes complex.

I tested that it also works for small values of Vb0. But if you already made the code work, it's fine.

请先登录，再进行评论。

eventfunction for stopping solver when output becomes complex doesn't work

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

类别

标签

Community Treasure Hunt

eventfunction for stopping solver when output becomes complex doesn't work

0 个评论 显示 -2更早的评论 隐藏 -2更早的评论

回答（2 个）

1 个评论 显示 -1更早的评论 隐藏 -1更早的评论

4 个评论 显示 2更早的评论 隐藏 2更早的评论

类别

标签

另请参阅

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论