With the center at $(0,0)$ and $\psi=0$, the equation of the ellipse is $$\frac{x^2}{a^2}+\frac{y^2}{b^2}=1\text{.}$$
So $$\frac{2x}{a^2}+\frac{2y}{b^2}\frac{dy}{dx}=0$$ and therefore $$\frac{dy}{dx}= -{\frac{xb^2}{ya^2}}\text{.}$$
The slope of the normal line is the negative reciprocal of this, so $$\tan(\phi)=\frac{ya^2}{xb^2}\text{.}$$
Meanwhile, $$\tan(\theta)=\frac{y}{x}\text{.}$$ So, eliminating $\frac{y}{x}$ from these two equations and clearing denominators, the relationship between $\phi$ and $\theta$ is: $$b^2\tan(\phi)=a^2\tan(\theta)$$
Matrix derivation
The original equation can be described in a more concise form as
$$\vec x A \vec x^T + \vec x \vec b^T + \vec b \vec x^T + f =0$$
or
$$\vec x A \vec x^T + 2\vec x \vec b^T + f =0$$
where $\vec b=(d,e)$ and $A=\begin{pmatrix}a&b\\b&c\end{pmatrix}$.
If the vector $\vec x_0$ represents the center then in a coordinate system with the origin in this point the equation will be
$$(\vec x-\vec x_0) A (\vec x-\vec x_0)^T + g =0$$
which can be rewritten as
$$\vec x A \vec x^T - 2\vec x A\vec x_0^T + \vec x_0 A \vec x_0^T +g.$$
Which means the we need to choose $\vec x_0$ in such way that $-2A\vec x_0^T=2\vec b^T$, i.e. $$A\vec x_0^T=-\vec b^T,$$
$$\begin{pmatrix}a&b\\b&c\end{pmatrix}
\begin{pmatrix}x\\y\end{pmatrix}=
\begin{pmatrix}d\\e\end{pmatrix}.$$
This is precisely the following system of linear equations written in the matrix form:
$$
\begin{align*}
ax + by &= -d\\
bx + cy &= -e
\end{align*}
$$
(To give credit where credit is due, I was shown this derivation by a colleague of mine who taught labs for the same subject as I did, but with another group of students.)
Derivatives
We can view the curve as a level set of the function
$$F(x,y)=ax^2+2bxy+cy^2+2dx+2ey.$$
The level sets of this function create the family of curves $ax^2+2bxy+cy^2+2dx+2ey+f=0$, where the parameter $f$ is changing.
For example, if $\delta>0$ we get concentric ellipses, like in this example - the plots are taken from WolframAlpha:
And a different plot:
If $\delta>0$, we get a family of hyperbolas. WolframAlpha:
And a different plot:
In the first case, the function $F$ has minimum at the center of each of the ellipse. In the second case there is a saddle point. In both cases, the center fulfills
$$\frac{\partial F}{\partial x}= \frac{\partial F}{\partial y} =0$$
which leads to the system of equations
\begin{align*}
ax + by +d&= 0\\
bx + cy +e&= 0
\end{align*}
which is precisely the linear system described above.
So far we have only drawn some pictures (which is good for intuition, but picture is not a proof.) If we want to give a more rigorous reason why this works we can use some known facts about quadratic forms.
We know from Principal axis theorem that there is a rotation and translation such that in the new coordinates the curve has
$$\lambda_1 x^2 + \lambda_2 x^2=const$$
Clearly $F(x,y)=\lambda_1 x^2 + \lambda_2 x^2$ has partial derivatives equal to zero at the origin. (And depending on the signs of $\lambda_{1,2}$ we can say whether it has minimum, maximum or saddle point there.)
Best Answer
The condition $\Delta/I < 0$ is necessary for the solution set to exist at all. For example, consider $x^2 + y^2 + 1 = 0$, which has no solutions; here $\Delta = 1$ and $I = 2$. I think the necessity of this condition can be shown using linear algebra; leave a comment if you'd like to see the proof.
Unless I'm missing something, $\Delta \neq 0$ is redundant, as it is implied by $\Delta/I < 0$.