Download The Kepler Problem | Astronomy Computational Problem Set 1 and more Exercises Astronomy in PDF only on Docsity!
MASSACHUSETTS INSTITUTE OF TECHNOLOGY
Physics Department
Earth, Atmospheric, and Planetary Sciences Department
Astronomy 8.282J12.402J February 28, 2006
Computational Problem Set 1
The Kepler Problem
Optional
Suggested Procedure:
The differential equation of motion for the Kepler problem, in vector form, is:
d
2 �r
dt
2
GM ˆr
r
2
GM�r
r
3
In Cartesian coordinates we have:
d
2 x
dt
2
GM x
r
3
d
2 y
dt
2
GM y
r
3
As four, firstorder, coupled equations:
dx
dt
= v x
dy
dt
= v y
dv x
dt
GM x
r
3
dvy
dt
GM y
r
3
Equations (4) through (7) can be solved by a simple 4thorder RungeKutta integration scheme.
A few pages of text describing this method, taken from the Numerical Recipes book, are attached.
In the language of the RK4 program (Numerical Recipes’ RungeKutta program):
Y = the four independent variables (x, y, vx, vy )
X = the time t
N = 4 (8 for the case of the twobody problem)
H = the time step Δt
DERIVS = user supplied analytic derivatives, [right hand sides of (4) through (7)]
YOUT = are the incremented Y variables after time step H
Or, you can use the canned integrators in mathematical packages such as Matlab, Mathematica,
IDL, etc.
Before carrying out any integrations, it is convenient to cast the equations in dimensionless
form. One possible choice of scaling parameters is:
[length] → a (semimajor axis)
[mass] → M
[speed] → Ω K a
[time] → 1 /Ω K
where Ω K ≡ GM/r
3 .
The time steps can be made either fixed (e.g., some small fraction of 1 /Ω K ), or variable,
depending on how rapidly the motion is changing. One possible variable timestep scheme is:
�r �r 1 − �r 2 Δt = ξ
or ξ
�v 1 − �v v 2
depending on whether you are integrating the oneor twobody problem. The choice of ξ basically
determines the approximate number of integration steps that will be taken around an orbit. A
choice of ξ = 0. 01 should produce rather accurate results.
For all three parts, you should display your results graphically.
Part 1
For the case of an infinite central mass and an orbiting test mass, integrate a circular orbit and
an eccentric orbit. Carry out the integration for several orbits, in both the circular and eccentric
orbit cases, to verify that the orbits close. Check that energy and angular momentum are conserved
around the orbit.
Part 2
Compute the orbits of two objects of comparable (but unequal) mass orbiting their common
center of mass. Carry out the calculations for both circular and eccentric orbits.
Part 3
Compute the precessing orbit of the mass in Part 1, if the potential is of the form
GM
φ = −
(r − Rg )
where R g ≡ 2 GM/c
2 , is the Schwarzschild radius associated with mass M. This is a postNewtonian
potential suggested by Paczynski (1980; A&A, 88, 23) that approximates General Relativistic dy
namics in strong gravitational fields. Treat every other aspect of the problem with simple Newtonian
mechanics.
Consider cases where a/R g = 5, 10, and 50. Make the orbit eccentric enough to be able to
discern the precession, but avoid periastron distances near 3 Rg.
Expand the potential (or the force) in a series involving terms in R g ; retain only the first
postNewtonian term. Use this form of the potential to find an analytic expression for the angular
precession per orbit. Compare the analytic and numerical results.
Finally, for fun, try to produce a circular orbit with r � 3 Rg , and with r � 3 Rg , and contrast
the results.
