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} $$


$$ 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.

  • $\begingroup$ To quote the accepted answer for this related post scicomp.stackexchange.com/questions/20197/…: "High order methods are not necessarily always more accurate than low order methods; they simply have a more rapid convergence rate. " $\endgroup$
    – Charlie S
    Commented Nov 4, 2020 at 20:48
  • $\begingroup$ @CharlieS so similar to how a higher order Maclaurin series is more accurate close to the origin, but usually more rapidly diverges for larger values? $\endgroup$
    – fibonatic
    Commented Nov 5, 2020 at 2:08
  • $\begingroup$ The error implied in $O(Δx^4)$ is not necessarily smaller than $O(Δx^2)$ for all values of $Δx$. There is an implied constant in front of those Big O's. However, when $Δx$ is sufficiently small, you would expect the error to drop much faster with the higher order scheme. For example, what function is bigger: $x^4$ or $10,000 x^2$? Depends on $x$! $\endgroup$
    – Charlie S
    Commented Nov 5, 2020 at 4:40
  • $\begingroup$ @fibonatic Does your stability analysis use Fourier expansion in the spatial coordinate? Does not look like it is done there. With Fourier decomposition the stability threshold could be interpreted in terms of some critical mode. $\endgroup$ Commented Nov 5, 2020 at 5:50
  • $\begingroup$ The basic idea is that the higher the order becomes, the larger the extremal eigenvalues from the corresponding FD matrices get. And then, you need smaller timesteps to enlarge the stability region of your chosen integration method (here explicit Euler), so that all those eigenvalues lie in the stability region. So you lose performance, but gain accuracy (in terms of polynomial orders that are integrated exactly). $\endgroup$
    – davidhigh
    Commented Nov 5, 2020 at 8:18


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.