I'm working on a programm to simulate the dynamics of a system of $n$ point charges using the solutions to the Liénard-Wiechert potential . For $n <= 2$ the forward euler method has served me well, but for bigger $n$ it quickly diverges.
For a given frame the accelaration is computed followingly ( $ Q_i = Q_1,...Q_n$): $$ \mathbf a(\mathbf r) = \frac q m( \mathbf E(\mathbf r, \mathbf x_i,\mathbf v_i,\mathbf a_i) + \mathbf v \times \mathbf B(\mathbf r ,\mathbf x_i,\mathbf v_i,\mathbf a_i)) $$ (Excuse my Latex knowledge, I dont know how to type multiple parameters) Since the acceleration depends not only on position but on velocity and acceleration itself, common methods, like leapfrog, do not apply (I believe). So I would like to know which method would be best suited for this task. Since I'm new to this field, advice on literature would also be much appreciated.
I'll post the code on request.
The reason the Euler method loses precision very rapidly (it doesn't actually diverge, but the solutions become worthless) is that for $n\geq 3$ solutions to the $n$-body problem generally exhibit chaotic behavior (that is, the behavior at time $T \gg t$ ids sensitive to exponentially small (in $T-t$) changes in the initial conditions at time $t$. This causes any algorithmic imprecision that would be polynomial in time step size to result in non-polynomial growth of the imprecision at large $T$. This is the "divergence" you observe.
There are some very special solutions for which this property does not hold, and there you can sensibly use the Euler method.
If you improve the order of the method (and you can develop an analogue of leapfrog or Runge-Kutta for this problem) this basic property of the solution still causes that approximation to lose precision at worse than a polynomial rate.