曲线的旋转、曲线的交集以及for循环的应用【浮点数比较问题】

Question

问题场景定义如下problem statement (Click here)

提供单次迭代中发生的信息here in this link

我应该得到一个小圆弧在两端接触圆（曲线 1 围绕点 A 旋转圆然后在点 B 旋转接触圆所以圆弧的两端接触圆圈）。当我手动执行此操作时，我得到了结果，通过使用 polyxpoly() 找到交点，将圆弧放在交点上并继续旋转它直到它接触到圆。 但是当我尝试使用 for 循环自动执行许多弧的过程时，我没有得到想要的结果。

由于发生了圆弧旋转，我给出了较小的旋转步长，但我得到了不同的结果。 有时for循环工作正常并给出期望的结果，有时在某些迭代中给出异常结果。我无法找出我的方法的问题。

结果如下click here to see the result

这些是显示预期结果与非预期结果之间差异的附加详细信息 Desired Result Vs Undesired Result for theta=500 & i=2 ,3 Desired Result Vs Undesired Result for theta=1800 & i=4

可以看出 theta = 500 & i = 2 弧线不旋转，因为它不旋转，所以 theta = 500 & i = 3 的结果相同但是当 theta = 1800 它给出不同的结果但是 在迭代 4 失败（对于 theta = 1800 i = 1,2,3 给出期望的结果但在 i = 4 失败） .检查结果图像。

参数theta可以在rotation()

我使用的代码如下，但 For-loop 根据 theta

对某些迭代给出异常结果

clc,clear
r = 10;                     % arc length
R = 55;                     % radius of a circle
aa = 60*pi/180;              % arc angle
ap = 0*pi/180;             % arc position angle
k=0;

t = linspace(0,pi/2);
[x,y] = pol2cart(t,R);      % circle data
t1 = linspace(0,aa)-aa/2+ap;
[x1,y1] = pol2cart(t1,r); % arc data
% Intersection of x axis (y=0) and quarter circle
x_axis_range=1:1000;
y_val_on_x_axis=zeros(1,length(x_axis_range));
[x1rc,y1rc]=polyxpoly(x_axis_range,y_val_on_x_axis,x,y);
xc_s(1)=x1rc;
yc_s(1)=y1rc;
% x1rc=55;
% y1rc=0;
for z=1:3
[x1_p,y1_p]= placement(x,y,x1,y1,x1rc,y1rc);
[x1_rotated,y1_rotated,xc,yc]=rotate(x,y,x1_p,y1_p,x1rc,y1rc);
%[x1_rotated,y1_rotated,xc,yc]=rotate(x,y,x1,y1,x1rc,y1rc);
x1rc=xc;
y1rc=yc;
end
% havent written the plot command

Placement()- placement()的作用是将新圆弧移动到原圆弧与圆的交点

function [x1_p,y1_p] = placement(x,y,x1,y1,x1rc,y1rc)

xcen=x1rc;
ycen=y1rc;

% shifting arc-lower-end to (t(1),0) or (
delx=xcen-x1(1); % Finding the x difference between arc-lower-end x-coordinate & x(1)[circle]
dely=ycen-y1(1); % Finding the y difference between arc-lower-end y-coordinate & y(1) [circle]

x1_p=x1+delx;
y1_p=y1+dely;
end

旋转function():

 function [x1_rotated,y1_rotated,xc,yc] = rotate(x,y,x1_p,y1_p,x1rc,y1rc)
    % [x1_p,y1_p]= placement(x,y,x1,y1,x1rc,y1rc);
    theta =linspace(0,pi,500);
    i=1;
    xc=[];
    yc=[];
    xcen=x1rc;
    ycen=y1rc;
    while isempty(xc)||xc==xcen && isempty(yc)||yc==ycen
    % create a matrix of these points, which will be useful in future calculations
    v = [x1_p;y1_p];
    % choose a point which will be the center of rotation
    x_center = xcen;
    y_center = ycen;
    % create a matrix which will be used later in calculations
    center = repmat([x_center; y_center], 1, length(x1_p));
    % define a theta degree counter-clockwise rotation matrix
    R = [cos(theta(i)) -sin(theta(i)); sin(theta(i)) cos(theta(i))];
    % do the rotation...
    s = v - center;     % shift points in the plane so that the center of rotation is at the origin
    so = R*s;           % apply the rotation about the origin
    vo = so + center;   % shift again so the origin goes back to the desired center of rotation
    % this can be done in one line as:
    % vo = R*(v - center) + center
    % pick out the vectors of rotated x- and y-data
    x1_rotated = vo(1,:);
    y1_rotated = vo(2,:);
    [xc,yc] = polyxpoly(x1_rotated,y1_rotated,x,y);
    [xc1,yc1] = polyxpoly(x1_p,y1_p,x,y);
    if length(xc)==2
    xc=xc(2);
    yc=yc(2);
    end
    i=i+1;
    end
    end

我尝试如下修改“if 语句”，但似乎在一定的迭代后它不起作用。 if语句新修改如下

   if length(xc)==2
   ubx=xc(1); lbx=xc(2);
   uby=yc(2);lby=yc(1);
   xc=xc(2);
   yc=yc(2)
    
  ub=[ubx uby];
  lb=[lbx-4 lby];
  end

然后我尝试 return ub 和 lb 并收到以下错误消息 输出参数“ub”（可能还有其他）不是在 for loop 的 第五次迭代 之后调用“旋转”时分配。。 if 语句 对 1,2,3,4 迭代工作正常

Answer 1

这是一个 floating point comparison problem。长话短说，3*(4/3 -1) == 1 和 sin(pi)==0 return 都是假的 (Avoiding Common Problems with Floating-Point Arithmetic).

将您的 while 条件更改为如下内容：

while (isempty(xc)||abs(xc-xcen)<=10*eps(xcen)) && (isempty(yc)||abs(yc-ycen)<=10*eps(ycen))

顺便说一句，您假设当 polyxpoly returns 2 点时，第二点是您的结果，这是不正确的。尝试找到距离圆弧起点最远的点：

if length(xc)>1
    [~, i] = max((xc-xcen).^2+(yc-ycen).^2);
    xc=xc(i);
    yc=yc(i);
end

曲线的旋转、曲线的交集以及for循环的应用【浮点数比较问题】

Rotation of curves , intersection of curve and for -loop application [Floating point comparison problem]]

math

matlab

geometry

for-loop

rotation