Celestial Mechanics: Solving Differential Equations for N-Body Systems

[rizumu]

I am getting stuck trying to solve a differential equation describing the gravitational relationship between two celestial bodies.

What I have so far is this:

a1 = G m2 (r2 - r1) / ||r2 - r1||^3

a2 = G m1 (r1 - r2) / ||r1 - r2||^3

where:
boldface denotes vectors in 3 dimensional Cartesian space
G for gravitational constant
m for mass
r for position
a for acceleration

I guess my question is where to go from here. I am assuming that the general solution will contain constants relating to initial position and initial velocity for both bodies. My goal is to arrive at a point where I can plot the trajectories of the two bodies using iterations. Essentially a system I can solve using numerical methods such as Euler's.

Anyway, any insight into how to solve this kind of problem would be great! Eventually I want to expand it to more than 2 bodies but I need to get comfortable with 2 first before taking the problem to N-bodies. This is for a numerical methods class, and possibly a final project/presentation.

Cheers!

 

[Petr Mugver]

The equations for the trajectories of the two bodies (not the equations of motion) can be solved exactly, without approximations. The key is to pass from the r1 r2 coordinates to the R r coordinates, where R are the coordinates of the center of mass of the system, and r are the relative coordinates of the two bodies:

R = (m1 r1 + m2 r2) / (m1 + m2)

r = r2 - r1

If you carry out the calculations, it turns out that R satisfies the equations for a free body of mass M = m1 + m2:

MR''=0

and r satisfies the equation of motion of a single body of reduced mass

m = m1 m2 / (m1 + m2)

in a gravitational central potential with center fixed at the origin:

r'' = -Gmr / |r|^3

So the whole problem is reduced to the one of a single body.

 

[phyzguy]

I think you are looking more for an iterative procedure that applies to more than two bodies. With more than two bodies, it can't be solved exactly. Given that you know the positions and velocities at time step n, and you have calculated the acceleration as above, you just time step forward to calculate the new positions and velocities. Then you calculate the new accelerations and iterate.

x_{n+1} = x_n +v_{xn}*\delta t; y_{n+1} = y_n +v_{yn}*\delta t; z_{n+1} = z_n +v_{zn}*\delta t; v_{x n+1} = v_{xn} +a_{xn}*\delta t; v_{y n+1} = v_{yn} +a_{yn}*\delta t; v_{z n+1} = v_{zn} +a_{zn}*\delta t