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

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?

Answer 1

[Walter Roberson]

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.

[Arne]

This gives the same results

Answer 2

[Torsten]

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.

[Arne]

This is the LoadConstants script:

rhol0 = 1.225;
L = 0.0065;
Tl0 =288.15;
R = 8.314;
g = 9.81;
Ml = 0.02896;
Pl0 =101325; 
Cw = 0.285; 

Vb0 can be any value between 10 and 10 000. Lambda can be 0.007, 0.0173 or 0.125. rhog0 is 0.169 and mb can be calculated via mb = (20*10^-3)*(36*pi*Vb0^2)^(1/3).

[Torsten]

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>

<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>

<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]

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]

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.