solx =

Are the ellipses in standard form (i.e. their axes parallel to x- and y-axis) ?

28 views (last 30 days)

Show older comments

I have two ellipse using the ellipse equation. I need to know whether the two ellipses are intersecting or not as a Boolean value. True if intersecting, and False if not intersecting. The usage of Solve function is not working in my matlab version. So please suggest me a numerical way to write this code. I have copied the code I am using for it.

ellipseTwo = '((x+0)^2)/2^2 + ((y-5)^2)/4^2 = 1'; ellipseOne = '((x+0)^2)/2^2 + (y-5)^2/4^2 = 1';

ezplot(ellipseOne, [-10, 10, -10, 10]);

hold on

ezplot(ellipseTwo, [-10, 10, -10, 10]);

Sathiya
on 4 Jul 2024

Bruno Luong
on 4 Jul 2024

Edited: Bruno Luong
on 4 Jul 2024

I don't think you get my idea. There is no approximation if tha aspect ratio of two ellipse are identical. By changing parameter, it becomes circle equation with new coordinates. This is NOT an approximation, but exact transformation.

If the aspect ratio is not identical (or radom tilted ellipses the problem become minimizing uadratic function under disc contraints; and ca, be solved with traust region algorithm. I won't go there yet since you have not confimed the asumption of aspect ratio.

Torsten
on 4 Jul 2024

Edited: Torsten
on 4 Jul 2024

I'm surprised that we get a division by zero by the below solution method if the aspect ratio of the two ellipses is equal.

But I think sigma18 and sigma19 have a common divisor. So the problem might vanish if you use

sol = simplify(solve(e1subs,y,'MaxDegree',4))

instead of

sol = solve(e1subs,y,'MaxDegree',4)

below.

If all solutions are complex, the two ellipses don't intersect. If not, you have intersections.

If you invest a little effort, use pencil and paper and deduce the coefficients of the polynomial e1subs of degree 4. Then you can use the numerical "roots" command to deduce the possible intersections (and thus the information whether the ellipses intersect or not).

syms x y

syms x1 x2 y1 y2 real

syms a1 b1 a2 b2 real positive

e1 = (x-x1)^2/a1^2 + (y-y1)^2/b1^2 == 1;

e2 = (x-x2)^2/a2^2 + (y-y2)^2/b2^2 == 1;

solx = solve(e1*a1^2-e2*a2^2,x)

e1subs = subs(e1,x,solx)

sol = solve(e1subs,y,'MaxDegree',4)

solsubs = subs(sol,[a1 a2 b1 b2 x1 x2 y1 y2],[4 4 2 2 0 -5 5 2])

solsubs(1:4)

Torsten
on 4 Jul 2024

Edited: Torsten
on 4 Jul 2024

There are still special cases for which the below code doesn't work, e.g. if Y1 = Y2.

% E1: (x-X1)^2/A1^2 + (y-Y1)^2/B1^2 = 1

% E2: (x-X2)^2/A2^2 + (y-Y2)^2/B2^2 = 1

%

% Numerical values for the ellipse parameters

A1 = 4;

A2 = 4;

B1 = 2;

B2 = 2;

X1 = 0;

X2 = -5;

Y1 = 5;

Y2 = 2;

%

syms x y

syms x1 x2 y1 y2 real

syms a1 b1 a2 b2 real positive

e1 = (x-x1)^2/a1^2 + (y-y1)^2/b1^2 == 1;

e2 = (x-x2)^2/a2^2 + (y-y2)^2/b2^2 == 1;

soly = solve(e1*b1^2-e2*b2^2,y);

e1subs = subs(e1,y,soly);

poly_x = fliplr(coeffs(lhs(e1subs)-rhs(e1subs),x));

sol_x = roots(subs(poly_x,[a1 a2 b1 b2 x1 x2 y1 y2],[A1 A2 B1 B2 X1 X2 Y1 Y2]));

poly_y = fliplr(coeffs(lhs(e1)-rhs(e1),y));

count = 0;

for i = 1:numel(sol_x)

if abs(imag(sol_x(i))) < 1e-8

sol_y = roots(subs(poly_y,[a1 a2 b1 b2 x1 x2 y1 y2 x],[A1 A2 B1 B2 X1 X2 Y1 Y2 sol_x(i)]));

for j = 1:2

if abs(imag(sol_y(j))) < 1e-8 & abs(subs(lhs(e2),[a1 a2 b1 b2 x1 x2 y1 y2 x y],[A1 A2 B1 B2 X1 X2 Y1 Y2 sol_x(i) sol_y(j)])-1) < 1e-8

count = count + 1;

SOL_X(count) = sol_x(i);

SOL_Y(count) = sol_y(j);

end

end

end

end

SOL_X

SOL_Y

if isempty(SOL_X)

disp('Ellipses don''t intersect')

else

disp('Ellipses intersect')

end

Aquatris
on 4 Jul 2024

Edited: Aquatris
on 4 Jul 2024

Here is one way where accuracy will depend on selecting the threshold correctly. You simply solve the equation for y given an x vector, and then find the distance between the elips 1 points and elips 2 points. Then find if any distance is less than a threshold

x = -10:0.01:10;

% solve the equation for y and get index of imaginary solutions

y1 = sqrt(-((x+5).^2/4^2-1)*2^2);idx1 = imag(y1) ==0;

y2 = sqrt(-((x+0).^2/4^2-1)*2^2);idx2 = imag(y2) ==0;

% form x-y points of the elipses

x1 = x(idx1) ;x1 = [x1 flip(x1)];

y1 = y1(idx1);y1 = [y1 -flip(y1)]+2;

x2 = x(idx2) ;x2 = [x2 flip(x2)];

y2 = y2(idx2);y2 = [y2 -flip(y2)]+5;

% loopt over points and calculate their distance, this might be vectorized

for i = 1:length(x1)

for j = 1:length(x2)

dist(i,j) = sqrt((x1(i)-x2(j))^2+(y1(i)-y2(j))^2);

end

end

% find indeces where distance is below a threshold

[a,b] = find(dist<2e-3);

% plot

figure(2)

plot(x1,y1,x2,y2,x1(a),y1(a),'mx',x2(b),y2(b),'ko')

axis([-10 10 -10 10])

Bruno Luong
on 4 Jul 2024

Edited: Bruno Luong
on 5 Jul 2024

% Ti and ci define two ellipses E1 and E2 of the form

% (x,y) st norm(Ti*Xi) <= 1 where Xi := [x;y]-ci

c1 = [0; 5]; % circle center #1

c2 = [-5; 2]; % circle center #2

T1 = [1/4 0;

0 1/2]; % Change variable matrix

T2 = T1;

T = T2/T1;

A = 2*(T'*T);

dc = c1-c2;

dc2 = T2*dc;

b = 2*T*dc2;

c = dc2'*dc2;

% TrustRegion.m File from FEX

% https://www.mathworks.com/matlabcentral/fileexchange/27596-least-square-with-2-norm-constraint

Xx = TrustRegion(A, b, 1);

% X = T1\Xx + c1

f = 1/2*Xx'*A*Xx+Xx'*b+c;

isintersect = f <= 1

Bruno Luong
on 5 Jul 2024

Edited: Bruno Luong
on 5 Jul 2024

Use FSOLVE to find numerical solution. Only xaimum one can be found

c1 = [0; 5]; % circle center #1

c2 = [-5; 2]; % circle center #2

T1 = [1/4 0;

0 1/2]; % Change variable matrix

T2 = T1;

[X, fval, exitflag] = fsolve(@(X) myfun(X,c1,T1,c2,T2), randn(2,1))

isintersect = norm(fval, Inf) < 1e-6 && exitflag >= 0

figure

hold on

h1=fimplicit(@(x,y) sum((T1*([x(:),y(:)]'-c1)).^2,1)-1,[-10 10]);

h2=fimplicit(@(x,y) sum((T2*([x(:),y(:)]'-c2)).^2,1)-1,[-10 10]);

h1.Color = 'k';

h2.Color = 'k';

plot(X(1),X(2),'ro','MarkerSize',10)

function F = myfun(X,c1,T1,c2,T2)

F = [sum((T1*(X-c1)).^2,1)-1;

sum((T2*(X-c2)).^2,1)-1];

end

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!