How to find common tangents between 2 ellipses ?

Hi everyone,

I am trying to find the common tangents between 2 ellipses by using solve function. I have f(x1) and g(x2) which are my 2 elipses (ellipse 1 the smallest and ellipse 2) equations and the 2 unknows (x1 and x2) and their derivative at each point which gives me the slope of their tangents at each point.

Then, my equations to solve are:

and here my Matlab code:

syms x1 x2

f_x1 = -sqrt(b1^2*(1-((x1-xC1)/a1)^2)) + yC1;
g_x2 = sqrt(b2^2*(1-((x2-xC2)/a2)^2)) + yC2;

f_x1_prime = b1^2*(x1-xC1)/(a1^2*sqrt(b1^2*(1-((x1-xC1)/a1)^2)));
g_x2_prime = -b2^2*(x2-xC2)/(a2^2*sqrt(b2^2*(1-((x2-xC2)/a2)^2)));

eq1 = (g_x2 - f_x1)/(x2 - x1) == f_x1_prime;
eq2 = f_x1_prime == g_x2_prime;

[S] = solve(eq1, eq2, x1, x2, 'ReturnConditions', true, 'Real', true);

assume(S.conditions)
restriction = [S.x1>xC1, S.x1<xC1+a1, S.x2<xC2, S.x2>xC2-a2];
solx = solve(restriction, S.parameters);

x1 = subs(S.x1, S.parameters, solx);
x2 = subs(S.x2, S.parameters, solx);

As result, I have something in S but the warning and no solution after having assume the condition and solve the parameter.

S = 

  struct with fields:

            x1: [1×1 sym]
            x2: [1×1 sym]
    parameters: [1×1 sym]
    conditions: [1×1 sym]

and the warning:

    Warning: Unable to find explicit solution. For options, see help. 
> In sym/solve (line 317)
In findTangentPoints (line 24)
In testTangent (line 11) 

Does anyone can help me to use this function solve or see any errors ? I am going to be crazy...

Thanks a lot everyone and have a nice evening!

[John D'Errico]

Note that as you have written f and g, you are restricting an ellipse to a specific branch of that curve. That happens as soon as you take a sqrt. In fact, a pair of non-intersecting ellipses will have FOUR common tangent lines.

Answer 1

[Matt J]

alpha1=30; alpha2=-45; %EDIT
a1=0.015; b1=0.020;  xC1=0.13; yC1=0.09;
a2=0.08;  b2=0.10;   xC2=0.27; yC2=0.11;

T1=eye(3); T1(1:2,end)=[-xC1;-yC1]; %translation matrix
T2=eye(3); T2(1:2,end)=[-xC2;-yC2];

E1=diag(1./[a1^2,b1^2 -1]);
E2=diag(1./[a2^2,b2^2,-1]);

R=@(t) [cosd(t), -sind(t) 0; sind(t) cosd(t) 0; 0 0 1]; %EDIT

C1=T1.'*R(alpha1).'*E1*R(alpha1)*T1; %homogeneous ellipse equation matrix
C2=T2.'*R(alpha2).'*E2*R(alpha2)*T2;


syms Lx Ly Lz  %line coefficients

equ1=[Lx,Ly,Lz]/C1*[Lx;Ly;Lz]==0; %tangent to ellipse 1
equ2=[Lx,Ly,Lz]/C2*[Lx;Ly;Lz]==0; %tangent to ellipse 2
equ3=Lx^2+Ly^2+Lz^2==1;       %resolve the scale ambiguity in the coefficients

sol=solve([equ1,equ2,equ3]);
 sol=double([sol.Lx,sol.Ly,sol.Lz]); %convert symbolic solutions to numeric doubles
 sol=sol(2:2:end,:);%discard redundant solutions

%%% Display the solutions
fimplicit(@(x,y) qform(x,y,C1));

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

hold on
fimplicit(@(x,y) qform(x,y,C2));
for i=1:4
 fimplicit(@(x,y) lform(x,y,sol(i,:))); %dumb warnings - ignore
end

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

hold off
axis equal, grid on
axis([0.1,0.4,0,.24])

function val=qform(x,y,Q)

 sz=size(x);
 xy=[x(:),y(:),x(:).^0];
 val=reshape( sum( (xy*Q).*xy ,2) ,sz);

end

function val=lform(x,y,L)

 sz=size(x);

 val=reshape( [x(:),y(:),x(:).^0]*L(:) ,sz);

end

[Youssef Kaddoura]

Thank you for such a valuable and helpful code. What about common tangents for two rotated ellipses by angles alpha1 and alpha2 accordingly.

[Matt J]

I've modified the code to handle rotated ellipses.

[Youssef Kaddoura]

I want to express my sincere appreciation for your quick and positive response.

[Youssef Kaddoura]

Hi Again, would be too much to ask you about the possibility of sharing the code for identifying the points of tangency.

[Matt J]

res=@(Z)Z(:,1:2)./Z(:,3);

T1=res(sol/C1); %tangent points lying on ellipse 1
T2=res(sol/C2); %tangent points lying on ellipse 2

[Youssef Kaddoura]

Thank you Matt for your quick and wonderful help. You are a blessing to all Matlab users!!!

[Varsha Radhakrishnan]

Hi ,

I am looking for a solution to same problem to find the common tangential points between 2 ellipsoids. Could someone help me out?

Thanks,

Varsha

[Matt J]

I am looking for a solution to same problem to find the common tangential points between 2 ellipsoids.

Interesting. What's the application?

[Varsha Radhakrishnan]

By finding the tangential points and use it to create a flight path.

[Matt J]

I would post as a new question since there are many new issues now that it's 3D.

For example, are we talking about tangent planes as opposed to tangent lines? Also, there may be infinite solutions in the 3D case. For example, in the case where the ellipsoids are two spheres of equal diameter, there is an entire tangent cylinder containing the two spheres. How do you plan to resolve those cases?

[Varsha Radhakrishnan]

Thanks Matt

Yes it will be tangent planes and the ellipsoids / spheres are of different diameter at different positions..for eg: one small ellipsoid and another one.....similar to the ellipse as in this previous question.

[Matt J]

and the ellipsoids / spheres are of different diameter at different positions..for eg: one small ellipsoid and another one.....similar to the ellipse as in this previous question.

I don't know if even that guarantees a discrete/finite set of solutions.

[Varsha Radhakrishnan]

@Matt J: Hi Matt, This is a similar figure representing my problem. I need to find the common tangential points passing through these ellipsoids.

I tried solving my equations using matlab..But there is something wrong with my codes. Could you please have a look or guide me by sharing some links to solve this.

a=1;b=2;c=3;d=4;e=5;f=6;g=7;h=8;%ellipsoid coefficients
x=-1.5:.01:1.5;y=x;%x,y axes
for i=1:max(size(x))
  for j=1:max(size(y))
      z(i,j)=real(sqrt((d-a*x(i)^2-b*y(j)^2)/c));%definition of z
  end
end
surf(x,y,z)
u=-1.5:.01:1.5;v=u;%u,v axes for other ellipsoid
for i=1:max(size(u))
    for j=1:max(size(v))
        w(i,j)=real(sqrt((h-e*x(i)^2-f*y(j)^2)/g));%definition of w
    end
end
surf(u,v,w)
zval=max(z(find(c*h*z==g*d*w)))%derived condition solution points maximum
xsol=find(a*h*x==d*e*u);xval=x(xsol)%x-position
ysol=find(b*h*y==f*d*v);yval=y(ysol)%y-position
wval=max(w(find(c*h*z==g*d*w)))%height at which tangent plane touches ellipsoid

Thanks in advance and have a nice day

[Matt J]

I can't understand what the code is doing, but as I said before, the question belongs in its own thread. Also, I am now even more certain that the problem is underdetermined. In the 3D case, you have three unknowns (a plane has 3 degrees of freedom) but only two equations (one for each ellipsoid).

[Salvatore Montella]

Can you please explain that line of the code:

"equ3=Lx^2+Ly^2+Lz^2==1; %resolve the scale ambiguity in the coefficients"

Why do you use it to solve the system?

I need these calculations to model the behaviour of a Tesla-shape superconducting accelerator.

[Matt J]

The coefficients of a line equation are unique only up to a scale factor. In other words x+2y=3 is the same as 2x+4y=6. We need equ3 to remove that ambiguity for the solver.

[Salvatore Montella]

Ok, but what is not clear to me is the reason why the sum of the squares of the coefficients is equal to 1. Can you explain it precisely, please? And is the line in the 3D space, right?

[Matt J]

but what is not clear to me is the reason why the sum of the squares of the coefficients is equal to 1

It's just an arbitrary choice (one of many possibilities) to resolve the scale ambiguity in the coefficients mentioned above.

And is the line in the 3D space, right?

No, the lines are 2D, as you can see from the plots above, but the equation for a 2D line Lx*x+Ly*y+Lz=0 contains three coefficients.

[Salvatore Montella]

Perfect! Thank you so much for your patience and your time. I appreciate your quick response and your code!

[Engineer_Park]

Excuse me sir..how can i find the circle tangent to the above two ellipses when the radius is known?

[Vojtech Zbíral]

Hello, Matt

Iam using the code you provided above with iterations. Basically based on the coordinations of the ellipses a force is generated. This force then moves with the ellipses. The code works fine when I do not use iterations, however when I do, I get this error:

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN. 

pointing to T2 matrix in this part:

res=@(Z)Z(:,1:2)./Z(:,3);
T1=res(sol/C1); %tangent points lying on ellipse 1
T2=res(sol/C2); %tangent points lying on ellipse 2

Please, do you have any idea on how to fix this problem?

kind regards.

[Matt J]

@Vojtech Zbíral Knowing nothing about your iterative process, I can only guess. It appears that your C2 matrix is singular, which means that it does not define an ellipse, but rather some degenerate conic. See also,

https://en.wikipedia.org/wiki/Matrix_representation_of_conic_sections#Classification

[Vojtech Zbíral]

Thank you so much for your quick response. I will look into it. Have a wonderful day.

[Niki Amini-Naieni]

Hi, thank you for this. I am wondering if you could explain and refer me to the resources to understand how you came up with these two equations:

equ1=[Lx,Ly,Lz]/C1*[Lx;Ly;Lz]==0; %tangent to ellipse 1

equ2=[Lx,Ly,Lz]/C2*[Lx;Ly;Lz]==0; %tangent to ellipse 2

[Matt J]

@Niki Amini-Naieni You would have to have studied some projective geometry for any brief explanation to make sense to you. Basically, given an ellipse, E, the equation coefficients of lines tangent to E form a dual ellipse D. Therefore, given two ellipses E1 and E2, the common tangents must therefore lie in the intersection of their duals D1 and D2.

This wiki page talks about deriving the dual ellipse when the given ellipse is non-rotated,

https://en.wikipedia.org/wiki/Dual_curve

The extension to rotated ellipses isn't that hard.

[Niki Amini-Naieni]

Great, thanks for the quick answer!