I was playing around with numerically solving the 1D wave equation with density and stiffness varying with position using central differencing methods and noticed that for certain discretization steps sizes $\Delta t$ and $\Delta x$ the solution could become unstable. Therefore, I thought maybe if I increased the order of the differencing method I could get away with a courser grid, but the opposite seems to be the case.
To make it easier to analyse I simplified the problem to constant material parameters, so the PDE considered can be reduced to
$$ \frac{\partial^2}{\partial t^2} u(t,x) = c^2 \frac{\partial^2}{\partial x^2} u(t,x), \tag{1} $$
with $c$ the wave speed and boundary conditions such that $u(t,x+L) = u(t,x)$, mainly because this was easiest to implement with both differencing methods. For discretization I used a uniform grid such that
$$ u^n_k = u(n\,\Delta t, k\,\Delta x). \tag{2} $$
For the integration with respect to time I used Verlet integration, so
$$ \frac{\partial^2}{\partial t^2} u(t,x) = \frac{u^{n+1}_k - 2\,u^n_k + u^{n-1}_k}{\Delta t^2} + O(\Delta t^2). \tag{3} $$
My initial approach used the same differencing method as $(3)$ for the spacial partial derivatives,
$$ \frac{\partial^2}{\partial x^2} u(t,x) = \frac{u^n_{k+1} - 2\,u^n_k + u^n_{k-1}}{\Delta x^2} + O(\Delta x^2). \tag{4} $$
Combining $(3)$ and $(4)$ therefore yields
$$ u^{n+1}_k = 2\,u^n_k - u^{n-1}_k + v^2 (u^n_{k+1} - 2\,u^n_k + u^n_{k-1}), \tag{5} $$
with
$$ v = \frac{c}{\Delta x / \Delta t}, $$
which can be seen as the wave speed normalized with respect to the discretization.
For $v>1$ the system becomes unstable, which also makes sense to me, since waves should propagate more than $\Delta x$ in one time step, but $(5)$ only considers neigboring grid positions.
For my second approach used a higher order differencing method for the spacial partial derivatives,
$$ \frac{\partial^2}{\partial x^2} u(t,x) = \frac{-u^n_{k+2} + 16\,u^n_{k+1} - 30\,u^n_k + 16\,u^n_{k-1} - u^n_{k-2}}{12\,\Delta x^2} + O(\Delta x^4). \tag{6} $$
Combining $(3)$ and $(6)$ therefore yields
$$ u^{n+1}_k = 2\,u^n_k - u^{n-1}_k + v^2 \frac{-u^n_{k+2} + 16\,u^n_{k+1} - 30\,u^n_k + 16\,u^n_{k-1} - u^n_{k-2}}{12}. \tag{7} $$
My initial guess would be that $(7)$ would be able to handle values for $v$ up to two, since it now sees two neigboring grid positions, so it can "see" waves coming from twice the distance. However, the opposite seems to be the case. Namely, now the system becomes unstable for $v > \sqrt{3/4} \approx 0.866$. And it is not obvious to me why $\sqrt{3/4}$, so I wonder if there is an intuitive explanation for this?
I determined the stability of each system by formulating it as $x[n+1] = A\,x[n]$, with
$$ x[n] = \begin{bmatrix}u^{n-1}_1 & \cdots & u^{n-1}_N & u^n_1 & \cdots & u^n_N\end{bmatrix}^\top $$
which is stable if $A$ has all it eigenvalues inside the unit disk.
