n-body problem

This article is about an ancient, classical problem in classical mechanics. For the problem in quantum mechanics, see Many-body problem.

The n-body problem is an ancient, classical problem of predicting the individual motions, and forces on same, of a group of celestial objects interacting with each other gravitationally. Solving this problem -- from the time of the Greeks and on -- has been motivated by the need to understand the motions of the Sun, planets and the visible stars. In the 20th century, globular-cluster star-structures became an important n-body problem too. The n-body problem in general relativity is considerably more difficult.

The classical physical problem can be informally stated as: given only the present positions and velocities of a group of celestial bodies, otherwise their orbits, predict their forces owning to their interactions; and as a consequence, predict their motions for all past and future times.

To this purpose the two-body problem has been completely solved. For n = 3, solutions exist for special cases, and also for the so-called Restricted 3-body Problem discussed below. A general, classical solution in terms of first integrals is known to be impossible. An exact theoretical solution for arbitrary n can be approximated via Taylor series, but in practice such an infinite series must be truncated, giving at best only an approximate solution. In addition, many solutions by numerical integration exist, but these too are approximate solutions (see below).

This page develops a general (statics), closed-form solution for determining the reactive forces for the general n-body problem case, owing to a concentrated applied body force (an external force and moment). Rigid-body (i.e., where A 's = 1, see below) applications include distribution of useful loads in a finite element model (FEM); 3D rigid-body rivet analyses; and Astronomy problems. Spatial (Astrodynamic) problems and the like are particular solutions to this new general n-body problem solution; and once the initial positions and forces of the bodies, {m_ζ}, ζ = 1, 2, ... N, are known, velocities and accelerations, i.e. their motions, may be determined. Determining equations and numerical values for Astronomy problems is beyond this page's scope (see Astrodynamics). Structural analyses or soil analyses elastic solutions (A 's < 1) are possible too. Again, the latter applications are beyond the scope of this page. Note especially: rigid-body solutions for elastic structures are not valid.

Informal version of the Newton n-body problem

Some textbooks state the n-body problem problem the following way: consider $n$ point-masses $m_{1}$, ..., $m_{n}$ in three-dimensional (physical) space. Suppose the attractive force experienced between each pair of particles is Newtonian (i.e., conform to Sir Isaac Newton's (1643-1727) Laws of Motion). Then determine the spatial position of each particle at every future (or past) moment of time if initial positions and velocities in space are specified for every particle at some instant, $t_{0}$ of time.

The above declaration is misleading because Newton, by using three orbital positions of a planet's orbit—obtained from John Flamsteed —was able to produce an equation predicting the planet's motion, i.e., gave its orbital properties. Having done so they soon discovered those equations did not predict correctly some orbits and Newton realized it was because of the gravitational interactive forces amongst all the planets effecting all their orbits. Thus the need and rise of the n-body problem in the early 17th century. So, strictly speaking, the attractive forces do not conform to Newton's Laws of Motion completely; they conform to Jean Le Rond D'Alembert's non-Newtonian 1st and 2nd Principles and to the n-body problem physics (and equations) developed below (see Planet Forces and Motions Within Our Solar System below). An aside: these mathematically undefined perturbations still exist even today and the planetary orbits have to be constantly updated.

The above goes right to the heart of the matter as to what exactly the n-body problem is all about. It is not sufficient to just specify the initial position and velocity to determine a planets orbit—the gravitational interactive forces have to be known too.

Newton does not say it directly but implies in his Principia the n-body problem is unsolvable^,. Newton said (ref Cohen) in his Principia, paragraph 21: "And hence it is that the attractive force is found in both bodies. The Sun attracts Jupiter and the other planets, Jupiter attracts its satellites and similarly the satellites act on one another. And although the actions of each of a pair of planets on the other can be distinguished from each other and can be considered as two actions by which each attracts the other, yet inasmuch as they are between the same, two bodies they are not two but a simple operation between two termini. Two bodies can be drawn to each other by the contraction of rope between them. The cause of the action is twofold, namely the disposition of each of the two bodies; the action is likewise twofold, insofar as it is upon two bodies; but insofar as it is between two bodies it is single and one"...Newton concluded via his 3rd Law that "according to this Law all bodies must attract each other."

Perhaps the best way to say what the n-body problem is, is to say what was initially stated: given only the present positions and velocities of a group of celestial bodies, otherwise their orbits, predict their forces owning to their interactions; and as a consequence, predict their motions for all past and future times.

King Oscar II Prize: Historical Perspective

The problem of finding the general solution of the n-body problem was considered very important and challenging. Indeed in the late 19th century King Oscar II of Sweden, advised by Gösta Mittag-Leffler, established a prize for anyone who could find the solution to the problem. The announcement was quite specific:

Given a system of arbitrarily many mass points that attract each according to Newton's law, under the assumption that no two points ever collide, try to find a representation of the coordinates of each point as a series in a variable that is some known function of time and for all of whose values the series converges uniformly.

In case the problem could not be solved, any other important contribution to classical mechanics would then be considered to be prize-worthy.

The prize was finally awarded to Poincaré, even though he did not solve the original problem. (The first version of his contribution even contained a serious error; for details see the article by Diacu). The version finally printed contained many important ideas which led to the development of chaos theory. The problem as stated originally was finally solved by Karl Fritiof Sundman for n = 3 (see below).

Classical Mathematical formulation

Since gravity is responsible Newton reasoned for the motions of planets and stars, he expressed their gravitational interactions in terms of his Law of Universal Gravitation, but soon found, and like everyone else since, it was not a general solution to the n-body problem (see below). However, his Law was the catalysis which started the modern quest for the solution to the n-body problem.

The Law of Universal Gravitation is expressed in scalar form as:

Scalar form wrong form!!!!!!!

${\vec {F}}={\frac {Gm_{1}m_{2}}{r^{3}}}{\vec {r}}\qquad (1)$

where F, a scalar, is, by Newton's 3rd Law, the force between the two masses; m₁ and m₂ are two mass bodies attracted to each other, assumed reduced to point-masses; and r, a scalar too, is the distance between the masses. Newton proved in his Principia a spherically-symmetric body can be modeled as a point-mass.

Notice the difference between the scalar form and the vector form:

${\vec {F}}={\frac {Gm_{1}m_{2}}{r^{3}}}{\vec {r}}\qquad (1)$

where both F and r are vectors giving both magnitude and direction in some pre-defined Euclidean coordinate system.

The constant G is special, it is the gravitational constant. G, the constant of proportionality, is the parameter making Newton's Law emperialistic (i.e., capable of being verified or disproved by observation or experiment). To this purpose see Lord Cavendish's experiment measuring of G. Then what is G physically? G allows mass (matter phenomenon) when combined with distance (spatial phenomenon) to be linked to force ("force" phenomenon). All three phenomena have different units of measurement. Newton's equation is empirical because G is empirical, but (m₁m₂/r³)r is not a force and so G(m₁m₂/r³)r is not equal to F_g (gravitational force) -- but only equivalent to it (and that only by assuming we know what "force" is: in this case it is a field force). His equation establishes only an equivalence relationship via G: the equation maps the mass and spatial phenomena -to- force. As will be explain below, force is a functional.

G allows Newton's 2nd Law and all others too, to not require incorporating any constants of proportionality.

Either the scalar or vector forms can be set equal to Newton's 2nd Law, F = m a , where a is the directed acceleration. Notice the absents of a constant of proportional. Continuing with this purpose of equating and using the vector form, let F _ij = f(m_i,m_j,r_ij), where i and j are integers, 1, 2, 3,... Otherwise, F is a function of mass and geometry. Then to employ Newton's Law Universal, let F₁₂ = (Gm₁m₂/r³) r be the force on m₁ owning to mass m₂; likewise, let F ₂₁ =(Gm₁m₂/r³)r be the force on m₂ owning to m₁. This way one arrives at two separate equations, and by Newton's 3rd Law, F ₁₂ = F ₂₁ (see the Two-body Problem below). Imagine if there were three masses involved: there would be six separate equations, F₁₂, F₂₁, F₁₃, F₃₁, F₂₃, F ₃₂. If four masses (or more) were involved there would be a mathematical mess.

Newton's Law of Universal Gravitation is generally considered applying only to celestial objects, but is equally useful for Earth science (calculating Ocean tides, calculating the 13-mile high ellipsoidal budge at the equator, calculating the fuel requirements verses weight for rockets, etc.). For those applications there is the "Laws of Weights". Bodies weigh most at the surface of the Earth of course. Below the surface the weight decreases as the distance to the center decreases linearly; above the surface the weight decrease as to the square of the distance increases (Newton's Law).

Newton thought in terms of a change in momenta ∆m''v; where as Euler, for changing parameters as a function of other variables, thought in calculus terms, and momenta became a rate of change between variables, like F = m''a = m d²r/dt². This was not just a change in notation started by G. W. Leibniz (1646-1716) but a change in concept. In calculus terms Newton's notation indicated a rate of change in variables—with variables wearing hat-like dots (for, say, a change in velocity with one dot), double dots for acceleration a. Newton's momentum became Euler's Eulerian acceleration with a = d²r/dt². In what follows Newton's notation is used.

Most textbooks present the following general solution approach for the n-body problem: it means in mathematical terms finding a global solution describing the $n$-body problem of celestial mechanics, which is an initial-value problem for ordinary, second order, differential equations. Given initial values for the positions $\mathbf {q} _{j}(0)$ and velocities ${\dot {\mathbf {q} }}_{j}(0)$ of n particles (j = 1,...,n) with $\mathbf {q} _{j}(0)\neq \mathbf {q} _{k}(0)$ for all mutually distinct j and k , find the solution of the second order system:

Eq. $\qquad (2)$

$m_{j}{\ddot {\mathbf {q} }}_{j}=G\sum \limits _{k\neq j}{\frac {m_{j}m_{k}(\mathbf {q} _{k}-\mathbf {q} _{j})}{|\mathbf {q} _{k}-\mathbf {q} _{j}|^{3}}},j=1,\ldots ,n$

where $m_{1},m_{2},\ldots m_{n}$ are constants representing the masses of n point-masses, $\mathbf {q} _{1},\mathbf {q} _{2},\ldots ,\mathbf {q} _{n}$, are 3-dimensional vector functions of the time variable t, describing the positions of the point masses, and G is the constants of proportionality; i.e., gravitational constant.

This equation is Newton's second law of motion on the (left side) combined with his Law of Universal Gravitation (right side). The left-hand side is the mass times acceleration for a particular j^th particle (F = ma in vector form); whereas the right-hand side is the sum of the forces on that particle relative to another particle. The forces are assumed here to be gravitational as given by Newton's law of universal gravitation; thus, they are proportional to the masses involved, and vary as the inverse square of the distance between the masses. The power in the denominator is three instead of two to balance the vector difference in the numerator, which is used to specify the direction of the force. An aside, Euler actually transformed Newton's 2nd Law equation — of changing momentum, assumed equivalent to net force - into an ordinary, second order differential equation.

A general, closed-form, n-body problem solution (i.,e., n-formulas) can not be derived employing Newton's Law of Universal Gravitation equation as a solution approach.

see General considerations for the Classical Solution below for further discussion.

Two-body problem

Main article: Two-body problem

The two-body problem, n = 2, was completely solved by Johann Bernoulli (1667-1748) by classical theory (and not by Newton) by assuming the main point-mass was fixed, is outlined here Consider then the motion of two bodies, say Sun-Earth, with the Sun fixed, then:

m₁a₁ = (Gm₁m₂/r³₁₂)(r₂ - r₁)…..Sun-to-Earth m₂a₂ = (Gm₂m₁/r³₂₁)(r₁ - r₂)…..Earth-to-Sun

The equation describing the motion of mass m₂ relative to mass m₁ is readily obtained from the differences between these two equations and after canceling common terms gives: α + (η/r³)r = 0, where

• α is the Eulerian acceleration d²r/dt²;
• r = r₂ - r₁ is the vector position of m₂ relative to m₁;
• and η = G(m₁ + m₂).

This is the fundamental differential equation for the two-body problem Bernoulli solved in 1734. Notice for this approach forces have to be determined first, then the equation of motion resolved. This differential equation has elliptic , or parabolic or hyperbolic solutions^{, ,}

It is incorrect to think of m₁ (the Sun) as fixed in space when applying Newton's Law of Universal Gravitation, and to do so leads to erroneous results. Dr. Clarence Cleminshaw calculated the approximate position of the Solar System's true barycenter, a result achieved mainly by combining only the masses of Jupiter and the Sun. Science Program stated in reference to his work: "The Sun contains 98 per cent of the mass in the solar system, with the superior planets beyond Mars accounting for most of the rest. On the average, the center of the mass of the Sun-Jupiter system, when the two most massive objects are considered alone, lies 462,000 miles from the Sun's center, or some 30,000 miles above the solar surface! Other large planets also influence the center of mass of the solar system, however. In 1951, for example, the systems' center of mass was not far from the Sun's center because Jupiter was on the opposite side from Saturn, Uranus and Neptune. In the late 1950s, when all four of these planets were on the same side of the Sun, the system's center of mass was more than 330,000 miles form the solar surface, Dr. C. H. Cleminshaw of Griffith Observatory in Los Angeles has calculated." This means the Sun wobbles and Sun-spots are possibly caused via the movement of the barycenter, owing to Jupiter's 11-year cycles, producing Sun-spots every 22 years. It further needs to be pointed out the total mass orbiting the Sun is probably equal to the Sun's own mass.

The Real Motion v.s. Kepler's Apparent Motion

The Sun wobbles as it rotates around the galactic center, dragging the Solar System and Earth along with it. What mathematician Kepler did in arriving at his three famous equations was curve-fit the apparent motions of the planets using Tycho Brahe's data, and not curve-fitting their true circular motions about the Sun (see Figure). Both Robert Hooke and Newton were well aware Newton's Law of Universal Gravitation did not hold for the forces associated with elliptical orbits. In fact, Newton's Universal Law doesn't account for the orbit of Mercury, the Asteroid Belt's gravitational behavior, or Saturn's Rings. Newton stated (in the 11th Section of the Principia) the main reason however for failing to predict the forces for elliptical orbits was his math model was for a body confined to a situation that "hardly exist in the real world," namely, the "motions of bodies attracted toward an unmoving center." Some present physics and astronomy textbooks don't emphasize the negative significance of Newton's assumption and end up teaching that his math model is in effect reality. It is to be understood the classical two-body problem solution above is a mathematical toy. See also Kepler's first law of planetary motion.

Newton conveniently fixed the Sun so he could do simple calculations (in effect he cheated), but all following after him have also made the same mistake by analytically fixing the Sun too. They perpetuated the mathematical toy (see Truesdell's Essays in the History of Mechanics referenced below). An aside: Newtonian physics doesn't include (among other things) relative motion and may be the root of the reason Newton "fixed" the Sun. The Sun's wobbling means the n-body problem's solution realistically is much more complicated than maybe previously thought (see below).

Three-body problem

Main article: Three-body problem

In the past not much was known about the n-body problem for n equal to or greater than three. The case for n = 3 has been the most studied. Many earlier attempts to understand the Three-body problem were quantitative, aiming at finding explicit solutions for special situations.

• In 1687 Isaac Newton published in the Principia the first steps in the study of the problem of the movements of three bodies subject to their mutual gravitational attractions, but his efforts resulted in verbal descriptions and geometrical sketches; see especially Book 1, Proposition 66 and its corollaries (Newton, 1687 and 1999 (transl.), see also Tisserand, 1894).
• In 1767 Euler found collinear motions, in which three bodies of any masses move proportionately along a fixed straight line. The circular restricted three-body problem is the special case in which two of the bodies are in circular orbits (approximated by the Sun-Earth-Moon system and many others).
• In 1772 Lagrange discovered two classes of periodic solution, each for three bodies of any masses. In one class, the bodies lie on a rotating straight line. In the other class, the bodies lie at the vertices of a rotating equilateral triangle. In either case, the paths of the bodies will be conic sections. Those solutions led to the study of central configurations , for which ${\ddot {q}}=kq$ for some constant k>0 .
• A major study of the Earth-Moon-Sun system was undertaken by Charles-Eugène Delaunay, who published two volumes on the topic, each of 900 pages in length, in 1860 and 1867. Among many other accomplishments, the work already hints at chaos, and clearly demonstrates the problem of so-called "small denominators" in perturbation theory.
• In 1917 Forest Ray Moulton published his now classic, An Introduction to Celestial Mechanics (see references) with its plot of the Three body restricted Solution (see below).

Motion of three particles under gravity, demonstrating chaotic behaviour

Specific solutions to the Three-body problem result in chaotic motion with no obvious sign of a repetitious path.

The restricted problem (both circular and elliptical) was worked on extensively by many famous mathematicians and physicists, most notably by Poincaré at the end of the 19th century. Poincaré's work on the Restricted Three-body Problem was the foundation of deterministic chaos theory. In the restricted problem, there exist five equilibrium points. Three are collinear with the masses (in the rotating frame) and are unstable. The remaining two are located on the third vertex of both equilateral triangles of which the two bodies are the first and second vertices. See Figure below.

This may be easier to visualize if one considers the more massive body (e.g., Sun) to be "stationary" in space, and the less massive body (e.g., Jupiter) to orbit around it, with the equilibrium points maintaining the 60 degree-spacing ahead of, and behind, the less massive body almost in its orbit (although in reality neither of the bodies are truly stationary, as they both orbit the center of mass of the whole system—about the barycenter). For sufficiently small mass ratio of the primaries, these triangular equilibrium points are stable, such that (nearly) massless particles will orbit about these points as they orbit around the larger primary (Sun). The five equilibrium points of the circular problem are known as the Lagrange points. See figure below:

Restricted 3-Body Problem

In the Restricted 3-Body Problem math model figure above (ref. Moulton), Lagrangian points L₄ and L₅ are where the Trojan planetoids resided; m₁ is the Sun and m₂ is Jupiter. L₂ is where the asteroid belt is. It has to be realized for this model, this whole Sun-Jupiter diagram is rotating about its barycenter. The h-circles and closed loops echo the electromagnetic fluxes issued from the Sun and Jupiter. (An aside, see Meirovitch, pages 414 and 413 for the Three-body Problem.)

The Restricted Three-body Problem assumes the mass of one of the bodies is negligible; For a discussion of the case where the negligible body is a satellite of the body of lesser mass, see Hill sphere; for binary systems, see Roche lobe; for another stable system, see Lagrangian point.

Practical n-body Problem Theoretical Solution

This solution occurred while trying to distribute useful applied loads (i.e., human body loads, luggage, attached extra hardware, etc.) to the grids (spatial points) belonging to surrounding basic structures of finite element models. How else would one distribute applied loads if not by a n-body problem type solution? The similarity for determining the reactive loads for a 3D-rivit analyses also became apparent.

Turns out any solution to the n-body problem requires using Jean Le Rond D'Alembert's 1st and 2nd Principles, is key.

Hookean Switches and Springs, Point-Masses and the Inertia Tensor

Index notation is used in this section. Indexes and scalars are plain italic text; vectors are in bold, italic text; Cartesian tensors are in plain italic text.

It is convenient to define a vector operator, a functional, to be employed in the solution of the n-body Problem.

Hookean Switches: Let there be a vector function, 0 ≤ A_ζi ≤ 1, where i := x, y, z or 1, 2, 3, are the component directions, ζ := 1,...N, is the point-masses' m_ζ coordinate number. The A_ζi's are functionally the feasible, translational reaction component's allowable degrees of freedom (DOF) at P_ζ⇒ (x_ζ, y_ζ, z_ζ), specified for any and all point-masses m_ζ, where the N point-masses m_ζ are defined in a right-handed, Euclidean coordinate system in File:Reals R.jpg³ with (x, y, z) coordinates. Coupling forces are reacted at these point-masses too and are subject to the same feasible translational constraints.

Applying Hookean switches A 's at the supports or reaction points for every feasible translational reaction, i.e., for each P_ζ, would be the same as showing those feasible or allowable translational vectors visually on a free-body diagram; but instead of showing by default all three DOF possible translational vectors, say ( P _ζ )_i, with i = 1,2,3. There must be at least three point-masses for stability, having a total of six DOF , meaning there must be at least six feasible reactions for these three points; i.e., six A_iζ's. (This latter requirement is the same as balancing an aircraft's FEM on three points and then applying Jean Le Rond D'Alembert's two Principles, resulting in no load on the three points.)

File:BOSCOVITCH Theorem.jpg

Radius Vectors and Massless Points

(Definition continues.) Let the radius vector be:

Eq. $\qquad (3)$

File:Radius Vector r.jpg

where i, j, k with the "hats" are the unit directional vectors, and the A 's are the Hookean switches or springs defined above. For example, for m_ζ at P_ζ ⇒ (x_ζ, y_ζ, z_ζ): if all A_ζi = 0 then r _ζ = 0; if all A_ζi = 1 then there are three reaction components at ζ. These conditions describe Hookean switches. All other values of A_ζi ≤ 1, but not zero, describe Hookean springs. For these conditions, if the reaction points' material is structural; i.e., elastic ( σ = E ε ), then it will be further assumed the structure will conform to a similar linear law, say Hooke's Law, ( F = k x ) of load-deformation proportionality, where k is the spring rate, and wherein the A 's are analogous to k (this assumption is not entirely necessary as F = f(x) can be any functional without its exact definition given).

(Definition continues.) Then for each point P_ζ there may be 0, or up to 3 DOF, meaning m_ζ may be, depending on A_ζi's components for any particular DOF in that direction, elastic or rigid body. If r _ζ is defined with Hookean Switches then r _ζ are said to be a Hookean vector.

The concept is: Hookean switches allow point-masses to respond, or not respond, as reaction forces in selected directions, by setting selected A_iζ's components to 0; and vice-verse. When all the feasible A_ζi = 1 the solution is a rigid body solution. In this way Hookean switches allows for the solution of the n-body Problem. Hookean springs allows a type of finite element model (FEM) to be constructed (i.e., 3D rivet analyses).

Let M := m₁ + m₂ + m₃ + ...+ m_N, for ζ = 1, 2, ..., N point-masses. Of course for each m_ζ there is a spatial point P_ζ ⇒ (x_ζ, y_ζ, z_ζ).

Let r _ζi at P_ζ, where ζ = 1,...N mass points, be a Hookean vector for each m_ζ in all (x, y, z) space, for a right-hand, rectangular Cartesian coordinate system. The A 's here are the components of A. Then the six massless barycenters (six naturally develop when using index notation) for {M} ≤ 3 are, by simple statics:

File:Centroids.jpg

Eq.s $\qquad (4)$

File:Centroid Set.jpg

If the reaction or support pattern is semi-symmetrical then the centroid set will be less than six, as, say when for example, x_ij = x_ji; if the pattern is entirely symmetrical then these six centroids will revert back to the usual three. Continuing, the massless centroidal mass moments of inertia for N masses, m_ς, again by simple statics are:

Eq.s $\qquad (5)$

File:Inertia Tensor.jpg

where, say x bar_xy, etc, are the barycenters calculated above.

But before that...qualification for a right-hand coordinate system arises because some USA Aerospace firms use the left-hand system, making their signs for the cross products positive. If the inertia tensor if symmetric and positive-definite (via the value of the its determinant), then the coordinate system is right-handed ― always check this, and if not then your signs are wrong.

Applied forces and moments may be Hookean too. These Hookean definitions (and concept) are missing from physics textbooks, but are somewhat similar in concept to grid-point releases in a finite elements analyses. Hookean points are not Lagrangian coordinates.

The Force

Interestingly, if empirical force laws and Euclidean geometry are true, and the Principle of Equilibrium is true of course, then any logical combination therein is also true. In all that follows is naught but simple physics (statics).

If, for example, when the vector s_i = Є _ijk r_j θ _k (analogous to s = rθ) from simple geometry is combined with D'Alembert's two Principles R = 0 and M = 0 , then a formula called The Force, can be developed giving the reactive forces for a pattern of points subject to an external resultant force F _i and moment M_0i. r_j may be a Hookean vector, and s_i covers all space. It will be shown The Force can be used to completely solve any n-body Problem for any n ≥ 3.

Let the radius vectors to the locations of the ζ = 1, 2, ...N reaction points or supports points be the set r ⇒{r_ζi)|i = 1, 2, 3} for m_ζ ≥ 3 mass-points, react to an applied resultant force F _i and moment ( M ₀)_i, i = 1, 2, 3 components, at File:C g s.jpg. Then reaction displacements (∆s_ζ)_i of any point (r _ζ)_i owing to angular rotations (∆θ_ζ)_i are:

Eq. $\qquad (6)$

Reaction Displacements: (∆s_ζ)_i= Є_ijkr_ζj∆θ_ζkA_ζk,

where File:Delta S 01.jpg is a vector. Assume only infinitesimal angular rotations (∆θ_ζ)_i are associated with reaction point motions File:Delta S 01.jpg in any instant of time. Further, since this is an instant of time (time is fixed) the solution becomes one of statics (i.e., quasi-steady). Remember, ζ := 1,...N, is the point-mass m_ζ coordinate number, where m ≥ 3.

The load center is defined to be at File:C g s.jpg, who coordinates are basic input. File:Centroids 01.jpg are use later to calculate the sum of the moments — are equal to the {P_ζ} coordinate set for ζ = 1, ... N reaction points, and are the centers of gravity of those point-masses. Real masses may be reduced to point-masses the same way Newton did it. In this analysis physical real mass is not a variable; only forces, reactions, moments, and geometry are real variables.

The Hookean File:Delta S 01.jpg equations in expanded matrix form are exactly analogous to M = r x F or bisor (second order tensor) M_i = Є_ijkr_jF_k.

File:Delta S 03.jpg

Assume infinitesimal moments (torques) M _ζi are used in lieu of ∆θ_ζi. Expressed mathematically:

File:M 01a.jpg

P _ζi are the three force components at the support or reaction points. File:Rbar 01.jpg_ζi are the distances from the support patterns centroid to the support points, which is a detail explained more fully below. This angular rotations -to- moments assumption will be incorporated below. Interestingly, this same assumption is employed in the displacement method in finite element analysis, but there in a different context. (Other assumption relationships may be employed too.) The left hand side expanded in matrix form is:

File:Moments 02a.jpg

The [f(File:Rbar 01.jpg)]_3,3 matrix is actually the inertia tensor of the reaction or support pattern, a mathematical detail given below. What follows next is to determine the set of distances File:Rbar 01.jpg_ςi from the centroids of the reaction or support pattern to the individual mass-points. To that purpose let:

File:Rbar Equation.jpg

The radius vector r_ςi are rays to point-masses m_ς at B; the radius vector r _{o i} connects to the tail of the position vector File:Rbar 01.jpg_ζi at A; and the position vector in turn connects to the point-masses at B; and closes the loop (see figure), where:

• File:Rbar 01.jpg_ζi, i = 1, 2, 3, are the position vectors from the six centroids of the reaction pattern -to- the reaction points or supports points (or planets, etc.). The three centroids really consist of the six centroids as developed in the Section directly above: File:Rbar Converted.jpg.
• r_0i are the radius vectors (rays from the origin to...) to the reaction or support pattern's mass-points at {P_ζi}.
• r_ζi are the distances from the origin to the reaction pattern's point-masses.

These equations in expanded form are:

File:Rbar Equation Expanded.jpg

Also note, m_ς can be the centroids of reduced point-masses, but usually are just a m_ζ 's coordinate.

Suffice it to say when the reactions P _ςi are multiplied by File:Rbar 01.jpg_ςi, theFile:Rbar 01.jpg_ςi takes the form of a skewed Hermitian matrix:

Eq. $\qquad (7)$

File:Hermitian.jpg

And that product times the inverse of the inertia tensor [I_ς]⁻¹, is equal to the internal moment M_ςi:; or re-written in terms of reactions:

File:Internal Loads.jpg

P _ζi are the individual coupling reactions from the reduced moments, plus the additional translational reactions distributed about the reaction points m_ζi owing to an externally applied force F _ζi and M_oi. The resultant internal moments M _ζi as coupling reactions at each point are equal to the sum of resultant external moments at the centroids of the reaction pattern, by equilibrium considerations; i.e., by D'Alembert's Second Principle (and not by Newton's Third Law, which would be incorrect reasoning).

File:Moment Transfer.jpg

Transferring a Force Produces a Moment

(Explanatory note.) Some explanation is in order before preceding: First, referring to the cantilever beam with the off-set load, if the load F₁ is moved to the cantilever's tip's end along the shaft's axis, then (to calculate loads at the root) if the original force, F₁, is reduplicated and transferred to the tip as F ₂, then it has to be balanced by an up-load, F₃, for the two to remain in equilibrium. This up-load F₃, combined with the original down-load F₁, is a couple and equivalent to a moment M (torque). Second, at the root of the cantilever the reaction forces at the attachment points (not shown) not only react with the vertical force F₂ in pure upward vertical shear divided by the number of attachment points, but by horizontal reaction couples in- and out- of the wall too. The new couple (moment M ) is reacted by only reaction couples at the root parallel to the wall. In this way (symbolically) the external moment M_external = M _internal = rF_reaction; i.e., the external moment(s) equal the in-plane and out-of-plane static reaction couples at the attachments in the wall (D'Alembert's Second Principle). Looking at the big picture, the total external forces and moments are equal to the sum of all the internal coupling reactions (moments) plus the direct shear reactions (translational reactions) divided by the number of attachments. D'Alembert's two Principles combined.

(Continuing.) Let the reactions P_ςi above owe to a known external applied force F_i and known applied external moment M_oi. The F_i when moved from its initial position to the centroids of the reaction or support pattern, creates moments and these moments are of the form [x_L]_3,3 F _i. The [x_L]_3,3 is developed now:

(rbar _p)_i = r _pi - r bar _i,

and r_pi ⇒File:C g s.jpg is the position of the applied force and moment. Then in expanded scalar form:

Eq. $\qquad (8)$

File:Hermitian Load.jpg

[x_L]_3,3 is a skewed Hermitian.

Putting all the above together results in an equation call The Force. The Force formula gives a set of balanced reactions at the reaction points or support points (Lb.). Its proof is by calculating equilibrium again, once those reactions have been found (see Equilibrium of Forces and Moments section below).

Eq. $\qquad (9)$

( P _ζ)_i = [{F_L}_i/ΣA_i + [x_s]_i[ I _s]⁻¹[ [x_L]_i{F _L}_i + {M ₀}_i]]{A_ζ}_i.

where i = x, y, z; ζ = 1, 2, ..., N.
This solution takes an external force and moment and moves it over to the centroids of the reactions or supporting pattern (system); and by incorporating the massless inertial properties of that pattern, beams the applied loads to those reaction points (or support points). Although somewhat messy algebraically, it's that simple. It's analogous to a 3D rigid body, rivet analysis. The above in expanded matrix form, except for the inverse of the massless inertia tensor [ I _ς]⁻¹, for the three orthogonal reaction components at a ζ-reaction point is:

Eq. $\qquad (10)$

File:The Force Eq.jpg

In the expanded version (except for [ I _ς]⁻¹ terms) one can clearly see the right side of the equation is composed of two parts: the translational forces and the coupling forces.

In the two above equations:

• ζ = 1, 2, 3, ... N point-masses m_ζ coordinate number; brackets [ ] and { } represent matrices and -1 is for inverse; and subscript L stands for Load position, s stands for supports;
• ( P _ζ )_i are the vector reactions at each support point P_i⇒ (x_i, y_i, z_i);
• ( F _L)_i are the three components of the applied force (lb.), transposed; ( M ₀)_i are the the three components of the applied moment (lb.-in.), transposed; F and M _o are basic inputs;
• The load centers at File:C g s.jpg is basic input (in.).
• [ I _ς]⁻¹ is the inverse of the massless inertia tensor of the reaction pattern (via Crout's Method), is a symmetrical 3x3 matrix (in.⁻²);
• [x_L]_i is the position of the applied force and moment, 3x3 matrix; [x_L]_i( F _L) are the resulting couples in the support pattern caused by moving the applied force over to the pattern's centroids (in.).
• [x_s]_i are the centroids of the reaction pattern, 3x3 matrix (in.);
• [x_s]_i[ I _ς]⁻¹[[x_L]_i{ F _L}_i^T + { M _o}_i^T]_3,3] are the unbalanced moments reacted as couples, 3x3 matrix (lb.);
• {( F _L)_i/(Σ(A_ζ)_i)}^T_1,3 are the symmetrical (translational) reactions (lb.).

The [x_s]_i[I_ς]⁻¹[x_L] term is a congruent transformation (Meirovitch) and for all physical entities having like congruent transformations, like stress, strain, they are always positive definite if the right-hand rule is used: a fundamental property of all matter.
{A_ζ}_i are the Hookean switches defined above. If a particular motion is planar then out-of-plane reaction forces (P _ζ)_i will be zero by default. ( P _ζ)_i has been computerized, see below.

Equilibrium of Forces and Moments

The sum of the forces for the entire support or reaction point set, {P_ζ} set, for ζ = 1, ... N reaction points, are:

Eq. $\qquad (11)$

File:Sum of Forces.jpg

The applied load center is at File:C g s.jpg and is basic input. The c.g.'s (centers of gravity) of the reaction points (the m_ζ's reduced to a support-point set of point-masses) are at File:Centroids 01.jpg, i = 1, 2, ...,N; are in most cases the same set as the {P_ζ} set for ζ = 1, ... N reaction points. Then the sum of the moments for the entire supports or reactions are:

Eq. $\qquad (12)$

Equilibrium Moments of Support Pattern

Equilibrium Moments of Support Pattern

Numerical Example

Example problem: a large video game monitor weighing 179 pounds is mounted into a cantilevered (from the floor) flying (airborne) console, is supported internally by the console by 20 support (attachment) points. The external applied force, has a load c.g. of File:C g s.jpg⇒(3.96442, 13.5625, 44.89365) inches, is close to these support points. The problem is to determine three rigid-body, unit-load sets for the entire support pattern. In this example, only (as an example) a unit load of 1.0 Lb. in the vertical direction is applied to determine its unit reactions. (Two more unit load runs would ordinarily be needed to determine sideways and forward unit reactions too. Otherwise, three sets of unit-load reactions are needed to obtain a complete loading profile for, say, a FEM.) Support coordinates in inches are given in Tables 1 and 2; from which the six support c.g.'s as well as the inertia components are calculated (only once) in the Tables below.

File:N-Body Table 01.jpg

Table 1, Calculations

File:N-Body Table 02.jpg

Table 2, Calculations

Extracting values from these tables gives:

File:N-Body Table 03.jpg

And the inertia tensor is:

small

small

And in matrix form:

File:Inertia 01.jpg

File:Inertia Inverse.jpg

The load transfer matrix is:

And the support transfer matrix which transfer the applied moment as couples to the support points, is:

File:Transfer Xs.jpg

All internal reactions at the N-support points are given by the following three equations:

File:Reaction Forces.jpg

Reaction moments equations:

File:Reaction Moments Eq.jpg

Planet Forces and Motions Within Our Solar System

If we define the Solar System's domains as having expanded as far as the gravitational fields of the Sun, Jupiter, and the other outer-planets, then the Solar System, beyond the planets Neptune and Pluto, extends perhaps a thousand times further than we have yet explored. But out there we know there are billions of small objects in that zone known as Oort Cloud. If all the known mass except the Sun — that is, the planets, asteroids, moons, Trojan planetoids, etc. — is subtracted from the Sun's known mass, then it is conjectured: the remaining mass of the Solar System would be equal to all that mass in the zone, such that, that mass plus the known body masses equals the Sun's mass. This is not only reasonable gravitationally, but application of the general n-body Problem indicates it is true too, because without this additional mass — the planets, etc., has too little mass to absorb the Sun's gravitational forces — result in higher interactive forces via the n-body solution calculations, than those via Newton's Law of Universal Gravitation. These higher forces are reduced to Newtonian forces by assuming the zone 's mass is almost equal to the Sun's.

Calculating Solar Forces: suppose high above the Solar plane a picture was taken by a satellite, instantly freezing all motion of all bodies in the Solar System. This would result in a quasi-steady (static) picture of fixed positions, zero velocities and accelerations. The positions become variables (coordinates) for input into the general solution of the n-body Problem; and the main masses are already known. With this knowledge, a suggested Solar System general n-body type solution may be effected by:

1. Calculate all forces of known bodies via Newton's Law of Universal Gravitation, whose values are to be used one at a time discretely as applied external forces; and as a comparisons of the n-body solution forces -to- Newtonian forces.
2. At the inner diameter of Oort's Cloud establish, say, a hundred (or more) mass-points uniformly distributed around the inside of its sphere, whose total mass sum is equal to the Sun's mass, less the known masses.
3. Use the quasi-steady positions of the known masses combined with the mass-points of Oort's Cloud mass-points, as input to the n-body program.
4. Run the n-body computer program first using the Sun as the known force; and the planets, etc., as unknown reactions; and record the position of the Solar System's barycenter.
5. Repeat a similar run for each known mass-point, with one mass-point selected as the known external applied force; (use Newton's Law to determine the applied force).
6. Repeat computer runs time-wise incrementally until the orbit of the Solar System's barycenter is known (shades of Cleminshaw).

If these forces are known then the eccentricities of the planet's orbits can be determined. This task is way to large for it to be presented here.

n-body particle simulations can be effected via computer. This program demonstrates body-to-body interactions, i.e., almost head-on collisions, almost always results in bodies orbiting bodies rather than collisions.

More General considerations for the Classical Solution

For every solution of the problem, not only applying an isometry or a time shift but also a reversal of time (unlike in the case of friction) gives a solution as well.

In the physical literature about the $n$-body problem ($n$ ≥ 3), sometimes reference is made to the impossibility of solving the $n$-body problem (via employing the above approach). However, care must be taken when discussing the 'impossibility' of a solution, as this refers only to the method of first integrals (compare the theorems by Abel and Galois about the impossibility of solving algebraic equations of degree five or higher by means of formulas only involving roots).

Employing Newton's generalized Law of Universal Gravitation formula above, the $n$-body problem contains 6$n$ variables, since each point particle is represented by three space (displacement) and three momentum components. First integrals (for ordinary differential equations) are functions that remain constant along any given solution of the system, the constant depending on the solution. In other words, integrals provide relations between the variables of the system, so each scalar integral would normally allow the reduction of the system's dimension by one unit. Of course, this reduction can take place only if the integral is an algebraic function not very complicated with respect to its variables. If the integral is transcendent the reduction cannot be performed.

The $n$-body problem has 10 independent algebraic integrals:

• three for the center of mass;
• three for the linear momentum;
• three for the angular momentum;
• one for the energy.

This allows the reduction of variables to 6$n$ − 10. The question at that time was whether there exist other integrals besides these 10. The answer was given in the negative in 1887 by Heinrich Bruns:

Theorem (First integrals of the $n$-body problem) The only linearly independent integrals of the $n$-body problem, which are algebraic with respect to $q$, $p$ and $t$ are the 10 described above.

(This theorem was later generalized by Poincaré). These results however do not imply there does not exist a general solution to the $n$-body problem, or that the perturbation series (Lindstedt series) diverges. Indeed Sundman provided such a solution by means of convergent series. (See Sundman's theorem for the 3-body problem).

(Classical) Power Series Solution

This is the most elemental way classically of solving the n-body Problem. The theoretical expression is often called "The n-body problem by Taylor series", which is an implementation of the Power series solution of differential equations.

We start by defining the differential equations system:

${\underset {i}{\overset {2}{x_{j}}}}(t)=G\sum _{k=0,k\neq i}^{n}{\frac {{\underset {k}{m}}\left({\underset {k}{\overset {0}{x_{j}}}}(t)-{\underset {i}{\overset {0}{x_{j}}}}(t)\right)}{\left(\left({\underset {k}{\overset {0}{x_{1}}}}(t)-{\underset {i}{\overset {0}{x_{1}}}}(t)\right){}^{2}+\left({\underset {k}{\overset {0}{x_{2}}}}(t)-{\underset {i}{\overset {0}{x_{2}}}}(t)\right){}^{2}+\left({\underset {k}{\overset {0}{x_{3}}}}(t)-{\underset {i}{\overset {0}{x_{3}}}}(t)\right){}^{2}\right){}^{3/2}}}$,

where (in ${\underset {i}{\overset {2}{x_{j}}}}(t)$) the upper index $2$ indicates the second derivative with respect to time $t$, $i$ represents the number of each body and $j$ the coordinate.

Because ${\underset {i}{\overset {0}{x_{j}}}}(t_{0})$ and ${\underset {i}{\overset {1}{x_{j}}}}(t_{0})$ are given as initial conditions, then every ${\underset {i}{\overset {2}{x_{j}}}}(t_{0})$ are known. Doing implicit derivation over every ${\underset {i}{\overset {2}{x_{j}}}}(t)$ results in ${\underset {i}{\overset {3}{x_{j}}}}(t)$ which at $t_{0}$ are known because each depends on known ${\underset {i}{\overset {k<3}{x_{j}}}}(t_{0})$ precalculated and given constants and then the Taylor series are constructed theoretically in such way, performing this process infinitely.

Numerical integration

Main article: N-body simulation

n-body problems can be solved by numerically integrating the differential equations of motion. Many different ways to do this to varying degrees of accuracy and speed exist.

The simplest integrator is the Euler method, but this is only first order. A second order method is leapfrog integration, but higher-order integration methods such as the Runge–Kutta methods can be employed. Symplectic integrators are often used for n-body problems.

Numerical integration has a time complexity of O(n²), but tree structured methods, such as Barnes-Hut simulation, can improve this to O(n log n), or even to O(n) such as with the fast multipole method.

Sundman's Theorem for the 3-body Problem

In 1912, the Finnish mathematician Karl Fritiof Sundman proved there existed a series solution in powers of $t^{1/3}$ for the 3-body Problem. This series is convergent for all real t, except initial data which correspond to zero angular momentum. However these initial data are not generic since they have Lebesgue measure zero.

An important issue in proving this result is the fact that the radius of convergence for this series is determined by the distance to the nearest singularity. Therefore it is necessary to study the possible singularities of the 3-body Problems. As it will be briefly discussed below, the only singularities in the 3-body Problem are:

1. binary collisions,
2. triple collisions.

Now collisions, whether binary or triple (in fact of arbitrary order), are somehow improbable since it has been shown they correspond to a set of initial data of measure zero. However there is no criterion known to be put on the initial state in order to avoid collisions for the corresponding solution. So Sundman's strategy consisted of the following steps:

1. He first was able, using an appropriate change of variables, to continue analytically the solution beyond the binary collision, in a process known as regularization.
2. He then proved triple collisions only occur when the angular momentum L vanishes. By restricting the initial data to $\mathbf {L} \neq 0$ he removed all real singularities from the transformed equations for the 3-body problem.
3. The next step consisted in showing that if $\mathbf {L} \neq 0$, then not only can there be no triple collision, but the system is strictly bounded away from a triple collision. This implies, by using Cauchy's existence theorem for differential equations, there are no complex singularities in a strip (depending on the value of L) in the complex plane centered around the real axis (shades of Kovalevskaya).
4. The last step is then to find a conformal transformation which maps this strip into the unit disc. For example if $s=t^{1/3}$ (the new variable after the regularization) and if $|\mathop {\text{Im}} \,s|\leq \beta$ then this map is given by:

$\sigma ={\frac {e^{\pi s/(2\beta )}-1}{e^{\pi s/(2\beta )}+1}}.$

This finishes the proof of Sundman's Theorem. Unfortunately the corresponding convergent series converges very slowly. That is, getting the value to any useful precision requires so many terms, that his solution is of little practical use.

A Generalized Sundman Global Solution

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 (see below).

Lastly, Sundman's result was generalized to the case of n > 3 bodies by Q. 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 )$.

Singularities of the n-body problem

There can be two types of singularities of the n-body problem:

• 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.)
• 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").

The latter ones are called Painlevé's conjecture (no-collisions singularities). Their existence has been conjectured for n > 3 by Painlevé (see Painlevé's conjecture). Examples of this behavior have been constructed by Xia and Gerver.

See also

• Einstein–Infeld–Hoffmann equations
• Few-body systems
• Natural units
• n-body choreography
• Virial theorem

References

• Heggie, Douglas and Hut, Piet: The Gravitational Million-Body Problem, A Multidisciplinary Approach to Star Cluster Dynamics, Cambridge University Press, 357 pages, 2003.
• Leimanis, E., and Minorsky, N.: Dynamics and Nonlinear Mechanics, Part I: Some Recent Advances in the Dynamics of Rigid Bodies and Celestial Mechanics (Leimanis), Part II: The Theory of Oscillations (Minorsky), John Wiley & Sons, Inc., 1958.
• Moulton, Forest Ray: An Introduction to Celestial Mechanics, Dover, 1970.
• Meirovitch, Leonard: Methods of Analytical Dynamics, McGraw-Hill Book Co., 1970.
• Brouwer, Dirk and Clemence, Gerald M.: Methods of Celestial Mechanics, Academic Press, 1961.
• Cohen, Bernard I.: "Newton's Discovery of Gravity," Scientific American, pp. 167–179, Vol. 244, No. 3, Mar. 1980.
• Cohen, Bernard I.: The Birth of a New Physics, Revised and Updated, W.W. Norton & Co., 1985.
• Science Program's “ The Nature of the Universe ,” booklet, published by Nelson Doubleday, Inc., in 1968:
• Bate, Roger R.; Mueller, Donald D.; and White, Jerry: Fundamentals of Astrodynamics, Dover, 1971.
• Batin, Richard H.: An Introduction to The Mathematics and Methods of Astrodynamics, AIAA, 1987.
• Gallian, Dave A. and Wilson, Henry E.: "The Integration of NASTRAN Into Helicopter Airframe Design/Analysis," American Helicopter Society Pub., Reprint No. 780, May 1973. This is a paper, not a book.
• Lindsay, Robert Bruce: Physical Mechanics, 3rd Ed., D. Van Nostrand Co., Inc., 1961.
• Gelman, Harry: Part I: The second orthogonality conditions in the theory of proper and improper rotations: Derivation of the conditions and of their main consequences, J. Res. NBS 72B (Math. Sci.)No. 3, 1968. Part II: The intrinsic vector; Part III: The Conjugacy Theorem,J. Res. NBS 72B (Math. Sci.) No. 2, 1969. A Note on the time dependence of the effective axis and angle of a rotation, J. Res. NBS 72B (Math. Sci.)No. 3&4, Oct. 1971. These papers are on the Internet.
• Meriam, J. L.: Engineering Mechanics, Volume 1 Statics, Volume 2 Dynamics, John Wiley & Sons, 1978.
• Quadling, Henley: "Gravitational N-Body Simulation: 16 bit DOS version," June 1994. nbody*.zip is available at the http://www.ftp.cica.indiana.edu: see external links.
• Korenev, G. V.: The Mechanics of Guided Bodies, CRC Press, 1967.
• Eisele, John A. and Mason, Robert M.: Applied Matrix and Tensor Analysis, John Wiley & Sons, 1970.
• Murray, Carl D. and Dermott, Stanley F.: Solar System Dynamics, Cambridge University Press, 606 pages, 2000.
• Bronowski, Jacob and Mazlish, Bruce: The Western Intellectual Tradition, from Leonardo to Hegel, Dorsey Press. 1986.
• Truesdell, Clifford: Essays in the History of Mechanics, Springer-Verlay, 1968.
• Hufbauer, Dr. Karl C. (History of Science): Exploring the Sun, Solar Science since Galileo, The Johns Hopkins University Press, 1991. This book was sponsored by the NASA History Office.
• Diacu, F.: The solution of the n-body problem, The Mathematical Intelligencer,1996,18,p. 66–70
• Mittag-Leffler, G.: The n-body problem (Price Announcement), Acta Matematica, 1885/1886,7
• Saari, D.: A visit to the Newtonian n-body problem via Elementary Complex Variables, American Mathematical Monthly, 1990, 89, 105–119
• Newton, I.: Philosophiae Naturalis Principia Mathematica, London, 1687: also English translation of 3rd (1726) edition by I. Bernard Cohen and Anne Whitman (Berkeley, CA, 1999).
• Wang, Qiudong (1991). "The global solution of the n-body problem". Celestial Mechanics and Dynamical Astronomy. 50 (1): 73–88. Bibcode:1991CeMDA..50...73W. ISSN 0923-2958. {{cite journal}}: |access-date= requires |url= (help)
• Sundman, K. F.: Memoire sur le probleme de trois corps, Acta Mathematica 36 (1912): 105–179.
• Tisserand, F-F.: Mecanique Celeste, tome III (Paris, 1894), ch.III, at p. 27.
• Hagihara, Y: Celestial Mechanics. (Vol I and Vol II pt 1 and Vol II pt 2.) MIT Press, 1970.
• Boccaletti, D. and Pucacco, G.: Theory of Orbits (two volumes). Springer-Verlag, 1998.
• Havel, Karel. N-Body Gravitational Problem: Unrestricted Solution (ISBN 978-09689120-5-8). Brampton: Grevyt Press, 2008. http://www.grevytpress.com
• Saari, D. G.; Hulkower, N. D. (1981). "On the Manifolds of Total Collapse Orbits and of Completely Parabolic Orbits for the n-Body Problem". Journal of Differential Equations. 41 (1): 27–43. doi:10.1016/0022-0396(81)90051-6.

External links

• Three-Body Problem at Scholarpedia
• More detailed information on the three-body problem
• Regular Keplerian motions in classical many-body systems
• Applet demonstrating chaos in restricted three-body problem
• Applets demonstrating many different three-body motions
• On the integration of the n-body equations
• Java applet simulating Solar System
• Java applet simulating a ring of bodies orbiting a large central mass
• Java applet simulating dust in the Solar System
• Java applet simulating a stable solution to the equi-mass 3-body problem
• Java applet simulating choreographies and other interesting n-body solutions
• A java applet to simulate the 3-d movement of set of particles under gravitational interaction
• Javascript Simulation of our Solar System
• The Lagrange Points - with links to the original papers of Euler and Lagrange, and to translations, with discussion
• http://ftp.math.utah.edu/pub/tex/bib/toc/sciam1980.html#242%281%29:January:1980
• http://www.ftp.cica.indiana.edu