newtonian mechanics - Kepler problem in time: how do two gravitationally attracted particles move?

Two particles with initial positions and velocities $r_1,v_1$ and $r_2,v_2$ are interacting by the inverse square law (with G=1), so that

$${d^2r_1\over dt^2} = - { m_2(r_1-r_2)\over |r_1-r_2|^3}$$ $${d^2r_2\over dt^2} = - { m_1(r_2-r_1)\over |r_1-r_2|^3}$$

(the inverse square law along the line of separation). What is the complete solution of these differential equations? What is the position of the two objects as a function of time?

After reading a lot on Wikipedia, I've come to the definition of center of mass and relative coordinates:

$$R(t) ~=~ \frac{m_1 r_1 + m_2 r_2}{m_1+m_2}$$

$$\ddot{r}(t) r(t)^2 ~=~ (m_1+m_2)G$$

Where $R$ is the center of mass, and $r$ is the displacement between the particles... Is this correct? How do I proceed to solve the differential equation?

Answer

Textbook solution

The first thing to do is to define the center of mass and relative coordinates:

$$R(t) = {m_1 r_1 + m_2 r_2 \over m_1 + m_2}$$

$$r(t) = r_2 - r_1$$

You invert this to find

$$r_1= R - {m_2\over M} r$$ $$r_2=R+ {m_1\over M} r$$

The equation of motion for R is trivial, since center of mass is a conservation law:

$${d^2 R \over dt^2 } = 0$$

and it is solved by $$R(t) = V_0 t + R_0$$ where $V_0$ and $R_0$ are the initial center of mass velocity and position respectively (which are calculated from the given initial conditions).

The nontrivial equation is for the relative coordinate: $$m{d^2 r \over dt^2} = - {m_1 m_2 r\over |r|^3}$$ where m is the reduced mass: ${1\over m} ={1\over m_1} + {1\over m_2}$.

Or: $${d^2 r \over dt^2} = - {M r\over |r|^3}$$

Where $M=m_1 + m_2$ is the total mass.

The problem is reduced to solving the Kepler motion in a 1/r potential. From now on, I will rescale time to make the mass parameter in the r equation 1.

You can choose the x axis to lie along the initial r, and the y-axis to lie along the component of the initial $\dot{r}$ perpendicular to the initial r. Another way of saying this is that you rotate the coordinates to make the angular momentum vector $r\times p$ where $p=m\dot{r}$ to lie along the z-axis. This rotation reduces the problem to a plane, and the rotation matrix columns are given by the normalized initial r (now along the x-axis), the component of the initial velocity perpendicular to r, normalized (along the y-axis), and normalized L along the z-axis.

You then use units to set the total over reduced mass to 1, and use polar coordinates in the x-y plane of the motion, and note that the angular momentum is constant:

$$r^2 {d\theta\over dt} = L$$

This tells you that equal areas are swept out by the r-vector in equal times. The equation of motion for r(t) (no longer a vector, now a scalar radial coordinate) is:

$${d^2 r \over dt^2} = {L^2 \over r^3} - {1 \over r^2}$$

Then you change time out for $\theta$, expressing everything in terms of $r(\theta)$, which you can do using the equal area law, whenever the anguar momentum is nonzero (if the initial angular momentum is zero, or very close to zero, this is a one-dimensional two-body problem which can be solved directly by more elementary means). The equation of motion for $r(\theta)$ simplifies when you make a coordinate transformation to $u={1\over r}$:

$${d^2 u \over d\theta^2} = C - u$$

Where C is some unimportant constant, and this is solved by

$$u(\theta) = {1\over A} (1 + a \cos(\theta -\theta_0))$$

Where A is the semi-major axis of the ellipse (if the orbit is an ellipse), $\theta_0$ determines the orientation in the x-y plane, and a is the eccentricity of the ellipse (if a<1), or determines the angle of the hyperbola (if a>1) or tells you the orbit is a parabola (a=1).

the only result you need is that

$$r(\theta) = {A\over 1+ a\cos(\theta-\theta_0)}$$

This gives you the solution of r as a function of $\theta$, which gives the shape of the orbit. This is where textbooks stop.

Finding $\theta$ as a function of $t$

But you then want the solution for $\theta$ as a function of time, to get the r and $\theta$ as functions of time. This is determined from the area law, conceptually:

$$r^2 {d\theta\over dt} = L$$

$${d\theta \over (1+ a \cos(\theta-\theta_0))^2 }= {L \over A^2} dt$$

and integrating this from time 0 to time t, tells you in principle what $\theta(t)$ is. The result can be written as:

$$F(a,\theta) = F(a,\theta_0) + {L\over A^2} t$$

Where $F(a,\theta)$ is the special function that gives you the area of a conic section of parameter a, in a wedge from the focus where one half-line is along the major axis, and the other half-line makes an angle $\theta$ with the first. This special function is not expressible in terms of elementary functions.

This function is defined by the integral above, and you can calculate it numerically using any numerical integration method. Finding this function, and inverting it, is the only difficult part of this problem. There are three limits which are necessary for perturbations:

• for a=0, $F(0,\theta) = \theta$


• for a=1, $F(1,\theta) = {y\over 4} + {y^3\over 12}$ where $y=r\sin{\theta} = {\sin(\theta)\over 1+\cos(\theta)}$


• for $a=\infty$, $F(a,\theta) \approx {1\over a} \tan(\theta)$

Each of these are elementary degenerations: the first is the circle, the second is the parabola, and the third is a straight line. The important thing is that each of these degenerations gives you x(t) and y(t) which are simple, and further, you can perturb around each of these three limits in a nice way. In the following, the parameter t is rescaled to absorb ${L\over A^2}$

• circle: $x(t) = \cos(t)$ $y(t) = \sin(t)$



• parabola: $y(t)= (\sqrt{1+36t^2} + 6t)^{1\over 3} - (\sqrt{1+36t^2} - 6t)^{1\over 3}$ $x(t) = {1\over 2} - {y^2\over 2}$


• line: $x(t) = {1\over a}$, $y(t) = t$

The line and circle are obvious, the parabola is found by inverting the cubic for y as a function of t using the cubic equation.

Near the circle, time is periodic with the orbital period, which is the area inside the ellipse $aA^2$ divided by the area sweep-rate $L/2$. So you have once-winding function from a circle to a circle, which can always be written as a Fourier series with a linear term, which is found from the power series of the integrand in a, integrated term by term. Near the straight-line hyperbola, you can similarly perturb in a series, and the only interesting degeneration is the parabola. Near the parabola, the perturbation theory is a little more complicated.