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 t0, 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.

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.

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.