# Anomaly in 2-body simulation error

I've generated a solution for a planet motion around the Sun (placed at the reference frame center). I wanted to test what would I get for the error of the method I'm using (Runge Kutta with 4 stages, with semimajor axis and Sun mass as distance and mass units, and time units evaluated in order to have $$G=1$$) by evaluating it with Richardson extrapolation, as described here, using specifically the perihelion of the orbit. I got this

What I actualy do is this (I'm using matlab)

p = 4
pers_001 = r_001(islocalmin(r_001));
pers_002 = r_002(islocalmin(r_002));

T_pers = (pers_002 - pers_001)./(2^(p+1) - 1);

figure;plot(1:length(T_pers),T_pers)
grid on


r_001 and r_002 are the position vector of the planet calculated with h/2 and h respectively, with h being the time-step used of $$0.001$$ and $$0.002$$. pers_001 and pers_002 takes the perihelion positions present in the two different position vectors. T_pers is the truncation error for every perihelion and it is presented in the picture above. This result seems good to me if I look at the order of the error, but that "modulated spikes thing" is something I didn't expect, what could be the reason? I mean I suppose the method's code is not wrong, since the general oscillation in perihelion position is of order $$10^{-7} A.U.$$, as can be seen from the plot below

EDIT

@Lutz Lehmann, yes the equation I'm using is precisely what you wrote, there's just the Sun mass at numerator also, but it's very close to one anyway. Plus I'm considering three dimensions, so for instance I use $$\ddot{\vec{r}} = M_{\odot}\vec{r}/\lvert r \rvert^3$$. I tried to verify what would have happened by changing $$h$$, and I actually got this (truncation errors)

The last picture is the same as the one before for the truncatin error, but the subplot ruins the details a bit. It seems actually something related to the aliasing since the amplitude of the effect decrease with $$h$$. So for some reason I periodically have a "bad" approximation of the perihelion (in respect of the other estimations).

• You're only showing one value per orbit. But surely your time step is smaller than one orbit. How do the curves look like with better resolution? Feb 4, 2021 at 5:24
• This could be artefacts of the discretization, an aliasing effect. I think you could observe similar if you determined the local extrema of the sine function in the same way. So the effect depends on how far away the time points $nh$ and $m(2h)$, where the local minimum of the computed sequence are found, are from the actual minimum, that is, the root of the radius velocity. And then also if $n=2m$, as $n=2m\pm1$ should be quite possible, and $n=2m\pm2$ can not be directly excluded. Feb 4, 2021 at 9:27
• Could you please write down the system you used? From what you wrote I can read off that you use something similar to $\ddot z=-\frac{z}{|z|^3}$ with $z(0)=(r_0,0)$ and $\dot z(0)=(0,v_0)$ with $r_0$ around $5$ by the second plot, and $v_0$ computed for some desired excentricity Feb 4, 2021 at 9:58
• @WolfgangBangerth Actually I'm showing just the perihelions, so there can't be more than one value per orbit. Unless I didn't understand what you mean.
– Zebx
Feb 4, 2021 at 13:40

The effect is due to quantization noise and/or aliasing, you are not computing with the true minima of the radius along the orbit, but with the closest sample point. This means that the sample point for the half-step integration can lie in-between. Close to the minimum this results in an $$O(e·dt^2)$$ ($$e$$ the eccentricity) difference between the computed radius minima just from the different time positions, in addition to the truncation errors of size $$O(dt^4)$$.

This all falls away if you use a better approximation of the true perihelion. (Or could be reduced by comparing positions at the same times points, even if they are not the local minima for both series.)

I have constructed a similar situation. The severity of the aliasing effect depends very much on the parameters of the orbit. By chance the first parameter I tried produced such a modulated oscillation as in the question, other parameters produced much smoother graphs.

First find a fixed-step solution of a slightly eccentric ellipse:

def fun(t,u):
x,y,vx,vy = u
r3 = (x*x+y*y)**1.5
return np.array([vx,vy,-x/r3,-y/r3])

r0, v0 = 1.0, 0.95 # v0=0.9 or 0.96 do not give short-period oscillations
u0 = [r0,0,0,v0]

t = np.arange(0,300,0.01)
u = RK4integrate(fun,t,u0)
x,y,vx,vy = u.T


Then plot the local minima

r = np.hypot(x,y)
per1 = [k for k,rk in enumerate(r[:-1]) if k>0 and r[k-1]>=rk and r[k+1]>=rk]
plt.plot(r[per1]); plt.grid(); plt.show()


Now the minima of the radius are at roots of the derivative of the radius function, that is, where $$r\dot r=x\dot x+y\dot y$$ changes its sign. Use a linear approximation to determine the root location between the sample points closer, then use one RK4 step to get the solution at this intermediate time with the same overall accuracy

dr = x*vx+y*vy

solve dr[k]*(1-s)+dr[k+1]*s = 0
return the radius of the RK4 step to this time
'''
s = dr[k]/(dr[k]-dr[k+1])
us = RK4integrate(fun,[t[k], (1-s)*t[k]+s*t[k+1]], u[k] )
xs,ys,vxs,vys = us[-1]
return np.hypot(xs,ys)
rmin = [ minrad(k) for k in range(len(x)-1) if dr[k]<=0 and dr[k+1]>=0]
plt.plot(rmin); plt.grid(); plt.show()


This returns the rather reduced plot

Note that the magnitude of the difference over the evolution is reduced from $$10^{-6}$$ in the first plot to $$10^{-9}$$ here.