Total energy oscillations with Runge-Kutta 4th order method - how to avoid them?

Question

Consider a free ideal pendulum, which obeys the equations:

$$\frac{d\varphi}{dt}=p \\ \frac{dp}{dt}=-\omega_{0}^{2}\sin\varphi$$

I am applying two 4th order Runge-Kutta schemes: the usual explicit one, and a special implicit scheme I developed.

Both methods are based on the Simpson's formula for a general 1st order ODE:

$$\frac{dy}{dt}=f(y)$$

$$y_{m+1}=y_{m}+\frac{\Delta t}{6}\left[f(y_{m})+4f\left(y_{m+1/2}\right)+f\left(y_{m+1}\right)\right]$$

The difference is in how we estimate $y_{m+1/2}$ and $y_{m+1}$ on the rign hand side.

In the explicit method we proceed with the recurrence:

$$r_1=f(y_{m}) \\ r_2=f\left(y_{m}+\frac{\Delta t}{2}r_1\right) \\ r_3=f\left(y_{m}+\frac{\Delta t}{2}r_2\right) \\ r_4=f\left(y_{m}+\Delta t r_3\right)$$

$$y_{m+1}=y_{m}+\frac{\Delta t}{6}\left[r_1+2r_2+2r_3+r_4\right]$$

In the implicit method I used a special fixed point iteration procedure, which I won't describe here.

---

Calculation parameters are as follows:

$$T=\frac{2\pi}{\omega_{0}}=2,\quad\Delta t=4\cdot10^{-3},\quad\varphi_{0}=5,\quad p_{0}=10$$

The results for the total energy are:

Where $n=4$ is the number of fixed point iterations I use for the implicit method.

---

While both methods are pretty accurate, there's an energy drift in the explicit method and spurious oscillations in both methods. The oscillations don't seem to be exactly equal to the pendulum period, but close. Maybe because I chose a large angle, the actual period differs from the theoretical value.

In any case, I would like to know how to modify the method to avoid these oscillations, as they can get quite annoying in more complicated, especially chaotic problems.

---

Edit: Changed the plot so it displays the energy change instead of absolute energy. Also improved the method.

Answer 1

The following is not an answer to the question, as that asks for an error free numerical method, which is impossible.

In the discussion it has been of interest how I would implement the implicit collocation method of order 4 based on the same points as RK4. This turned out to be relatively short,

from numba import jit

@jit()
def rk4_step(f,t,y,h):
    k1=f(t,y)
    k2=f(t+h/2,y+h*k1/2)
    k3=f(t+h/2,y+h*k2/2)
    k2=(k2+k3)/2
    k3=f(t+h,y+h*k3)
    return y+h*(k1+4*k2+k3)/6,k1,k2,k3

@jit()
def coll4_step(f,t,y,h,ncorr):
    y3,k1,k2,k3 = rk4_step(f,t,y,h)
    for it in range(ncorr):
        y2 = y+h*(5*k1+8*k2-k3)/24
        k2 = f(t+h/2,y2)
        k3 = f(t+h,y3)
        y3 = y+h*(k1+4*k2+k3)/6
    return y3

@jit()
def odeint4(f,t,y0,ncorr):
    y = np.empty((len(t),)+y0.shape)
    y[0] = y0
    for k in range(1,len(t)):
        h=t[k]-t[k-1]
        y[k] = coll4_step(f,t[k-1],y[k-1],h,ncorr)
    return y

The explicit predictor method is standard RK4, so calling with ncorr=0 corrector steps gives the RK4 base case. The implicit method is the collocation method for the nodes $c_i=0,\frac12,1$. It proceeds by constructing a cubic polynomial $p(s)\approx y(t_n+sh)$ and satisfying $p'(c_i)=hk_i$ where $$ k_i=y'(t_n+c_ih)=f(t_n+c_ih,p(c_i))$. $$ This gives via interpolation $$ p'(s)=h(k_1(s-1)(2s-1)-4k_2s(s-1)+k_3s(2s-1)) $$ The polynomial can now be reconstructed via integration, which concludes the corrector step $$ p(s)=y_n+\frac{h}{6}(k_1(4s^3-9s^2+6s)-4k_2(2s^3-3s^2)+k_3(4s^3-3s^2) $$ The only values we need from that polynomial are the point in the middle $$ p(\tfrac12)=y_n+\frac{h}{24}(5 k_1+8k_2-k_3). $$ and the next point in the numerical trajectory $$ y_{n+1}=p(1)=y_n+\frac{h}{6}(k_1+4k_2+k_3), $$ which is the Simpson rule as expected.

Called with the task parameters

from numpy import sin, cos, pi

w0 = pi
u0 = np.array([5.,10.])
t = np.arange(0,20,4e-3)

@jit()
def sys(t,u): y,v = u; return np.array([v,-w0**2*sin(y)])

@jit()
def energy(u): y,v = u; return v**2/2+w0**2*(1-cos(y))

E0 = energy(u0)

%matplotlib inline
plt.figure(figsize=(8,5))
for k in range(4):
    print(f"num corrections = {k}", end="\r")
    u = odeint4(sys,t,u0,k)
    dE = energy(u.T)-E0
    plt.plot(t,dE,label=f"ncorr={k}")
    
plt.grid(); plt.legend(); plt.show()

results in the plot of the energy oscillations

As the 3/8 method in general has half the error of the Simpson method, this is repeated in the implicit methods based on them, as is readily visible.

To get the same reduced error with the Simpson scheme, one would need to divide the step size by $\sqrt[4]2=1.1892$, which increases the number of function evaluations by the same factor. This indicates that the graphs of error vs. function evaluations would look similar for both methods.