Furthest point in direction ellipsoid with Newton's method I'm trying to find the furthest point on the surface of an ellipsoid in a given direction. To do this I figured the simplest method would be to use Newton's method.
$\vec{x}(\theta,\phi)=r*\begin{bmatrix}
a*\cos(\phi)*\cos(\theta)\\ 
b*\cos(\phi)*\sin(\theta)\\ 
c*\sin(\phi)
\end{bmatrix}$
let $\vec{d}=\begin{bmatrix}
x\\ 
y\\ 
z
\end{bmatrix}$ be the direction we are trying to find the furthest point in.
Then it must me true that the furthest point is $\max\limits_{\phi,\ \theta} \vec{x}(\theta,\phi) \cdot \vec{d}$
From there it should just be a matter of maximizing the dot product of the dot product mentioned above.
$f(\theta, \phi)=\vec{x}(\theta,\phi) \cdot \vec{d}= r*(a*x*\cos(\phi)*\cos(\theta) + b*y*\cos(\phi)*\sin(\theta) + c*z*\sin(\phi))$
$F(\theta,\phi)=\frac{\partial f}{\partial \theta}=r*(b*y*\cos(\phi)*\cos(\theta)-a*x*\cos(\phi)*\sin(\theta))=0$
$G(\theta,\phi)=\frac{\partial f}{\partial \phi}=-r*(a*x*\cos(\theta)*\sin(\phi)+b*y*\sin(\theta)*\sin(\phi)-c*z*\cos(\phi)=0$
Find the jacobian:
$J(\theta,\phi)=\begin{bmatrix}
\frac{\partial F}{\partial \theta} & \frac{\partial F}{\partial \phi}\\ 
\frac{\partial G}{\partial \theta} & \frac{\partial G}{\partial \phi}
\end{bmatrix}$
Where:
$\frac{\partial F}{\partial \theta}=-r*(b*y*\cos(\phi)*\sin(\theta)+a*x*\cos(\phi)*\cos(\theta))$
$\frac{\partial F}{\partial \phi}=-r*(b*y*\sin(\phi)*\cos(\theta)-a*x*\sin(\phi)*\sin(\theta))$
$\frac{\partial G}{\partial \theta}=r*(a*x*\sin(\theta)*\sin(\phi)-b*y*\cos(\theta)*\sin(\phi))$
$\frac{\partial G}{\partial \phi}=-r*(a*x*\cos(\theta)*\cos(\phi)+b*y*\sin(\theta)*\cos(\phi)-+c*z*\sin(\phi))$
Finally, using newton's method:
$\begin{bmatrix}
\theta_{n+1}\\ 
\phi_{n+1}
\end{bmatrix}
=\begin{bmatrix}
\theta_{n}\\ 
\phi_{n}
\end{bmatrix}-J^{-1}(\theta_{n},\phi_{n})\begin{bmatrix}
F(\theta_{n},\phi_{n})\\ 
G(\theta_{n},\phi_{n})
\end{bmatrix}$
Some example inputs that give me the wrong result:
$\vec{d}=\begin{bmatrix}
-0.170260981\\ 
0.141882509\\ 
-0.975131035
\end{bmatrix}$
a=5, b=2.5, c=2.5 
To get the initial guess at the angles (which I can confirm gives a good initial guess):
theta = atan2(ay, bx)
phi = atan2(z, c*sqrt((x/a)^2 + (y/b)^2))
This leaves me with an output of (after 10 iterations): $\begin{bmatrix}
1.63307858\\ 
-0.340220749\\ 
2.33827138
\end{bmatrix}$
Which we can see isn't even a point in the same direction as the initial direction!?
$\begin{bmatrix}
1.63307858\\ 
-0.340220749\\ 
2.33827138
\end{bmatrix} \cdot \begin{bmatrix}
-0.170260981\\ 
0.141882509\\ 
-0.975131035
\end{bmatrix}=-2.6064419746399$
That being said it seemingly works for some directions, but not others. I must be approaching this wrong or have messed up the algebra. Any help or suggestions?
 A: I'm not willing to slog through your algebra, but another method you can use is to say that at the farthest point, the normal vector to the ellipsoid (which is pretty easy to write down) must be parallel to $d$. This gives you three equations to solve, and you're done (except that you'll also find the NEAREST point). 
Note that the "farthest point" may not be well defined -- if your ellipsoid is, for instance, a sphere, then from the origin, for any direction, there are two equally distance points on the surface.
As yet another alternative: consider your ellipsoid as the result of applying a matrix transformation to all points of the unit sphere; in your case, the matrix is 
$$
M = \begin{bmatrix}
ra & 0 & 0 \\0 & rb & 0 \\ 0 & 0 & rc
\end{bmatrix}.
$$
You seek a point $p$ of the unit sphere (i.e., $\|p\| = 1$) with 
$$
f(p) = \langle d, Mp \rangle
$$
as large as possible. Well, 
\begin{align}
f(p) &= d^t (Mp) & \text{by definition of inner product}\\
&= (d^t M) p & \text{associativity of matrix multiply} \\
&= (M^t d)^t p & \text{property of transpose}
\end{align}
If you let $d' = M^t d$, then you're seeking the point of the unit sphere that maximizes $d' \cdot p$, and that point is simply $p_0 = \frac{1}{\|d'\|}d'$. You then multiply this point $p_0$ by the matrix $M$ to get $Mp_0$, which is your desired maximum-distance point. (Although $-Mp_0$ also works). 
