From: Dave Rusin <rusin@math.niu.edu>
Subject: Re:  distance to ellipsoid
Date: Mon, 30 Aug 1999 17:44:15 -0500 (CDT)
Newsgroups: [missing]
To: aniem@gub.ruhr-uni-bochum.de
Keywords: Lagrange multipliers: nearest point on an ellipsoid to a given point

>Dear Dave Rusin,
>I have read your text
>96/normals
>titled '' How many Surface Normals ?'' and I hope that my 
>question will be just a special case.
>
>I am looking for a distance between an ellipsoid 
>F(x,y,z) = (x/a)^2 + (y/b)^2 + (z/c)^2 - 1 = 0 and a
>point  P = (x0,y0,z0) inside the ellipsoid but 
>not (0,0,0). 
>Elimination of x1, y1, z1 from these equations leads to a nasty
>polynomial equation  which  Mathematica does not solve. 
>Where can I find an analytical solution?

You asked for the distance between P = (x0,y0,z0)  and the ellipsoid F=0
where  F(x,y,z) = (x/a)^2 + (y/b)^2 + (z/c)^2 - 1. There is really no
trick here: you are seeking to minimize  D = (x-x0)^2 + (y-y0)^2 + (z-z0)^2
on the set where  F = 0, and so just need to solve the simultaneous
equations  { x/a^2 = lambda * (x-x0), y/b^2 = lambda * (y-y0), 
z/c^2 = lambda * (z-z0), (x/a)^2 + (y/b)^2 + (z/c)^2 - 1 = 0}. It is
trivial to eliminate  lambda; then, two of the equations can be used
to solve for, say, x and y in terms of  z  and the constants. This reduces
the situation to one equation in one unknown  (z)  which I find to be
of degree 6. I attach this polynomial below for your convenience.
Mathematica cannot solve this equation symbolically (the Galois group is
generically  Sym(6),  which means in particular that the equation cannot
be solved using only the extraction of n-th roots), but there's no reason
it cannot be solved numerically.

If you think about points withing an ellipsoid shaped like a long cigar
or thin plate, it is clear that there should be at least two points where
the distance from  P  to the surface will attain a local minimum, and at
least two where it will attain a local maximum. I suppose I ought to be
able to visualize a third pair of points separated by a line segment
more or less perpendicular to the first two, where the distance function
again attains a critical point; for some reason I'm not seeing them in
my head, but their existence would justify the fact that the problem
appears to have 6 solutions in general.

dave

The polynomial to be solved is

z^6 *
 (b-c)^2*(b+c)^2*(a-c)^2*(a+c)^2
+z^5 *
 2*c^2*z0*(b-c)*(b+c)*(a-c)*(a+c)*(b^2+a^2-2*c^2)
+z^4 *
 -c^2*(-2*a^4*b^2*c^2-6*c^6*z0^2+6*c^4*z0^2*b^2+a^4*b^4-a^4*c^2*z0^2+c^4*b^4-2*
 c^6*b^2-2*a^2*c^2*b^4-4*c^2*a^2*z0^2*b^2-2*a^2*c^6+6*c^4*a^2*z0^2+c^8+a^4*c^4+
 4*a^2*c^4*b^2+2*c^2*a^2*x0^2*b^2-c^4*y0^2*b^2-a^4*y0^2*b^2-c^2*z0^2*b^4-b^4*x0
 ^2*a^2-c^4*a^2*x0^2+2*c^2*b^2*y0^2*a^2)
+z^3 *
 -2*c^4*z0*(-c^2*b^4+c^2*y0^2*b^2+2*c^4*z0^2-4*a^2*b^2*c^2+3*b^2*c^4-2*c^6-c^2*
 z0^2*b^2-b^2*x0^2*a^2+3*c^4*a^2-a^2*y0^2*b^2+b^2*a^4-a^4*c^2-c^2*a^2*z0^2+c^2*
 a^2*x0^2+b^4*a^2)
+z^2 *
 -c^6*z0^2*(4*a^2*b^2-6*a^2*c^2-c^2*z0^2-6*b^2*c^2-y0^2*b^2-x0^2*a^2+6*c^4+a^4+
 b^4)
+z^1 *
 -2*c^8*z0^3*(b^2+a^2-2*c^2)
+z^0 *
 -c^10*z0^4