# Solve implicit ODE numerically in orbit simulation

I'm trying to plot the orbit of a compact binary star system where general relativistic effects become important. I'm using post-Newtonian approximation and I want to solve the orbit numerically based on the equation below:

\begin{aligned} \ddot{\vec{r}}_{i} &=& \sum_{j\neq i}^{ } \frac{\mu_{j}(\vec{r}_{j}-\vec{r}_{i})}{r_{ij}^{3}}\left[1-\frac{4}{c^{2}}\sum_{k\neq i}^{ }\frac{\mu_{k}}{r_{ik}}-\frac{1}{c^{2}}\sum_{k\neq j}^{ }\frac{\mu_{k}}{r_{jk}}+\frac{\vec{v}_{i}^{2}}{c^2}+2\frac{\vec{v}_{j}^{2}}{c^2}-\frac{4}{c^2}\vec{v}_{i} \cdot \vec{v}_{j} - \frac{3}{2c^{2}}\left(\frac{(\vec{r}_{i}-\vec{v}_{j})\cdot\vec{v}_{j}}{r_{ij}}\right)^{2}+\frac{1}{2c^{2}}(\vec{r}_{j}-\vec{r}_{i})\cdot\ddot{\vec{r}}_{j}\right] \\&& + \frac{1}{c^2}\sum_{j\neq i}^{ }\frac{\mu_{j}}{r_{ij}^{3}}\left[(\vec{r}_{i}-\vec{r}_{j})\cdot(3\vec{v}_{i}-3\vec{v}_{j})\right](\vec{v}_{i}-\vec{v}_{j}) + \frac{7}{2c^{2}}\sum_{j\neq i}^{ }\frac{\mu_{j}\ddot{\vec{r}}_{j}}{r_{ij}}, \end{aligned} where $\mu_{i}=Gm_{i}$ and $G$ is gravitational constant, $\vec{r}_{i}$,$\vec{v}_{i}$ and $\ddot{\vec{r}}_{i}$ are the position vector, velocity vector and acceleration vector of star $i$, and $r_{ij} = |\vec{r}_{i}-\vec{r}_{j}|$.

and this is listed as Eq. 1 in this paper: http://arxiv.org/pdf/0709.1160v2.pdf

As you can see on the right hand side of this equation there is $\ddot{\vec{r}}_{j}$ term which makes finding $\ddot{\vec{r}}_{i}$ not straightforward. There are two ways I could think of:

1. In the code, just plug in the value of $\ddot{\vec{r}}_{j}$ that was obtained at one time-step ago, but I assume this would give large errors after a long simulation time.
2. Use only Newtonian part as the value of $\ddot{\vec{r}}_{j}$ to compute $\ddot{\vec{r}}_{i}$ but this gave me a very messy plot, and to improve it I used the value of $\ddot{\vec{r}}_{i}$ that I just calculated to be the new value of $\ddot{\vec{r}}_{j}$ and use this to obtain a new $\ddot{\vec{r}}_{i}$, although the plot I got is indeed better now but it's still far from accurate. Maybe I could do the same thing recursively but I don't know how much accurate it can get and if it's worth the computation time.
3. Convert this to explicit expression of $\ddot{\vec{r}}_{i}$ but this expression is too complicated to do so.

Any suggestions on how to numerically solve this $\ddot{\vec{r}}_{i}$ correctly?

(I posted a similar question in Astronomy but here might be a more suitable place to ask this question)

• Cross-posted: astronomy.stackexchange.com/q/11338 Please avoid posting the same question on multiple SE sites, it is considered a breach of etiquette: meta.stackexchange.com/q/64068 For a question where the main issue is solving a set of equations numerically, scicomp.SE is likely the appropriate place to post it. – Kirill Jul 22 '15 at 0:51

(1) Using the previous value of $\ddot{r}_j$ is like adding an error term to the r.h.s. of your equation of magnitude $\mathit{const}\times(\ddot{r}_j(t+\delta t) - \ddot{r}_j(t))$, meaning your scheme will only be first-order correct, regardless of whether the integration method you use has a higher order.
(2) I believe what you are describing is similar to applying the iteration method to solve the equation in the unknowns $\ddot{r}_i$. But what you're describing differs from how the iteration method would normally be applied because normally you would iterate until convergence.
(3) This is the best approach here. The system of equations is linear in the unknowns $\ddot{r}_i$. If I denote by $\ddot R$ the vector whose components are the components of all $\ddot{r}_1,\ldots,\ddot{r}_n$, what you can do is rewrite the system of equations in the form $$\ddot R = f(r, v) + A(r, v)\ddot R,$$ where $A(r,v)$ is a simple, but not diagonal (because of the $(r_j-r_i)\cdot \ddot{r}_j$ term), matrix, and then just solve this as a system of linear equations. It only needs to be linear in the unknowns, after all, and you do not need to solve the linear system by paper and pencil, only construct it and hand it to a linear solver (which would use LU/GE, for example).