How to actually plot the Kepler Problem?

Question

As far as I understand the Kepler problem tries to predict the motion of two bodies under the influence of gravity.

Here the trajectory of each particle is chosen to be:

$$\vec{r}_1 = \vec{R}+\dfrac{m_2}{M}\,\vec{r} \qquad \vec{r}_2 = \vec{R}-\dfrac{m_1}{M}\,\vec{r}$$

where $\vec{R} = \dfrac{m_1\,\vec{r}_1+m_2\,\vec{r}_2}{M}$ is the the common barycentre, $\vec{r} = \vec{r}_1 -\vec{r}_2$ the relative Distance and $M = m_1+m_2$ the total mass.

Now on the way to the solution some simplifications have been made. For instance motion is just locked to the $x$-$y$-plane. Subsequently angular momentum $\vec{L}$ is constant and just non-zero in the $z$-combonent: $L_z = x\,\dot{y}-\dot{y}\,x $. Energy is defined as usual: $E = \dfrac{1}{2}\,\mu\vec{\dot{r}}^2+\dfrac{1}{2}\,M\,\vec{\dot{R}}^2+V$, $\quad\mu = \dfrac{m_1\,m_2}{M}$ being the reduced mass and $V$ the potential. Respectively using polar coordinates: $x = r\,\cos\varphi \quad y = r\,\sin\varphi$ the equations simplifies further more: $L_z = \mu\,r^2\,\dot{\varphi} \quad E_r = \dfrac{1}{2}\,\mu\dot{r}^2+\dfrac{{L_z}^2}{2\,\mu\,r^2}+V(r)$.

At this point using differentials the equation of motion becomes:

$$\varphi-\varphi_0 = \int_{r_0}^r \dfrac{L_z}{r^2\,\sqrt{(E_r-V(r))\,2\,\mu-\dfrac{{L_z}^2}{r^2}}}\,\mathrm{dr}$$

During the Kepler problem the potential is predefined: $V(r) = -\dfrac{\gamma}{r}$. In this case the integral can be solved and rearranged:

$$r(\varphi) = \dfrac{\rho}{1+\varepsilon\cos{\varphi}}\quad \text{where} \quad \varepsilon = \sqrt{\dfrac{2\,E_r\,{L_z}^2}{\gamma^2\,\mu}+1},\quad \rho = \dfrac{{L_z}^2}{\mu\,\gamma}$$

Now here is what the actual question is. How do I plot this thing? Apparently using the equations above it follows: $$\vec{r_1} = \left(\begin{array}{c}x_1 \\y_1\end{array}\right) = R+\dfrac{m_2}{M}\,\left(\begin{array}{c}r(\varphi)\,\cos\varphi \\r(\varphi)\,\sin\varphi\end{array}\right)\qquad \text{$\vec{r}_2$ by analogy}$$

But how about initial conditions? For example if the trajectory of one object moving is set to be $\vec{r}_a = \left(\begin{array}{c}a_1+t^2 \\ a_2+t\end{array}\right)$ and the other $\vec{r}_b = \left(\begin{array}{c}b_1-2\,t \\ b_2-1\end{array}\right)$ it seems not matching the formula for $R_1$, $L_z$ and $E_r$ thus $\vec{r}_1$ and $\vec{r}_2$ the way it was clarified earlier. Probably because time is now important. In general I do not understand how $\quad \vec{r}_1, \quad \vec{r}_2\quad $ are calculated at all, because therefore $\vec{R}$ needs to be determined, but it depends on $\vec{r}_{1,2}$ itself.

Answer 1

Two body equation of motions

\begin{align*} & \ddot{r} -r{\dot\varphi }^{2}+{\frac { \left( m+M \right) G}{{r}^{2}}}=0\\ &\ddot\varphi +{2\, \frac {\dot\varphi \,{\dot r}}{r}}= 0 \end{align*}

solve those differential equations with the initial conditions \begin{align*} &\textbf{vis-viva -equation}\\\\ &v=\sqrt{G\,(M+m)\,\left(\frac{2}{r}-\frac{1}{a}\right)}\\ &\text{with}\\ &r=\frac{a\,(1-e^2)}{1+e\,\cos(\varphi)}\\ &\text{the value of r at the perigee}\\ &r_0=r(\varphi=0)=\frac{a\,(1-e^2)}{1+e}=a\,(1+e)\\ &v_0=v(r=r_0)=\sqrt {G \left( M+m \right) \left( 2\,{\frac {1+e}{a \left( 1-{e}^{2} \right) }}-\frac{1}{a} \right) } \\\\ &\text{hence}\\ &r(0)=r_{0}\,,\dot r(0)=0\\\\ &\varphi(0)=0\,,\dot \varphi (0)=\frac{v_0}{r_0} \end{align*}

you obtain the position of the bodies

\begin{align*} &x_1=\frac{\mu}{m}\,r(t)\\ &x_2=-\frac{\mu}{m}\,r(t)\\ &\text{with}\\ &\mu=\frac{m\,M}{m+M} \end{align*}

Answer 2

I'm not sure I understand the question... why are you trying to enforce that $\vec{r}_1$ and $\vec{r}_2$ are these arbitrary functions of time? You already have expressions for $\vec{r}_1$ and $\vec{r}_2$, and your initial conditions should enforce what the conserved quantities are.

If you specifically would like to plot the time-dependence, you should take the formula:

$E_r=\frac{1}{2}\mu\dot{r}^2+\frac{L_z^2}{2\mu r^2}+V(r)$

You can then manipulate this to

$\frac{dr}{dt}=\sqrt{2\left(\frac{E_r}{\mu}-\frac{L_z^2}{2\mu^2r^2}-\frac{V(r)}{\mu}\right)}$

Then the time can be related to the function $r(\varphi)$ (and thus $\varphi$) by the following integral, which can be done numerically:

$t=\int dt=\int\frac{dr'}{\sqrt{2\left(\frac{E_r}{\mu}-\frac{L_z^2}{2\mu^2r'^2}-\frac{V(r')}{\mu}\right)}}$