3
$\begingroup$

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

enter image description here

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

p = 4
r_001 = radius(J_001,1);
r_002 = radius(J_002,1);
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

enter image description here

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)

enter image description here

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).

Improve this question
$\endgroup$
4
  • $\begingroup$ 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? $\endgroup$
    – Wolfgang Bangerth
    Commented Feb 4, 2021 at 5:24
  • 1
    $\begingroup$ 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. $\endgroup$
    – Lutz Lehmann
    Commented Feb 4, 2021 at 9:27
  • $\begingroup$ 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 $\endgroup$
    – Lutz Lehmann
    Commented Feb 4, 2021 at 9:58
  • $\begingroup$ @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. $\endgroup$
    – Zebx
    Commented Feb 4, 2021 at 13:40

1 Answer 1

Reset to default
3
$\begingroup$

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()

enter image description here

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

def minrad(k):
    '''minimal radius of segment k
    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

enter image description here

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

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.