The standard, linear wave equation in 1D for a function \( u(x,t) \) reads $$ \begin{equation} \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{63} \end{equation} $$ where \( c \) is the constant wave velocity of the physical medium \( [0,L] \). The equation can also be more compactly written as $$ \begin{equation} u_{tt} = c^2u_{xx} + f,\quad x\in (0, L),\ t\in (0,T], \tag{64} \end{equation} $$ Centered, second-order finite differences are a natural choice for discretizing the derivatives, leading to $$ \begin{equation} [D_t D_t u = c^2 D_xD_x u + f]^n_i \tag{65} \tp \end{equation} $$
Inserting the exact solution \( \uex(x,t) \) in (65) makes this function fulfill the equation if we add the term \( R \): $$ \begin{equation} [D_t D_t \uex = c^2 D_xD_x \uex + f + R]^n_i \tag{68} \end{equation} $$
Our purpose is to calculate the truncation error \( R \). From (15)-(16) 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 (68) 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) (64), the first, third, and fifth terms cancel out, and we are left with $$ \begin{equation} 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{67} \end{equation} $$ showing that the scheme (65) is of second order in the time and space mesh spacing.
Can we add correction terms to the PDE and increase the order of \( R^n_i \) in (67)? The starting point is $$ \begin{equation} [D_t D_t \uex = c^2 D_xD_x \uex + f + C + R]^n_i \tag{68} \end{equation} $$ From the previous analysis we simply get (67) again, but now with \( C \): $$ \begin{equation} 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{69} \end{equation} $$ 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 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}\tp$$ 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}\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.
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 $$ \begin{equation} u_{tt} = (\lambda u_x)_x\tp \tag{70} \end{equation} $$ The discrete counterpart to this equation, using arithmetic mean for \( \lambda \) and centered differences, reads $$ \begin{equation} [D_t D_t u = D_x \overline{\lambda}^{x}D_x u]^n_i\tp \tag{71} \end{equation} $$ The truncation error is the residual \( R \) in the equation $$ \begin{equation} [D_t D_t \uex = D_x \overline{\lambda}^{x}D_x \uex + R]^n_i\tp \tag{72} \end{equation} $$ The difficulty in the present 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: $$ \begin{equation} [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{73} \end{equation} $$ With the aid of (3)-(4) and (19)-(20) 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 (73) 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 (3)-(4) 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.
The two-dimensional extension of (63) takes the form $$ \begin{equation} \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{74} \end{equation} $$ 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 (74) can be written $$ \begin{equation} 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{75} \end{equation} $$ in 2D, while the 3D version reads $$ \begin{equation} u_{tt} = c^2(u_{xx} + u_{yy} + u_{zz}) + f(x,y,z,t), \tag{76} \end{equation} $$ 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 (15)-(16) yields the scheme $$ \begin{equation} [D_t D_t u = c^2(D_xD_x u + D_yD_y u) + f]^n_{i,j,k} \tp \end{equation} $$ The truncation error is found from $$ \begin{equation} [D_t D_t \uex = c^2(D_xD_x \uex + D_yD_y \uex) + f + R]^n \tp \end{equation} $$ 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} +\\ &\quad \Oof{\Delta t^4,\Delta x^4,\Delta y^4,\Delta z^4}\nonumber \tp \end{align} $$