So, just finished reading the Three Body Problem by Cixin Liu (very good novel), and was curious about the two-body problem. Of course, the two-body scenario is much simpler and more predictable than the two-body problem, so I was hoping to find a general solution to the two-body problem but haven't found one online, just solutions of specific cases. I decided to try to go at it on my own for the 2-dimensional system, but I got stuck.
I started by saying the first planet, Planet $x$, had a 2-dimensional position vector $\vec{x}$ and a mass $m_x$, and the second (Planet $y$) had a position vector of $\vec{y}$ and mass $m_y$. With the generalization that the two bodies are point masses with the only force being gravity, and assuming that both the initial positions/velocities and the masses of the two bodies are known, I came up with these two second-order differential equations:
$$\frac{d^2\vec{x}}{dt^2}=\frac{Gm_y}{||\vec{x}-\vec{y}||^3}(\vec{y}-\vec{x}) \\
\frac{d^2\vec{y}}{dt^2}=\frac{Gm_x}{||\vec{x}-\vec{y}||^3}(\vec{x}-\vec{y}) \tag{1}$$
I got these from these classical physics equations:
$$\underbrace{F_{y \rightarrow x}}_\text{force of gravity on x} = \underbrace{\frac{G m_y m_x}{||\vec{x}-\vec{y}||^2}}_\text{gravitational force}=\underbrace{m_x\frac{d^2\vec{x}}{dt^2}}_\text{mass and acceleration}$$
Then, I saw that there was a lot similar in the terms $\frac{d^2\vec{x}}{dt^2}$ and $\frac{d^2\vec{y}}{dt^2}$, so I used substitution to get this relationship:
$$m_x\frac{d^2\vec{x}}{dt^2} + m_y\frac{d^2\vec{y}}{dt^2} = 0$$
And that makes sense with the equal and opposite forces. Then I integrated twice:
$$m_x\vec{x} + m_y\vec{y} = \vec{c}t + \vec{d} \tag{2}$$
… where $\vec{c}$ and $\vec{d}$ are vector integration constants (I realized that they'd need to be vectors and not just constants). Using initial conditions, this means:
$$\vec{c} = m_x\vec{x}'(0) + m_y\vec{y}'(0) \\
\vec{d} = m_x\vec{x}(0) + m_y\vec{y}(0) \tag{3}$$
I solved $(3)$ for $\vec{y}$ (and later for $\vec{x}$) to substitute into $\frac{d^2\vec{x}}{dt^2}$ and $\frac{d^2\vec{y}}{dt^2}$ and uncouple the two differential equations in $(1)$:
$$\frac{d^2\vec{x}}{dt^2}=\frac{-Gm_y^3}{||(m_x+m_y)\vec{x}-\vec{c}t-\vec{d}||^3}((m_x+m_y)\vec{x}-\vec{c}t-\vec{d}) \\
\frac{d^2\vec{y}}{dt^2}=\frac{-Gm_x^3}{||(m_x+m_y)\vec{y}-\vec{c}t-\vec{d}||^3}((m_x+m_y)\vec{y}-\vec{c}t-\vec{d}) \tag{4}$$
But this is where I got stuck. I coded something in Matlab that uses $(4)$ in the ode45 function to see how things move at least (code below), but I definitely have not been able to find any methods to solve either equation in $(4)$. Does such a closed form even exist? Of course that would mean that $\vec{x}$ and $\vec{y}$ aren't always functions, but I was wondering if maybe there's a parameterized solution or something like that. Thanks!
Here's the code for my demonstration, it's kinda janky:
% %%%%%%%%% Set stuff up %%%%%%%%%
% planet x
x_init = [6 9]';
dx_init = [0 -0.001]';
mx = 6900;
% planet y
y_init = [-6 -9]';
dy_init = [0 0]';
my = 420000;
% other things
tlength = 100000; % diffeq soln time length (bigger means longer duration)
scale = 300; % scale of velocity vectors
speed = 5; % speed of animation (bigger means faster)
% %%%%%%%%% Animate Motion %%%%%%%%%
plotPlanets([x_init dx_init y_init dy_init],[mx my tlength],scale,speed)
%% Functions
function [T,D] = getPositions(vecs, cons)
% vecs := matrix containing info for the planet in this format:
% [[x0] [dx0] [y0] [dy0]]
% cons := vector with relevant constants in this format:
% [[mass of x] [mass of y] [time length]]
syms x1(t) x2(t) dx1(t) dx2(t) y1(t) y2(t) dy1(t) dy2(t)
funcx = formula([x1 x2 dx1 dx2]);
funcy = formula([y1 y2 dy1 dy2]);
% Unpack
x_init = vecs(:,1);
dx_init = vecs(:,2);
y_init = vecs(:,3);
dy_init = vecs(:,4);
mx = cons(1);
my = cons(2);
tlength = cons(3);
c = mx*dx_init + my*dy_init;
d = mx*x_init + my*y_init;
eqnsx = getFormulas([c d],[mx my],[funcx t]);
eqnsy = getFormulas([c d],[my mx],[funcy t]);
eqns = [eqnsx eqnsy];
vars = [funcx funcy];
init = [x_init' dx_init' y_init' dy_init'];
F = odeFunction(eqns,vars);
[T, D] = ode45(F,[0 tlength],init);
end
function eqns = getFormulas(vecs, mass, funcs)
% vecs := matrix containing info for the planet in this format:
% [[c] [d]]
% mass := vector with relevant constants in this format:
% [[own mass] [other mass]]
% funcs := symbolic functions in the format [z1 z2 dz1 dz2 t].
% Set up
G = 6.67408*10^(-11); % Gravitational constant
% Unpack
c = vecs(:,1);
d = vecs(:,2);
mz = mass(1);
mw = mass(2);
% fromz is vector between planets
fromz = formula((mz+mw)*[funcs(1);funcs(2)] - c*funcs(5) - d);
eqn1 = -G*mw^3/norm(fromz)^3*fromz(1);
eqn2 = -G*mw^3/norm(fromz)^3*fromz(2);
eqns = formula([funcs(3:4) eqn1 eqn2]);
end
function plotPlanets(vecs, cons, scale, speed)
[T,D] = getPositions(vecs,cons);
Dx = D(:,1:4);
Dy = D(:,5:8);
xmin = min([min(Dx(:,1)) min(Dy(:,1))]);
xmax = max([max(Dx(:,1)) max(Dy(:,1))]);
ymin = min([min(Dx(:,2)) min(Dy(:,2))]);
ymax = max([max(Dx(:,2)) max(Dy(:,2))]);
tlast = T(1);
figure(1)
hold on
axis(sizeUp(xmin,xmax,ymin,ymax))
pathx = plot(1);
px = plot(Dx(1,1),Dx(1,2),"bo");
vx = quiver(Dx(1,1),Dx(1,2),scale*Dx(1,3),...
scale*Dx(1,4),"Color",[0 0 0]);
pathy = plot(1);
py = plot(Dy(1,1),Dy(1,2),"ro");
vy = quiver(Dy(1,1),Dy(1,2),scale*Dy(1,3),...
scale*Dy(1,4),"Color",[0 0 0]);
pause(0.5)
for i = 2:size(D,1)
if ~any(findobj(allchild(0), 'flat', 'type', 'figure') == 1)
break
end
delete(pathx);
delete(px);
delete(vx);
delete(pathy);
delete(py);
delete(vy);
pathx = plot(Dx(1:i,1),Dx(1:i,2),"b");
px = plot(Dx(i,1),Dx(i,2),"bo");
vx = quiver(Dx(i,1),Dx(i,2),scale*Dx(i,3),...
scale*Dx(i,4),"Color",[0 0 0]);
pathy = plot(Dy(1:i,1),Dy(1:i,2),"r");
py = plot(Dy(i,1),Dy(i,2),"ro");
vy = quiver(Dy(i,1),Dy(i,2),scale*Dy(i,3),...
scale*Dy(i,4),"Color",[0 0 0]);
pause((T(i)-tlast)/speed/1000)
tlast = T(i);
end
if any(findobj(allchild(0), 'flat', 'type', 'figure') == 1)
hold off
end
end
function sizes = sizeUp(xmin, xmax, ymin, ymax)
if xmax - xmin > ymax - ymin
len = (xmax - xmin)/2;
ymid = sum([ymax ymin])/2;
ymin = ymid - len;
ymax = ymid + len;
elseif xmax - xmin < ymax - ymin
len = (ymax - ymin)/2;
xmid = sum([xmax xmin])/2;
xmin = xmid - len;
xmax = xmid + len;
end
pad = abs(len*0.1);
sizes = [xmin-pad, xmax+pad, ymin-pad, ymax+pad];
end
Best Answer
The easiest way to solve the two-body problem is to reduce it into two independent one-body problems by introducing the center of mass.
This effectively decouples the motion of the two bodies into a (1) center of mass motion (= 1st one-body problem) and a (2) displacement vector motion (= 2nd one-body problem).
For detailed info, you may wish to refer to this Wikipedia link: https://en.wikipedia.org/wiki/Two-body_problem#Center_of_mass_motion_(1st_one-body_problem)