N体问题

Mission: Solve the "$N$-Body Problem"!

Calling all physicists, astronomers, mathematicians, and aspiring students! We are given an extremely tough mission to accomplish: solve the $N$-Body Problem. A problem that is worthy of being called one of the most difficult, yet important trials that we have faced. The great Issac Newton left this problem in his legacy to all of us, yet after many years none have successfully solved it. In this lesson, we will go through understanding the problem and its significance, a couple of variants of the problem (and one or two respective solutions) where $n$ has a clear numerical value (We have one other variant included as a special case that will be mentioned briefly), and attempts to solve the problem with no numerical restriction.

WARNING:

1.) This will be a considerably long lesson, so please make sure you have time to cover all the material. It is preferable that for any lesson, that you go through it in one sitting so that there is less of a risk of forgetting information.

2.) This lesson will involve levels of mathematics and physics that are at a very high level (which you can problably tel from the fact that this problem hasn't been solved yet), so please be forewarned that it may be very difficult to understand this lesson unless you have previously studied physics or mathematics with relation to astronomy (relativity is included in this criteria).

Stage 1: Understanding the Mission Statement

It is imperative that before we get into this mission any further, that we begin by understanding what we need to accomplish. So now we ask ourselves what is the $n$-body problem?

The $n$-body problem is a probem that asks us to create a mathematical model that can predict the individual motions of a group of celestial objects interacting with each other gravitationally. Over its many years of existence, solving this problem has been motivated by the desire to understand the motions of the Sun, Moon, planets and the visible stars.

If we are to get a more informal problem statement, it would be as follows:

Given the quasi-steady orbital properties (instantaneous position, velocity and time) of a group of celestial bodies, predict their interactive forces; and consequently, predict their true orbital motions for all future times. (Wikipedia)

Or in other words,

Given $N$ particles which interact with one another only via Newtonian inverse-square gravity, and some intitial configuration, predict the time-evolution of each particle’s position and velocity. (Senchyna, Peter 2013)

Stage 2: Taking the First Steps by Analyzing Special Cases/Variants of the $N$-Body Problem

Case 1.) Two Bodies

The two-body problem can be simply defined as determining the motion of two bodies orbiting one another under their mutual gravitational attraction. The following presents one of the solutions to this problem by reducing the two body problem into an equivalent one-body problem and then solved as a special case of motion in a potential field.

(Source: Drexel University)

The equations of motion for two bodies of masses $m_1$ and $m_2$, positions $x_1$ and $x_2$, and velocities $v_1$ and $v_2$ moving under gravity are

(1) $$m_1a_1 = \frac{Gm_1m_2}{r^3_{12}} (x_2 ? x_1)$$

(2) $$m_2a_2 = \frac{Gm_1m_2}{r^3_{12}} (x_1 ? x_2)$$

where $a = \ddot{x}$ and $r_12 = \mid x_1 ? x_2 \mid$. Using these equations, we can note that $m_1a_1 + m_2a_2 = 0$, so after integrating once, $m_1v_1 + m_2v_2 =$ constant and the center of mass of the bodies moves at a constant velocity. Dividing Equation (1) by $m_1$, Equation (2) by $m_2$, and subtracting, we get the equation of motion for the relative separation $r = x_2 ? x_1$:

(3) $$\ddot{r} = - \frac{GMr}{r^3_{12}}$$

where $M = m_1 + m_2$. Thus the relative motion in the two-body problem is identical to motion of a test particle in the potential field of a mass $M$. Momentum in a central potential is planar, since angular momentum is conserved, and in polar coordinates $(r, \theta)$ in that plane, the equations of motion are

(4) $$\ddot{r} - r\dot{\theta}^2 = F(r)$$

(5) $$r^2\dot{\theta}= L = \text{constant}$$

where $F(r) = -\varphi '(r)$ and the second expression simply states conservation of angular momentum. We can use the angular momentum integral (5) to simplify the radial equation (4) by eliminating the transverse motion $\dot{\theta}$:

(6) $$\ddot{r} = \frac{L^2}{r^3} + F(r) = \frac{d\varphi_{\text{eff}}}{dr}$$

where the effective potential $\varphi_{\text{eff}} = \varphi + \frac{1}{2} \frac{L^2}{r^2}$ combines the external and the centrifugal forces intoa single quantity. We can also use Equation (5) to eliminate time as a variable in Equation (4) by writing

(7) $$\frac{d}{dt} = \frac{L}{r^2} \frac{d}{d\theta}$$

so the radial equation becomes

(8) $$\frac{L^2}{r^2} \frac{d}{d\theta} \left(\frac{1}{r^2} \frac{dr}{d\theta} \right) - \frac{L^2}{r^3} = F(r).$$

We now make the key substitution $u = \frac{1}{r}$, so

(9) $$\frac{dr}{d\theta}= -\frac{1}{u^2} \frac{du}{d\theta} = - r^2\frac{du}{d\theta}$$

and the radial equation (8) becomes

(10) $$\frac{d^2u}{d\theta ^2} + u = - \frac{F(\frac{1}{u})}{L^2u^2}.$$

The first integral of this equation, obtained by multiplying by $\frac{du}{d\theta}$ and integrating with respect to $\theta$ (and noting that $\frac{d\varphi}{dr} = u^2\frac{d\varphi}{du}$) is

(11) $$\frac{1}{2} \frac{du}{d\theta}^2 + \frac{1}{2}u^2 + \frac{\varphi}{L^2} = \text{constant}.$$

We call this the radial energy equation. Since

(12) $$\frac{du}{d\theta} = \frac{1}{Lu^2} \frac{du}{dt} = - \frac{1}{L} \frac{dr}{dt} = -\frac{v_r}{L}$$

where $v_r$ is the radial velocity, and since the transverse velocity is $v\theta = r\dot{\theta} = \frac{L}{r} = Lu$, the constant on the right side of Equation (11) is just $\frac{E}{L^2}$, where

(13) $$E = \frac{1}{2} v^2_r + \frac{1}{2} v^2 \theta + \varphi (r)$$

is the energy per unit mass. The equations become particularly simple for the Kepler problem describing motion under gravity, since

(14) $$F(r) = - \frac{GM}{r^2} = ? GMu^2$$

and Equation (8) becomes

(15) $$\frac{d^2u}{d\theta^2} + u = \frac{GM}{L^2}$$

where the solution to Equation (15) is

(16) $$u = C \cos(\theta ? \theta_0) + \frac{GM}{L^2}$$

where $C$ and $\theta_0$ are constants. Writing

(17) $$e = \frac{CL^2}{GM}$$

(18) $$a = \frac{L^2}{GM(1 - e^2)}$$

So Equation (16) becomes

(19) $$r [1 + e \cos(\theta - \theta_0)] = a(1 ? e^2)$$

which is the standard equation of an ellipse of semi-major axis $a$ and eccentricity $e$, with the mass $M$ at one focus and the long axis oriented along the direction $\theta - \theta_0$. For $e < 1$, $r$ is bounded and the orbit is a closed ellipse. For $e \ge 1$ it is possible for $r \to \infty$ and the orbit is open a hyperbola $(e > 1)$ or a parabola $(e = 1)$. By substituting Equation (16) into the radial energy equation (11) and using equations (17) and (18), it is easily verified that the total energy (per unit mass) of the motion is

(20) $$E = -\frac{GM}{2a}$$

The shape of the ellipse is completely specified by the geometric parameters $a$ and $e$, or, equivalently, by the dynamical parameters $E$ and $L$. Note that a hyperbolic solution with $e > 1$ corresponds to $a < 0$ and $E > 0$, so the dividing line between bound and unbound motion is $E = 0$. We can use Equation (20) to eliminate $C$ from Equation (17) to find

(21) $$e^2 = 1 + \frac{2EL^2}{G^2M^2}$$

It's not possible to write down a closed-form analytic expression for $r$ as a function of time $t$, but we can derive a parametric solution to the problem, as follows. Let's assume $E < 0$. A similarderivation holds for the unbound case. Writing

(22) $$E = \frac{1}{2}v^2_r + \frac{1}{2} v^2_\theta - \frac{GM}{r} = \frac{1}{2} v^2_r + \frac{1}{2} \frac{L^2}{r^2} - \frac{GM}{r}$$

we can rearrange to find

(23) $$v^2_r = 2\left(E + \frac{GM}{r} \right) - \frac{L^2}{r^2}$$

is equal to

(24) $$- \frac{2E}{r^2} \left(- r^2 - \frac{GM}{2E} r + \frac{L^2}{2E}\right)$$

But the solutions to this quadratic in $r$ are known, since $v_r = 0$ at the two turning points of the orbit, corresponding to pericenter, $\theta - \theta_0$, $r = rp = a(1?e)$, and apocenter, $\theta = \theta_0 \pm \pi$, $r = ra = a(1 +e)$. The quadratic must therefore factorize as $(ra ? r)(r ? rp)$, and we have

(25) $$v^2_r = - \frac{2E}{r^2}(ra ? r)(r ? rp)$$.

Since $v_r = \frac{dr}{dt}$, we can say

(26) $$t = \int \frac{dr}{v_r}$$

is equal to

(27) $$\frac{1}{\sqrt{-2E}} \int \frac{r dr}{\sqrt{(ra ? r)(r ? rp)}}$$

We can do the integral by using the substitution $r = a(1 - e \cos \eta)$. Then $dr = ae \sin \eta d\eta$ and $ra - r = ae(1 +\cos \eta)$, $r ? rp = ae(1 ? \cos \eta)$ so the denominator is $ae \sin \eta$. The result is

(28) $$t = \frac{a}{\sqrt{-2E}} \int (1 ? e \cos \eta) d\eta$$

is equal to

(29) $$\frac{a}{\sqrt{-2E}} (\eta ? e \sin \eta)$$

By the defintion of $\eta$, $v_r < 0$ for $\eta < 0$ and $v_r > 0$ for $\eta > 0$, so this result is correct over the entire orbit. Writing $E = \frac{-GM}{2a}$, the leading coefficient is $\sqrt{\frac{a^3}{GM}} \equiv \Omega^{-1}$, where $\Omega$ is the mean motion. Thus we have a complete parametric solution

(30) $$\Omega t = \eta - e \sin \eta$$

(31) $$r = a(1 ? e \cos \eta)$$

Equation (30) is called Kepler’s equation. It must be solved numerically for $\eta$ given $t$. The parameter $\eta$ is called the eccentric anomaly. Once $r$ is known from Equation (31), $\theta$ can be determined from Equation (19). Finally, the period of the motion is

(32) $$T = \frac{2\pi}{\Omega} = 2\pi \left(\frac{a^3}{GM} \right)^{1/2}$$

which is Kepler’s Third Law. $\blacksquare$

Whew, that was long. With this solution, we explored a lot of topics and ideas at once so if you need a break then this is a fine place to stop for a bit. Just make sure that you remember the basic idea of the solution, and [more importnatly] the meaning of the $N$-Body Problem before we move on.

Case 2: Three Body Problem

The hree body problem holds the same concept as the two body problem, however instead of $n = 2$ we have $n = 3$. While it is true that this makes the problem even harder, our friend Euler comes up with a phenomenal solution that is simple and exact.

(Source: Frank, Juhan) Note: The numbers for denoting equations start from (7) as this soution was midway through the lecture published by Prof. Juhan Frank.

If all particles are collinear, all the vectors $s_i$, $x_i$ and $G$ are proportional to one another. Without loss of generality let’s suppose that $m_2$ lies in between the other two masses. Then $s_3$ points from $m_1$ to $m_2$, $s_1$ points in the same direction and sense as $s_3$ from $m_2$ to $m_3$, and $s_2$ points back from $m_3$ to $m_1$. Therefore, we can write

(7) $$s_1 = \lambda s_3, s_2 = -(1 + \lambda)s_3$$

where $\lambda$ is a positive scalar. Expressing everything in the equations of motion in terms of $s_3$ and $\lambda$, one obtains after some algebra a fifth degree polynomial in $\lambda$ with one single positive real root which is a function of the three masses; and, $s_3$ obeys a two-body equation of the form

(8) $$\ddot{s}_3 = - \frac{m_2 + m_3(1 + \lambda)^{-2}}{m_2 + m_3(1 + \lambda)} \frac{GMs^3}{s^3_3}$$

Thus the particles move along confocal ellipses of the same eccentricity (i.e. similar ellipses) and the same orbital period around the common center of mass, always lined up and separated by distances obeying eq. (7).

Case 3: The Planetary Problem

The planetary problem is the $n$-body problem in the case that one of the masses is much larger than all the others. A prototypical example of a planetary problem is the Sun-Jupiter-Saturn system, where the mass of the Sun is about $1000$ times bigger than the masses of Jupiter or Saturn. An approximate solution to the problem is to decompose it into $N - 1$ pairs of star-planet Kepler problems, treating interactions among the planets as perturbations. Perturbative approximation works well as long as there are no orbital resonances in the system, that is none of the ratios of unperturbed Kepler frequencies is a rational number. Resonances appear as small denominators in the expansion.

Stage 3: So Close yet So Far (Sundman and Wang)

Sundman shows that one way of solving the classical $n$ - body problem is "the $n$ - body problem by Taylor series", which is an implementation of the Power series solution of differential equations.

We start by defining the system of differential equations:

$$\frac{d^2\mathbf{x}_i(t)}{dt^2}=G \sum_{k=1,k\neq i}^n \frac{m_k \left(\mathbf{x}_k(t)-\mathbf{x}_i(t)\right)}{\left|\mathbf{x}_k(t)-\mathbf{x}_i(t)\right|^{3}}$$

As $x_i (t = t_0)$ and $\frac{dx_i(t)}{dt_{t=t_0}}$ are given as initial conditions, every $\frac{d^2x_i(t)}{dt^2}$ are known. Differentiating $\frac{d^2x_i(t)}{dt^2}$ results in $\frac{d^3x_i(t)}{dt^3}$ which at $t_0$ which is also known, and the Taylor series is constructed iteratively.

In order to generalize Sundman's result for the case $n > 3$ (or $n = 3$ and $c = 0$) one has to face two obstacles:

1.) As it has been shown by Siegel, collisions which involve more than two bodies cannot be regularized analytically, hence Sundman's regularization cannot be generalized.

2.) The structure of singularities is more complicated in this case: other types of singularities may occur:

a.) Collisions of two or more bodies, but for which $q(t)$ (the bodies' positions) remains finite. (In this mathematical sense, a "collision" means that two point-like bodies have identical positions in space.)

b.) Singularities in which a collision does not occur, but $q(t)$ does not remain finite. In this scenario, bodies diverge to infinity in a finite time, while at the same time tending towards zero separation (an imaginary collision occurs "at infinity").

Lastly, Sundman's result was generalized to the case of $n > 3$ bodies by Wang in the 1990s. Since the structure of singularities is more complicated, Wang had to leave out completely the questions of singularities. The central point of his approach is to transform, in an appropriate manner, the equations to a new system, such that the interval of existence for the solutions of this new system is $[0,\infty)$.

Concluding Remarks

Whew, that was a difficult and long endeavor. We didn't get to solve the problem at hand, but we got to understand it, and see the different attempts of solving the problem. If one was to solve this, our pathway to understanding outer space would be much clearer. I hope that this was not too difficult to understand, and that you enjoyed this great exploration. Keep exploring on expii!