三角形两顶点在椭圆上滑动时第三点的轨迹

        受到《趣题：三角形两顶点在直线上滑动时第三点的轨迹》(趣题：三角形两顶点在直线上滑动时第三点的轨迹 | Matrix67: The Aha Moments)的启发，决定研究一下三角形两顶点在椭圆上滑动时第三点的轨迹。△ABC和△ABD的A、B两点在椭圆x^2/a^2+y^2/b^2=1(a>b)上滑动，C和D两点关于线段AB对称，其中C点更靠近椭圆中心。C点不一定在椭圆内部，D点也不一定在椭圆外部。C、D两点轨迹的参数方程比较容易得出，先假设A点的坐标(Ax,Ay)，Ax=acosθ，Ay=bsinθ，然后求以A为圆心、AB长度为半径的圆与椭圆的交点(求解四次方程)，即得B点坐标。
    ■|AB|较小时，四次方程只有2个实根，可以得到2个B点；    
    ■|AB|较大时，四次方程可能出现4个实根，还可能有重根。(比如|AB|=2b，A位于椭圆短轴端点，B可以位于椭圆左侧、右侧以及另一个短轴端点，Bx=0就是二重根)
    ■|AB|=2a，则A、B分别位于2个长轴端点，三角形无法滑动。
    ■|AB|>2a，则A、B无法同时位于椭圆上。

       第一次解四次方程后，从实根中随意选出一个作为第一个B点的横坐标。那么如何从之后解出的2个或4个实根中选出合适的根作为下一个B点的横坐标呢？最直接的想法是要尽可能地靠近前一个B点。四次方程有2个实根时，在(Bx(1),±b*sqrt(1-Bx(1)^2/a^2))，(Bx(2),±b*sqrt(1-Bx(2)^2/a^2))这4个点中，只有2个点与A点的距离是|AB|,用
abs((NewBx-Ax)^2+(b*sqrt(1-NewBx^2/a^2)-Ay)^2-AB^2)<1e-8
这一条件就能选出2个正确的B点，在这2个正确的点中找出与前一个B点距离最近的，就是所需的下一个B点。不直接写成
(NewBx-Ax)^2+(b*sqrt(1-NewBx^2/a^2)-Ay)^2==AB^2
是因为浮点运算的误差可能导致并不严格相等(误差在1e-15量级)。四次方程有4个实根时，共有(Bx(1),±By(1)),(Bx(2),±By(2)),(Bx(3),±By(3)),(Bx(4),±By(4))8个可能的点，先从前4个可能点中选出一个距离最小的，再从后4个可能点中选出一个距离最小的，最后从这两个距离“最小”的点中选出一个距离更小的。

       提供一个有4个实根的数值例子,a=2,b=1,|AB|=1.9,
A(2cos(1.31),sin(1.31))=A(0.515700,0.966185),B1(-1.369414,0.728818),B2(1.814578,-0.420507),B3(1.436365,-0.695855),B4(0.868872,-0.900703)。

MATLAB代码

clc;close all;
tic
global a b AB
a=2; b=1; 
AB=1.90;BC=0.7;AC=1.3;
cos_BAC=(AB^2+AC^2-BC^2)/(2*AB*AC); %余弦定理
sin_BAC=sqrt(1-cos_BAC^2);
N=629;
Bpoint=zeros(2,N);  Cpoint=zeros(2,N);  Dpoint=zeros(2,N);

int1=1;
QA=int1/100;
Ax=a*cos(QA);
Ay=b*sin(QA);
T=Ax^2+Ay^2+b^2-AB^2;
root4=roots([(b^2/a^2 - 1)^2, 4*Ax*(b^2/a^2 - 1), 4*Ax^2 - 2*T*(b^2/a^2 - 1)...
    + (4*Ay^2*b^2)/a^2, -4*Ax*T,- 4*Ay^2*b^2 + T^2]);
Bx=zeros(2,1);
[Bx(1),Bx(2)]=SelectRealRoot(root4);
Bpoint(1,int1)=Bx(1);
if abs((Bx(1)-Ax)^2+(b*sqrt(1-Bx(1)^2/a^2)-Ay)^2-AB^2)<1e-8
    Bpoint(2,int1)=b*sqrt(1-Bx(1)^2/a^2); 
else
    Bpoint(2,int1)=-b*sqrt(1-Bx(1)^2/a^2);
end

ABvector=[Bpoint(1,int1)-Ax,Bpoint(2,int1)-Ay];
ACvector=ABvector*AC/AB*[cos_BAC,sin_BAC;-sin_BAC,cos_BAC]; %实现向量的旋转
ADvector=ABvector*AC/AB*[cos_BAC,-sin_BAC;sin_BAC,cos_BAC];
if ((Ax+ACvector(1))^2+(Ay+ACvector(2))^2)>((Ax+ADvector(1))^2+(Ay+ADvector(2))^2)
    %假定C点更靠近椭圆中心，这与选择第一个B点的坐标一样，具有一定的随机性。
    TempVector=ACvector;
    ACvector=ADvector;
    ADvector=TempVector;
end
Cpoint(:,int1)=[Ax+ACvector(1);Ay+ACvector(2)];
Dpoint(:,int1)=[Ax+ADvector(1);Ay+ADvector(2)];

for int1=2:N
     int1 %显示int1监测进度
% QA=0.01:0.01:2*pi
QA=int1/100;
Ax=a*cos(QA);
Ay=b*sin(QA) ;
T=Ax^2+Ay^2+b^2-AB^2;
root4=roots([(b^2/a^2 - 1)^2, 4*Ax*(b^2/a^2 - 1), 4*Ax^2 - 2*T*(b^2/a^2 - 1)...
           + (4*Ay^2*b^2)/a^2, -4*Ax*T, - 4*Ay^2*b^2 + T^2]);

if abs(imag(root4(1)))<1e-8 && abs(imag(root4(2)))<1e-8  && ...
   abs(imag(root4(3)))<1e-8 && abs(imag(root4(4)))<1e-8  %如果有4个实根
        root4
   [Bx(1),~]=...
   DetermineBpoint(root4(1),root4(2),Bpoint(1,int1-1),Bpoint(2,int1-1),Ax,Ay); 
   [Bx(2),~]=...
   DetermineBpoint(root4(3),root4(4),Bpoint(1,int1-1),Bpoint(2,int1-1),Ax,Ay); 
else %如果有2个实根
   [Bx(1),Bx(2)]=SelectRealRoot(root4);
end
   [Bpoint(1,int1),Bpoint(2,int1)]=...
    DetermineBpoint(Bx(1),Bx(2),Bpoint(1,int1-1),Bpoint(2,int1-1),Ax,Ay); 

    ABvector=[Bpoint(1,int1)-Ax,Bpoint(2,int1)-Ay];
    ACvector=ABvector*AC/AB*[cos_BAC,sin_BAC;-sin_BAC,cos_BAC]; %实现向量的旋转
    ADvector=ABvector*AC/AB*[cos_BAC,-sin_BAC;sin_BAC,cos_BAC];
    if ((Ax+ACvector(1))^2+(Ay+ACvector(2))^2)>((Ax+ADvector(1))^2+(Ay+ADvector(2))^2)
        %因为C点更靠近椭圆中心
        TempVector=ACvector;
        ACvector=ADvector;
        ADvector=TempVector;
    end
    Cpoint(:,int1)=[Ax+ACvector(1);Ay+ACvector(2)];
    Dpoint(:,int1)=[Ax+ADvector(1);Ay+ADvector(2)];
end

LWidth=1;
for int1=1:30:N
    QA=int1/100;
    Ax=a*cos(QA);
    Ay=b*sin(QA); 
    ColorArray=[mod(5*int1,N)/N,mod(10*int1,N)/N,mod(15*int1,N)/N];
    line([Ax,Bpoint(1,int1)],[Ay,Bpoint(2,int1)],'color',ColorArray,'LineWidth',LWidth); hold on
    line([Ax,Cpoint(1,int1)],[Ay,Cpoint(2,int1)],'color',ColorArray,'LineWidth',LWidth); 
    line([Cpoint(1,int1),Bpoint(1,int1)],[Cpoint(2,int1),Bpoint(2,int1)],'color',ColorArray,'LineWidth',LWidth); 
    line([Ax,Dpoint(1,int1)],[Ay,Dpoint(2,int1)],'color',ColorArray,'LineWidth',LWidth); 
    line([Dpoint(1,int1),Bpoint(1,int1)],[Dpoint(2,int1),Bpoint(2,int1)],'color',ColorArray,'LineWidth',LWidth); 
end

plotEllipse=plot(a*cos(0:0.01:2*pi),b*sin(0:0.01:2*pi),'LineWidth',2,'color',[0,0,0]);hold on
plotC=plot(Cpoint(1,:),Cpoint(2,:),'LineWidth',2,'color',[1,0,0]);
plotD=plot(Dpoint(1,:),Dpoint(2,:),'LineWidth',2,'color',[0,0,1]);
legend([plotEllipse,plotC,plotD],{'椭圆x^2/a^2+y^2/b^2=1','C点轨迹','D点轨迹'})
title(strcat('a=',num2str(a),', b=',num2str(b),'; AB=',num2str(AB),...
    ', BC=',num2str(BC),', AC=',num2str(AC),'; 作者:磁悬浮青蛙呱呱呱'))
h=max([AB,AC,BC])+0.1;
axis([-(a+h),a+h,-(b+h),b+h]);axis equal
hold off
toc 

function [Bx1,Bx2]=SelectRealRoot(root4)
int3=1;
Bx=zeros(2,1);
    for int2=1:4   %从4个根中挑选出2个实根
        if abs(imag(root4(int2)))<1e-8
           Bx(int3)= root4(int2);
           int3=int3+1;
        end
    end
Bx1=Bx(1);
Bx2=Bx(2);
end

function [Bx,By]=DetermineBpoint(NewBx1,NewBx2,PreviousBx,PreviousBy,Ax,Ay)
global a b AB
BTemp=zeros(2,2);
int2=1;

    if  abs((NewBx1-Ax)^2+(b*sqrt(1-NewBx1^2/a^2)-Ay)^2-AB^2)<1e-8
        BTemp(:,int2)=[NewBx1;b*sqrt(1-NewBx1^2/a^2)];
        int2=int2+1;
    end
    if  abs((NewBx1-Ax)^2+(-b*sqrt(1-NewBx1^2/a^2)-Ay)^2-AB^2)<1e-8
        BTemp(:,int2)=[NewBx1;-b*sqrt(1-NewBx1^2/a^2)];
        int2=int2+1;
    end    
    if  abs((NewBx2-Ax)^2+(b*sqrt(1-NewBx2^2/a^2)-Ay)^2-AB^2)<1e-8
        BTemp(:,int2)=[NewBx2;b*sqrt(1-NewBx2^2/a^2)];
        int2=int2+1;
    end    
    if  abs((NewBx2-Ax)^2+(-b*sqrt(1-NewBx2^2/a^2)-Ay)^2-AB^2)<1e-8
        BTemp(:,int2)=[NewBx2;-b*sqrt(1-NewBx2^2/a^2)];
    end
    
    if norm(BTemp(:,1)-[PreviousBx;PreviousBy])<norm(BTemp(:,2)-[PreviousBx;PreviousBy])
        Bx=BTemp(1,1); By=BTemp(2,1);  return;
    else
        Bx=BTemp(1,2); By=BTemp(2,2);
    end
end

       下面给出的C、D两点的轨迹图，主要是0<|AB|<2b情况下的，其它情形可自行研究。观察下面这些轨迹图，当|AC|=|BC|时，C点和D点的轨迹存在两根相互垂直的对称轴。|AC|≠|BC|时，C点和D点的轨迹不存在对称轴，而是一个中心对称图形。

等边三角形

等腰三角形

锐角不等腰三角形

直角三角形

钝角三角形