Three-Body Problem Problem ( Astronomers )

[choppu]

Three-Body Problem Problem ("Astronomers")

Hi,


I am having quite a few problems with the three-body problem. First of all "my" Three-body problem consists of

a static sun (mass 200000 * 10²⁴kg) in the middle of the programme, and two non static planets that revolve around the sun.


I believe this version of the problem is called the "Astronomers" version.


I am not able to find out the differential equation for my version, and would be glad if somebody could point it out to me (preferably with explanations...).


Thanks

choppu

 

[ideasrule]

Writing out the differential equations is easy: just write out Newton's second law for all three objects in vector form. Solving them analytically is impossible, as discussed in Wikipedia's article on the n-body problem:

http://en.wikipedia.org/wiki/N-body_problem

 

[Chronos]

The n-body problem is incredibly complicated. Even n=3 has not been completely solved. Fortunately, good approximations have be derived. Even large n's can be reasonably approximated with high end computers.

 

[D H]

The n-body problem is incredibly complicated. Even n=3 has not been completely solved.

Not quite true. The three-body problem can be deemed fully solved in the sense that an infinite series representation does exist. Karl Sundman developed this series in the early 1900s. He was awarded with the Pontécoulant prize for this work. However (http://scienceworld.wolfram.com/physics/RestrictedThree-BodyProblem.html)

Sundman found a uniformly convergent infinite series involving known functions that "solves" the restricted three-body problem in the whole plane (once singularities are removed through the process of regularization). Since such global regularizations are available for this problem, the restricted problem of three bodies can be considered to be complete "solved." However, this "solution" does not address issues of stability, allowed regions of motion, and so on, and so is of limited practical utility (Szebehely 1967, p. 42). Furthermore, an unreasonably large number of terms (of order 10^8,000,000) of Sundman's series are required into attain anything like the accuracy required for astronomical observations.​

Fortunately, good approximations have be derived. Even large n's can be reasonably approximated with high end computers.

It doesn't take that high-end a computer. The typical desktop computer or laptop has plenty of oomph. Think of it this way: A 21st century programmable calculator has a lot more computing power than the computers of the 1960s that were used to design trajectories for vehicles that took people to the Moon and took probes to other planets.

Where a lot of computing power is needed is in grinding through the myriad of observations of the planets made over the course of hundreds of years and eclipses made over the course of thousands of years to determine a best-fit state for all the planets at some epoch time. This is a boundary value problem from Hades. Numerically propagating a given state over time is far a simpler problem.

 

[twofish-quant]

Also there are restricted three body problems whose solutions are very well known. One common problem for which simple solutions do exist is the problem of two bodies and a spaceship.

 

[D H]

Only if the spaceship is at one of the libration points. The general problem of a spaceship subject to the gravitational force of two bodies is the restricted three body problem. The restricted three body problem, like the unrestricted three body body, is insoluble in terms of elementary functions.

 

[choppu]

Hi guys,

Writing out the differential equations is easy: just write out Newton's second law for all three objects in vector form. Solving them analytically is impossible, as discussed in Wikipedia's article on the n-body problem:

I do not quite understand how I can write out F=ma for all three objects in vector form. Can you maybe demonstrate it to me with an example please?

 

[Matterwave]

F=ma and F=-GMm/r^2 is needed (in vector form).

Now, for one object, you have

$$F_1=m_1a_1=\frac{-Gm_1m_2}{r_{12}^2}-\frac{Gm_1m_3}{r_{13}^2}$$
(Put the vectors in there, I'm too lazy atm)

So, that one object is affected by the gravitational forces on it due to the other two objects.

You get 3 such equations.

 

[choppu]

Ok, thanks for explaining.



I still do not understand how I can split up these equations to gain two (six?) first order differential equations.


Can anybody point it out to me please?

 

[D H]

One common way to express Newton's law of gravitation vectorially is

$$\vec F_{i,j} = \frac {Gm_im_j}{||\vec r_j - \vec r_i||^3}(\vec r_j - \vec r_i)$$

where

• $$\vec F_{i,j}$$ is the gravitational force exerted on body i by body j.
• $$G$$ is the universal gravitational coefficient.
• $$m_i$$ and $$m_j$$ are the masses of bodies i and j.
• $$\vec r_i$$ and $$\vec r_j$$ are the position vectors of bodies i and j expressed in some inertial coordinates and with respect to the origin of some inertial reference frame, which is typically the ensemble center of mass.

Acceleration rather than force is needed to propagate state. On combining the above with Newton's 2nd law the mass of body i cancels out. The component of the acceleration on body i as induced by body j is thus

$$\vec a_i = \frac {\mu_j}{||\vec r_j - \vec r_i||^3}(\vec r_j - \vec r_i)$$

Here, $$\mu_j$$ is the "standard gravitational coefficient" for body j. Conceptually, this is just the product of the universal gravitational coefficient G and the mass of body j. However, solar system astrophysicists *never* use this product, G*m_j. The reason is twofold. First, the standard gravitational coefficients for many of the bodies in the solar system are observable quantities. The masses of the Sun and the planets are not directly observable. The second reason is precision. The uncertainty in G is a lousy 1 part in of 10,000. The standard gravitational parameters for the Sun and the eight planets are known to much greater precision (e.g., 1 part in 10¹⁰ for the Sun). See [thread=271501]this thread[/thread] for details.

To find the total acceleration of body i, just sum the accelerations induced by all of the other bodies in the simulation:

$$\vec a_i = \sum_{j\ne i} \frac {\mu_j}{||\vec r_i - \vec r_j||^3}(\vec r_i - \vec r_j)$$

 

[Chronos]

Apologies, I do not see a simple n=3 solution in this mix. I confess to being dense.

 

[choppu]

Hi and thanks D H,



I was attempting to use a numerical solution on this problem, however I do not see how to begin. How to attempt to solve these two differential equations numerically (by example eulers method) ?


All you got is the equation for the accelerations...

 

[D H]

All you got is the equation for the accelerations...

What else do you need? Remember, acceleration is the derivative of velocity, and velocity is the derivative of position.

 

[choppu]

So i need to integrate the equations two times?

 

[D H]

So i need to integrate the equations two times?

You will definitely need to integrate N*6 scalar states, where N is the number of bodies. If you want to look at this as integrating twice, then the answer is "yes".

However, a second-order differential equation can always be reduced to a first-order differential equation. For example, Newton's second law is a second order differential equation. Ignoring the vectorial nature of things for a bit, Newton's second law for a constant-mass object subject to some external force is

$$\ddot x(t) = \frac{F(t)}{m}$$

To convert this to a first-order differential equation, construct a state vector u that comprises the object's position x and velocity v:

$$\vec u \equiv \bmatrix x \\ v \endbmatrix$$

The derivative of this vector is

$$\frac{d}{dt}\vec u(t) = \bmatrix u_2(t) \\ F(t)/m \endbmatrix$$

where u₂ is just the second element (the velocity element) of the vector u.

Bottom line: Augmenting the state with the velocity changed the initial second order differential equation to a first order differential equation.¹ This technique also works for vectors as well as scalars. Instead of adding 1 element to the state, you simply add n elements to the state -- the time derivatives of each of the n elements of the vector. This is exactly what you are going to need to do if you want to use a Runge-Kutta method to propagate the state of the system.

To answer the question posed in your last post: If you envision the position & velocity vectors for all N bodies in your simulation as forming one big vector with N*6 elements, you will only have to integrate once.²There is a problem with this augmented state approach. Suppose you use this approach and you use Euler's method to perform the integration. Your simulated solar system won't hang together for long. This basic use of Euler's method fails to conserve energy. There is a simple solution variously called symplectic Euler, Euler–Cromer, semi-implicit Euler (and others). Calculate all of the accelerations (see the second footnote). Then for each body, integrate the velocity vector using the calculated acceleration for that body and then integrate the position vector using the integrated velocity.³ This simple change makes a huge difference. Energy is conserved and orbits are bounded (more-or-less).⁴

==============

¹ This technique works for any order differential equation. For example, see http://en.wikipedia.org/wiki/Ordinary_differential_equation#Reduction_to_a_first_order_system.

² One advantage of this "big honkin state vector" approach is that it ensures that the accelerations are time synchronous. If you integrate each body separately you will need to ensure that you calculate all of the accelerations before you propagate any of the positions.

³ Using the Euler method to integrate the augmented state essentially does this in reverse order: It integrates the position vectors using the current velocity vectors and then integrates the velocity vectors (but using the accelerations calculated from the position vectors prior to the update).

⁴ If you want to know more about these kinds of techniques, google the terms "symplectic integrator" and "geometric integrator". Just as basic Euler method is a starting point for numerical integration in general, the symplectic Euler method is a starting point for integration methods that conserve some quantity of interest (e.g., energy) in some dynamic system.

 

[choppu]

Thank you so much D H, your post is really helpful.

However I still don't understand how I can figure out the x and y coordinates... as in the "current" x and y position. Because all there is now is that first order Differential equation. How do I figure out the X and Y coordinates? Is there a way to split them up or some other way?


Thanks

 

[Redbelly98]

You'll need to supply initial values for x, y, v_x, v_y for each planet. Then use a numerical procedure (eg. Runge-Kutta) to work out how x and y evolve over time for each planet.