33.4 Planetary Motion

The gravitational interaction between a planet and the sun is described by the inverse square central force law.

We assume that the planet is much lighter than the sun and so we can imagine that the sun is fixed in the center of the system, and the planet moves about it. Actually we can avoid this assumption by looking at motion relative to the Center of Mass of the system, but will not bother to do so.

We therefore consider the sun to be at the origin, with coordinates (0, 0, 0).

We choose our coordinates so that the planet is at the position (1, 0, 0) at initial time t₀, and that its motion at that time is in the xy plane, by which we mean the plane obeying z = 0.

Because the acceleration that the planet experiences always points toward the sun (origin) the planet will never leave that plane and we can ignore the z component of everything entirely.

It is an empirical fact that all the planets move to a first approximation in the same plane, so that even taking the effects of attractions from other planets into consideration can be described in two dimensions.

Physicists attack this problem by defining conserved quantities (energy, momentum and angular momentum), and using the values and properties of these to characterize the motion.

Our approach would have been considered an impossible brute force approach to the problem in the past, but now it is quite easy to implement, and provides a refreshing supplement to the standard treatment of the subject about which we recommend any standard text on mechanics and motion under gravity.

The actual behavior of planets was carefully observed by astronomers over centuries and was crisply summarized in Kepler's three laws, which are as follows:

1. The motion of planets and other bodies subject to the same force is in orbits that are "conic sections": ellipses or hyperbolae or in very special circumstances parabolas (all with the sun as a focus), or straight lines.

2. The area swept out per unit time in any orbit is constant.

3. There is a certain specific relation between the period of an elliptical orbit and a measure of its radius, which relation we will not discuss further.

On Conic Sections, Conservation of Angular momentum and Kepler's Second Law

We will confine ourselves here to showing how to integrate the equations of motion numerically on a spreadsheet for this problem, and how to chart the results. By doing so you can look at x or y or r as a function of time or at the orbit of the system.

We shall see that it is not significantly more difficult than handling a single first order differential equation. It differs from same in that it is a second order equation, and we have two dependent variables, x and y that will be functions of time t.

The differential equation that we will solve is

subject to initial conditions r(0) = (1, 0, 0) and where you choose p and q.

We use units for which MG is 1, for convenience, but you need not do this, nor need you start at (1, 0, 0).

To solve it we will devote column A to the variable t, B to x, C to y,

All we have to do is to insert initial conditions and conditions for change once in each of these and copy down and we will have our solution.

What then do we put in the various columns?

I would start the action in row 11 leaving the first 10 rows for notes, constants (which you can change later) initial conditions and your basic time interval d.

The initial conditions for each variable can then be entered into row 11 in the appropriate columns.

We begin by giving the simplest way to do all this:

Column A: the time column a11 = t₀, a12 = a11+b2. a13 = 2*a12-a11, copy a13 down.

Column B: the x column b11 = x₀, b12 = b11 + e11*($b$2), copy this down and into column C.

Column C: the y column c11 = y₀, the rest comes from column B.

Column D: the r column d11 =sqrt(b11*bb11+c11*c11), copy it down.

Column E: the x dot column: e11 = dx/dt(t₀), e12 = e11 + ($b$2)*g11, copy down and into column F.

Column F: the y dot column: f11 = dy/dt(t₀), the rest comes from column E.

Column G: the x double dot column g11 = - b11*($d11^3), copy down and into column H.

And that is it.

To see the orbit chart highlight columns B and C and do an xy scatter chart on them.

Exercise 33.2 Do this for different initial conditions and look at the orbit you get (make each column as long as you can. Vary d to get different accuracy).

Is greater accuracy possible?

The procedure described above is like a left hand rule, and is similarly inaccurate. With only slightly more effort you can get trapezoid quality accuracy. (And of course extrapolation is possible if you really want it.)

How?

Change the entries in columns B and C by changing b12 to

b12 = b11 + (e11+e12)*($b$2)/2

and copying it down and to the right into C.

Similarly change E and F by changing e13 to

e13 = e12 + ($b$2)*(2*g12-g11)

and copying it down and to the right into F.

You can do even better by making more complex changes.

To learn more you should study numerical methods.

Can this method fail?

Yes it will.

When the planet gets too close to the sun, the second derivative will get so large that the various changes in a single time interval will be very badly approximated. This will cause a big change in energy and the planet will jump to a very different orbit.

Use of a smaller time interval can alleviate this problem.