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.

  • 1
    $\begingroup$ What is the problem? Can't you olot plot r$_i$ as a function of $\phi$, which a function of t? $\endgroup$
    – my2cts
    Nov 7, 2021 at 11:41

2 Answers 2


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*}

enter image description here


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


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)}}$


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.