Elliptical Motion in Classical Mechanics
Consider elliptical motion in three-dimensional space, typically that of a planet orbiting a large star such as the Sun. In classical mechanics, the energy (Hamiltonian) of elliptical motion is given by
\[ \small E = H(p,q) = \frac{p^2}{2m} – \frac{GmM}{r}. \]
Assume that one object (the fixed star) has a very large mass compared to the other object (the planet), and that the object (the fixed star) does not move from its origin. In reality, there is no need to assume that one of the masses is very large, and in problems that only consider the relative motion of two objects, it is acceptable to assume one of the coordinates is the origin without loss of generality. Let us find the equation of motion for the other object (planet) at this time. In the theory of relativity, this would correspond to the problem known as the Schwarzschild solution. Since the calculations would be lengthy, we will cover elliptical motion in the theory of relativity in a separate article.
For simplicity, we assume that the direction of rotation of the elliptical motion is in the \(\small xy\)-plane and that \(\small z=0\). If we calculate the acceleration, it becomes
\[ \small \begin{align*} &-\frac{\partial H}{\partial x}=F_x = m\frac{d^2x}{dt^2} = -\frac{GmM}{r^2}\frac{x}{r} \\ &-\frac{\partial H}{\partial y}=F_y = m\frac{d^2y}{dt^2} = -\frac{GmM}{r^2}\frac{y}{r} \\ &-\frac{\partial H}{\partial z}=F_z = m\frac{d^2z}{dt^2} = -\frac{GmM}{r^2}\frac{z}{r}=0. \end{align*} \]
Calculating the cross product of the coordinates and the force gives us
\[ \small q \times F = \left[ \begin{array}{c} yF_z-zF_y \\ zF_x-xF_z \\ xF_y-yF_x \end{array} \right] = \left[ \begin{array}{c} 0 \\ 0 \\ 0 \end{array} \right]. \]
This corresponds to the time derivative of angular momentum \(\small q\times p\), which is 0, and we can see that angular momentum is a constant (a conserved quantity). Calculating the cross product of force and angular momentum gives us
\[ \small \begin{align*} F \times l \; & = \left[ \begin{array}{c} F_yl_z -F_z l_y \\ F_zl_x -F_x l_z \\ F_xl_y -F_y l_x \end{array} \right] = -\frac{GmM}{r^3}\left[ \begin{array}{c} yl_z-zl_y \\ zl_x-x l_z \\ xl_y-y l_x \end{array} \right] \\ & = -\frac{GmM}{r^3}\left[ \begin{array}{c} x(xp_x+yp_y+zp_z)-(x^2+y^2+z^2)p_x \\ y(xp_x+yp_y+zp_z)-(x^2+y^2+z^2)p_y \\ z(xp_x+yp_y+zp_z)-(x^2+y^2+z^2)p_z \end{array} \right], \end{align*} \]
and the right-hand side corresponds to
\[ \small Gm^2M\frac{d}{dt}\left(\frac{q}{r} \right) = -Gm^2M\frac{qr\frac{dr}{dt}-\frac{dq}{dt}r^2}{r^3}, \quad q = x,y,z, \]
but since the angular momentum is a constant, we can integrate both sides to get
\[ \small p \times l = Gm^2M\left[ \begin{array}{c} \frac{x}{r} + \Lambda_x \\ \frac{y}{r} + \Lambda_y \\ \frac{z}{r} + \Lambda_z \end{array} \right]. \]
Where \(\small \Lambda\) is the integral constant.
Since
\[ \small \left[ \begin{array}{c} \Lambda_x \\ \Lambda_y \\ \Lambda_z \end{array} \right] = p\times l – \frac{Gm^2M}{r}\left[ \begin{array}{c} x \\ y \\ z \end{array} \right], \]
To be more specific, the calculation is:
\[ \small \left[ \begin{array}{c} \Lambda_x \\ \Lambda_y \\ \Lambda_z \end{array} \right] = \left[ \begin{array}{c} x(p_x^2+p_y^2+p_z^2)- p_x(xp_x + yp_y + zp_z) \\ y(p_x^2+p_y^2+p_z^2)- p_y(xp_x + yp_y + zp_z) \\ z(p_x^2+p_y^2+p_z^2)- p_y(xp_x + yp_y + zp_z) \end{array} \right]-\frac{Gm^2M}{r}\left[ \begin{array}{c} x \\ y \\ z \end{array} \right]. \]
Since \(\small \Lambda_x,\Lambda_y,\Lambda_z\) are constants, \(\small |\Lambda|^2= \Lambda^2_x+\Lambda^2_y+\Lambda^2_z\) is also a constant and we get
\[ \small \begin{align*} \Lambda^2\;&=G^2m^4M^2+\left(p_x^2+p_y^2+p_z^2- 2\frac{Gm^2M}{r}\right)(l_x^2+l_y^2+l_z^2) \\ & =G^2m^4M^2+2mEL^2 \end{align*} \]
by calculating it honestly. where \(\small E\) is the energy and \(\small L^2=l_x^2+l_y^2+l_z^2\) is the orbital angular momentum squared. By finding the angle from the dot product of the coordinates and the integral constant vector, we get
\[ \small \cos \omega = \frac{q_x\Lambda_x+q_y\Lambda_y+q_z\Lambda_z}{|\Lambda|r} = \frac{l_x^2+l_y^2+l_z^2}{|\Lambda|r}-\frac{Gm^2M}{|\Lambda|}, \]
and by calculating the orbital radius \(\small r\), we get
\[ \small r = \frac{\frac{l_x^2+l_y^2+l_z^2}{Gm^2M}}{1+\frac{|\Lambda|}{Gm^2M}\cos \omega} = \frac{\frac{L^2}{Gm^2M}}{1+\sqrt{1+\frac{2EL^2}{G^2m^3M^2}}\cos\omega}. \]
By defining
\[ \small \begin{align*} &a= -\frac{1}{2}\frac{GmM}{E} \\ &\epsilon = \sqrt{1+\frac{2EL^2}{G^2m^3M^2}}, \end{align*} \]
this equation can be transformed into
\[ \small r = \frac{a(1-\epsilon^2)}{1+\epsilon \cos \omega}. \]
\(\small a\) is called the semimajor axis (the longer radius of the ellipse), and \(\small \epsilon\) is called the eccentricity. By inverting the equation,
\[ \small \begin{align*} &E = -\frac{1}{2}\frac{GmM}{a} \\ &L^2 = Gm^2Ma(1-\epsilon^2) \end{align*} \]
holds true.
To identify the equations of motion, we fix the direction of rotation and assume \(\small z=0\). The coordinate in this case can be expressed as:
\[ \small \begin{align*} &x(t) = r(t)\cos(\omega (t)) \\ &y(t) = r(t)\sin(\omega (t)) \\ &z(t) = 0. \end{align*} \]
Note that unlike circular motion, \(\small r(t)\) is not a constant. Find the velocity and acceleration of this motion. For simplicity, let denote the notation:
\[ \small \dot{r} = \frac{dr(t)}{dt}, \quad \ddot{r}=\frac{d^2r(t)}{dt^2}, \quad \dot{\omega}=\frac{d\omega(t)}{dt}, \quad \ddot{\omega}=\frac{d^2\omega(t)}{dt^2}. \]
By calculating, we get
\[ \small \begin{align*} &\frac{dx}{dt} = \dot{r} \cos \omega(t) – r\dot{\omega} \sin \omega(t) \\ &\frac{dy}{dt} = \dot{r} \sin \omega(t) + r\dot{\omega} \cos \omega(t) \\ &\frac{d^2x}{dt^2} = (\ddot{r} – r\dot{\omega}^2) \cos \omega(t) – (r\ddot{\omega}+2\dot{r}\dot{\omega}) \sin \omega(t) \\ &\frac{d^2y}{dt^2} = (\ddot{r} – r\dot{\omega}^2) \sin \omega(t) + (r\ddot{\omega}+2\dot{r}\dot{\omega}) \cos \omega(t). \end{align*} \]
Since
\[ \small \begin{align*} &m\frac{d^2x}{dt^2} = -\frac{GmM}{r^2}\frac{x}{r} = -\frac{GmM}{r^2}\cos \omega(t) \\ &m\frac{d^2y}{dt^2} = -\frac{GmM}{r^2}\frac{y}{r} = -\frac{GmM}{r^2}\sin \omega(t) \end{align*} \]
holds, a correspondence relationship:
\[ \small \begin{align*} &m(\ddot{r} – r\dot{\omega}^2) = m\frac{d^2r}{dt^2} -mr \left(\frac{d\omega}{dt} \right)^2 = -\frac{GmM}{r^2} \\ &r\ddot{\omega}+2\dot{r}\dot{\omega} = r\frac{d^2\omega}{dt^2}+2\frac{dr}{dt}\frac{d\omega}{dt} =\frac{1}{r}\frac{d}{dt}\left(r^2 \frac{d\omega}{dt} \right) = 0 \end{align*} \]
can be found. These equations are known as the law of constant areal velocity (Kepler’s second law). The law of constant areal velocity is equivalent to the law of conservation of angular momentum in Newtonian mechanics. Calculating the acceleration with respect to the orbit radius in the usual way gives us
\[ \small \begin{align*} m\frac{d^2r}{dt^2} \; & = m\frac{d}{dt}\left(\frac{x\frac{dx}{dt}+y\frac{dy}{dt}}{r} \right) \\ &=\frac{m}{r}\left(x\frac{d^2x}{dt^2}+y\frac{d^2y}{dt^2}+\left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2 -\frac{\left(x\frac{dx}{dt}+y\frac{dy}{dt}\right)^2}{r^2} \right) \\ &=\frac{1}{r}\left(-\frac{GmM}{r}+m\frac{(x^2+y^2)\left(\left(\frac{dx}{dt}\right)^2+\left(\frac{dy}{dt}\right)^2\right)-\left(x\frac{dx}{dt}+y\frac{dy}{dt}\right)^2}{r^2} \right) \\ & = \frac{1}{r}\left(-\frac{GmM}{r}+\frac{\left(xp_y-yp_x\right)^2}{mr^2} \right). \end{align*} \]
The second term is the squared orbital angular momentum \(\small L^2\) and we can obtain
\[ \small m\frac{d^2r}{dt^2}-\frac{L^2}{mr^3} = -\frac{GmM}{r^2}. \]
Since
\[ \small m\frac{d^2r}{dt^2} -mr \left(\frac{d\omega}{dt} \right)^2 = -\frac{GmM}{r^2}, \]
holds, , we get
\[ \small \frac{d\omega(t)}{dt} = \frac{\sqrt{L^2}}{mr^2} = \frac{|L|}{mr^2} \]
by comparing the two equations.
Therefore, the final equation of motion is obtained as:
\[ \small \begin{align*} &r(t) = \frac{a(1-\epsilon^2)}{1+\epsilon \cos \omega(t)} \\ &x(t) = r(t)\cos \omega(t) \\ &y(t) = r(t)\sin \omega(t) \\ &z(t) = 0 \\ &\frac{d\omega}{dt} = \frac{|L|}{mr^2}. \end{align*} \]
Note that the angular momentum is calculated as:
\[ \small \begin{align*} &l_x = yp_z – zp_y = 0 \\ &l_y = zp_x – xp_z = 0 \\ &l_z = xp_y – yp_x = m\dot{\omega} r^2 \\ &L^2 = l_x^2+l_y^2+l_z^2 = m^2\dot{\omega}^2r^4. \end{align*} \]
The energy is calculated as:
\[ \small E = \frac{p^2}{2m} – \frac{GmM}{r}=\frac{L^2}{2mr^2} – \frac{GmM}{r}=-\frac{1}{2}\frac{GmM}{a} \]
from the orbital radius.
In addition, the area of the ellipse described by this elliptical motion is given by
\[ \small S = \pi a b = \pi a^2 \sqrt{1-\epsilon^2} \]
where the semimajor axis \(\small a\) and the semiminor axis \(\small b = a\sqrt{1-\epsilon^2}\). Since the areal velocity is:
\[ \small \frac{dS}{dt} = \frac{1}{2}r^2\dot{\omega} = \frac{|L|}{2m}, \]
the period of this elliptical motion \(\small T\) (the time it takes to complete one rotation and return to the original coordinates) can be calculated as:
\[ \small T = \frac{S}{dS/dt} = \frac{2m\pi a^2\sqrt{1-\epsilon^2}}{\sqrt{Gm^2Ma(1-\epsilon^2)}}=\frac{2\pi}{\sqrt{GM}}a^{3/2}. \]
It was the formula known as Kepler’s third law.
Elliptical Motion Rotating in Any Direction
Just like uniform circular motion, let us extend elliptical motion to rotation in any direction. All we need to do is multiply it by a rotational coordinate transformation matrix:
\[ \small \left[ \begin{array}{ccc} 1 & 0 & 0 \\ 0 & \cos \phi & \sin \phi \\ 0 & -\sin \phi & \cos \phi \end{array} \right] \left[ \begin{array}{ccc} \cos \theta & 0 & \sin \theta \\ 0 & 1 & 0 \\ -\sin \theta & 0 & \cos \theta \end{array} \right] \left[ \begin{array}{ccc} \cos \psi & \sin \psi & 0 \\ -\sin \psi & \cos \psi & 0 \\ 0 & 0 & 1 \end{array} \right]. \]
These represent rotations around the \(\small x,y,z\) axes, respectively. Around the \(\small z\) axis, all that changes is the value of the phase (just replace \(\small \omega (t)\) with \(\small \omega (t) + \psi)\), i.e. the initial coordinates, so if we ignore this and calculate, the result can be transformed into:
\[ \small \begin{align*} &x(t) = r(t)\cos\theta \cos(\omega (t))\cos \phi + r(t)\sin \phi \sin(\omega (t)) \\ &y(t) = r(t)\cos\theta \cos(\omega (t))(-\sin \phi) + r(t)\cos \phi \sin(\omega (t)) \\ &z(t) = r(t)\sin\theta \cos(\omega (t)). \end{align*} \]
The coordinate axes have been swapped so that the original equation results when \(\small \phi=\theta=0\). By changing the values of \(\small \theta,\phi\), it is possible to express elliptical motion rotating in any direction.
Find the velocity, momentum, and angular momentum. Note that unlike uniform circular motion, \(\small r,\omega\) is not a constant. The velocity can be calculated as:
\[ \small \begin{align*} &v_x =(\dot{r}(t)\cos(\omega (t)) -r(t)\dot{\omega}(t)\sin(\omega (t)))\cos\theta\cos \phi \\ & \qquad + (\dot{r}(t)\sin(\omega (t))+r(t)\dot{\omega}(t)\cos(\omega (t)))\sin \phi \\ &v_y = (\dot{r}(t) \cos(\omega (t))-r(t)\dot{\omega}(t)\sin(\omega (t)))\cos\theta(-\sin \phi) \\ & \qquad + (\dot{r}(t)\sin(\omega (t))+r(t)\dot{\omega}(t)\cos(\omega (t)))\cos \phi \\ &v_z = (\dot{r}(t)\cos(\omega (t))-r(t)\dot{\omega}(t)\sin(\omega (t)))\sin\theta \\ &v = \sqrt{v_x^2+v_y^2+v_z^2} = \sqrt{r^2(t) \dot{\omega}^2(t)+\dot{r}^2(t)}. \end{align*} \]
The momentum can be calculated using \(\small p=mv\). Angular momentum is calculated as :
\[ \small \begin{align*} &l_x = yp_z – zp_y = m\dot{\omega}(t) r^2(t)\sin\theta\cos \phi \\ &l_y = zp_x – xp_z = m\dot{\omega}(t) r^2(t)\sin\theta\sin \phi \\ &l_z = xp_y – yp_x = m\dot{\omega}(t) r^2(t)\cos\theta \\ &L^2 = l_x^2+l_y^2+l_z^2 = m^2\dot{\omega}^2(t)r^4(t). \end{align*} \]
As mentioned in the previous section, in elliptical motion, angular momentum is a constant, so
\[ \small r^2(t)\dot{\omega}(t) = \text{const.} \]
holds.
Remaining Issues
You might think that we have now calculated all the equations of motion for elliptical motion, but in fact you will notice that we have not calculated the specific functional form of \(\small \omega(t)\) or \(\small r(t)\) which is derived from it. It cannot be written as:
\[ \small \omega(t) = \bar{\omega} t \]
as in uniform circular motion. Generally, in elliptical motion, the closer to the focal point the faster the speed tends to be, and the farther from the focal point the slower the speed tends to be. It is necessary to find a functional form for \(\small \omega(t)\) that satisfies this property.
\(\small \omega(t)\) satisfies differential equation:
\[ \small \frac{d\omega(t)}{dt} = \frac{|L|}{mr^2(t)}, \]
and since
\[ \small \begin{align*} &r(t) = \frac{a(1-\epsilon^2)}{1+\epsilon \cos \omega(t)} \\ &L^2 = Gm^2Ma(1-\epsilon^2) \ \end{align*} \]
hold,
\[ \small \frac{d\omega(t)}{dt} = \sqrt{GM}\frac{(1+\epsilon \cos \omega(t))^2}{a^{3/2}(1-\epsilon^2)^{3/2}} \]
can be obtained. All we need to do is solve this differential equation, but is it possible to find a solution to this differential equation? Well, in this day and age of IT advances, we could just use a computer to do the calculations…This needs to be determined from a slightly different equation, and the solution can be found by solving an algebraic equation called the Kepler equation. As mentioned in the article on Bessel functions, we will use Bessel functions in this case. We will consider Kepler’s equations in more detail in the next article.
Comments