Prove that the locus of the vertices of the right circular cones that pass through the ellipse $\frac{x^2}{a^2}+\frac{y^2}{b^2}=1, z=0$ is $\frac{x^2}{a^2-b^2}-\frac{z^2}{b^2}=1, y=0$ or $\frac{y^2}{a^2-b^2}+\frac{z^2}{a^2}=-1, z=0$.
EDIT:
Here $\frac{x^2}{a^2}+\frac{y^2}{b^2}=1, z=0$ is the base of the cone. Let the vertex of the cone be $(x_1,y_1,z_1)$. Let the generator be $$\frac{x-x_1}{l}=\frac{y-y_1}{m}=\frac{z-z_1}{n}.$$
Then $z=0$ implies any point on ellipse be $(x_1-lz_1/n, y_1-mz_1/n,0)$. It lies on the ellipse $\frac{x^2}{a^2}+\frac{y^2}{b^2}=1$ then the point satisfies the ellipse we get $$\frac{(nx_1-lz_1)^2}{a^2}+\frac{(ny_1-mz_1)^2}{a^2}=n^2$$
Eliminating $l,m,n$ we get the equation of the cone as $$\frac{1}{a^2}(zx_1-xz_1)^2+\frac{1}{b^2}(zy_1-yz_1)^2=(z-z_1)^2.$$ How to get the locus of vertex $(x_1,y_1,z_1)$?
Edit 2
How to get the locus of the vertices in the given two forms using purely mathematical way? I not able to solve the problem. Please help.

It appears that your solution has stalled because you haven't incorporated the fact that the cone is circular. I'll start from scratch.
As you did, we take the cone's generator to pass through $(x_1, y_1, z_1)$ (with $z_1 \neq 0$), and to have unit direction vector $(\ell, m, n)$, so that we can write
$$\frac{x-x_1}{\ell} = \frac{y-y_1}{m} = \frac{z-z_1}{n} = t \tag{1}$$
As the vector has unit length, we have $$1 = \ell^2 + m^2 + n^2 \quad\to\quad (x-x_1)^2 + (y-y_1)^2 + (z-z_1)^2 = t^2 \tag{2}$$ To ensure circularity, we impose the condition that the generator makes a constant (non-right) angle with some unit axis vector, $(p,q,r)$. (We may assume $r\neq 0$, so that the axis is not parallel to the $xy$-plane; otherwise, the intersection with that plane would be a hyperbola.) For this, we need only assume that the dot product of the vectors is some (non-zero) constant: $$(p, q, r)\cdot (\ell,m,n) = k \neq 0 \quad\to\quad p (x-x_1) + q (y-y_1) + r (z-z_1) = kt \tag{3}$$
Eliminating the parameter $t$ from $(1)$ and $(2)$, by expressing $k^2t^2$ in two ways, gives the equation of a circular cone. $$k^2\left(\;(x-x_1)^2 + (y-y_1)^2 + (z-z_1)^2\;\right) = \left(\;p (x-x_1) + q (y-y_1) + r (z-z_1)\;\right)^2 \tag{4}$$
The intersection of this cone with the plane $z=0$ has this equation: $$\begin{align} 0 &= x^2 (k^2 - p^2) + y^2 ( k^2 - q^2 ) - 2 x y p q \\ &-2 x (\; (k^2-p^2) x_1 - p q y_1 - p r z_1\;) \\ &-2 y (\; (k^2-q^2) y_1 - p q x_1 - q r z_1\;) \\ &+ (k^2-p^2) x_1^2 + (k^2-q^2)y_1^2 + (k^2-r^2)z_1^2 - 2 p q x_1 y_1 \\ &- 2 p r x_1 z_1 - 2 q r y_1 z_1 \tag{5} \end{align}$$
For this to match the target ellipse, the equation must be a scalar multiple, say, with factor $u \neq 0$, of $x^2 b^2 + y^2 a^2 - a^2 b^2 = 0$. That is: $$\begin{align} k^2 - p^2 &= u b^2 \tag{6.1}\\ k^2 - q^2 &= u a^2 \tag{6.2}\\ p q &= 0 \tag{6.3}\\ ub^2 x_1 &= p r z_1 \tag{6.4}\\ ua^2 y_1 &= q r z_1 \tag{6.5}\\ ub^2 x_1^2 + ua^2 y_1^2 - (k^2-r^2)z_1^2 &= u a^2 b^2 \tag{6.6} \end{align}$$ where subsequent equations conveniently incorporate substitutions from previous ones.
From here, eliminating $u$, $k$, $p$, $q$ from the system $(6.x)$ results in the target relations, but the process isn't particularly fun. Things simplify a bit by observing from $(6.3)$ that one of $p$ and $q$ vanishes; by $(6.4)$ or $(6.5)$, this implies that $x_1$ or $y_1$, respectively, also vanishes. (This is how we get the separate loci.) Even so, we can make some interesting progress before splitting into cases.
For instance, it's easy enough to solve $(6.1)$ and $(6.2)$ to get (for $a\neq b$) $$u = \frac{p^2- q^2}{a^2-b^2} \qquad\qquad k^2 = \frac{a^2 p^2 - b^2 q^2}{a^2-b^2} \tag{7}$$ Then, from $(6.4)$ and $(6.5)$, we calculate thusly, $$( p^2 - q^2 )r^2 z_1^2 = u^2 (b^4 x_1^2 - a^4 y_1^2) = \frac{( p^2 - q^2 )^2}{(a^2-b^2)^2}\,(b^4 x_1^2 - a^4 y_1^2) \tag{8}$$ so that, dividing-through by $p^2-q^2$ (which safely assumes that one of $p$ and $q$ is non-zero when $a\neq b$), we have $$r^2 z_1^2 = \frac{u(b^4x_1^2-a^4y_1^2)}{a^2-b^2} \tag{9}$$
Therefore, we may re-write $(6.6)$ ... $$\begin{align} ub^2 x_1^2 + ua^2 y_1^2 + r^2 z_1^2 - u a^2 b^2 &= k^2z_1^2 \tag{10.1}\\[6pt] ub^2 x_1^2 + ua^2 y_1^2 + \frac{u(b^4x_1^2-a^4y_1^2)}{a^2-b^2} - u a^2 b^2 &= k^2z_1^2 \tag{10.2}\\[6pt] u a^2 b^2\left( \frac{x_1^2 - y_1^2}{a^2-b^2} - 1 \right) &= k^2z_1^2 \tag{10.3}\\[6pt] \frac{x_1^2 - y_1^2}{a^2-b^2} - \frac{k^2z_1^2}{ua^2b^2} &= 1 \tag{10.4} \\[6pt] \frac{x_1^2 - y_1^2}{a^2-b^2} - \frac{z_1^2}{a^2b^2}\frac{a^2p^2-b^2q^2}{p^2-q^2} &= 1 \tag{10.5} \end{align}$$ Finally, we can invoke the cases,
Given how elegant the geometric solution is, I suspect there's a more-clever route through the algebra here, but this is the best I have at the moment.