#  Global space-time Trefftz DG schemes for the time-dependent linear
nonhomogeneous and anisotropic wave equation ††thanks: The author was
supported by Shandong Provincial Natural Science Foundation under the grant
ZR2020MA046.

LONG YUAN  College of Mathematics and Systems Science, Shandong University of
Science and Technology, Qingdao 266590, China (yuanlong@lsec.cc.ac.cn).

###### Abstract

In this paper we are concerned with Trefftz discretizations of the time-
dependent linear wave equation in anisotropic media in arbitrary space
dimensional domains $\Omega\subset\mathbb{R}^{d}~{}(d\in\mathbb{N})$. We
propose two variants of the Trefftz DG method, define novel plane wave basis
functions based on rigorous choices of scaling transformations and coordinate
transformations, and prove that the corresponding approximate solutions
possess optimal-order error estimates with respect to the meshwidth $h$ and
the condition number of the coefficient matrices, respectively. Furthermore,
we propose the global Trefftz DG method combined with local DG methods to
solve the time-dependent linear nonhomogeneous wave equation in anisotropic
media. In particular, the error analysis holds for the (nonhomogeneous)
Dirichlet, Neumann, and mixed boundary conditions from the original PDEs. The
numerical results verify the validity of the theoretical results, and show
that the resulting approximate solutions possess high accuracy.

###### keywords:

time-dependent wave equation, nonhomogeneous, anisotropic, local discontinuous
Galerkin, Trefftz method, error estimates.

###### AMS:

65N30, 65N55.

##  1 Introduction

The idea at the heart of Trefftz method, which are named after the seminal
work [19] of E. Trefftz, is to choose the Trefftz approximation functions from
a class of piecewise solutions of the same governing partial differential
equation (PDE) without boundary conditions. Trefftz methods turned out to be
particularly effective, and popular, for wave propagation problems in time-
harmonic regime at medium and high frequencies, where the oscillatory nature
of the solutions makes standard methods computationally too expensive; see the
recent survey [7] and references therein. The Trefftz method has an important
advantage over Lagrange finite elements for discretization of the Helmholtz
equation and time-harmonic Maxwell equations [12, 13, 16, 14, 20, 21, 36]: to
achieve the same accuracy, relatively smaller degrees of freedom are enough in
the plane wave-type methods owing to the particular choice of the basis
functions that (may approximately) satisfy the considered PDE without boundary
conditions.

Much work has been devoted to Trefftz discontinuous Galerkin (DG) methods for
time-dependent linear isotropic wave phenomena, see in particular [1, 6, 7,
22, 24, 25, 27, 31]. A Trefftz interior-penalty formulation for the second
order wave equation is proposed and studied in [2]. Another Trefftz DG method
for time-dependent electromagnetic problems has been analysed in [24, 25] in
one space dimension, and then it has been extended to three-dimensional time
dependent Maxwell's equations in [6, 7, 22]. The stabilized space-time
discretizations for linear wave equations were proposed by [3, 18, 19, 34].
The work [18, 19] investigates the error bounds of DG approximations for
smooth solutions and numerical examples in one spatial dimension were
reported. The work [3] gives numerical analysis results of space-time
discretization schemes for evolution equations by the assumption on the
maximal regularity of the analytic solution in the spatial and the temporal
variable. Moreover, the recent work on explicit marching-type space-time
schemes, which are so-called tent-pitched space-time meshes, are fully
investigated by [8, 10, 11, 31], where the PDE are explicitly evolved from the
``bottom" to the ``top" of the space-time cylinder element by element.

Recently, a class of space-time DG discretizations of the linear isotropic
acoustic wave equation in two space dimensions in polygonal domains occupied
by possibly heterogeneous media are addressed in [1]. The consistency analysis
of the space-time discretization are generalized to non-Trefftz discrete
spaces, and the realistic setting of solutions exhibiting spatial point
singularities is allowed, where high $h-$convergence rates of the
approximations generated by the space-time scheme with local corner mesh
refinement on the spatial domain still hold. Moreover, the consistency error
bounds in mesh-dependent norm holds true without any time-step size
constraint.

Anisotropy can result from periodic layering of fine layers [4], preferential
alignment of fractures and cracks [33]. Anisotropy may greatly influence
seismic wave propagation, seismic data acquisition and subsequent data
analysis and processing procedures [35]. It is therefore important to design
accurate and efficient numerical methods for modeling wave propagation in
anisotropic media. An acoustic wave equation for anisotropic media in [23] is
introduced to describe a wave type that propagates at speeds slower than the
P-wave for a positive anisotropy parameter. An improved rotated staggered-grid
finite-difference method in [9] with fourth-order temporal accuracy has been
developed to solve elastic-wave modeling in anisotropic media, where the
symmetry axes of anisotropy are not aligned with the coordinate axes. A
compensated-amplitude vertical transverse isotropic least-squares reverse time
migration method in [32] is adopted to correct the anisotropy effect and
compensate amplitude attenuation. [37] develops a closed expression of Riemann
solvers for the discontinuous Galerkin time domain method, applied to wave
propagation modeling in distinct anisotropic material properties. Recently,
the PWDG methods [36] have been developed to solve Helmholtz equation and
time-harmonic Maxwell equations in three-dimensional anisotropic media.

In this paper we construct a global space-time Trefftz DG scheme for the
linear anisotropic wave equation in arbitrary space dimensional domains
$\Omega\subset\mathbb{R}^{d}~{}(d\in\mathbb{N})$. In order to build better
convergence results, we have to carefully define plane wave basis functions by
rigorous choices of the scaling transformations and the coordinate
transformations. We prove that the approximate solutions generated by the
proposed method possess satisfactory and optimal-order error estimates with
respect to meshwidth $h$ and the condition number $\rho$ of the coefficient
matrices, respectively. Besides, we propose an alternative standard Trefftz DG
method with almost the same computational cost and with the same convergence
order with respect to $h$ and $\rho$ in Section 6.

Numerical results indicates that, the consistency error bounds in the mesh-
dependent $|\cdot|_{\text{DG}}-$norm and mesh-independent
$L^{2}(\Omega\times\\{T\\})-$norm for Trefftz DG methods are optimal with
respect to $h$ and $\rho$, respectively. Moreover, numerical experiments in
Section 8.3 show that the approximations generated by the Trefftz DG methods
are clearly more accurate than that generated by the high-order DG finite
element method.

Since Trefftz basis functions on each element are solutions of the homogeneous
wave equations without boundary conditions, the Trefftz methods can not be
directly applied to discretizations of the nonhomogeneous wave equations.
Motivated by the coupled discontinuous Galerkin formulation developed in [17],
we develop the global Trefftz DG method combined with overlapping local DG
method. Numerical results indicates that, the consistency error bounds in the
mesh-dependent $|\cdot|_{\text{DG}}-$norm and mesh-independent
$L^{2}(\Omega\times\\{T\\})-$norm for Trefftz DG method combined with
overlapping local DG are optimal with respect to $h$ and $\rho$, respectively,
in the presence of the nonhomogeneous source and the anisotropic media.
Besides, we propose another alternative to define nonhomogeneous local
problems in each nonoverlapping time slab, and solve it by the space-time DG
method. The resulting residue problem on the global solution domain is still
solved by the Trefftz DG method. We call the new method as `` Trefftz DG
method combined with nonoverlapping local DG". Numerical results in Section
8.4 indicates that, Trefftz DG method combined with overlapping local DG is
comparable to Trefftz DG method combined with overlapping local DG.

Comparing against the space-time DG method of [1] for isotropic wave equations
which employs the piecewise-polynomial discrete space, our global (resp.
local) discrete space is available for the Trefftz (resp. DG finite element)
space in arbitrary space dimensions, and the Dirichlet, Neumann, and mixed
boundary conditions on space-time domain boundary are considered. In
particular, there is no constraint on the data of homogeneous Neumann boundary
conditions from the original initial boundary value problem (IBVP), thus the
analysis holds for the nonhomogeneous Neumann boundary conditions from the
original IBVP, which has been also verified by the numerical tests in Section
8.3.

Moreover, in comparison to most existing Trefftz methods [22, 24, 27] for the
isotropic wave equations, the proposed variational formulation with three
relaxation parameters is applied to the anisotropic model with Dirichlet,
Neumann, and mixed boundary conditions; combined with the local DG method on
auxiliary smooth subdomains, the new method can generate the approximations
with the spectral convergence orders for the nonhomogeneous case; the space-
time domain partition ${\cal T}_{h}$ is obtained as the tensor product of
space and time mesh grids ${\cal T}_{h_{{\bf x}}}^{\bf x}$ and ${\cal
T}_{h_{t}}^{t}$, the space-time mesh $\hat{\cal T}_{\hat{h}}$ of the
transformed space-time domain $\hat{Q}$ satisfies quasi-uniform assumption,
and the space mesh grid $\hat{\cal T}_{\hat{h}_{\hat{\bf x}}}^{\hat{\bf x}}$
of $\hat{\Omega}$ is shape regular and quasi-uniform.

The paper is organized as follows: In Section 2, we state the initial boundary
value problem for the acoustic wave equation in both first- and second-order
formulation. Section 3 describes the proposed method for nonhomogeneous PDEs.
In Section 4, we explain how to discretize the resulting variational problems.
Section 5 provides the desired error estimates for the approximate solutions.
In Section 6, we propose an alternative standard Trefftz DG method. In order
to solve the nonhomogeneous and anisotropic model, we develop a global Trefftz
DG method combined with overlapping local DG method for the nonhomogeneous in
Section 7. Finally, we report some numerical results to confirm the
effectiveness of the proposed method .

##  2 Considered model

We consider the first order acoustic wave IBVP posed on a space-time domain
$Q=\Omega\times I$, where $\Omega\subset\mathbb{R}^{d}~{}(d\in\mathbb{N})$ is
an open bounded Lipschitz polytope and $I=(0,T),T>0$. ${\bf n}^{x}_{\Omega}$
is an outward-pointing unit normal vector on $\partial\Omega$. The boundary of
$\Omega$ denoted by $\Gamma$, is divided in two parts, with mutually disjoint
interiors, denoted $\Gamma_{D}$ or $\Gamma_{N}$ corresponding to Dirichlet and
Neumann boundary conditions, respectively. The model reads as

(1) $\left\\{\begin{aligned} &A^{\frac{1}{2}}\nabla
v+\frac{\partial{\bm{\sigma}}}{\partial t}={\bf 0}&\text{in}\quad Q,\\\
&\nabla\cdot A^{\frac{1}{2}}{\bm{\sigma}}+c^{-2}\frac{\partial v}{\partial
t}=0&\text{in}\quad Q,\\\
&v(\cdot,0)=v_{0},\quad{\bm{\sigma}}(\cdot,0)={\bm{\sigma}}_{0}&\text{on}\quad\Omega,\\\
&v=g_{D}&\text{on}\quad\Gamma_{D}\times[0,T],\\\
&A^{\frac{1}{2}}{\bm{\sigma}}\cdot{\bf
n}^{x}_{\Omega}=g_{N}&\text{on}\quad\Gamma_{N}\times[0,T].\end{aligned}\right.$

Here $v_{0},{\bm{\sigma}}_{0},g_{D},g_{N}$ are the given source data, $c$ is
the wave speed, which is constant in the whole space domain $\Omega$ and
independent of time $t$. $A$ is positive definite matrix independent of ${\bf
x}$ and $t$. The gradient $\nabla$ and divergence $\nabla\cdot$ operators are
meant in the space variable ${\bf x}$ only.

If there exists a scalar field $U_{0}$ such that
${\bm{\sigma}}_{0}=-A^{\frac{1}{2}}\nabla U_{0}$, then IBVP (1) is equivalent
to the following second order scalar wave equation, by setting
$v=\frac{\partial U}{\partial t}$ and ${\bm{\sigma}}=-A^{\frac{1}{2}}\nabla
U$,

(2) $\left\\{\begin{aligned} &-\nabla\cdot(A\nabla
U)+c^{-2}\frac{\partial^{2}U}{\partial t^{2}}=0&\text{in}\quad Q,\\\
&\frac{\partial U}{\partial t}(\cdot,0)=v_{0},\quad
U(\cdot,0)=U_{0}&\text{on}\quad\Omega,\\\ &\frac{\partial U}{\partial
t}=g_{D}&\text{on}\quad\Gamma_{D}\times[0,T],\\\ &-A\nabla U\cdot{\bf
n}^{x}_{\Omega}=g_{N}&\text{on}\quad\Gamma_{N}\times[0,T].\end{aligned}\right.$

Let the time domain $(0,T)$ be divided into $N\in\mathbb{N}$ intervals
$I_{n}(1\leq n\leq N)$ composing a partition ${\cal T}_{h_{t}}^{t}$, with

$I_{n}=(t_{n-1},t_{n}),~{}~{}h_{n}=t_{n}-t_{n-1}=|I_{n}|,~{}~{}h_{t}=\mathop{\text{missing}}{max}\limits_{1\leq
n\leq N}h_{n}.$

Let us introduce the following notation for the time slabs and the partial
cylinders, respectively,

$D_{n}=\Omega\times I_{n},~{}~{}Q_{n}=\Omega\times(0,t_{n}),~{}~{}1\leq n\leq
N.$

For each $1\leq n\leq N$, we introduce a same polygonal finite element mesh
${\cal T}_{h_{{\bf x}}}^{\bf x}=\\{K_{\bf x}\\}$ of the spatial domain
$\Omega$ with

$h_{K_{\bf x}}=\text{diam}K_{\bf x},~{}~{}h_{{\bf
x}}=\mathop{\text{missing}}{max}\limits_{K_{\bf x}\in{\cal T}_{h_{\bf
x}}^{x}}h_{K_{\bf x}}.$

Then the space-time domain $Q=\Omega\times(0,T)$ can be partitioned with a
finite element mesh ${\cal T}_{h}$ given by

${\cal T}_{h}=\\{K=K_{\bf x}\times I_{n},K_{\bf x}\in{\cal T}_{h_{{\bf
x}}}^{\bf x},1\leq n\leq N\\}.$

Here ${\cal T}_{h}$ is a tensor product mesh. Besides, we define the time-
truncated mesh

${\cal T}_{h}(Q_{n})=\\{K\in{\cal T}_{h},K\subset Q_{n}\\},~{}1\leq n\leq N$

and

${\cal T}_{h}(D_{n})=\\{K\in{\cal T}_{h},K\subset D_{n}\\},~{}1\leq n\leq N.$

Assume the space-time grid ${\cal T}_{h}$ satisfies the assumptions presented
in Section 4 of [1]: on an internal face $F=\partial K_{1}\bigcap\partial
K_{2}$, either

(3) $\left\\{\begin{aligned} &{\bf n}_{F}^{\bf x}=0&\text{and}~{}F~{}\text{is
called ``space-like" face, or}\\\ &{\bf n}_{F}^{t}=0&\text{and}~{}F~{}\text{is
called ``time-like" face},\\\ \end{aligned}\right.$

where $({\bf n}_{F}^{\bf x},{\bf n}_{F}^{t})$ is a unit vector of the face
$F$. On space-like faces, by convention, we choose ${\bf n}_{F}^{t}>0$, which
means that the unit normal vector $({\bf n}_{F}^{\bf x},{\bf n}_{F}^{t})$
points towards future time. Moreover, all time-like faces are of the form
$F=F_{\bf x}\times F_{t}$ with $h_{F_{\bf x}}=|F_{\bf x}|$ and
$h_{F_{t}}=|F_{t}|$; we recall that $F_{t}=I_{n},1\leq n\leq N$. Finally, We
denote the outward-pointing unit normal vector on $\partial K$ by $({\bf
n}_{K}^{\bf x},{n}_{K}^{t})$ .

We denote by $\mathcal{F}_{h}=\bigcup\limits_{K\in{\cal T}_{h}}\partial K$ the
skeleton of the mesh, by $\mathcal{F}_{h}^{\text{space}}$ the union of the
internal space-like faces, and by $\mathcal{F}_{h}^{\text{time}}$ the union of
the internal time-like faces, respectively. Set
$\mathcal{F}_{h}^{0}=\Omega\times\\{t=0\\}$,
$\mathcal{F}_{h}^{T}=\Omega\times\\{t=T\\}$,
$\mathcal{F}_{h}^{D}=\Gamma_{D}\times[0,T]$ and
$\mathcal{F}_{h}^{N}=\Gamma_{N}\times[0,T]$.

Let $w$ and ${\bm{\tau}}$ be a piecewise smooth function and vector field on
${\cal T}_{h}$, respectively. On $F=\partial K_{1}\bigcap\partial K_{2}$, we
define

(4) $\displaystyle\text{the
averages:}~{}~{}\\{\\{w\\}\\}:=\frac{w_{|K_{1}}+w_{|K_{2}}}{2},~{}~{}\\{\\{{\bm{\tau}}\\}\\}:=\frac{{\bm{\tau}}_{|K_{1}}+{\bm{\tau}}_{|K_{2}}}{2},$
(5) $\displaystyle\text{space normal jumps:}~{}\llbracket{w}\rrbracket_{\bf
N}:=w_{|K_{1}}{\bf n}^{\bf x}_{K_{1}}+w_{|K_{2}}{\bf n}^{\bf
x}_{K_{2}},~{}\llbracket{\bm{\tau}}\rrbracket_{\bf
N}={\bm{\tau}}_{|K_{1}}\cdot{\bf n}^{\bf
x}_{K_{1}}+{\bm{\tau}}_{|K_{2}}\cdot{\bf n}^{\bf x}_{K_{2}},$ (6)
$\displaystyle\text{time full jumps:}~{}\llbracket
w\rrbracket_{t}:=w_{|K_{1}}{n}^{t}_{K_{1}}+w_{|K_{2}}{n}^{t}_{K_{2}}=(w^{-}-w^{+}){n}^{t}_{F},$
(7) $\displaystyle\text{time full
jumps:}~{}\llbracket{\bm{\tau}}\rrbracket_{t}:={\bm{\tau}}_{|K_{1}}{n}^{t}_{K_{1}}+{\bm{\tau}}_{|K_{2}}{n}^{t}_{K_{2}}=({\bm{\tau}}^{-}-{\bm{\tau}}^{+}){n}^{t}_{F}.$

Here $w^{-}$ and $w^{+}$ denote the traces of the function $w$ from the
adjacent elements at lower and higher times, respectively, and similarly for
${\bm{\tau}}^{\pm}$.

##  3 The variational formulation

Set local Trefftz space:

(9) $\displaystyle{\bf T}(K)$ $\displaystyle=$
$\displaystyle\bigg{\\{}(w,{\bm{\tau}})\in
H^{1}(K)^{1+d}s.t.~{}{\bm{\tau}}|_{\partial K}\in L^{2}(\partial
K)^{d},~{}\frac{\partial w}{\partial t},\nabla\cdot{\bm{\tau}}\in L^{2}(K),$
$\displaystyle\frac{\partial\bm{\tau}}{\partial t},\nabla w\in
L^{2}(K)^{d},~{}A^{\frac{1}{2}}\nabla w+\frac{\partial{\bm{\tau}}}{\partial
t}={\bf 0},\nabla\cdot A^{\frac{1}{2}}{\bm{\tau}}+c^{-2}\frac{\partial
w}{\partial t}=0\bigg{\\}}~{}\quad\forall K\in{\cal T}_{h},$

and set global Trefftz space:

(10) $\displaystyle{\bf T}({\cal T}_{h})=\bigg{\\{}(w,{\bm{\tau}})\in
L^{2}(Q)^{1+d}~{}s.t.~{}(w|_{K},{\bm{\tau}}|_{K})\in{\bf T}(K)~{}\forall
K\in{\cal T}_{h}\bigg{\\}}.$

To derive the Trefftz-DG variational formulation, we multiply the first two
equation of (1) with test fields ${\bm{\tau}}$ and $w$ and integrated by parts
on each $K\in{\cal T}_{h}$:

(11) $\displaystyle-\int_{K}$ $\displaystyle\bigg{(}v(\nabla\cdot
A^{\frac{1}{2}}{\bm{\tau}}+c^{-2}\frac{\partial w}{\partial
t})+{\bm{\sigma}}\cdot(A^{\frac{1}{2}}\nabla
w+\frac{\partial\bm{\tau}}{\partial t})\bigg{)}dV$ (12) $\displaystyle+$
$\displaystyle\int_{\partial
K}\bigg{(}(vA^{\frac{1}{2}}{\bm{\tau}}+wA^{\frac{1}{2}}{\bm{\sigma}})\cdot{\bf
n}^{x}_{K}+({\bm{\sigma}}\cdot{\bm{\tau}}+c^{-2}vw)n^{t}_{K}\bigg{)}dS=\int_{K}f~{}w~{}dV,~{}\forall(w,{\bm{\tau}})\in{\bf
T}({\cal T}_{h}).$

Replacing the traces of $v$ and ${\bm{\sigma}}$ on the mesh skeleton by the
single-valued numerical fluxes $\check{v}$ and $\check{\bm{\sigma}}$, we have

(13) $\int_{\partial K}\bigg{(}\check{v}(A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf
n}^{x}_{K}+c^{-2}wn^{t}_{K})+\check{\bm{\sigma}}\cdot(wA^{\frac{1}{2}}{\bf
n}^{x}_{K}+{\bm{\tau}}n^{t}_{K})\bigg{)}dS=\int_{K}f~{}w~{}dV.$

Define the numerical fluxes as follows.

$\check{v}=\left\\{\begin{aligned} &v^{-}\\\ &v\\\ &v_{0}\\\
&\\{\\{v\\}\\}+\beta\llbracket A^{\frac{1}{2}}{\bm{\sigma}}\rrbracket_{N}\\\
&g_{D}\\\ &v+\beta(A^{\frac{1}{2}}{\bm{\sigma}}\cdot{\bf
n}^{x}_{\Omega}-g_{N})\end{aligned}\right.\check{\bm{\sigma}}=\left\\{\begin{aligned}
&{\bm{\sigma}}^{-}&\quad{\text{o}n}~{}{\mathcal{F}}_{h}^{\text{space}},\\\
&{\bm{\sigma}}&\quad{\text{o}n}~{}{\mathcal{F}}_{h}^{T},\\\
&{\bm{\sigma}}_{0}&\quad{\text{o}n}~{}{\mathcal{F}}_{h}^{0},\\\
&\\{\\{\bm{\sigma}\\}\\}+\alpha A^{\delta}\llbracket
v\rrbracket_{N}&\quad{\text{o}n}~{}{\mathcal{F}}_{h}^{\text{time}},\\\
&\bm{\sigma}+\alpha(v-g_{D})A^{\delta}{\bf
n}^{x}_{\Omega}&\quad{\text{o}n}~{}{\mathcal{F}}_{h}^{D},\\\
&g_{N}A^{-\frac{1}{2}}{\bf
n}^{x}_{\Omega}&\quad{\text{o}n}~{}{\mathcal{F}}_{h}^{N}.\end{aligned}\right.$

The stabilization parameters $\alpha\in
L^{\infty}({\mathcal{F}}_{h}^{\text{time}}\bigcup{\mathcal{F}}_{h}^{D})$,
$\beta\in
L^{\infty}({\mathcal{F}}_{h}^{\text{time}}\bigcup{\mathcal{F}}_{h}^{N})$ are
positive constant on each time-like face. $\delta\in
L^{\infty}({\mathcal{F}}_{h}^{\text{time}}\bigcup{\mathcal{F}}_{h}^{D})$ is
constant on each time-like face, and its best choice will be given in Lemma 3.

By summing the elemental DG equation over the element $K\in{\bf T}({\cal
T}_{h})$ and using the defined fluxes, we can obtain the Trefftz-DG
variational formulation: Find $(v,{\bm{\sigma}})\in{\bf T}({\cal T}_{h})$ such
that

(14)
$\mathcal{A}(v,{\bm{\sigma}};w,{\bm{\tau}})=\ell(w,{\bm{\tau}})\quad\forall(w,{\bm{\tau}})\in{\bf
T}({\cal T}_{h}),$

where

(15) $\displaystyle\mathcal{A}(v,{\bm{\sigma}};w,{\bm{\tau}})$
$\displaystyle=$
$\displaystyle\int_{\mathcal{F}_{h}^{\text{space}}}\big{(}c^{-2}v^{-}\llbracket
w\rrbracket_{t}+{\bm{\sigma}}^{-}\cdot\llbracket{\bm{\tau}}\rrbracket_{t}\big{)}~{}d{\bf
x}+\int_{\mathcal{F}_{h}^{\text{time}}}\bigg{(}\\{\\{v\\}\\}\llbracket{A^{\frac{1}{2}}\bm{\tau}}\rrbracket_{\bf
N}+\\{\\{{\bm{\sigma}}\\}\\}\cdot A^{\frac{1}{2}}\llbracket w\rrbracket_{\bf
N}$ (16) $\displaystyle+$ $\displaystyle\alpha A^{\delta}\llbracket
v\rrbracket_{\bf N}\cdot A^{\frac{1}{2}}\llbracket w\rrbracket_{\bf
N}+\beta\llbracket A^{\frac{1}{2}}{\bm{\sigma}}\rrbracket_{\bf N}\llbracket
A^{\frac{1}{2}}{\bm{\tau}}\rrbracket_{\bf N}\bigg{)}~{}dS$ (17)
$\displaystyle+$
$\displaystyle\int_{\mathcal{F}_{h}^{T}}\big{(}c^{-2}vw+{\bm{\sigma}}\cdot{\bm{\tau}}\big{)}~{}d{\bf
x}+\int_{\mathcal{F}_{h}^{D}}({\bm{\sigma}}\cdot wA^{\frac{1}{2}}{\bf
n}^{x}_{\Omega}+\alpha vwA^{\delta}{\bf n}^{x}_{\Omega}\cdot
A^{\frac{1}{2}}{\bf n}^{x}_{\Omega})~{}dS$ (18) $\displaystyle+$
$\displaystyle\int_{\mathcal{F}_{h}^{N}}\bigg{(}v(A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf
n}^{x}_{\Omega})+\beta(A^{\frac{1}{2}}{\bm{\sigma}}\cdot{\bf
n}^{x}_{\Omega})(A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf
n}^{x}_{\Omega})\bigg{)}~{}dS,$

and

(19) $\displaystyle\ell(w,{\bm{\tau}})$ $\displaystyle=$
$\displaystyle\int_{\mathcal{F}_{h}^{0}}\big{(}c^{-2}v_{0}w+{\bm{\sigma}}_{0}\cdot{\bm{\tau}}\big{)}~{}d{\bf
x}+\int_{\mathcal{F}_{h}^{N}}g_{N}\bigg{(}\beta
A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf n}^{x}_{\Omega}-w\bigg{)}~{}dS$ (20)
$\displaystyle+$ $\displaystyle\int_{\mathcal{F}_{h}^{D}}\alpha
g_{D}wA^{\delta}{\bf n}^{x}_{\Omega}\cdot A^{\frac{1}{2}}{\bf
n}^{x}_{\Omega}~{}dS-\int_{\mathcal{F}_{h}^{D}}g_{D}A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf
n}^{x}_{\Omega}~{}dS.$

######  Remark 3.1.

The choice of numerical fluxes $(\check{v},\check{\bm{\sigma}})$ defined on
``time-like" interfaces and boundary faces is such that the Trefftz DG
formulation is consistent; namely, if $(v,\bm{\sigma})\in H^{1}(Q)$ solves
(1), then it satisfies (14). In particular, the flux $\check{\bm{\sigma}}$
satisfies $A^{\frac{1}{2}}\check{\bm{\sigma}}\cdot{\bf n}^{x}_{\Omega}=g_{N}$
on ${\mathcal{F}}_{h}^{N}$ coinciding with the Neumann boundary condition
satisfied by the exact solution $(v,\bm{\sigma})$ of the IBVP (1).

##  4 Discretization of the variational problems

The proposed Trefftz DG method for (1) depends on two transformations.

###  4.1 A coordinate transformation and a scaled transformation

Since $A$ is positive definite matrix, there exists an orthogonal matrix $P$
and a diagonal positive definite matrix
$\Lambda=\text{diag}(\lambda_{1},\lambda_{2},\cdots,\lambda_{d})$ such that
$A=P^{T}\Lambda P$, where $\lambda_{i}\leq\lambda_{i+1}(1\leq i\leq d-1)$ and
the superscript $T$ denotes matrix transposition. Set
$\lambda_{\text{min}}=\lambda_{1}$ and $\lambda_{\text{max}}=\lambda_{d}$. Of
course, we can assume that $\text{det}(P)=1$. Define a coordinate
transformation:

(21) $\hat{\bf x}=\Lambda^{-\frac{1}{2}}P{\bf x}\xlongequal{\Delta}S{\bf
x},\quad S=\Lambda^{-\frac{1}{2}}P.$

Under the coordinate transformation (21), let $\hat{\Omega}$ and $\hat{Q}$
denote the images of $\Omega$ and $Q$, respectively, and denote by $\hat{\cal
T}_{{\hat{h}}_{\hat{\bf x}}}^{\hat{\bf x}}=\\{\hat{K}_{\hat{\bf x}}\\}$ the
transformed finite element mesh of the spatial domain $\hat{\Omega}$ with

$\hat{h}_{\hat{K}_{\hat{\bf x}}}=\text{diam}\hat{K}_{\hat{\bf
x}},~{}~{}\hat{h}_{\hat{\bf
x}}=\mathop{\text{missing}}{max}\limits_{\hat{K}_{\hat{\bf x}}\in\hat{\cal
T}_{{\hat{h}}_{\hat{\bf x}}}^{\hat{x}}}\hat{h}_{\hat{K}_{\hat{\bf x}}}.$

Furthermore, the transformed space-time domain
$\hat{Q}=\hat{\Omega}\times(0,T)$ can be partitioned with a finite element
mesh $\hat{\cal T}_{\hat{h}}$ given by

$\hat{\cal T}_{\hat{h}}=\\{\hat{K}=\hat{K}_{\hat{\bf x}}\times
I_{n},\hat{K}_{\hat{\bf x}}\in\hat{\cal T}_{\hat{h}_{\hat{\bf x}}}^{\hat{\bf
x}},1\leq n\leq N\\}.$

Assume that the space-time mesh satisfies the condition

(22) $\hat{h}=\mathop{\text{missing}}{max}\big{\\{}\hat{h}_{\hat{\bf
x}},\mathop{\text{missing}}{max}\limits_{n^{{}^{\prime}}=1,\cdots,n}ch_{n^{{}^{\prime}}}\big{\\}}\leq\hat{\rho}\mathop{\text{missing}}{min}\big{\\{}\hat{h}_{\hat{\bf
x}},\mathop{\text{missing}}{min}\limits_{n^{{}^{\prime}}=1,\cdots,n}ch_{n^{{}^{\prime}}}\big{\\}},$

for each discrete time $t_{n}$ and some $\hat{\rho}>1$.

We denote by $\hat{\mathcal{F}}_{\hat{h}}=\bigcup\limits_{\hat{K}\in{\hat{\cal
T}}_{\hat{h}}}\partial\hat{K}$ the skeleton of the mesh, by
$\hat{\mathcal{F}}_{\hat{h}}^{\text{space}}$ the union of the internal space-
like faces, and by $\hat{\mathcal{F}}_{\hat{h}}^{\text{time}}$ the union of
the internal time-like faces, respectively. We use $\hat{\bf n}$ to denote the
unit outer normal vector on the boundary of each element $\hat{K}_{\hat{\bf
x}}$. Denote by $\hat{\Gamma}_{D}$ and $\hat{\Gamma}_{N}$ the images of
$\Gamma_{D}$ and $\Gamma_{N}$ under the coordinate transformation (21),
respectively. Set
$\hat{\mathcal{F}}_{\hat{h}}^{0}=\hat{\Omega}\times\\{t=0\\}$,
$\hat{\mathcal{F}}_{\hat{h}}^{T}=\hat{\Omega}\times\\{t=T\\}$,
$\hat{\mathcal{F}}_{\hat{h}}^{D}=\hat{\Gamma}_{D}\times[0,T]$ and
$\hat{\mathcal{F}}_{\hat{h}}^{N}=\hat{\Gamma}_{N}\times[0,T]$.

Denote by $\nabla_{h}$ and $\hat{\nabla}_{\hat{h}}$ the element application of
the spacial gradient operator $\nabla=(\frac{\partial}{\partial
x_{1}}~{}\frac{\partial}{\partial x_{2}}~{}\cdots~{}\frac{\partial}{\partial
x_{d}})^{T}$ and
$\hat{\nabla}=(\frac{\partial}{\partial\hat{x}_{1}}~{}\frac{\partial}{\partial\hat{x}_{2}}~{}\cdots~{}\frac{\partial}{\partial\hat{x}_{d}})^{T}$,
respectively. Define the spacial Laplace operator $\hat{\triangle}$ on
$\hat{\Omega}$ by
$\hat{\triangle}=\frac{\partial^{2}}{\partial\hat{x}_{1}^{2}}+\frac{\partial^{2}}{\partial\hat{x}_{2}^{2}}+\cdots+\frac{\partial^{2}}{\partial\hat{x}_{d}^{2}}$.

Define the scaled fields $(\hat{v},\hat{\bm{\sigma}})$ as

(23) $\hat{v}=v,~{}~{}\hat{\bm{\sigma}}=P{\bm{\sigma}}.$

By some patient calculation in ``Appendix", we can obtain the following
relationships transforming ``anisotropic" into ``isotropic":

(24) $A^{\frac{1}{2}}\nabla
v=P^{T}\hat{\nabla}\hat{v},~{}~{}\text{and}~{}~{}\nabla\cdot(A^{\frac{1}{2}}{\bm{\sigma}})=\hat{\nabla}\cdot\hat{\bm{\sigma}}.$

Thus the anisotropic wave equation (1) is transformed into the isotropic wave
equation:

(25) $\left\\{\begin{aligned}
&\hat{\nabla}\hat{v}+\frac{\partial\hat{\bm{\sigma}}}{\partial t}={\bf
0}&\text{in}\quad\hat{Q},\\\
&\hat{\nabla}\cdot\hat{\bm{\sigma}}+c^{-2}\frac{\partial\hat{v}}{\partial
t}=\hat{f}=f&\text{in}\quad\hat{Q}.\\\ \end{aligned}\right.$

Conversely, if $(\hat{v},\hat{\bm{\sigma}})$ satisfies the isotropic isotropic
wave equation (25), $(v,{\bm{\sigma}})$ defined by the inverse scaled
transformation of (23) and the coordinate transformation (21) satisfies the
original anisotropic wave equation (1).

###  4.2 Anisotropic Trefftz basis function spaces

In order to derive a finite dimensional Trefftz space ${\bf V}_{h}({\cal
T}_{h})\subset{\bf T}({\cal T}_{h})$ satisfying the original anisotropic wave
equation (1), we first give the definition of a discretized Trefftz space
$\hat{\bf V}_{\hat{h}}(\hat{\cal T}_{\hat{h}})$ satisfying isotropic wave
equation (25). We refer the reader to [27, Remark 13] for a detailed
construction of $\hat{\bf V}_{\hat{h}}(\hat{\cal T}_{\hat{h}})$ by evolving in
time polynomial initial conditions.

Assuming that the first order problem (1) is derived from the second order
problem (2), and define the polynomial Trefftz space for the second order
problem from (25):

(26)
$\hat{\mathbb{U}}^{p}(\hat{K})=\\{\hat{U}\in\mathbb{P}^{p}(\hat{K})~{}s.t.~{}-\Delta\hat{U}+c^{-2}\frac{\partial^{2}\hat{U}}{\partial
t^{2}}=0\\},$

where the subscript $p$ is related to the dimension of the local spaces.
Denote some multi-index notation for ${\bm{\alpha}}\in\mathbb{N}_{0}^{d}$ by
$|{\bm{\alpha}}|=\alpha_{1}+\cdots+\alpha_{d}$,
$D^{\bm{\alpha}}\phi=\frac{\partial^{|{\bm{\alpha}}|}\phi}{\partial\hat{x}_{1}^{\alpha_{1}}\cdots\partial\hat{x}_{d}^{\alpha_{d}}}$,
${\bf
x}^{\bm{\alpha}}=\hat{x}_{1}^{\alpha_{1}}\cdots\hat{x}_{d}^{\alpha_{d}}$, and
for a space-time field $\phi$, by
$D^{\alpha_{t},\bm{\alpha}}\phi=\frac{\partial^{\alpha_{t}+|{\bm{\alpha}}|}\phi}{\partial
t^{\alpha_{t}}\partial\hat{x}_{1}^{\alpha_{1}}\cdots\partial\hat{x}_{d}^{\alpha_{d}}}$.
If the polynomial

$\hat{U}(t,\hat{\bf
x})=\sum\limits_{k\in\mathbb{N}_{0},{\bm{\alpha}}\in\mathbb{N}_{0}^{d},k+|{\bm{\alpha}}|\leq
p}a_{k,{\bm{\alpha}}}t^{k}\hat{\bf x}^{\bm{\alpha}}$

with $a_{k,{\bm{\alpha}}}\in\mathbb{R}$ satisfies the second order wave
equation in (26), then the coefficients $a_{k,{\bm{\alpha}}}$ satisfy the
recurrence

$a_{k,{\bm{\alpha}}}=\frac{c^{2}}{k(k-1)}\sum_{m=1}^{d}(\alpha_{m}+2)(\alpha_{m}+1)a_{k-2,{\bm{\alpha}}+2{\bf
e}_{m}},$

where ${\bf e}_{m}$ is the $m$th row of the identity matrix of order $d$.

In order to start the recursion, one can start by choosing polynomial basis
functions $\\{\tilde{b}_{1},\cdots,\tilde{b}_{C_{p+d}^{d}}\\}$ for the space
$\mathbb{P}^{p}(\mathbb{R}^{d})$ for $k=0$ and
$\\{\tilde{\tilde{b}}_{1},\cdots,\tilde{\tilde{b}}_{C_{p-1+d}^{d}}\\}$ for the
space $\mathbb{P}^{p-1}(\mathbb{R}^{d})$ for $k=1$. Then a basis for
$\hat{\mathbb{U}}^{p}(\hat{K})$ can be defined such that either
$U(\cdot,0)=\tilde{b}_{j}$ and $\frac{\partial U}{\partial t}(\cdot,0)=0$, or
$U(\cdot,0)=0$ and $\frac{\partial U}{\partial
t}(\cdot,0)=\tilde{\tilde{b}}_{j}$ for some $j$. It leads to the dimension of
$\hat{\mathbb{U}}^{p}(\hat{K})$ as

$\text{dim}\hat{\mathbb{U}}^{p}(\hat{K})=C_{p+d}^{d}+C_{p-1+d}^{d},$

where $C_{k}^{j}=\frac{k!}{j!(k-j)!}$ for $j\leq k\in\mathbb{N}_{0}$.

We denote by $\hat{b}_{j}(1\leq j\leq\text{dim}\hat{\mathbb{U}}^{p}(K))$ the
basis functions of the space $\hat{\mathbb{U}}^{p}(\hat{K})$. Then, a Trefftz
space $\hat{\mathbb{W}}^{p}(\hat{K})$ for the first order system can be
derived from

$\hat{\mathbb{W}}^{p}(\hat{K})=\text{span}\bigg{\\{}(\frac{\partial\hat{b}_{j}}{\partial
t},-\hat{\nabla}\hat{b}_{j}),~{}\hat{b}_{j}\in\hat{\mathbb{U}}^{p+1}(\hat{K}),~{}1\leq
j\leq\text{dim}\hat{\mathbb{U}}^{p+1}(\hat{K})\bigg{\\}}.$

Since the constants in $\hat{\mathbb{U}}^{p+1}(\hat{K})$ have no contribution
to $\hat{\mathbb{W}}^{p}(\hat{K})$, the dimension of
$\hat{\mathbb{W}}^{p}(\hat{K})$ equals
$\text{dim}\hat{\mathbb{U}}^{p+1}(\hat{K})-1$. Meanwhile we have the isotropic
Trefftz space defined on $\hat{\cal T}_{\hat{h}}$

$\hat{\bf V}_{\hat{h}}(\hat{\cal
T}_{\hat{h}})=\prod\limits_{\hat{K}\in\hat{\cal
T}_{\hat{h}}}\hat{\mathbb{W}}^{p}(\hat{K}).$

By the coordinate transformation (21) and (23), we get the anisotropic Trefftz
space ${\mathbb{W}}^{p}(K)$ on $K\in{\cal T}_{h}$ for the first two equations
of first order system (1)

(27)
${\mathbb{W}}^{p}(K)=\text{span}\bigg{\\{}(\frac{\partial\hat{b}_{j}}{\partial
t},-P^{T}\hat{\nabla}\hat{b}_{j}),~{}\hat{b}_{j}\in\hat{\mathbb{U}}^{p+1}(\hat{K}),~{}K=S^{-1}\hat{K},~{}1\leq
j\leq\text{dim}\hat{\mathbb{U}}^{p+1}(\hat{K})\bigg{\\}}.$

Furthermore, we get the anisotropic Trefftz space defined on ${\cal T}_{h}$

(28) ${\bf V}_{h}({\cal T}_{h})=\prod\limits_{K\in{\cal
T}_{h}}{\mathbb{W}}^{p}(K).$

Then, we can obtain the discretized Trefftz-DG variational formulation
corresponding to (14): Find $(v_{h},{\bm{\sigma}}_{h})\in{\bf V}_{h}({\cal
T}_{h})$ such that

(29)
$\mathcal{A}(v_{h},{\bm{\sigma}}_{h};w,{\bm{\tau}})=\ell(w,{\bm{\tau}})\quad\forall(w,{\bm{\tau}})\in{\bf
V}_{h}({\cal T}_{h}).$

##  5 Error estimates

In this Section we derive error estimates of approximations generated by the
global Trefftz DG method.

###  5.1 The required partition

In order to derive the desired error estimates of the approximate solutions,
we require that the partition must satisfy some assumptions. In this part we
introduce a kind of particular triangulation such that these assumptions can
be met.

We adopt a non-regularity triangulation ${\cal T}_{h_{{\bf x}}}^{\bf x}$ for
the three-dimensional domain $\Omega$ as follows (see Figure 1).

Mesh Generation Algorithm:

Step 1. Determine the transformed domain $\hat{\Omega}$ under the coordinate
transformation (21).

Step 2. Decompose $\hat{\Omega}$ into polyhedron elements
$\\{\hat{K}_{\hat{\bf x}}\\}$ such that $\hat{\cal T}_{\hat{h}_{\hat{\bf
x}}}^{\hat{\bf x}}$ is shape regular and quasi-uniform in the usual manner.

Step 3. Determine the triangulation ${\cal T}_{h_{{\bf x}}}^{\bf x}=\\{K_{\bf
x}\\}$ of $\Omega$ by using the inverse transformation of (21) acting on the
elements of $\hat{\cal T}_{\hat{h}_{\hat{\bf x}}}^{\hat{\bf x}}$.

Fig. 1: Mesh generation.

Under the proposed mesh triangulation, we introduce the two important
geometric properties, which are the direct generalization to multidimensional
space cases of Lemma 2.2 and 2.3 of [36], respectively.

######  Lemma 1.

For the proposed triangulation, we have

(30) $c_{0}||\Lambda^{\frac{1}{2}}||^{-1}h_{{\bf x}}\leq\hat{h}_{\hat{\bf
x}}\leq C_{0}||\Lambda^{\frac{1}{2}}||^{-1}h_{{\bf
x}},~{}~{}\text{and}~{}~{}c_{0}||\Lambda^{\frac{1}{2}}||^{-1}h\leq\hat{h}\leq
C_{0}||\Lambda^{\frac{1}{2}}||^{-1}h,$

where $c_{0}$ and $C_{0}$ denote two constants independent of $A$.

The next Lemma gives a relation between the areas of two bounded hyperplanes
based on the coordinate transformation (21).

######  Lemma 2.

For the proposed triangulation, denote by $\Gamma$ a hyperplane in
$\mathbb{R}^{d}$ which belongs to $\mathcal{F}_{h}$, and by $\hat{\Gamma}$ the
correspondingly transformed hyperplane which belongs to
$\hat{\mathcal{F}}_{\hat{h}}$. Then we have

(31)
$\frac{|\Gamma|}{|\hat{\Gamma}|}\leq\text{det}(\Lambda^{\frac{1}{2}})\lambda_{1}^{-\frac{1}{2}},$

where $|f|$ denotes the area of a bounded hyperplane $f$ in the
$d$-dimensional space.

For the simplicity of notation, let $\rho$ denote the condition number
$\text{cond}(A)$ of the anisotropic matrix $A$. Then
$\rho=\text{cond}(\Lambda)=\text{cond}^{2}(\Lambda^{\frac{1}{2}})=\text{cond}^{2}(S)$.
Without losing generality, we set $||A||=1$. Equivalently, the original model
(1) can be normalized such that $||A||=1$.

###  5.2 The transformation stability with respect to mesh-dependent norms

For the global Trefftz DG space ${\bf V}_{h}({\cal T}_{h})$, we define the
following DG norms (see [27]):

(32) $\displaystyle|(w,{\bm{\tau}})|_{\text{DG}(Q)}^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{2}||c^{-1}\llbracket
w\rrbracket_{t}||^{2}_{L^{2}(\mathcal{F}_{h}^{\text{space}})}+\frac{1}{2}||\llbracket{\bm{\tau}}\rrbracket_{t}||^{2}_{{L^{2}(\mathcal{F}_{h}^{\text{space}})}^{d}}+\frac{1}{2}||c^{-1}w||^{2}_{L^{2}(\mathcal{F}_{h}^{0}\cup\mathcal{F}_{h}^{T})}$
(33) $\displaystyle+$
$\displaystyle\frac{1}{2}||{\bm{\tau}}||^{2}_{L^{2}(\mathcal{F}_{h}^{0}\cup\mathcal{F}_{h}^{T})^{d}}+||\alpha^{\frac{1}{2}}A^{\frac{1}{4}+\frac{\delta}{2}}\llbracket
w\rrbracket_{\bf
N}||^{2}_{{L^{2}(\mathcal{F}_{h}^{\text{time}})}^{d}}+||\beta^{\frac{1}{2}}\llbracket
A^{\frac{1}{2}}{\bm{\tau}}\rrbracket_{\bf
N}||^{2}_{{L^{2}(\mathcal{F}_{h}^{\text{time}})}}$ (34) $\displaystyle+$
$\displaystyle\big{|}\big{|}\alpha^{\frac{1}{2}}wA^{\frac{1}{4}+\frac{\delta}{2}}{\bf
n}^{x}_{\Omega}\big{|}\big{|}^{2}_{L^{2}(\mathcal{F}_{h}^{D})^{d}}+||\beta^{\frac{1}{2}}A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf
n}^{x}_{\Omega}||^{2}_{{L^{2}(\mathcal{F}_{h}^{N})}},$

and

(35) $\displaystyle|(w,{\bm{\tau}})|_{\text{DG}(Q)^{+}}^{2}$ $\displaystyle=$
$\displaystyle|(w,{\bm{\tau}})|_{\text{DG}(Q)}^{2}+2||c^{-1}w^{-}||^{2}_{L^{2}(\mathcal{F}_{h}^{\text{space}})}+2||{\bm{\tau}}^{-}||^{2}_{{L^{2}(\mathcal{F}_{h}^{\text{space}})}^{d}}$
(36) $\displaystyle+$
$\displaystyle||\alpha^{-\frac{1}{2}}A^{\frac{1}{4}-\frac{\delta}{2}}\\{\\{\bm{\tau}\\}\\}||^{2}_{{L^{2}(\mathcal{F}_{h}^{\text{time}})^{d}}}+||\alpha^{-\frac{1}{2}}A^{\frac{1}{4}-\frac{\delta}{2}}\bm{\tau}||^{2}_{{L^{2}(\mathcal{F}_{h}^{D})^{d}}}$
(37) $\displaystyle+$
$\displaystyle||\beta^{-\frac{1}{2}}\\{\\{w\\}\\}||^{2}_{{L^{2}(\mathcal{F}_{h}^{\text{time}})}}+||\beta^{-\frac{1}{2}}w||^{2}_{{L^{2}(\mathcal{F}_{h}^{N})}}.$

In addition, we endow the space $\hat{\bf V}_{\hat{h}}(\hat{\cal
T}_{\hat{h}})$ with the norm

(38) $\displaystyle|(\hat{w},\hat{\bm{\tau}})|_{\text{DG}(\hat{Q})}^{2}$
$\displaystyle=$
$\displaystyle\frac{1}{2}||c^{-1}\llbracket\hat{w}\rrbracket_{t}||^{2}_{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{space}})}+\frac{1}{2}||\llbracket\hat{\bm{\tau}}\rrbracket_{t}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{space}})}^{d}}+\frac{1}{2}||c^{-1}\hat{w}||^{2}_{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{0}\cup\hat{\mathcal{F}}_{\hat{h}}^{T})}$
(39) $\displaystyle+$
$\displaystyle\frac{1}{2}||\hat{\bm{\tau}}||^{2}_{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{0}\cup\hat{\mathcal{F}}_{\hat{h}}^{T})^{d}}+||\alpha^{\frac{1}{2}}\llbracket\hat{w}\rrbracket_{\bf
N}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}})}^{d}}+||\beta^{\frac{1}{2}}\llbracket\hat{\bm{\tau}}\rrbracket_{\bf
N}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}})}}$ (40)
$\displaystyle+$
$\displaystyle||\alpha^{\frac{1}{2}}\hat{w}||^{2}_{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{D})}+||\beta^{\frac{1}{2}}\hat{\bm{\tau}}\cdot\hat{\bf
n}^{\hat{x}}_{\hat{\Omega}}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{N})}},$

and the augmented norm

(41) $\displaystyle|(\hat{w},\hat{\bm{\tau}})|_{\text{DG}(\hat{Q})^{+}}^{2}$
$\displaystyle=$
$\displaystyle|(\hat{w},\hat{\bm{\tau}})|_{\text{DG}(\hat{Q})}^{2}+2||c^{-1}\hat{w}^{-}||^{2}_{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{space}})}+2||\hat{\bm{\tau}}^{-}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{space}})}^{d}}$
(42) $\displaystyle+$
$\displaystyle||\alpha^{-\frac{1}{2}}\\{\\{\hat{\bm{\tau}}\\}\\}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}})^{d}}}+|\alpha^{-\frac{1}{2}}\hat{\bm{\tau}}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{D})^{d}}}$
(43) $\displaystyle+$
$\displaystyle||\beta^{-\frac{1}{2}}\\{\\{\hat{w}\\}\\}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}})}}|+||\beta^{-\frac{1}{2}}\hat{w}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{N})}}.$

Moreover, we denote the truncated seminorms
$|(w,{\bm{\tau}})|_{\text{DG}(Q_{n})}:=|(w\chi_{Q_{n}},{\bm{\tau}}\chi_{Q_{n}})|_{\text{DG}(Q)}$
and
$|(w,{\bm{\tau}})|_{\text{DG}(Q_{n})^{+}}:=|(w\chi_{Q_{n}},{\bm{\tau}}\chi_{Q_{n}})|_{\text{DG}(Q)^{+}}$,
respectively, where $\chi_{Q_{n}}$ is the characteristic function of $Q_{n}$.
The same definition is applied to
$|(\hat{w},\hat{\bm{\tau}})|_{\text{DG}(\hat{Q}_{n})}$ and
$|(\hat{w},\hat{\bm{\tau}})|_{\text{DG}(\hat{Q}_{n})^{+}}$. Similarly we
define the truncated bilinear form
$\mathcal{A}_{\text{DG}(Q_{n})}(v,\bm{\sigma};w,\bm{\tau})=\mathcal{A}(v\chi_{Q_{n}},\bm{\sigma}\chi_{Q_{n}};w\chi_{Q_{n}},\bm{\tau}\chi_{Q_{n}})$.

The following Lemma states the transformation stability with respect to mesh-
dependent norms, which indicates that, in order to obtain the optimal
stability estimates and error estimates with respect to $\rho$, the best
choice of $\delta$ is set to be $\frac{1}{2}$ throughout the rest paper.

######  Lemma 3.

For $(w_{h},{\bm{\tau}}_{h})\in{\bf V}_{h}({\cal T}_{h})$ and
$\delta=\frac{1}{2}$, we have

(44)
$\displaystyle\begin{split}|(w_{h},{\bm{\tau}}_{h})|_{\text{DG}(Q)}&\leq~{}\text{det}(\Lambda^{\frac{1}{4}})\lambda_{\text{min}}^{-\frac{1}{4}}~{}|(\hat{w}_{\hat{h}},\hat{\bm{\tau}}_{\hat{h}})|_{\text{DG}(\hat{Q})},\\\
|(w_{h},{\bm{\tau}}_{h})|_{\text{DG}(Q)^{+}}&\leq~{}\text{det}(\Lambda^{\frac{1}{4}})\lambda_{\text{min}}^{-\frac{1}{4}}~{}|(\hat{w}_{\hat{h}},\hat{\bm{\tau}}_{\hat{h}})|_{\text{DG}(\hat{Q})^{+}}.\end{split}$

Proof. We divide the proof into two steps.

Step 1: To estimate the terms of $|(w_{h},{\bm{\tau}}_{h})|_{\text{DG}(Q)}$ on
$\mathcal{F}_{h}^{\text{space}}\cup\mathcal{F}_{h}^{0}\cup\mathcal{F}_{h}^{T}$.

By the coordinate transformation (21) and direct calculation, we obtain

$\displaystyle\begin{split}&\frac{1}{2}||c^{-1}\llbracket
w\rrbracket_{t}||^{2}_{L^{2}(\mathcal{F}_{h}^{\text{space}})}+\frac{1}{2}||\llbracket{\bm{\tau}}\rrbracket_{t}||^{2}_{{L^{2}(\mathcal{F}_{h}^{\text{space}})}^{d}}+\frac{1}{2}||c^{-1}w||^{2}_{L^{2}(\mathcal{F}_{h}^{0}\cup\mathcal{F}_{h}^{T})}+\frac{1}{2}||{\bm{\tau}}||^{2}_{L^{2}(\mathcal{F}_{h}^{0}\cup\mathcal{F}_{h}^{T})^{d}}&\\\
&\leq
C\text{det}(S^{-1})~{}\bigg{(}\frac{1}{2}||c^{-1}\llbracket\hat{w}\rrbracket_{t}||^{2}_{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{space}})}+\frac{1}{2}||\llbracket\hat{\bm{\tau}}\rrbracket_{t}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{space}})}^{d}}+\frac{1}{2}||c^{-1}\hat{w}||^{2}_{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{0}\cup\hat{\mathcal{F}}_{\hat{h}}^{T})}+\frac{1}{2}||\hat{\bm{\tau}}||^{2}_{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{0}\cup\hat{\mathcal{F}}_{\hat{h}}^{T})^{d}}\bigg{)}.\end{split}$

Step 2: To estimate the terms of $|(w_{h},{\bm{\tau}}_{h})|_{\text{DG}(Q)}$ on
$\mathcal{F}_{h}^{\text{time}}\cup\mathcal{F}_{h}^{D}\cup\mathcal{F}_{h}^{N}$.
It is easy to check that (here ${\bf n}^{t}_{K}=0$)

(45) ${\bf n}^{\bf x}_{K}=|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{K}|~{}P^{T}\Lambda^{\frac{-T}{2}}\hat{\bf n}^{\hat{\bf x}}_{\hat{K}}.$

Combining

(46) $\displaystyle\begin{split}&A^{\frac{1}{4}+\frac{\delta}{2}}\llbracket
w\rrbracket_{\bf N}=|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{K}|~{}P^{T}\Lambda^{\frac{\delta}{2}-\frac{1}{4}}\llbracket\hat{w}\rrbracket_{\bf
N},~{}~{}A^{\frac{1}{4}+\frac{\delta}{2}}{\bf
n}^{x}_{\Omega}=|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{\Omega}|P^{T}\Lambda^{\frac{\delta}{2}-\frac{1}{4}}\hat{\bf
n}^{\hat{x}}_{\hat{\Omega}},\\\ &\llbracket
A^{\frac{1}{2}}{\bm{\tau}}\rrbracket_{\bf N}=|\Lambda^{\frac{1}{2}}P{\bf
n}^{\bf x}_{K}|~{}\llbracket\hat{\bm{\tau}}\rrbracket_{\bf
N},~{}~{}A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf
n}^{x}_{\Omega}=|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{\Omega}|~{}\hat{\bm{\tau}}\cdot\hat{\bf
n}^{\hat{x}}_{\hat{\Omega}},\end{split}$

with (31), we get

(47)
$\displaystyle\begin{split}&||\alpha^{\frac{1}{2}}A^{\frac{1}{4}+\frac{\delta}{2}}\llbracket
w\rrbracket_{\bf
N}||_{{L^{2}(\mathcal{F}_{h}^{\text{time}})}^{d}}\leq\text{det}(\Lambda^{\frac{1}{4}})\lambda_{\text{min}}^{-\frac{1}{4}}~{}||\Lambda^{\frac{\delta}{2}-\frac{1}{4}}||\cdot||\alpha^{\frac{1}{2}}\llbracket\hat{w}\rrbracket_{\bf
N}||_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}})}^{d}},\\\
&||\alpha^{-\frac{1}{2}}A^{\frac{1}{4}-\frac{\delta}{2}}\\{\\{\bm{\tau}\\}\\}||_{{L^{2}(\mathcal{F}_{h}^{\text{time}})^{d}}}\leq~{}\text{det}(\Lambda^{\frac{1}{4}})\lambda_{\text{min}}^{-\frac{1}{4}}~{}||\Lambda^{\frac{1}{4}-\frac{\delta}{2}}||\cdot||\\{\\{\hat{\bm{\tau}}\\}\\}||_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}})^{d}}},\\\
&||\beta^{\frac{1}{2}}\llbracket A^{\frac{1}{2}}{\bm{\tau}}\rrbracket_{\bf
N}||_{{L^{2}(\mathcal{F}_{h}^{\text{time}})}}\leq\text{det}(\Lambda^{\frac{1}{4}})\lambda_{\text{min}}^{-\frac{1}{4}}~{}||\beta^{\frac{1}{2}}\llbracket\hat{\bm{\tau}}\rrbracket_{\bf
N}||_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}})}},\\\
&\big{|}\big{|}\alpha^{\frac{1}{2}}wA^{\frac{1}{4}+\frac{\delta}{2}}{\bf
n}^{x}_{\Omega}\big{|}\big{|}^{2}_{L^{2}(\mathcal{F}_{h}^{D})^{d}}\leq\text{det}(\Lambda^{\frac{1}{4}})\lambda_{\text{min}}^{-\frac{1}{4}}~{}||\Lambda^{\frac{\delta}{2}-\frac{1}{4}}||\cdot||\hat{w}||_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{D})}},\\\
&||\beta^{\frac{1}{2}}A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf
n}^{x}_{\Omega}||^{2}_{{L^{2}(\mathcal{F}_{h}^{N})}}\leq\text{det}(\Lambda^{\frac{1}{4}})\lambda_{\text{min}}^{-\frac{1}{4}}~{}||\beta^{\frac{1}{2}}\hat{\bm{\tau}}\cdot\hat{\bf
n}^{\hat{x}}_{\hat{\Omega}}||_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{N})}},\\\
&||\beta^{-\frac{1}{2}}\\{\\{w\\}\\}||_{{L^{2}(\mathcal{F}_{h}^{\text{time}})}}\leq\text{det}(\Lambda^{\frac{1}{4}})\lambda_{\text{min}}^{-\frac{1}{4}}~{}||\beta^{-\frac{1}{2}}\\{\\{\hat{w}\\}\\}||^{2}_{{L^{2}(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}})}}.\end{split}$

Combining the two steps with the best choice $\delta=\frac{1}{2}$ yields the
desired results (44). $\Box$

###  5.3 Error estimates of Trefftz discontinuous Galerkin approximations

We prove existence and uniqueness of the Trefftz DG solution of (29) and the
bilinear form in (15) admits the following upper bounds. Throughout this
paper, $C$ denotes a generic positive constant that may have different values
in different occurrences, where $C$ depends on the mesh of $\hat{\cal
T}_{\hat{h}}$ and the shape of the elements.

######  Lemma 4.

There exists a unique solution $(v_{h},~{}{\bm{\sigma}}_{h})$ to (29);
moreover, for $\forall(v,{\bm{\sigma}}),~{}(w,{\bm{\tau}})\in{\bf V}_{h}({\cal
T}_{h})$ we have

(48)
$|(w,{\bm{\tau}})|_{\text{DG}(Q)}^{2}=\mathcal{A}((w,{\bm{\tau}});(w,{\bm{\tau}})),$
(49) $|\mathcal{A}((v,{\bm{\sigma}}),(w,{\bm{\tau}}))|\leq
2|(v,{\bm{\sigma}})|_{\text{DG}(Q)^{+}}|(w,{\bm{\tau}})|_{\text{DG}(Q)},$

and

(50) $|\mathcal{A}((v,{\bm{\sigma}}),(w,{\bm{\tau}}))|\leq
2|(v,{\bm{\sigma}})|_{\text{DG}(Q)}~{}|(w,{\bm{\tau}})|_{\text{DG}(Q)^{+}}.$

Proof. Provided that $(v,{\bm{\sigma}}),~{}(w,{\bm{\tau}})\in{\bf V}_{h}({\cal
T}_{h})$, local integration by parts permits us to rewrite the bilinear form
$\mathcal{A}((v,{\bm{\sigma}}),~{}(w,{\bm{\tau}}))$ as

(51)
$\displaystyle\mathcal{A}((v,{\bm{\sigma}}),~{}(w,{\bm{\tau}}))=-\int_{\mathcal{F}_{h}^{\text{space}}}\big{(}c^{-2}\llbracket
v\rrbracket_{t}w^{+}+\llbracket{\bm{\sigma}}\rrbracket_{t}\cdot{\bm{\tau}}^{+}\big{)}~{}d{\bf
x}+\int_{\mathcal{F}_{h}^{0}}\big{(}c^{-2}vw+{\bm{\sigma}}\cdot{\bm{\tau}}\big{)}~{}d{\bf
x}$ (52)
$\displaystyle+\int_{\mathcal{F}_{h}^{\text{time}}}\big{(}-A^{\frac{1}{2}}\llbracket
v\rrbracket_{\bf N}\cdot\\{\\{{\bm{\tau}}\\}\\}-\llbracket
A^{\frac{1}{2}}{\bm{\sigma}}\rrbracket_{\bf N}\\{\\{w\\}\\}+\alpha
A^{\delta}\llbracket v\rrbracket_{\bf N}\cdot A^{\frac{1}{2}}\llbracket
w\rrbracket_{\bf N}+\beta\llbracket
A^{\frac{1}{2}}{\bm{\sigma}}\rrbracket_{\bf N}\llbracket
A^{\frac{1}{2}}{\bm{\tau}}\rrbracket_{\bf N}\big{)}~{}dS$ (53)
$\displaystyle+\int_{\mathcal{F}_{h}^{D}}(-vA^{\frac{1}{2}}{\bf
n}^{x}_{\Omega}\cdot{\bm{\tau}}+\alpha vwA^{\delta}{\bf n}^{x}_{\Omega}\cdot
A^{\frac{1}{2}}{\bf n}^{x}_{\Omega})~{}dS$ (54)
$\displaystyle+\int_{\mathcal{F}_{h}^{N}}(-(A^{\frac{1}{2}}{\bm{\sigma}}\cdot{\bf
n}^{x}_{\Omega})w+\beta(A^{\frac{1}{2}}{\bm{\sigma}}\cdot{\bf
n}^{x}_{\Omega})(A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf n}^{x}_{\Omega}))~{}dS.$

By taking $(v,{\bm{\sigma}})=(w,{\bm{\tau}})$ and summing the two expressions
given in (15) and (51), we obtain (48).

$|\cdot|_{\text{DG}(Q)}$ and $|\cdot|_{\text{DG}(Q)^{+}}$ are only seminorms
on broken Sobolev spaces defined on the mesh ${\cal T}_{h}(Q)$, but are norms
on ${\bf T}({\cal T}_{h}(Q))$: indeed $|(w,{\bm{\tau}})|_{\text{DG}}=0$ for
$(w,{\bm{\tau}})\in{\bf T}({\cal T}_{h}(Q))$ implies that $(w,{\bm{\tau}})$ is
solution of the homogeneous IBVP (1) with zero initial and boundary
conditions, so $(w,{\bm{\tau}})={\bf 0}$ by the well-posedness of the IBVP
itself (see [27, Section 5.1] and [24, Lemma 4.1]. Thus the variational
formulation (14) has a unique discrete solution. Existence of the solution
follows from linearity of the problem and finite dimensionality.

By applying the Cauchy-Schwarz inequality to (15) and (51), respectively, we
obatin (49) and (50). $\Box$

######  Theorem 5.

Assume that $f\in\prod\limits_{K\in{\cal T}_{h}}C_{0}^{\infty}(K^{\ast})$,
$s=\text{min}\\{p,k_{t}-1,k_{\bf x}-1\\}$ and that the IBVP solution
$(v,{\bm{\sigma}})\in C^{k_{t}-1}(I;H^{k_{\bf x}+1}(\Omega))\times
C^{k_{t}}(I;H^{k_{\bf x}}(\Omega)^{d})$. Then we have, for $1\leq n\leq N$,

(55) $\displaystyle\frac{1}{2}$
$\displaystyle||c^{-1}(v-v_{h})||_{L^{2}(\Omega\times\\{t_{n}\\})}+\frac{1}{2}||\bm{\sigma}-\bm{\sigma}_{h}||_{L^{2}(\Omega\times\\{t_{n}\\})^{d}}\leq|(v,{\bm{\sigma}})-(v_{h},{\bm{\sigma}}_{h})|_{\text{DG}(Q_{n})^{+}}$
$\displaystyle\leq$ $\displaystyle
C\rho^{\frac{1}{4}}~{}h^{s+\frac{1}{2}}|(v,{\bm{\sigma}})|_{H^{s+1}(Q_{n})^{1+d}}.$

Proof. By (48), (14) and (29) applied to $Q_{n}$, we obtain, for
$\forall(w_{h},{\bm{\tau}}_{h})\in{\bf V}_{h}({\cal T}_{h}(Q_{n}))$,

(56)
$|(v,{\bm{\sigma}})-(v_{h},{\bm{\sigma}}_{h})|_{\text{DG}(Q_{n})}^{2}=\mathcal{A}_{\text{DG}(Q_{n})}((v,{\bm{\sigma}})-(v_{h},{\bm{\sigma}}_{h});(v,{\bm{\sigma}})-(w_{h},{\bm{\tau}}_{h})).$

Taking into account (50), we get the abstract error estimate:

(57)
$|(v,{\bm{\sigma}})-(v_{h},{\bm{\sigma}}_{h})|_{\text{DG}(Q_{n})}\leq\mathop{\text{inf}}\limits_{(w_{h},{\bm{\tau}}_{h})\in{\bf
V}_{h}({\cal
T}_{h}(Q_{n}))}|(v,{\bm{\sigma}})-(w_{h},{\bm{\tau}}_{h})|_{\text{DG}(Q_{n})^{+}}.$

By the existing approximation result Corollary 4 of [27], there exists
$\hat{Q}_{\hat{h}}(\hat{v},\hat{\bm{\sigma}})=(\hat{Q}_{\hat{h}}\hat{v},\hat{Q}_{\hat{h}}\hat{\bm{\sigma}})\in\hat{\bf
V}_{\hat{h}}(\hat{\cal T}_{\hat{h}}(\hat{Q}_{n}))$ such that,

(58)
$|(\hat{v},\hat{\bm{\sigma}})-\hat{Q}_{\hat{h}}(\hat{v},\hat{\bm{\sigma}})|_{H^{j}(\hat{K})^{1+d}}\leq
C\hat{h}^{s+1-j}|(\hat{v},\hat{\bm{\sigma}})|_{H^{s+1}(\hat{K})^{1+d}}.$

Using the inverse transformation of (23), set

(59)
$Q_{h}(v,\bm{\sigma})=(\hat{Q}_{\hat{h}}\hat{v},P^{T}\hat{Q}_{\hat{h}}\hat{\bm{\sigma}}).$

Using (57), (59) and (44), we obtain

(60)
$\displaystyle|(v,{\bm{\sigma}})-(v_{h},{\bm{\sigma}}_{h})|_{\text{DG}(Q_{n})}\leq|(v,{\bm{\sigma}})-Q_{h}(v,\bm{\sigma})|_{\text{DG}(Q_{n})^{+}}$
(61)
$\displaystyle\leq~{}\text{det}(\Lambda^{\frac{1}{4}})\lambda_{\text{min}}^{-\frac{1}{4}}~{}|(\hat{v},\hat{\bm{\sigma}})-\hat{Q}_{\hat{h}}(\hat{v},\hat{\bm{\sigma}})|_{\text{DG}(\hat{Q}_{n})^{+}}.$

By the transformation stability of Lemma 3, we only need to derive a bound of
$(\hat{v},\hat{\bm{\sigma}})-\hat{Q}_{\hat{h}}(\hat{v},\hat{\bm{\sigma}})$ in
terms of elementwise sums of traces, tracking the dependence on spatial and
temporal meshsizes. Taking into account the defintion (41), we get, for
$(\hat{w},\hat{\bm{\tau}})\in\hat{\bf T}(\hat{\cal
T}_{\hat{h}}(\hat{Q}_{n}))$,

(62)
$\displaystyle|(\hat{w},\hat{\bm{\tau}})|_{\text{DG}(\hat{Q}_{n})^{+}}^{2}$
$\displaystyle\leq$ $\displaystyle\sum\limits_{\hat{K}=\hat{K}_{\hat{x}}\times
I_{n^{\prime}}\in\hat{\cal
T}_{\hat{h}}(\hat{Q}_{n})}\bigg{[}||c^{-1}\hat{w}||^{2}_{L^{2}(\hat{K}_{\hat{x}}\times\\{t_{n^{\prime}-1},t_{n^{\prime}}\\})}+||\hat{\bm{\tau}}||^{2}_{L^{2}(\hat{K}_{\hat{x}}\times\\{t_{n^{\prime}-1},t_{n^{\prime}}\\})^{d}}$
(63) $\displaystyle+$
$\displaystyle\sum\limits_{F\in\partial\hat{K}\cap\hat{Q}_{n}\cap(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}}\cup\hat{\mathcal{F}}_{\hat{h}}^{D})}||\hat{w}\hat{\bf
n}^{\hat{x}}_{\hat{F}}||_{L^{2}(\hat{F})^{d}}^{2}+\sum\limits_{\hat{F}\in\partial\hat{K}\cap\hat{Q}_{n}\cap(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}}\cup\hat{\mathcal{F}}_{\hat{h}}^{N})}||\hat{w}||_{L^{2}(\hat{F})}^{2}$
(64) $\displaystyle+$
$\displaystyle\sum\limits_{\hat{F}\in\partial\hat{K}\cap\hat{Q}_{n}\cap(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}}\cup\hat{\mathcal{F}}_{\hat{h}}^{N})}||\hat{\bm{\tau}}\cdot\hat{\bf
n}^{\hat{x}}_{\hat{F}}||_{L^{2}(\hat{F})}^{2}+\sum\limits_{\hat{F}\in\partial\hat{K}\cap\hat{Q}_{n}\cap(\hat{\mathcal{F}}_{\hat{h}}^{\text{time}}\cup\hat{\mathcal{F}}_{\hat{h}}^{D})}||\hat{\bm{\tau}}||_{L^{2}(\hat{F})^{d}}^{2}\bigg{]}.$

By the standard weighted trace inequality applied in the time and space
directions independently (see [2, Sec 1.6.6]), the following bound holds true:

(65)
$\displaystyle|(\hat{w},\hat{\bm{\tau}})|_{\text{DG}(\hat{Q}_{n})^{+}}^{2}$
$\displaystyle\leq$ $\displaystyle
C\sum\limits_{\hat{K}=\hat{K}_{\hat{x}}\times I_{n^{\prime}}\in\hat{\cal
T}_{\hat{h}}(\hat{Q}_{n})}\bigg{[}\hat{h}^{-1}_{n^{\prime}}\bigg{(}||c^{-1}\hat{w}||^{2}_{L^{2}(\hat{K})}+||\hat{\bm{\tau}}||^{2}_{L^{2}(\hat{K})^{d}}\bigg{)}$
(66) $\displaystyle+$
$\displaystyle\hat{h}_{n^{\prime}}\bigg{(}|c^{-1}\hat{w}|^{2}_{H^{1}(I_{n^{\prime}};L^{2}(\hat{K}_{\hat{x}}))}+|\hat{\bm{\tau}}|^{2}_{H^{1}(I_{n^{\prime}};L^{2}(\hat{K}_{\hat{x}})^{d})}\bigg{)}$
(67) $\displaystyle+$
$\displaystyle\hat{h}_{\hat{K}_{\hat{x}}}^{-1}\bigg{(}||\hat{w}||^{2}_{L^{2}(\hat{K})}+||\hat{w}||^{2}_{L^{2}(\hat{K})}\bigg{)}+\hat{h}_{\hat{K}_{\hat{x}}}\bigg{(}|\hat{w}|^{2}_{L^{2}(I_{n^{\prime}};H^{1}(\hat{K}_{\hat{x}}))}+|\hat{w}|^{2}_{L^{2}(I_{n^{\prime}};H^{1}(\hat{K}_{\hat{x}}))}\bigg{)}$
(68) $\displaystyle+$
$\displaystyle\hat{h}_{\hat{K}_{\hat{x}}}^{-1}||\hat{\bm{\tau}}||^{2}_{L^{2}(\hat{K})^{d}}+\hat{h}_{\hat{K}_{\hat{x}}}|\hat{\bm{\tau}}|^{2}_{L^{2}(I_{n^{\prime}};H^{1}(\hat{K}_{\hat{x}})^{d})}\bigg{]}.$

Using (60), (65), (58), (30) and the scaling argument, we obtain

(69)
$\displaystyle|(v,{\bm{\sigma}})-(v_{h},{\bm{\sigma}}_{h})|_{\text{DG}(Q_{n})}\leq~{}\text{det}(\Lambda^{\frac{1}{4}})\lambda_{\text{min}}^{-\frac{1}{4}}\sum\limits_{\hat{K}=\hat{K}_{\hat{x}}\times
I_{n^{\prime}}\in\hat{\cal
T}_{\hat{h}}(\hat{Q}_{n})}\bigg{[}\hat{h}^{-\frac{1}{2}}||(\hat{v}^{(2)},\hat{\bm{\sigma}}^{(2)})-\hat{Q}_{\hat{h}}(\hat{v}^{(2)},\hat{\bm{\sigma}}^{(2)})||_{L^{2}(\hat{K})^{1+d}}$
(70)
$\displaystyle+\hat{h}^{\frac{1}{2}}|(\hat{v}^{(2)},\hat{\bm{\sigma}}^{(2)})-\hat{Q}_{\hat{h}}(\hat{v}^{(2)},\hat{\bm{\sigma}}^{(2)})|_{H^{1}(\hat{K})^{1+d}}+\hat{h}^{-\frac{1}{2}}\big{|}\big{|}\big{(}(\hat{v}^{(2)},\hat{\bm{\sigma}}^{(2)})-\hat{Q}_{\hat{h}}(\hat{v}^{(2)},\hat{\bm{\sigma}}^{(2)})\big{)}\big{|}\big{|}_{L^{2}(\hat{K})^{1+d}}$
(71)
$\displaystyle+\hat{h}\big{|}\big{(}(\hat{v}^{(2)},\hat{\bm{\sigma}}^{(2)})-\hat{Q}_{\hat{h}}(\hat{v}^{(2)},\hat{\bm{\sigma}}^{(2)})\big{)}\big{|}_{H^{1}(\hat{K})^{1+d}}^{2}\bigg{]}\leq
C\text{det}(\Lambda^{\frac{1}{4}})~{}\lambda_{\text{min}}^{-\frac{1}{4}}~{}\hat{h}^{s+\frac{1}{2}}\sum\limits_{\hat{K}=\hat{K}_{\hat{x}}\times
I_{n^{\prime}}\in\hat{\cal
T}_{\hat{h}}(\hat{Q}_{n})}|(\hat{v},\hat{\bm{\sigma}})|_{H^{s+1}(\hat{K})^{1+d}}$
(72) $\displaystyle\leq
C\rho^{\frac{1}{4}}~{}h^{s+{\frac{1}{2}}}|(v,{\bm{\sigma}})|_{H^{s+1}(Q_{n})^{1+d}}.\quad\Box$

##  6 A standard Trefftz DG method

A natural idea is to apply the standard Trefftz DG method to the isotropic
wave equation (25) derived by the coordinate transformation $S$ and the scaled
transformation (23), and then use the image of the resulting approximation
under the inverse transformation $S^{-1}$ and the inverse scaled
transformation (23) as the desired approximation of $(v,{\bm{\sigma}})$. We
will give the detailed derivation of variational formulation in this Section.

Let $(\hat{v},\hat{\bm{\sigma}})$ denote the analytic solution of the equation
(25) with the transformed boundary and initial conditions from the original
boundary and initial conditions (1):

(73) $\left\\{\begin{aligned}
&\hat{v}(\cdot,0)=v_{0},\quad\hat{\bm{\sigma}}(\cdot,0)=P{\bm{\sigma}}_{0}&\text{on}\quad\hat{\Omega},\\\
&\hat{v}=g_{D}&\text{on}\quad\hat{\Gamma}_{D}\times[0,T],\\\
&|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{\Omega}|\hat{\bm{\sigma}}\cdot\hat{\bf
n}^{x}_{\hat{\Omega}}=g_{N}&\text{on}\quad\hat{\Gamma}_{N}\times[0,T].\end{aligned}\right.$

Define the numerical fluxes as follows.

$\breve{v}=\left\\{\begin{aligned} &\hat{v}^{-}\\\ &\hat{v}\\\ &v_{0}\\\
&\\{\\{\hat{v}\\}\\}+\beta\llbracket\hat{\bm{\sigma}}\rrbracket_{N}\\\
&g_{D}\\\ &\hat{v}+\beta(|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{\Omega}|\hat{\bm{\sigma}}\cdot\hat{\bf
n}^{x}_{\hat{\Omega}}-g_{N})\end{aligned}\right.\breve{\bm{\sigma}}=\left\\{\begin{aligned}
&\hat{\bm{\sigma}}^{-}&\quad{\text{o}n}~{}\hat{\mathcal{F}}_{\hat{h}}^{\text{space}},\\\
&\hat{\bm{\sigma}}&\quad{\text{o}n}~{}\hat{\mathcal{F}}_{\hat{h}}^{T},\\\
&P{\bm{\sigma}}_{0}&\quad{\text{o}n}~{}\hat{\mathcal{F}}_{\hat{h}}^{0},\\\
&\\{\\{\hat{\bm{\sigma}}\\}\\}+\alpha\llbracket\hat{v}\rrbracket_{N}&\quad{\text{o}n}~{}\hat{\mathcal{F}}_{\hat{h}}^{\text{time}},\\\
&\hat{\bm{\sigma}}+\alpha(\hat{v}-g_{D})\hat{\bf
n}^{x}_{\hat{\Omega}}&\quad{\text{o}n}~{}\hat{\mathcal{F}}_{\hat{h}}^{D},\\\
&g_{N}|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf x}_{\Omega}|^{-1}\hat{\bf
n}^{x}_{\hat{\Omega}}&\quad{\text{o}n}~{}\hat{\mathcal{F}}_{\hat{h}}^{N}.\end{aligned}\right.$

Then $(\hat{v},\hat{\bm{\sigma}})\in\hat{\bf T}(\hat{\cal T}_{\hat{h}})$
satisfies

(74)
$\hat{\mathcal{A}}(\hat{v},\hat{\bm{\sigma}};\hat{w},\hat{\bm{\tau}})=\hat{\ell}(\hat{w},\hat{\bm{\tau}})\quad\forall(\hat{w},\hat{\bm{\tau}})\in\hat{\bf
T}(\hat{\cal T}_{\hat{h}}),$

where

(75) $\displaystyle\hat{\mathcal{A}}(v,{\bm{\sigma}};w,{\bm{\tau}})$
$\displaystyle=$
$\displaystyle\int_{\hat{\mathcal{F}}_{\hat{h}}^{\text{space}}}\big{(}c^{-2}v^{-}\llbracket
w\rrbracket_{t}+{\bm{\sigma}}^{-}\cdot\llbracket{\bm{\tau}}\rrbracket_{t}\big{)}~{}d\hat{\bf
x}+\int_{\hat{\mathcal{F}}_{\hat{h}}^{\text{time}}}\bigg{(}\\{\\{v\\}\\}\llbracket{\bm{\tau}}\rrbracket_{\bf
N}+\\{\\{{\bm{\sigma}}\\}\\}\cdot\llbracket w\rrbracket_{\bf N}$ (76)
$\displaystyle+$ $\displaystyle\alpha\llbracket v\rrbracket_{\bf
N}\cdot\llbracket w\rrbracket_{\bf
N}+\beta\llbracket{\bm{\sigma}}\rrbracket_{\bf
N}\llbracket{\bm{\tau}}\rrbracket_{\bf N}\bigg{)}~{}dS$ (77) $\displaystyle+$
$\displaystyle\int_{\hat{\mathcal{F}}_{\hat{h}}^{T}}\big{(}c^{-2}vw+{\bm{\sigma}}\cdot{\bm{\tau}}\big{)}~{}d\hat{\bf
x}+\int_{\hat{\mathcal{F}}_{\hat{h}}^{D}}({\bm{\sigma}}\cdot\hat{\bf
n}^{x}_{\hat{\Omega}}w+\alpha vw)~{}dS$ (78) $\displaystyle+$
$\displaystyle\int_{\hat{\mathcal{F}}_{\hat{h}}^{N}}\bigg{(}v({\bm{\tau}}\cdot\hat{\bf
n}^{x}_{\hat{\Omega}})+\beta|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{\Omega}|({\bm{\sigma}}\cdot\hat{\bf
n}^{x}_{\hat{\Omega}})({\bm{\tau}}\cdot\hat{\bf
n}^{x}_{\hat{\Omega}})\bigg{)}~{}dS,$

and

(79) $\displaystyle\hat{\ell}(w,{\bm{\tau}})$ $\displaystyle=$
$\displaystyle\int_{\hat{\mathcal{F}}_{\hat{h}}^{0}}\big{(}c^{-2}v_{0}w+P{\bm{\sigma}}_{0}\cdot{\bm{\tau}}\big{)}~{}d\hat{\bf
x}+\int_{\hat{\mathcal{F}}_{\hat{h}}^{N}}g_{N}\bigg{(}\beta{\bm{\tau}}\cdot\hat{\bf
n}^{x}_{\hat{\Omega}}-|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{\Omega}|^{-1}w\bigg{)}~{}dS$ (80) $\displaystyle+$
$\displaystyle\int_{\hat{\mathcal{F}}_{\hat{h}}^{D}}\alpha
g_{D}w~{}dS-\int_{\hat{\mathcal{F}}_{\hat{h}}^{D}}g_{D}{\bm{\tau}}\cdot\hat{\bf
n}^{x}_{\hat{\Omega}}~{}dS.$

######  Remark 6.1.

The choice of numerical fluxes $(\breve{v},\breve{\bm{\sigma}})$ defined on
``time-like" interfaces and boundary faces is such that the Trefftz DG
formulation is consistent; namely, if $(\hat{v},\hat{\bm{\sigma}})\in
H^{1}(\hat{Q})$ solves (25) with boundary and initial conditions (73), then it
satisfies (74). In particular, the flux $\breve{\bm{\sigma}}$ satisfies
$|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{\Omega}|\breve{\bm{\sigma}}\cdot\hat{\bf n}^{x}_{\hat{\Omega}}=g_{N}$
coinciding with the Neumann boundary condition satisfied by the exact solution
$(\hat{v},\hat{\bm{\sigma}})$.

Denote by $(\hat{v}_{\hat{h}},\hat{\bm{\sigma}}_{\hat{h}})$ the discrete
approximation of $(\hat{v},\hat{\bm{\sigma}})$, and let $(\tilde{v}_{h}({\bf
x}),\tilde{\bm{\sigma}}_{h}({\bf x}))$ denote the image of
$(\hat{v}_{\hat{h}},\hat{\bm{\sigma}}_{\hat{h}})$ under the inverse scaled
transformation of (23) and the coordinate transformation (21). By the
definition (15) of the sesquilinear form $\mathcal{A}(\cdot;\cdot)$, the
definition (75) of $\hat{\mathcal{A}}(\cdot;\cdot)$, (45), and (46), we have

(81) $\left\\{\begin{aligned} &c^{-2}v^{-}\llbracket w\rrbracket_{t}&=&\quad
c^{-2}\hat{v}^{-}\llbracket\hat{w}\rrbracket_{t}\\\
&{\bm{\sigma}}^{-}\cdot\llbracket{\bm{\tau}}\rrbracket_{t}&=&\quad\hat{\bm{\sigma}}^{-}\cdot\llbracket\hat{\bm{\tau}}\rrbracket_{t}\end{aligned}\right.\quad{\text{o}n}~{}\mathcal{F}_{h}^{\text{space}},$

and

(82) $\left\\{\begin{aligned}
&\\{\\{v\\}\\}\llbracket{A^{\frac{1}{2}}\bm{\tau}}\rrbracket_{\bf
N}&=&\quad|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{K}|~{}\\{\\{\hat{v}\\}\\}\llbracket\hat{\bm{\tau}}\rrbracket_{\bf N}\\\
&\\{\\{{\bm{\sigma}}\\}\\}\cdot\llbracket A^{\frac{1}{2}}w\rrbracket_{\bf
N}&=&\quad|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{K}|~{}\\{\\{\hat{\bm{\sigma}}\\}\\}~{}\llbracket\hat{w}\rrbracket_{\bf
N}\\\ &\alpha\llbracket A^{\frac{1}{2}}v\rrbracket_{\bf N}\cdot\llbracket
A^{\frac{1}{2}}w\rrbracket_{\bf N}&=&\quad\alpha|\Lambda^{\frac{1}{2}}P{\bf
n}^{\bf x}_{K}|^{2}~{}\llbracket\hat{v}\rrbracket_{\bf
N}\cdot\llbracket\hat{w}\rrbracket_{\bf N}\\\ &\beta\llbracket
A^{\frac{1}{2}}{\bm{\sigma}}\rrbracket_{\bf N}\llbracket
A^{\frac{1}{2}}{\bm{\tau}}\rrbracket_{\bf
N}&=&\quad\beta|\Lambda^{\frac{1}{2}}P{\bf n}^{\bf
x}_{K}|^{2}~{}\llbracket\hat{\bm{\sigma}}\rrbracket_{\bf
N}\llbracket\hat{\bm{\tau}}\rrbracket_{\bf
N}\end{aligned}\right.\quad{\text{o}n}~{}\mathcal{F}_{h}^{\text{time}}.$

Thus we can see that, no matter how the stabilization parameters $\alpha$ and
$\beta$ defined on $\mathcal{F}_{h}^{\text{time}}$ are chosen, the
sesquilinear form $\mathcal{A}(\cdot;\cdot)$ for the anisotropic case can not
coincide with $\mathcal{\hat{A}}_{\hat{h}}(\cdot;\cdot)$ in the sense of
proportionality. Thus the proposed approximation $({v}_{h}({\bf
x}),{\bm{\sigma}}_{h}({\bf x}))$ is different from the image
$(\tilde{v}_{h}({\bf x}),\tilde{\bm{\sigma}}_{h}({\bf x}))$. However, we still
build the desired error estimates of $(\tilde{v}_{h}({\bf
x}),\tilde{\bm{\sigma}}_{h}({\bf x}))$ as follows.

######  Theorem 6.

Assume that $f\in\prod\limits_{K\in{\cal T}_{h}}C_{0}^{\infty}(K^{\ast})$,
$s=\text{min}\\{p,k_{t}-1,k_{\bf x}-1\\}$ and that the IBVP solution
$(v,{\bm{\sigma}})\in C^{k_{t}-1}(I;H^{k_{\bf x}+1}(\Omega))\times
C^{k_{t}}(I;H^{k_{\bf x}}(\Omega)^{d})$. Then we have, for $1\leq n\leq N$,

(83) $\displaystyle\frac{1}{2}$
$\displaystyle||c^{-1}(v-\tilde{v}_{h})||_{L^{2}(\Omega\times\\{t_{n}\\})}+\frac{1}{2}||\bm{\sigma}-\tilde{\bm{\sigma}}_{h}||_{L^{2}(\Omega\times\\{t_{n}\\})^{d}}\leq|(v,{\bm{\sigma}})-(\tilde{v}_{h},\tilde{\bm{\sigma}}_{h})|_{\text{DG}(Q_{n})^{+}}$
$\displaystyle\leq$ $\displaystyle
C\rho^{\frac{1}{4}}~{}h^{s+\frac{1}{2}}|(v,{\bm{\sigma}})|_{H^{s+1}(Q_{n})^{1+d}}.$

Proof. By Lemma 3 and the scaling argument, we have

(84)
$|(v,{\bm{\sigma}})-(\tilde{v}_{h},\tilde{\bm{\sigma}}_{h})|_{\text{DG}(Q_{n})^{+}}\leq\text{det}(\Lambda^{\frac{1}{4}})~{}\lambda_{\text{min}}^{-\frac{1}{4}}|(\hat{v},\hat{\bm{\sigma}})-(\hat{v}_{\hat{h}},\hat{\bm{\sigma}}_{\hat{h}})|_{\text{DG}(\hat{Q}_{n})^{+}}.$

By (30), (58) and (62)-(69) for the isotropic case, we get

(85)
$\displaystyle|(v,{\bm{\sigma}})-(\tilde{v}_{h},\tilde{\bm{\sigma}}_{h})|_{\text{DG}(Q_{n})^{+}}\leq\text{det}(\Lambda^{\frac{1}{4}})~{}\lambda_{\text{min}}^{-\frac{1}{4}}~{}\hat{h}^{s+\frac{1}{2}}\sum\limits_{\hat{K}=\hat{K}_{\hat{x}}\times
I_{n^{\prime}}\in\hat{\cal
T}_{\hat{h}}(\hat{Q}_{n})}|(\hat{v},\hat{\bm{\sigma}})|_{H^{s+1}(\hat{K})^{1+d}}$
(86) $\displaystyle\leq
C\rho^{\frac{1}{4}}~{}h^{s+{\frac{1}{2}}}|(v,{\bm{\sigma}})|_{H^{s+1}(Q_{n})^{1+d}}.\quad\Box$

######  Remark 6.2.

By Theorems 5 and 6, we can see that the proposed method and the standard
Trefftz DG method have the same convergence order with respect to $h$ and
$\rho$, and almost have the same computational cost. Besides, we believe that
the orders of the condition number $\rho$ in the error estimates are optimal
since the transformation stability estimates seem sharp.

######  Remark 6.3.

We address that, compared with the valid error analysis of [1] only for the
homogeneous Neumann boundary conditions from the original IBVP and part valid
numerical results for nonhomogeneous Neumann boundary conditions, there is no
constraint for our proposed methods introduced in Sections 5 and 6 on the data
of homogeneous Neumann boundary conditions from the original IBVP, owing to
the use of trace estimate in (65) instead of the employed inverse trace
estimate as in [1, Section 6.1].

##  7 A nonhomogeneous model

The model reads as

(87) $\left\\{\begin{aligned} &A^{\frac{1}{2}}\nabla
v+\frac{\partial{\bm{\sigma}}}{\partial t}={\bf 0}&\text{in}\quad Q,\\\
&\nabla\cdot A^{\frac{1}{2}}{\bm{\sigma}}+c^{-2}\frac{\partial v}{\partial
t}=f&\text{in}\quad Q,\\\
&v(\cdot,0)=v_{0},\quad{\bm{\sigma}}(\cdot,0)={\bm{\sigma}}_{0}&\text{on}\quad\Omega,\\\
&v=g_{D}&\text{on}\quad\Gamma_{D}\times[0,T],\\\
&A^{\frac{1}{2}}{\bm{\sigma}}\cdot{\bf
n}^{x}_{\Omega}=g_{N}&\text{on}\quad\Gamma_{N}\times[0,T].\end{aligned}\right.$

In the framework of the global Trefftz DG method combined with overlapping
local DG method, we decompose the solution $(v,~{}{\bm{\sigma}})$ of the
problem $(\ref{nonmodel})$ into
$(v,~{}{\bm{\sigma}})=(v^{(1)},~{}{\bm{\sigma}}^{(1)})+(v^{(2)},~{}{\bm{\sigma}}^{(2)})$,
where $(v^{(1)},~{}{\bm{\sigma}}^{(1)})$ is a particular local solution of the
first two equations (87) on each fictitious domain with homogeneous boundary
and initial conditions, and $(v^{(2)},~{}{\bm{\sigma}}^{(2)})$ satisfies the
locally homogeneous wave equation.

Similarly to the derivation of (14), we can obtain the Trefftz-DG variational
formulation: Find $(v^{(2)},{\bm{\sigma}}^{(2)})\in{\bf T}({\cal T}_{h})$ such
that

(88)
$\mathcal{A}(v^{(2)},{\bm{\sigma}}^{(2)};w,{\bm{\tau}})=\tilde{\ell}(w,{\bm{\tau}})-\mathcal{A}(v^{(1)},{\bm{\sigma}}^{(1)};w,{\bm{\tau}})\quad\forall(w,{\bm{\tau}})\in{\bf
T}({\cal T}_{h}),$

where $\mathcal{A}(\cdot;\cdot)$ is defined as (15), and

(89) $\displaystyle\tilde{\ell}(w,{\bm{\tau}})$ $\displaystyle=$
$\displaystyle\int_{Q}fw~{}dV+\int_{\mathcal{F}_{h}^{0}}\big{(}c^{-2}v_{0}w+{\bm{\sigma}}_{0}\cdot{\bm{\tau}}\big{)}~{}d{\bf
x}$ (90) $\displaystyle+$ $\displaystyle\int_{\mathcal{F}_{h}^{D}}\alpha
g_{D}wA^{\delta}{\bf n}^{x}_{\Omega}\cdot A^{\frac{1}{2}}{\bf
n}^{x}_{\Omega}~{}dS-\int_{\mathcal{F}_{h}^{D}}g_{D}A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf
n}^{x}_{\Omega}~{}dS$ (91) $\displaystyle+$
$\displaystyle\int_{\mathcal{F}_{h}^{N}}g_{N}\bigg{(}\beta
A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf n}^{x}_{\Omega}-w\bigg{)}~{}dS.$

###  7.1 Nonhomogeneous local problems

For each space-time element $K=K_{\bf x}\times I_{n}\in{\cal T}_{h},K_{\bf
x}\in{\cal T}_{h_{{\bf x}}}^{\bf x}$, let $K_{\bf x}^{\ast}$ be a fictitious
domain that contains $K_{\bf x}$ as its subdomain. Set the fictitious domain
$K^{\ast}=K_{\bf x}^{\ast}\times I_{n}$,
$\mathcal{F}_{K^{\ast}}^{t_{n}}=K_{\bf x}^{\ast}\times\\{t=t_{n}\\}$, and
$\mathcal{F}_{K^{\ast},n}^{D}=\partial K_{\bf x}^{\ast}\times I_{n}$.

Define ${\bf V}^{(1)}_{K^{\ast}}=H^{1}(I_{n};H^{2,2}(K_{\bf x}^{\ast}))\times
H^{1}(I_{n};H^{1,1}(K_{\bf x}^{\ast})^{d})$. The particular solution
$(v^{(1)},~{}{\bm{\sigma}}^{(1)})\in(L^{2}(Q))^{1+d}$ is defined as
$(v^{(1)},~{}{\bm{\sigma}}^{(1)})\big{|}_{K}=(v^{(1)}_{K},~{}{\bm{\sigma}}^{(1)}_{K})=(v^{(1)}_{K^{\ast}},~{}{\bm{\sigma}}^{(1)}_{K^{\ast}})\big{|}_{K}$,
where $(v^{(1)}_{K^{\ast}},~{}{\bm{\sigma}}^{(1)}_{K^{\ast}})\in{\bf
V}^{(1)}_{K^{\ast}}$ satisfies the nonhomogeneous local acoustic equation on
the fictitious domain $K^{\ast}$:

(92) $\left\\{\begin{aligned} &A^{\frac{1}{2}}\nabla
v_{K^{\ast}}^{(1)}+\frac{\partial{\bm{\sigma}}_{K^{\ast}}^{(1)}}{\partial
t}={\bf 0}&\text{in}\quad K^{\ast},\\\ &\nabla\cdot
A^{\frac{1}{2}}{\bm{\sigma}}_{K^{\ast}}^{(1)}+c^{-2}\frac{\partial
v_{K^{\ast}}^{(1)}}{\partial t}=f&\text{in}\quad K^{\ast},\\\
&v_{K^{\ast}}^{(1)}=0,~{}~{}{\bm{\sigma}}_{K^{\ast}}^{(1)}={\bf
0}&\quad\quad\quad\text{on}\quad K_{\bf
x}^{\ast}\times\\{t=t_{n-1}\\},\hfill\\\
&v_{K^{\ast}}^{(1)}=0&\text{on}\quad\partial K_{\bf x}^{\ast}\times I_{n}.\\\
\end{aligned}\right.$

Similarly to the derivation of (11)-(13), we can obtain

(93) $\displaystyle-\int_{K^{\ast}}$ $\displaystyle\bigg{(}v(\nabla\cdot
A^{\frac{1}{2}}{\bm{\tau}}+c^{-2}\frac{\partial w}{\partial
t})+{\bm{\sigma}}\cdot(A^{\frac{1}{2}}\nabla
w+\frac{\partial\bm{\tau}}{\partial
t})\bigg{)}dV+\int_{\partial{K^{\ast}}}\bigg{(}\check{v}(A^{\frac{1}{2}}{\bm{\tau}}\cdot{\bf
n}^{x}_{K}+c^{-2}wn^{t}_{K})$ (94) $\displaystyle+$
$\displaystyle\check{\bm{\sigma}}\cdot(wA^{\frac{1}{2}}{\bf
n}^{x}_{K}+{\bm{\tau}}n^{t}_{K})\bigg{)}dS=\int_{K^{\ast}}f~{}w~{}dV,~{}\forall(w,{\bm{\tau}})\in{\bf
V}^{(1)}_{K^{\ast}}.$

Define the numerical fluxes as follows.

$\check{v}=\left\\{\begin{aligned} &v\\\ &0\\\ &0\\\
\end{aligned}\right.\quad\quad\check{\bm{\sigma}}=\left\\{\begin{aligned}
&{\bm{\sigma}}&\quad{\text{o}n}~{}\mathcal{F}_{K^{\ast}}^{t_{n}},\\\ &{\bf
0}&\quad{\text{o}n}~{}\mathcal{F}_{K^{\ast}}^{t_{n-1}},\\\ &\bm{\sigma}+\alpha
vA^{\delta}{\bf
n}^{x}_{\Omega}&\quad{\text{o}n}~{}\mathcal{F}_{K^{\ast},n}^{D}.\\\
\end{aligned}\right.$

Using (93) and the defined fluxes, the variational problem of (92) is to: Find
$(v_{K^{\ast}}^{(1)},{\bm{\sigma}}_{K^{\ast}}^{(1)})\in{\bf
V}^{(1)}_{K^{\ast}}$ such that

(95)
$\mathcal{A}^{(1)}(v_{K^{\ast}}^{(1)},{\bm{\sigma}}_{K^{\ast}}^{(1)};w,{\bm{\tau}})=\ell^{(1)}(w,{\bm{\tau}})\quad\forall(w,{\bm{\tau}})\in{\bf
V}^{(1)}_{K^{\ast}},$

where

(96)
$\displaystyle\mathcal{A}^{(1)}(v,{\bm{\sigma}};w,{\bm{\tau}})=-\int_{K^{\ast}}\bigg{(}{\bm{\sigma}}\cdot(\frac{\partial{\bm{\tau}}}{\partial
t}+A^{\frac{1}{2}}\nabla w)+v(\nabla\cdot
A^{\frac{1}{2}}{\bm{\tau}}+c^{-2}\frac{\partial w}{\partial t})\bigg{)}~{}dV$
(97)
$\displaystyle+\int_{\mathcal{F}_{K^{\ast}}^{t_{n}}}\big{(}c^{-2}vw+{\bm{\sigma}}\cdot{\bm{\tau}}\big{)}~{}d{\bf
x}+\int_{\mathcal{F}_{K^{\ast},n}^{D}}\bigg{(}{\bm{\sigma}}\cdot
wA^{\frac{1}{2}}{\bf n}^{x}_{K^{\ast}}+\alpha vwA^{\delta}{\bf
n}^{x}_{K^{\ast}}\cdot A^{\frac{1}{2}}{\bf n}^{x}_{K^{\ast}}\bigg{)}~{}dS,$

and

$\ell^{(1)}(w,{\bm{\tau}})=\int_{K^{\ast}}fw~{}dV.$

###  7.2 Discretization of the variational problems

We decompose the discrete solution $(v_{h},~{}{\bm{\sigma}}_{h})$ of the
problem $(\ref{model})$ into
$(v_{h},~{}{\bm{\sigma}}_{h})=(v_{h}^{(1)},~{}{\bm{\sigma}}_{h}^{(1)})+(v_{h}^{(2)},~{}{\bm{\sigma}}_{h}^{(2)})$,
where $(v_{h}^{(1)},~{}{\bm{\sigma}}_{h}^{(1)})$ defined later is the discrete
solution of continuous variational formulations (95), and
$(v_{h}^{(2)},~{}{\bm{\sigma}}_{h}^{(2)})\in{\bf V}_{h}({\cal T}_{h})$ is the
discrete solution of continuous Trefftz DG variational formulation (88);
namely, find $(v_{h}^{(2)},~{}{\bm{\sigma}}_{h}^{(2)})\in{\bf V}_{h}({\cal
T}_{h})$ such that

(98)
$\mathcal{A}(v_{h}^{(2)},{\bm{\sigma}}_{h}^{(2)};w,{\bm{\tau}})=\tilde{\ell}(w,{\bm{\tau}})-\mathcal{A}(v_{h}^{(1)},{\bm{\sigma}}_{h}^{(1)};w,{\bm{\tau}})\quad\forall(w,{\bm{\tau}})\in{\bf
V}_{h}({\cal T}_{h}).$

Let $\mathbb{Q}_{q}(K^{\ast})$ denote the set of polynomials of the same
degree $q$ in each of the $n+1$ variables. Define ${\bf
V}^{(1)}_{h,K^{\ast}}=\mathbb{Q}_{q}(K^{\ast})^{1+d}$ and ${\bf
V}^{(1)}_{h,K}=\mathbb{Q}_{q}(K)^{1+d}$.

Then a discretized version of the continuous variational problem (95): Find
$(v^{(1)}_{h,K^{\ast}},~{}{\bm{\sigma}}^{(1)}_{h,K^{\ast}})\in{\bf
V}^{(1)}_{h,K^{\ast}}$ such that

(99)
$\mathcal{A}^{(1)}(v_{h,K^{\ast}}^{(1)},{\bm{\sigma}}_{h,K^{\ast}}^{(1)};w,{\bm{\tau}})=\ell^{(1)}(w,{\bm{\tau}})\quad\forall(w,{\bm{\tau}})\in{\bf
V}^{(1)}_{h,K^{\ast}}.$

Define $(v_{h}^{(1)},~{}{\bm{\sigma}}_{h}^{(1)})\in\prod_{K\in{\cal
T}_{h}}{\bf V}^{(1)}_{h,K}$ by
$(v_{h}^{(1)},~{}{\bm{\sigma}}_{h}^{(1)})|_{K}=(v^{(1)}_{h,K^{\ast}},~{}{\bm{\sigma}}^{(1)}_{h,K^{\ast}})|_{K}$.

Determine each local fictitious domain $K_{\bf x}^{\ast}$ by using the inverse
transformation of (21) acting on $\hat{K}_{\hat{\bf x}}^{\ast}$. A natural way
is to choose $\hat{K}_{\hat{\bf x}}^{\ast}$ as the geometric sphere, e.g. the
disc for the two-dimensional case and the sphere for the three-dimensional
case, whose radius and center are denoted by $r_{\hat{K}_{\hat{\bf x}}}$ and
$O_{\hat{K}_{\hat{\bf x}}}$, respectively. Notice that the center and the
radius can be calculated easily. Then the variational problems (99) can be
solved easily by using the polar coordinate transformation for the calculation
of the involved integrations. We would like to emphasize that the discrete
problems (99) are local and independent each other for $K\in{\cal T}_{h}$, so
they can be explicitly solved in parallel and the cost is small for low
dimensional space case.

######  Remark 7.1.

We would like to point out, there is no constraint for our proposed method for
the nonhomogeneous model on the data of homogeneous Neumann boundary
conditions from the original IBVP, owing to the fact that the analytic
solution of the artificially constructed local IBVP (92) automatically
satisfies the homogeneous Neumann boundary condition (Note that
$\mathcal{F}_{K^{\ast},n}^{N}=\emptyset$), which is necessary in the
theoretical error analysis of the DG method [1].

######  Remark 7.2.

As stated in [17], if the nonhomogeneous local problem is defined on each
element $K$, which is a non-smooth domain, then the analytic solution has only
low regularity even if the analytic solution of the original problem defined
on the global solution domain is smooth enough. Surprisingly, if we choose
another alternative to define nonhomogeneous local problems on each
nonoverlapping time slab $D_{n}(1\leq n\leq N)$, the combined DG scheme can
also preserve the same orders of convergence as for the local smooth
overlapping case (see Section 8.4).

##  8 Numerical experiments

In this Section, we apply the proposed methods to solve the wave propagation
in anisotropic media, and we report numerical results to verify the efficiency
of the method.

The wave speed is fixed at $c=1$. As described in Section 4, we choose the
same number $p$ of basis functions for every elements $\Omega_{k}$, and
consider the following choice of numerical fluxes for the proposed methods:
the constant parameters $\alpha=\beta=1$. Meanwhile, we compare numerical
performances among the Trefftz methods and the high-order DG method introduced
in [1]. The new proposed assumption on the shape regularity of polyhedral
meshes $\hat{{\cal T}}_{\hat{h}_{\hat{\bf x}}}$ is employed for the high-order
DG method. For the convenience of statement, we call the method (14) as
``Method-I", and (74) described in Section 6 as ``Method-II".

In our tests, we estimate the convergence orders of the approximations by
measuring the relative errors in $L^{2}(\Omega)$-norm at a given time $T$, and
the errors in $|\cdot|_{\text{DG}}-$norm, respectively. All of the
computations have been done in MATLAB, and the system matrix was computed by
numerical integration. ``DOFs" represents the number of degree of freedoms
equal to the elements multiplied by the number of basis functions per element.

###  8.1 Homogeneous case

We use uniform mesh with $h_{{\bf x}}\approx h_{t}=2^{-l},l\in\mathbb{N}$. In
the Tables 8.1.1-8.1.2, the collum labelled ``Error" shows the numerical
errors in relative $L^{2}(\Omega\times\\{T\\})-$norm and
$|\cdot|_{\text{DG}}-$norm: if $p=1$, it is given for the mesh level $l=4$,
else if $p\in\\{2,3\\}$ for $l=3$.

####  8.1.1 2D space case with Neumann boundary conditions

We choose the space-time domain $Q=(0,1)^{2}\times(0,1)$, and set the
anisotropic matrix
$A=\left(\begin{array}[]{cc}\lambda_{1}^{2}a^{2}+\lambda_{2}^{2}b^{2}&ab(\lambda_{2}-\lambda_{1})\\\
ab(\lambda_{2}-\lambda_{1})&\lambda_{1}^{2}b^{2}+\lambda_{2}^{2}a^{2}\\\
\end{array}\right)$, where $a=\frac{1}{\sqrt{2}},~{}b=\frac{1}{\sqrt{2}}$, and
$0<\lambda_{1}<1,\lambda_{2}=1$. Consider the exact smooth solution

(100) $\displaystyle\hat{U}(\hat{\bf x},t)$ $\displaystyle=$
$\displaystyle\text{sin}(\pi\hat{x}_{1})\text{sin}(\pi\hat{x}_{2})\text{sin}(\sqrt{2}\pi
t),$ (101) $\displaystyle v$ $\displaystyle=$
$\displaystyle\hat{v}=\frac{\partial\hat{U}}{\partial
t},~{}{\bm{\sigma}}=P^{T}\hat{\bm{\sigma}}=-P^{T}\hat{\nabla}\hat{U},$

to the IBVP (1), with nonhomogeneous Neumann boundary conditions.

The convergence rates with respect to $h$ are given in the Table 8.1.1.

Table 1:

Convergence rates of the TDG schemes w.r.t. $h$.  $v_{h}$ $\bm{\sigma}_{h}$
$|\cdot|_{\text{DG}}$ Method-I $\rho$ Error Rate Error Rate Error Rate $p=1$ 2
4.70e-2 2.30 8.27e-3 2.41 2.02e-1 1.49 4 4.31e-2 2.31 8.25e-3 2.39 1.93e-1
1.49 16 4.60e-2 2.30 8.45e-3 2.33 1.94e-1 1.48 $p=2$ 2 1.16e-2 3.38 2.56e-3
3.17 5.49e-2 2.51 4 1.04e-2 3.31 2.51e-3 3.16 5.03e-2 2.52 16 1.18e-2 3.25
2.86e-3 3.14 5.75e-2 2.58 $p=3$ 2 7.55e-4 4.12 2.51e-4 3.90 3.83e-3 3.59 4
6.77e-4 4.10 2.39e-4 3.92 3.41e-3 3.59 16 7.96e-4 4.06 2.87e-4 3.96 3.33e-3
3.55 Method-II $\rho$ Error Rate Error Rate Error Rate $p=1$ 2 4.85e-2 2.30
8.68e-3 2.43 2.01e-1 1.49 4 4.52e-2 2.30 8.88e-3 2.43 1.91e-1 1.49 16 5.02e-2
2.31 9.24e-3 2.41 1.91e-1 1.48 $p=2$ 2 1.20e-2 3.44 2.58e-3 3.18 5.47e-2 2.52
4 1.10e-2 3.40 2.53e-3 3.18 4.99e-2 2.52 16 1.27e-2 3.43 2.83e-3 3.21 5.17e-2
2.51 $p=3$ 2 7.69e-4 4.13 2.51e-4 3.89 3.82e-3 3.49 4 7.02e-4 4.13 2.38e-4
3.89 3.39e-3 3.50 16 8.37e-4 4.11 2.84e-4 3.89 3.59e-3 3.49

We can obtain that $||v-v_{h}||_{L^{2}(\Omega\times
T)}\approx||\bm{\sigma}-\bm{\sigma}_{h}||_{L^{2}(\Omega\times
T)^{2}}=\mathcal{O}(h^{p+1})$. The last column showing the experimental
convergence rates of the errors measured in $|\cdot|_{\text{DG}}-$norm
indicates that the estimates of Theorems 5 and 6 are sharp.

The convergence rates with respect to $\rho$ are given in the Table 8.1.1.

Table 2:

Convergence rates of the TDG schemes w.r.t. $\rho$.  $v_{h}$ $\bm{\sigma}_{h}$
$|\cdot|_{\text{DG}}$ Method-I $\rho$ Error Rate Error Rate Error Rate $p=1$
$32$ 4.61e-2 8.12e-3 1.92e-1 64 4.84e-2 0.0702 8.53e-3 0.0711 1.96e-1 0.0297
128 5.01e-2 0.0498 8.89e-3 0.0596 1.97e-1 0.0073 $p=2$ $32$ 1.15e-2 2.94e-3
5.16e-2 64 1.28e-2 0.1545 3.35e-3 0.1883 5.43e-2 0.0736 128 1.36e-2 0.0875
3.77e-3 0.1704 5.62e-2 0.0496 $p=3$ 8 6.87e-4 2.54e-4 3.37e-3 16 7.96e-4
0.2125 2.87e-4 0.1762 3.61e-3 0.0993 $32$ 7.89e-4 -0.0127 2.96e-4 0.0445
3.54e-3 -0.0282 Method-II $\rho$ Error Rate Error Rate Error Rate $p=1$ $32$
5.06e-2 8.95e-3 1.89e-1 64 5.35e-2 0.0804 9.52e-3 0.0891 1.92e-1 0.0227 128
5.58e-2 0.0607 9.99e-3 0.0695 1.94e-1 0.0150 $p=2$ $32$ 1.25e-2 2.88e-3
5.08e-2 64 1.40e-2 0.1635 3.06e-3 0.0875 5.34e-2 0.0720 128 1.49e-2 0.0899
3.21e-3 0.0690 5.52e-2 0.0478 $p=3$ 8 7.20e-4 2.52e-4 3.35e-3 16 8.37e-4
0.2172 2.84e-4 0.1725 3.59e-3 0.0998 $32$ 8.36e-4 -0.0017 2.94e-4 0.0499
3.52e-3 -0.0284

We can obtain that $||v-v_{h}||_{L^{2}(\Omega\times
T)},||\bm{\sigma}-\bm{\sigma}_{h}||_{L^{2}(\Omega\times
T)^{2}},|(v,{\bm{\sigma}})-(v_{h},{\bm{\sigma}}_{h})|_{\text{DG}^{+}}\lesssim
C(v,\bm{\sigma})\mathcal{O}(\rho^{\frac{1}{4}})$, which indicates that the
estimates of Theorems 5 and 6 are sharp.

####  8.1.2 3D space case with Neumann boundary conditions

We choose the space-time domain $Q=(0,1)^{2}\times(0,1)$, and set the
anisotropic matrix
$A=\left(\begin{array}[]{ccc}\lambda_{1}^{2}a^{2}+\lambda_{2}^{2}b^{2}&ab(\lambda_{2}-\lambda_{1})&0\\\
ab(\lambda_{2}-\lambda_{1})&\lambda_{1}^{2}b^{2}+\lambda_{2}^{2}a^{2}&0\\\
0&0&1\\\ \end{array}\right)$, where
$a=\frac{1}{\sqrt{2}},~{}b=\frac{1}{\sqrt{2}}$, and
$0<\lambda_{1},\lambda_{2}<1$. Choose the space-time domain
$Q=(0,1)^{3}\times(0,1)$, and consider the exact smooth solution

(102) $\displaystyle\hat{U}(\hat{\bf x},t)$ $\displaystyle=$
$\displaystyle\text{sin}(\pi\hat{x}_{1})\text{sin}(\pi\hat{x}_{2})\text{sin}(\pi\hat{x}_{3})\text{sin}(\sqrt{3}\pi
t),$ (103) $\displaystyle v$ $\displaystyle=$
$\displaystyle\hat{v}=\frac{\partial\hat{U}}{\partial
t},~{}{\bm{\sigma}}=P^{T}\hat{\bm{\sigma}}=-P^{T}\hat{\nabla}\hat{U},$

to the IBVP (1), with nonhomogeneous Neumann boundary conditions.

The convergence rates with respect to $h$ are given in the Table 8.1.2. The
convergence rates with respect to $\rho$ are given in the Table 8.1.2.

Table 3:

Convergence rates of the TDG schemes w.r.t. $h$.  $v_{h}$ $\bm{\sigma}_{h}$
$|\cdot|_{\text{DG}}$ ($p,\rho$) Method Error Rate Error Rate Error Rate (1,2)
-I 2.41e-2 2.35 7.19e-3 2.53 1.68e-1 1.57 -II 2.39e-2 2.37 7.45e-3 2.54
1.62e-1 1.59 (2,2) -I 1.03e-2 3.41 7.54e-3 3.26 8.16e-2 2.51 -II 1.08e-2 3.41
7.78e-3 3.25 8.02e-2 2.53 (3,2) -I 7.76e-4 4.12 8.45e-4 4.05 8.62e-3 3.52 -II
7.94e-4 4.11 8.75e-4 4.02 8.85e-3 3.48

Table 4:

Convergence rates of the TDG schemes w.r.t. $\rho$.  $v_{h}$ $\bm{\sigma}_{h}$
$|\cdot|_{\text{DG}}$ Method-I $\rho$ Error Rate Error Rate Error Rate $p=1$
$4$ 1.53e-2 7.83e-3 1.05e-1 8 1.48e-2 -0.0479 8.02e-3 0.0346 1.06e-1 0.0137 16
1.48e-2 0 8.58e-3 0.0974 1.09e-1 0.0403 $p=2$ $4$ 9.67e-3 4.97e-3 4.63e-2 8
9.48e-3 -0.0286 4.78e-3 -0.0562 4.63e-2 0 16 9.73e-3 0.0376 4.88e-3 0.0299
4.59e-2 -0.0125 $p=3$ $4$ 9.18e-4 7.07e-4 5.34e-3 8 8.98e-4 -0.0318 7.03e-4
-0.0082 5.34e-3 0 16 9.57e-4 0.0918 7.46e-4 0.0857 5.45e-3 0.0294 Method-II
$\rho$ Error Rate Error Rate Error Rate $p=1$ $4$ 1.52e-2 8.14e-3 1.03e-1 8
1.45e-2 -0.0680 8.67e-3 0.0910 1.03e-1 0 16 1.44e-2 -0.0100 9.66e-3 0.1560
1.01e-1 -0.0283 $p=2$ $4$ 1.04e-2 4.86e-3 4.56e-2 8 1.06e-2 0.0275 4.66e-3
-0.0606 4.59e-2 0.0095 16 1.13e-2 0.0923 4.75e-3 0.0276 4.59e-2 0 $p=3$ $4$
9.92e-4 7.07e-4 5.34e-3 8 1.01e-3 0.0259 7.07e-4 0 5.36e-3 0.0054 16 1.10e-3
0.1231 7.42e-4 0.0697 5.51e-3 0.0398

We can obtain that $||v-v_{h}||_{L^{2}(\Omega\times
T)}\approx||\bm{\sigma}-\bm{\sigma}_{h}||_{L^{2}(\Omega\times
T)^{2}}=\mathcal{O}(h^{p+1})$, and $||v-v_{h}||_{L^{2}(\Omega\times
T)},||\bm{\sigma}-\bm{\sigma}_{h}||_{L^{2}(\Omega\times
T)^{2}},|(v,{\bm{\sigma}})-(v_{h},{\bm{\sigma}}_{h})|_{\text{DG}^{+}}\lesssim
C(v,\bm{\sigma})\mathcal{O}(\rho^{\frac{1}{4}})$, which indicate that the
estimates of Theorems 5 and 6 are sharp.

###  8.2 Nonhomogeneous case

We use uniform mesh with $h_{{\bf x}}\approx h_{t}=2^{-l},l\in\mathbb{N}$.

####  8.2.1 1D space case for smooth solution with Dirichlet boundary
conditions

We choose the space-time domain $Q=(0,1)\times(0,1)$. Consider the exact
smooth solution

(104) $\displaystyle U(x,t)$ $\displaystyle=$ $\displaystyle\text{sin}(\pi
x)\text{sin}(\sqrt{2}\pi t),$ (105) $\displaystyle v$ $\displaystyle=$
$\displaystyle\partial_{t}U,~{}\sigma=-\partial_{x}U,$

to the IBVP (1), with $A=1$, nonhomogeneous soure $f$:

$f=-\pi^{2}\text{sin}(\pi x)\text{sin}(\sqrt{2}\pi t),$

and homogeneous Dirichlet boundary conditions.

The convergence rates are given in the Table 8.2.1 for different choices of
$p$ and $q$.

Table 5:

Convergence rates of the combined space-time DG scheme w.r.t. $h$.  $v_{h}$
$\sigma_{h}$ $|\cdot|_{\text{DG}}$ ($p,q$) $h$ Error Rate Error Rate Error
Rate $(1,1)$ 1/8 5.32e-2 4.63e-2 3.64e-1 1/16 1.22e-2 2.24 8.70e-3 2.24
1.31e-1 1.47 1/32 2.98e-3 2.11 1.77e-3 2.11 4.66e-2 1.49 $(1,2)$ 1/8 5.19e-2
4.78e-2 3.30e-1 1/16 1.21e-2 2.22 9.08e-3 2.22 1.21e-1 1.45 1/32 2.98e-3 2.11
1.83e-3 2.11 4.31e-2 1.49 $(2,1)$ 1/8 9.55e-3 9.11e-3 1.51e-1 1/16 2.10e-3
2.33 1.58e-3 2.33 5.36e-2 1.50 1/32 5.04e-4 2.15 3.10e-4 2.15 1.90e-2 1.50
$(2,2)$ 1/4 2.70e-2 2.95e-2 1.07e-1 1/8 3.02e-3 3.20 3.17e-3 3.20 1.99e-2 2.52
1/16 3.50e-4 3.13 3.62e-4 3.13 3.57e-3 2.48 $(2,3)$ 1/4 2.72e-2 3.02e-2
9.52e-2 1/8 3.04e-3 3.21 3.24e-3 3.21 1.73e-2 2.46 1/16 3.55e-4 3.12 3.72e-4
3.12 3.72e-3 2.52 $(3,2)$ 1/4 4.57e-3 8.12e-3 6.10e-2 1/8 6.26e-4 3.21 7.99e-4
3.22 1.08e-2 2.50 1/16 8.45e-5 2.99 9.61e-5 2.99 1.91 2.50 $(3,3)$ 1/2 5.90e-2
3.59e-2 1.10e-1 1/4 2.69e-3 4.06 3.14e-3 4.04 9.02e-3 3.61 1/8 1.29e-4 4.28
1.68e-4 4.29 7.69e-4 3.55

From the Table 8.2.1, we obtain that, as typical for DG methods, in the case
of a regular enough solution, we observe the convergence orders of the errors
in $L^{2}(\Omega\times\\{T\\})-$norm with the rate
$\mathcal{O}(h^{\text{min}\\{p+1,q+1\\}})$, and that the convergence rates of
errors in $|\cdot|_{\text{DG}}-$norm are
$O(h^{\text{min}\\{p+\frac{1}{2},q+\frac{1}{2}\\}})$, which support
convergence rate optimality on the uniform refined meshes of the combined
numerical DG scheme just as for the homogeneous case. Furthermore, the Table
8.2.1 shows that for the $L^{2}(\Omega\times\\{T\\})$ norms and
$|\cdot|_{\text{DG}}$ norms of errors generated by the proposed method, the
choice $p=q+1$ is preferable compared with the choice $p=q$.

####  8.2.2 2D space case for smooth solution with Neumann boundary conditions

We choose the space-time domain $Q=(0,1)^{2}\times(0,1)$, and set the
anisotropic matrix
$A=\left(\begin{array}[]{cc}\lambda_{1}^{2}a^{2}+\lambda_{2}^{2}b^{2}&ab(\lambda_{2}-\lambda_{1})\\\
ab(\lambda_{2}-\lambda_{1})&\lambda_{1}^{2}b^{2}+\lambda_{2}^{2}a^{2}\\\
\end{array}\right)$, where $a=\frac{1}{\sqrt{2}},~{}b=\frac{1}{\sqrt{2}}$, and
$0<\lambda_{1}<1,\lambda_{2}=1$. Consider the exact smooth solution

(106) $\displaystyle U({\bf x},t)$ $\displaystyle=$
$\displaystyle\text{sin}(\pi x_{1})\text{sin}(\pi x_{2})\text{sin}(\sqrt{3}\pi
t),$ (107) $\displaystyle v$ $\displaystyle=$ $\displaystyle\frac{\partial
U}{\partial t},~{}{\bm{\sigma}}=-A^{\frac{1}{2}}\nabla U,$

to the IBVP (1), with nonhomogeneous soure $f$ and nonhomogeneous Neumann
boundary conditions.

In the Tables 8.2.2-8.2.2, the collum labelled ``Error" shows the numerical
errors in relative $L^{2}(\Omega\times\\{T\\})-$norm and
$|\cdot|_{\text{DG}}-$norm: if $p=\in\\{1,2\\}$, it is given for the mesh
level $l=4$, else if $p\in\\{3,4\\}$ for $l=3$. The convergence rates with
respect to $h$ are given in the Table 8.2.2 for different choices of $p=q+1$.

Table 6:

Convergence rates of the combined space-time DG scheme w.r.t. $h$.  $v_{h}$
$\bm{\sigma}_{h}$ $|\cdot|_{\text{DG}}$ $(p,q)$ Error Rate Error Rate Error
Rate $(2,1)$ 1.38e-2 3.32 1.23e-2 3.13 1.56e-1 1.87 $(3,2)$ 1.43e-3 3.95
2.11e-3 3.42 2.22e-2 2.53 $(4,3)$ 5.55e-5 4.80 6.24e-5 5.07 8.49e-4 3.67

From the Table 8.2.2, we observe the convergence orders of the errors in
$L^{2}(\Omega\times\\{T\\})-$norm between $\mathcal{O}(h^{p+\frac{1}{2}})$ and
$\mathcal{O}(h^{p+1})$. Besides, the convergence rates of errors in
$|\cdot|_{\text{DG}}-$norm are
$O(h^{\text{min}\\{p+\frac{1}{2},q+\frac{1}{2}\\}})$, which support
convergence rate optimality on the uniform refined meshes of the combined
numerical DG scheme just as for the homogeneous case.

The convergence rates with respect to $\rho$ are given in the Table 8.2.2.

Table 7:

Convergence rates of the TDG schemes w.r.t. $\rho$.  $v_{h}$ $\bm{\sigma}_{h}$
$|\cdot|_{\text{DG}}$ ($p,q$) $\rho$ Error Rate Error Rate Error Rate $(1,0)$
$16$ 3.59e-2 7.42e-2 6.58e-1 32 5.55e-2 0.6285 8.32e-2 0.1652 5.68e-1 -0.2122
64 5.51e-2 -0.0104 9.87e-2 0.2465 5.30e-1 -0.0999 $(2,1)$ $16$ 2.63e-2 3.40e-2
3.86e-1 32 2.93e-2 0.1558 4.01e-2 0.2381 4.01e-1 0.0550 64 3.27e-2 0.1584
4.68e-2 0.2229 4.24e-1 0.0805 $(3,2)$ $2$ 3.25e-3 8.01e-3 9.81e-2 4 3.77e-3
0.2141 8.97e-3 0.1633 1.04e-1 0.0843 8 3.84e-3 0.0265 9.26e-3 0.0459 1.05e-1
0.0138

We can see that $||v-v_{h}||_{L^{2}(\Omega\times
T)},||\bm{\sigma}-\bm{\sigma}_{h}||_{L^{2}(\Omega\times
T)^{2}},|(v,{\bm{\sigma}})-(v_{h},{\bm{\sigma}}_{h})|_{\text{DG}^{+}}\lesssim
C(v,\bm{\sigma})\mathcal{O}(\rho^{\frac{1}{4}})$, which coincides with the
optimal convergence rates of the errors with respect to $\rho$ indicated by
the Theorem 5 for the homogeneous case.

####  8.2.3 3D space case for smooth solution with Neumann boundary conditions

We choose the space-time domain $Q=(0,1)^{2}\times(0,1)$, and set the
anisotropic matrix
$A=\left(\begin{array}[]{ccc}\lambda_{1}^{2}a^{2}+\lambda_{2}^{2}b^{2}&ab(\lambda_{2}-\lambda_{1})&0\\\
ab(\lambda_{2}-\lambda_{1})&\lambda_{1}^{2}b^{2}+\lambda_{2}^{2}a^{2}&0\\\
0&0&1\\\ \end{array}\right)$, where
$a=\frac{1}{\sqrt{2}},~{}b=\frac{1}{\sqrt{2}}$, and
$0<\lambda_{1},\lambda_{2}<1$. Choose the space-time domain
$Q=(0,1)^{3}\times(0,1)$, and consider the exact smooth solution

(108) $\displaystyle U({\bf x},t)$ $\displaystyle=$
$\displaystyle\text{sin}(\pi x_{1})\text{sin}(\pi x_{2})\text{sin}(\pi
x_{3})\text{sin}(2\pi t),$ (109) $\displaystyle v$ $\displaystyle=$
$\displaystyle\frac{\partial U}{\partial
t},~{}{\bm{\sigma}}=-A^{\frac{1}{2}}\nabla U,$

to the IBVP (1), with nonhomogeneous soure $f$ and nonhomogeneous Neumann
boundary conditions.

The convergence rates with respect to $h$ are given in the Table 8.2.3 for
different choices of $p=q+1$. The collum labelled ``Error" shows the numerical
errors in relative $L^{2}(\Omega\times\\{T\\})-$norm and
$|\cdot|_{\text{DG}}-$norm: if $p=1$, it is given for the mesh level $l=4$,
else if $p\in\\{2,3\\}$ for $l=3$.

Table 8:

Convergence rates of the combined space-time DG scheme w.r.t. $h$.  $v_{h}$
$\bm{\sigma}_{h}$ $|\cdot|_{\text{DG}}$ ($p,q$) Error Rate Error Rate Error
Rate $(1,0)$ 6.96e-2 1.85 1.20e-1 1.21 2.07e-1 1.51 $(2,1)$ 2.10e-2 2.71
5.23e-2 2.07 2.22e-1 1.47 $(3,2)$ 2.97e-3 3.53 6.26e-3 3.17 2.87e-2 2.55

In the Table 8.2.3, we can see that, the convergence orders of the errors in
$L^{2}(\Omega\times\\{T\\})-$norm are between $\mathcal{O}(h^{p})$ and
$\mathcal{O}(h^{p+1})$. Besides, the convergence rates of errors in
$|\cdot|_{\text{DG}}-$norm are
$O(h^{\text{min}\\{p+\frac{1}{2},q+\frac{1}{2}\\}})$, which support
convergence rate optimality on the uniform refined meshes of the combined
numerical DG scheme just as for the homogeneous case.

The convergence rates with respect to $\rho$ are given in the Table 8.2.3.

Table 9:

Convergence rates of the TDG schemes w.r.t. $\rho$.  $v_{h}$ $\bm{\sigma}_{h}$
$|\cdot|_{\text{DG}}$ ($p,q$) $\rho$ Error Rate Error Rate Error Rate $(1,0)$
$4$ 7.47e-2 7.41e-2 1.13e-1 8 7.25e-2 -0.0431 7.86e-2 0.0851 1.09e-1 -0.0520
16 7.07e-2 -0.0362 8.11e-2 0.0452 1.06e-1 -0.0403 $(2,1)$ $4$ 2.39e-2 5.69e-2
2.54e-1 8 2.66e-2 0.1544 6.56e-2 0.2053 2.72e-1 0.0988 16 3.01e-2 0.1783
7.73e-2 0.2368 2.98e-1 0.1317 $(3,2)$ $4$ 3.03e-3 7.24e-3 3.01e-2 8 3.49e-2
0.2039 8.25e-3 0.1884 3.34e-2 0.1501 16 4.05e-2 0.2147 9.42e-3 0.1913 3.69e-2
0.1438

We can see that $||v-v_{h}||_{L^{2}(\Omega\times
T)},||\bm{\sigma}-\bm{\sigma}_{h}||_{L^{2}(\Omega\times
T)^{2}},|(v,{\bm{\sigma}})-(v_{h},{\bm{\sigma}}_{h})|_{\text{DG}^{+}}\lesssim
C(v,\bm{\sigma})\mathcal{O}(\rho^{\frac{1}{4}})$, which coincides with the
optimal convergence rates of the errors with respect to $\rho$ indicated by
the Theorem 5 for the homogeneous case.

###  8.3 Comparisons of TDG and DG methods for 2D space case with
nonhomogeneous sources

We choose the space-time domain $Q=(0,1)^{2}\times(0,1)$. Consider the exact
smooth solution

(110) $\displaystyle u({\bf x},t)$ $\displaystyle=$
$\displaystyle\text{sin}(\pi x_{1})\text{sin}(\pi x_{2})\text{sin}(\sqrt{3}\pi
t),$ (111) $\displaystyle v$ $\displaystyle=$ $\displaystyle\frac{\partial
U}{\partial t},~{}\bm{\sigma}=-A^{\frac{1}{2}}\nabla U,$

to the IBVP (1), with nonhomogeneous soure $f$. We use uniform mesh with
$h_{{\bf x}}\approx h_{t}=2^{-l},l\in\mathbb{N}$.

####  8.3.1 The case of Dirichlet boundary conditions

We would like to compare the errors of the approximations generated by the
proposed combined DG method of Section 7 and the high-order DG method of [1].
We choose $\Gamma_{D}=\partial\Omega$. The convergence rates are given in the
Table 8.3.1 for different choices of $p=q+1$.

Table 10:

Convergence rates of the combined space-time DG scheme w.r.t. $h$.  $v_{h}$
$\bm{\sigma}_{h}$ $|\cdot|_{\text{DG}}$ $(p,q)$ Method DOFs Error Rate Error
Rate Error Rate $(2,1)$ TDG 7680 1.30e-2 2.79 6.33e-3 3.61 1.57e-1 1.86 DG
12288 2.18e-2 2.89 4.16e-2 2.66 3.88e-1 1.48 $(3,2)$ TDG 5184 2.96e-3 3.48
3.08e-3 3.93 3.23e-2 2.77 DG 5184 8.50e-3 3.41 6.65e-3 4.00 1.46e-1 2.52
$(4,3)$ TDG 7560 2.06e-4 4.64 2.86e-4 4.99 2.45e-3 3.70 DG 12288 6.05e-4 4.01
8.28e-4 4.14 1.42e-2 3.53

We can see from Table 8.3.1 that, the convergence orders of errors generated
by the global TDG method in $L^{2}(\Omega\times\\{T\\})-$norm are between
$\mathcal{O}(h^{p+\frac{1}{2}})$ and $\mathcal{O}(h^{p+1})$. Besides, the
convergence rates of errors in $|\cdot|_{\text{DG}}-$norm are
$O(h^{\text{min}\\{p+\frac{1}{2},q+\frac{1}{2}\\}})$, which support
convergence rate optimality on the uniform refined meshes of the global TDG
method just as for the homogeneous case. Moreover, the approximations
generated by the TDG are more accurate than those generated by the high-order
DG method, even if the DOFs of the DG method are significantly higher than
that of the TDG method.

####  8.3.2 The case of mixed boundary conditions

We choose $\Gamma_{D}=\\{x_{1}=0,1\\}\times[0,1]$ and
$\Gamma_{N}=[0,1]\times\\{x_{2}=0,1\\}$. The convergence rates are given in
the Table 8.3.2 for different choices of $p=q+1$.

Table 11:

Convergence rates of the combined space-time DG scheme w.r.t. $h$.  $v_{h}$
$\bm{\sigma}_{h}$ $|\cdot|_{\text{DG}}$ $(p,q)$ Method DOFs Error Rate Error
Rate Error Rate $(2,1)$ TDG 7680 1.47e-2 3.45 8.04e-3 2.98 1.60e-1 1.86 DG – –
– – – – – $(3,2)$ TDG 5184 2.81e-3 3.55 2.21e-3 3.74 2.13e-2 2.65 DG – – – – –
– – $(4,3)$ TDG 7560 6.66e-4 4.53 4.94e-4 4.92 8.05e-3 3.56 DG – – – – – – –

It can be seen from Table 8.3.2 that, the TDG scheme still works, but high-
order DG scheme fails (here ``–" represents that the corresponding numerical
method has no accuracy) for the case of nonhomogeneous Neumann boundary
conditions, which verifies the existing theoretical and numerical conclusions,
see Remarks 6.3, 7.1 and [1, Sections 3 and 6].

####  8.3.3 The case of Neumann boundary conditions

We choose $\Gamma_{N}=\partial\Omega$. The convergence rates are given in the
Table 8.3.3 for different choices of $p=q+1$.

Table 12:

Convergence rates of the combined space-time DG scheme w.r.t. $h$.  $v_{h}$
$\bm{\sigma}_{h}$ $|\cdot|_{\text{DG}}$ $(p,q)$ Method DOFs Error Rate Error
Rate Error Rate $(2,1)$ TDG 7680 1.38e-2 3.32 1.23e-2 3.13 1.56e-1 1.87 DG – –
– – – – – $(3,2)$ TDG 5184 1.43e-3 3.95 2.11e-3 3.42 2.22e-2 2.53 DG – – – – –
– – $(4,3)$ TDG 7560 5.55e-5 4.80 6.24e-5 5.07 8.49e-4 3.67 DG – – – – – – –

The conclusion coincides with the above Section.

###  8.4 Discussion on the size of local nonhomogeneous problems

Consider the exact solution of the Section 8.2.2. The convergence rates are
given in the Table 8.4 for different space size of local nonhomogeneous
problems and different choices of $p=q+1$.

Table 13:

Convergence rates of the combined space-time DG scheme w.r.t. $h$.  $v_{h}$
$\bm{\sigma}_{h}$ $|\cdot|_{\text{DG}}$ ($p,q$)
$(\text{size}_{x_{1}},\text{size}_{x_{2}})$ Error Rate Error Rate Error Rate
$(2,1)$ (1 1) 1.38e-2 3.32 1.23e-2 3.13 1.56e-1 1.87 $(h^{-1}_{{\bf
x}},h^{-1}_{{\bf x}})$ 1.45e-2 3.29 1.09e-2 3.28 1.67e-1 1.95 $(3,2)$ (1 1)
1.43e-3 3.95 2.11e-3 3.42 2.22e-2 2.53 $(h^{-1}_{{\bf x}},h^{-1}_{{\bf x}})$
3.12e-3 3.23 3.00e-3 3.15 3.45e-2 2.53 $(4,3)$ (1 1) 5.55e-5 4.80 6.24e-5 5.07
8.49e-4 3.67 $(h^{-1}_{{\bf x}},h^{-1}_{{\bf x}})$ 4.77e-5 5.02 5.97e-5 5.13
7.43e-4 3.86

Here $(\text{size}_{x_{1}},\text{size}_{x_{2}})=(1,1)$ represents the involved
number of spacial elements employed by the overlapping local DG variational
formulation (92) in each coordinate direction. For the case of
$(\text{size}_{x_{1}},\text{size}_{x_{2}})=(h^{-1}_{{\bf x}},h^{-1}_{{\bf
x}})$, it represents that the nonoverlapping variant of local DG variational
formulation (92) are defined on each nonoverlapping time slab $D_{n}(1\leq
n\leq N)$. It can be seen that, if nonhomogeneous local problems are defined
on $D_{n}$, the combined DG scheme can also preserve the same orders of
convergence as for the case of
$(\text{size}_{x_{1}},\text{size}_{x_{2}})=(1,1)$, i.e. the local smooth
overlapping case (92).

##  9 Future Work

The presently proposed space-time Trefftz DG discretizations will be further
investigated for the case of piecewise anisotropic media; namely, $A$ is
piecewise positive definite matrix. Besides, the Trefftz numerical methods for
the anisotropic wave equations with the Robin boundary conditions will also be
studied.

##  10 Appendix: the derivation of Eq. (24)

For convenience, we use ${\bf p}_{1},{\bf p}_{2},\cdots,{\bf p}_{d}$ to denote
the column vectors of $P$, and use ${\bf q}_{1}^{T},{\bf
q}_{2}^{T},\cdots,{\bf q}_{d}^{T}$ to denote the row vectors of $P$. Then each
of these vectors is a unit vector, and ${\bf p}_{1},{\bf p}_{2},\cdots,{\bf
p}_{d}$ (and ${\bf q}_{1},{\bf q}_{2},\cdots,{\bf q}_{d}$) are orthogonal each
other.

For $1\leq i\leq d$, by the coordinate transformation (21), we have

$\frac{\partial\hat{v}}{\partial
x_{i}}=\sum_{k=1}^{d}\frac{\partial\hat{v}}{\partial\hat{x}_{k}}\frac{\partial\hat{x}_{k}}{\partial
x_{i}}=\sum_{k=1}^{d}\frac{\partial\hat{v}}{\partial\hat{x}_{k}}\lambda_{k}^{-\frac{1}{2}}q_{ik},$

which yields

$\nabla\hat{v}=P^{T}\Lambda^{-\frac{1}{2}}\hat{\nabla}\hat{v}.$

Combining it with (23), yields

(112) $A^{\frac{1}{2}}\nabla
v=A^{\frac{1}{2}}\nabla\hat{v}=A^{\frac{1}{2}}P^{T}\Lambda^{-\frac{1}{2}}\hat{\nabla}\hat{v}=P^{T}\hat{\nabla}\hat{v}.$

Next, by (23), it holds that

$A^{\frac{1}{2}}{\bm{\sigma}}=P^{T}\Lambda^{\frac{1}{2}}\hat{\bm{\sigma}}=\sum_{i=1}^{d}\lambda_{i}^{\frac{1}{2}}{\bf
q}_{i}\hat{\sigma}_{i}.$

By the chain rule, we obtain

(113) $\displaystyle\nabla\cdot\bigg{(}A^{\frac{1}{2}}{\bm{\sigma}}\bigg{)}$
$\displaystyle=$
$\displaystyle\nabla\cdot\bigg{(}\sum_{i=1}^{d}\lambda_{i}^{\frac{1}{2}}{\bf
q}_{i}\hat{\sigma}_{i}\bigg{)}=\sum_{j=1}^{d}\sum_{i=1}^{d}\lambda_{i}^{\frac{1}{2}}q_{ji}\frac{\partial\hat{\sigma}_{i}}{\partial
x_{j}}$ (114) $\displaystyle=$
$\displaystyle\sum_{j=1}^{d}\sum_{i=1}^{d}\lambda_{i}^{\frac{1}{2}}q_{ji}\bigg{(}\sum_{k=1}^{d}\frac{\partial\hat{\sigma}_{i}}{\partial\hat{x}_{k}}\frac{{\partial\hat{x}_{k}}}{{\partial
x_{j}}}\bigg{)}=\sum_{j=1}^{d}\sum_{i=1}^{d}\lambda_{i}^{\frac{1}{2}}q_{ji}\bigg{(}\sum_{k=1}^{d}\frac{\partial\hat{\sigma}_{i}}{\partial\hat{x}_{k}}\lambda_{k}^{-\frac{1}{2}}q_{jk}\bigg{)}$
(115) $\displaystyle=$
$\displaystyle\sum_{i=1}^{d}\sum_{k=1}^{d}\lambda_{i}^{\frac{1}{2}}\lambda_{k}^{-\frac{1}{2}}\frac{\partial\hat{\sigma}_{i}}{\partial\hat{x}_{k}}\bigg{(}\sum_{j=1}^{d}q_{ji}q_{jk}\bigg{)}=\sum_{i=k=1}^{d}\frac{\partial\hat{\sigma}_{i}}{\partial\hat{x}_{k}}=\hat{\nabla}\cdot\hat{\bm{\sigma}}.$

## References

  * [1] P. Bansal, A. Moiola, I. Perugia and C. Schwab, Space-time discontinuous Galerkin approximation of acoustic waves with point singularities, IMA J. Numer. Anal., 41 (2021), pp. 2056-2109.
  * [2] S. Brenner and L. Scott, Mathematical theory of finite element methods, 3rd ed., Texts Appl. Math., Springer-Verlag, New York, 2007.
  * [3] A. Cangiani, Z. Dong and E. Georgoulis, $hp-$Version space-time discontinuous Galerkin methods for parabolic problems on prismatic meshes, SIAM J. Sci. Comput., 39(2017), pp. A1251-A1279.
  * [4] J. Carcione, Wave Fields in Real Media: Wave Propagation in Anisotropic, Anelastic, Porous and Electromagnetic Media, 2nd edition, Elsevier, Ams-terdam, Netherlands, 2007.
  * [5] O. Cessenat, B. Despres, Application of an ultra weak variational formulation of elliptic PDEs to the two-dimensional Helmholtz problem. SIAM J. Numer. Anal. 35(1998), pp. 255-299.
  * [6] H. Egger, F. Kretzschmar, S. Schnepp, I. Tsukerman and T. Weiland, Transparent boundary conditions for a discontinuous Galerkin Trefftz method, Appl. Math. Comput., 267(2015), pp. 42-55.
  * [7] H. Egger, F. Kretzschmar, S. Schnepp and T. Weiland, A space-time discontinuous Galerkin Trefftz method for time dependent Maxwell's equations, SIAM J. Sci. Comput., 37(2015), pp. B689-B711.
  * [8] R. Falk and G. Richter, Explicit finite element methods for symmetric hyperbolic equations, SIAM J. Numer. Anal., 36 (1999), pp. 935-952.
  * [9] L. Gao and L. Huang, An improved rotated staggered-grid finite-difference method with fourth-order temporal accuracy for elastic-wave modeling in anisotropic media, J. Comput. Phys., 350(2017), pp. 361-386.
  * [10] J. Gopalakrishnan, M. Hochsteger, J. Schöberl, C. Wintersteiger, An explicit mapped tent pitching scheme for Maxwell equations, in Spectral and High Order Methods for Partial Differential Equations: ICOSAHOM 2018, S. J. Sherwin, D. Moxey, J. Peiró, P. E. Vincent, and C. Schwab, eds., vol. 134 of Lecture Notes in Computational Science and Engineering, 2020, pp. 359-369.
  * [11] J. Gopalakrishnan, J. Schöberl, and C. Wintersteiger, Mapped tent pitching schemes for hyperbolic systems, SIAM J. Sci. Comput., 39 (2017), pp. B1043-B1063.
  * [12] R. Hiptmair, A. Moiola, and I. Perugia, Plane wave discontinuous Galerkin methods for the 2D Helmholtz equation: analysis of the $p$-version. SIAM J. Numer. Anal., 49(2011), pp. 264-284.
  * [13] R. Hiptmair, A. Moiola, and I. Perugia, Error analysis of Trefftz-discontinuous Galerkin methods for the time-harmonic Maxwell equations, Math. Comp., 82(2013), pp. 247-268.
  * [14] R. Hiptmair, A. Moiola, and I. Perugia, A Survey of Trefftz Methods for the Helmholtz Equation, Lect. Notes Comput. Sci. Eng., Springer, 2016. pp. 237-278.
  * [15] C. Howartha, P. Childs b and A. Moiola, Implementation of an interior point source in the ultra weak variational formulation through source extraction, J. Comput. Appl. Math., 271(2014), 295-306.
  * [16] Q. Hu and L. Yuan, A Plane Wave Least-Squares Method for Time-Harmonic Maxwell's Equations in Absorbing Media, SIAM J. Sci. Comput., 36(2014), pp. A1911-A1936.
  * [17] Q. Hu and L. Yuan, A Plane wave method combined with local spectral elements for nonhomogeneous Helmholtz equation and time-harmonic Maxwell equations, Adv. Comput. Math., 44(2018), pp. 245-275.
  * [18] T. Hughes and G. Hulbert, Space-time finite element methods for elastodynamics: formulations and error estimates, Comput. Methods Appl. Mech. Engrg., 66(1988), pp. 339-363.
  * [19] G. Hulbert and T. Hughes, Space-time finite element methods for second-order hyperbolic equations, Comput. Methods Appl. Mech. Engrg., 84 (1990), pp. 327-348.
  * [20] T. Huttunen, M. Malinen, P. Monk, Solving Maxwell's equations using the ultra weak variational formulation, J. Comput. Phys., 223 (2007), pp. 731-758.
  * [21] S. Kapita, P. Monk and T. Warburton, Residual-based adaptivity and PWDG methods for the helmholtz equation, SIAM J. Sci. Comput., 37(2015), pp. A1525-A1553.
  * [22] F. Kretzschmar, The discontinuous Galerkin Trefftz method, Ph.D. thesis, Technische Universität Darmstadt, http://tuprints.ulb.tu-darmstadt.de/5166/.
  * [23] T. Alkhalifah, An acoustic wave equation for anisotropic media, Geophysics, 65(2000), pp. 1239-1250.
  * [24] F. Kretzschmar, A. Moiola, I. Perugia and S. Schnepp, A priori error analysis of space-time Trefftz discontinuous Galerkin methods for wave problems, IMA J. Numer. Anal., 36(2016), pp. 1599-1635.
  * [25] F. Kretzschmar, S. Schnepp, I. Tsukerman and T. Weiland, Discontinuous Galerkin methods with Trefftz approximations, J. Comput. Appl. Math., 270(2014), pp. 211-222.
  * [26] J. Lions, E. Magenes, Non-homogeneous boundary value problems and applications, vol. I. Springer, New York: (Translated from the French by P, p. 181. Kenneth, Die Grundlehren der mathematischen Wissenschaften, Band (1972)).
  * [27] A. Moiola and I. Perugia, A space-time Trefftz discontinuous Galerkin method for the acoustic wave equation in first-order formulation, Numer. Math., 138(2018), pp. 389-435.
  * [28] P. Monk and G. Richter, A Discontinuous Galerkin Method for Linear Symmetric Hyperbolic Systems in Inhomogeneous Media, J. Sci. Comput., 22(2005), pp. 443-477.
  * [29] F. Müller, Numerical analysis of finite element methods for second order wave equations in polygons, PhD thesis, ETH Zurich, 2017. Available at https://www.research-collection.ethz.ch/handle/20.500.11850/167502.
  * [30] F. Müller, D. Schötzau, and C. Schwab, Discontinuous Galerkin methods for acoustic wave propagation in polygons, J. Sci. Comput., 77(2018), pp. 1909-1935.
  * [31] I. Perugia, J. Schöberl, P. Stocker and C. Wintersteiger, Tent pitching and Trefftz-DG method for the acoustic wave equation, Comput. Math. Appl., 79 (2020), pp. 2987-3000.
  * [32] Y. Qu, J. Huang, Z. Li, Z. Guan and J. Li, Attenuation compensation in anisotropic least-squares reverse time migration, Geophysics, 82(2017), pp. S411-S423.
  * [33] C. Sayers, Seismic characterization of reservoirs containing multiple fracture sets, Geophys. Prospect., 57(2009), pp. 187-192.
  * [34] O. Steinbach and M. Zank, A Stabilized Space Time Finite Element Method for the Wave Equation, Springer International Publishing, Cham, 2019, pp. 341-370.
  * [35] I. Tsvankin, Seismic Signatures and Analysis of Reflection Data in Anisotropic Media, 2nd edition, Elsevier Science Publ. Co., Inc., 2005.
  * [36] L. Yuan and Q. Hu, Plane wave discontinuous Galerkin methods for the Helmholtz equation and Maxwell equations in Anisotropic Media, Comput. Math. Appl., 97( 2021), pp. 355-374.
  * [37] Q. Zhan, M. Zhuang, Y. Mao, Q. Liu, Unified Riemann solution for multi-physics coupling: Anisotropicporoelastic/elastic/fluid interfaces, J. Comput. Phys., 402(2020), ID: 108961.