Why is the central difference method dispersing my solution?

Question

I am solving numerically the ODE $\ddot x(t)=-c\dot x(t) -\sin(x(t))+F\cdot \cos(\omega t), \;\dot x(0)=x(0)=0$ for $t\in [0,20\pi]$ on an $N=2000$ dimensional grid. I am working on Python, and I replaced the time derivatives by the finite difference operators \begin{align} \dot x(t)&=\frac{x(t+dt)-x(t)}{dt}, & \ddot x(t) &=\frac{x(t+dt)-2x(t)+x(t-dt)}{dt^2} \end{align} Thus, the discretized ODE can be solved by the extrapolation scheme $$x(t+dt) =\frac 1{1+c\cdot dt}\Bigr( x(t)\bigr(2+c\cdot dt\bigr) +x(t-dt) + dt^2\bigr( -\sin(x(t))+ F\cdot \cos(\omega t)\bigr) \Bigr) $$ I performed this scheme and the result was great. To compare, I also reduced the 2nd order ODE to a system of two 1st order ODE's \begin{align} \begin{cases} \dot x=y\\ \dot y(t) =-cy(t)-\sin x(t)+Fcos(\omega t) \end{cases} \end{align} with initial values $(x(0),y(0))=(0,0)$. Solving this system with a Forward-Euler scheme, yields the a solution that starts similar to the first scheme, but is not quite the same. Note that using the Forward-Euler scheme is the same as solving for $x(t+dt)$ after replacing the derivatives by the forward difference operator $f'(t)=\frac1 {dt}(f(t+dt)-f(t)).$ However, if we rather replace the derivatives with the central difference operator $$f(t)=\frac{f(t+dt)-f(t-dt)}{2\cdot dt}, $$ we obtain the following extrapolation scheme: \begin{align} \begin{cases} x(t+dt) = x(t-dt)+2\cdot dt \cdot y(t)\\ y(t+dt) = y(t-dt) + 2\cdot dt \bigr( - c y(t) -\sin x(t)+ F\cdot \cos(\omega t) \bigr) \end{cases} \end{align} Performing this scheme with values $c=0.05, \; \omega = 0.7, \; F = 0.4$ yields a horrible solution. I do not know what is going on, here are the plots for the solutions obtained in each scheme.

Also, here is the code I wrote

import matplotlib
import matplotlib.pyplot as plt


tmax = 20 * np.pi
tmin = 0
n = 1000
dt = (tmax-tmin)/(n-1)
t = np.linspace(0,(n+1)*dt,n)
c = 0.05
w = 0.7
F = 0.4
x = np.zeros(n)
y = np.zeros(n)
for i in range(1,n-1):
    """
    First method
    """
    x[i+1] = (x[i] * (2+ c*dt) - x[i-1] + (dt ** 2)  * (-np.sin(x[i]) \
        + F * np.cos(w*t[i])) )/(1+c*dt)
plt.scatter(t,x,s = 0.5)
plt.title('Second order ODE scheme')
plt.xlabel('t')
plt.ylabel('x')
plt.show()
plt.close()    
for i in range(1,n-1): 
    """
    Second method
    """      
    x[i+1] = x[i] + dt * y[i]
    y[i+1] = y[i] + dt * ( - c * y[i] - np.sin(x[i])\
        + F * np.cos(w * t[i]))
plt.scatter(t,x,s = 0.5)
plt.xlabel('t')
plt.ylabel('x')
plt.title('First order ODE system (Forward)')
plt.show()
plt.close()    
for i in range(1,n-1):
    """
    Third method
    """
    x[i+1] = x[i-1] +   2*dt * y[i]
    y[i+1] = y[i-1] +   2*dt * ( - c * y[i] - np.sin(x[i]) \
        + F * np.cos(w * t[i]))           
plt.scatter(t,x,s = 0.5)
plt.xlabel('t')
plt.ylabel('x')
plt.title('First order ODE system (Central)')
plt.show()
plt.close()

Answer 1

I'll write the equation short as $$\ddot x(t)+c\dot x(t)=a(t,x(t))$$ to separate the "easy" linear parts from the non-linear and forcing terms.

On the first method

The claimed order of the method is two, the implemented order is one. This is due to the implementation of the first derivative as a one-sided difference quotient. Changing that to the central difference quotient makes the whole method second order, $$ \frac{x(t+dt)-2x(t)+x(t-dt)}{dt^2}+c\frac{x(t+dt)-x(t-dt)}{2dt}=a(t,x(t)) $$ gives

for i in range(1,n-1):
    x[i+1] = ( 2*x[i] - x[i-1]*(1-0.5*c*dt) +  dt**2 * a(t[i],x[i]) )/(1+0.5*c*dt)

On consistent initial data

Multi-step methods need a bootstrap to fill an initial segment of values of the same error order or a little more precise. For the correct implementation of second order multi-step methods, this initial segment consists of he first two states. The required bootstrap step should at least be second order accurate, that is, with an error of $O(dt^3)$.⁽¹⁾ Else this initial error will dominate all the later second-order truncation errors, potentially reducing the effective order to one.

_{(1) : using the convention that for methods with global error order $p$ the local truncation errors, not converted to unit-step errors, are of size $O((dt)^{p+1})$}

So with the stated initial conditions $x(0)=\dot x(0)=0$ one gets that the first step should be $$x(dt)=\frac12\ddot x(0)\,dt^2+O(dt^3)=\frac12F\,dt^2+O(dt^3).$$ Alternatively, to get the implemented behavior $x(0)=x(dt)=0$, one can force $x(dt)=0+O(dt^3)$ by setting $\dot x(0)+\frac12\ddot x(0)dt=0$, giving $$y(0)=\dot x(0)=-\frac{F\,dt}{2-c\,dt}$$

One of these variants should be consistently followed, for this answer it will be the second one. To compare the results for different step sizes, one should follow the first convention.

On the second method

It makes no sense to compare methods of different orders with the same step size for a chaotic dynamical process like the forced and dampened pendulum. So either use a good number of Euler sub-steps with a smaller step size for every step of the second order methods, or use an explicit second order method like Heun's method

for i in range(0,n-1): 
    kx1, ky1 = y[i],  - c * y[i] + a(t[i],x[i])
    x[i+1], y[i+1] = x[i]+dt*kx1, y[i]+dt*ky1
    kx2, ky2 = y[i+1],  - c * y[i+1] + a(t[i+1],x[i+1])
    x[i+1], y[i+1] = x[i]+dt*0.5*(kx1+kx2), y[i]+dt*0.5*(ky1+ky2)

Note that the iteration has to start with the initial values at i=0.

On the third method

The central Euler method, aka Nyström method, is weakly stable, that is, its stability region is the segment $[-i,i]$ on the imaginary axis. In its error formula it has an oscillating term, $e_k=c_k+(-1)^kd_k$, where the sequences $c_k$ and $d_k$ are relatively smooth and both potentially exponentially growing (see CE1,CE2). $d_k$ is mainly proportional to the error in the first step, so again it is important to keep that small.

The failure to keep the first-step error small leads to the observed (early) dispersal of the numerical solution. As you keep the first two steps as zero from the start, you get a first order error $O(dt^2)$ in the x and a zero-order error $O(dt)$ in the velocity y component. These get exponentially magnified in time, resulting in the visual separation of the odd and even indexed sub-sequences of the numerical solution sequence.

Re-using the Heun values from the (modified) second method is sufficient to keep these two sub-sequences together over the integration interval. It is to be expected that over longer intervals the separation still occurs, as the initial error is not zero and floating-point errors also accumulate and magnify in the coefficients of the alternating error terms.

Visually, the methods now give the same result.

The values stay sufficiently close together for this range, but a certain deviation between the methods is visible

   x[:4]                                         x[-4:]

[0.         0.         0.00157828 0.00471908] [-5.71998607 -5.66266182 -5.60624037 -5.55089183]
[0.         0.         0.00157984 0.00472373] [-5.73204895 -5.67600568 -5.62081879 -5.56665659]
[0.         0.         0.00158141 0.00472833] [-5.66988533 -5.61833293 -5.54752403 -5.50001843]