Numerical errors in highly eccentric motion in gravitational field

Question

I am trying to solve the following problem:

A test particle at coordinate $\vec{r_0}=(1,0)$ is thrown with initial velocity $\vec{v_0}=(0,v)$ and is influenced by a gravitational field $\vec{F}=\frac{-\vec{r}}{r^3}$ (an $M=1$ mass at $(0,0)$) where $v$ can have the values: 1 (circ motion) , 0.1 or 0.01 (for eliptical motion). I need to solve this problem using RK-4 method for $T_{sol}=2\pi$ and use a timestep to reach a desired accuracy $\frac{dE}{E}$ (relative error of the energy at the end of the calculation and the starting,analytical energy), for all 3 values.

I already prepared the equations and the code, and I get the desired accuracy for v=1 with 60 points for the solution. However, so far I couldn't reach the same accuracy for v=0.1 or 0.01, and I have tried as many as 25,000 points. I see that in every solution, the motion begins as a fall towards the mass until it reaches some minimum radius, and then the trajectory gets ruined due to some numeric error, and the relative error in the end is not close at all to the desired number. Is there any way to estimate how many points I will need for the solution to prevent that from happening?

Answer 1

The Jacobian of the force, and thus the Lipschitz constant of the system, has a size of about $|r|^{-3}$. For $v=0.1$ the periapsis is $|r|=0.005$, giving a Lipschitz constant of about $L=10^{7}$. The stability of the RK4 method requires $Lh<2.5$, thus $h<h_s\simeq 2.5·10^{-7}$. For results approaching "qualitatively correct" with a fixed step-size method one would need step sizes smaller $h_s/10$. With a period of $T=2.222$ one would need about $10^8$ points per orbit period.

A step-size controller will keep the step-size times dominant eigenvalue of the Jacobian at the boundary of the stability region, that is, $h\approx h_s$, unless the error tolerance used is smaller than $c·h_s^4\sim c·10^{-27}$, where $c\sim |r|^{-6}\sim 10^{14}$ is the leading coefficient of the error expansion. Then the step size will be even smaller, increasing the run-time even further. Care must be paid to avoid the accumulation of floating-point errors in the summation of many small quantities.

For $v=0.01$ the apoapsis is at radius $|r|=0.0001$, which makes the stiff situation explained above even more untractable, $L=10^{12}$ and so on.