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)$. 935
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 numpy as np 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+ cdt) - x[i-1] + (dt 2) * (-np.sin(x[i]) \ + F * np.cos(wt[i])) )/(1+cdt) 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. 509