I am currently studying the effects of group velocity on the finite difference solution of the wave equation. Most of what I learned is from this source. I understand that high frequency components in the initial and boundary conditions can cause numerical dispersion and anisotropy in the 2D case. While the dispersion problem can be compensated to some degree with an appropriate temporal and spatial step size when dealing with a one-dimensional wave equation, compensating for dispersion (and anisotropy) in the solution of the 2D wave equation seems very hard. I have a question on how to do that.

Let's consider the 2D wave equation:

$$ u_{tt}=c^2 \Big ( u_{xx} + u_{yy} \Big ) \tag 1 $$

The first thing that comes to my mind is to improve the accuracy of the numerical solution. The temporal accuracy can probably be improved by using a fourth-order Runge-Kuta algorithm instead of a standard three-point central difference approximation for the second-order time derivative. However, I have never seen this in practice. So, my question is, can this be done and does this help in eliminating numerical dispersion and anisotropy?

  • 1
    $\begingroup$ Welcome to Scicomp.SE! As part of the guidelines for this site, you should ask (only) one clear question per post. $\endgroup$ Commented Feb 5 at 3:36
  • $\begingroup$ Thank you for the notice, I didn't know. I will keep that in mind next time. $\endgroup$
    – Amilox Lex
    Commented Feb 5 at 13:12
  • 1
    $\begingroup$ You can edit your question to make it follow the guidelines. I think you're more likely to get good answers if you do so. $\endgroup$ Commented Feb 5 at 14:39

1 Answer 1


Question 1

The temporal accuracy can probably be improved by using a fourth order Runge-Kuta algorithm instead of a standard three-point central difference approximation for the second order time derivative. However, I have never seen this in practice. So, my first question is, can this be done and does this help in eliminating numerical dispersion and anisotropy?

I think the anisotropy is mostly a function of how well you have discretised things in space, i.e. the spatial discretisation step size and to what degree your discretisation of $\Delta$ is rotationally invariant. I am not sure about dispersion. But I have provided a solution that you may use as a reference in question 4.

Question 2

The spatial accuracy can be improved by using a five-point central difference approximation instead of a three-point one. However, this requires knowing additional boundary conditions because instead of one ghost node, there are two if we are dealing with Neumann type boundary conditions. Usually, those additional boundary conditions are not known. My second question is, can we use a three-point spatial central difference approximation to solve the PDE at the boundaries, and then use a five-point central difference approximation to solve the PDE in the interior of the numerical grid?

You can use a three point stencil near the boundaries, but I think that your interpretation that you need higher order boundary conditions is wrong. E.g. assume that $x_{n+1}$ is a point on the boundary and $x_n$ is on the inside. Then with a three point stencil you would have $d_{xx} u(x_n) \approx \frac{u(x_{n-1}) - 2u(x_n) + u(x_{n+1})}{h^2}$ and you can figure out $u(x_{n+1})$ from e.g. a reflecting boundary condition $d_x u(x_{n+1}) = 0$, then $\frac{u(x_{n+1}) - u(x_n)}{h} = 0$ and you can substitute $u(x_{n+1}) = u(x_n)$ in order to get the reduced stencil: $d_{xx} u(x_n) \approx \frac{u(x_{n-1}) - u(x_n)}{h^2}$. I am pretty sure you can do the same thing with a five point stencil by just using first derivatives (using second derivatives would in fact be wrong afaik): $$d_{xx} u(x_n) \approx \frac{1}{12h^2}(-u(x_{n-2}) + 16 u(x_{n-1}) -30u(x_n) + 16 u(x_{n+1}) - u(x_{n+2})).$$ Say you still have the reflecting boundary conditions $d_x u(x_n) = 0$, but now you could try to make two realizations of those: $$d_x u(x_n) \approx \frac{u(x_{n+1})-u(x_{n})}{h} = 0, \quad d_x u(x_n) \approx \frac{u(x_{n+2}) - u(x_{n})}{2h} = 0.$$ From here you can take both to be the same as $u(x_n)$ and then you get the reduced stencil: $$d_{xx} u(x_n) \approx \frac{1}{12h^2}(-u(x_{n-2}) + 16 u(x_{n-1}) -15u(x_n)).$$

I don't like that the second term is larger in magnitude here though, so I would probably just stick to the three point stencil. In fact in Wesseling's book you can also find something to that effect (it was in the context of multigrid though - that interpolation operators ought to be modified to lower order ones near boundaries).

Question 3

My third question is, is it ok to use a three point spatial central difference approximation for one variable, and a five point spatial central difference approximation for the second variable?

Yes it is technically ok, you would just have a higher order of consistency in one direction. I have no idea why you would want to do so, however, unless you know that you need a higher consistency in some direction. At this point you may want to look at hp-FEM especially if you want a spatially adaptive irregular mesh instead of just along one direction.

Question 4

And finally, what are some other methods we can use to compensate for numerical anisotropy and dispersion when solving the 2D wave equation with the finite difference method?

If you want to compute the solution without discretising time you can use the following approach:

$$\partial^2_t u = c^2\Delta u \implies \begin{cases} \partial_t u = v, \\ \partial_t v = c^2 \Delta u. \end{cases}$$

With this substitution $v = \partial_t u$ you can rewrite the problem as a vector-valued problem with a single $\partial_t$:

$$\partial_t \begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} 0 & id \\ c^2\Delta & 0 \end{bmatrix} \begin{bmatrix} u \\ v \end{bmatrix}.$$

To solve a vector problem $d_t w = Mw$ (equivalently $d_t w - Mw = 0$) you want to bring it to a form $d_t(Qw) = 0$ so that you may integrate both sides and get $Q(t)w(t) - Q(0)w(0) = 0$ and thus $w(t) = (Q(t))^{-1}Q(0)w(0)$. So it pays off to study $d_t (Qw) = Q w_t + Q_t w = 0$, specifically you can multiply $d_t w - Mw = 0$ by $Q$ in order to get to that form: $Q d_t w - QM w = 0$. Then it follows that you need $Q_t = -QM$, this holds if $Q(t) = \exp(-tM)$ (note that $\exp(-tM)$ and $M$ commute). With this you get that $w(t) = (\exp(-tM))^{-1}\exp(-0M)w(0) = \exp(tM)w(0)$. You can now apply this to the original problem to get:

$$\begin{bmatrix} u(t) \\ v(t) \end{bmatrix} = \left(\exp\left(t\begin{bmatrix} 0 & id \\ c^2\Delta & 0 \end{bmatrix}\right) \begin{bmatrix} u(\cdot) \\ v(\cdot) \end{bmatrix}\right)(0).$$

You can now use the series expansion of the exponent:

$$\exp(tM) = I + tM + \frac{t^2M^2}{2} + \frac{t^3M^3}{3!} + \ldots.$$

With a bit of linear algebra you can arrive at the following:

$$exp\left(tM\right) = \begin{bmatrix} id & 0 \\ 0 & id \end{bmatrix} + t\begin{bmatrix} 0 & id \\ c^2\Delta & 0 \end{bmatrix} + \frac{t^2}{2} \begin{bmatrix} c^2\Delta & 0 \\ 0 & c^2\Delta\end{bmatrix} + \frac{t^3}{3!}\begin{bmatrix} 0 & c^2\Delta \\ (c^2\Delta)^2 & 0 \end{bmatrix} + \ldots.$$

Consequently for $u$ you get the following solution: \begin{align} u(t) &= \left((t\sqrt{c^2\Delta})^0 + \frac{(t\sqrt{c^2\Delta})^2}{2} + \frac{(t\sqrt{c^2\Delta})^4}{4!}+\ldots\right)u(0) \\ &+\frac{1}{\sqrt{c^2\Delta}}\left((t\sqrt{c^2\Delta}) + \frac{(t\sqrt{c^2\Delta})^3}{3!} + \frac{(t\sqrt{c^2\Delta})^5}{5!}+\ldots\right)\partial_t u(0). \end{align}

You can use the fact that $\exp(ix) = cos x + i\sin x$ to rewrite these series as a cosine and a sine:

\begin{align} u(t) &= \left((it\sqrt{c^2\Delta})^0 - \frac{(it\sqrt{c^2\Delta})^2}{2} + \frac{(it\sqrt{c^2\Delta})^4}{4!}-\ldots\right)u(0) \\ &+\frac{1}{i\sqrt{c^2\Delta}}\left((it\sqrt{c^2\Delta}) -\frac{(it\sqrt{c^2\Delta})^3}{3!} + \frac{(it\sqrt{c^2\Delta})^5}{5!}-\ldots\right)\partial_t u(0) \\ &= \cos(it\sqrt{c^2\Delta})u(0) +(i\sqrt{c^2\Delta})^{-}\sin(it\sqrt{c^2\Delta})\partial_t u(0) \\ &= \cos(t\sqrt{-c^2\Delta})u(0) +(\sqrt{-c^2\Delta})^{-}\sin(t\sqrt{-c^2\Delta})\partial_t u(0). \end{align}

If you discretise $-\Delta$ with $L$ and $c^2$ with a diagonal matrix $C^2$ then your problem boils down to computing $$u(t) = \cos(t\sqrt{C^2L}) u(0) + (\sqrt{C^2L})^{-}\sin(t\sqrt{C^2L})d_t u(0),$$ where $A^{-}$ is a generalised inverse (e.g. the Moore-Penrose pseudoinverse). If you know the eigendecomposition $C^2L = PDP^{-1}$ then you can compute the above matrix functions by applying them to the eigenvalues: $$f(C^2L) = P f(D) P^{-1} = P\begin{bmatrix} f(d_1) & & \\ & \ddots & \\ && f(d_n)\end{bmatrix}P^{-1}.$$ Notably, if you have a discretisation on a grid then you know the eigendecomposition of $L$ (e.g. diagonalization by using DST/DCT/DFT depending on the BCs), and if $C^2$ is $c^2 I$, then you can compute the solution directly. In fact you could directly use the spectral discretisation of $-\Delta$ in this case which would yield in the spatial domain a full stencil, but it's diagonal in the DST/DCT/DFT domain.

If $c$ is space-variant but not zero anywhere then you can perform the following trick: $$\partial^2_t u = c^2\Delta u, \,\, u = cw \implies \partial^2_t w = c \Delta cw \implies -c\Delta c \approx CLC.$$ Now at least $CLC$ is symmetric unlike $C^2L$ (and it has the same eigenvalues). I can't think of a way to easily compute the eigenvalues and eigenvectors knowing the ones of $L$ however. But at least if your matrix is smaller you can compute those with a Lanczos iteration considering the matrix is symmetric. You could also try methods that estimate matrix functions without directly computing the eigendecomposition, for example you can look into Limited-memory polynomial methods for large-scale matrix functions.


Not the answer you're looking for? Browse other questions tagged or ask your own question.