While loop stops short of meeting the condition

5 次查看（过去 30 天）
Milos Krsmanovic
Milos Krsmanovic 2022-10-12
评论： Milos Krsmanovic 2022-10-19
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.
  2 个评论
Ghazwan
Ghazwan 2022-10-12
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
Milos Krsmanovic 2022-10-19
Thank you for your reply @Ghazwan. I'm still trying to figure out what the issue is, nevertheless I do appreciate your input.

请先登录，再进行评论。

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

采纳的回答

Mathieu NOE
Mathieu NOE 2022-10-12
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
  1 个评论
Milos Krsmanovic
Milos Krsmanovic 2022-10-19
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Arithmetic Operations 的更多信息

标签

产品


版本

R2019b

Community Treasure Hunt

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

Start Hunting!

Translated by