Solving the Three-body problem numerically

Question

I want to create a program in $Mathematica$ that solves numerically the Three-body problem by Euler-Lagrange's equations. I was searching some methods to sucessfully do it. So I found a way to solve Two-body problem in http://www.maths.usyd.edu.au/u/joachimw/thesis.pdf (page 12). I found, also, how to minimize the number of generalized coordinates using relative-position vectors $\vec{s_i}=\vec{r_j}-\vec{r_k}$ and center of mass $\vec{r_G}$, as you can see it in the following scheme scanned from Goldstein, 3rd edition:

With this, I reached the Lagrangin to a Three-body system:

$$\mathcal{L}=\frac{1}{2}M\dot{r}_{G}^2+\frac{1}{2}\frac{m_1m_2}{M}|\vec{\dot{r}_1}-\vec{\dot{r}_2}|^2+\frac{1}{2}\frac{m_2m_3}{M}|\vec{\dot{r}_2}-\vec{\dot{r}_3}|^2+\frac{1}{2}\frac{m_1m_3}{M}|\vec{\dot{r}_1}-\vec{\dot{r}_3}|^2+2G \left[\frac{m_1m_2}{|\vec{{r}_1}-\vec{{r}_2}|}+\frac{m_2m_3}{|\vec{{r}_2}-\vec{{r}_3}|}+\frac{m_1m_3}{|\vec{{r}_1}-\vec{{r}_3}|}\right]$$

I want to use Euler-Lagrange's equations:

$$\frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial{\dot{q}_i}}\right)-\frac{\partial \mathcal{L}}{\partial q_i}=0$$

where $q_i=r_G,s_1,s_2,s_3$.

Finding Lagrange's equations for $r_G$ was really easy, because $\frac{\partial \mathcal{L}}{\partial \dot{r}_{G}}$ is a conserved quantity ($r_G$ doesn't appear explicitly in $\mathcal{L}$). But finding Lagrange's equations for $s_i$, $i=1,2,3$ was a little bit confusing.

I have a doubt: is it true that $$\frac{d}{dt}\left(|\vec{{r}_j}-\vec{{r}_k}|\right)=|\vec{\dot{r}_j}-\vec{\dot{r}_k}|$$ or more specifically, $$\dot{s_i}=|\vec{\dot{r}_j}-\vec{\dot{r}_k}|~?$$ Because if it doesn't (and that's what I think) it wasn't useful using relative positions vectors. It would require new generalized coordinates as the $x,y,z$ components of $\vec{s_i}$.

What would be your approach to solve this problem? Could this be a good way to do it?

Answer 1

Well, I've done some calculus to your problem.

The problem, in fact, is that $|\vec{r_i}-\vec{r_j}|$ everywhere. Also the temporal derivative of that is pain. Of course you can write it without that angle you commented if you write the derivatives in terms of escalar products.

However, I'm not going to follow that. You have a system of 3 bodies, isolated, so energy is conserved and you can write the Hamiltonian of the system as:

$$H = \sum_i \frac{p_{x_i}^2+p_{y_i}^2+p_{z_i}^2}{2m_i} - \sum_{i,j>i}\frac{Gm_im_j}{|\vec{r_i}-\vec{r_j}|}$$

Note the sum is $j$ is with $j>i$ to avoid counting the same interaction twice. You could use the condition $i\neq j$ and after divide by 2, but doing this you avoid unneccesary iterations. It looks like the problem is still there. However, note that I have used as generlized coordinates simply the cartesian coordinates of each particle and not the center of mass positions.

Now let's apply Hamilton equation. Suposse you want to calculate the $\dot{p}_x$ of the particle $k$. Then you must solve:

$$\dot{p}_{x_k} = -\frac{\partial H}{\partial x_k} = \frac{\partial}{\partial x_k}\sum_{i,j>i}\frac{Gm_im_j}{|\vec{r_i}-\vec{r_j}|} = \sum_{j>k}Gm_im_j\frac{\partial}{\partial x_k}\frac{1}{|\vec{r_k}-\vec{r_j}|}$$

Note that in the last step I've eliminated the sum in $i$, because if $k\neq i$ then the derivative is 0. Thit last derivative is easy to calculate. You can see, writing the expression of $|\vec{r_k}-\vec{r_j}|$ with coordinates, as I indicated in the comments, that this derivative is:

$$\frac{\partial}{\partial x_k}\frac{1}{|\vec{r_k}-\vec{r_j}|} = \frac{x_j-x_k}{|\vec{r_k}-\vec{r_j}|} $$

This derivative is easier to evaluate than the one you have, because this is partial in the coordinates, but you have a total derivative in time. That's why Hamilton is better in this case. At the end, you'll have the following system of equations:

$$\dot{p}_{q_i}= \sum_{j>i}Gm_im_j \frac{q_j-q_i}{|\vec{r_i}-\vec{r_j}|}$$ $$\dot{q}_i = \frac{p_{q_i}}{m}$$

With $q\equiv x,y,z$. Now you only have to solve this system of 1st order equations with your favorite method (Euler, Runge Kutta, etc). Note that the sums can be to an arbitrary number of particles $N$.