An adaptive multiresolution interior penalty discontinuous Galerkin method for
wave equations in second order form

Juntao Huang 111Department of Mathematics, Michigan State University, East
Lansing, MI 48824, USA. E-mail<EMAIL_ADDRESS>Yuan Liu222Department of
Mathematics, Statistics and Physics, Wichita State University, Wichita, KS
67260, USA. E-mail<EMAIL_ADDRESS>Research supported in part by a
grant from the Simons Foundation (426993, Yuan Liu).  Wei Guo 333Department of
Mathematics and Statistics, Texas Tech University, Lubbock, TX, 70409. E-mail:
<EMAIL_ADDRESS>Research is supported by NSF grant DMS-1830838 Zhanjing
Tao 444School of Mathematics, Jilin University, Changchun, Jilin 130012,
China<EMAIL_ADDRESS>Corresponding author Yingda Cheng 555Department of
Mathematics, Department of Computational Mathematics, Science and Engineering,
Michigan State University, East Lansing, MI 48824, USA. E-mail:
<EMAIL_ADDRESS>Research is supported by NSF grants DMS-1453661 and
DMS-1720023

Abstract

In this paper, we propose a class of adaptive multiresolution (also called
adaptive sparse grid) discontinuous Galerkin (DG) methods for simulating
scalar wave equations in second order form in space. The two key ingredients
of the schemes include an interior penalty DG formulation in the adaptive
function space and two classes of multiwavelets for achieving multiresolution.
In particular, the orthonormal Alpert’s multiwavelets are used to express the
DG solution in terms of a hierarchical structure, and the interpolatory
multiwavelets are further introduced to enhance computational efficiency in
the presence of variable wave speed or nonlinear source. Some theoretical
results on stability and accuracy of the proposed method are presented.
Benchmark numerical tests in 2D and 3D are provided to validate the
performance of the method.

Key Words: Sparse grid; Multiresolution; Interior Penalty Discontinuous
Galerkin Method; Wave Equation; Adaptivity.

##  1 Introduction

Wave propagation, governed by the wave equation, is ubiquitous in science and
engineering, such as sound waves, light waves, and water waves propagating in
acoustics, electromagnetics and geoscience. Designing efficient and robust
numerical methods to solve the wave equation is of fundamental and practical
importance in those applications. The goal of this work is to design a class
of numerical solvers that are adaptive, high order accurate, and more
importantly, work efficiently in high dimensions. In particular, we develop a
class of adaptive multiresolution (also called adaptive sparse grid)
discontinuous Galerkin (DG) method for the following model second-order wave
equation

$u_{tt}=\nabla\cdot(c^{2}(\mathbf{x})\nabla u)+f$ (1.1)

on the bounded domain $\Omega=[0,1]^{d}$ in arbitrary $d$ dimensions, subject
to initial conditions

$u(\mathbf{x},0)=u_{0}(\mathbf{x}),\qquad
u_{t}(\mathbf{x},0)=v_{0}(\mathbf{x}).$ (1.2)

We assume that the wave speed $c(\mathbf{x})$ is piecewise smooth and bounded
below and above uniformly, i.e., $0<C_{*}\leq c^{2}(\mathbf{x})\leq
C^{*}<\infty$ . For simplicity, we only consider periodic, Dirichlet and
Neumann boundary conditions in this paper. Extensions to more complicated
domains and other types of boundary conditions will be considered in the
future work.

A vast amount of numerical methods have been developed in the literature on
the numerical approximations of the wave equation, including finite difference
discretization [23, 16, 24, 36], spectral and spectral element discretization
[18, 33, 38] and finite element discretization [29, 1], to name a few. As a
special class of finite element discretization, the DG methods [32, 15] have
become very popular recently in approximating partial differential equations
(PDEs) due to their distinguished advantages in handling geometry, boundary
conditions and accommodating adaptivity. In the context of the wave
simulations, DG methods have been successfully developed for simulating wave
equations in first-order form [25, 30, 40], second-order form [20, 41, 13, 3],
and with $hp$-adaptivity [17]. In this paper, we utilize the symmetric
interior penalty DG (IPDG) method [4] for wave equation in second order form
[20], though our framework can work with other types of DG schemes.

Adaptivity is crucial for efficient simulations of the wave equation due to
the multiscale nature of the solution structures. The well-known adaptive mesh
refinement (AMR) [7, 6] adjusts the computational grid adaptively to track
small scale features of the underlying problems, improving computational
efficiency significantly. AMR has been incorporated in various software
framework and packages to simulate wave propagation with great success [9,
11]. In contrast, this paper considers adaptive simulations in the
multiresolution sense. The main idea of multiresolution analysis (MRA) [31] is
to explore mesh hierarchy, which induces nested polynomial approximation
spaces to accelerate the computation and in the mean time circumvents the need
for _a posteriori_ error indicators. MRA is also the foundation of sparse grid
methods [10], which is known as a popular dimension reduction technique for
solving high dimensional problems. As a continuation of our previous research
for adaptive multiresolution (also called adaptive sparse grid) DG methods
[22, 27] for first order equations, this paper develops an adaptive
multiresolution IPDG solver for 2D and 3D scalar wave equations (1.1). In
particular, we employ the Alpert’s multiwavelets as the DG bases in the IPDG
formulation, following the approach proposed in [39, 21, 22] for linear
equations, together with the interpolatory multiwavelets for efficiently
computing variable wave speed problems as done in [27] for nonlinear
hyperbolic conservation laws. We refer the readers to [27] for more details on
the background of adaptive multiresolution DG methods [12, 26]. It is worth
noting that a fast matrix-vector multiplication algorithm [34, 42] is
essential for efficient implementation of the method with varying wave speed.
We conducts error analysis for the semi-discrete formulation for the scheme
with and without interpolations. First, when the sparse grid piecewise
polynomial space of degree $k$ is employed in the IPDG formulation as in [39],
the newly proposed method converges with order $k$ and a polylogarithmic
factor in the energy norm for sufficiently smooth problems with constant
coefficients. Second, in the case of smooth problems with variable
coefficients, the proposed interpolatory technique ensures a high order local
truncation error and hence preserve the original accuracy of the scheme given
sufficient high order accuracy of interpolation. Numerical experiments in 2D
and 3D verify the accuracy of the methods. In particular, the adaptive scheme
is demonstrated to capture the fine scale structure presented in inhomogeneous
media.

The rest of the paper is organized as follows. In Section 2, we review
Alpert’s and interpolatory multiwavelets. Section 3 describes the numerical
schemes with details on some theoretical results and implementations. Section
4 contains numerical examples. In Section 5, we make conclusions and discuss
future work. Appendix collects detailed formulas of interpolatory
multiwavelets used in this paper.

##  2 MRA and multiwavelets

In this section, we first review the $L^{2}$ orthonormal Alpert’s
multiwavelets [2] and the sparse grid DG finite element space [39, 21]. Next,
we review the interpolatory multiwavelets proposed in [37], which has been
used for the calculation of nonlinear conservation laws in [27].

###  2.1 Alpert’s multiwavelets

In this subsection, we review the construction of sparse grid DG finite
element space based on Alpert’s multiwavelets [2]. For a unit domain
$\Omega=[0,1]$ in 1D, we define a set of nested grids, where the $n$-th level
grid $\Omega_{n}$ consists of $2^{n}$ uniform cells

$I_{n}^{j}=(2^{-n}j,2^{-n}(j+1)],\quad j=0,\ldots,2^{n}-1$

for $n\geq 0.$ For notational convenience, we also denote $I_{-1}=[0,1].$ The
piecewise polynomial space of degree at most $k\geq 1$ on the $n$-th level
grid $\Omega_{n}$ for $n\geq 0$ is denoted by

$V_{n}^{k}:=\\{v:v\in P^{k}(I_{n}^{j}),\,\forall\,j=0,\ldots,2^{n}-1\\}.$
(2.1)

Because of the nested structure

$V_{0}^{k}\subset V_{1}^{k}\subset V_{2}^{k}\subset V_{3}^{k}\subset\cdots,$

we define the multiwavelet subspace $W_{n}^{k}$, $n=1,2,\ldots$ as the
orthogonal complement of $V_{n-1}^{k}$ in $V_{n}^{k}$ with respect to the
$L^{2}$ inner product on $[0,1]$, i.e.,

$V_{n-1}^{k}\oplus W_{n}^{k}=V_{n}^{k},\quad W_{n}^{k}\perp V_{n-1}^{k}.$

Denote $W_{0}^{k}:=V_{0}^{k}$, we have $V_{n}^{k}=\bigoplus_{0\leq l\leq
n}W_{l}^{k}$. A set of orthonormal basis can be defined on $W_{l}^{k}$ as
follows. When $l=0$, the basis $v^{0}_{i,0}(x)$, $i=0,\ldots,k$ are the
normalized shifted Legendre polynomials in $[0,1]$. When $l>0$, the Alpert’s
orthonormal multiwavelets are employed [2] as the bases and denoted by

$v^{j}_{i,l}(x),\quad i=0,\ldots,k,\quad j=0,\ldots,2^{l-1}-1.$

We then follow a tensor-product approach to construct the hierarchical finite
element space in multi-dimensional space. Denote
$\mathbf{l}=(l_{1},\cdots,l_{d})\in\mathbb{N}_{0}^{d}$ as the mesh level in a
multivariate sense, where $\mathbb{N}_{0}$ denotes the set of nonnegative
integers, we can define the tensor-product mesh grid
$\Omega_{\mathbf{l}}=\Omega_{l_{1}}\otimes\cdots\otimes\Omega_{l_{d}}$ and the
corresponding mesh size $h_{\mathbf{l}}=(h_{l_{1}},\cdots,h_{l_{d}}).$ Based
on the grid $\Omega_{\mathbf{l}}$, we denote
$I_{\mathbf{l}}^{\mathbf{j}}=\\{\mathbf{x}:x_{m}\in(h_{m}j_{m},h_{m}(j_{m}+1)),m=1,\cdots,d\\}$
as an elementary cell, and

${\bf V}_{\mathbf{l}}^{k}:=\\{{\bf v}:{\bf v}\in
Q^{k}(I^{\mathbf{j}}_{\mathbf{l}}),\,\,\mathbf{0}\leq\mathbf{j}\leq
2^{\mathbf{l}}-\mathbf{1}\\}=V_{l_{1},x_{1}}^{k}\times\cdots\times
V_{l_{d},x_{d}}^{k}$

as the tensor-product piecewise polynomial space, where
$Q^{k}(I^{\mathbf{j}}_{\mathbf{l}})$ represents the collection of polynomials
of degree up to $k$ in each dimension on cell $I^{\mathbf{j}}_{\mathbf{l}}$.
If we use equal mesh refinement of size $h_{N}=2^{-N}$ in each coordinate
direction, the grid and space will be denoted by $\Omega_{N}$ and ${\bf
V}_{N}^{k}$, respectively. Based on a tensor-product construction, the multi-
dimensional increment space can be defined as

$\mathbf{W}_{\mathbf{l}}^{k}=W_{l_{1},x_{1}}^{k}\times\cdots\times
W_{l_{d},x_{d}}^{k}.$

The basis functions in multi-dimensions are defined as

$v^{\mathbf{j}}_{\mathbf{i},\mathbf{l}}(\mathbf{x}):=\prod_{m=1}^{d}v^{j_{m}}_{i_{m},l_{m}}(x_{m}),$
(2.2)

for $\mathbf{l}\in\mathbb{N}_{0}^{d}$, $\mathbf{j}\in
B_{\mathbf{l}}:=\\{\mathbf{j}\in\mathbb{N}_{0}^{d}:\,\mathbf{0}\leq\mathbf{j}\leq\max(2^{\mathbf{l}-\mathbf{1}}-\mathbf{1},\mathbf{0})\\}$
and $\mathbf{1}\leq\mathbf{i}\leq\mathbf{k}+\mathbf{1}$. The orthonormality of
the bases can be easily verified.

Using the notation of

$|\mathbf{l}|_{1}:=\sum_{m=1}^{d}l_{m},\qquad|\mathbf{l}|_{\infty}:=\max_{1\leq
m\leq d}l_{m}.$

and the same component-wise arithmetic operations and relations as defined in
[39], we reach the decomposition

${\bf V}_{N}^{k}=\bigoplus_{\begin{subarray}{c}|\mathbf{l}|_{\infty}\leq N\\\
\mathbf{l}\in\mathbb{N}_{0}^{d}\end{subarray}}\mathbf{W}_{\mathbf{l}}^{k}.$
(2.3)

On the other hand, a standard choice of sparse grid space [39, 21] is

$\hat{{\bf V}}_{N}^{k}=\bigoplus_{\begin{subarray}{c}|\mathbf{l}|_{1}\leq N\\\
\mathbf{l}\in\mathbb{N}_{0}^{d}\end{subarray}}\mathbf{W}_{\mathbf{l}}^{k}\subset{\bf
V}_{N}^{k}.$ (2.4)

We skip the discussions on the details with regard to the property of the
space, but refer the readers to [39, 21]. In Section 3, we will describe the
adaptive scheme which adapts a subspace of ${\bf V}_{N}^{k}$ according to the
numerical solution, hence offering more flexibility and efficiency.

###  2.2 Interpolatory multiwavelets

Alpert’s multiwavelets described in Section 2.1 are associated with the
$L^{2}$ projection operator. The idea of interpolatory multiwavelet bases [37]
is based on interpolation operators and is essential for the computation of
variable coefficient problems. In this work, only Lagrange interpolation is
considered, while we note that Hermite interpolation can be used. The details
are provided below.

We define the set of interpolation points on the interval $I=[0,1]$ at mesh
level 0 by $X_{0}=\\{x_{i}\\}_{i=0}^{M}\subset I$. Here, the number of points
in $X_{0}$ is $(M+1)$. We defer the discussion of the relations between $M$
and $k$ to Section 3.2.

The interpolation points at mesh level $n\geq 1$, $X_{n}$ can be obtained
correspondingly as

$X_{n}=\\{x_{i,n}^{j}:=2^{-n}(x_{i}+j),\quad i=0,\dots,M,\quad
j=0,\dots,2^{n}-1\\}.$

We require the points to be nested, i.e.

$X_{0}\subset X_{1}\subset X_{2}\subset X_{3}\subset\cdots.$ (2.5)

This can be achieved by requiring $X_{0}\subset X_{1}$.

Given the nodes, we define the basis functions on the zeroth level grid as
Lagrange interpolation polynomials of degree $\leq M$ which satisfy the
property:

$\phi_{i}(x_{i^{\prime}})=\delta_{ii^{\prime}},$

for $i,i^{\prime}=0,\dots,M$. It is easy to see that
$\textrm{span}\\{\phi_{i},i=0,\dots,M\\}=V_{0}^{M}.$ With the basis function
at mesh level zero, we can define the basis functions at mesh level $n\geq 1$:

$\phi_{i,n}^{j}:=\phi_{i}(2^{n}x-j),\quad i=0,\dots,M,\quad j=0,\dots,2^{n}-1$

which form a complete basis set for $V_{n}^{M}.$

We now introduce the hierarchical representations and the interpolatory
multiwavelets. Define $\tilde{X}_{0}:=X_{0}$ and
$\tilde{X}_{n}:=X_{n}\backslash X_{n-1}$ for $n\geq 1$, then we have the
decomposition

$X_{n}=\tilde{X}_{0}\cup\tilde{X}_{1}\cup\cdots\cup\tilde{X}_{n}.$

Denote the points in $\tilde{X}_{1}$ by
$\tilde{X}_{1}=\\{\tilde{x}_{i}\\}_{i=0}^{M}$. Then the points in
$\tilde{X}_{n}$ for $n\geq 1$ can be represented by

$\tilde{X}_{n}=\\{\tilde{x}_{i,n}^{j}:=2^{-(n-1)}(\tilde{x}_{i}+j),\quad
i=0,\dots,M,\quad j=0,\dots,2^{n-1}-1\\}.$

For notational convenience, we let $\tilde{W}_{0}^{M}:=V_{0}^{M}.$ The
increment function space $\tilde{W}_{n}^{M}$ for $n\geq 1$ is introduced as a
function space that satisfies

$V_{n}^{M}=V_{n-1}^{M}\oplus\tilde{W}_{n}^{M},$ (2.6)

and is defined through the multiwavelets $\psi_{i}\in V_{1}^{M}$ that
satisfies

$\psi_{i}(x_{i^{\prime}})=0,\quad\psi_{i}(\tilde{x}_{i^{\prime}})=\delta_{i,i^{\prime}},$

for $i,i^{\prime}=0,\dots,M$. Then $\tilde{W}_{n}^{M}$ is given by

$\tilde{W}_{n}^{M}=\textrm{span}\\{\psi_{i,n}^{j},\quad i=0,\dots,M,\quad
j=0,\dots,2^{n-1}-1\\}$

where $\psi_{i,n}^{j}(x):=\psi_{i}(2^{n-1}x-j)$.

The multi-dimensional construction follows similar lines as in Section 2.1. We
let

$\tilde{\mathbf{W}}_{\mathbf{l}}^{M}=\tilde{W}_{l_{1},x_{1}}^{M}\times\cdots\times\tilde{W}_{l_{d},x_{d}}^{M},$

then

${\bf V}_{N}^{M}=\bigoplus_{\begin{subarray}{c}|\mathbf{l}|_{\infty}\leq N\\\
\mathbf{l}\in\mathbb{N}_{0}^{d}\end{subarray}}\tilde{\mathbf{W}}_{\mathbf{l}}^{M},$

while the sparse grid approximation space is

$\hat{{\bf V}}_{N}^{M}=\bigoplus_{\begin{subarray}{c}|\mathbf{l}|_{1}\leq N\\\
\mathbf{l}\in\mathbb{N}_{0}^{d}\end{subarray}}\tilde{\mathbf{W}}_{\mathbf{l}}^{M}.$

Note that the construction by Alpert’s multiwavelet and the interpolatory
multiwavelet gives the same sparse grid space because of the same nested
structure. Finally, the interpolation operator in multidimension is defined as
$\mathcal{I}^{M}_{N}:C(\Omega)\rightarrow\mathbf{V}^{M}_{N}$:

$\displaystyle\mathcal{I}^{M}_{N}[f](\mathbf{x})=\sum_{\begin{subarray}{c}\left|\mathbf{n}\right|_{\infty}\leq
N\\\
\mathbf{0}\leq\mathbf{j}\leq\max(2^{\mathbf{n}-1}-\mathbf{1},\mathbf{0})\\\
\mathbf{0}\leq\mathbf{i}\leq\mathbf{M}\end{subarray}}b^{\mathbf{j}}_{\mathbf{i},\mathbf{n}}\psi^{\mathbf{j}}_{\mathbf{i},\mathbf{n}}(\mathbf{x}),$

where the multi-dimensional basis functions
$\psi^{\mathbf{j}}_{\mathbf{i},\mathbf{n}}(\mathbf{x})$ are defined in the
same approach as (2.2) by tensor products:

$\psi^{\mathbf{j}}_{\mathbf{i},\mathbf{n}}(\mathbf{x}):=\prod_{m=1}^{d}\psi^{j_{m}}_{i_{m},n_{m}}(x_{m}).$
(2.7)

For the sparse grid space $\hat{{\bf V}}_{N}^{M}$ or any adaptively chosen
subspace of ${\bf V}_{N}^{M},$ the interpolation operator, which is denoted by
$\mathcal{I}_{h}$ in later sections, can be defined accordingly, by taking
only multiwavelet basis functions that belong to that space. For completeness,
we collect the detailed formulas of the interpolation points and the
associated interpolatory multiwavelets used in this work in the Appendix.

##  3 Adaptive multiresolution DG scheme

In this section, we construct our numerical schemes for $d$-dimensional wave
equation (1.1). We start by reviewing the semi-discrete IPDG formulation and
its properties in Section 3.1. For variable wave speed, schemes with
multiresolution interpolation are described in Section 3.2. Time stepping,
adaptivity and fast implementations are discussed in Sections 3.3 and 3.4.

###  3.1 Semi-discrete scheme

We use the IPDG formulation [20] for solving (1.1). Namely, we look for
$u_{h}\in{\bf V},$ such that for any test function $v\in{\bf V}$,

$\displaystyle\int_{\Omega}(u_{h})_{tt}v\ d\textbf{x}+B(u_{h},v)=L(v).$ (3.1)

where the bilinear form is defined as

$\displaystyle B(u_{h},v)$ $\displaystyle=\int_{\Omega}c^{2}\nabla
u_{h}\cdot\nabla v\ d\textbf{x}-\sum_{e\in\Gamma}\int_{e}\\{c^{2}\nabla
u_{h}\\}\cdot[v]\ ds-\sum_{e\in\Gamma}\int_{e}\\{c^{2}\nabla v\\}\cdot[u_{h}]\
ds$ (3.2)
$\displaystyle+\sum_{e\in\Gamma}\frac{\sigma}{h_{N}}\int_{e}[u_{h}]\cdot[v]\
ds$

and

$L(v)=\int_{\Omega}fv\ d\textbf{x}$ (3.3)

for periodic or homogeneous Dirichlet boundary condition, and

$L(v)=\int_{\Omega}fvd\textbf{x}+\sum_{e\in\Gamma_{D}}\int_{e}(-c^{2}\nabla
v\cdot\textbf{n}+\frac{\sigma}{h_{N}}v)g_{D}ds+\sum_{e\in\Gamma_{N}}\int_{e}g_{N}vds$
(3.4)

for Dirichlet and Neumann boundary conditions
$u(\mathbf{x},t)|_{\mathbf{x}\in\Gamma_{D}}=g_{D}$ and $\nabla
u(\mathbf{x},t)\cdot\textbf{n}|_{\mathbf{x}\in\Gamma_{N}}=g_{N}$. $\Gamma$ is
the union of the boundaries for all the elements in the partition
$\Omega_{N}$, and $\sigma$ is the penalty parameter depending on the dimension
$d.$ The average and jump are defined as,

$\displaystyle[q]=q^{-}\textbf{n}^{-}+q^{+}\textbf{n}^{+},\qquad$
$\displaystyle\\{q\\}=\frac{1}{2}(q^{-}+q^{+}),$
$\displaystyle[\textbf{q}]=\textbf{q}^{-}\cdot\textbf{n}^{-}+\textbf{q}^{+}\cdot\textbf{n}^{+},\qquad$
$\displaystyle\\{\textbf{q}\\}=\frac{1}{2}(\textbf{q}^{-}+\textbf{q}^{+}).$
(3.5)

where n is the unit normal. ‘-’ and ‘+’ represent that the directions of the
vector point to interior and exterior at $e$ respectively. If $e$ is part of
the boundary, then we let $[q]=q\textbf{n}$ (n is the outward unit normal) and
$\\{\textbf{q}\\}=\textbf{q}$.

Depending on the choice of space ${\bf V}$, various IPDG methods with distinct
properties are obtained. If ${\bf V}={\bf V}^{k}_{N},$ we recover the IPDG
scheme in [20] on tensor-product meshes. If ${\bf V}=\hat{{\bf V}}^{k}_{N},$
then we obtain the sparse grid IPDG method. If ${\bf V}$ is chosen adaptively
as described in Section 3.3, we have the adaptive multiresolution scheme. Note
that besides the IPDG formulation, other DG formulations can be used as well,
such as the local DG method [13] and the energy-based DG method [3]. The main
novelty of this work is the choice of the multiresolution polynomial space
which is not tied specifically to the weak formulation in use.

For completeness, we now review some properties of the semi-discrete IPDG
scheme (3.1). Define the discrete energy of wave propagation by

$\displaystyle E_{h}(t):=\frac{1}{2}\left\|\frac{\partial u_{h}}{\partial
t}\right\|^{2}+\frac{1}{2}B(u_{h},u_{h}),$ (3.6)

Then the stability inherently holds true since the bilinear form
$B(\cdot,\cdot)$ is symmetric and coercive:

######  Theorem 3.1 (Energy stability [20]).

The discrete energy (3.6) is conserved by semi-discrete DG scheme (3.1)-(3.3)
when $f=0$ with periodic boundary condition for arbitrary choice of space
including ${\bf V}={\bf V}^{k}_{N}$ and ${\bf V}=\hat{{\bf V}}^{k}_{N}$.

We then review some results in the error estimates [20], and extend it to the
sparse grid method with ${\bf V}=\hat{{\bf V}}^{k}_{N}$ based on the
approximation properties of the space $\hat{{\bf V}}^{k}_{N}$ in [21]. We use
$\|\cdot\|$ to represent the standard $L^{2}$ norm on $\Omega$ or
$\Omega_{N},$ $\|\cdot\|_{L^{2}(\Gamma)}$ to represent the $L^{2}$ norm on the
collection of the cell interfaces of the mesh $\Omega_{N}:\Gamma,$ and define
the energy norm of a function $v\in H^{2}(\Omega_{N})$ as

$\displaystyle|||v|||^{2}:=\int_{\Omega}|\nabla
v|^{2}\,d\mathbf{x}\,+\sum_{\begin{subarray}{c}e\in\Gamma\end{subarray}}h_{N}\int_{e}\left\\{\frac{\partial
v}{\partial{\bf
n}}\right\\}^{2}\,ds\,+\sum_{\begin{subarray}{c}e\in\Gamma\end{subarray}}\frac{1}{h_{N}}\int_{e}[v]^{2}\,ds.$
(3.7)

Some basic properties of the bilinear operator $B(\cdot,\cdot)$ are listed
below.

######  Lemma 3.2 (Boundedness [4, 5]).

There exists a positive constant $C_{b}$, depending only on $C^{*},\sigma$,
such that

$|B(w,v)|\leq C_{b}|||w|||\cdot|||v|||,\quad\forall\,w,v\in
H^{2}(\Omega_{N}).$

######  Lemma 3.3 (Coercivity [4, 5]).

When $\sigma$ is taken large enough, there exists a positive constant $C_{s}$
depending only on $C_{*}$, such that

$B(v,v)\geq C_{s}|||v|||^{2},\quad\forall\,v\in\hat{{\bf V}}_{N}^{k}.$

Then we arrive at the following error estimate.

######  Theorem 3.4 (Error estimate in energy norm for sparse grid IPDG
method).

Let $u$ be the solution of (1.1)-(1.2) satisfying $u\in
L^{\infty}(0,T;\mathcal{H}^{p+1}(\Omega))$, $u_{t}\in
L^{\infty}(0,T;\mathcal{H}^{p+1}(\Omega))$, $u_{tt}\in
L^{1}(0,T;\mathcal{H}^{p}(\Omega)).$ $u_{h}$ is the semi-discrete DG solution
obtained by (3.1)-(3.3) with ${\bf V}=\hat{{\bf V}}^{k}_{N}$ and the initial
condition $u_{h}(0)=\mathbf{P}u_{0}$ and $(u_{h})_{t}(0)=\mathbf{P}v_{0}$,
where $\mathbf{P}$ denotes the $L^{2}$ projection of a function onto the space
$\hat{{\bf V}}^{k}_{N}$. Then for $k\geq 1$ and any $1\leq
q\leq\min\\{p,k\\}$, the error $e=u_{h}-u$ satisfies the estimation

$\displaystyle||e_{t}||_{L^{\infty}(0,T;L^{2}(\Omega))}$
$\displaystyle+\sup_{t\in[0,T]}|||e|||\leq C(||e_{t}(0)||+|||e(0)|||)$ (3.8)
$\displaystyle+C\left|\log_{2}h_{N}\right|^{d}h_{N}^{q}\left(|u|_{L^{\infty}(0,T;\mathcal{H}^{q+1}(\Omega))}+T|u_{t}|_{L^{\infty}(0,T;\mathcal{H}^{q+1}(\Omega))}+|u_{tt}|_{L^{\infty}(0,T;\mathcal{H}^{q}(\Omega))}\right)$

where the dimension $d\geq 2.$ $|\cdot|_{\mathcal{H}^{q+1}(\Omega)}$ denotes
mixed derivative norm of a function and was defined in [21]. Here and below,
$C$ denotes a generic constant that does not depend on $h_{N}$ or the solution
$u.$

###### Proof.

Following [20], we let ${\bf V}(h)=H^{1}(\Omega)+\hat{{\bf V}}^{k}_{N},$ and
for any $v\in{\bf V}(h),$ we define the lifted function
$\mathcal{L}_{c}(v)\in(\hat{{\bf V}}^{k}_{N})^{d}$ by requiring

$\displaystyle\int_{\Omega}\mathcal{L}_{c}(v)\cdot
wd\textbf{x}=\sum_{e\in\Gamma}\int_{e}[v]\cdot\\{c^{2}w\\}ds,\qquad
w\in(\hat{{\bf V}}^{k}_{N})^{d}.$ (3.9)

Using similar arguments as in Lemma 4.3 in [20], we conclude the lifting
operator $\mathcal{L}_{c}$ exists and is stable in the DG norm. Then the
auxiliary bilinear form can be introduced as

$\displaystyle\hat{B}(u,v)$ $\displaystyle=\int_{\Omega}c^{2}\nabla
u\cdot\nabla v\ d\textbf{x}-\int_{\Omega}\mathcal{L}_{c}(u)\cdot\nabla{v}\
ds-\int_{\Omega}\mathcal{L}_{c}(v)\cdot\nabla{u}\ ds$ (3.10)
$\displaystyle+\sum_{e\in\Gamma}\frac{\sigma}{h_{N}}\int_{e}[u]\cdot[v]\ ds.$

$\hat{B}(u,v)$ can be viewed as an extension of the wave operator and bilinear
form $B(u,v)$ to the space ${\bf V}(h)\times{\bf V}(h)$, since

$\displaystyle\hat{B}(u,v)=B(u,v)\quad\text{on}\quad\hat{{\bf
V}}^{k}_{N}\times\hat{{\bf V}}^{k}_{N},$ (3.11)
$\displaystyle\hat{B}(u,v)=\int_{\Omega}c^{2}\nabla u\cdot\nabla
vd\textbf{x}-\int_{\partial\Omega}(c^{2}\nabla u)v\cdot{\bf
n}ds\quad\text{on}\quad H^{1}(\Omega)\times H^{1}(\Omega).$ (3.12)

Moreover, it can be verified that

$\displaystyle\hat{B}(u,v)\leq C_{b}|||u|||\cdot|||v|||,$ (3.13)
$\displaystyle\hat{B}(u,u)\geq C_{s}|||u|||^{2}.$

Similar to Lemma 4.5 in [20], $e$ satisfies the equation

$\displaystyle(e_{tt},v)+\hat{B}(e,v)=r_{h}(u,v),\quad\forall v\in\hat{{\bf
V}}^{k}_{N}$ (3.14)

where

$\displaystyle r_{h}(u,v)=\sum_{e\in\Gamma}\int_{e}[v]\cdot\\{c^{2}\nabla
u-c^{2}\mathbf{P}(\nabla u)\\}ds.$ (3.15)

Therefore, we will have

$\displaystyle\frac{1}{2}\frac{d}{dt}[||e_{t}||^{2}+\hat{B}(e,e)]$
$\displaystyle=(e_{tt},e_{t})+\hat{B}(e,e_{t})$ (3.16)
$\displaystyle=(e_{tt},(u-\mathbf{P}u)_{t})+\hat{B}(e,(u-\mathbf{P}u)_{t})+r_{h}(u,(\mathbf{P}u-u_{h})_{t}).$

Integrating (3.16) over $[0,s]$ for any $s\in[0,T]$ yields

$\displaystyle\frac{1}{2}||e_{t}(s)||^{2}$
$\displaystyle+\frac{1}{2}\hat{B}(e(s),e(s))=\frac{1}{2}||e_{t}(0)||^{2}+\frac{1}{2}\hat{B}(e(0),e(0))+\int^{s}_{0}(e_{tt},(u-\mathbf{P}u)_{t})dt$
(3.17)
$\displaystyle+\int^{s}_{0}\hat{B}(e,(u-\mathbf{P}u)_{t})dt+\int^{s}_{0}r_{h}(u,(\mathbf{P}u-u_{h})_{t})dt.$

Because

$\displaystyle\int^{s}_{0}(e_{tt},(u-\mathbf{P}u)_{t})dt=-\int^{s}_{0}(e_{t},(u-\mathbf{P}u)_{tt})dt+[(e_{t},(u-\mathbf{P}u)_{t})]^{t=s}_{t=0},$
(3.18)

and the inequalities (3.13) hold, together with Holder’s inequalities, we will
have

$\displaystyle\frac{1}{2}||e_{t}(s)||^{2}+\frac{1}{2}C_{s}|||e(s)|||^{2}$
$\displaystyle\leq\frac{1}{2}||e_{t}(0)||^{2}+\frac{1}{2}C_{b}|||e(0)|||^{2}$
(3.19)
$\displaystyle+||e_{t}||_{L^{\infty}(0,T;L^{2}(\Omega))}(||(u-\mathbf{P}u)_{tt}||_{L^{1}(0,T;L^{2}(\Omega))}+2||(u-\mathbf{P}u)_{t}||_{L^{\infty}(0,T;L^{2}(\Omega))})$
$\displaystyle+C_{b}T|||e|||\cdot|||(u-\mathbf{P}u)_{t}|||$
$\displaystyle+\left|\int^{T}_{0}r_{h}(u,(\mathbf{P}u-u_{h})_{t})dt\right|.$

Since the inequality (3.19) holds for any $s\in[0,T]$, taking the maximum on
$[0,T]$ will result in

$\displaystyle||e_{t}||^{2}_{L^{\infty}(0,T;L^{2}(\Omega))}+C_{s}||e||^{2}_{L^{\infty}(0,T;{\bf
V}(h))}\leq||e_{t}(0)||^{2}+C_{b}|||e(0)|||^{2}+T_{1}+T_{2}+T_{3}$ (3.20)

where the short-hand notation $||e||_{L^{\infty}(0,T;{\bf
V}(h))}:=\sup_{t\in[0,T]}|||e|||$ is introduced, and

$\displaystyle
T_{1}=2||e_{t}||_{L^{\infty}(0,T;L^{2}(\Omega))}(||(u-\mathbf{P}u)_{tt}||_{L^{1}(0,T;L^{2}(\Omega))}+2||(u-\mathbf{P}u)_{t}||_{L^{\infty}(0,T;L^{2}(\Omega))})$
(3.21) $\displaystyle T_{2}=2C_{b}T|||e|||\cdot|||(u-\mathbf{P}u)_{t}|||$
$\displaystyle
T_{3}=2\left|\int^{T}_{0}r_{h}(u,(\mathbf{P}u-u_{h})_{t})dt\right|.$

Using the geometric-arithmetic mean inequality, and Lemma 3.2 in [21], we
conclude

$\displaystyle T_{1}$
$\displaystyle\leq\frac{1}{2}||e_{t}||^{2}_{L^{\infty}(0,T;L^{2}(\Omega))}+2(||(u-\mathbf{P}u)_{tt}||_{L^{1}(0,T;L^{2}(\Omega))}+2||(u-\mathbf{P}u)_{t}||_{L^{\infty}(0,T;L^{2}(\Omega))})^{2}$
(3.22)
$\displaystyle\leq\frac{1}{2}||e_{t}||^{2}_{L^{\infty}(0,T;L^{2}(\Omega))}+4||(u-\mathbf{P}u)_{tt}||_{L^{1}(0,T;L^{2}(\Omega))}^{2}+16||(u-\mathbf{P}u)_{t}||_{L^{\infty}(0,T;L^{2}(\Omega))}^{2}$
$\displaystyle\leq\frac{1}{2}||e_{t}||^{2}_{L^{\infty}(0,T;L^{2}(\Omega))}+C\left|\log_{2}h_{N}\right|^{2d}h_{N}^{2q}(|u_{tt}|^{2}_{L^{\infty}(0,T;\mathcal{H}^{q}(\Omega))}+h_{N}^{2}|u_{t}|^{2}_{L^{2}(0,T;\mathcal{H}^{q+1}(\Omega))}).$

Similarly

$\displaystyle T_{2}$
$\displaystyle\leq\frac{1}{4}C_{s}|||e|||^{2}+4\frac{C^{2}_{b}}{C_{s}}T^{2}|||(u-\mathbf{P}u)_{t}|||^{2}$
(3.23) $\displaystyle\leq\frac{1}{4}C_{s}||e||_{L^{\infty}(0,T;{\bf
V}(h))}^{2}+CT^{2}|\log_{2}h_{N}|^{2d}h_{N}^{2q+2}|u_{t}|^{2}_{L^{2}(0,T;\mathcal{H}^{q+1}(\Omega))}.$

We then start to bound the term $T_{3}$. From (3.15), we can derive

$\displaystyle|r_{h}(u,v)|$
$\displaystyle=|\sum_{e\in\Gamma}\int_{e}[v]\cdot\\{c^{2}\nabla
u-c^{2}\mathbf{P}(\nabla u)\\}ds|$ (3.24)
$\displaystyle\leq(\sum_{e\in\Gamma}\int_{e}\frac{\sigma}{h_{N}}[v]^{2}ds)^{\frac{1}{2}}(\sum_{e\in\Gamma}\int_{e}\frac{h_{N}}{\sigma}|c^{2}\nabla
u-c^{2}\mathbf{P}(\nabla u)|^{2}ds)^{\frac{1}{2}}$ $\displaystyle\leq
C|||v|||(\sum_{K\in\Omega_{N}}h_{N}||\nabla u-\mathbf{P}(\nabla u)||_{\partial
K}^{2})^{\frac{1}{2}},$

with $h_{N}=\frac{1}{2^{N}}$, using trace inequality and Lemma 3.2 in [21], we
have

$\displaystyle|r_{h}(u,v)|\leq C||v||_{L^{\infty}(0,T;{\bf
V}(h))}\cdot\left|\log_{2}h_{N}\right|^{d}h_{N}^{q}|u|_{L^{\infty}(0,T;\mathcal{H}^{q+1}(\Omega))}.$
(3.25)

Therefore,

$\displaystyle\left|\int^{T}_{0}r_{h}(u,v_{t})\right|$
$\displaystyle=\left|-\int^{T}_{0}r_{h}(u_{t},v)dt+r_{h}(u,v)|^{t=T}_{t=0}\right|$
(3.26) $\displaystyle\leq CT||v||_{L^{\infty}(0,T;{\bf
V}(h))}\left|\log_{2}h_{N}\right|^{d}h_{N}^{q}|u_{t}|_{L^{\infty}(0,T;\mathcal{H}^{q+1}(\Omega))}$
$\displaystyle+2C||v||_{L^{\infty}(0,T;{\bf
V}(h))}\left|\log_{2}h_{N}\right|^{d}h_{N}^{q}|u|_{L^{\infty}(0,T;\mathcal{H}^{q+1}(\Omega))}.$

Denote
$\mathcal{R}=T|u_{t}|_{L^{\infty}(0,T;\mathcal{H}^{q+1}(\Omega))}+2|u|_{L^{\infty}(0,T;\mathcal{H}^{q+1}(\Omega))}$,
we will have

$\displaystyle T_{3}$ $\displaystyle\leq
2C\mathcal{R}\left|\log_{2}h_{N}\right|^{d}h_{N}^{q}||\mathbf{P}u-u_{h}||_{L^{\infty}(0,T;{\bf
V}(h))}$ (3.27) $\displaystyle\leq
2C\mathcal{R}\left|\log_{2}h_{N}\right|^{d}h_{N}^{q}\left[||e||_{L^{\infty}(0,T;{\bf
V}(h))}+||u-\mathbf{P}u||_{L^{\infty}(0,T;{\bf V}(h))}\right]$
$\displaystyle\leq\frac{1}{4}C_{s}||e||^{2}_{L^{\infty}(0,T;{\bf
V}(h))}+C\left|\log_{2}h_{N}\right|^{2d}h_{N}^{2q}\left[|u|^{2}_{L^{\infty}(0,T;H^{q+1}(\Omega))}+\mathcal{R}^{2}\right].$

Together with (3.20) and the estimates for $T_{1}$, $T_{2}$ and $T_{3}$, we
arrive at the estimate

$\displaystyle\frac{1}{2}||e_{t}||^{2}_{\infty}$
$\displaystyle+\frac{1}{2}C_{s}\sup_{t\in[0,T]}|||e|||^{2}\leq||e_{t}(0)||^{2}+C|||e(0)|||^{2}$
(3.28)
$\displaystyle+C\left|\log_{2}h_{N}\right|^{2d}h_{N}^{2q}(|u_{tt}|^{2}_{L^{\infty}(0,T;\mathcal{H}^{q}(\Omega))}+T^{2}|u_{t}|^{2}_{L^{\infty}(0,T;\mathcal{H}^{q+1}(\Omega))}+|u|^{2}_{L^{\infty}(0,T;\mathcal{H}^{q+1}(\Omega))}),$

and this completes the proof. ∎

###  3.2 Semi-discrete scheme with multiresolution interpolation

To treat variable coefficient case, we follow the idea in [34, 27] and
interpolate the functions $c^{2}u_{h}$ and $c^{2}\nabla u_{h}$ (or
$(c^{2})^{-}\nabla u_{h}$ and $(c^{2})^{+}\nabla u_{h}$ in the case when
$c^{2}(\mathbf{x})$ contains discontinuity on the cell interfaces of
$\Omega_{N}$) by using the multiresolution Lagrange interpolation discussed in
Section 2.2. For simplicity of discussion, we only focus on the homogeneous
Dirichlet boundary condition with no source term. However, similar results can
be established for mixed boundary conditions and also with source terms.

We first assume $c=c(\mathbf{x})$ is continuous. In this case, we can
reformulate (3.1) into an equivalent form

$\displaystyle B(u_{h},v)$ $\displaystyle=\int_{\Omega}c^{2}\nabla
u_{h}\cdot\nabla v\ d\textbf{x}-\sum_{e\in\Gamma}\int_{e}\\{c^{2}\nabla
u_{h}\\}\cdot[v]\ ds-\sum_{e\in\Gamma}\int_{e}\\{\nabla v\\}\cdot[c^{2}u_{h}]\
ds$
$\displaystyle+\sum_{e\in\Gamma}\frac{\sigma}{h_{N}}\int_{e}[u_{h}]\cdot[v]\
ds$ (3.29)

then the scheme is implemented by the modified operator with interpolation

$\displaystyle\tilde{B}(u_{h},v)$
$\displaystyle=\int_{\Omega}\mathcal{I}_{h}(c^{2}\nabla u_{h})\cdot\nabla v\
d\textbf{x}-\sum_{e\in\Gamma}\int_{e}\\{\mathcal{I}_{h}(c^{2}\nabla
u_{h})\\}\cdot[v]\ ds-\sum_{e\in\Gamma}\int_{e}\\{\nabla
v\\}\cdot[\mathcal{I}_{h}(c^{2}u_{h})]\ ds$
$\displaystyle+\sum_{e\in\Gamma}\frac{\sigma}{h_{N}}\int_{e}[u_{h}]\cdot[v]\
ds$ (3.30)

Here, $\mathcal{I}_{h}(\cdot)$ denote the interpolation operator defined in
Section 2.2 with interpolation parameter $M$ to be specified later.

If $c=c(\mathbf{x})$ is discontinuous along the cell interface, then some
special care has to be taken for the third term
$\sum_{e\in\Gamma}\int_{e}\\{c^{2}\nabla v\\}\cdot[u_{h}]\ ds$. We first
reformulate it into another form:

$\displaystyle\\{c^{2}\nabla v\\}\cdot[u_{h}]$ $\displaystyle={}$
$\displaystyle\frac{1}{2}\left((c^{2}\nabla v)^{-}+(c^{2}\nabla
v)^{+}\right)\cdot(u_{h}^{-}\textbf{n}^{-}+u_{h}^{+}\textbf{n}^{+})$
$\displaystyle={}$ $\displaystyle\frac{1}{2}\left((c^{2}\nabla v)^{-}\cdot
u_{h}^{-}\textbf{n}^{-}+(c^{2}\nabla v)^{-}\cdot
u_{h}^{+}\textbf{n}^{+}+(c^{2}\nabla v)^{+}\cdot
u_{h}^{-}\textbf{n}^{-}+(c^{2}\nabla v)^{+}\cdot
u_{h}^{+}\textbf{n}^{+}\right)$ $\displaystyle={}$
$\displaystyle\frac{1}{2}\left((c^{2})^{-}u_{h}^{-}\textbf{n}^{-}\cdot(\nabla
v)^{-}+(c^{2})^{-}u_{h}^{+}\textbf{n}^{+}\cdot(\nabla
v)^{-}+(c^{2})^{+}u_{h}^{-}\textbf{n}^{-}\cdot(\nabla
v)^{+}+(c^{2})^{+}u_{h}^{+}\textbf{n}^{+}\cdot(\nabla v)^{+}\right)$
$\displaystyle={}$
$\displaystyle\frac{1}{2}((c^{2})^{-}u_{h}^{-}\textbf{n}^{-}+(c^{2})^{-}u_{h}^{+}\textbf{n}^{+})\cdot(\nabla
v)^{-}+\frac{1}{2}((c^{2})^{+}u_{h}^{-}\textbf{n}^{-}+(c^{2})^{+}u_{h}^{+}\textbf{n}^{+})\cdot(\nabla
v)^{+}$ $\displaystyle={}$
$\displaystyle\frac{1}{2}[(c^{2})^{-}u_{h}]\cdot(\nabla
v)^{-}+\frac{1}{2}[(c^{2})^{+}u_{h}]\cdot(\nabla v)^{+}.$

Here
$[(c^{2})^{-}u_{h}]:=((c^{2})^{-}u_{h}^{-}\textbf{n}^{-}+(c^{2})^{-}u_{h}^{+}\textbf{n}^{+})$
and
$[(c^{2})^{+}u_{h}]:=((c^{2})^{+}u_{h}^{-}\textbf{n}^{-}+(c^{2})^{+}u_{h}^{+}\textbf{n}^{+})$.
Now the bilinear form (3.1) is rewriten into

$\displaystyle B(u_{h},v)$ $\displaystyle=\int_{\Omega}c^{2}\nabla
u_{h}\cdot\nabla v\ d\textbf{x}-\sum_{e\in\Gamma}\int_{e}\\{c^{2}\nabla
u_{h}\\}\cdot[v]\ ds$
$\displaystyle-\sum_{e\in\Gamma}\int_{e}\left(\frac{1}{2}[(c^{2})^{-}u_{h}]\cdot(\nabla
v)^{-}+\frac{1}{2}[(c^{2})^{+}u_{h}]\cdot(\nabla v)^{+}\right)\
ds+\sum_{e\in\Gamma}\frac{\sigma}{h_{N}}\int_{e}[u_{h}]\cdot[v]\ ds$

and then the interpolation operator is performed on $c^{2}\nabla u_{h}$,
$(c^{2})^{-}u_{h}$ and also $(c^{2})^{+}u_{h},$ which gives:

$\displaystyle\tilde{B}(u_{h},v)$
$\displaystyle=\int_{\Omega}\mathcal{I}_{h}(c^{2}\nabla u_{h})\cdot\nabla v\
d\textbf{x}-\sum_{e\in\Gamma}\int_{e}\\{\mathcal{I}_{h}(c^{2}\nabla
u_{h})\\}\cdot[v]\ ds$ (3.31)
$\displaystyle-\sum_{e\in\Gamma}\int_{e}\left(\frac{1}{2}[\mathcal{I}_{h}((c^{2})^{-}u_{h})]\cdot(\nabla
v)^{-}+\frac{1}{2}[\mathcal{I}_{h}((c^{2})^{+}u_{h})]\cdot(\nabla
v)^{+}\right)\
ds+\sum_{e\in\Gamma}\frac{\sigma}{h_{N}}\int_{e}[u_{h}]\cdot[v]\ ds.$

Following [14], we can now write the DG scheme with interpolation (3.31) into
the semi-discrete form as

$\frac{d^{2}u_{h}}{dt^{2}}=L_{h}(u_{h}),$ (3.32)

where $L_{h}(u)$ is an operator onto ${\bf V}$ which is a discrete
approximation of $-\nabla\cdot(c^{2}(\mathbf{x})\nabla u)$ and satisfies

$\displaystyle\sum_{K\in{\Omega_{N}}}\int_{K}L_{h}(u_{h})v_{h}\ d\textbf{x}$
$\displaystyle=-\int_{\Omega}\mathcal{I}_{h}(c^{2}\nabla u_{h})\cdot\nabla
v_{h}\ d\textbf{x}+\sum_{e\in\Gamma}\int_{e}\\{\mathcal{I}_{h}(c^{2}\nabla
u_{h})\\}\cdot[v_{h}]\ ds$
$\displaystyle+\sum_{e\in\Gamma}\int_{e}\left(\frac{1}{2}[\mathcal{I}_{h}((c^{2})^{-}u_{h})]\cdot(\nabla
v_{h})^{-}+\frac{1}{2}[\mathcal{I}_{h}((c^{2})^{+}u_{h})]\cdot(\nabla
v_{h})^{+}\right)\ ds$
$\displaystyle-\sum_{e\in\Gamma}\frac{\sigma}{h_{N}}\int_{e}[u_{h}]\cdot[v_{h}]\
ds$ (3.33)

for any $v_{h}\in{\bf V}$.

To preserve the accuracy of the original DG scheme, interpolation operator
$\mathcal{I}_{h}(\cdot)$ needs to reach certain accuracy. Using similar
techniques as in [14, 28], we have the following proposition on local
truncation error of the sparse grid method with ${\bf V}=\hat{{\bf
V}}^{k}_{N}.$ We only discuss the case when $c(\mathbf{x})$ is discontinuous,
since the similar approach can be applied when $c(\mathbf{x})$ is continuous.

######  Proposition 3.5 (Local truncation error analysis).

If the interpolation operator $\mathcal{I}_{h}$ in (3.31) has the accuracy of
order $\left|\log_{2}h_{N}\right|^{d}h_{N}^{k+3}$ for sufficiently smooth
functions, then the local truncation error of the semi-discrete DG scheme with
interpolation (3.31) is of order $\left|\log_{2}h_{N}\right|^{d}h_{N}^{k+1}$.
To be more precise, for sufficiently smooth function $u$, the sparse grid DG
method with interpolation (3.31) has the truncation error:

$\left\|L_{h}(u)+\nabla\cdot(c^{2}(\mathbf{x})\nabla
u)\right\|_{L^{2}(\Omega)}\leq C\left|\log_{2}h_{N}\right|^{d}h_{N}^{k+1}.$
(3.34)

Here, we use $C$ to denote any generic constant that may depend on the
solution $u$ and $c(\mathbf{x}),$ but does not depend on $N.$

###### Proof.

We denote the standard $L^{2}$ projection operator onto the sparse grid DG
finite element space by $\mathbf{P}$, then

$\left\|L_{h}(u)+\nabla\cdot(c^{2}(\mathbf{x})\nabla u)\right\|\leq
e_{1}+e_{2},$ (3.35)

where

$e_{1}:=\left\|L_{h}(u)+\mathbf{P}(\nabla\cdot(c^{2}(\mathbf{x})\nabla
u))\right\|,$

and

$e_{2}:=\left\|\mathbf{P}(\nabla\cdot(c^{2}(\mathbf{x})\nabla
u)-\nabla\cdot(c^{2}(\mathbf{x})\nabla u)\right\|$

The estimate for $e_{2}$ can be obtained by projection properties [21]:

$e_{2}\leq C\left|\log_{2}h_{N}\right|^{d}h_{N}^{k+1}.$ (3.36)

To estimate $e_{1},$ we consider any test function $v_{h}$ in DG space, and
obtain

$\displaystyle\sum_{K\in{\Omega_{N}}}\int_{K}(L_{h}(u)+\mathbf{P}(\nabla\cdot(c^{2}(\mathbf{x})\nabla
u))v_{h}\
d\textbf{x}=\sum_{K\in{\Omega_{N}}}\int_{K}(L_{h}(u)+\nabla\cdot(c^{2}(\mathbf{x})\nabla
u))v_{h}\ d\textbf{x}$ $\displaystyle={}$
$\displaystyle-\int_{\Omega}\mathcal{I}_{h}(c^{2}\nabla u)\cdot\nabla v_{h}\
d\textbf{x}+\sum_{e\in\Gamma}\int_{e}\\{\mathcal{I}_{h}(c^{2}\nabla
u)\\}\cdot[v_{h}]\ ds$
$\displaystyle+\sum_{e\in\Gamma}\int_{e}\left(\frac{1}{2}[\mathcal{I}_{h}((c^{2})^{-}u)]\cdot(\nabla
v_{h})^{-}+\frac{1}{2}[\mathcal{I}_{h}((c^{2})^{+}u)]\cdot(\nabla
v_{h})^{+}\right)\
ds-\sum_{e\in\Gamma}\frac{\sigma}{h_{N}}\int_{e}[u]\cdot[v_{h}]\ ds$
$\displaystyle+\int_{\Omega}c^{2}\nabla u\cdot\nabla v_{h}\
d\textbf{x}-\sum_{e\in\Gamma}\int_{e}\\{c^{2}\nabla u\\}\cdot[v_{h}]\ ds$
$\displaystyle-\sum_{e\in\Gamma}\int_{e}\left(\frac{1}{2}[(c^{2})^{-}u]\cdot(\nabla
v_{h})^{-}+\frac{1}{2}[(c^{2})^{+}u]\cdot(\nabla v_{h})^{+}\right)\
ds+\sum_{e\in\Gamma}\frac{\sigma}{h_{N}}\int_{e}[u]\cdot[v_{h}]\ ds$
$\displaystyle={}$
$\displaystyle-\int_{\Omega}\left(\mathcal{I}_{h}(c^{2}\nabla u)-c^{2}\nabla
u\right)\cdot\nabla v_{h}\
d\textbf{x}+\sum_{e\in\Gamma}\int_{e}\\{\mathcal{I}_{h}(c^{2}\nabla
u)-c^{2}\nabla u\\}\cdot[v_{h}]\ ds$
$\displaystyle+\sum_{e\in\Gamma}\int_{e}\left(\frac{1}{2}[\mathcal{I}_{h}((c^{2})^{-}u)-(c^{2})^{-}u]\cdot(\nabla
v_{h})^{-}+\frac{1}{2}[\mathcal{I}_{h}((c^{2})^{+}u)-(c^{2})^{+}u]\cdot(\nabla
v_{h})^{+}\right)\ ds$ $\displaystyle\leq{}$
$\displaystyle\left\|\mathcal{I}_{h}[c^{2}\nabla u]-c^{2}\nabla
u\right\|\left\|\nabla v_{h}\right\|+\left\|\mathcal{I}_{h}[c^{2}\nabla
u]-c^{2}\nabla
u\right\|_{L^{2}(\Gamma_{h})}\left\|v_{h}\right\|_{L^{2}(\Gamma_{h})}$
$\displaystyle+\left\|\mathcal{I}_{h}[c^{2}\nabla u]-c^{2}\nabla
u\right\|_{L^{2}(\Gamma_{h})}\left\|\nabla v_{h}\right\|_{L^{2}(\Gamma_{h})}$
$\displaystyle\leq{}$ $\displaystyle
C\left|\log_{2}h_{N}\right|^{d}h_{N}^{k+3}h_{N}^{-1}\left\|v_{h}\right\|+Ch_{N}^{-\frac{1}{2}}\left|\log_{2}h_{N}\right|^{d}h_{N}^{k+3}h_{N}^{-\frac{1}{2}}\left\|v_{h}\right\|+Ch_{N}^{-\frac{1}{2}}\left|\log_{2}h_{N}\right|^{d}h_{N}^{k+3}h_{N}^{-\frac{3}{2}}\left\|v_{h}\right\|$
$\displaystyle={}$ $\displaystyle
C\left|\log_{2}h_{N}\right|^{d}h_{N}^{k+1}\left\|v_{h}\right\|.$

Here, we have used the multiplicative trace inequality and the inverse
inequality, see e.g. Lemma 2.1 and Lemma 2.3 in [28]. By taking
$v_{h}=(L_{h}(u)+\mathbf{P}(\nabla\cdot(c^{2}(\mathbf{x})\nabla u))$ in the
inequality above, we have

$e_{1}=\left\|L_{h}(u)+\mathbf{P}(\nabla\cdot(c^{2}(\mathbf{x})\nabla
u))\right\|\leq C\left|\log_{2}h_{N}\right|^{d}h_{N}^{k+1}.$ (3.37)

Combining (3.37) and (3.36), we have the estimate for the truncation error
(3.34). ∎

######  Remark 3.6.

The proposition above indicates that, to preserve the order of the original
scheme, we should use $M\geq k+2.$ For example, if we take piecewise linear
polynomials for the DG space, then it is required to apply cubic interpolation
operator to treat the nonlinear terms. From our numerical tests, however, this
seems that it is not a necessary condition. To reach the desired convergence
rate, one only needs to take $M\geq k$.

In Proposition 3.5, we only estimate the local truncation error, and this is
far from a rigorous error estimate that takes into account stability. Unlike
the scheme with the symmetric bilinear form $B(u_{h},v)$ as in Theorem 3.1,
the symmetry is lost in the interpolated bilinear form $\tilde{B}(u_{h},v).$
Hence, energy stability is not automatic. In numerical experiments, we observe
that the sparse grid DG method with Lagrange interpolation with only inner
interface points is unstable for polynomials of high degrees (see the
numerical results in Table 4.5 in Section 4). With the interpolation points at
the _interface_ , the sparse grid DG scheme is stable and yields satisfactory
convergence rate (see Table 4.6 in Section 4).

###  3.3 Time stepping and adaptivity

For time discretizations, we first write the second order semi-discrete scheme
(3.31)

$(u_{h})_{tt}=L_{h}(u_{h})$ (3.38)

into a first order system

$\displaystyle(u_{h})_{t}$ $\displaystyle=w_{h},$ $\displaystyle(w_{h})_{t}$
$\displaystyle=L_{h}(u_{h}),$

and then appy the standard Runge-Kutta scheme. The reason why we use the one-
step RK method instead of the multistep method is that the maximum allowed
time step size from the CFL restriction may change with the adaptive mesh in
different time steps. This would result in additional computational cost in
extrapolation or interpolation between different time steps for the multistep
methods.

The adaptive scheme uses the procedure developed in [8, 22] to determine the
space ${\bf V}$ that dynamically evolves over time. The method is very similar
to those in [8, 22], and the details are omitted for brevity. The main
difference is that we need to keep track of two sets of basis functions
corresponding to the same adaptive space are involved [27]. Another difference
is that the refinement and the coarsening criteria are determined by the
$L^{2}$ norms of both $u_{h}$ and $w_{h}$, which are both important for
predicting solution profiles for wave equations. There are some cases which
start with a zero displacement $u$ but a non-zero velocity $u_{t}$. If we only
take the norms of $u_{h}$ as an indicator, the adaptive procedure will result
in poor resolutions. Only by considering the norms of both $u_{h}$ and
$w_{h}$, one can capture the profiles well.

###  3.4 Fast algorithms

We now describe the fast matrix-vector multiplication algorithm, which is
essential for efficient implementation of our schemes. Because the
multiwavelet bases are global, the evaluation of the residual yields denser
matrix than those obtained by standard local bases. Efficient implementations
are therefore essential to ensure that the computational cost is on par with
element-wise implementation of traditional DG schemes. This issue has been
also discussed in our work for conservation laws [27], which extends the fast
matrix-vector multiplication in [34, 42] to adaptive index set.

Following [34, 27], we consider matrix-vector multiplication in multi-
dimensions in an abstract framework.

$f_{\bm{n}}=\sum_{H(\bm{n}^{\prime})\leq
0}f^{\prime}_{\bm{n}^{\prime}}t_{n_{1}^{\prime},n_{1}}^{(1)}t_{n_{2}^{\prime},n_{2}}^{(2)}\cdots
t_{n_{d}^{\prime},n_{d}}^{(d)},\quad H(\bm{n})\leq 0,$ (3.39)

where $\bm{n}=({n}_{1},{n}_{2},\dots,{n}_{d})$ and
$\bm{n}^{\prime}=({n}_{1}^{\prime},{n}_{2}^{\prime},\dots,{n}_{d}^{\prime})$
can be thought of as the level of the mesh, and
$t_{n_{1}^{\prime},n_{1}}^{(i)}=T^{(i)}_{n_{1}^{\prime},n_{1}}$ represents the
calculations in the $i$-th dimension. It is assumed that the constraint
function
$H=H(\bm{n}^{\prime})=H({n}_{1}^{\prime},{n}_{2}^{\prime},\dots,{n}_{d}^{\prime})$
is non-decreasing with respect to each variable. This holds true for sparse
grid (by taking $H(\bm{n}^{\prime})=|\bm{n}^{\prime}|_{1}$) and also for
adaptive multiresolution method.

One can compute the sum (3.39) dimension-by-dimension, i.e. we first perform
the transformation in the $x_{1}$ dimension:

$g^{(1)}_{(n_{1},n_{2}^{\prime},\dots,n_{d}^{\prime})}=\sum_{H(n_{1}^{\prime},n_{2}^{\prime},\dots,n_{d}^{\prime})\leq
0}f^{\prime}_{(n_{1}^{\prime},n_{2}^{\prime},\dots,n_{d}^{\prime})}t_{n_{1}^{\prime},n_{1}}^{(1)},$
(3.40)

and then in the $x_{2}$ dimension:

$g^{(2)}_{(n_{1},n_{2},n_{3}^{\prime}\dots,n_{d}^{\prime})}=\sum_{H(n_{1},n_{2}^{\prime},\dots,n_{d}^{\prime})\leq
0}g^{(1)}_{(n_{1},n_{2}^{\prime},\dots,n_{d}^{\prime})}t_{n_{2}^{\prime},n_{2}}^{(2)},$
(3.41)

and all the way up to $x_{d}$ dimension:

$f_{(n_{1},n_{2},n_{3}\dots,n_{d})}=\sum_{H(n_{1},n_{2},\dots,n_{d-1},n_{d}^{\prime})\leq
0}g^{(d-1)}_{(n_{1},n_{2},\dots,n_{d-1},n_{d}^{\prime})}t_{n_{d}^{\prime},n_{d}}^{(d)}.$
(3.42)

It can be proved that (3.40)-(3.42) is equivalent to the original summation
(3.39), if assuming that, for some integer $1\leq k\leq d$, $T^{(i)}$ for
$i=1,\dots,k-1$ are strictly lower triangular and $T^{(i)}$ for
$i=k+1,\dots,d$ are upper triangular (or $T^{(i)}$ for $i=1,\dots,k-1$ are
lower triangular and $T^{(i)}$ for $i=k+1,\dots,d$ are strictly upper
triangular) [34]. Here, $T^{(i)}$ denotes the $i$-th transformation matrix.
When such properties for $T^{(i)}$ matrices are not true, one can perform
$L+U$ split and (3.39) becomes:

$f_{\bm{n}}=\sum_{H(\bm{n}^{\prime})\leq
0}f^{\prime}_{\bm{n}^{\prime}}(l_{n_{1}^{\prime},n_{1}}^{(1)}+u_{n_{1}^{\prime},n_{1}}^{(1)})(l_{n_{1}^{\prime},n_{1}}^{(1)}+u_{n_{1}^{\prime},n_{1}}^{(1)})\cdots(l_{n_{d-1}^{\prime},n_{d-1}}^{(d-1)}+u_{n_{d-1}^{\prime},n_{d-1}}^{(d-1)})t_{n_{d}^{\prime},n_{d}}^{(d)},$
(3.43)

where there are in total $2^{d-1}$ terms that can be computed dimension-by-
dimension. The overall computational cost is $\mathcal{O}(2^{d-1}\cdot
DoF\cdot N)$ if the cost of one-dimensional transform is log-linear, i.e.,
$\mathcal{O}(\mathcal{N}\log\mathcal{N})$ where $\mathcal{N}$ denotes the DoF
in one-dimension [34]. This assumption holds true for our sparse grid DG
scheme.

We apply this fast matrix-vector multiplication in several parts of our
algorithm. We will discuss the details about initialization, which is the
procedure to project the initial value onto the DG finite element space
represented by multiwavelet bases. When the given initial value is separable,
i.e.,

$u(x)=\prod_{i=1}^{d}u_{i}(x_{i}),$ (3.44)

one just need to project each 1D function $u_{i}=u_{i}(x_{i})$ for
$i=1,\dots,d$ onto 1D multiwavelet bases and then we can easily get the
projection of $u(x)$. This approach naturally extends to the case when the
initial value is a summation of separable functions:

$u(x)=\sum_{j=1}^{n}(\prod_{i=1}^{d}u_{j,i}(x_{i})).$ (3.45)

However, if the function is non-separable, direct evaluation of $L^{2}$
projection would result in very large computational cost if using numerical
quadratures in multi-dimensions. Rather, we propose to apply the adaptive
multiresolution interpolation introduced in [37] with appropriate error
tolerance and approximate the solution using a collection of interpolation
basis functions:

$u_{h}(\mathbf{x})=\sum_{\begin{subarray}{c}(\mathbf{l},\mathbf{j})\in G,\\\
\mathbf{1}\leq\mathbf{i}\leq\mathbf{k}+\mathbf{1}\end{subarray}}{b}^{\mathbf{j}}_{\mathbf{i},\mathbf{l}}\psi^{\mathbf{j}}_{\mathbf{i},\mathbf{l}}(\mathbf{x})$
(3.46)

with $G$ the index of all active elements. Next we use the fast matrix-vector
multiplication (3.39) to transform coefficients of interpolation basis
$\\{{b}^{\mathbf{j}}_{\mathbf{i},\mathbf{l}}\\}$ to coefficients of Alpert’s
basis $\\{{c}^{\mathbf{j}}_{\mathbf{i},\mathbf{l}}\\}$:

$u_{h}(\mathbf{x})=\sum_{\begin{subarray}{c}(\mathbf{l},\mathbf{j})\in G,\\\
\mathbf{1}\leq\mathbf{i}\leq\mathbf{k}+\mathbf{1}\end{subarray}}{c}^{\mathbf{j}}_{\mathbf{i},\mathbf{l}}v^{\mathbf{j}}_{\mathbf{i},\mathbf{l}}(\mathbf{x})$
(3.47)

In (3.39), $f_{\bm{n}}$ and $f^{\prime}_{\bm{n}^{\prime}}$ represent
$\\{{b}^{\mathbf{j}}_{\mathbf{i},\mathbf{l}}\\}$ and
$\\{{c}^{\mathbf{j}}_{\mathbf{i},\mathbf{l}}\\}$, respectively. The matrix
$t_{n_{i}^{\prime},n_{i}}^{(i)}$ is the product of 1D interpolation basis and
1D Alpert’s basis.

We apply similar approach in the multiresolution interpolation and the
evaluation of the right hand side of the weak formulation (3.31). We refer
readers to [27] for details.

##  4 Numerical examples

In this section, we perform numerical experiments to validate the performance
of our scheme. We consider 2D or 3D problems with computational domain being
$[0,1]^{d}$ with $d=2,3$. The CFL number is taken to be $0.1$ in 2D and $0.05$
in 3D. The penalty parameter $\sigma$ is taken to be 10 in 2D and 30 in 3D,
unless otherwise stated. For the accuracy test, we check the convergence order
for $P^{k}$ DG with $k=1,2,3$ coupled with RK time discretization. In
particular, for $k=1,2,$ we use the second and the third-order strong
stability preserving Runge-Kutta method [35, 19], and for $k=3,$ we use the
classical RK4 methods. All adaptive calculations are obtained by $k=3$ and RK4
time stepping. In the adaptive scheme, we take $\eta=\epsilon/10$.
$\textrm{DoF=dim}({\bf V}^{k})$ refers to the number of Alperts’ multiwavelets
basis functions in the adaptive grids. The maximum mesh level $N$ is taken to
be 8, unless otherwise stated.

######  Example 4.1 (wave equation with constant coefficient).

In this example, consider the $d$-dimensional wave equation with a constant
coefficient

$u_{tt}=\sum_{i=1}^{d}u_{x_{i}x_{i}}\\\ $ (4.1)

on the domain $[0,1]^{d}$. We take the exact solution to be

$u(\mathbf{x},t)=\sin(a\sqrt{d}\pi t)\prod\limits_{i=1}^{d}\cos(a\pi x_{i})$

with $a$ being a constant and various types of boundary conditions.

  1. (a)

We take $a=2$ and $d=2,3$ with periodic boundary conditions.

  2. (b)

We take $a=1$ and $d=2,3$ and incorporate Dirichlet boundary condition in the
$x_{1}$-direction and Neumann boundary in other directions.

Note that in this example, (3.1) is implemented with no interpolation because
$c$ is a constant. To output the $L^{2}$-error between the numerical solution
$u_{h}$ and the exact solution $u(x),$ we use the fact that

$\displaystyle\int_{\Omega}(u_{h}-u)^{2}dx=\int_{\Omega}u_{h}^{2}dx-2\int_{\Omega}u_{h}udx+\int_{\Omega}u^{2}dx.$

The first term $\int_{\Omega}u_{h}^{2}dx$ can be easily computed with the aid
of the orthonormality of the Alpert’s basis functions. The second term
$\int_{\Omega}u_{h}udx$ can be computed by the same fast approach as the
initial projection, which has been explained in detail in Section 3.4. The
third term $\int_{\Omega}u^{2}dx$ can be computed analytically.

The numerical results obtained by sparse grid DG method are presented in Table
4.1 for case (a) and in Table 4.2 for case (b). For both cases, the
convergence order is slightly bigger than $k$ but smaller than $k+1$, which is
higher than the predicted rate in Theorem 3.4, but similar to the results for
linear transport equation in [21]. The numerical results with adaptive method
are shown in Tables 4.3 and 4.4. Similar to [22], we measure the convergence
rates with respect to DoF: $R_{\textrm{DoF}}$ and $\epsilon:R_{\epsilon}.$ We
can clearly observe the effectiveness of the adaptive algorithm, i.e.
$R_{\epsilon}$ is close to 1. The convergence order $R_{\textrm{DoF}}$ is
bigger than $\frac{k+1}{d},$ which is the rate obtained by an optimally
convergent non-adaptive scheme.

Table 4.1: Example 4.1(a): wave equation with constant coefficients, periodic boundary conditions, sparse grid DG, $k=1,2,3$, $d=2,3$. $t=0.1$. | $N$ | $k=1$ | $N$ | $k=2$ | $N$ |  $k=3$
---|---|---|---|---|---|---
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$d=2$ | 5 | 5.90e-03 | - | 5 | 1.96e-04 | - | 3 | 2.80e-04 | -
6 | 1.69e-03 | 1.81 | 6 | 3.03e-05 | 2.69 | 4 | 1.80e-05 | 3.96
7 | 4.66e-04 | 1.86 | 7 | 4.43e-06 | 2.77 | 5 | 1.48e-06 | 3.61
8 | 1.23e-04 | 1.92 | 8 | 6.21e-07 | 2.83 | 6 | 1.10e-07 | 3.76
| $N$ | $k=1$ | $N$ | $k=2$ | $N$ |  $k=3$
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$d=3$ | 5 | 1.58e-02 | - | 5 | 7.38e-04 | - | 3 | 4.50e-04 | -
6 | 8.66e-03 | 0.87 | 6 | 1.68e-04 | 2.14 | 4 | 8.02e-05 | 2.49
7 | 2.42e-03 | 1.84 | 7 | 3.03e-05 | 2.47 | 5 | 5.24e-06 | 3.94
8 | 8.41e-04 | 1.53 | 8 | 5.30e-06 | 2.51 | 6 | 5.30e-07 | 3.30
Table 4.2: Example 4.1(b): wave equation with constant coefficients, Dirichlet and Neumann boundary conditions, sparse grid DG, $k=1,2,3$, $d=2,3$. $t=0.1$. | $N$ | $k=1$ | $N$ | $k=2$ | $N$ |  $k=3$
---|---|---|---|---|---|---
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$d=2$ | 3 | 3.30e-03 | - | 3 | 1.21e-04 | - | 1 | 4.51e-04 | -
4 | 1.21e-03 | 1.44 | 4 | 1.79e-05 | 2.75 | 2 | 6.34e-05 | 2.83
5 | 2.94e-04 | 2.04 | 5 | 2.43e-06 | 2.89 | 3 | 8.30e-06 | 2.93
6 | 8.15e-05 | 1.85 | 6 | 3.41e-07 | 2.83 | 4 | 7.64e-07 | 3.44
| $N$ | $k=1$ | $N$ | $k=2$ | $N$ |  $k=3$
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$d=3$ | 3 | 2.15e-02 | - | 3 | 2.14e-04 | - | 1 | 4.44e-04 | -
4 | 7.06e-03 | 1.61 | 4 | 3.39e-05 | 2.66 | 2 | 3.66e-05 | 3.60
5 | 2.04e-03 | 1.79 | 5 | 5.13e-06 | 2.73 | 3 | 2.39e-06 | 3.94
6 | 5.27e-04 | 1.95 | 6 | 1.08e-06 | 2.24 | 4 | 1.23e-07 | 4.28
Table 4.3: Example 4.1(a): wave equation with constant coefficients, periodic boundary conditions, adaptive sparse grid DG, $k=3$. | $\epsilon$ | DoF |  $L^{2}$-error | $R_{\textrm{DoF}}$ | $R_{\epsilon}$
---|---|---|---|---|---
$d=2$ | 1e-1 | 128 | 1.25e-3 | - | -
1e-2 | 320 | 2.80e-4 | 1.63 | 0.65
1e-3 | 1088 | 2.61e-5 | 1.94 | 1.03
1e-4 | 1536 | 4.02e-6 | 5.43 | 0.81
$d=3$ | 1e-1 | 896 | 2.16e-3 | - | -
1e-2 | 2432 | 4.74e-4 | 1.52 | 0.66
1e-3 | 5888 | 8.30e-5 | 1.97 | 0.76
1e-4 | 28160 | 6.69e-6 | 1.61 | 1.09
Table 4.4: Example 4.1(b): wave equation with constant coefficients, Dirichlet and Neumann boundary conditions, adaptive sparse grid DG, $k=3$. | $\epsilon$ | DoF |  $L^{2}$-error | $R_{\textrm{DoF}}$ | $R_{\epsilon}$
---|---|---|---|---|---
$d=2$ | 1e-1 | 32 | 1.11e-3 | - | -
1e-2 | 112 | 6.35e-5 | 2.29 | 1.24
1e-3 | 208 | 1.70e-5 | 2.13 | 0.57
1e-4 | 384 | 2.35e-6 | 3.22 | 0.86
$d=3$ | 1e-1 | 64 | 1.29e-3 | - | -
1e-2 | 640 | 3.71e-5 | 1.54 | 1.54
1e-3 | 1280 | 2.41e-5 | 0.62 | 0.19
1e-4 | 2368 | 3.12e-6 | 3.32 | 0.89

######  Example 4.2 (wave equation with smooth variable coefficient).

This example tests wave equation with smooth variable coefficient

$u_{tt}-\nabla\cdot(c^{2}(\mathbf{x})\nabla u)=f,$ (4.2)

on the computational domain $[0,1]^{d}$ with $d=2,3$ and periodic boundary
conditions.

For 2D case, we take

$c^{2}(x_{1},x_{2})=(\cos(2\pi x_{1})\cos(2\pi x_{2})+2)/3,$ (4.3)

and the corresponding source term $f=f(x_{1},x_{2},t)$ such that the exact
solution is

$u=\sin(\pi t)\sin(2\pi x_{1})\cos(2\pi x_{2}).$ (4.4)

For 3D case, we take

$c^{2}(x_{1},x_{2},x_{3})=(\sin(2\pi x_{1})\sin(2\pi x_{2})\cos(2\pi
x_{3})+2)/3,$ (4.5)

and the corresponding source term $f=f(x_{1},x_{2},x_{3},t)$ such that the
exact solution is

$u=\sin(\pi t)\sin(2\pi x_{1})\cos(2\pi x_{2})\cos(2\pi x_{3}).$ (4.6)

This problem needs to invoke the fast interpolation methods to handle the
variable coefficient. We first compare different choices of interpolation
points. We use the inner interpolation points in Table 4.5 for 2D. The
interpolation points and the basis functions are listed in the Appendix. When
$k=1$, the convergence order seems satisfactory. However, for $k=2,3,$ the
scheme is unstable. If we use Lagrange interpolation with the interface
points, one will observe good convergence rate for $M\geq k$, as shown in
Table 4.6. We also find that the error is almost the same for $M=k+1$ and
$M=k+2$, and both much smaller than $M=k$. Therefore, in applications, we
recommend taking $M=k+1$ for accuracy considerations. Notice that this is a
more relaxed condition from what is indicated by the local truncation analysis
Proposition 3.5. We also remark that for nonlinear conservation laws in [27],
Lagrange interpolation is unstable even with interface points, and Hermite
interpolation has to be employed. However, for all numerical examples in this
paper for linear wave equations with variable coefficients, Lagrange
interpolation with interface points yields a stable scheme, and we choose to
use this instead of Hermite interpolation due to its easier implementation.

For 3D cases, to save space, we only show numerical results with interface
interpolation points in Table 4.7, in which good convergence rate is also
observed. The result using adaptive method with $k=3$ and $M=4$ are presented
in Table 4.8 for both 2D and 3D, and the conclusions are similar to the
constant coefficient case.

Table 4.5: Example 4.2: wave equation with smooth variable coefficients in 2D, sparse grid DG, Lagrange interpolation with inner interpolation points, $k=1,2,3$. $t=0.1$. | $N$ | $M=1$ | $M=2$ |  $M=3$
---|---|---|---|---
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$k=1$ | 3 | 2.52e-02 | - | 2.55e-02 | - | 2.52e-02 | -
4 | 1.68e-02 | 0.59 | 1.63e-02 | 0.64 | 1.63e-02 | 0.63
5 | 3.67e-03 | 2.19 | 3.37e-03 | 2.28 | 3.36e-03 | 2.28
6 | 9.43e-04 | 1.96 | 9.62e-04 | 1.81 | 8.33e-04 | 2.01
| $N$ | $M=2$ | $M=3$ |  $M=4$
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$k=2$ | 3 | 1.64e-02 | - | 3.95e-03 | - | 4.73e-02 | -
4 | 1.03e-02 | 0.67 | 7.89e-04 | 2.32 | 1.46e-01 | -1.63
5 | 3.57e-03 | 1.53 | 7.80e-04 | 0.02 | 3.17e+02 | -11.08
6 | 1.13e-02 | -1.66 | 1.60e-01 | -7.68 | 6.27e+10 | -27.56
| $N$ | $M=3$ | $M=4$ |  $M=5$
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$k=3$ | 3 | 5.28e-03 | - | 2.14e+03 | - | 5.13e+03 | -
4 | 5.55e-02 | -3.40 | 1.11e+10 | -22.31 | 2.23e+09 | -18.73
5 | 1.93e+03 | -15.09 | 1.39e+24 | -46.83 | 4.78e+21 | -40.97
6 | 8.15e+20 | -58.55 | 9.83e+56 | -109.12 | 1.40e+51 | -97.88
Table 4.6: Example 4.2: wave equation with smooth variable coefficients in 2D, sparse grid DG, Lagrange interpolation with interface interpolation points, $k=1,2,3$. $t=0.1$. | $N$ | $M=2$ | $M=3$ |  $M=4$
---|---|---|---|---
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$k=1$ | 3 | 2.52e-02 | - | 2.52e-02 | - | 2.52e-02 | -
4 | 1.65e-02 | 0.61 | 1.63e-02 | 0.63 | 1.63e-02 | 0.63
5 | 3.52e-03 | 2.23 | 3.36e-03 | 2.28 | 3.36e-03 | 2.28
6 | 9.52e-04 | 1.89 | 8.30e-04 | 2.02 | 8.27e-04 | 2.02
| $N$ | $M=2$ | $M=3$ |  $M=4$
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$k=2$ | 3 | 2.69e-03 | - | 2.08e-03 | - | 2.08e-03 | -
4 | 5.24e-04 | 2.36 | 4.38e-04 | 2.25 | 4.37e-04 | 2.25
5 | 1.25e-04 | 2.07 | 7.58e-05 | 2.53 | 7.58e-05 | 2.53
6 | 1.64e-05 | 2.93 | 1.16e-05 | 2.71 | 1.16e-05 | 2.71
| $N$ | $M=3$ | $M=4$ |  $M=5$
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$k=3$ | 3 | 2.92e-04 | - | 9.28e-05 | - | 8.75e-05 | -
4 | 2.66e-05 | 3.46 | 1.05e-05 | 3.15 | 1.03e-05 | 3.09
5 | 3.04e-06 | 3.13 | 7.80e-07 | 3.74 | 7.68e-07 | 3.74
6 | 1.83e-07 | 4.05 | 5.10e-08 | 3.94 | 5.03e-08 | 3.93
Table 4.7: Example 4.2: wave equation with smooth variable coefficients in 3D, sparse grid DG, Lagrange interpolation with interface interpolation points, $k=1,2,3$. $t=0.1$. | $N$ | $M=1$ | $M=2$ |  $M=3$
---|---|---|---|---
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$k=1$ | 3 | 1.17e-01 | - | 1.17e-01 | - | 1.17e-01 | -
4 | 2.20e-02 | 2.41 | 2.20e-02 | 2.41 | 2.20e-02 | 2.41
5 | 1.74e-02 | 0.34 | 1.71e-02 | 0.36 | 1.71e-02 | 0.36
6 | 4.65e-03 | 1.90 | 4.52e-03 | 1.92 | 4.51e-03 | 1.92
| $N$ | $M=2$ | $M=3$ |  $M=4$
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$k=2$ | 4 | 2.96e-03 | - | 1.58e-03 | - | 1.58e-03 | -
5 | 7.78e-04 | 1.93 | 3.28e-04 | 2.27 | 3.27e-04 | 2.27
6 | 2.93e-04 | 1.41 | 6.58e-05 | 2.32 | 6.58e-05 | 2.32
7 | 3.88e-05 | 2.92 | 1.15e-05 | 2.52 | 1.15e-05 | 2.52
| $N$ | $M=3$ | $M=4$ |  $M=5$
$L^{2}$-error | order |  $L^{2}$-error | order |  $L^{2}$-error | order
$k=3$ | 3 | 8.96e-04 | - | 3.88e-04 | - | 3.17e-04 | -
4 | 2.05e-04 | 2.13 | 3.58e-05 | 3.44 | 2.19e-05 | 3.85
5 | 4.87e-05 | 2.07 | 3.27e-06 | 3.45 | 3.01e-06 | 2.86
6 | 5.80e-06 | 3.07 | 2.55e-07 | 3.68 | 2.31e-07 | 3.71
Table 4.8: Example 4.2, wave equation with smooth variable coefficient, adaptive sparse grid DG, 2D and 3D. $k=3$, $M=4$, $t=0.1$. | $\epsilon$ | DoF |  $L^{2}$-error | $R_{\textrm{DoF}}$ | $R_{\epsilon}$
---|---|---|---|---|---
$d=2$ | 1e-1 | 96 | 1.66e-3 | - | -
1e-2 | 224 | 3.03e-4 | 2.00 | 0.74
1e-3 | 672 | 2.78e-5 | 2.18 | 1.04
1e-4 | 1088 | 3.17e-6 | 4.50 | 0.94
$d=3$ | 1e-1 | 576 | 2.11e-3 | - | -
1e-2 | 1152 | 5.26e-4 | 2.00 | 0.60
1e-3 | 3584 | 8.73e-5 | 1.58 | 0.78
1e-4 | 8704 | 1.26e-5 | 2.18 | 0.84

######  Example 4.3 (wave equation with discontinuous coefficients).

In this example, we consider wave equation with discontinuous coefficients.
The jump of the coefficient aligns with the cell interface on the fine mesh
$\Omega_{N}.$

For 2D case, the domain $\Omega=[0,1]^{2}$ is composed of two subdomains
$\Omega_{1}=[\frac{1}{4},\frac{3}{4}]\times[0,1]$ and
$\Omega_{2}=\Omega\backslash\Omega_{1}$. The coefficient $c^{2}$ is a constant
in each subdomain:

$\displaystyle c^{2}=\begin{cases}1,\qquad\text{in}\quad\Omega_{1},\\\
{\frac{5}{37}},\qquad\text{in}\quad\Omega_{2}.\end{cases}$ (4.7)

Periodic boundary conditions are imposed in both $x_{1}$\- and $x_{2}$\-
directions. With this setup, the exact solution is a standing wave

$\displaystyle u=\begin{cases}\sin(\sqrt{20}\pi t)\cos(4\pi x_{1})\cos(2\pi
x_{2}),\qquad\text{in}\quad\Omega_{1},\\\ \sin(\sqrt{20}\pi t)\cos(12\pi
x_{1})\cos(2\pi x_{2}),\qquad\text{in}\quad\Omega_{2}.\end{cases}$ (4.8)

For 3D case, $\Omega_{1}=[\frac{1}{4},\frac{3}{4}]\times[0,1]\times[0,1]$ and
$\Omega_{2}=\Omega\backslash\Omega_{1}$

$\displaystyle c^{2}=\begin{cases}1,\qquad\text{in}\quad\Omega_{1},\\\
{\frac{3}{19}},\qquad\text{in}\quad\Omega_{2}.\end{cases}$ (4.9)

Periodic boundary conditions are imposed in all directions. With this setup,
the exact solution is a standing wave

$\displaystyle u=\begin{cases}\sin(\sqrt{24}\pi t)\cos(4\pi x_{1})\cos(2\pi
x_{2})\cos(2\pi x_{3}),\qquad\text{in}\quad\Omega_{1},\\\ \sin(\sqrt{24}\pi
t)\cos(12\pi x_{1})\cos(2\pi x_{2})\cos(2\pi
x_{3}),\qquad\text{in}\quad\Omega_{2}.\end{cases}$ (4.10)

Since the solution is only piecewise smooth, the sparse grid DG method is not
expected to have good convergence rate. Therefore, we only show the
convergence result obtained by the adaptive method in Table 4.9 for both 2D
and 3D. In addition, the adaptive result with the parameter $N=8$ and
$\epsilon=1\times 10^{-4}$ in 2D is shown in Fig. 4.1. There are fewer DoFs in
the $x_{1}$ direction since the solution is smooth in that direction, and as
expected, there are more DoFs located in the subdomain $\Omega_{1}$ than that
in $\Omega\backslash\Omega_{1}$.

Table 4.9: Example 4.3. discontinuous coefficient, adaptive sparse grid DG, 2D and 3D. $t=0.01$. | $\epsilon$ | DoF |  $L^{2}$-error | $R_{\textrm{DoF}}$ | $R_{\epsilon}$
---|---|---|---|---|---
$d=2$ | 1e-1 | 480 | 2.93e-4 | - | -
1e-2 | 1088 | 8.43e-5 | 1.52 | 0.54
1e-3 | 2240 | 8.46e-6 | 3.18 | 1.00
1e-4 | 4224 | 1.05e-6 | 3.29 | 0.91
$d=3$ | 1e-1 | 2304 | 5.59e-4 | - | -
1e-2 | 7040 | 1.28e-4 | 1.32 | 0.64
1e-3 | 18176 | 1.65e-5 | 2.16 | 0.89
1e-4 | 41472 | 1.55e-6 | 2.87 | 1.03

(a) numerical solution

(b) centers of active elements

Figure 4.1: Example 4.3: Discontinuous coefficient in 2D at $t=0.1$. Adaptive
sparse grid DG with $N=8$ and $\epsilon=10^{-4}$.

######  Example 4.4 (Expanding wave in homogeneous medium).

We consider the wave equation with constant wave speed $c=1$ on the
computational domain $\Omega=[0,1]^{d}$. The homogeneous Neumann boundary
conditions are used on all boundaries. The initial condition is taken as

$u(\mathbf{x},0)=0,\quad u_{t}(\mathbf{x},0)=100e^{-500r^{2}}$ (4.11)

with $r=(\sum_{i=1}^{d}x_{i}^{2})^{1/2}$ being the radius.

For small time $t$ (before the wave front touch the outside boundary), the
exact solution in 2D can be represented by an integral which is derived by
Hadamard’s method of descent:

$u(x_{1},x_{2},t)=\frac{1}{2\pi}\iint_{\rho<t}\frac{100e^{-500(y_{1}^{2}+y_{2}^{2})}}{\sqrt{t^{2}-\rho^{2}}}dy_{1}dy_{2}$
(4.12)

with $\rho:=\sqrt{(y_{1}-x_{1})^{2}+(y_{2}-x_{2})^{2}}$ and then computed by
using numerical integrations with sufficiently small error tolerance. For
$d=3$, there exists the analytic solution:

$u(\mathbf{x},t)=\frac{1}{20r}\left(e^{-500(t-r)^{2}}-e^{-500(t+r)^{2}}\right).$
(4.13)

The numerical results for 2D and 3D are presented in Fig. 4.2 and Fig. 4.3. In
both cases, our numerical solutions coincide with the exact solutions quite
well. The $L^{\infty}$ errors between the numerical and the exact solutions at
$t=0.5$ are $9.06\times 10^{-5}$ and $6.79\times 10^{-4}$ for 2D and 3D, which
are both in the same magnitude as the adaptive parameter $\epsilon=1\times
10^{-4}$. This indicates that our adaptive algorithm controls the error really
well. The DoFs are 14896 and 188672 for 2D and 3D. It can be also observed
that the active elements in 3D are more “sparse” than 2D. This is a numerical
evidence that the Huyghens principle only holds for wave equations in odd
dimensions.

(a) numerical solution

(b) 1D cut along diagonal $x_{1}=x_{2}$

(c) error between exact and numerical solutions

(d) centers of active elements

Figure 4.2: Example 4.4: Expanding wave in homogeneous medium in 2D at
$t=0.5$. Adaptive sparse grid DG. $N=7$ and $\epsilon=10^{-4}$.

(a) numerical solution cut in 2D along $x_{3}=0$

(b) 1D cut along $x_{1}=x_{2}$ and $x_{3}=0$

(c) centers of active elements in 3D

(d) centers of active elements on $x_{3}=0.5$

Figure 4.3: Example 4.4: Expanding wave in homogeneous medium in 3D at
$t=0.5$. Adaptive sparse grid DG. $N=7$ and $\epsilon=10^{-4}$.

######  Example 4.5 (Isotropic wave propagation in heterogeneous media).

We consider the wave equation with discontinuous coefficient on the
computational domain $\Omega=[0,1]^{d}$ for $d=2$ and 3 [13].

$\displaystyle c^{2}=\begin{cases}\frac{1}{4},\qquad\text{if}\quad
0.35\leq{x}_{1}\leq 0.65,\\\ 1,\qquad\text{otherwise}.\end{cases}$ (4.14)

Note that the jump in material coefficient is not aligned with the cell
interface on $\Omega_{N}.$

For both 2D and 3D case, the initial conditions are taken as

$\displaystyle u(\mathbf{x},0)=0,\quad u_{t}(\mathbf{x},0)=100e^{-500r^{2}}.$
(4.15)

with $r=\left(\sum_{i=1}^{d}(x_{i}-\frac{1}{2})^{2}\right)^{\frac{1}{2}}$. The
zero Dirichlet boundary conditions are used.

The profiles and centers of active elements obtained by the adaptive scheme
are shown in Fig. 4.4 for 2D and Fig. 4.5 for 3D. We see that the wave fronts
propagate at different speeds in these two media and our adaptive scheme
capture this phenomenon and obtain comparable results to the literature [13].

(a) solution profile at $t=0.1$

(b) centers of active elements at $t=0.1$

(c) solution profile at $t=0.3$

(d) centers of active elements at $t=0.3$

Figure 4.4: Example 4.5: Isotropic wave propagation within heterogeneous media
in 2D at $t=0.1$ and $t=0.3$. Adaptive sparse grid DG. $N=8$ and
$\epsilon=10^{-4}$. Left: solution profile; right: centers of active elements.

(a) solution profile at $t=0.1$ (cut in 2D on $z=0.5$)

(b) centers of active elements at $t=0.1$

(c) solution profile at $t=0.3$ (cut in 2D on $z=0.5$)

(d) centers of active elements at $t=0.3$

Figure 4.5: Example 4.5: Isotropic wave propagation within heterogeneous media
in 3D at $t=0.1$ and $t=0.3$. Adaptive sparse grid DG. $N=7$ and
$\epsilon=10^{-4}$. Left: solution profile; right: centers of active elements.

##  5 Conclusion

In this paper, we develop an adaptive multiresolution DG scheme for wave
equations in second order form in multi-dimensions. Our method can achieve
similar computational complexity as the sparse grid DG method for smooth
solutions like those proposed for equations in first order form [21, 22].
Extensive numerical tests in 2D and 3D verify the accuracy and robustness of
the adaptive schemes for smooth and piecewise smooth wave propagation speed.
Though the formulation is based on IPDG scheme for scalar wave equation, it
can be extended to other DG method and other wave applications. Future work
includes extensions to other boundary conditions and investigation on
stability of schemes with interpolation. In an effort for promoting
reproducible research, the code generating the results in this paper can be
found at the github link: https://github.com/JuntaoHuang/adaptive-
multiresolution-DG.

## Acknowledgment

We would like to thank Daniel Appelö for discussions on numerical examples of
wave propagation, Qi Tang and Kai Huang for the assistance and discussion in
code implementation.

## References

  * [1] M. Ainsworth, P. Monk, and W. Muniz.  Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second-order wave equation.  Journal of Scientific Computing, 27(1-3):5–40, 2006.
  * [2] B. Alpert.  A class of bases in $L$2 for the sparse representation of integral operators.  SIAM Journal on Mathematical Analysis, 24(1):246–262, 1993.
  * [3] D. Appelö and T. Hagstrom.  A new discontinuous Galerkin formulation for wave equations in second-order form.  SIAM Journal on Numerical Analysis, 53(6):2705–2726, 2015.
  * [4] D. N. Arnold.  An interior penalty finite element method with discontinuous elements.  SIAM Journal on Numerical Analysis, 19(4):742–760, 1982.
  * [5] D. N. Arnold, F. Brezzi, B. Cockburn, and L. D. Marini.  Unified analysis of discontinuous Galerkin methods for elliptic problems.  SIAM Journal on Numerical Analysis, 39(5):1749–1779, 2002.
  * [6] M. J. Berger and P. Colella.  Local adaptive mesh refinement for shock hydrodynamics.  Journal of Computational Physics, 82(1):64–84, 1989.
  * [7] M. J. Berger and J. Oliger.  Adaptive mesh refinement for hyperbolic partial differential equations.  Journal of Computational Physics, 53(3):484 – 512, 1984.
  * [8] O. Bokanowski, J. Garcke, M. Griebel, and I. Klompmaker.  An adaptive sparse grid semi-Lagrangian scheme for first order Hamilton-Jacobi Bellman equations.  Journal of Scientific Computing, 55(3):575–605, 2013.
  * [9] D. L. Brown, G. S. Chesshire, W. D. Henshaw, and D. J. Quinlan.  Overture: An object-oriented software system for solving partial differential equations in serial and parallel environments.  Technical report, Los Alamos National Lab., NM (United States), 1997.
  * [10] H.-J. Bungartz and M. Griebel.  Sparse Grids.  Acta Numerica, 13:147–269, 2004.
  * [11] C. Burstedde, L. C. Wilcox, and O. Ghattas.  p4est: Scalable algorithms for parallel adaptive mesh refinement on forests of octrees.  SIAM Journal on Scientific Computing, 33(3):1103–1133, 2011.
  * [12] J. L. D. Calle, P. R. B. Devloo, and S. M. Gomes.  Wavelets and adaptive grids for the discontinuous Galerkin method.  Numerical Algorithms, 39(1-3):143–154, July 2005.
  * [13] C.-S. Chou, C.-W. Shu, and Y. Xing.  Optimal energy conserving local discontinuous Galerkin methods for second-order wave equation in heterogeneous media.  Journal of Computational Physics, 272:88–107, 2014.
  * [14] B. Cockburn, S. Hou, and C.-W. Shu.  The Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. IV. The multidimensional case.  Mathematics of Computation, 54(190):545–581, 1990.
  * [15] B. Cockburn, G. Karniadakis, and C.-W. Shu.  The development of discontinuous Galerkin methods.  In B. Cockburn, G. Karniadakis, and C.-W. Shu, editors, Discontinuous Galerkin methods: theory, computation and applications, volume 11, pages 3–50. Springer, 2000.
  * [16] G. C. Cohen.  Higher-Order Numerical Methods for Transient Wave Equations.  Springer, 2002.
  * [17] V. Etienne, E. Chaljub, J. Virieux, and N. Glinsky.  An hp-adaptive discontinuous Galerkin finite-element method for 3-d elastic wave modelling.  Geophysical Journal International, 183(2):941–962, 2010.
  * [18] D. Gottlieb and S. A. Orszag.  Numerical analysis of spectral methods: theory and applications.  SIAM, 1977.
  * [19] S. Gottlieb, C.-W. Shu, and E. Tadmor.  Strong stability preserving high order time discretization methods.  SIAM review, 43:89–112, 2001.
  * [20] M. J. Grote, A. Schneebeli, and D. Schötzau.  Discontinuous Galerkin finite element method for the wave equation.  SIAM Journal on Numerical Analysis, 44(6):2408–2431, 2006.
  * [21] W. Guo and Y. Cheng.  A sparse grid discontinuous Galerkin method for high-dimensional transport equations and its application to kinetic simulations.  SIAM Journal on Scientific Computing, 38(6):A3381–A3409, 2016.
  * [22] W. Guo and Y. Cheng.  An adaptive multiresolution discontinuous Galerkin method for time-dependent transport equations in multidimensions.  SIAM Journal on Scientific Computing, 39(6):A2962–A2992, 2017.
  * [23] B. Gustafsson, H.-O. Kreiss, and J. Oliger.  Time Dependent Problems and Difference Methods, volume 24.  John Wiley & Sons, 1995.
  * [24] W. D. Henshaw.  A high-order accurate parallel solver for Maxwell’s equations on overlapping grids.  SIAM Journal on Scientific Computing, 28(5):1730–1765, 2006.
  * [25] J. Hesthaven and T. Warburton.  Nodal high-order methods on unstructured grids: I. time-domain solution of Maxwell’s equations.  Journal of Computational Physics, 181:186–221, 2002.
  * [26] N. Hovhannisyan, S. Müller, and R. Schäfer.  Adaptive multiresolution discontinuous Galerkin schemes for conservation laws.  Mathematics of Computation, 83(285):113–151, 2014.
  * [27] J. Huang and Y. Cheng.  An adaptive multiresolution discontinuous Galerkin method with artificial viscosity for scalar hyperbolic conservation laws in multidimensions.  arXiv preprint arXiv:1906.00829, 2019.
  * [28] J. Huang and C.-W. Shu.  Error estimates to smooth solutions of semi-discrete discontinuous Galerkin methods with quadrature rules for scalar conservation laws.  Numerical Methods for Partial Differential Equations, 33(2):467–488, 2017.
  * [29] P. Joly.  Variational methods for time-dependent wave propagation problems., volume 31 of Topics in computational wave propagation, pages 201–264.  Springer, Berlin, 2003.
  * [30] M. Käser and M. Dumbser.  An arbitrary high-order discontinuous Galerkin method for elastic waves on unstructured meshes—i. the two-dimensional isotropic case with external source terms.  Geophysical Journal International, 166(2):855–877, 2006.
  * [31] S. Mallat.  A wavelet tour of signal processing.  Elsevier, 1999.
  * [32] W. Reed and T. Hill.  Tiangular mesh methods for the neutron transport equation.  Technical report, Los Alamos National Laboratory, Los Alamos, NM, 1973\.
  * [33] G. Seriani and E. Priolo.  Spectral element method for acoustic wave simulation in heterogeneous media.  Finite Elements in Analysis and Design, 16(3-4):337–348, 1994.
  * [34] J. Shen and H. Yu.  Efficient spectral sparse grid methods and applications to high-dimensional elliptic problems.  SIAM Journal on Scientific Computing, 32(6):3228–3250, 2010.
  * [35] C.-W. Shu and S. Osher.  Efficient implementation of essentially non-oscillatory shock-capturing schemes.  Journal of Computational Physics, 77(2):439–471, 1988.
  * [36] B. Sjögreen and N. A. Petersson.  A fourth order accurate finite difference scheme for the elastic wave equation in second order formulation.  Journal of Scientific Computing, 52(1):17–48, 2012.
  * [37] Z. Tao, Y. Jiang, and Y. Cheng.  An adaptive high-order piecewise polynomial based sparse grid collocation method with applications.  arXiv preprint arXiv:1912.03982, 2019.
  * [38] J. Tromp, D. Komatitsch, and Q. Liu.  Spectral-element and adjoint methods in seismology.  Communications in Computational Physics, 3(1):1–32, 2008.
  * [39] Z. Wang, Q. Tang, W. Guo, and Y. Cheng.  Sparse grid discontinuous Galerkin methods for high-dimensional elliptic equations.  Journal of Computational Physics, 314:244–263, 2016.
  * [40] L. C. Wilcox, G. Stadler, C. Burstedde, and O. Ghattas.  A high-order discontinuous Galerkin method for wave propagation through coupled elastic–acoustic media.  Journal of Computational Physics, 229(24):9373–9396, 2010.
  * [41] Y. Xing, C.-S. Chou, and C.-W. Shu.  Energy conserving local discontinuous Galerkin methods for wave propagation problems.  Inverse Problems & Imaging, 7(3), 2013.
  * [42] A. Zeiser.  Fast matrix-vector multiplication in the sparse-grid Galerkin method.  Journal of Scientific Computing, 47(3):328–346, 2011.

##  Appendix A Interpolation basis functions

For completeness of the paper, we present details of the multiresolution
interpolation basis functions, which are first introduced in [37]. We will
first focus on the case in which the interpolation points are imposed in the
inner domain, as implemented in Table 4.5. Then we discuss the case in which
the points includes the cell interface points. Here, we only discuss the case
when $M=4$ and $M=5$. For $M=1,2,3$, we refer readers to the appendix in [27].

The basis functions in $\tilde{W}_{1}$ are piecewise polynomials on
$I_{l}:=(0,\frac{1}{2})$ and $I_{r}:=(\frac{1}{2},1)$. Note that the functions
may be discontinuous at the interface $x=1/2$, thus $I_{l}$ and $I_{r}$ are
both defined to be open intervals. The basis functions in $\tilde{W}_{1}$ in
this paper are all supported on one half interval $I_{l}$ or $I_{r}$ and
vanish on the other half. For simplicity, we will only declare the function on
its support. For example, $\psi_{0}(x)|_{I_{r}}$ gives the definition of
$\psi_{0}$ on $I_{r}$ and indicates that $\psi_{0}$ vanishes on $I_{l}$.

###  A.1 Interpolation points in the inner domain

####  A.1.1 $M=4$

The interpolation points are

$\tilde{X}_{0}=\\{\frac{1}{6},\frac{7}{24},\frac{1}{3},\frac{7}{12},\frac{2}{3}\\},\quad\tilde{X}_{1}=\\{\frac{1}{12},\frac{7}{48},\frac{31}{48},\frac{19}{24},\frac{5}{6}\\}.$

The basis functions in $\tilde{W}_{0}^{4}$ and $\tilde{W}_{1}^{4}$ are

$\displaystyle\begin{array}[]{ll}\phi_{0}(x)=\frac{4}{45}(3x-2)(3x-1)(12x-7)(24x-7),\\\
\phi_{1}(x)=-\frac{512}{189}(3x-2)(3x-1)(6x-1)(12x-7),\\\
\phi_{2}(x)=\frac{1}{3}(3x-2)(6x-1)(12x-7)(24x-7),\\\
\phi_{3}(x)=-\frac{32}{105}(3x-2)(3x-1)(6x-1)(24x-7),\\\
\phi_{4}(x)=\frac{1}{27}(3x-1)(6x-1)(12x-7)(24x-7).\end{array}$

and

$\displaystyle\begin{array}[]{l}\psi_{0}(x)|_{I_{l}}=\frac{8}{45}(3x-1)(6x-1)(24x-7)(48x-7),\\\
\psi_{1}(x)|_{I_{l}}=-\frac{1024}{189}(3x-1)(6x-1)(12x-1)(24x-7),\\\
\psi_{2}(x)|_{I_{r}}=-\frac{1024}{189}(3x-2)(6x-5)(12x-7)(24x-19),\\\
\psi_{3}(x)|_{I_{r}}=-\frac{64}{105}(3x-2)(6x-5)(12x-7)(48x-31),\\\
\psi_{4}(x)|_{I_{r}}=\frac{2}{27}(3x-2)(12x-7)(24x-19)(48x-31)\end{array}$

####  A.1.2 $M=5$

The interpolation points are

$\tilde{X}_{0}=\\{\frac{1}{12},\frac{1}{6},\frac{7}{24},\frac{1}{3},\frac{7}{12},\frac{2}{3}\\},\quad\tilde{X}_{1}=\\{\frac{7}{48},\frac{1}{24},\frac{31}{48},\frac{19}{24},\frac{5}{6},\frac{13}{24}\\}.$

The basis functions in $\tilde{W}_{0}^{5}$ and $\tilde{W}_{1}^{5}$ are

$\displaystyle\begin{array}[]{ll}\phi_{0}(x)=\frac{1}{315}(-16)(3x-2)(3x-1)(6x-1)(12x-7)(24x-7),\\\
\phi_{1}(x)=\frac{4}{45}(3x-2)(3x-1)(12x-7)(12x-1)(24x-7),\\\
\phi_{2}(x)=-\frac{1024}{945}(3x-2)(3x-1)(6x-1)(12x-7)(12x-1),\\\
\phi_{3}(x)=\frac{1}{9}(3x-2)(6x-1)(12x-7)(12x-1)(24x-7),\\\
\phi_{4}(x)=-\frac{16}{315}(3x-2)(3x-1)(6x-1)(12x-1)(24x-7),\\\
\phi_{5}(x)=\frac{1}{189}(3x-1)(6x-1)(12x-7)(12x-1)(24x-7),\end{array}$

and

$\displaystyle\begin{array}[]{l}\psi_{0}(x)|_{I_{l}}=-\frac{2048}{945}(3x-1)(6x-1)(12x-1)(24x-7)(24x-1),\\\
\psi_{1}(x)|_{I_{l}}=-\frac{32}{315}(3x-1)(6x-1)(12x-1)(24x-7)(48x-7),\\\
\psi_{2}(x)|_{I_{r}}=-\frac{2048}{945}(3x-2)(6x-5)(12x-7)(24x-19)(24x-13),\\\
\psi_{3}(x)|_{I_{r}}=-\frac{32}{315}(3x-2)(6x-5)(12x-7)(24x-13)(48x-31),\\\
\psi_{4}(x)|_{I_{r}}=\frac{2}{189}(3x-2)(12x-7)(24x-19)(24x-13)(48x-31),\\\
\psi_{5}(x)|_{I_{r}}=-\frac{32}{315}(3x-2)(6x-5)(12x-7)(24x-19)(48x-31)\end{array}$

###  A.2 Interpolation points with the interface points

####  A.2.1 $M=4$

The interpolation points are

$\tilde{X}_{0}=\\{0^{+},\left(\frac{1}{4}\right)^{-},\left(\frac{1}{2}\right)^{-},\left(\frac{3}{4}\right)^{-},1^{-}\\},\quad\tilde{X}_{1}=\\{\left(\frac{1}{8}\right)^{-},\left(\frac{3}{8}\right)^{-},\left(\frac{1}{2}\right)^{+},\left(\frac{5}{8}\right)^{-},\left(\frac{7}{8}\right)^{-}\\}.$

The basis functions in $\tilde{W}_{0}^{4}$ and $\tilde{W}_{1}^{4}$ are

$\displaystyle\begin{array}[]{ll}\phi_{0}(x)=\frac{1}{3}(x-1)(2x-1)(4x-3)(4x-1),\\\
\phi_{1}(x)=-\frac{16}{3}(x-1)x(2x-1)(4x-3),\\\
\phi_{2}(x)=4(x-1)x(4x-3)(4x-1),\\\
\phi_{3}(x)=-\frac{16}{3}(x-1)x(2x-1)(4x-1),\\\
\phi_{4}(x)=\frac{1}{3}x(2x-1)(4x-3)(4x-1).\end{array}$

and

$\displaystyle\begin{array}[]{l}\psi_{0}(x)|_{I_{l}}=-\frac{32}{3}x(2x-1)(4x-1)(8x-3),\\\
\psi_{1}(x)|_{I_{l}}=-\frac{32}{3}x(2x-1)(4x-1)(8x-1),\\\
\psi_{2}(x)|_{I_{r}}=\frac{2}{3}(x-1)(4x-3)(8x-7)(8x-5),\\\
\psi_{3}(x)|_{I_{r}}=-\frac{32}{3}(x-1)(2x-1)(4x-3)(8x-7),\\\
\psi_{4}(x)|_{I_{r}}=-\frac{32}{3}(-32)(x-1)(2x-1)(4x-3)(8x-5)\end{array}$

####  A.2.2 $M=5$

The interpolation points are

$\tilde{X}_{0}=\\{0^{+},\frac{1}{5},\frac{2}{5},\frac{3}{5},\frac{4}{5},1^{-}\\},\quad\tilde{X}_{1}=\\{\frac{1}{10},\frac{3}{10},\left(\frac{1}{2}\right)^{-},\left(\frac{1}{2}\right)^{+},\frac{7}{10},\frac{9}{10}\\}.$

The basis functions in $\tilde{W}_{0}^{5}$ and $\tilde{W}_{1}^{5}$ are

$\displaystyle\begin{array}[]{ll}\phi_{0}(x)=-\frac{1}{24}(x-1)(5x-4)(5x-3)(5x-2)(5x-1),\\\
\phi_{1}(x)=\frac{25}{24}(x-1)x(5x-4)(5x-3)(5x-2),\\\
\phi_{2}(x)=-\frac{25}{12}(x-1)x(5x-4)(5x-3)(5x-1),\\\
\phi_{3}(x)=\frac{25}{12}(x-1)x(5x-4)(5x-2)(5x-1),\\\
\phi_{4}(x)=-\frac{25}{24}(x-1)x(5x-3)(5x-2)(5x-1),\\\
\phi_{5}(x)=\frac{1}{24}x(5x-4)(5x-3)(5x-2)(5x-1),\end{array}$

and

$\displaystyle\begin{array}[]{l}\psi_{0}(x)|_{I_{l}}=\frac{25}{3}x(2x-1)(5x-2)(5x-1)(10x-3),\\\
\psi_{1}(x)|_{I_{l}}=\frac{50}{3}x(2x-1)(5x-2)(5x-1)(10x-1),\\\
\psi_{2}(x)|_{I_{r}}=\frac{1}{3}x(5x-2)(5x-1)(10x-3)(10x-1),\\\
\psi_{3}(x)|_{I_{r}}=-\frac{1}{3}(x-1)(5x-4)(5x-3)(10x-9)(10x-7),\\\
\psi_{4}(x)|_{I_{r}}=-\frac{50}{3}(x-1)(2x-1)(5x-4)(5x-3)(10x-9),\\\
\psi_{5}(x)|_{I_{r}}=-\frac{25}{3}(x-1)(2x-1)(5x-4)(5x-3)(10x-7).\end{array}$