受到《趣题:三角形两顶点在直线上滑动时第三点的轨迹》(趣题:三角形两顶点在直线上滑动时第三点的轨迹 | 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点的轨迹不存在对称轴,而是一个中心对称图形。
等边三角形等腰三角形
锐角不等腰三角形
直角三角形
钝角三角形