Three point particles(m1,m2,m0) of mass 6, 8, 10 are located at
(1,3) the green dotted one ...
(-2,-1) the red one
( 1,-1) the blue one
These initial positions were denoted by a green circle in the last
frame of the movie.
The method used to solve the nonlinear system of equations will be the Runge-Kutta
4th or higher. In matlab there is a 4-5 ode solver that does adaptive step size.
The code comes from L.F.Shampine et. al. and probably
follows the general derivation from Taylor series expansion. In what follows y0dotx
represents first derivative of y with respect to x evaluated at x0 and y0dotxdotx
is the second derivative of y with respect to x evaluated at x=x0.
The Taylor expansion of y(x + h) is
y(x0 + h) = y(x0) +
ydotx(x0)*h + (1/2!)*ydotxdotx(x0) * h^2 + (1/3!)*ydotxdotxdotx(x0) * h^3
+ (1/4!)*ydotxdotxdotxdotx(x0) * h^4 + ...and so on.
We want to find a numerical solution to ydotx=f(x,y) as starting point with
( initial condition y(x0)=y0 ) so differentiate and substitute into the Taylor series
obtaining intermediate equation
ydotxdotx(x) = fdotx(x,y) + ydotx(x)*fdoty(x,y).
Substituting these into the Taylor expansion we get
y(x0+h)=y(x0) + f(x0,y0)*h + [ fdotx(x0,y0) + fdoty(x0,y0)*f(x0,y0)]*h^2 + higher terms
We want to write the right hand side of the equation as:
y(x0 + h)=y(x0) + A*h*f(x0,y0) + B*h*f[x0+C*h,y0 + D*h*f(x0,y0)]
and then find values of A,B,C,D so that this is true. Using the Taylor
series expansion of a function of two variables about x=x0 + a and y = y0 + b
f(x0+a,y0+b)=f(x0,y0) + a*fdotx(x0,y0) + b*fdoty(x0,y0) +....
Letting a=C*h and b=D*h*f(x0,y0) then
f[x0+C*h,y0 +D*h*f(x0,y0)]=f(x0,y0) + C*h*fdotx(x0,y0) + D*h*f(x0,y0)*fdoty(x0,y0)+...
substituting into the right hand side of y(x0+h) equation:
y(x0+h)=y(x0) + (A+B)*h*f(x0,y0) +B*h*h*[C*fdotx(x0,y0) +D*f(x0,y0)*fdoty(x0,y0)]
Compare coefficients, one choice is A=1/2=B , C=1=D and this gives the second-order
formula Runge-Kutta as
y(x0+h)=y(x0)+(1/2)*h*f(x0,y0)+(1/2)*h*f[x0+h,y0+h*f(x0,y0)]
This can be written as y(x0+h) = y(x0) + (1/2)*(u1+u2)
u1=h*f(x0,y0) and u2=h*f(x0+h,y0+u1).
If the first five terms of the original Taylor series are used we get fourth order:
y(x0+h)=y(x0)+ (1/6)*(w1 + 2*w2 + 2*w3 + w4) + ..terms of order h^5 neglected
w1=h*f(x0,y0)
w2=h*f(x0+h/2,y0 + w1/2)
w3=h*f(x0+h/2,y0+w2/2)
w4=h*f(x0+h,y0+w3)
The input to ode45 is the functions (contained in function orbit3.m)
in vector form that describes the equations
to be solved , the time span of interest, and the initial conditions.
The general three body problem in space involves 3 times 6 (velocity and
acceleration in 3 dimensions)equations but restricting
the center of mass to lie at the origin and conservation of momentum
cuts the equations down by 6 so the planar solution needs 12 equations(or less
as will be mentioned later these lead to some messes).
For the three body problem the initial conditions are 12
component vector consisting of the x,y initial
positions and the x,y initial velocities and these four are repeated three
times to get twelve. The output of ode45 [t,Y] is a twelve component row
vector for each time in the timespan considered. The time is a colomn vector.
The order of the row vector
output from ode45 is the same order as the IV (initial values for starting off the
calculation):x0,y0,x0dott,y0dott,x1,y1,x1dott,y1dott,x2,y2,x2dott,y2dott
Absorbing the gravitational constant into the masses, the x component
of the force exerted by m1 onto m0 divided by the mass of m0 is
x0dottdott=m1*(x1-x0)/r301 with
r301=(sqrt((x1-x0)^2 + (y1-y0)^2))^3
The y component of the force is m1*(y1-y0)/r301
Of course there is a force on m0 from m2 and it needs to be added to
produce the x component of the force divided by m0 or the x component of the
acceleration of x0 is
x0dottdott= m1*x10 -m2*x02 with the definitions
x10=(x1-x0)/r301
x02=(x0-x2)/r302
r302=(sqrt((x0-x2)*(x0-x2)+(y0-y2)*(y0-y2)))^3
The y component of the force on x0 divided by the mass of m0 or the
y component of the acceleration of m0 is
y0dottdott=m1*y10 - m2*y02
y10=(y1-y0)/r301
y02=(y0-y2)/r302
The x and y components of the acceleration of x1 are
x1dottdott=m2*x21 -m0*x10 with the following definitions
r312=(sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)))^3
x21=(x2-x1)/r312 the y component is
y1dottdott=m2*y21 - m0*y10 with the definition
y21=(y2-y1)/r312
The x and y components of the acceleration of x2 are
x2dottdott=m0*x02 - m1*x21 with previous definitions
y2dottdott=m0*y02 - m1*y21
There are six equations for the acceleration components of three bodies in two
dimensions, and these equations are coupled non-linearly.
The input of orbit3.m is the output of last call to ode45. The output of orbit3.m
is a twelve component colomn vector with the components in the order:
x0dott,y0dott,x0dottdott,y0dottdott,x1dott,y1dott,x1dottdott,y1dottdott,
x2dott,y2dott,x2dottdott,y2dottdott.
forming the quantities x10=(x1-x0)/r301 and y10 in similar manner
the twelve components of the output of the function orbit3 are
f(1)=Y(3); % x0dot
f(2)=Y(4); %y0dot
f(3)=m1*x10 -m2*x02 ; % x0dotdot
f(4)=m1*y10 -m2*y02 ; % y0dotdot
f(5)=Y(7); %x1dot
f(6)=Y(8) ;%y1dot
f(7)=m2*x21 -m0*x10 ;% x1dotdot
f(8)=m2*y21 - m0*y10; % y1dotdot
f(9)=Y(11); %x2dot
f(10)=Y(12);%y2dot
f(11)=m0*x02 - m1*x21 ;% x2dotdot
f(12)=m0*y02 - m1*y21; % y2dotdot
All programs are listed at the end of the page.
The matlab program ode45 has a default absolute error 1e-6 and absolute relative error
1e-3. If the initial velocity of m1 instead of being zero is altered to
be (0,-1e-5) changes in the orbit can be seen.
Even with m1 initial velocity (0,-1e-6)
a visible orbit change can be seen.
m1 initial velocity (0,-1e-7) produces no visible orbit change
as you would expect from 1e-7 being less
tha 1e-6 (the default error setting of ode45 program).
ThreePlay.m calls ode45 (matlab's numerical differential equation solver) with function orbit3.m and generates low resolution jpeg files that are run through the javascript animator. You can see this crude display by pressing start and scrolling down a bit. This still needs some work to find the bug that doesn't allow the animation to stop. At the present time ode45 solves for masses m0=5, m1=3, m2=4 subject to gravitational forces initial positions and velocities are in IV=[1,-1,0,0,1,3,0,0,-2,-1,0,0] and the time is only from 0 to 20. The default (.000001) accuracy is used. Using masses m0=10,m1=6,m2=8 we were able to increase the length of the timespan. D.Gruntz and J.Waldvogel in Chapter 4 [1]Orbits in the Planar Three-Body Problem page 41 get tighter loops in the orbits of all three from time 10 to 20. They use 10 million times the accuracy. At t=15.8299 there is a near collision of m0 and m2. This affects the accuracy from then on. The final evolution of all three orbits is incorrect.
Introducing new variables such that singularities due to
all possible binary collisions are regularized leads to more believable orbits.
Considering conservation of energy and angular momentum reduces the equations required to
eight instead of twelve. New variables are introduced that depend on the
relative positions of the three bodies.
From before we have
(s1)^3=r302=(sqrt((x0-x2)*(x0-x2)+(y0-y2)*(y0-y2)))^3
(s2)^3=r301=(sqrt((x0-x1)*(x0-x1)+(y0-y1)*(y0-y1)))^3
(s0)^3=r312=(sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)))^3
and define xs1=x0-x2 ys1=y0-y2 as the x and y components of coordinate s1
also xs2=x1-x0 ys2=y1-y0
ans xs0=x2-x1 ys0=y2-y1
m=m0 + m1 + m2
define xG=(xs1/r302) + ( xs2/r301) + (xs0/r312)
similarl for yG so that the equations of motion look more symmetrical
xs1dottdott=-(m*xs1/(s1)^3 +m1* xG is the second derivative of the x component of s1
ys1dottdott=-(m*ys1/(s1)^3 +m1* yG is the second derivative of the y component of s1
xs2dottdott=-(m*xs2/(s2)^3 +m2* xG is the second derivative of the x component of s2
ys2dottdott=-(m*ys2/(s2)^3 +m2* yG is the second derivative of the y component of s2
xs0dottdott=-(m*xs0/(s0)^3 +m0* xG is the second derivative of the x component of s0
ys0dottdott= -(m*ys0/(s0)^3 +m0*yG is the second derivative of the y component of s0
Using conservation of energy and momentum the equations of motion are alot harder to
write down but can be done; this is most useful at the points of near collisions
as you can tell from the movie (time around 16 is critical).
After all their work (transforming coordinates at near collisions,etc)
, Gruntz and Waldvogel conclude that their accuracy
was only 1e-6 and eventually the lightest mass does escape leaving behind a binary.
Forces proportional to distance are dual to
forces proportional to inverse distance squared.
Z -> z ^ 2 transforms orbits of F to those of f [2,3]. After they made these changes and computed they found that there are some periodicities. The three body problem has been solved exactly by K.F.Sundman, Memoire sur le probleme des trois corps, Acta Math 36,1912, pages 105-179. and T.Levi-Civita,Sur la regularisation du probleme des trois corps, Acta Math. 42,1920,pp99-144. It would be nice to plug these equations into a computer and see what
obits these series solutions predict for the Pythagorean three body and others. Real three bodies are not point masses so this problem is not too realisitic.
[1]W.Gander and J.Hrebicek, Solving Problems in Scientific Computing Using Maple and Matlab isbn 3-540-58746-2
[2] The Kasner-Arnold theorem from page 247 of Visual Complex Variables Associated with each power law F~R^A is precisely one power law f~r^a that is dual in the sense that orbits of the former are mapped to orbits of the latter by Z -> z ^ m and the relationship between the forces and the mapping are (A + 3)*(a + 3) = 4 and m=(A + 3)/2. In general, positive energy orbits in either the attractive or repulsive field F ~ R ^ A map to attractive orbits in the dual field while negative energy orbits map to repulsive ones. However if -3 < A < -1 (e.g. gravity) then these roles are reversed. In all cases, zero energy orbits map to force-free rectilinear orbits.
[3] http://www.math.psu.edu/hall/newton/alternative.pdf
%threePlay.m global m0 m1 m2 THE CODE IS DOWN FOR REJUVENATION
