Find the point $P$ on an ellipse such that $\overline{AP} + \overline{BP}$ is minimum for given points $A$ and $B$

Question

An ellipse in $3D$ space is specified in parametric form as follows:

$ E(t) = V_0 + V_1 \cos t + V_2 \sin t $

In addition, two points in space $A$ and $B$ are given. What I would like to find is the point $P$ on the ellipse $E(t)$ that will minimize the distance sum $\overline{AP} + \overline{BP} $

Initial remarks:

First, note that if points $A$ an $B$ lie in the same plane as the given ellipse, then point $P$ is the reflection point where a ray from point $A$ landing on the given ellipse at $P$ gets reflected into a ray passing through point $B$.

It is possible to just evaluate the given sum of distances as a function of the parameter $t$, and then explicitly differentiate this function, and set its derivative to zero, and solve the resulting equation numerically using Newton's method for the value of $t$.

Apart from that, I'd like to know if there are alternative ways to view this problem and its solution method.

Answer 1

Given a point $P = E(t)$ for some $t$, it lies on a unique ellipsoid (of revolution) whose foci are $A$ and $B$. It was shown here that the equation of this ellipsoid is

$ (p - C)^T (a^2 I - {UU}^T ) (p - C) = a^2 (a^2 - c^2) \tag{1}$

where $C = \frac{1}{2}(A + B)$, $U = \frac{1}{2}(B - A) $ and $c = \| U \| $.

This gives the first relation that the minimizing point $P$ must satisfy. There are basically two unknowns here, which are $t$ and $a^2$. At the minimizing $t$, we have $\dfrac{d}{dt}(a^2) = 0$. Using this fact and differentiating $(1)$ with respect to $t$, gives us

$ (p - C)^T (a^2 I - {UU}^T) \dot{p} = 0 \tag{2}$

where $\dot{p} = \dfrac{d p }{dt } $.

We can solve $(1), (2)$ using the Newton-Raphson multivariate method.

Depending on the initial guess, we will get either the minimizing $t$ or the maximizing $t$ and the corresponding value of $a^2$.

Here's the MATLAB code that implements the solution of these two nonlinear equations. In this example, I took $V_0 = [2, 3, 7]^T, V_1=[1, 0.5, 0.5]^T , V_2 = [-2, 6, 2]^T, A = [5,10,2]^T, B = [1,-5,6]^T $.

function F = find_f(y)
t = y(1);
u = y(2);
v0 = [2;3;7];
v1 = [1; .5; .5];
v2 = [-2; 6; 2];
A = [5; 10; 2];
B = [1;-5;6];
C = 0.5(A + B);
U = 0.5(B - A);
c = norm(U);
Q = ueye(3) - UU.';
p = v0 + v1cos(t)+ v2sin(t);
pp = - v1sin(t) + v2cos(t);
% (p - C)^T (u I - UU^T) (p - C) = u ( u - c^2)
F(1) = (p - C).'* Q * (p - C ) - u * (u - c^2);
% (p - C)^T (u I - UU^T ) (p') = 0
F(2) = (p - C).'* Q * pp;
end
clc
v0 = [2;3;7];
v1 = [1; .5; .5];
v2 = [-2; 6; 2];
A = [5; 10; 2];
B = [1;-5;6];
C = 0.5(A + B);
U = 0.5(B - A);
c = norm(U);
mindist = 1000;
maxdist = 0;
tmin = 0;
tmax = 0;
for t = 0:0.1:6.28
    p = v0 + v1cos(t)+ v2sin(t);
    dist = norm(p - A) + norm(p - B);
    if dist < mindist
        mindist = dist;
        tmin = t;
    end
    if dist > maxdist
        maxdist = dist;
        tmax = t;
    end
end
u1 = (0.5mindist)^2
u2 = (0.5maxdist)^2
tmin
tmax
y0 = [tmin; u1];
fun = @find_f;
[x, fval] = fsolve( fun , y0)
t = x(1);
u = x(2);
u - c^2
dist = 2*sqrt(u)
y0 = [tmax; u2];
[x, fval ] = fsolve(fun, y0)
t = x(1);
u = x(2);
u - c^2
dist = 2*sqrt(u)

---

The output for the example included in the above code is shown below. From the generated output, we find that the minimum round trip distance is $16.0835$ and the maximum round trip distance is $23.0181$.

---

u1 =
64.6787
u2 =
132.3988
tmin =
3.8000



tmax =
1.5000



Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
<stopping criteria details>
x =
3.7792

64.6696
fval =
1.0e-11 *
-0.0735    0.1222
ans =
0.4196



dist =
16.0835
Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
<stopping criteria details>
x =
1.5354

132.4584
fval =
1.0e-09 *
-0.9659   -0.0780
ans =
68.2084
dist =
23.0181

---

Answer 2

Given $p = (x,y,z)'$ eliminating $t$ in the equation

$$ p = v_0+v_1\cos t+v_2\sin t $$

with the given values

$$ \cases{ v_0 = (2,3,7)'\\ v_1 = (1,0.5,0.5)'\\ v_2 = (-2,6,2)\\ p_A = (5,10,2)'\\ p_B=(1,-5,6)'} $$

we obtain equivalently

$$ \cases{ 3 y + 36 + 2 x - 7 z = 0\\ 28 x z + 20 z^2 - 336z + 1404 - 264 x + 17 x^2=0 } $$

Note that $3 y + 36 + 2 x - 7 z = 0$ is the plane containing the ellipse. Now eliminating $z$ we arrive at

$$ 92 x y - 304y + 20 y^2 + 1116 - 856 x + 145 x^2 = 0 $$

this is the ellipse equation referred to it's plane.

Now eliminating $z$ in

$$ \cases{ d = \|p-p_A\|+\|p-p_B\|\\ 3 y + 36 + 2 x - 7 z = 0} $$

we obtain the distance $d$ as a function of points pertaining to the ellipse plane or

$$ 573049 - 15838 d^2 + 49 d^4 - 60560 x + 1048 d^2 x + 1600 x^2 - 212 d^2 x^2 - 281604 y + 788 d^2 y + 14880 x y - 48 d^2 x y + 34596 y^2 - 232 d^2 y^2 = 0 $$

now at minimum/maximum the two curves

$$ \cases{ 92 x y - 304y + 20 y^2 + 1116 - 856 x + 145 x^2 = 0\\ 573049 - 15838 d^2 + 49 d^4 - 60560 x + 1048 d^2 x + 1600 x^2 - 212 d^2 x^2 - 281604 y + 788 d^2 y + 14880 x y - 48 d^2 x y + 34596 y^2 - 232 d^2 y^2 = 0 } $$

are tangent. Now eliminating $y$ between then, we arrive at a polynomial which contains a double root at the tangency point so,

$$ P_d(x)=k_0(x-x_t)^2(x^2 + c_1 x + c_0)\ \ \ \forall \{x,d\} $$

or equivalently

$\alpha_i(d) = 0,\ \ \{i = 0,\dots,4\}$ in $P_d(x)-k_0(x-x_t)^2(x^2 + c_1 x + c_0)=\sum_{i=0}^4 \alpha_i(d)x^i$

Finally, using Groebner bases we can eliminate $\{k_0,c_1,c_0,x_t\}$ from the polynomial system $\alpha_i(d) = 0$ thus obtaining a polynomial in $d$ which contains our solutions.

Follows a MATHEMATICA script showing this process.

v0 = {2, 3, 7};
v1 = {1, 0.5, 0.5};
v2 = {-2, 6, 2};
pA = {5, 10, 2};
pB = {1, -5, 6};
p = {x, y, z};
equs = p - (v0 + v1 Cos[t] + v2 Sin[t]);
equsxy = Quiet@Eliminate[equs == 0, t];
planeE = equsxy[[1]];
ellxy = Quiet@Eliminate[equs == 0, {t, z}];
distxy = Quiet@Eliminate[{Sqrt[(p - pA) . (p - pA)] + Sqrt[(p - pB) . (p - pB)] - d == 0, planeE}, z];
equ0 = First[Eliminate[{ellxy, distxy}, y]] - k0 (x - xt)^2(x^2 + c1 x + c0);
coefk0 = Solve[D[equ0, x, x, x, x]/4! == 0, k0][[1]];
equ00 = equ0 /. coefk0;
rels = CoefficientList[equ00, x];
vars = Variables[rels]
Length[rels]
pold = Eliminate[rels == 0, {c0, c1, xt}];
dist = Solve[pold, d, Reals] // N

thus obtaining

$$ \cases{ d_{min}= 16.083482892850036\\ d_{max}= 23.018111746188083 } $$