Python coding: Use of Verlet method
Let's use the Verlet method to calculate the orbit of the Earth around the Sun. The equations of motion for the position r = (x,y) of the planet in its orbital plane are the same as those for any orbiting body and are d^2x/dt^2 = -GMx/x^3 and d^2y/dt^2 = -GMy/y^3.
That is, the general form is:
d^2x/dt^2 = -GMx/x^3
d^2y/dt^2 = -GMy/y^3
The orbit is not perfectly circular, with the planet being sometimes closer to and sometimes further from the Sun. When it is at its closest point, or perihelion, it is moving precisely tangentially (i.e., perpendicular to the line between itself and the Sun) and it has a distance of 1.4710 * 10^11 m from the Sun and a linear velocity of 3.0287 * 10^4 m/s.
a) Write a program to calculate the orbit of the Earth using the Verlet method with a time-step of h = 1 hour. Make a plot of the orbit, showing several complete revolutions about the Sun. The orbit should be very slightly, but visibly, non-circular.
b) The gravitational potential energy of the Earth is -GMm/r, where m = 5.9722 * 10^24 kg is the mass of the planet, and its kinetic energy is (1/2)mu^2 as usual. Modify your program to calculate both of these quantities at each step, along with their sum (which is the total energy), and make a plot showing all three as a function of time on the same axes.
c) Now plot the total energy alone without the others, and you should be able to see a slight variation over the course of an orbit. However, energy should always return to its starting value at the end of each complete orbit, as energy is conserved in the long term.
Hints: Use the following equation for the Verlet method:
v(t+h) = v(t) - h * f(r(t), t) / 2
and, repeated application of:
r(t+h) = r(t) + h * v(t+h) + k * h^2 * f(r(t+h), t+h)
v(t+h) = v(t) + h * (f(r(t), t) + f(r(t+h), t+h)) / 2