Solving $\textbf{r}''(t)=\frac{GM}{(r(t))^3}\textbf{r}(t)$

Solution 1:

Your equation is wrong (you forgot the minus sign). It should instead be $$\ddot{\mathbf{r}}=\frac{-GM}{r^3}\mathbf{r}$$ Or, $$\ddot{\mathbf{r}}=\frac{-GM}{r^2}\hat{\mathbf{r}} ~~~~~(1)$$ With that out of the way, let's get started. It is well known that in polar coordinates, $$\ddot{\mathbf{r}} =\left(\ddot{r} -r\dot{\theta }^{2}\right)\hat{\mathbf{r}} +( r\ddot{\theta } +2\dot{r}\dot{\theta })\widehat{\mathbf{e}_{\theta }}$$ However, since the gravitational force is entirely radial, in our case this is simply $$\ddot{\mathbf{r}} =\left(\ddot{r} -r\dot{\theta }^{2}\right)\hat{\mathbf{r}}$$ Then, returning to (1), $$\ddot{r} -r\dot{\theta }^{2}=\frac{-GM}{r^2}$$ Now here's the clever part: We'll make a change of variable $u \equiv \frac{1}{r}$. Therefore $$\dot{r}=\frac{\mathrm{d}r}{\mathrm{d}t}=\frac{-1}{u^2}\frac{\mathrm{d}u}{\mathrm{d}t}=\frac{-1}{u^2}\frac{\mathrm{d}u}{\mathrm{d}\theta}\frac{\mathrm{d}\theta}{\mathrm{d}t}=\frac{-1}{u^2}\frac{\mathrm{d}u}{\mathrm{d}\theta}\dot{\theta}$$ We now recall the fact that $ r\ddot{\theta } +2\dot{r}\dot{\theta }=0$ and also observe that $$r\ddot{\theta } +2\dot{r}\dot{\theta }=\frac{1}{r}\frac{\mathrm{d}}{\mathrm{d}t}(r^2\dot{\theta})=0$$ Therefore $r^2\dot{\theta}$ is a conserved quantity, and we make the definition $h \equiv r^2\dot{\theta}=\frac{1}{u^2}\dot{\theta}$. Therefore, $$\dot{r}=\frac{-1}{u^2}\frac{\mathrm{d}u}{\mathrm{d}\theta}\dot{\theta}=\frac{-1}{u^2}\frac{\mathrm{d}u}{\mathrm{d}\theta}\frac{h}{r^2}=-h\frac{\mathrm{d}u}{\mathrm{d}\theta}$$ So, $$\ddot{r}=\frac{\mathrm{d}}{\mathrm{d}t}\left(-h\frac{\mathrm{d}u}{\mathrm{d}\theta}\right)=-h^2u^2\frac{\mathrm{d}^2u}{\mathrm{d}\theta^2}$$ Then returning to (1) again, we can substitute in our new variables to obtain $$-h^2u^2\frac{\mathrm{d}^2u}{\mathrm{d}\theta^2}-\frac{1}{u}\left(hu^2\right)^2=-GMu^2$$ $$\frac{\mathrm{d}^2u}{\mathrm{d}\theta^2}+u=\frac{GM}{h^2}$$ This is a second order linear constant coefficient ODE for $u(\theta)$ which has the following solution: $$u(\theta)=A\sin(\theta)+B\cos(\theta)+\frac{GM}{h^2}=\frac{1}{r(\theta)}$$ And if we let $r(\theta)=r_{\text{max}}$ at $\theta=0$ and $r'(\theta)=0$ at $\theta=0$, then we find that $A=0$ and $B=\frac{1}{r_{\text{max}}}-\frac{GM}{h^2}$ Thus $$r(\theta)=\frac{1}{B\cos(\theta)+\frac{GM}{h^2}}$$ If we define $C \equiv \frac{h^2B}{GM} = \frac{h^2}{GMr_{\text{max}}}-1$, $$r(\theta)=\frac{h^2/GM}{1+C\cos(\theta)}$$ More work is needed to express this in terms of the semi minor axis $a$, the semi major axis $b$, and the eccentricity $\epsilon$, but I hope this is satisfactory enough. Finding $r(\theta)$ instead of $r(t)$ is considerably more useful, so I would recommend to stick to this sort of formulation. Please ask if anything is unclear!

THE PARAMETRIC FORM (if you're interested)

The parametric form of the ellipse is $$\mathbf{r}(t)=a\cos(t)\hat{\mathbf{i}}+b\sin(t)\hat{\mathbf{j}}$$ Which in polar coordinates is $$\mathbf{r}(t)=\sqrt{a^2\cos^2(t)+b^2\sin^2(t)}\hat{\mathbf{r}}+\operatorname{atan2}\left(b\sin(t),a\cos(t)\right)\widehat{\mathbf{e}_{\theta}}$$ Information on $\operatorname{atan2}$ can be found here

EDIT: A more general analysis:

Suppose our equation is $$\ddot{\mathbf{r}}=\frac{D}{r^3}\mathbf{r}$$ Where $D \in \mathbb{R}$ is some positive or negative constant. Then by similar analysis as before we reach $$r(\theta)=\frac{1}{B\cos(\theta)-\frac{D}{h^2}}$$ The solutions behave rather differently based on the sign of $D$. If it is negative, we get elliptical orbits as before. If $D$ is positive, however, the "orbits" (bit of an abuse of terminology, as they aren't periodic) are either parabolic or hyperbolic, depending on $h$, $r(0)$, and $r'(0)$.