Context: My current work involves ellipsoidal cones - the intersection of a cone and ellipsoid where the origin of each is at $(0,0,0)$ and where I say it is in 'canonical position' when the cone axis is the $(0,0,1)$ vector and the ellipsoid is rotated accordingly.
For visualization, the pictures
show two orientations of the same example where the ellipsoid radii are $(36.8, 23.2, 15.99)$, the elliptical cone's half angles are $(23.67, 52.8)$ degrees, and the cone's position within the ellipsoid is rotated as desired (I can provide exact details if needed). The triangle empirically approximates the $z$-maximizer point.
The red contour represents the $P(t)$ function for $360$ values of $t \in [0, 2\pi]$ and $P$ is the $z \ge 0$ intersection of the line through $f(t)$ on the ellipse and the ellipsoid itself in canonical ellipsoidal cone position. $f(t)$ is the elliptic curve formed by the intersection of the plane $z=k$ with the cone, where for simplicity, $k=1$. The points from $f(t)$ form the lines connected from the cone origin that I intersect with the ellipsoid. The intersection of the line and the ellipsoid is done via the standard quadratic equation method described by Paul Bourke (http://www.paulbourke.net/geometry/circlesphere/) where of course $f(t)$ is first rotated into canonical (axis-aligned) ellipsoid position.
Problem To Solve: I want to analytically compute $min/max$ $z(t)$ where $t$ is parameter $\in [0, 2*\pi]$ on the ellipse formed by intersecting the cone with the plane $z=k$. For simplicity I use $k=1$.
The standard methodology for solving my problem analytically is to create an expression for $z(t)$, differentiate w.r.t. $t$, set to $0$, and solve for $t$. Of course, I can solve this numerically through a 'line search' on $t$ but would like to do better.
My Approach: Briefly, the Bourke's line-ellipsoid intersection method is this:
Find points $P (x,y,z)$ on a line defined by two points $P_1 (x_1,y_1,z_1)$ and $P_2 (x_2,y_2,z_2)$ which is described by $P = P_1 + \mu(P_2 - P_1)$ or in each coordinate $x(\mu) = x_1 + \mu(x_2 - x_1)$, $y(\mu) = y_1 + \mu(y_2 - y_1)$,and $z(\mu) = z_1 + \mu(z_2 - z_1)$.
An axis-aligned ellipsoid centered at $P_3 (x_3,y_3,z_3)$ with radii $a$, $b$, and $c$ is described by
$(\frac{x - x_3}{a})^2 + (\frac{y - y_3}{b})^2 + (\frac{z - z_3}{c})^2 = 1^2$.
Substituting the equation of the line into the sphere gives a quadratic equation of the form
$A\mu^2 + B\mu + C = 0$
Because $P_1$ and $P_3$ are both $(0,0,0)$, this simplifies dramatically to
$A\mu^2 + C =0$ where
$A = [(bcx_2)^2 + (acy_2)^2 + (abz_2)^2]$ and
$C = -(abc)^2$. I can provide the algebra if requested.
Now, note that $A$ is really $A(t)$ because $P_2$, the point on the cone's ellipse at $z=1$ is really $P_2(t)$.
$P_2(t) = [r_{x(1)}cos(t), r_{y(1)}sin(t), 1]$ where $z=1$ and $r_{x(1)}$ and $r_{y(1)}$ are the radii of the cone ellipse on the plane $z=1$.
Note that the positive/negative roots of the quadratic, $\mu$, are actually the $z(t)$-values I will need in my solution...so solving for $z(t)=\mu(t)$ would answer my question.
Finally, note that to put $P_2$ in axis-aligned ellipsoid space, apply a rotation matrix, $R*P_2$, and substitute this into $A(t)$.
Now we finally have a means of expressing the quadratic as a function of $t$.
Again, I can add the ugly algebra, but I end up with an equation $A(t)\mu^2 + C = 0$ or
$A(t)z(t)^2 + C = 0$
that will contain a collection of $cos^2(t)$, $sin^2(t)$, and $cos(t)sin(t)$ terms that I will then need to differentiate.
$A(t) = [(bc(R_{11}cos(t) + R_{12}sin(t) + R_{13}))^2 + (ac(R_{21}cos(t) + R_{22}sin(t) + R_{23}))^2 + (ab(R_{31}cos(t) + R_{32}sin(t) + R_{33}))^2]$
where the $R_{ij}$ values are the elements of the rotation matrix $R$ and
$z(t) = [R_{31}r_{x}cos(t) + R_{32}r_{y}sin(t) + R_{33}]$
Question for the Community Do I have the right approach? How do I get rid of the $cos(t)sin(t)$ terms in the derivative when solving for $t$? I can handle a derivative with $cos^2$ and $sin^2$ but not one with $cos(t)sin(t)$ terms.
Thank you!
I would approach this problem by using the Lagrange multiplier method. Let the vector $k = [0,0,1]^T$, and let $ r = [x,y,z]^T $, then you cone and ellipsoid equations are
$ r^T Q_c r = 0 $ (for the cone), and
$r^T Q_e r = 1 $ (for the ellipsoid)
And you want to maximize $z = k^T r $, so define
$ G(r, \lambda_1, \lambda_2) = k^T r + \lambda_1( r^T Q_c r) + \lambda_2 ( r^T Q_e r - 1) $
Take the gradient of $G$ and set it to zero,
$ \nabla G = k + \lambda_1 (2 Q_c r) + \lambda_2 (2 Q_e r ) = \mathbf{0} $
Since vectors $i = [1,0,0]^T$ and $j = [0,1,0]^T$ satisfy $i^T k = 0 $ and $j^T k = 0 $ , then we now have two equations
$ i^T \nabla G = \lambda_1 (2 i^T Q_c r) + \lambda_2 (2 i^T Q_e r) = 0 $
$ j^T \nabla G = \lambda_1 (2 j^T Q_c r) + \lambda_2 (2 j^T Q_e r ) = 0 $
Since $\lambda_1$ and $\lambda_2$ cannot be zero, this implies that the determinant of this homogenous system of equations must be zero, i.e.
$ (2 i^T Q_c r) (2 j^T Q_e r) - (2 j^T Q_c r) (2 i^T Q_e r ) = 0 $
Simplifying and re-arraning, this becomes
$ r^T ( Q_c (i j^T - j i^T ) Q_e ) r = 0 $
Note that
$ i j^T - j i^T = \begin{bmatrix} 0 && 1 && 0 \\ 0 && 0 && 0 \\ 0 && 0 && 0 \end{bmatrix} - \begin{bmatrix} 0 && 0 && 0 \\ 1 && 0 && 0 \\ 0 && 0 && 0 \end{bmatrix} = \begin{bmatrix} 0 && 1 && 0 \\ -1 && 0 && 0 \\ 0 && 0 && 0 \end{bmatrix} $
If we call this matrix $E$, then our three equations that we have to solve are
$r^T Q_c r = 0 \tag{1} $
$r^T Q_e r = 1 \tag{2} $
$r^T Q_c E Q_e r = 0 \tag{3} $
And these can be solved using Mathematica or Sage
I've worked an example in this Sage Worksheet. It contains the plot of an ellipsoid and a cone, their intersection curve for $z \gt 0$, and the critical points which are solutions of $(1),(2),(3)$.