While loop stops short of meeting the condition

I have a loop function which is supposed to find the x=0 intercept of two trigonometric functions x=x(t). The two functions alter depending on whether the x is greater or lower than 0. Although I'm using the same while loop for both x(t)>0 and x(t)<0 domains, the while loop unexpectedly stops short of reaching x(t)=0 for only x(t)<0 intervals of the function. That is, only one half of the same while loop works properly.

Please forgive me for not uploading the complete code - I have half a dozen .m files, and also there is an issue of intellectual property (I'm not the owner, just a collaborator). I will try to explain the best I can what the problem is.

The subfunction and the while loop are as follows:

function[x] = get_x(f_dn, X, phi, wndn, wddn, zetadn, ttemp)
    if f_dn == 'u'
        x = X*cos(wndn*ttemp-phi);
    elseif f_dn == 'd'
        x = X*exp(-zetadn*wndn*ttemp)*cos(wddn*ttemp-phi);
    else fprintf('Incorrect f_up/f_dn entry (should be U or D).');
    end
end

ttemp = 1.05*t; % t is the lower limit of the interval, where x(t)=0,...
% already provided as an argument
xapproach = 0.001;

if rem(m,2) == 0 % even interval
    xup = xapproach;
    xdn = -xapproach;
else
    xup = xapproach;
    xdn = -xapproach;    
end

while (abs(x) > xapproach)
    if rem(m,2) == 0 % even interval
        if x > xup
            ttemp = ttemp + abs(x)*ttemp/m;
        elseif x < xdn
            ttemp = ttemp - abs(x)*ttemp/m;
        end
    elseif rem(m,2) ~= 0 % odd interval
        if x < xup
            ttemp = ttemp + abs(x)*ttemp/m;
        elseif x > xdn
            ttemp = ttemp - abs(x)*ttemp/m;
        end
    else
        disp('Some error.');
    end
    x = get_x(f_dn, X, phi, wndn, wddn, zetadn, ttemp)
end

Arguments such as X, phi, etc. are inherited from other functions, defined in other .m files.

The issue is visible from the plot:

As you can tell, each odd (red) interval starting from the third one onwards ends prematurely, which in turn affects the following even (blue) interval's start. At the same time, even (blue) intervals complete just fine, and I am able to find the x(t)=0 upper limits of those intervals just fine. But both of these are the part of the same while loop. I'm really confused as to why is this happening.

Let us look at a typical command window print of x for the odd (red) interval:

x = -0.1280
x = -0.2757
x = -0.5566
x = -0.8112
x = -0.1367
x = 0.0759
x = -0.0463
x = 0.0270
x = -0.0162
x = 0.0096
x = -0.0057
x = 0.0034
x = -0.0020
x = 0.0012

So from what I can tell my function starts right, it sets a point on the odd (red) interval so that x(t)<0, then starts moving the point right along the red line. If it overshoots, it corrects and moves the point to the left by an increasingly smaller increment. However, it ultimately suddenly stops at x = 0.0012, where the condition abs(x) > xapproach is still ans = 1. Increasing the precision of xapproach produces the same result. As noted, the very same while loop finds the upper limit x(t)=0 of the blue interval just fine, as you can alo see from the plot.

Any ideas on what I can investigate in more detail to understand what the issue is? Quite a few similar questions here had an issue of long precision, but I do not see that being an issue here (again, the same while loop, etc.) Thank you in advance.

[Ghazwan]

Given the provided information, it looks like it is an issue of rounding the double numbers when it is approaching zero. A possible solution is that you equate the outputs when reaching these points. If not, increasing the number of calculating points would reduce the gap.

[Milos Krsmanovic]

Thank you for your reply @Ghazwan. I'm still trying to figure out what the issue is, nevertheless I do appreciate your input.

Answer 1

[Mathieu NOE]

hello @Milos Krsmanovic

why not using this code ? as you see the zero (or any other threshold value) points are obtained with linear interpolation, so better accuracy as using nearest data point.

hope it helps

% dummy data
n = 100;
x= 25*(0:n-1)/n;
y1 = -cos(x).*exp(-0.1*x);

threshold = 0; % "zero" crossing or any other value here

[t0_pos1,s0_pos1,t0_neg1,s0_neg1]= crossing_V7(y1,x,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% ind => time index (samples)
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"


figure(1)
plot(x,y1,'b.-',t0_pos1,s0_pos1,'*r',t0_neg1,s0_neg1,'*g','linewidth',2,'markersize',12);grid on
legend('signal','signal positive slope crossing points','signal negative slope crossing points');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%


% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if ~isempty(ind)

    if strcmp(imeth,'linear')
        % linear interpolation of crossing
        for ii=1:length(t0)
            %if abs(S(ind(ii))) >= eps(S(ind(ii)))    % MATLAB V7 et +
            if abs(S(ind(ii))) >= eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)

                % interpolate only when data point is not already zero
                NUM = (t(ind(ii)+1) - t(ind(ii)));
                DEN = (S(ind(ii)+1) - S(ind(ii)));
                slope =  NUM / DEN;
                slope_sign(ii) = sign(slope);
                t0(ii) = t0(ii) - S(ind(ii)) * slope;
                s0(ii) = level;
            end
        end
    end

    % extract the positive slope crossing points 
    ind_pos = find(sign(slope_sign)>0);
    t0_pos = t0(ind_pos);
    s0_pos = s0(ind_pos);

    % extract the negative slope crossing points 
    ind_neg = find(sign(slope_sign)<0);
    t0_neg = t0(ind_neg);
    s0_neg = s0(ind_neg);

else

    % empty output
    ind_pos = [];
    t0_pos = [];
    s0_pos = [];
    ind_neg = [];
    t0_neg = [];
    s0_neg = [];
    
end

end

[Milos Krsmanovic]

Thank you for the reply and for the code sample @Mathieu NOE. I'm still very interested in finding out what the issue with my code is, exactly, as I need to reuse the same algorithm elsewhere, where your solution would not apply. I will mark your answer as a solution so that others can refer to it in the future.