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.pdfAs 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: 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.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.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)
Solve implicit ODE numerically in orbit simulation
ode;simulation;computational physics
(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.This is a very conventional idea described in numerical analysis textbooks in relation to nonlinear root finding.(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).Because the system is linear, the usual methods for solving nonlinear equations like iteration method or Newton's method aren't really necessary. If you do end up using an iterative method, then what you're describing is close to how it's usually done: you iterate until convergence, using either direct iteration, Newton's method, or something similar, and the guess used for the first iteration is something like the value from previous time step. There is already a lot of research literature on this in the context of differential-algebraic equations (DAE).