Planetary Motion

Physical Model

The motion of planets around the Sun is a classic problem in physics that beautifully demonstrates Newton’s laws of motion and universal gravitation. While this motion might seem very different from a spring pendulum, the mathematical description is surprisingly similar. The main difference is that instead of a spring force, we now have gravity as our central force.

In both cases, we describe the motion using two coordinates: the distance from the center (\(r\)) and an angle (\(\theta\)). For planetary motion, \(r\) is the distance between the planet and the Sun, and \(\theta\) describes the planet’s angular position in its orbit.

The equations of motion contain two parts: The first equation describes the radial acceleration (\(\ddot{r}\)), and the second describes the angular acceleration (\(\ddot{\theta}\)):

\[\begin{eqnarray} \ddot{r}&=&r\dot{\theta}^2-\frac{G\, M}{r^2}\\ \ddot{\theta}&=&-\frac{1}{r}2\dot{r}\dot{\theta} \end{eqnarray}\]

In the first equation, the term \(r\dot{\theta}^2\) represents the centrifugal effect - the tendency of the planet to move away from the Sun due to its orbital motion. The term \(-\frac{G\, M}{r^2}\) is Newton’s gravitational force (divided by mass), pulling the planet toward the Sun. \(G\) is the gravitational constant, and \(M\) is the mass of the Sun.

The second equation describes how the angular motion changes. The term \(2\dot{r}\dot{\theta}\) represents the coupling between radial and angular motion - as the planet moves closer to or farther from the Sun, its angular velocity must change to conserve angular momentum, similar to how an ice skater spins faster when pulling their arms in.

The solution to these equations gives us the famous orbital equation:

\[ r(\theta)=\frac{p}{1+\epsilon \cos(\theta)} \]

This is the equation of a conic section, where \(p\) and \(\epsilon\) determine the orbit’s shape. The parameter \(p\) is related to the angular momentum \(L\):

\[\begin{equation} p=\frac{L^2}{G M m^2} \end{equation}\]

The eccentricity \(\epsilon\) depends on both the energy \(E\) and angular momentum \(L\):

\[\begin{equation} \epsilon=\sqrt{1+\frac{2\frac{E}{m}\frac{L^2}{m^2}}{G^2M^2}} \end{equation}\]

When \(0 < \epsilon < 1\), we get an elliptical orbit (the case for all planets in our solar system). A perfect circle occurs when \(\epsilon = 0\). These equations, derived by Newton and studied by Kepler, explain not only planetary orbits but also the paths of comets, artificial satellites, and many other celestial objects.

Numerical Solution

For the numerical solution, we will use the solve_ivp function from the scipy.integrate module. This function solves ordinary differential equations (ODEs) given an initial state, time span, and optional time points for evaluation. We will integrate the equations of motion for a planet with an initial radius \(r_0\), radial velocity \(v_0\), angle \(\theta_0\), and angular velocity \(\omega_0\).

Initial Parameters: Planets

We first define the initial parameters for the planet’s motion. We set the mass of the Sun \(M=1\) and the mass of the planet \(m=1\) for simplicity. The gravitational constant \(G\) is set to \(4\pi^2\) to simplify the equations. We also define the initial radius \(r_0\), radial velocity \(v_0\), angle \(\theta_0\), and angular velocity \(\omega_0\).

Solution: Planets

We now solve the equations of motion using a higher accuracy solver method.

Plotting: Planets

We can now plot the numerical solution of the planetary motion. The black dashed line represents the analytical solution for the orbit of the planet. The black solid line shows the numerical solution obtained by integrating the equations of motion.

Trajectory

Energy

Finally, we can plot the total energy of the planet as a function of time. The total energy is the sum of the kinetic and potential energy of the planet. We can see that the energy is conserved over time, as expected for a system with no external forces.