$$\newcommand{\half}{\frac{1}{2}} \newcommand{\halfi}{{1/2}} \newcommand{\tp}{\thinspace .} \newcommand{\uex}{{u_{\small\mbox{e}}}} \newcommand{\uexd}[1]{{u_{\small\mbox{e}, #1}}} \newcommand{\vex}{{v_{\small\mbox{e}}}} \newcommand{\Aex}{{A_{\small\mbox{e}}}} \newcommand{\E}[1]{\hbox{E}\lbrack #1 \rbrack} \newcommand{\Var}[1]{\hbox{Var}\lbrack #1 \rbrack} \newcommand{\xpoint}{\boldsymbol{x}} \newcommand{\normalvec}{\boldsymbol{n}} \newcommand{\Oof}[1]{\mathcal{O}(#1)} \newcommand{\x}{\boldsymbol{x}} \renewcommand{\u}{\boldsymbol{u}} \renewcommand{\v}{\boldsymbol{v}} \newcommand{\acc}{\boldsymbol{a}} \newcommand{\rpos}{\boldsymbol{r}} \newcommand{\e}{\boldsymbol{e}} \newcommand{\f}{\boldsymbol{f}} \newcommand{\F}{\boldsymbol{F}} \newcommand{\stress}{\boldsymbol{\sigma}} \newcommand{\I}{\boldsymbol{I}} \newcommand{\T}{\boldsymbol{T}} \newcommand{\q}{\boldsymbol{q}} \newcommand{\g}{\boldsymbol{g}} \newcommand{\dfc}{\alpha} % diffusion coefficient \newcommand{\ii}{\boldsymbol{i}} \newcommand{\jj}{\boldsymbol{j}} \newcommand{\kk}{\boldsymbol{k}} \newcommand{\ir}{\boldsymbol{i}_r} \newcommand{\ith}{\boldsymbol{i}_{\theta}} \newcommand{\Ix}{\mathcal{I}_x} \newcommand{\Iy}{\mathcal{I}_y} \newcommand{\Iz}{\mathcal{I}_z} \newcommand{\It}{\mathcal{I}_t} \newcommand{\setb}[1]{#1^0} % set begin \newcommand{\sete}[1]{#1^{-1}} % set end \newcommand{\setl}[1]{#1^-} \newcommand{\setr}[1]{#1^+} \newcommand{\seti}[1]{#1^i} \newcommand{\stepzero}{*} \newcommand{\stephalf}{***} \newcommand{\stepone}{**} \newcommand{\baspsi}{\psi} \newcommand{\dx}{\, \mathrm{d}x} \newcommand{\ds}{\, \mathrm{d}s} \newcommand{\Real}{\mathbb{R}} \newcommand{\Integer}{\mathbb{Z}}$$

# Wave equations

## Linear wave equation in 1D

The standard, linear wave equation in 1D for a function $$u(x,t)$$ reads $$$$\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], \tag{7.68}$$$$ where $$c$$ is the constant wave velocity of the physical medium in $$[0,L]$$. The equation can also be more compactly written as $$$$u_{tt} = c^2u_{xx} + f,\quad x\in (0, L),\ t\in (0,T], \tag{7.69}$$$$ Centered, second-order finite differences are a natural choice for discretizing the derivatives, leading to $$$$[D_t D_t u = c^2 D_xD_x u + f]^n_i \tag{7.70} \tp$$$$

Inserting the exact solution $$\uex(x,t)$$ in (7.70) makes this function fulfill the equation if we add the term $$R$$: $$$$[D_t D_t \uex = c^2 D_xD_x \uex + f + R]^n_i \tag{7.71}$$$$

Our purpose is to calculate the truncation error $$R$$. From (7.17)-(7.18) we have that $$[D_t D_t\uex]_i^n = \uexd{tt}(x_i,t_n) + \frac{1}{12}\uexd{tttt}(x_i,t_n)\Delta t^2 + \Oof{\Delta t^4},$$ when we use a notation taking into account that $$\uex$$ is a function of two variables and that derivatives must be partial derivatives. The notation $$\uexd{tt}$$ means $$\partial^2\uex /\partial t^2$$.

The same formula may also be applied to the $$x$$-derivative term: $$[D_xD_x\uex]_i^n = \uexd{xx}(x_i,t_n) + \frac{1}{12}\uexd{xxxx}(x_i,t_n)\Delta x^2 + \Oof{\Delta x^4},$$ Equation (7.71) now becomes \begin{align*} \uexd{tt} + \frac{1}{12}\uexd{tttt}(x_i,t_n)\Delta t^2 &= c^2\uexd{xx} + c^2\frac{1}{12}\uexd{xxxx}(x_i,t_n)\Delta x^2 + f(x_i,t_n) + \\ & \quad \Oof{\Delta t^4,\Delta x^4} + R^n_i \tp \end{align*} Because $$\uex$$ fulfills the partial differential equation (PDE) (7.69), the first, third, and fifth term cancel out, and we are left with $$$$R^n_i = \frac{1}{12}\uexd{tttt}(x_i,t_n)\Delta t^2 - c^2\frac{1}{12}\uexd{xxxx}(x_i,t_n)\Delta x^2 + \Oof{\Delta t^4,\Delta x^4}, \tag{7.72}$$$$ showing that the scheme (7.70) 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 (7.72)? The starting point is $$$$[D_t D_t \uex = c^2 D_xD_x \uex + f + C + R]^n_i \tag{7.73}$$$$ From the previous analysis we simply get (7.72) again, but now with $$C$$: $$$$R^n_i + C_i^n = \frac{1}{12}\uexd{tttt}(x_i,t_n)\Delta t^2 - c^2\frac{1}{12}\uexd{xxxx}(x_i,t_n)\Delta x^2 + \Oof{\Delta t^4,\Delta x^4}\tp \tag{7.74}$$$$ The idea is to let $$C_i^n$$ cancel the $$\Delta t^2$$ and $$\Delta x^2$$ terms to make $$R^n_i = \Oof{\Delta t^4,\Delta x^4}$$: $$C_i^n = \frac{1}{12}\uexd{tttt}(x_i,t_n)\Delta t^2 - c^2\frac{1}{12}\uexd{xxxx}(x_i,t_n)\Delta x^2\tp$$ 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 the operator equality $$\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}\tp$$ Assuming $$u$$ is smooth enough, so that $$u_{xxtt}=u_{ttxx}$$, these relations lead to $$C = \frac{1}{12}((c^2\Delta t^2 - \Delta x^2)u_{xx})_{tt}\tp$$ 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\tp$$ Writing out $$[D_xD_xD_tD_t u]^n_i$$ as $$[D_xD_x (D_tD_t u)]^n_i$$ gives \begin{align*} \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{align*} 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

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 $$$$u_{tt} = (\lambda u_x)_x\tp \tag{7.75}$$$$ The discrete counterpart to this equation, using arithmetic mean for $$\lambda$$ and centered differences, reads $$$$[D_t D_t u = D_x \overline{\lambda}^{x}D_x u]^n_i\tp \tag{7.76}$$$$ The truncation error is the residual $$R$$ in the equation $$$$[D_t D_t \uex = D_x \overline{\lambda}^{x}D_x \uex + R]^n_i\tp \tag{7.77}$$$$ The difficulty with (7.77) is how to compute the truncation error of the term $$[D_x \overline{\lambda}^{x}D_x \uex]^n_i$$.

We start by writing out the outer operator: $$$$[D_x \overline{\lambda}^{x}D_x \uex]^n_i = \frac{1}{\Delta x}\left( [\overline{\lambda}^{x}D_x \uex]^n_{i+\half} - [\overline{\lambda}^{x}D_x \uex]^n_{i-\half} \right). \tag{7.78}$$$$ With the aid of (7.5)-(7.6) and (7.21)-(7.22) we have \begin{align*} \lbrack D_x \uex \rbrack^n_{i+\half} & = \uexd{x}(x_{i+\half},t_n) + \frac{1}{24}\uexd{xxx}(x_{i+\half},t_n)\Delta x^2 + \Oof{\Delta x^4},\\ \lbrack\overline{\lambda}^{x}\rbrack_{i+\half} &= \lambda(x_{i+\half}) + \frac{1}{8}\lambda''(x_{i+\half})\Delta x^2 + \Oof{\Delta x^4},\\ [\overline{\lambda}^{x}D_x \uex]^n_{i+\half} &= (\lambda(x_{i+\half}) + \frac{1}{8}\lambda''(x_{i+\half})\Delta x^2 + \Oof{\Delta x^4})\times\\ &\quad (\uexd{x}(x_{i+\half},t_n) + \frac{1}{24}\uexd{xxx}(x_{i+\half},t_n)\Delta x^2 + \Oof{\Delta x^4})\\ &= \lambda(x_{i+\half})\uexd{x}(x_{i+\half},t_n) + \lambda(x_{i+\half}) \frac{1}{24}\uexd{xxx}(x_{i+\half},t_n)\Delta x^2 + \\ &\quad \uexd{x}(x_{i+\half}) \frac{1}{8}\lambda''(x_{i+\half})\Delta x^2 +\Oof{\Delta x^4}\\ &= [\lambda \uexd{x}]^n_{i+\half} + G^n_{i+\half}\Delta x^2 +\Oof{\Delta x^4}, \end{align*} where we have introduced the short form $$G^n_{i+\half} = (\frac{1}{24}\uexd{xxx}(x_{i+\half},t_n)\lambda((x_{i+\half}) + \uexd{x}(x_{i+\half},t_n) \frac{1}{8}\lambda''(x_{i+\half}))\Delta x^2\tp$$ Similarly, we find that $$\lbrack\overline{\lambda}^{x}D_x \uex\rbrack^n_{i-\half} = [\lambda \uexd{x}]^n_{i-\half} + G^n_{i-\half}\Delta x^2 +\Oof{\Delta x^4}\tp$$ Inserting these expressions in the outer operator (7.78) results in \begin{align*} \lbrack D_x \overline{\lambda}^{x}D_x \uex \rbrack^n_i &= \frac{1}{\Delta x}( [\overline{\lambda}^{x}D_x \uex]^n_{i+\half} - [\overline{\lambda}^{x}D_x \uex]^n_{i-\half} )\\ &= \frac{1}{\Delta x}( [\lambda \uexd{x}]^n_{i+\half} + G^n_{i+\half}\Delta x^2 - [\lambda \uexd{x}]^n_{i-\half} - G^n_{i-\half}\Delta x^2 + \Oof{\Delta x^4} )\\ &= [D_x \lambda \uexd{x}]^n_i + [D_x G]^n_i\Delta x^2 + \Oof{\Delta x^4}\tp \end{align*} The reason for $$\Oof{\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 (7.5)-(7.6) to express the $$D_x$$ operator in $$[D_x \lambda \uexd{x}]^n_i$$ as a derivative and a truncation error: $$[D_x \lambda \uexd{x}]^n_i = \frac{\partial}{\partial x}\lambda(x_i)\uexd{x}(x_i,t_n) + \frac{1}{24}(\lambda\uexd{x})_{xxx}(x_i,t_n)\Delta x^2 + \Oof{\Delta x^4}\tp$$ 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 + \Oof{\Delta x^4}\tp$$

There will be a number of terms with the $$\Delta x^2$$ factor. We lump these now into $$\Oof{\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 \uex]^n_i = \frac{\partial}{\partial x} \lambda(x_i)\uexd{x}(x_i,t_n) + \Oof{\Delta x^2}\tp$$ After having treated the $$[D_tD_t\uex]^n_i$$ term as well, we achieve $$R^n_i = \Oof{\Delta x^2} + \frac{1}{12}\uexd{tttt}(x_i,t_n)\Delta t^2 \tp$$ 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.

## Linear wave equation in 2D/3D

The two-dimensional extension of (7.68) takes the form $$$$\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], \tag{7.79}$$$$ 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 (7.79) can be written $$$$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], \tag{7.80}$$$$ in 2D, while the 3D version reads $$$$u_{tt} = c^2(u_{xx} + u_{yy} + u_{zz}) + f(x,y,z,t), \tag{7.81}$$$$ 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 (7.17)-(7.18) yields the scheme $$$$[D_t D_t u = c^2(D_xD_x u + D_yD_y u) + f]^n_{i,j,k} \tp \tag{7.82}$$$$ The truncation error is found from $$$$[D_t D_t \uex = c^2(D_xD_x \uex + D_yD_y \uex) + f + R]^n \tp \tag{7.83}$$$$ 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 \begin{align} R^n_{i,j,k} & = [\frac{1}{12}\uexd{tttt}\Delta t^2 - c^2\frac{1}{12}\left( \uexd{xxxx}\Delta x^2 + \uexd{yyyy}\Delta x^2 + \uexd{zzzz}\Delta z^2\right)]^n_{i,j,k} + \tag{7.84}\\ &\quad \Oof{\Delta t^4,\Delta x^4,\Delta y^4,\Delta z^4}\nonumber \tp \end{align}