效果是由于量化噪声和/或混叠造成的，您不是使用沿轨道的半径的真实最小值进行计算，而是使用最近的采样点进行计算。这意味着半步积分的采样点可以介于两者之间。接近最小值，这会导致$O(e·dt^2)$($e$除了大小的截断误差之外，仅来自不同时间位置的计算半径最小值之间的偏心度差$O(dt^4)$.

如果您使用更好的真实近日点近似值，这一切都会消失。（或者可以通过比较同一时间点的位置来减少，即使它们不是两个系列的局部最小值。）

我构建了一个类似的情况。混叠效应的严重程度很大程度上取决于轨道的参数。碰巧我尝试的第一个参数产生了问题中的调制振荡，其他参数产生了更平滑的图形。

首先求一个稍微偏心的椭圆的定步解：

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

然后绘制局部最小值

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

现在半径的最小值在半径函数的导数的根处，即$r\dot r=x\dot x+y\dot y$改变它的符号。使用线性近似来确定样本点之间的根位置更接近，然后使用一个 RK4 步长得到这个中间时间的解，总体精度相同

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

这将返回相当减少的情节

请注意，演化过程中差异的大小从$10^{-6}$在第一个情节中$10^{-9}$这里。