Elementary Solution of Kepler Problem (and a few other problems)

Abstract

We present a simple method to obtain the solution of a few orbital problems: the Kepler problem, the modified Kepler problem by the addition of an inverse square potential and linear force.

Keywords
classical mechanics; Kepler problem; celestial mechanics

1. Introduction

There is a handful of non trivial problems in classical mechanics one can solve exactly, such as the one-dimensional harmonic oscillator (1dHO) and the Kepler problem[¹, ², ³, ⁴, ⁵, ⁶, ⁷, ⁸, ⁹]. The latter one, that is, to find the orbit of a point mass under the action of a inverse square force, is particularly important for its historical role as well as applications in celestial mechanics.

Besides their important applications in the real world, these problems also serve as theoretical playground where one can apply different methods an techniques, with varying degrees of sophistication, as can be easily learned from analytical mechanics books [¹⁰, ¹¹].

In this note we present yet another solution by means of a very simple trick, using minimal technical background, avoiding cumbersome manipulations or involved mathematical calculations, making the approach suitable for beginners.

Initially we solve the 1dHO as a way to illustrate the method. Then we move to the Kepler problem, followed by the study of the Kepler problem modified by the addition of a $1/r^2$ potential and, at the end we discuss the orbit of a point mass under a quadratic potential which is, surprisingly, a bit more involved than the Kepler problem.

2. Warmup Problem: 1d Harmonic Oscillator

Consider a $1d$ harmonic oscillator, a point mass $m$ attached to a spring of constant $k$. This is a conservative system and we know its total mechanical energy, $E$, is conserved and is given by

(1) $$E=\frac{m}{2}{\stackrel{.}{x}}^2+\frac{k}{2}{x}^2.$$

We can rewrite equation (1) as

(2) $${\left(\sqrt{\frac{m}{2E}}\stackrel{.}{x}\right)}^2+{\left(\sqrt{\frac{k}{2E}}x\right)}^2=1.$$

Whenever the sum of the squares of two real quantities adds up to $1$ we can parametrize them by trigonometric functions, namely, we can write each term as

(3) $$\sqrt{{\displaystyle \frac{m}{2E}}}\stackrel{.}{x}=\sin \varphi ,$$

(4) $$\sqrt{{\displaystyle \frac{k}{2E}}}x=\cos \varphi .$$

where $\varphi$ is a function of time, to be determined. Taking the time derivative of equation (4) and using equation (3), we obtain

(5) $$\sqrt{\frac{k}{2E}}\stackrel{.}{x}=-\stackrel{.}{\varphi }\sin \varphi =-\stackrel{.}{\varphi }\sqrt{\frac{m}{2E}}\stackrel{.}{x},$$

which implies

(6) $$\stackrel{.}{\varphi }=-\sqrt{\frac{k}{m}}=-\omega ,$$

where we introduced $\omega =\sqrt{k/m}$. Therefore $\varphi =-\omega t+{\varphi }_0$ and the solutions for $x$ is

(7) $$x=\sqrt{\frac{2E}{k}}\cos (\omega t-{\varphi }_0).$$

This is the complete solution of the $1d$ harmonic oscillator. Notice that the amplitude is already expressed in terms of the total energy of the system. The only adjustable parameter is the phase ${\varphi }_0$, which can always be set to zero, by demanding our clock starts at the moment the particle is most distant from the origin.

3. The Kepler Problem

Explaining the motion of celestial bodies, in particular the orbits of comets and planets in the solar system, is one of the greatest triumphs of Newtonian mechanics. This is a problem that deserves to be studied carefully and from as many different angles as possible.

Let us treat the Kepler problem, that is, the motion of a point particle of mass $m$ orbiting a planet of mass $M$ under the central potential $V(r)=-GMm/r$.

Initially, we note that the angular momentum $\overrightarrow{L}=\overrightarrow{r}\times \overrightarrow{p}$ is conserved for central forces, for

(8) $$\frac{d}{dt}\overrightarrow{L}=\frac{d}{dt}\overrightarrow{r}\times \overrightarrow{p}+\overrightarrow{r}\times \frac{d}{dt}\overrightarrow{p}=\overrightarrow{v}\times \overrightarrow{p}+\overrightarrow{r}\times \overrightarrow{F}=0$$

since the velocity $\overrightarrow{v}$ is parallel to the momentum $\overrightarrow{p}$ and the position vector $\overrightarrow{r}$ is parallel to the central force $\overrightarrow{F}=F(r)\widehat{r}$. This general fact implies that the motion of a particle under a central force takes place in a plane.

We can, therefore, use polar coordinates to describe the motion of the planet, with $\overrightarrow{r}=r\widehat{r}$, and its velocity $\overrightarrow{v}=\stackrel{.}{r}\widehat{r}+r\stackrel{.}{\theta }\widehat{\theta }$. Therefore, the (conserved) mechanical energy $E$ is given by

(9) $$E=\frac{1}{2}m{v}^2+V(r)=\frac{1}{2}m({\stackrel{.}{r}}^2+r^2{\stackrel{.}{\theta }}^2)-\frac{GMm}{r}.$$

The angular momentum is given by $\overrightarrow{L}=\overrightarrow{r}\times \overrightarrow{p}=m{r}^2\stackrel{.}{\theta }\widehat{z}=L\widehat{z}$, where $L=m{r}^2\stackrel{.}{\theta }$ is a constant. Using that $\stackrel{.}{\theta }=L/m{r}^2$, the energy can be written as

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

We complete the squares, obtaining

(11) $$E=\frac{m}{2}{\stackrel{.}{r}}^2+\frac{L^2}{2m}{\left(\frac{1}{r}-\frac{GM{m}^2}{L^2}\right)}^2-\frac{G^2{M}^2{m}^3}{2L^2}.$$

Introducing $K\equiv E+G^2{M}^2{m}^3/2L^2$, we have

(12) $$K={\left(\sqrt{\frac{m}{2}}\stackrel{.}{r}\right)}^2+{\left(\frac{L}{\sqrt{2m}}\left(\frac{1}{r}-\frac{GM{m}^2}{L^2}\right)\right)}^2.$$

Since $K$ is the sum of two squares, it is a positive constant, unless $\stackrel{.}{r}=0$ and $1/r=GM{m}^2/L^2$, which corresponds to circular orbits. Therefore we can rewrite (12) as

(13) $${\left(\sqrt{\frac{m}{2K}}\stackrel{.}{r}\right)}^2+{\left(\frac{L}{\sqrt{2mK}}\left(\frac{1}{r}-\frac{GM{m}^2}{L^2}\right)\right)}^2=1.$$

Equation (13) allows us to write each of the squared terms as

(14) $$\sqrt{{\displaystyle \frac{m}{2K}}}\stackrel{.}{r}=\sin \varphi ,$$

(15) $${\displaystyle \frac{L}{\sqrt{2mK}}}\left({\displaystyle \frac{1}{r}}-{\displaystyle \frac{GM{m}^2}{L^2}}\right)=\cos \varphi$$

where $\varphi$ is a function of time, to be determined. Taking the time derivative of the second expression in equation (15) and equation (14), we obtain

(16) $$\frac{\stackrel{.}{r}}{r^2}=\frac{\sqrt{2mK}}{L}\stackrel{.}{\varphi }\sin \varphi =\frac{m}{L}\stackrel{.}{\varphi }\stackrel{.}{r},$$

implying $L=m{r}^2\stackrel{.}{\varphi }$. But since $L=m{r}^2\stackrel{.}{\theta }$, we conclude that $\varphi =\theta +{\varphi }_0$. Without loss of generality we can take ${\varphi }_0=0$, and solve for $r$ in equation (15),

(17) $$r=\frac{p}{1+e\cos \theta },$$

where $p=L^2/GM{m}^2$ and $e=\sqrt{2L^{2}K/G^2{M}^2{m}^3}$ is the eccentricity of the orbit. For completeness, we show that this is the equation of conic sections in polar coordinates.

We can rewrite equation (17) as

(18) $$r+er\cos \theta =p,$$

which implies

(19) $$\sqrt{x^2+y^2}+ex=p.$$

From this we readily derive that

(20) $$x^2\left(\frac{1-e^2}{e^2}\right)+\frac{2p}{e}x+\frac{y^2}{e^2}=\frac{p^2}{e^2},$$

that is, for it is an ellipse, for $e=1$ a parabola and for $e>1$ a hyperbole.

4. The Kepler Problem as a Harmonic Oscillator

We can rewrite Equation (12) in such a way that we see a direct connection with the one dimensional harmonic oscillator. Using $L=m{r}^2\stackrel{.}{\theta }$ is a constant, we have

(21) $$K={\displaystyle \frac{L^2}{L^2}}{\left(\sqrt{{\displaystyle \frac{m}{2}}}\stackrel{.}{r}\right)}^2+{\left({\displaystyle \frac{L}{\sqrt{2m}}}\left({\displaystyle \frac{1}{r}}-{\displaystyle \frac{GM{m}^2}{L^2}}\right)\right)}^2$$ $$={\left(L\sqrt{{\displaystyle \frac{m}{2}}}{\displaystyle \frac{1}{m{r}^2\stackrel{.}{\theta }}}\stackrel{.}{r}\right)}^2+{\left({\displaystyle \frac{L}{\sqrt{2m}}}\left({\displaystyle \frac{1}{r}}-{\displaystyle \frac{GM{m}^2}{L^2}}\right)\right)}^2$$ $$={\left({\displaystyle \frac{L}{\sqrt{2m}}}{\displaystyle \frac{1}{r^2}}{\displaystyle \frac{dr}{d\theta }}\right)}^2+{\left({\displaystyle \frac{L}{\sqrt{2m}}}\left({\displaystyle \frac{1}{r}}-{\displaystyle \frac{GM{m}^2}{L^2}}\right)\right)}^2.$$

Finally, introducing $u=1/r$, we have

(22) $${\left(\frac{du}{d\theta }\right)}^2+{\left(u-\frac{GM{m}^2}{L^2}\right)}^2=\frac{2mK}{L^2},$$

which is quadratic in $u$ and $\stackrel{.}{u}$, being, formally, the same energy function as the one-dimensional harmonic oscillator.

5. A Related Problem

Careful observations of the orbit of Mercury revealed it not to be an ellipse, as expected from a simple model of spherical celestial bodies [¹²]. During it’s orbit, Mercury’s perihelion, the point nearest to the sun, precesses by a small amount. This led the french mathematician Urbain Le Verrier to conjecture the existence of a small planet, Vulcan, as it was called, with an orbit between Mercury and the sun. A first correction to the gravitational potential is proportional to $1/r^2$. We can use the same method to solve this problem, namely, consider the potential

(23) $$V(r)=-\frac{GMm}{r}+\frac{\alpha }{r^2}$$

where $\alpha$ is a constant. In this case, the mechanical energy is given by

(24) $$E=\frac{1}{2}m{\stackrel{.}{r}}^2+\frac{L^2}{2m{r}^2}-\frac{GMm}{r}+\frac{\alpha }{r^2}=\frac{1}{2}m{\stackrel{.}{r}}^2+\frac{{\tilde{L}}^2}{2m{r}^2}-\frac{GMm}{r}$$

where ${\tilde{L}}^2=L^2+2m\alpha$. Following the calculations in section ³, in particular, equation (16) we find $\tilde{L}=m{r}^2\stackrel{.}{\varphi }$, and therefore

(25) $$L{\left(1+\frac{2m\alpha }{L^2}\right)}^{1/2}=m{r}^2\stackrel{.}{\varphi }$$

from which we derive $\varphi =\theta \sqrt{1+2m\alpha /L^2}$, that is, the angular variable $\theta$ is multiplied by a constant, which depends on the total angular momentum. The solution is, similarly to the Kepler problem,

(26) $$r=\frac{p}{1+e\cos (\theta \sqrt{1+2m\alpha /L^2})}$$

where $p={\tilde{L}}^2/GM{m}^2$ and $e=\sqrt{2{\tilde{L}}^{2}K/G^2{M}^2{m}^3}$, with $K\equiv E+G^2{M}^2{m}^3/2{\tilde{L}}^2$ is the excentricity of the orbit. Notice that the perihelion happens at ${\theta }_n\sqrt{1+2m\alpha /L^2}=2n\pi$, $n=0,1,\mathrm{\dots }$. The angular separation between two perihelions is $\mathrm{\Delta }{\theta }_n=2\pi /\sqrt{1+2m\alpha /L^2}$, which shows the perihelion precession explicitly.

6. Orbit Under Linear Force

We can use the same technique to solve the problem of the orbit of a point mass $m$ under the influence of a linear force, that is, a quadratic potential. The energy is given by

(27) $$E=\frac{1}{2}m{\stackrel{.}{r}}^2+\frac{1}{2}m{r}^2{\stackrel{.}{\theta }}^2+\frac{1}{2}k{r}^2,$$

where $k>0$. The case can be dealt in the same way or even obtained by taking $k\to -k$ in the final result. Since this is a central potential, the angular momentum $L=m{r}^2\stackrel{.}{\theta }$ is conserved, and we have

(28) $$E=\frac{1}{2}m{\stackrel{.}{r}}^2+\frac{L^2}{2m{r}^2}+\frac{1}{2}k{r}^2$$

We can complete the squares, like in the previous examples,

(29) $$E={\displaystyle \frac{1}{2}}m{\stackrel{.}{r}}^2+{\displaystyle \frac{L^2}{2m{r}^2}}+{\displaystyle \frac{1}{2}}k{r}^2$$ $$={\displaystyle \frac{1}{2}}m{\stackrel{.}{r}}^2+{\left({\displaystyle \frac{L}{\sqrt{2m}}}{\displaystyle \frac{1}{r}}\pm \sqrt{{\displaystyle \frac{k}{2}}}r\right)}^2\mp L\sqrt{{\displaystyle \frac{k}{m}}}.$$

Introducing $\omega =\sqrt{k/m}$ and $K_{\pm }=E\pm L\omega$, we have

(30) $${\left(\sqrt{\frac{m}{2K_{\pm }}}\stackrel{.}{r}\right)}^2+{\left(\frac{L}{\sqrt{2m{K}_{\pm }}}\frac{1}{r}\pm \sqrt{\frac{k}{2K_{\pm }}}r\right)}^2=1.$$

We can, therefore, write the terms inside the parentheses as $\sin {\varphi }_{\pm }$ and $\cos {\varphi }_{\pm }$, where ${\varphi }_{\pm }$ is a function yet to be determined,

(31) $$\sqrt{{\displaystyle \frac{m}{2K_{\pm }}}}\stackrel{.}{r}=\sin {\varphi }_{\pm };$$

(32) $${\displaystyle \frac{L}{\sqrt{2m{K}_{\pm }}}}{\displaystyle \frac{1}{r}}\pm \sqrt{{\displaystyle \frac{k}{2K_{\pm }}}}r=\cos {\varphi }_{\pm }.$$

Taking the time derivative of the second equation and using the first equation, we obtain

(33) $$-{\displaystyle \frac{L}{\sqrt{2m{K}_{\pm }}}}{\displaystyle \frac{1}{r^2}}\stackrel{.}{r}\pm \sqrt{{\displaystyle \frac{k}{2K_{\pm }}}}\stackrel{.}{r}$$ $$=-{\stackrel{.}{\varphi }}_{\pm }\sin {\varphi }_{\pm }=-{\stackrel{.}{\varphi }}_{\pm }\sqrt{{\displaystyle \frac{m}{2K_{\pm }}}}\stackrel{.}{r},$$

that is,

(34) $$-\frac{L}{m{r}^2}\pm \omega =-{\stackrel{.}{\varphi }}_{\pm }.$$

Using $L=m{r}^2\stackrel{.}{\theta }$ and integrating, we obtain

(35) $${\varphi }_{\pm }=\theta \mp \omega t\pm {\alpha }_{\pm }$$

where ${\alpha }_{\pm }$ are integration constants. We can choose $\theta =0$ at $t=0$. Moreover, we can also choose $\theta =0$ as the point in the orbit when $\stackrel{.}{r}=0$ (closest or farthest from the center of force), which corresponds to ${\varphi }_{\pm }=0$. All this together implies we can take ${\alpha }_{\pm }=0$ with the appropriate choice of coordinates and initial moment.

After some slightly tedious but straightforward algebra, we find

(36) $$r_{\pm }^2=\frac{2L^2{\omega }^2+(E^2-L^2{\omega }^2)({\sin }^{2}2\theta \pm {\sin }^{2}2\theta )}{2Ek+2k\sqrt{E^2-L^2{\omega }^2}\cos 2\theta }.$$

Let us analyze the solutions, starting with $r_{+}$,

(37) $$r_{+}^2={\displaystyle \frac{L^2{\omega }^2+(E^2-L^2{\omega }^2){\sin }^{2}2\theta }{k(E+\sqrt{E^2-L^2{\omega }^2}\cos 2\theta )}}={\displaystyle \frac{1}{k}}(E-\sqrt{E^2-L^2{\omega }^2}\cos 2\theta ).$$

Plugging Equation (37) in the energy expression for the linear orbit does not give a constant, therefore being a spurious solution which must be discarded.

A simpler way to convince oneself that Equation (37) is not a solution of our problem is by noticing a peculiar behaviour, namely, for $E>3L\omega$, the curvature at $\theta =0$, say, is negative, contradicting the fact that this is a central force: the motion can not be ‘curved’ towards outside the trajectory. One qualitative way to see that is by inspecting the formal case $L=0$. Here the polar equation describes two circles tangent at the origin, in an ‘8-shaped’ figure. Increasing $L$ it is natural to expect the negative curvature effect.

This leaves us with the ‘minus’ solution for the physical motion of the particle

(38) $$r(\theta )=\frac{L\omega }{\sqrt{Ek+\sqrt{E^2-L^2{\omega }^2}k\cos 2\theta }}.$$

7. Conclusions

We obtained the complete solutions of the $1d$ harmonic oscillator and the central potential problems that depend on the radial coordinate as $1/r$, $1/r^2$ and $r^2$, using minimal technical background. The parameters of the solutions expressed directly in terms of the physical quantities, such as energy and angular momentum. We hope students, specially beginners, find this an appealing approach to such beautiful problems.

Acknowledgments

The author would like to thank Reinaldo de Melo e Souza and two anonymous referees, for helping improve this paper considerably.

References

• ^[1]
  K. Symon, Mechanics (Addison-Wesley, Manlo Park, 1971), 3rd ed.
• ^[2]
  A.P. French, Newtonian Mechanics (Thomas Nelson & Sons LTD, London, 1971), 3rd ed.
• ^[3]
  D. Morin, Introduction to Classical Mechanics (Cambridge University Press, Cambridge, 2008).
• ^[4]
  A.C. Chen, Am. J. Phys. 55 , 250 (1987).
• ^[5]
  E.D. Davis, Am. J. Phys. 79 , 1246 (2011).
• ^[6]
  E. Vogt, Am. J. Phys. 64 , 392 (1996).
• ^[7]
  J. Milnor, Amer. Math. Montlhy 83 , 353 (1983).
• ^[8]
  A. Simha, Am. J. Phys. 89 , 1009 (2021).
• ^[9]
  W.G. Unruh, arXiv1803.06770 (2018).
• ^[10]
  A. Fasano and S. Marmi, Analytical Mechanics (Oxford University Press, Oxford, 2006).
• ^[11]
  N. Lemos, Analytical Mechanics (Cambridge University Press, Cambridge, 2018).
• ^[12]
  J.D. Wells, arXiv1106.1568v1 (2011).

Publication Dates

• Publication in this collection
  05 Jan 2022
• Date of issue
  2022

History

• Received
  12 Oct 2021
• Reviewed
  18 Nov 2021
• Accepted
  22 Nov 2021