Wave equations
Linear wave equation in 1D
The standard, linear wave equation in 1D for a function \(u(x,t)\) reads
\[\tag{721}
\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 in \([0,L]\).
The equation can also be more compactly written as
\[\tag{722}
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
\[ \begin{align}\begin{aligned}\tag{723}
[D_t D_t u = c^2 D_xD_x u + f]^n_i\\ {\thinspace .}\end{aligned}\end{align} \]
Inserting the exact solution \({u_{\small\mbox{e}}}(x,t)\) in (723)
makes this function fulfill the equation if we add the
term \(R\):
\[\tag{724}
[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 (670)-(671) 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 (724) now becomes
\[\begin{split}\begin{align*}
{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{align*}\end{split}\]
Because \({u_{\small\mbox{e}}}\) fulfills the partial differential equation (PDE)
(722), the first, third, and fifth term cancel out,
and we are left with
\[\tag{725}
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 (723) 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 (725)? The starting point is
\[\tag{726}
[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 (725)
again, but now with \(C\):
\[\tag{727}
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 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}{\thinspace .}\]
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}{\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}\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*}\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
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
\[\tag{728}
u_{tt} = (\lambda u_x)_x{\thinspace .}\]
The discrete counterpart to this equation, using arithmetic mean for
\(\lambda\) and centered differences, reads
\[\tag{729}
[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
\[\tag{730}
[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 with (730)
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:
\[\tag{731}
[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 (658)-(659)
and (674)-(675) we have
\[\begin{split}\begin{align*}
\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{align*}\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 (731)
results in
\[\begin{split}\begin{align*}
\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{align*}\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 (658)-(659)
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 (721) takes the form
\[\tag{732}
\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
(732) can be written
\[\tag{733}
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
\[\tag{734}
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 (670)-(671)
yields the scheme
\[\tag{735}
[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
\[\tag{736}
[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
\[\tag{737}
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 .}\]
