To derive the polar coordinates of an ellipse, I believe it is simpler to start from a geometric perspective rather than an algebraic one. If you begin with the equation
$$
\frac{x^2}{a^2}+\frac{y^2}{b^2}=1
$$
you will encounter a quadratic equation. While this approach will ultimately yield the correct result, it involves quite a bit of computation. A geometric approach, on the other hand, simplifies the process.
Let us consider the polar form centered at the left focus of a horizontal ellipse.

Define the major axis length of the ellipse as $2a$. Let $P$ be an arbitrary point on the ellipse, and let $F_1$ and $F_2$ be the foci. Denote $PF_1 = r$ and $PF_2 = 2a - r$. Define the angles $\angle PF_1F_2 = \theta$ and $\angle PF_2F_1 = \phi$.
By the Law of Sines, we have
$$
\frac{r}{\sin\phi}=\frac{2a-r}{\sin\theta}
$$
Rearranging gives
$$
\sin\phi=\frac{r\sin\theta}{2a-r}
$$
Using the identity $\cos^2\phi+\sin^2\phi=1$, we obtain
$$
\cos^2\phi+\frac{r^2\sin^2\theta}{(2a-r)^2}=1
$$
which simplifies to
$$
\cos\phi=\sqrt{1-\frac{r^2\sin^2\theta}{(2a-r)^2}}
$$
Next, we use the equation
$$
r\cos\theta+(2a-r)\cos\phi=2c
$$
Substituting $\cos\phi$ with $\sqrt{1-\frac{r^2\sin^2\theta}{(2a-r)^2}}$, we get
$$
r\cos\theta+(2a-r)(\sqrt{1-\frac{r^2\sin^2\theta}{(2a-r)^2}})=2c
$$
Rearranging,
$$
2c-r\cos\theta=(2a-r)(\sqrt{1-\frac{r^2\sin^2\theta}{(2a-r)^2}})
$$
Squaring both sides, we obtain
$$
(2c-r\cos\theta)^2=(2a-r)^2-r^2\sin^2\theta
$$
$$
\implies r^2\cos^2\theta-4cr\cos\theta+4c^2=4a^2-4ar+r^2-r^2\sin^2\theta
$$
$$
\implies r^2(\cos^2\theta+\sin^2\theta-1)+r(-4c\cos\theta+4a)=4a^2-4c^2
$$
It is satisfying to see that the term $\cos^2\theta+\sin^2\theta-1=0$, eliminating $r^2$. This simplifies to
$$
r(-4c\cos\theta+4a)=4a^2-4c^2
$$
$$
\implies r=\frac{4a^2-4c^2}{-4c\cos\theta+4a} =\frac{a^2-c^2}{a-c\cos\theta}
$$
Using $c=a\epsilon$, we express this in terms of the eccentricity:
$$
r=\frac{a^2-c^2}{a-c\cos\theta}=\frac{a^2-a^2\epsilon^2}{a-a\epsilon\cos\theta}=\frac{a(1-\epsilon^2)}{1-\epsilon\cos\theta}
$$
This is the form when centered at the left focus. If instead we center at the right focus and follow the same steps, we obtain
$$
r=\frac{a(1-\epsilon^2)}{1+\epsilon\cos\theta}
$$
This geometric approach also extends naturally to hyperbolas and parabolas. Its main advantage is that it avoids solving a complicated quadratic equation. I hope this helps!