Axes of Symmetry for a General Ellipse

Let $f(x,y)=Ax^2+Bxy+Cy^2+Dx+Ey+F$. This function represents a quadratic surface, specifically an elliptic paraboloid in our case of $B^2<4AC$. The level curves of this surface form a family of concentric ellipses. The center of this family coincides with an extremum of the surface, so we can find it by solving the equation $\nabla F=0$. This yields $$\mathbf p=\left({BE-2CD\over4AC-B^2},{BD-2AE\over4AC-B^2}\right)$$ for the center of the family of ellipses. Alternatively, one can look for a translation that makes the linear terms vanish, which leads to the same equations, but with a bit more work.

Translating the paraboloid clearly doesn’t change its orientation, so we need examine only the quadratic part $Q(x,y)=Ax^2+Bxy+Cy^2$, which is unchanged by translations, to determine the directions of the principal axes. A fairly simple way to do this is to compute the eigenvectors of the matrix $\bigl(\begin{smallmatrix}A&B/2\\B/2&C\end{smallmatrix}\bigr)$ of this quadratic form. It’s a symmetric real matrix, so its eigenspaces are orthogonal and correspond to the ellipses’ axes: the minor axis to the eigenspace of the greater of the two eigenvalues, the major axis to that of the lesser. A bit of work produces the eigenvalues $$\frac12\left(A+C\pm\sqrt{(A-C)^2+B^2}\right)$$ with corresponding eigenvectors $$\left(A-C\pm\sqrt{(A-C)^2+B^2},B\right).$$You can verify that these vectors are indeed orthogonal.

Another way to compute these direction vectors is to find the extrema of $Q$ restricted to the unit circle. As with the eigenvalues, the major axis will lie in the direction of the local minima, the minor axis in the direction of the maxima. This can be solved via Lagrange multipliers (you will find that they are equal to the eigenvalues of $Q$) or by pulling $Q$ back by $\phi:t\mapsto(\cos t,\sin t)$, though this last method might lead to “messy” trigonometric manipulations.

With the center of the ellipse and direction vectors of the axes in hand, it’s a simple matter to construct the equations of the axis lines.

The principal axes are:

$$\left( y-\frac{2AE-BD}{B^2-4AC} \right)= \frac{C-A \color{red}{\pm} \sqrt{(A-C)^2+B^2}}{B} \left( x-\frac{2CD-BE}{B^2-4AC} \right)$$

Also refer to another answer.

Idea for solution using multivariable calculus:

(1) Find the max an min of the functions $H(x,y) = x$ , $V(x,y) = y$ with the restriction $Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$ (easy with Lagrange multipliers). This gives a "bounding box" for the ellipse. The center of the box will be the center of the ellipse.

(2) Find the max and min distances of $(x,y)$ to the center of the ellipse with the restriction $Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0$.