两个椭圆的公切线求法（Matlab） |
您所在的位置：网站首页 › 直线与圆相切的切点怎么求 › 两个椭圆的公切线求法（Matlab） |
两个椭圆的公切线求法（Matlab）
|
文章目录
问题描述问题求解椭圆生成求解过程切线方程切点坐标特殊情况处理
测试程序代码
问题描述
已知两个椭圆e1、e2的一般方程，求与两个椭圆相切的所有公切线方程和切点坐标（数值解）。 椭圆一般方程定义如下： A x 2 + B x y + C y 2 + D x + E y + F = 0 Ax^2 + Bxy + Cy^2 +Dx + Ey + F= 0 Ax2+Bxy+Cy2+Dx+Ey+F=0 问题求解使用Matlab 2015a进行求解。 两个椭圆的公切线存在以下情况： 两个椭圆重合，无数条公切线，这个不做讨论其中一个椭圆被另一个椭圆完全包含，不存在公切线两个椭圆内切，只有一个切点，存在1条公切线两个椭圆内切，只有两个切点，存在2条公切线两个椭圆相交，只有两个交点，存在2条公切线两个椭圆相交，只有三个交点，存在3条公切线两个椭圆相交，只有四个交点，存在4条公切线两个椭圆外切，存在3条公切线两个椭圆相离，存在4条公切线 椭圆生成首先需生成两个椭圆，通过定义椭圆的半长轴 a a a，半短轴 b b b，中心坐标 ( x 0 , y 0 ) (x_0,y_0) (x0,y0)，长轴倾角（逆时针为正） θ \theta θ，根据维基百科Ellipse，得到椭圆的一般方程的系数为： A = a 2 s i n ( θ ) 2 + b 2 c o s ( θ ) 2 A = a^2sin(\theta)^2 + b^2cos(\theta)^2 A=a2sin(θ)2+b2cos(θ)2 B = 2 ( b 2 − a 2 ) ∗ s i n ( θ ) c o s ( θ ) B = 2(b^2-a^2)*sin(\theta)cos(\theta) B=2(b2−a2)∗sin(θ)cos(θ) C = a 2 c o s ( θ ) 2 + b 2 s i n ( θ ) 2 C = a^2cos(\theta)^2+b^2sin(\theta)^2 C=a2cos(θ)2+b2sin(θ)2 D = − 2 A x 0 − B y 0 D = -2Ax_0-By_0 D=−2Ax0−By0 E = − B x 0 − 2 C y 0 E = -Bx_0-2Cy_0 E=−Bx0−2Cy0 F = A x 0 2 + B x 0 y 0 + C y 0 2 − a 2 b 2 F = Ax_0^2+Bx_0y_0+Cy_0^2-a^2b^2 F=Ax02+Bx0y0+Cy02−a2b2 求解过程 切线方程令公切线方程为 y = k x + b y=kx+b y=kx+b 联立以下公式: A 1 x 2 + B 1 x y + C 1 y 2 + D 1 x + E 1 y + F 1 = 0 A_1x^2 + B_1xy + C_1y^2 +D_1x + E_1y + F_1= 0 A1x2+B1xy+C1y2+D1x+E1y+F1=0 A 2 x 2 + B 2 x y + C 2 y 2 + D 2 x + E 2 y + F 2 = 0 A_2x^2 + B_2xy + C_2y^2 +D_2x + E_2y + F_2= 0 A2x2+B2xy+C2y2+D2x+E2y+F2=0 y = k x + b y=kx+b y=kx+b 化简得 ( A 1 + B 1 k + C 1 k 2 ) x 2 + ( B 1 b + 2 C 1 k b + D 1 + E 1 k ) x + ( C 1 b 2 + E 1 b + F 1 ) = 0 (A_1+B_1k+C_1k^2)x^2+(B_1b+2C_1kb+D_1+E_1k)x+(C_1b^2+E_1b+F_1)=0 (A1+B1k+C1k2)x2+(B1b+2C1kb+D1+E1k)x+(C1b2+E1b+F1)=0 ( A 2 + B 2 k + C 2 k 2 ) x 2 + ( B 2 b + 2 C 2 k b + D 2 + E 2 k ) x + ( C 2 b 2 + E 2 b + F 2 ) = 0 (A_2+B_2k+C_2k^2)x^2+(B_2b+2C_2kb+D_2+E_2k)x+(C_2b^2+E_2b+F_2)=0 (A2+B2k+C2k2)x2+(B2b+2C2kb+D2+E2k)x+(C2b2+E2b+F2)=0 当以上两式同时具有重根时，表示直线 y = k x + b y=kx+b y=kx+b同时与两椭圆相切， 一元二次方程 a x 2 + b x + c = 0 ax^2+bx+c=0 ax2+bx+c=0的重根判别式 Δ = b 2 − 4 a c = 0 \Delta=b^2-4ac=0 Δ=b2−4ac=0，重根 x 1 , 2 = − b / ( 2 a ) x_{1,2}=-b/(2a) x1,2=−b/(2a) 联立两个判别式： ( B 1 b + 2 C 1 k b + D 1 + E 1 k ) 2 − 4 ( A 1 + B 1 k + C 1 k 2 ) ( C 1 b 2 + E 1 b + F 1 ) = 0 (B_1b+2C_1kb+D_1+E_1k)^2-4(A_1+B_1k+C_1k^2)(C_1b^2+E_1b+F_1) = 0 (B1b+2C1kb+D1+E1k)2−4(A1+B1k+C1k2)(C1b2+E1b+F1)=0 ( B 2 b + 2 C 2 k b + D 2 + E 2 k ) 2 − 4 ( A 2 + B 2 k + C 2 k 2 ) ( C 2 b 2 + E 2 b + F 2 ) = 0 (B_2b+2C_2kb+D_2+E_2k)^2-4(A_2+B_2k+C_2k^2)(C_2b^2+E_2b+F_2) = 0 (B2b+2C2kb+D2+E2k)2−4(A2+B2k+C2k2)(C2b2+E2b+F2)=0 采用matlab中的solve函数进行求解。 syms k b eq1 = (B1*b+2*C1*k*b+D1+E1*k)^2-4*(A1+B1*k+C1*k^2)*(C1*b^2+E1*b+F1); eq2 = (B2*b+2*C2*k*b+D2+E2*k)^2-4*(A2+B2*k+C2*k^2)*(C2*b^2+E2*b+F2); s=solve(eq1,eq2,'k','b'); k = double(s.k); b = double(s.b); 切点坐标求出 k k k和 b b b后，再利用重根公式计算切点的横坐标， aa1=A1+B1*k+C1*k.^2; bb1=B1*b+2*C1*k.*b+D1+E1*k; x1 = -bb1./(2*aa1); y1 = k.*x1+b; aa2=A2+B2*k+C2*k.^2; bb2=B2*b+2*C2*k.*b+D2+E2*k; x2 = -bb2./(2*aa2); y2 = k.*x2+b;然后代入 y = k x + b y=kx+b y=kx+b中，得出纵坐标。 特殊情况处理当切线斜率无穷大时， k k k不存在，因此需要作为特殊情况进行处理，令此时的直线方程为: x = t x=t x=t 代入两椭圆方程中, A 1 t 2 + B 1 t y + C 1 y 2 + D 1 t + E 1 y + F 1 = 0 A_1t^2 + B_1ty + C_1y^2 +D_1t + E_1y + F_1= 0 A1t2+B1ty+C1y2+D1t+E1y+F1=0 A 2 t 2 + B 2 t y + C 2 y 2 + D 2 t + E 2 y + F 2 = 0 A_2t^2 + B_2ty + C_2y^2 +D_2t + E_2y + F_2= 0 A2t2+B2ty+C2y2+D2t+E2y+F2=0 同理，当以上两式同时具有重根时，表示直线 x = t x=t x=t同时与两椭圆相切。根据重根条件，得, ( B 1 t + E 1 ) 2 − 4 C 1 ( A 1 t 2 + D 1 t + F 1 ) = 0 (B_1t+E_1)^2-4C_1(A_1t^2+D_1t+F_1) = 0 (B1t+E1)2−4C1(A1t2+D1t+F1)=0 ( B 2 t + E 2 ) 2 − 4 C 2 ( A 2 t 2 + D 2 t + F 2 ) = 0 (B_2t+E_2)^2-4C_2(A_2t^2+D_2t+F_2) = 0 (B2t+E2)2−4C2(A2t2+D2t+F2)=0 切线横坐标 t t t和切点坐标求解方式同样采用solve函数，以及重根公式计算， syms t eq1 = (B1*t+E1)^2-4*C1*(A1*t^2+D1*t+F1); eq2 = (B2*t+E2)^2-4*C2*(A2*t^2+D2*t+F2); s1=solve(eq1,'t'); s2=solve(eq2,'t'); b1 = double(s1); b2 = double(s2); kinf = 0; if( norm(b1-b2) == 0) || (norm(b1- flipud(b2)) == 0) sx1 = b1; sy1 = -(B1*b1+E1)/(2*C1); sx2 = b1; sy2 = -(B2*b1+E2)/(2*C2); kinf = 1; end 测试matlab求解程序和测试代码附在文后。 测试图片如下: FindCommonTangentofTwoEllipse.m FindCommonTangentofTwoEllipseTest.m function [] = FindCommonTangentofTwoEllipse(ellipse1, ellipse2) % Function: find common tangent line for two separate ellipses % Author: Yang Yang % Date: 2019-02-18 % Usage: input two ellipses paramters, [halfMajor, halfMinor, centerX, centerY, theta(deg)] %% generate ellipses % ellipse 1 [ A, B, C, D, E, F] = generateEllipse(ellipse1(1), ellipse1(2), ellipse1(3), ellipse1(4), ellipse1(5)); % check the ellipse equation coefficients by calcuating center xc=(B*E-2*C*D)/(4*A*C-B^2); yc=(B*D-2*A*E)/(4*A*C-B^2); center1 = [xc,yc]; eps1= [ A B C D E F]'; % ellipse 2 [ A, B, C, D, E, F] = generateEllipse(ellipse2(1), ellipse2(2), ellipse2(3), ellipse2(4), ellipse2(5)); % check the ellipse equation coefficients by calcuating center xc=(B*E-2*C*D)/(4*A*C-B^2); yc=(B*D-2*A*E)/(4*A*C-B^2); center2 = [xc,yc]; eps2= [ A B C D E F]'; %draw ellipse1 and ellipse2 figure; A=eps1; a=A(1); b=A(2); c=A(3); d=A(4); e=A(5); f=A(6); eq0=@(x,y) a*x^2 + b*x*y + c*y^2 +d*x + e*y + f; hold on; h=ezplot(eq0,[-250,250,-180,180]); % ellipse 1 set(h,'Color','r'); plot(center1(1), center1(2), '+', 'Color', 'r'); A=eps2; a=A(1); b=A(2); c=A(3); d=A(4); e=A(5); f=A(6); eq0=@(x,y) a*x^2 + b*x*y + c*y^2 +d*x + e*y + f; h=ezplot(eq0,[-250,250,-180,180]); % ellipse 2 set(h,'Color','g'); plot(center2(1), center2(2), '+', 'Color', 'g'); %% solve the common tangent line eps1 = eps1/eps1(6); % F normalization eps2 = eps2/eps2(6); % F normalization A1=eps1(1);B1=eps1(2);C1=eps1(3);D1=eps1(4);E1=eps1(5);F1=eps1(6); A2=eps2(1);B2=eps2(2);C2=eps2(3);D2=eps2(4);E2=eps2(5);F2=eps2(6); % check condition: slope of the line is infinite, i.e. k=inf % deal with this special case, let x=t. syms t eq1 = (B1*t+E1)^2-4*C1*(A1*t^2+D1*t+F1); eq2 = (B2*t+E2)^2-4*C2*(A2*t^2+D2*t+F2); s1=solve(eq1,'t'); s2=solve(eq2,'t'); b1 = double(s1); b2 = double(s2); kinf = 0; if( norm(b1-b2) == 0) || (norm(b1- flipud(b2)) == 0) sx1 = b1; sy1 = -(B1*b1+E1)/(2*C1); sx2 = b1; sy2 = -(B2*b1+E2)/(2*C2); kinf = 1; end % deal with general case syms k b % eq1 = (B1*b+2*C1*k*b+D1+E1*k)^2-4*(A1+B1*k+C1*k^2)*(C1*b^2+E1*b+F1); % eq2 = (B2*b+2*C2*k*b+D2+E2*k)^2-4*(A2+B2*k+C2*k^2)*(C2*b^2+E2*b+F2); % work better when expand the above equations, the expansion process is % given below % syms k b A1 B1 C1 D1 E1 F1 % ss = (B1*b+2*C1*k*b+D1+E1*k)^2-4*(A1+B1*k+C1*k^2)*(C1*b^2+E1*b+F1); % eq = expand(ss); eq1 = B1^2*b^2 + 2*B1*D1*b - 2*B1*E1*b*k - 4*F1*B1*k + D1^2 + 2*D1*E1*k + ... 4*C1*D1*b*k + E1^2*k^2 - 4*A1*E1*b - 4*A1*C1*b^2 - 4*C1*F1*k^2 - 4*A1*F1; eq2 = B2^2*b^2 + 2*B2*D2*b - 2*B2*E2*b*k - 4*F2*B2*k + D2^2 + 2*D2*E2*k + ... 4*C2*D2*b*k + E2^2*k^2 - 4*A2*E2*b - 4*A2*C2*b^2 - 4*C2*F2*k^2 - 4*A2*F2; % first attempt, use default calculation precision of Matlab s=solve(eq1,eq2,'k','b'); k = double(s.k); b = double(s.b); % check solution, only real solutions are accepted realFlag = 1; for i=1:size(k,1) if (~isreal(k(i)) ) realFlag = 0; end end % if first attempt fails, increase the calculation precision and try again. if((size(k,1) == 0) || realFlag == 0) disp('Try again, with higher precision!'); digits(20); s=solve(vpa(eq1),vpa(eq2),'k','b'); k = double(s.k); b = double(s.b); end % sometimes the soultions above have a very large k, which is not good for % calculating point of tangency, it should be neglected if(kinf) temp = find(abs(k) > 1e14 ); % hack for nearly infinity line slope k(temp) = []; b(temp) = []; end % calculate point of tangency aa1=A1+B1*k+C1*k.^2; bb1=B1*b+2*C1*k.*b+D1+E1*k; x1 = -bb1./(2*aa1); y1 = k.*x1+b; if(kinf) % deal with infinity k x1 = [sx1;x1]; y1 = [sy1;y1]; end aa2=A2+B2*k+C2*k.^2; bb2=B2*b+2*C2*k.*b+D2+E2*k; x2 = -bb2./(2*aa2); y2 = k.*x2+b; if(kinf) % deal with infinity k x2 = [sx2;x2]; y2 = [sy2;y2]; end % draw point of tangency and common tangent line eum = ['r', 'g', 'b', 'c']; for i=1:size(x1,1) plot(x1(i), y1(i), '*', 'Color',eum(i)); plot(x2(i), y2(i), '*', 'Color',eum(i)); line([x1(i), x2(i)], [y1(i), y2(i)], 'Color',eum(i)); end hold off; end function [A, B, C, D, E, F] = generateEllipse(halfMajor, halfMinor, centerX, centerY, theta) % use a general ellipse equation found in: % https://en.wikipedia.org/wiki/Ellipse a = halfMajor; b = halfMinor; x0 = centerX; y0 = centerY; t = deg2rad(theta); A = a^2*sin(t)^2 + b^2*cos(t)^2; B = 2*(b^2-a^2)*sin(t)*cos(t); C = a^2*cos(t)^2+b^2*sin(t)^2; D = -2*A*x0-B*y0; E = -B*x0-2*C*y0; F = A*x0^2+B*x0*y0+C*y0^2-a^2*b^2; end % Description: test code for FindCommonTangentofTwoEllipse.m % Author: Yang Yang % Date: 2019-02-20 % clean up clc; clear; close all; % test case 1, infinity slope ep1 = [30, 10, 0, 50, 45]; ep2 = [30, 10, 0, -50, 45]; FindCommonTangentofTwoEllipse(ep1, ep2); % test case 2, zero slope ep1 = [50, 20, 70, 0, -10]; ep2 = [50, 20, -70, 0, -10]; FindCommonTangentofTwoEllipse(ep1, ep2); % test case 3, common case ep1 = [30, 20, 50, 10, -10]; ep2 = [45, 20, -50, 20, 30]; FindCommonTangentofTwoEllipse(ep1, ep2); % test case 4, random case a=10+(50-10)*rand(1); b=10+(50-10)*rand(1); xc = 60+(100-60)*rand(1); yc = 60+(100-60)*rand(1); theta = 180*rand(1); ep1 = [a, b, xc, yc, theta]; a=10+(50-10)*rand(1); b=10+(50-10)*rand(1); xc = -(60+(100-60)*rand(1)); yc = 60+(100-60)*rand(1); theta = 180*rand(1); ep2 = [a, b, xc, yc, theta]; save ep1-sep.mat ep1; % save data for debug save ep2-sep.mat ep2; % save data for debug FindCommonTangentofTwoEllipse(ep1, ep2); |
【本文地址】
今日新闻 |
推荐新闻 |