How did people get ellipses of Newton's equations of motion and gravitation? [duplicate]

It has become somewhat general knowledge that from Newton's laws of motion and his law of universal gravitation it is possible to verify that the trajectories of the planets are ellipses, i remember reading in the feynman's lectures on physics Vol.1 that one could use sucessive linear approximations on newton's laws to get a trajectory that visualy approximates an ellipse, and as the time intervals taken get smaller and smaller the approximations improve.

Nonetheless is there any other way to get something like the general equation of the ellipse or some parametrization of an ellipse out of newton's equations in a more rigourous way, instead of the usual way of iterative approximations? I imagine that there are probably problems in trying to solve sistems of many bodies, given the existence of the three body problem and its generalization the many body problem, so considering only two planets is there a way to obtain ellipses out of Newton's equations in a more mathematically rigorous way if for so to be accomplished it is needed other developments in classical mechanics such as variational principles i would be okay with that, but only if its truly needed.

I have realized that such a demonstration may be too long to be done in here, so i would gladly accept a reference to some book or site in wich there is such a demonstration.

marked as duplicate by sammy gerbil, Kyle Kanos, Qmechanic♦Jul 28 '17 at 10:45

Start from the total energy of the System (kinetic energy formulated in polar coordinates)

$E = \frac{1}{2}m(\dot{r}^2+r^2 \dot{\phi}^2) - \frac{GMm}{r}$

where $E= const.$, $m,M$ are the masses of the planets and $G$ is gravitational constant. If you would derive this Expression by time, you would obtain Newton's equation of Motion. Moreover, angular momentum $L$ is a conserved quantity, i.e.

$L = mrv_\phi = mr^2\dot{\phi}$.

Combining yields

$E = \frac{1}{2}m(\dot{r}^2+ \frac{L^2}{m^2r^2}) - \frac{GMm}{r}$

or

$\dot{r} = \frac{dr}{dt} = \sqrt{\frac{2E}{m}+\frac{2GM}{r}-\frac{L^2}{m^2r^2}}$

which can be integrated to

$\int \frac{dr}{\sqrt{\frac{2E}{m}+\frac{2GM}{r}-\frac{L^2}{m^2r^2}}}=\int dt = t+c;c=const.$

Now you can solve this integral by performing the Substitution $r = \frac{1}{u}$. After some calculus you will obtain an equation of the form

$u = \frac{1}{r} = A+Bcos(\Omega t)$

for parameters $A,B,\omega$. These Parameters have such values that an equation of an Ellipse is obtained.

Radial forces conserve angular momentum so motion is planar and $r^2\dot{\theta}$ is a constant, say $C$. From $\hat{\mathbf{r}}=\cos\theta\mathbf{i}+\sin\theta\mathbf{j},\,\hat{\mathbf{\theta}}=-\sin\theta\mathbf{i}+\cos\theta\mathbf{j}$ we have $\dot{\hat{\mathbf{r}}}=\dot{\theta}\hat{\mathbf{\theta}},\,\dot{\hat{\mathbf{\theta}}}=-\dot{\theta}\hat{\mathbf{r}}$ so $\dot{\mathbf{r}}=\dot{r}\hat{\mathbf{r}}+r\dot{\theta}\hat{\mathbf{\theta}}$ and $\ddot{\mathbf{r}}=(\ddot{r}-r\dot{\theta}^2)\hat{\mathbf{r}}+r^{-1}\frac{d}{dt}(r^2\dot{\theta})\hat{\mathbf{\theta}}$, but the last term vanishes by angular momentum conservation (as expected, since the acceleration should be radial). From $\dot{\theta}=Cu^2$ with $u=r^{-1}$ we have $\frac{d}{dt}=Cu^2\frac{d}{d\theta}$ so $\frac{dr}{dt}=-C\frac{du}{d\theta},\,\frac{d^2r}{dt^2}=-C^2u^2\frac{d^2u}{d\theta^2},\,\ddot{r}-r\dot{\theta}^2=-C^2u^2(\frac{d^2 u}{d\theta^2}+u)$. We can first restate the radial part of Newton's second law as the Binet equation, which for inverse-square forces simplifies to $\frac{d^2 u}{d\theta^2}+u=\frac{1}{\ell}$ (say). Rotating the axes allows us to set $u=\frac{1+e\cos\theta}{\ell}$ for $e\ge 0$, a conic section. Closed stable orbits require $e<1$, i.e. an ellipse