# Solving the Three-body problem numerically

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?

• Are the Euler-Lagrange equations a requeriment? To solve this kind of differential equations, Hamilton equations are easier because they're 1st order equations. By the way, this year I computed numerically a Solar System using a Verlet algorithm. It is fast and it gives an acceptable accuracy. Second thing. About your doubt: I recommend you to write explicitly the form of $|\vec{r_j}-\vec{r_k}|=\sqrt{(x_j-x_k)^2+(y_j-y_k)^2+(z_j-z_k)^2}$, then derive, then re-order the result in terms of $\vec{\dot{r_j}}$ and $\vec{\dot{r_k}}$. – VictorSeven May 23 '15 at 13:20
• I could try to solve this problem with Hamilton equations, thank you for your advice. Relatively to my doubt, I just found that $\frac{d}{dt}\left(|\vec{r_i}-\vec{r_j}|\right)=|\vec{\dot{r}_i}-\vec{\dot{r}_j}|\cos(\theta)\neq|\vec{\dot{r}_i}-\vec{\dot{r}_j}|$, where $\theta=\angle(\vec{\dot{r}_i}-\vec{\dot{r}_j},\vec{{r}_i}-\vec{{r}_j})$. I don't think that I could express $\cos(\theta)$ in terms of $r_i,r_j,\dot{r}_i$ or $\dot{r}_j$. – Élio Pereira May 23 '15 at 14:10
• But we can't simplify the problem using Hamilton equations, because generalized coordinates would be the same as in by Lagrangian equations. So, we still have to use the $x,y,z$ components of $\vec{s_i}$, $i=1,2,3$. – Élio Pereira May 23 '15 at 15:28

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$.