Truncation errors in wave equations

Linear wave equation in 1D

The standard, linear wave equation in 1D for a function \(u(x,t)\) reads

(1)\[ \frac{\partial^2 u}{\partial t^2} = c^2\frac{\partial^2 u}{\partial x^2} + f(x,t),\quad x\in (0, L),\ t\in (0,T],\]

where \(c\) is the constant wave velocity of the physical medium \([0,L]\). The equation can also be more compactly written as

(2)\[ u_{tt} = c^2u_{xx} + f,\quad x\in (0, L),\ t\in (0,T],\]

Centered, second-order finite differences are a natural choice for discretizing the derivatives, leading to

(3)\[ [D_t D_t u = c^2 D_xD_x u + f]^n_i\]\[ {\thinspace .}\]

Inserting the exact solution \({u_{\small\mbox{e}}}(x,t)\) in (3) makes this function fulfill the equation if we add the term \(R\):

(4)\[ [D_t D_t {u_{\small\mbox{e}}} = c^2 D_xD_x {u_{\small\mbox{e}}} + f + R]^n_i\]

Our purpose is to calculate the truncation error \(R\). From (1.15)-(1.16) we have that

\[[D_t D_t{u_{\small\mbox{e}}}]_i^n = {u_{\small\mbox{e}, tt}}(x_i,t_n) + \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 + {\mathcal{O}(\Delta t^4)},\]

when we use a notation taking into account that \({u_{\small\mbox{e}}}\) is a function of two variables and that derivatives must be partial derivatives. The notation \({u_{\small\mbox{e}, tt}}\) means \(\partial^2{u_{\small\mbox{e}}} /\partial t^2\).

The same formula may also be applied to the \(x\)-derivative term:

\[[D_xD_x{u_{\small\mbox{e}}}]_i^n = {u_{\small\mbox{e}, xx}}(x_i,t_n) + \frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2 + {\mathcal{O}(\Delta x^4)},\]

Equation (6) now becomes

\[\begin{split}{u_{\small\mbox{e}, tt}} + \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 &= c^2{u_{\small\mbox{e}, xx}} + c^2\frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2 + f(x_i,t_n) + \\ & \quad {\mathcal{O}(\Delta t^4,\Delta x^4)} + R^n_i {\thinspace .}\end{split}\]

Because \({u_{\small\mbox{e}}}\) fulfills the partial differential equation (PDE) (2), the first, third, and fifth terms cancel out, and we are left with

(5)\[ R^n_i = \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 - c^2\frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2 + {\mathcal{O}(\Delta t^4,\Delta x^4)},\]

showing that the scheme (3) is of second order in the time and space mesh spacing.

Finding correction terms

Can we add correction terms to the PDE and increase the order of \(R^n_i\) in (5)? The starting point is

(6)\[ [D_t D_t {u_{\small\mbox{e}}} = c^2 D_xD_x {u_{\small\mbox{e}}} + f + C + R]^n_i\]

From the previous analysis we simply get (5) again, but now with \(C\):

(7)\[ R^n_i + C_i^n = \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 - c^2\frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2 + {\mathcal{O}(\Delta t^4,\Delta x^4)}{\thinspace .}\]

The idea is to let \(C_i^n\) cancel the \(\Delta t^2\) and \(\Delta x^2\) terms to make \(R^n_i = {\mathcal{O}(\Delta t^4,\Delta x^4)}\):

\[C_i^n = \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 - c^2\frac{1}{12}{u_{\small\mbox{e}, xxxx}}(x_i,t_n)\Delta x^2{\thinspace .}\]

Essentially, it means that we add a new term

\[C = \frac{1}{12}\left( u_{tttt}\Delta t^2 - c^2u_{xxxx}\Delta x^2\right),\]

to the right-hand side of the PDE. We must either discretize these 4th-order derivatives directly or rewrite them in terms of lower-order derivatives with the aid of the PDE. The latter approach is more feasible. From the PDE we have that

\[\frac{\partial^2}{\partial t^2} = c^2\frac{\partial^2}{\partial x^2},\]

so

\[u_{tttt} = c^2u_{xxtt},\quad u_{xxxx} = c^{-2}u_{ttxx}{\thinspace .}\]

Assuming \(u\) is smooth enough that \(u_{xxtt}=u_{ttxx}\), these relations lead to

\[C = \frac{1}{12}((c^2\Delta t^2 - \Delta x^2)u_{xx})_{tt}{\thinspace .}\]

A natural discretization is

\[C^n_i = \frac{1}{12}((c^2\Delta t^2 - \Delta x^2) [D_xD_xD_tD_t u]^n_i{\thinspace .}\]

Writing out \([D_xD_xD_tD_t u]^n_i\) as \([D_xD_x (D_tD_t u)]^n_i\) gives

\[\begin{split}\frac{1}{\Delta t^2}\biggl( &\frac{u^{n+1}_{i+1} - 2u^{n}_{i+1} + u^{n-1}_{i+1}}{\Delta x^2} -2\\ &\frac{u^{n+1}_{i} - 2u^{n}_{i} + u^{n-1}_{i}}{\Delta x^2} + &\frac{u^{n+1}_{i-1} - 2u^{n}_{i-1} + u^{n-1}_{i-1}}{\Delta x^2} \biggr)\end{split}\]

Now the unknown values \(u^{n+1}_{i+1}\), \(u^{n+1}_{i}\), and \(u^{n+1}_{i-1}\) are coupled, and we must solve a tridiagonal system to find them. This is in principle straightforward, but it results in an implicit finite difference schemes, while we had a convenient explicit scheme without the correction terms.

Extension to variable coefficients (2)

Now we address the variable coefficient version of the linear 1D wave equation,

\[\frac{\partial^2 u}{\partial t^2} = \frac{\partial}{\partial x} \left( \lambda(x)\frac{\partial u}{\partial x}\right),\]

or written more compactly as

(8)\[ u_{tt} = (\lambda u_x)_x{\thinspace .}\]

The discrete counterpart to this equation, using arithmetic mean for \(\lambda\) and centered differences, reads

(9)\[ [D_t D_t u = D_x \overline{\lambda}^{x}D_x u]^n_i{\thinspace .}\]

The truncation error is the residual \(R\) in the equation

(10)\[ [D_t D_t {u_{\small\mbox{e}}} = D_x \overline{\lambda}^{x}D_x {u_{\small\mbox{e}}} + R]^n_i{\thinspace .}\]

The difficulty in the present is how to compute the truncation error of the term \([D_x \overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_i\).

We start by writing out the outer operator:

(11)\[ [D_x \overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_i = \frac{1}{\Delta x}\left( [\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_{i+\frac{1}{2}} - [\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_{i-\frac{1}{2}} \right).\]

With the aid of (1.3)-(1.4) and (1.19)-(1.20) we have

\[\begin{split}\lbrack D_x {u_{\small\mbox{e}}} \rbrack^n_{i+\frac{1}{2}} & = {u_{\small\mbox{e}, x}}(x_{i+\frac{1}{2}},t_n) + \frac{1}{24}{u_{\small\mbox{e}, xxx}}(x_{i+\frac{1}{2}},t_n)\Delta x^2 + {\mathcal{O}(\Delta x^4)},\\ \lbrack\overline{\lambda}^{x}\rbrack_{i+\frac{1}{2}} &= \lambda(x_{i+\frac{1}{2}}) + \frac{1}{8}\lambda''(x_{i+\frac{1}{2}})\Delta x^2 + {\mathcal{O}(\Delta x^4)},\\ [\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_{i+\frac{1}{2}} &= (\lambda(x_{i+\frac{1}{2}}) + \frac{1}{8}\lambda''(x_{i+\frac{1}{2}})\Delta x^2 + {\mathcal{O}(\Delta x^4)})\times\\ &\quad ({u_{\small\mbox{e}, x}}(x_{i+\frac{1}{2}},t_n) + \frac{1}{24}{u_{\small\mbox{e}, xxx}}(x_{i+\frac{1}{2}},t_n)\Delta x^2 + {\mathcal{O}(\Delta x^4)})\\ &= \lambda(x_{i+\frac{1}{2}}){u_{\small\mbox{e}, x}}(x_{i+\frac{1}{2}},t_n) + \lambda(x_{i+\frac{1}{2}}) \frac{1}{24}{u_{\small\mbox{e}, xxx}}(x_{i+\frac{1}{2}},t_n)\Delta x^2 + \\ &\quad {u_{\small\mbox{e}, x}}(x_{i+\frac{1}{2}}) \frac{1}{8}\lambda''(x_{i+\frac{1}{2}})\Delta x^2 +{\mathcal{O}(\Delta x^4)}\\ &= [\lambda {u_{\small\mbox{e}, x}}]^n_{i+\frac{1}{2}} + G^n_{i+\frac{1}{2}}\Delta x^2 +{\mathcal{O}(\Delta x^4)},\end{split}\]

where we have introduced the short form

\[G^n_{i+\frac{1}{2}} = (\frac{1}{24}{u_{\small\mbox{e}, xxx}}(x_{i+\frac{1}{2}},t_n)\lambda((x_{i+\frac{1}{2}}) + {u_{\small\mbox{e}, x}}(x_{i+\frac{1}{2}},t_n) \frac{1}{8}\lambda''(x_{i+\frac{1}{2}}))\Delta x^2{\thinspace .}\]

Similarly, we find that

\[\lbrack\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}\rbrack^n_{i-\frac{1}{2}} = [\lambda {u_{\small\mbox{e}, x}}]^n_{i-\frac{1}{2}} + G^n_{i-\frac{1}{2}}\Delta x^2 +{\mathcal{O}(\Delta x^4)}{\thinspace .}\]

Inserting these expressions in the outer operator (11) results in

\[\begin{split}\lbrack D_x \overline{\lambda}^{x}D_x {u_{\small\mbox{e}}} \rbrack^n_i &= \frac{1}{\Delta x}( [\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_{i+\frac{1}{2}} - [\overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_{i-\frac{1}{2}} )\\ &= \frac{1}{\Delta x}( [\lambda {u_{\small\mbox{e}, x}}]^n_{i+\frac{1}{2}} + G^n_{i+\frac{1}{2}}\Delta x^2 - [\lambda {u_{\small\mbox{e}, x}}]^n_{i-\frac{1}{2}} - G^n_{i-\frac{1}{2}}\Delta x^2 + {\mathcal{O}(\Delta x^4)} )\\ &= [D_x \lambda {u_{\small\mbox{e}, x}}]^n_i + [D_x G]^n_i\Delta x^2 + {\mathcal{O}(\Delta x^4)}{\thinspace .}\end{split}\]

The reason for \({\mathcal{O}(\Delta x^4)}\) in the remainder is that there are coefficients in front of this term, say \(H\Delta x^4\), and the subtraction and division by \(\Delta x\) results in \([D_x H]^n_i\Delta x^4\).

We can now use (1.3)-(1.4) to express the \(D_x\) operator in \([D_x \lambda {u_{\small\mbox{e}, x}}]^n_i\) as a derivative and a truncation error:

\[[D_x \lambda {u_{\small\mbox{e}, x}}]^n_i = \frac{\partial}{\partial x}\lambda(x_i){u_{\small\mbox{e}, x}}(x_i,t_n) + \frac{1}{24}(\lambda{u_{\small\mbox{e}, x}})_{xxx}(x_i,t_n)\Delta x^2 + {\mathcal{O}(\Delta x^4)}{\thinspace .}\]

Expressions like \([D_x G]^n_i\Delta x^2\) can be treated in an identical way,

\[[D_x G]^n_i\Delta x^2 = G_x(x_i,t_n)\Delta x^2 + \frac{1}{24}G_{xxx}(x_i,t_n)\Delta x^4 + {\mathcal{O}(\Delta x^4)}{\thinspace .}\]

There will be a number of terms with the \(\Delta x^2\) factor. We lump these now into \({\mathcal{O}(\Delta x^2)}\). The result of the truncation error analysis of the spatial derivative is therefore summarized as

\[[D_x \overline{\lambda}^{x}D_x {u_{\small\mbox{e}}}]^n_i = \frac{\partial}{\partial x} \lambda(x_i){u_{\small\mbox{e}, x}}(x_i,t_n) + {\mathcal{O}(\Delta x^2)}{\thinspace .}\]

After having treated the \([D_tD_t{u_{\small\mbox{e}}}]^n_i\) term as well, we achieve

\[R^n_i = {\mathcal{O}(\Delta x^2)} + \frac{1}{12}{u_{\small\mbox{e}, tttt}}(x_i,t_n)\Delta t^2 {\thinspace .}\]

The main conclusion is that the scheme is of second-order in time and space also in this variable coefficient case. The key ingredients for second order are the centered differences and the arithmetic mean for \(\lambda\): all those building blocks feature second-order accuracy.

1D wave equation on a staggered mesh

Linear wave equation in 2D/3D

The two-dimensional extension of (1) takes the form

(12)\[ \frac{\partial^2 u}{\partial t^2} = c^2\left(\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) + f(x,y,t),\quad (x,y)\in (0, L)\times (0,H),\ t\in (0,T],\]

where now \(c(x,y)\) is the constant wave velocity of the physical medium \([0,L]\times [0,H]\). In the compact notation, the PDE (12) can be written

(13)\[ u_{tt} = c^2(u_{xx} + u_{yy}) + f(x,y,t),\quad (x,y)\in (0, L)\times (0,H), \ t\in (0,T],\]

in 2D, while the 3D version reads

(14)\[ u_{tt} = c^2(u_{xx} + u_{yy} + u_{zz}) + f(x,y,z,t),\]

for \((x,y,z)\in (0, L)\times (0,H)\times (0,B)\) and \(t\in (0,T]\).

Approximating the second-order derivatives by the standard formulas (1.15)-(1.16) yields the scheme

\[[D_t D_t u = c^2(D_xD_x u + D_yD_y u) + f]^n_{i,j,k} {\thinspace .}\]

The truncation error is found from

\[[D_t D_t {u_{\small\mbox{e}}} = c^2(D_xD_x {u_{\small\mbox{e}}} + D_yD_y {u_{\small\mbox{e}}}) + f + R]^n {\thinspace .}\]

The calculations from the 1D case can be repeated to the terms in the \(y\) and \(z\) directions. Collecting terms that fulfill the PDE, we end up with

\[R^n_{i,j,k} = [\frac{1}{12}{u_{\small\mbox{e}, tttt}}\Delta t^2 - c^2\frac{1}{12}\left( {u_{\small\mbox{e}, xxxx}}\Delta x^2 + {u_{\small\mbox{e}, yyyy}}\Delta x^2 + {u_{\small\mbox{e}, zzzz}}\Delta z^2\right)]^n_{i,j,k} +\]

\[\quad {\mathcal{O}(\Delta t^4,\Delta x^4,\Delta y^4,\Delta z^4)}\nonumber {\thinspace .}\]