# Why is the central difference method dispersing my solution?

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

• Interesting. If you increase the number of points to n=10000, the solutions seem to be very similar to each other until about t=28. Then the Forward Euler solution becomes different than the other two. I'm not sure, though, why the central difference scheme seems to diverge to two different solutions. They seem to approximate the second-order one well until about t=46 for n=10000. I suggest plotting everything in the same figure and add a legend. Sep 13 at 21:11
• Is it possible that there is an actual bifurcation point there for the underlying differential equation? Sep 14 at 0:40
• Here is a GIF of the $n$ dependence on the congruence of solutions for different methods, I added Runge-Kutta:<"Click to enlarge."> Sep 14 at 4:33
• I haven't looked carefully at the parameters you're using, but in certain regimes the problem you are studying is chaotic, so you should expect numerical solutions to diverge strongly after some time. See e.g. galileoandeinstein.phys.virginia.edu/7010/… Sep 14 at 5:47
• Code nitpick: Your construction of t and dt is strange. To get n segments with the correct time step use t, dt = np.linspace(tmin, tmax, n+1, retstep=True). The present mismatch of dt and the time step in t could also introduce subtle errors in the result. /// In the multi-step method you use a first-order method to compute the first step in a second-order method. This initial error will dominate all the later second-order truncation errors, reducing the effective order to one. Sep 14 at 7:27

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)$$.(1) 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]