Cauchy–Schwarz inequality, we infer

$\displaystyle\left|\mathfrak{T}_{3}\right|=\sum_{e\in\mathcal{E}_{h}^{I}}\int_{e}{\boldsymbol{\beta}}\\{\\!\\!\\{v\\}\\!\\!\\}\cdot[\\![w_{h}]\\!]$
$\displaystyle\leq\left(2\sum_{e\in{\mathcal{E}}_{h}^{I}}\int_{e}\left|{\boldsymbol{\beta}}\cdot{{\mathbf{n}}}_{e}\right|\\{\\!\\!\\{v\\}\\!\\!\\}^{2}\right)^{\frac{1}{2}}\
\left(\frac{1}{2}\sum_{e\in{\mathcal{E}}_{h}^{I}}\int_{e}\left|{\boldsymbol{\beta}}\cdot{{\mathbf{n}}}_{e}\right|[\\![w_{h}]\\!]^{2}\right)^{\frac{1}{2}}$
$\displaystyle\leq\left(\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(\Omega)}\sum_{T\in{\mathcal{T}}_{h}}\left\|v\right\|^{2}_{L^{2}(\partial
T)}\right)^{\frac{1}{2}}\
{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|w_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{ar}$
$\displaystyle\leq\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(\Omega)}^{\frac{1}{2}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{ar,*}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|w_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{ar},$

where in the second step we use the formula
$2\\{\\!\\!\\{v\\}\\!\\!\\}^{2}=\frac{1}{2}(v_{1}+v_{2})^{2}\leq
v_{1}^{2}+v_{2}^{2}$ that follows from the Young’s inequality. Collecting all
the above bounds yields the thesis with
$M_{ar}=\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(\Omega)}k_{min}^{-\frac{1}{2}}+\left\|\sigma\right\|_{L^{\infty}(\Omega)}+2+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(\Omega)}^{\frac{1}{2}}$.
∎

The following lemma establishes the boundedness of the diffusion-advection-
reaction bilinear form $a^{dar}_{h}$ on $V_{*h}\times V_{h}$ using the norms
$\left\|\cdot\right\|_{dar,*}$ and $\left\|\cdot\right\|_{dar}$.

######  Lemma 4.3.3 (Boundedness).

There exists $M>0$, independent of $h$, such that

$a_{h}^{\text{dar}}(v,w_{h})\leq
M{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{dar,*}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|w_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{dar}\quad\forall(v,w_{h})\in
V_{*h}\times V_{h},$

where $M=M_{d}+M_{ar}$, with $M_{d}$ as in Lemma 4.3.1 and $M_{ar}$ as in
Lemma 4.3.2.

###### Proof.

Let $(v,w_{h})\in V_{*h}\times V_{h}$, Lemma 4.3.1 states

$a_{h}^{d}(v,w_{h})\leq
M_{d}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{d,*}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|w_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{d},$

and Lemma 4.3.2 yields

$a_{h}^{ar}(v,w_{h})\leq
M_{ar}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{ar,*}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|w_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{dar}.$

For the norm definitions, from these two bounds we derive the assertion with
$M=M_{d}+M_{ar}$. ∎

###  4.4 Error analysis

The aim of this section is to infer a convergence rate in $h$ for the Galerkin
error $u-u_{h}$ measured in the $\left\|\cdot\right\|_{dar}$-norm. Since
consistency, discrete coercivity and boundedness are satisfied, we can apply
Theorem 1.1.3 to the DG method (3.8) for any discrete space
$V_{h}\subseteq\mathbb{P}^{p}_{d}({\mathcal{T}}_{h})$.

######  Theorem 4.4.1 (Error estimate).

Let $u\in V_{*h}$ solve (1.15) and let $u_{h}$ solve (3.8) with the penalty
parameter $\gamma$ as in Lemma 4.2.1. Then,

${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-u_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{dar}\leq\left(1+\frac{M}{\alpha}\right)\inf_{v_{h}\in
V_{h}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{dar,*},$
(4.6)

with $\alpha$ as in Lemma 4.2.3 and $M$ as in Lemma 4.3.3.

In order to deduce the $h$-convergence, given the quasi-optimality inequality
(4.6), we need to estimate the right-hand side $\inf_{v_{h}\in
V_{h}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{dar,*}$.
The estimation of this quantity depends on the approximation properties of the
discrete space $V_{h}$.

If we choose as discrete space $V_{h}$ in the DG scheme (3.8) the broken space
of degree-$p$ polynomials $\mathbb{P}^{p}_{d}({\mathcal{T}}_{h})$ we obtain
the standard DG method. The Bramble-Hilbert Lemma states the approximation
properties that are achieved using the full polynomial space [5, Lemma 4.3.8].

######  Lemma 4.4.2 (Bramble-Hilbert).

Let ${\mathcal{T}}_{\mathcal{H}}$ be a mesh sequence with the star-shaped
property 2.1.12. Then, for all $h\in\mathcal{H}$, all $T\in{\mathcal{T}}_{h}$
and for all polynomial degrees $p$, there exist a polynomial
$v\in{P}^{p}_{d}(T)$ and $C_{opt}>0$, independent of $T$ and $h$, such that,
for all $q\in\\{0,\dots,p+1\\}$ and all $u\in H^{q}(T)$, there holds

$\left|u-v\right|_{H^{m}(T)}\leq
C_{opt}h_{T}^{q-m}\left|v\right|_{H^{q}(T)}\quad\forall m\in\\{0,\dots,q\\}.$
(4.7)

The next theorem establishes the convergence rate in $h$ for the error
$u-u_{h}$ in the $\left\|\cdot\right\|_{dar}$-norm. From the quasi-optimality
inequality (4.6) we deduce the convergence rate for the DG method, using the
the error bound (4.7). We recall that $V_{*}:=H^{1}_{\Gamma_{D}}(\Omega)\cap
H^{2}({\mathcal{T}}_{h})$.

######  Theorem 4.4.3 (Convergence rate).

Let $p\in\mathbb{N}$ and let $u\in V_{*}\cap H^{p+1}({\mathcal{T}}_{h})$ solve
(1.6)-(1.8). Let $u_{h}$ solve (3.8) with discrete space
$V_{h}=\mathbb{P}_{d}^{p}({\mathcal{T}}_{h})$ and the penalty parameter
$\gamma$ as in Lemma 4.2.1. Then, the following convergence rate holds

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-u_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{dar}\leq$
$\displaystyle\left(1+\frac{M}{\alpha}\right)C_{opt}\Biggl{(}\sum_{T\in{\mathcal{T}}_{h}}\left[\left(1+C_{tr}\right)\left\|{\bm{K}}\right\|_{L^{\infty}(T)}+2C_{g}\gamma
C_{tr}\right.$ (4.8)
$\displaystyle\left.+C_{tr}\left(\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}+1\right)h_{T}+h_{T}^{2}\right]h_{T}^{2p}\left|u\right|^{2}_{H^{p+1}(T)}\Biggr{)}^{\frac{1}{2}}.$

with $\alpha$ as in Lemma 4.2.3 and $M$ as in Lemma 4.3.3.

###### Proof.

Since we have the quasi-optimality inequality (4.6), we now focus on
estimating the quantity $\inf_{v_{h}\in
V_{h}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{dar,*}$.
We recall the definition of the $\left\|\cdot\right\|_{dar,*}$-norm, for $v\in
V_{*h}$,

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}=$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\int_{T}{\bm{K}}\nabla v\cdot\nabla
v+\sum_{e\in{\mathcal{E}}_{h}^{I}\cup{\mathcal{E}}_{h}^{D}}\frac{\gamma}{h_{e}}\int_{e}[\\![v]\\!]^{2}+\frac{1}{2}\sum_{e\in{\mathcal{E}}_{h}}\int_{e}\left|{\boldsymbol{\beta}}\cdot{{\mathbf{n}}}_{e}\right|[\\![v]\\!]^{2}$
$\displaystyle+\left\|v\right\|^{2}_{L^{2}(\Omega)}+\sum_{T\in{\mathcal{T}}_{h}}h_{T}\left\|{\bm{K}}^{\frac{1}{2}}\nabla
v\cdot{{\mathbf{n}}}_{T}\right\|^{2}_{L^{2}(\partial
T)}+\sum_{T\in{\mathcal{T}}_{h}}\left\|v\right\|^{2}_{L^{2}(\partial T)}.$

We recall the assumption that the mesh is graded (2.2) and observe that, from
the Young’s inequality, $[\\![v]\\!]^{2}=(v_{1}-v_{2})^{2}\leq
2(v_{1}^{2}+v_{2}^{2})$ on an internal facet $e=\partial T_{1}\cap\partial
T_{2}$ and $[\\![v]\\!]^{2}=v^{2}$ on a boundary facet $e$. Using these two
facts in the second and third terms of the norm, we rearrange the sums over
mesh elements as sums over elements and obtain the bound:

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}\leq$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\biggl{(}\left\|{\bm{K}}\right\|_{L^{\infty}(T)}\left\|\nabla
v\right\|^{2}_{L^{2}(T)}+2C_{g}\frac{\gamma}{h_{T}}\left\|v\right\|^{2}_{L^{2}(\partial
T)}+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}\left\|v\right\|^{2}_{L^{2}(\partial
T)}$
$\displaystyle+\left\|v\right\|^{2}_{L^{2}(T)}+\left\|{\bm{K}}\right\|_{L^{\infty}(T)}h_{T}\left\|\nabla
v\right\|^{2}_{L^{2}(\partial T)}+\left\|v\right\|^{2}_{L^{2}(\partial
T)}\biggr{)}.$

Owing to the discrete trace inequality (2.8), we get

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}\leq$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\biggl{(}\left\|{\bm{K}}\right\|_{L^{\infty}(T)}\left\|\nabla
v\right\|^{2}_{L^{2}(T)}+2C_{g}\frac{\gamma}{h_{T}}C_{tr}h_{T}^{-1}\left\|v\right\|^{2}_{L^{2}(T)}$
$\displaystyle+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}C_{tr}h_{T}^{-1}\left\|v\right\|^{2}_{L^{2}(T)}+\left\|v\right\|^{2}_{L^{2}(T)}$
$\displaystyle+\left\|{\bm{K}}\right\|_{L^{\infty}(T)}h_{T}C_{tr}h_{T}^{-1}\left\|\nabla
v\right\|^{2}_{L^{2}(T)}+C_{tr}h_{T}^{-1}\left\|v\right\|^{2}_{L^{2}(T)}\biggr{)}.$

Considering the quantity of our interest we have

$\displaystyle\inf_{v_{h}\in
V_{h}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}\leq$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\inf_{v_{h}\in
V_{h}}\biggl{[}\left(1+C_{tr}\right)\left\|{\bm{K}}\right\|_{L^{\infty}(T)}\left\|\nabla(u-v_{h})\right\|^{2}_{L^{2}(T)}$
$\displaystyle+\left(2C_{g}\gamma
h_{T}^{-2}C_{tr}+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}C_{tr}h_{T}^{-1}+1+C_{tr}h_{T}^{-1}\right)\left\|u-v_{h}\right\|^{2}_{L^{2}(T)}\biggr{]}.$

The optimal polynomial approximation properties of the Bramble-Hilbert Lemma
(4.7) yield

$\displaystyle\inf_{v_{h}\in
V_{h}}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}\leq$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\biggl{[}\left(1+C_{tr}\right)\left\|{\bm{K}}\right\|_{L^{\infty}(T)}C_{opt}^{2}h_{T}^{2p}\left|u\right|^{2}_{H^{p+1}(T)}$
$\displaystyle+\biggl{(}2C_{g}\gamma
h_{T}^{-2}C_{tr}+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}C_{tr}h_{T}^{-1}+1+C_{tr}h_{T}^{-1}\biggr{)}\times$
$\displaystyle
C_{opt}^{2}h_{T}^{2(p+1)}\left|u\right|^{2}_{H^{p+1}(T)}\biggr{]}.$
$\displaystyle=$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\biggl{[}\left(1+C_{tr}\right)\left\|{\bm{K}}\right\|_{L^{\infty}(T)}+2C_{g}\gamma
C_{tr}$
$\displaystyle+C_{tr}\left(\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}+1\right)h_{T}+h_{T}^{2}\biggr{]}C_{opt}^{2}h_{T}^{2p}\left|u\right|^{2}_{H^{p+1}(T)}.$

By combining this bound with the inequality (4.6) we infer the assertion. ∎

##  Chapter 5 Quasi-Trefftz space

In this chapter we define the quasi-Trefftz space, prove its high-order
approximation properties and describe an algorithm for the construction of the
basis functions. First, in Section 5.1 we fix some notation and give the
definitions. In Section 5.2 we prove approximation bounds for the local quasi-
Trefftz space and then in Section 5.3 we prove the $h$-convergence of the
quasi-Trefftz DG method for the diffusion-advection-reaction problem. In
Section 5.4 we define a family of basis functions and we explicitly construct
them in Section 5.5.

###  5.1 Definitions and notation

We call multi-indices the vectors
${\boldsymbol{i}}=(i_{1},\ldots,i_{d})\in\mathbb{N}^{d}$. We define their
length $|{\boldsymbol{i}}|:=i_{1}+\cdots+i_{d}$ and we establish a partial
order ${\boldsymbol{i}}\leq{\boldsymbol{j}}$ if $i_{k}\leq j_{k}$ for all
$k\in\\{1,\dots,d\\}$. The factorial and the binomial coefficient are defined,
respectively, as

${\boldsymbol{i}}!:=i_{1}!\cdots
i_{d}!\quad\text{and}\quad\binom{{\boldsymbol{i}}}{{\boldsymbol{j}}}:=\frac{{\boldsymbol{i}}!}{{\boldsymbol{j}}!({\boldsymbol{i}}-{\boldsymbol{j}})!}=\binom{i_{1}}{j_{1}}\cdots\binom{i_{d}}{j_{d}}.$

We use standard multi-index notation for partial derivatives and monomials:

$D^{\boldsymbol{i}}f:=\partial_{x_{1}}^{i_{1}}\cdots\partial_{x_{d}}^{i_{d}}f,\quad\text{and}\quad{{\mathbf{x}}}^{\boldsymbol{i}}=x_{1}^{i_{1}}\cdots
x_{d}^{i_{d}}.$

We recall the following formula

$\sum_{{\bm{i}}\in\mathbb{N}^{d},\left|{\bm{i}}\right|=k}\frac{1}{{\bm{i}}!}=\frac{d^{k}}{k!}\quad\forall
k\in\mathbb{N},$ (5.1)

obtained from the multinomial theorem

$(x_{1}+\dots+x_{d})^{k}=\sum_{{\bm{i}}\in\mathbb{N}^{d},\left|{\bm{i}}\right|=k}\frac{k!}{{\bm{i}}!}{{\mathbf{x}}}^{{\bm{i}}}\quad\forall{{\mathbf{x}}}\in\mathbb{R}^{d},k\in\mathbb{N},$

by choosing ${{\mathbf{x}}}=(1,\dots,1)$ [27, p. 198]. We use the canonical
basis of $\mathbb{R}^{d}$, namely
$\\{{{\mathbf{e}}}_{k}\in\mathbb{R}^{d},1\leq k\leq
d,({{\mathbf{e}}}_{k})_{l}=\delta_{kl}\\}$. We recall the Leibniz product rule
for multi-indices:

$D^{\boldsymbol{i}}(fg)=\sum_{\begin{subarray}{c}{\boldsymbol{j}}\in\mathbb{N}^{d},\;{\boldsymbol{j}}\leq{\boldsymbol{i}}\end{subarray}}\binom{{\boldsymbol{i}}}{{\boldsymbol{j}}}D^{\boldsymbol{j}}fD^{{\boldsymbol{i}}-{\boldsymbol{j}}}g.$
(5.2)

Let $p\in\mathbb{N}$ , let $\Upsilon\subset\mathbb{R}^{d}$ be an open,
bounded, Lipschitz set and let ${{\mathbf{z}}}\in\Upsilon$.

######  Definition 5.1.1 (Taylor polynomial).

Given $v\in C^{p}(\Upsilon)$, the Taylor polynomial of order $p+1$ of $v$,
centred at ${{\mathbf{z}}}\in\Upsilon$, is defined as

${\mathsf{T}}^{p+1}_{{{\mathbf{z}}}}[v]({{\mathbf{x}}}):=\sum_{|{\boldsymbol{j}}|\leq
p}\frac{1}{{\boldsymbol{j}}!}D^{\boldsymbol{j}}v({{\mathbf{z}}})({{\mathbf{x}}}-{{\mathbf{z}}})^{{\boldsymbol{j}}}.$

Notice that ${\mathsf{T}}^{p+1}_{{{\mathbf{z}}}}[v]$ is a polynomial of degree
at most $p$. We highlight the following property: for every multi-index
${\boldsymbol{i}}$

$\begin{split}D^{\boldsymbol{i}}{\mathsf{T}}^{p+1}_{{{\mathbf{z}}}}[v]({{\mathbf{x}}})&=\sum_{\begin{subarray}{c}\left|{\bm{j}}\right|\leq
p\\\
{\bm{j}}\geq{\bm{i}}\end{subarray}}\frac{1}{{\boldsymbol{j}}!}D^{{\boldsymbol{j}}}v({{\mathbf{z}}})\frac{{\boldsymbol{j}}!}{({\boldsymbol{j}}-{\boldsymbol{i}})!}({{\mathbf{x}}}-{{\mathbf{z}}})^{{\boldsymbol{j}}-{\boldsymbol{i}}}\\\
&=\sum_{\left|\bm{\gamma}\right|\leq
p-\left|{\bm{i}}\right|}\frac{1}{\bm{\gamma}!}D^{\bm{\gamma}+{\bm{i}}}v({{\mathbf{z}}})({{\mathbf{x}}}-{{\mathbf{z}}})^{\bm{\gamma}}={\mathsf{T}}^{p+1-\left|{\bm{i}}\right|}_{{{\mathbf{z}}}}[D^{{\bm{i}}}v]({{\mathbf{x}}}).\end{split}$
(5.3)

From the first step, evaluated at ${{\mathbf{z}}}$, it follows that

$D^{\boldsymbol{i}}{\mathsf{T}}^{p+1}_{{{\mathbf{z}}}}[v]({{\mathbf{z}}})=\left\\{\begin{array}[]{ll}D^{\boldsymbol{i}}v({{\mathbf{z}}})&\text{
if }|{\boldsymbol{i}}|\leq p,\\\ 0&\text{ if
}|{\boldsymbol{i}}|>p.\end{array}\right.$ (5.4)

We also recall the Lagrange’s form of the Taylor remainder [9, Cor. 3.19]: if
$v\in C^{p+1}(\Upsilon)$ and the segment $S$ with extremes ${{\mathbf{z}}}$
and ${{\mathbf{x}}}$ is contained in $\Upsilon$, then exists
${{\mathbf{x}}}_{*}\in S$ such that

$v({{\mathbf{x}}})-{\mathsf{T}}^{p+1}_{{{\mathbf{z}}}}[v]({{\mathbf{x}}})=\sum_{|{\boldsymbol{j}}|=p+1}\frac{1}{{\boldsymbol{j}}!}D^{\boldsymbol{j}}v({{\mathbf{x}}}_{*})({{\mathbf{x}}}-{{\mathbf{z}}})^{{\boldsymbol{j}}}.$
(5.5)

If $\Upsilon$ is a mesh element $T\in{\mathcal{T}}_{h}$, which is the case of
our interest, we denote the centre of the Taylor expansion as
${{\mathbf{x}}}_{T}\in T$ and the Taylor polynomial as
${\mathsf{T}}^{p+1}_{T}[v]({{\mathbf{x}}})={\mathsf{T}}^{p+1}_{{{\mathbf{x}}}_{T}}[v]({{\mathbf{x}}})$.

For each $T\in{\mathcal{T}}_{h}$, for $m\in\mathbb{N}$ we use the standard
$C^{m}$ norms and seminorms denoted by

$\left\|v\right\|_{C^{0}(T)}:=\sup_{{{\mathbf{x}}}\in
T}|v({{\mathbf{x}}})|,\qquad\left|v\right|_{C^{m}(T)}:=\max_{|{\boldsymbol{i}}|=m}\left\|D^{\boldsymbol{i}}v\right\|_{C^{0}(T)}.$

We define the global space $C^{m}({\mathcal{T}}_{h}):=\\{v\in
L^{2}(\Omega)\mid v_{|_{T}}\in C^{m}(T)\quad\forall T\in\mathcal{T}_{h}\\}$
and the seminorms

$\left|v\right|_{C^{m}{({\mathcal{T}}_{h})}}:=\max_{T\in{\mathcal{T}}_{h}}\left|v_{|_{T}}\right|_{C^{m}(T)}.$

We now define the quasi-Trefftz space and prove its optimal approximation
properties for a general linear PDE with varying coefficients. This is
possible since no special property of the diffusion-advection-reaction
equation will be used in the proof of Theorem 5.2.1, which provides local
best-approximation bounds. Then, we will focus on the homogeneous diffusion-
advection-reaction equation of our model problem.

We consider a linear partial differential operator of order $m\in\mathbb{N}$
denoted as

${\mathcal{M}}:=\sum_{|\bm{j}|\leq m}\alpha_{\bm{j}}D^{\bm{j}},$ (5.6)

with coefficients $\alpha_{{\bm{j}}}:\Omega\to\mathbb{R}$ for
$\left|{\bm{j}}\right|\leq m$. We introduce local Trefftz and quasi-Trefftz
spaces.

######  Definition 5.1.2 (Local Trefftz space).

Let $p\in\mathbb{N}$, let $T\in{\mathcal{T}}_{h}$ and assume that, for all
$\left|{\bm{j}}\right|\leq m$, the coefficients $\alpha_{{\bm{j}}}$ of
${\mathcal{M}}$ are constant inside the mesh element $T$. We define the
polynomial Trefftz space for the homogeneous equation ${\mathcal{M}}u=0$ on
the mesh element $T$ as

$\mathbb{T}^{p}(T):=\\{v\in\mathbb{P}_{d}^{p}(T)\mid{\mathcal{M}}v=0\text{ in
}T\\}.$

This is the space of degree-$p$ polynomials such that are exact solutions of
${\mathcal{M}}u=0$ in $T$. We observe that, by definition, $\mathbb{T}^{p}(T)$
is a subspace of the full polynomial space $\mathbb{P}_{d}^{p}(T)$.

######  Remark 5.1.1.

When all the terms in the differential operator ${\mathcal{M}}$ are
derivatives of the same order, such as, for example, in the Laplace’s equation
and in the wave equation, the polynomial Trefftz space $\mathbb{T}^{p}(T)$
offers the same orders of $h$-convergence as the full polynomial space
$\mathbb{P}^{p}_{d}(T)$ [28, Lemma 1]. On the other hand, if the differential
operator ${\mathcal{M}}$ includes derivatives of different orders, in general,
the convergence rates for the Trefftz subspace of degree-$p$ polynomials are
not the same as those for full polynomial space of the same degree [15]. For
example, in [15] for the linear time-dependent Schrödinger equation the same
accuracy as for full polynomials of degree at most $p$ is reached using
Trefftz polynomials of degree $2p$. Moreover, when the derivatives are of
different orders, the polynomial Trefftz space defined above can also be
trivial, such as in the case of the Helmholtz’s equation, which does not admit
polynomial solutions. The quasi-Trefftz space, instead, is always rich enough
to give the same approximation bounds as the full polynomial space of the same
degree, as we will see below in Theorem 5.2.1. This suggests that also for
problems with piecewise-constant coefficients could be a good idea to use
quasi-Trefftz methods.

######  Definition 5.1.3 (Local quasi-Trefftz space).

Let $p\in\mathbb{N}$, let $T\in{\mathcal{T}}_{h}$ and let
${{\mathbf{x}}}_{T}\in T$. Assume that the coefficients $\alpha_{\bm{j}}\in
C^{\max\\{p-m,0\\}}(T)$ for all $\left|{\bm{j}}\right|\leq m$. We define the
quasi-Trefftz space for the homogeneous equation ${\mathcal{M}}u=0$ on the
mesh element $T$ as

${\mathbb{Q\\!T}}^{p}(T):=\\{v\in\mathbb{P}^{p}_{d}(T)\mid
D^{\bm{i}}{\mathcal{M}}v(\mathbf{x}_{T})=0,\forall\bm{i}\in\mathbb{N}^{d},|\bm{i}|\leq
p-m\\}.$ (5.7)

This is the space of degree-$p$ polynomials such that the Taylor polynomial of
order $p-m+1$ of their image by the differential operator ${\mathcal{M}}$
vanishes at ${{\mathbf{x}}}_{T}$. Notice that, by definition,
${\mathbb{Q\\!T}}^{p}(T)$ is a subspace of the full polynomial space
$\mathbb{P}_{d}^{p}(T)$ and that, for $p<m$, ${\mathbb{Q\\!T}}^{p}(T)$
coincides with $\mathbb{P}_{d}^{p}(T)$.

We notice that, in general,
${\mathbb{Q\\!T}}^{p}(T)\not\subseteq{\mathbb{Q\\!T}}^{p+1}(T)$. In order to
show this, we can consider, for example, as linear partial differential
operator ${\mathcal{M}}$ the second-order diffusion-advection-reaction
operator ${\mathcal{L}}$ with ${\bm{K}}({{\mathbf{x}}})=Id$,
${\boldsymbol{\beta}}({{\mathbf{x}}})=(1,\cdots,1)^{T}$ and
$\sigma({{\mathbf{x}}})=\frac{2}{x_{1}^{2}+1}$. If we consider the point
${{\mathbf{x}}}_{T}=\bm{0}$ and the function
$v({{\mathbf{x}}})=x_{1}^{2}+1\in\mathbb{P}^{2}_{d}(T)$, then,
${\mathcal{L}}v({{\mathbf{x}}})=-2+2x_{1}+2$, so
${\mathcal{L}}v({{\mathbf{x}}}_{T})=0$. Hence $v\in{\mathbb{Q\\!T}}^{2}(T)$ ,
but $\partial_{x_{1}}{\mathcal{L}}v({{\mathbf{x}}})=2$, implying that
$v\in{\mathbb{Q\\!T}}^{2}(T)\setminus{\mathbb{Q\\!T}}^{3}(T)$. This means
that, in general, for increasing polynomial degree $p$ the quasi-Trefftz
spaces are not nested.

For $\left|{\bm{i}}\right|\leq p-m$, using the definition of ${\mathcal{M}}$
(5.6) and the Leibniz formula (5.2), we obtain, for any $v\in
C^{\left|{\bm{i}}\right|+m}(T)$,

$\displaystyle D^{\boldsymbol{i}}{\mathcal{M}}v({{\mathbf{x}}}_{T})$
$\displaystyle=\sum_{\left|{\bm{j}}\right|\leq
m}D^{\bm{i}}\left(\alpha_{{\bm{j}}}({{\mathbf{x}}}_{T})D^{\bm{j}}v({{\mathbf{x}}}_{T})\right)$
(5.8) $\displaystyle=\sum_{\left|{\bm{j}}\right|\leq
m}\sum_{{\bm{r}}\leq{\bm{i}}}\binom{{\bm{i}}}{{\bm{r}}}D^{\bm{r}}\alpha_{{\bm{j}}}({{\mathbf{x}}}_{T})D^{{\bm{i}}-{\bm{r}}+{\bm{j}}}v({{\mathbf{x}}}_{T}).$

From (5.8), the space ${\mathbb{Q\\!T}}^{p}(T)$ is well-defined if, for all
$\left|{\bm{j}}\right|\leq m$, $\alpha_{{\bm{j}}}\in C^{\max\\{p-m,0\\}}$ in a
neighbourhood of ${{\mathbf{x}}}_{T}$.

The global quasi-Trefftz space on the whole mesh is defined as follows.

######  Definition 5.1.4 (Global quasi-Trefftz space).

For $p\in\mathbb{N}$, we define the global quasi-Trefftz space for the
homogeneous equation ${\mathcal{M}}u=0$ as

${\mathbb{Q\\!T}}^{p}(\mathcal{T}_{h}):=\\{v\in L^{2}(\Omega)\mid
v_{|_{T}}\in{\mathbb{Q\\!T}}^{p}(T)\quad\forall T\in\mathcal{T}_{h}\\}.$ (5.9)

We want to choose this space as discrete space $V_{h}$ in the DG method (3.8).

###  5.2 Approximation properties of the quasi-Trefftz space

The following theorem states the approximation properties of the local quasi-
Trefftz spaces. We prove that Taylor polynomials of smooth solutions of the
homogeneous linear PDE considered are quasi-Trefftz functions. This implies
that we can approximate smooth solutions in ${\mathbb{Q\\!T}}^{p}(T)$ with
optimal convergence rates with respect to the mesh size, meaning that the
orders of $h$-convergence for quasi-Trefftz spaces are the same as those for
full polynomial spaces.

######  Theorem 5.2.1.

Let $p\in\mathbb{N}$ and $T\in{\mathcal{T}}_{h}$. Assume that
$\alpha_{\bm{j}}\in C^{\max\\{p-m,0\\}}(T)$ for all $\left|{\bm{j}}\right|\leq
m$ and that $u\in C^{p+1}(T)$ satisfies ${\mathcal{M}}u=0$ in $T$.

Then the Taylor polynomial
$\mathsf{T}^{p+1}_{T}[u]\in{\mathbb{Q\\!T}}^{p}(T)$.

Moreover, if $T$ is star-shaped with respect to $\mathbf{x}_{T}$, then, for
all $q\in\mathbb{N}$ such that $q\leq p$,

$\inf_{v\in{\mathbb{Q\\!T}}^{p}(T)}\left|u-v\right|_{C^{q}(T)}\leq\frac{d^{p+1-q}}{(p+1-q)!}h_{T}^{p+1-q}\left|u\right|_{C^{p+1}(T)}.$
(5.10)

###### Proof.

First we prove that $\mathsf{T}^{p+1}_{T}[u]\in{\mathbb{Q\\!T}}^{p}(T)$. By
definition, we already know that
$\mathsf{T}^{p+1}_{T}[u]\in\mathbb{P}^{p}_{d}(T)$. For all
$|{\boldsymbol{i}}|\leq p-m$, we have

$\displaystyle
D^{\boldsymbol{i}}{\mathcal{M}}\mathsf{T}^{p+1}_{T}[u]({{\mathbf{x}}}_{T})$
$\displaystyle=\sum_{\left|{\bm{j}}\right|\leq
m}\sum_{{\bm{r}}\leq{\bm{i}}}\binom{{\bm{i}}}{{\bm{r}}}D^{\bm{r}}\alpha_{{\bm{j}}}({{\mathbf{x}}}_{T})D^{{\bm{i}}-{\bm{r}}+{\bm{j}}}\mathsf{T}^{p+1}_{T}[u]({{\mathbf{x}}}_{T})$
$\displaystyle=\sum_{\left|{\bm{j}}\right|\leq
m}\sum_{{\bm{r}}\leq{\bm{i}}}\binom{{\bm{i}}}{{\bm{r}}}D^{\bm{r}}\alpha_{{\bm{j}}}({{\mathbf{x}}}_{T})D^{{\bm{i}}-{\bm{r}}+{\bm{j}}}u({{\mathbf{x}}}_{T})$
$\displaystyle=D^{\boldsymbol{i}}{\mathcal{M}}u({{\mathbf{x}}}_{T}).$

The first step is the identity (5.8) with $v=\mathsf{T}^{p+1}_{T}[u]\in
C^{\left|{\bm{i}}\right|+m}(T)$, while the second step follows from the
property (5.4) with partial derivatives of order at most equal to
$\left|{\bm{i}}\right|+m\leq p$, and the third one is (5.8) again with $v=u\in
C^{\left|{\bm{i}}\right|+m}(T)$. Since $u$ is solution of ${\mathcal{M}}u=0$
in $T$, the last term $D^{\boldsymbol{i}}{\mathcal{M}}u({{\mathbf{x}}}_{T})$
vanishes, implying that
$D^{\boldsymbol{i}}{\mathcal{M}}\mathsf{T}^{p+1}_{T}[u](\mathbf{x}_{T})=0$ for
all $|{\boldsymbol{i}}|\leq p-m$. This shows that the Taylor polynomial
$\mathsf{T}^{p+1}_{T}[u]$ belongs to the quasi-Trefftz space
${\mathbb{Q\\!T}}^{p}(T)$ defined in (5.7).

In order to prove the inequality (5.10), fix $q\in\mathbb{N}$ such that $q\leq
p$. Using the fact that $\mathsf{T}^{p+1}_{T}[u]\in{\mathbb{Q\\!T}}^{p}(T)$,
the best approximation error is estimated by the error of the Taylor
polynomial

$\inf_{v\in{\mathbb{Q\\!T}}^{p}(T)}\left|u-v\right|_{C^{q}(T)}\leq\left|u-\mathsf{T}^{p+1}_{T}[u]\right|_{C^{q}(T)}.$

Using the $\left|\cdot\right|_{C^{q}}$-norm definition, the linearity of the
derivatives and the identity
$D^{\boldsymbol{i}}{\mathsf{T}}^{p+1}_{T}[u]={\mathsf{T}}^{p+1-|{\boldsymbol{i}}|}_{T}[D^{\boldsymbol{i}}u]$
for $|{\boldsymbol{i}}|=q\leq p$ from (5.3), we obtain

$\displaystyle\left|u-\mathsf{T}^{p+1}_{T}[u]\right|_{C^{q}(T)}$
$\displaystyle=\max_{{\boldsymbol{i}}\in\mathbb{N}^{d},\;|{\boldsymbol{i}}|=q}\left\|D^{\boldsymbol{i}}(u-\mathsf{T}^{p+1}_{T}[u])\right\|_{C^{0}(T)}$
$\displaystyle=\max_{{\boldsymbol{i}}\in\mathbb{N}^{d},\;|{\boldsymbol{i}}|=q}\left\|D^{\boldsymbol{i}}u-\mathsf{T}^{p+1-q}_{T}[D^{\boldsymbol{i}}u]\right\|_{C^{0}(T)}.$

Since $u\in C^{p+1}(T)$ and $T$ is star-shaped with respect to
$\mathbf{x}_{T}$, we can consider the Lagrange’s form of the Taylor remainder
(5.5) and estimate the last term as follows:

$\displaystyle\max_{{\boldsymbol{i}}\in\mathbb{N}^{d},\;|{\boldsymbol{i}}|=q}\left\|D^{\boldsymbol{i}}u-\mathsf{T}^{p+1-q}_{T}[D^{\boldsymbol{i}}u]\right\|_{C^{0}(T)}$
$\displaystyle\leq\max_{{\boldsymbol{i}}\in\mathbb{N}^{d},\;|{\boldsymbol{i}}|=q}\sum_{\begin{subarray}{c}|{\boldsymbol{j}}|=p+1-q\end{subarray}}\frac{1}{{\boldsymbol{j}}!}\sup_{{{\mathbf{x}}}\in
T}\left|D^{{\boldsymbol{i}}+{\boldsymbol{j}}}u({{\mathbf{x}}}_{*})({{\mathbf{x}}}-{{\mathbf{x}}}_{T})^{\boldsymbol{j}}\right|.$
$\displaystyle\leq\frac{d^{p+1-q}}{(p+1-q)!}h_{T}^{p+1-q}\left|u\right|_{C^{p+1}(T)},$

where we use the formula
$\sum_{|{\boldsymbol{j}}|=p+1-q}\frac{1}{{\boldsymbol{j}}!}=\frac{d^{p+1-q}}{(p+1-q)!}$,
obtained from (5.1) with $k=p+1-q$. From the chain of inequalities, the thesis
follows. ∎

###  5.3 Convergence analysis of the quasi-Trefftz DG method

####  5.3.1 Quasi-Trefftz space for the homogenous diffusion-advection-
reaction equation

We now focus on the diffusion-advection-reaction operator ${\mathcal{L}}$ with
variable coefficients ${\bm{K}}:\Omega\to\mathbb{R}^{d\times d}$,
${\boldsymbol{\beta}}:\Omega\to\mathbb{R}^{d}$ and
$\sigma:\Omega\to\mathbb{R}$:

$\mathcal{L}v=\text{div}\left(-\bm{K}\nabla v+\bm{\beta}v\right)+\sigma v.$

Since ${\mathcal{L}}$ is a linear partial differential operator of the second
order, the quasi-Trefftz space (5.7) for the homogeneous diffusion-advection-
reaction equation on a mesh element $T\in\mathcal{T}_{h}$ is

${\mathbb{Q\\!T}}^{p}(T)=\\{v\in\mathbb{P}^{p}_{d}(T)\mid
D^{\bm{i}}(\mathcal{L}v)(\mathbf{x}_{T})=0,\forall\bm{i}\in\mathbb{N}^{d},|\bm{i}|\leq
p-2\\},\quad p\in\mathbb{N}.$ (5.11)

As seen before, this space is well-defined if the PDE coefficients are
sufficiently smooth. The operator ${\mathcal{L}}$ can be written in the form
(5.6) with $m=2$, for which we have shown that the coefficients
$\alpha_{\bm{j}}$ need to belong to $C^{\max\\{p-2,0\\}}(T)$. In order to
deduce the regularity assumptions of ${\bm{K}}$, ${\boldsymbol{\beta}}$ and
$\sigma$, under which the quasi-Trefftz space is well-defined, we expand the
operator ${\mathcal{L}}$, using linearity and Leibniz formula (5.2):

$\displaystyle\mathcal{L}v$
$\displaystyle=\sum_{j=1}^{d}D^{\bm{e}_{j}}(-\bm{K}\nabla
v+\bm{\beta}v)_{j}+\sigma v$
$\displaystyle=\sum_{j=1}^{d}\biggl{(}-\sum_{m=1}^{d}D^{\bm{e}_{j}}\big{(}\bm{K}_{jm}D^{\bm{e}_{m}}v\big{)}+D^{\bm{e}_{j}}(\bm{\beta}_{j}v)\biggl{)}+\sigma
v$
$\displaystyle=\sum_{j=1}^{d}\sum_{\begin{subarray}{c}\bm{\ell}\in\mathbb{N}^{d},\;\bm{\ell}\leq\bm{e}_{j}\end{subarray}}\binom{\bm{e}_{j}}{\bm{\ell}}\biggl{(}-\sum_{m=1}^{d}D^{\bm{\ell}}\bm{K}_{jm}D^{\bm{e}_{j}-\bm{\ell}+\bm{e}_{m}}v+D^{\bm{\ell}}\bm{\beta}_{j}D^{\bm{e}_{j}-\bm{\ell}}v\biggl{)}+\sigma
v$
$\displaystyle=\sum_{j=1}^{d}\biggl{[}-\sum_{m=1}^{d}\bigg{(}\bm{K}_{jm}D^{\bm{e}_{j}+\bm{e}_{m}}v+D^{\bm{e}_{j}}\bm{K}_{jm}D^{\bm{e}_{m}}v\bigg{)}+\bm{\beta}_{j}D^{\bm{e}_{j}}v+D^{\bm{e}_{j}}{\boldsymbol{\beta}}_{j}v\biggl{]}+\sigma
v.$

It follows that $D^{\bm{\ell}}\bm{K}_{jm}$, $D^{\bm{\ell}}\bm{\beta}_{j}$ and
$\sigma$ need to belong to $C^{\max\\{p-2,0\\}}(T)$ for all
$\bm{\ell}\leq\bm{e}_{j}$ and for all $j,m=1,\dots,d$. Hence, the quasi-
Trefftz space ${\mathbb{Q\\!T}}^{p}(T)$ for the diffusion-advection-reaction
equation is well-defined if ${\bm{K}}_{jm}\in C^{\max\\{p-1,1\\}}(T)$,
${\boldsymbol{\beta}}_{j}\in C^{\max\\{p-1,1\\}}(T)$ and $\sigma\in
C^{\max\\{p-2,0\\}}(T)$ for all $j,m=1,\dots,d$.

Using linearity and Leibniz formula (5.2), we also derive the following
expression of the partial derivatives of ${\mathcal{L}}v$ which will be useful
for the construction of the quasi-Trefftz basis functions:

$\displaystyle D^{\bm{i}}\mathcal{L}v=$
$\displaystyle-\sum_{j=1}^{d}\sum_{m=1}^{d}D^{\bm{i}+\bm{e}_{j}}(\bm{K}_{jm}D^{\bm{e}_{m}}v)+\sum_{j=1}^{d}D^{{\bm{i}}+\bm{e}_{j}}(\bm{\beta}_{j}v)+D^{{\bm{i}}}\sigma
v$ $\displaystyle=$
$\displaystyle-\sum_{j=1}^{d}\sum_{m=1}^{d}\sum_{\bm{\ell}\leq\bm{i}+\bm{e}_{j}}\binom{\bm{i}+\bm{e}_{j}}{\bm{\ell}}D^{\bm{\ell}}\bm{K}_{jm}D^{\bm{i}+\bm{e}_{j}-\bm{\ell}+\bm{e}_{m}}v$
(5.12)
$\displaystyle+\sum_{j=1}^{d}\sum_{\bm{\ell}\leq\bm{i}+\bm{e}_{j}}\binom{\bm{i}+\bm{e}_{j}}{\bm{\ell}}D^{\bm{\ell}}\bm{\beta}_{j}D^{\bm{i}+\bm{e}_{j}-\bm{\ell}}v+\sum_{\bm{\ell}\leq\bm{i}}\binom{\bm{i}}{\bm{\ell}}D^{\bm{\ell}}\sigma
D^{\bm{i}-\bm{\ell}}v.$

####  5.3.2 h-convergence of the quasi-Trefftz DG method

We consider only homogenous diffusion-advection-reaction equation, i.e, with
no source term $f$, since the quasi-Trefftz space defined in (5.11) is suited
for this kind of problems. Choosing the global quasi-Trefftz space of degree
$p$ defined in (5.9) as discrete space $V_{h}$ in the DG scheme (3.8) with
$f=0$, we obtain the quasi-Trefftz DG method.

The next theorem establishes the convergence rate in $h$ for the error
$u-u_{h}$ in the $\left\|\cdot\right\|_{dar}$-norm. From the quasi-optimality
inequality (4.6) we deduce the convergence rate for the quasi-Trefftz DG
method, using the the error bound (5.10). We recall that
$V_{*}:=H^{1}_{\Gamma_{D}}(\Omega)\cap H^{2}({\mathcal{T}}_{h})$.

######  Theorem 5.3.1 (Quasi-Trefftz DG convergence rate).

Let $p\in\mathbb{N}$ and let $u\in V_{*}\cap C^{p+1}({\mathcal{T}}_{h})$ solve
(1.6)-(1.8) with ${\bm{K}}_{jm}\in C^{\max\\{p-1,1\\}}({\mathcal{T}}_{h})$,
${\boldsymbol{\beta}}_{j}\in C^{\max\\{p-1,1\\}}({\mathcal{T}}_{h})$ and
$\sigma\in C^{\max\\{p-2,0\\}}({\mathcal{T}}_{h})$ for all $j,m=1,\dots,d$ and
$f=0$. Let $u_{h}$ solve (3.8) with discrete space
$V_{h}={\mathbb{Q\\!T}}^{p}({\mathcal{T}}_{h})$ and the penalty parameter
$\gamma$ as in Lemma 4.2.1. Assume that each mesh element $T$ is star-shaped
with respect to $\mathbf{x}_{T}$. Then, the following convergence rate holds

$\begin{split}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-u_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{dar}\leq&\left(1+\frac{M}{\alpha}\right)\frac{d^{p}}{p!}\left[\sum_{T\in{\mathcal{T}}_{h}}\left|T\right|\biggl{(}\left(1+\eta\right)\left\|{\bm{K}}\right\|_{L^{\infty}(T)}\right.\\\
&\left.+\frac{d^{2}}{(p+1)^{2}}\left(2C_{g}\gamma\eta+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}\eta
h_{T}+h_{T}^{2}+\eta
h_{T}\right)\biggr{)}h_{T}^{2p}\left|u\right|^{2}_{C^{p+1}(T)}\right]^{\frac{1}{2}},\end{split}$
(5.13)

with $\alpha$ as in Lemma 4.2.3 and $M$ as in Lemma 4.3.3.

###### Proof.

Since we have the quasi-optimality inequality (4.6), we now focus on
estimating the quantity
$\inf_{v_{h}\in{\mathbb{Q\\!T}}^{p}({\mathcal{T}}_{h})}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{dar,*}$.
We recall, for $v\in V_{*h}=V_{*}+{\mathbb{Q\\!T}}^{p}({\mathcal{T}}_{h})$,
the definition of the $\left\|\cdot\right\|_{dar,*}$-norm defined in Section
4.1

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}=$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\int_{T}{\bm{K}}\nabla v\cdot\nabla
v+\sum_{e\in{\mathcal{E}}_{h}^{I}\cup{\mathcal{E}}_{h}^{D}}\frac{\gamma}{h_{e}}\int_{e}[\\![v]\\!]^{2}+\frac{1}{2}\sum_{e\in{\mathcal{E}}_{h}}\int_{e}\left|{\boldsymbol{\beta}}\cdot{{\mathbf{n}}}_{e}\right|[\\![v]\\!]^{2}$
$\displaystyle+\left\|v\right\|^{2}_{L^{2}(\Omega)}+\sum_{T\in{\mathcal{T}}_{h}}h_{T}\left\|{\bm{K}}^{\frac{1}{2}}\nabla
v\cdot{{\mathbf{n}}}_{T}\right\|^{2}_{L^{2}(\partial
T)}+\sum_{T\in{\mathcal{T}}_{h}}\left\|v\right\|^{2}_{L^{2}(\partial T)}.$

We recall the assumption (2.2): the mesh ${\mathcal{T}}_{h}$ is graded, i.e.,
there exists a constant $C_{g}>0$ such that, for all $T\in{\mathcal{T}}_{h}$
and all $e\in{\mathcal{F}}_{T}$,

$C_{g}h_{e}\geq h_{T}.$

Observe that, from the Young’s inequality,
$[\\![v]\\!]^{2}=(v_{1}-v_{2})^{2}\leq 2(v_{1}^{2}+v_{2}^{2})$ on an internal
facet $e=\partial T_{1}\cap\partial T_{2}$ and $[\\![v]\\!]^{2}=v^{2}$ on a
boundary facet $e$. Using these two facts in the second and third terms of the
norm, we rearrange the sums over parts of the mesh skeleton as sums over
elements and obtain the bound:

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}\leq$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\biggl{(}\left\|{\bm{K}}\right\|_{L^{\infty}(T)}\left\|\nabla
v\right\|^{2}_{L^{2}(T)}+2C_{g}\frac{\gamma}{h_{T}}\left\|v\right\|^{2}_{L^{2}(\partial
T)}+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}\left\|v\right\|^{2}_{L^{2}(\partial
T)}$
$\displaystyle+\left\|v\right\|^{2}_{L^{2}(T)}+\left\|{\bm{K}}\right\|_{L^{\infty}(T)}h_{T}\left\|\nabla
v\right\|^{2}_{L^{2}(\partial T)}+\left\|v\right\|^{2}_{L^{2}(\partial
T)}\biggr{)}.$

Next, we use the definition of the $\left\|\cdot\right\|_{C^{m}}$-norms and
obtain

$\displaystyle{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|v\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}\leq$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\biggl{(}\left\|{\bm{K}}\right\|_{L^{\infty}(T)}\left|T\right|\left\|\nabla
v\right\|^{2}_{C^{0}(T)}+2C_{g}\frac{\gamma}{h_{T}}\left|\partial
T\right|\left\|v\right\|^{2}_{C^{0}(T)}$
$\displaystyle+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}\left|\partial
T\right|\left\|v\right\|^{2}_{C^{0}(T)}+\left|T\right|\left\|v\right\|^{2}_{C^{0}(T)}$
$\displaystyle+\left|\partial
T\right|\left\|v\right\|^{2}_{C^{0}(T)}+\left\|{\bm{K}}\right\|_{L^{\infty}(T)}\left|\partial
T\right|h_{T}\left\|\nabla v\right\|^{2}_{C^{0}(T)}\biggr{)}.$

Considering the quantity of our interest we have

$\displaystyle\inf_{v_{h}\in{\mathbb{Q\\!T}}^{p}({\mathcal{T}}_{h})}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}\leq$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\inf_{v_{h}\in{\mathbb{Q\\!T}}^{p}(T)}\biggl{[}\left(\left|T\right|+\left|\partial
T\right|h_{T}\right)\left\|{\bm{K}}\right\|_{L^{\infty}(T)}\left\|\nabla(u-v_{h})\right\|^{2}_{C^{0}(T)}$
$\displaystyle+\left(2C_{g}\frac{\gamma}{h_{T}}\left|\partial
T\right|+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}\left|\partial
T\right|+\left|T\right|+\left|\partial
T\right|\right)\left\|u-v_{h}\right\|^{2}_{C^{0}(T)}\biggr{]}.$

Using the Quasi-Trefftz approximation properties (5.10) we get

$\displaystyle\inf_{v_{h}\in{\mathbb{Q\\!T}}^{p}({\mathcal{T}}_{h})}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}\leq$
$\displaystyle\sum_{T\in{\mathcal{T}}_{h}}\left[\left(\left|T\right|+\left|\partial
T\right|h_{T}\right)\left\|{\bm{K}}\right\|_{L^{\infty}(T)}\frac{d^{2p}}{(p!)^{2}}h_{T}^{2p}\left|u\right|^{2}_{C^{p+1}(T)}\right.$
$\displaystyle\left.+\left(2C_{g}\frac{\gamma}{h_{T}}\left|\partial
T\right|+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}\left|\partial
T\right|+\left|T\right|+\left|\partial T\right|\right)\times\right.$
$\displaystyle\left.\frac{d^{2(p+1)}}{((p+1)!)^{2}}h_{T}^{2(p+1)}\left|u\right|^{2}_{C^{p+1}(T)}\right].$

We now use the assumption (2.3): exists $\eta>0$ such that
$h_{T}\left|\partial T\right|\leq\eta\left|T\right|$ for all
$T\in{\mathcal{T}}_{h}$, and we infer that

$\displaystyle\inf_{v_{h}\in{\mathbb{Q\\!T}}^{p}({\mathcal{T}}_{h})}{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|u-v_{h}\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}^{2}_{dar,*}\leq$
$\displaystyle\frac{d^{2p}}{(p!)^{2}}\sum_{T\in{\mathcal{T}}_{h}}\biggl{[}\left(1+\eta\right)\left|T\right|\left\|{\bm{K}}\right\|_{L^{\infty}(T)}+\frac{d^{2}}{(p+1)^{2}}\times$
$\displaystyle\left(2C_{g}\frac{\gamma}{h_{T}}\eta
h_{T}^{-1}+\left\|{\boldsymbol{\beta}}\right\|_{L^{\infty}(T)}\eta
h_{T}^{-1}+1+\eta
h_{T}^{-1}\right)\left|T\right|h_{T}^{2}\biggr{]}h_{T}^{2p}\left|u\right|^{2}_{C^{p+1}(T)}.$

Combining this bound with the inequality (4.6) yields the assertion. ∎

We remark that we use the $C^{m}$-norms in the convergence analysis of the
quasi-Trefftz DG method, instead of the classical Sobolev norms $H^{m}$. We
observe that in Theorem 5.2.1 and in Theorem 5.3.1 the solution $u$ needs to
belong to $C^{p+1}({\mathcal{T}}_{h})$, which is a stronger regularity
assumption than the one required in Theorem 4.4.3, i.e., $u\in
H^{p+1}({\mathcal{T}}_{h})$. For the Trefftz space the analysis has been
extended to the case of less regular solutions using the fact that “averaged
Taylor polynomials” [28, p. 421] of solutions are Trefftz functions. However,
we cannot use this argument since, in general, they are not quasi-Trefftz
functions and to our knowledge a study using Sobolev norms for the quasi-
Trefftz space is still missing [23, Remark 4.7]. Apart from this difference,
we have optimal convergence rates in the $\left\|\cdot\right\|_{dar}$-norm.

###  5.4 Basis functions

In this section, we define, for any $p\in\mathbb{N}$ and
$T\in{\mathcal{T}}_{h}$, a family of bases for the quasi-Trefftz space
${\mathbb{Q\\!T}}^{p}(T)$ for the homogeneous diffusion-advection-reaction
equation and derive an explicit formula to compute them.

Since each function $v\in{\mathbb{Q\\!T}}^{p}(T)$ is a degree-$p$ polynomial,
$v$ can be expressed as a linear combination of scaled monomials centred in
${{\mathbf{x}}}_{T}\in T$:

$v({{\mathbf{x}}})=\sum_{{\boldsymbol{k}}\in\mathbb{N}^{d},|{\boldsymbol{k}}|\leq
p}a_{\boldsymbol{k}}\left(\frac{\mathbf{x}-\mathbf{x}_{T}}{h_{T}}\right)^{{\boldsymbol{k}}},$
(5.14)

where the coefficients
$\\{a_{\boldsymbol{k}}\\}_{{\boldsymbol{k}}\in\mathbb{N}^{d},|{\boldsymbol{k}}|\leq
p}$ are such that

$a_{\boldsymbol{k}}=\frac{h_{T}^{\left|\bm{k}\right|}}{{\boldsymbol{k}}!}D^{\boldsymbol{k}}v({{\mathbf{x}}}_{T}).$

This is obtained by considering the derivatives of $v$

$D^{{\bm{i}}}v({{\mathbf{x}}})=\sum_{\begin{subarray}{c}{\boldsymbol{k}}\in\mathbb{N}^{d},\left|{\boldsymbol{k}}\right|\leq
p\\\
{\boldsymbol{k}}\geq{\bm{i}}\end{subarray}}a_{\boldsymbol{k}}\left(\frac{\mathbf{x}-\mathbf{x}_{T}}{h_{T}}\right)^{{\boldsymbol{k}}-{\bm{i}}}\frac{{\boldsymbol{k}}!}{({\boldsymbol{k}}-{\bm{i}})!}\frac{1}{h_{T}^{\left|{\bm{i}}\right|}},$

and evaluating them at ${{\mathbf{x}}}_{T}$. The construction of
$v\in{\mathbb{Q\\!T}}^{p}(T)$ as in (5.14) is equivalent to finding the set of
coefficients
$\\{a_{\boldsymbol{k}}\\}_{{\boldsymbol{k}}\in\mathbb{N}^{d},|{\boldsymbol{k}}|\leq
p}$ such that $v$ satisfies $D^{\bm{i}}\mathcal{L}v({{\mathbf{x}}}_{T})=0$ for
all $\left|{\bm{i}}\right|\leq p-2$. Evaluating (5.3.1) in
${{\mathbf{x}}}_{T}$ and using the relation
$D^{\bm{k}}v({{\mathbf{x}}}_{T})=\frac{a_{\bm{k}}\bm{k}!}{h_{T}^{\left|\bm{k}\right|}}$,
we obtain:

$\displaystyle D^{\bm{i}}\mathcal{L}v({{\mathbf{x}}}_{T})=$
$\displaystyle-\sum_{j=1}^{d}\sum_{m=1}^{d}\sum_{\bm{\ell}\leq\bm{i}+\bm{e}_{j}}\binom{\bm{i}+\bm{e}_{j}}{\bm{\ell}}D^{\bm{\ell}}\bm{K}_{jm}({{\mathbf{x}}}_{T})\frac{a_{\bm{i}+\bm{e}_{j}-\bm{\ell}+\bm{e}_{m}}(\bm{i}+\bm{e}_{j}-\bm{\ell}+\bm{e}_{m})!}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|+2}}$
$\displaystyle+\sum_{j=1}^{d}\sum_{\bm{\ell}\leq\bm{i}+\bm{e}_{j}}\binom{\bm{i}+\bm{e}_{j}}{\bm{\ell}}D^{\bm{\ell}}\bm{\beta}_{j}({{\mathbf{x}}}_{T})\frac{a_{\bm{i}+\bm{e}_{j}-\bm{\ell}}(\bm{i}+\bm{e}_{j}-\bm{\ell})!}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|+1}}$
$\displaystyle+\sum_{\bm{\ell}\leq\bm{i}}\binom{\bm{i}}{\bm{\ell}}D^{\bm{\ell}}\sigma({{\mathbf{x}}}_{T})\frac{a_{\bm{i}-\bm{\ell}}(\bm{i}-\bm{\ell})!}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|}}.$

Using the definition of the binomial coefficients, this can be rewritten as

$\displaystyle D^{\bm{i}}\mathcal{L}v({{\mathbf{x}}}_{T})=$
$\displaystyle-\sum_{j=1}^{d}\sum_{m=1}^{d}\sum_{\bm{\ell}\leq\bm{i}+\bm{e}_{j}}\frac{(\bm{i}+\bm{e}_{j})!}{\bm{\ell}!}D^{\bm{\ell}}\bm{K}_{jm}({{\mathbf{x}}}_{T})\frac{a_{\bm{i}+\bm{e}_{j}-\bm{\ell}+\bm{e}_{m}}(i_{m}+(\bm{e}_{j})_{m}-\ell_{m}+1)}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|+2}}$
$\displaystyle+\sum_{j=1}^{d}\sum_{\bm{\ell}\leq\bm{i}+\bm{e}_{j}}\frac{(\bm{i}+\bm{e}_{j})!}{\bm{\ell}!}D^{\bm{\ell}}\bm{\beta}_{j}({{\mathbf{x}}}_{T})\frac{a_{\bm{i}+\bm{e}_{j}-\bm{\ell}}}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|+1}}$
(5.15)
$\displaystyle+\sum_{\bm{\ell}\leq\bm{i}}\frac{{\bm{i}}!}{\bm{\ell}!}D^{\bm{\ell}}\sigma({{\mathbf{x}}}_{T})\frac{a_{\bm{i}-\bm{\ell}}}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|}}.$

We recall the assumption (1.11):

${\bm{K}}_{11}({{\mathbf{x}}})\geq k_{min}>0\quad\text{for
all}\quad{{\mathbf{x}}}\in\Omega.$

From the conditions $D^{\bm{i}}\mathcal{L}v({{\mathbf{x}}}_{T})=0$ for all
$\left|{\bm{i}}\right|\leq p-2$, we can compute the element with
$\bm{\ell}=\bm{0}$, $j=1$ and $m=1$ in the first sum of (5.4), dividing by its
non-zero coefficient ${\bm{K}}_{11}({{\mathbf{x}}}_{T})>0$:

$\displaystyle a_{{\bm{i}}+2{{\mathbf{e}}}_{1}}=$
$\displaystyle\frac{{h_{T}^{\left|{\bm{i}}\right|+2}}}{{\bm{K}}_{11}({{\mathbf{x}}}_{T})({\bm{i}}+2{{\mathbf{e}}}_{1})!}\times$
$\displaystyle\Biggl{(}-\sum_{\begin{subarray}{c}j,m=1,\dots,d\\\
\bm{\ell}\leq\bm{i}+\bm{e}_{j}\\\
(j,m,\bm{\ell})\neq(1,1,\bm{0})\end{subarray}}\frac{(\bm{i}+\bm{e}_{j})!}{\bm{\ell}!}D^{\bm{\ell}}\bm{K}_{jm}({{\mathbf{x}}}_{T})\frac{a_{\bm{i}+\bm{e}_{j}-\bm{\ell}+\bm{e}_{m}}(i_{m}+(\bm{e}_{j})_{m}-\ell_{m}+1)}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|+2}}$
$\displaystyle+\sum_{\begin{subarray}{c}j=1,\dots,d\\\
\bm{\ell}\leq\bm{i}+\bm{e}_{j}\end{subarray}}\frac{(\bm{i}+\bm{e}_{j})!}{\bm{\ell}!}D^{\bm{\ell}}\bm{\beta}_{j}({{\mathbf{x}}}_{T})\frac{a_{\bm{i}+\bm{e}_{j}-\bm{\ell}}}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|+1}}$
(5.16)
$\displaystyle+\sum_{\bm{\ell}\leq\bm{i}}\frac{{\bm{i}}!}{\bm{\ell}!}D^{\bm{\ell}}\sigma({{\mathbf{x}}}_{T})\frac{a_{\bm{i}-\bm{\ell}}}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|}}\Biggr{)}.$

We have found this recurrence relation (5.4) between the coefficients
$a_{{\boldsymbol{k}}}$, which will be used to compute the coefficients
iteratively. We will describe an algorithm in Section 5.5 that at every steps
relies only on values already computed and is initialized by the coefficients
$a_{{\boldsymbol{k}}}$ with $k_{1}\in\\{0,1\\}$. It is clear from (5.4) that
the recursive formula is based only on the partial derivatives at
${{\mathbf{x}}}_{T}$ of the PDE coefficients ${\bm{K}}$,
${\boldsymbol{\beta}}$ and $\sigma$. We note that the equations (5.4) could
also be written as a square linear system with $a_{{\boldsymbol{k}}}$ with
$k_{1}\neq 0,1$ as unknowns and $a_{{\boldsymbol{k}}}$ with
$k_{1}\in\\{0,1\\}$ as right-hand side. However, we will not follow this path
but a recursive implementation.

We use the relation (5.4) to construct a basis for the local quasi-Trefftz
space. We recall the dimension of the space of polynomials of maximum degree
$p$ in $d$ variables defined in 2.2.7:

$S_{d,p}=\text{dim}(\mathbb{P}^{p}(\mathbb{R}^{d}))=\binom{p+d}{d}.$

To initialize the recursion, we choose two $2$ polynomial bases:

$\left\\{\widehat{f}_{s}\right\\}_{s=1,\dots,S_{d-1,p}}\text{ basis for
}\mathbb{P}^{p}(\mathbb{R}^{d-1})\quad\text{ and
}\quad\left\\{\widetilde{f}_{s}\right\\}_{s=1,\dots,S_{d-1,p-1}}\text{ basis
for }\mathbb{P}^{p-1}(\mathbb{R}^{d-1}).$

Their total cardinality is

$N_{d,p}:=S_{d-1,p}+S_{d-1,p-1}=\binom{p+d-1}{d-1}+\binom{p+d-2}{d-1}=\frac{(p+d-2)!(2p+d-1)}{(d-1)!p!}.$

We define the following set of $N_{d,p}$ elements of
${\mathbb{Q\\!T}}^{p}(T)$:

$\displaystyle\mathcal{B}_{T}^{p}:=\left\\{\begin{array}[]{ll}b_{J}\in{\mathbb{Q\\!T}}^{p}(T):\\\
b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)=\widehat{f}_{J}\;\text{ and
}\;\partial_{x_{1}}b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)=0\text{ for }J\leq
S_{d-1,p},\\\ b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)=0\;\text{ and
}\;\partial_{x_{1}}b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)=\widetilde{f}_{J-S_{d-1,p}}\text{
for }J\geq S_{d-1,p}+1.\end{array}\right\\}_{J=1,\dots,N_{p,d},}$

where $b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)$ denotes the restriction of
$b_{J}$ to the hyperplane $\\{x_{1}=({{\mathbf{x}}}_{T})_{1}\\}$.
Equivalently, the elements $b_{J}$ are degree-$p$ polynomials that satisfy the
following system of equations

$\displaystyle\begin{cases}D^{\boldsymbol{i}}{\mathcal{L}}b_{J}({{\mathbf{x}}}_{T})=0&\quad
J=1,\ldots,N_{d,p},\;{\boldsymbol{i}}\in\mathbb{N}^{d},\;|{\boldsymbol{i}}|\leq
p-m,\\\ b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)=\widehat{f}_{J}\text{ and
}\partial_{x_{1}}b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)=0&\quad J=1,\ldots
S_{d-1,p},\\\ b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)=0\text{ and
}\partial_{x_{1}}b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)=\widetilde{f}_{J-S_{d-1,p}}&\quad
J=S_{d-1,p}+1,\ldots,N_{d,p}.\end{cases}$ (5.17)

We observe that the second and third set of conditions in (5.17) determine the
coefficients $a_{{\boldsymbol{k}}}$ with $k_{1}\in\\{0,1\\}$ of the monomial
expansion (5.14) of an element $b_{J}$. In the next proposition we show that
the elements $b_{J}$ are well-defined and that they constitute a basis of
${\mathbb{Q\\!T}}^{p}(T)$.

######  Proposition 5.4.1.

For any $p\in\mathbb{N}$ and $T\in{\mathcal{T}}_{h}$, the set
$\mathcal{B}^{p}_{T}$ constitutes a basis for the space
${\mathbb{Q\\!T}}^{p}(T)$.

###### Proof.

We need to verify that the elements $b_{J}$ are well defined and that
$\mathcal{B}_{T}^{p}$ is a set of linearly independent generators. The
algorithms described in Section 5.5 show that the coefficients
$a_{{\boldsymbol{k}}}$ of an element $b_{J}$ are uniquely determined once the
coefficients $a_{{\boldsymbol{k}}}$ with $k_{1}\in\\{0,1\\}$ are assigned.
This implies that the $N_{d,p}$ polynomials $b_{J}$ are well and uniquely
defined, since the second and third set of conditions in (5.17) fix the
coefficients $a_{0,k_{2},\dots,k_{d}}$ and $a_{1,k_{2},\dots,k_{d}}$. We also
notice that, since (5.4) holds for any $v\in{\mathbb{Q\\!T}}^{p}(T)$, we can
apply the algorithms described in Section 5.5 to construct $v$. From this we
deduce that each $v\in{\mathbb{Q\\!T}}^{p}(T)$ is uniquely determined by the
coefficients $a_{0,k_{2},\dots,k_{d}}$ and $a_{1,k_{2},\dots,k_{d}}$ and hence
by its restriction $v(({{\mathbf{x}}}_{T})_{1},\cdot)$ and the restriction of
its derivative $\partial_{x_{1}}v(({{\mathbf{x}}}_{T})_{1},\cdot)$.

For any $v\in{\mathbb{Q\\!T}}^{p}(T)$, since its restriction
$v\left(({{\mathbf{x}}}_{T})_{1},\cdot\right)$ is a polynomial of degree $p$
and the restriction of its derivative
$\partial_{x_{1}}v(({{\mathbf{x}}}_{T})_{1},\cdot)$ is a polynomial of degree
$p-1$, there exist some coefficients
$\\{\lambda_{J}\\}_{J=1}^{N_{p,d}}\subseteq\mathbb{R}$ such that

$\displaystyle v(\left({{\mathbf{x}}}_{T})_{1},\cdot\right)$
$\displaystyle=\sum_{J=1}^{S_{d-1,p}}\lambda_{J}\widehat{f}_{J}=\sum_{J=1}^{S_{d-1,p}}\lambda_{J}b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot),$
$\displaystyle\partial_{x_{1}}v(({{\mathbf{x}}}_{T})_{1},\cdot)$
$\displaystyle=\sum_{J=S_{d-1,p}+1}^{N_{d,p}}\lambda_{J}\widetilde{f}_{J-S_{d-1,p}}=\sum_{J=S_{d-1,p}+1}^{N_{d,p}}\lambda_{J}\partial_{x_{1}}b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot).$

Hence, $v=\sum_{J=1}^{N_{d,p}}\lambda_{J}b_{J}$ is a linear combination of
$b_{J}$.

The fact that the polynomials $\\{b_{J}\\}_{J=1,\dots,N_{d,p}}$ are linearly
independent is due to their restrictions to $x_{1}=({{\mathbf{x}}}_{T})_{1}$.
Assume that $\sum_{J=1}^{N_{d,p}}c_{J}b_{J}=0$ for some coefficients
$\\{c_{J}\\}_{J=1}^{N_{d,p}}\subseteq\mathbb{R}$. Then, restricting to
$x_{1}=({{\mathbf{x}}}_{T})_{1}$, we obtain

$\displaystyle\sum_{J=1}^{S_{d-1,p}}c_{J}b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)=$
$\displaystyle\sum_{J=1}^{S_{d-1,p}}c_{J}\widehat{f}_{J}=0,$
$\displaystyle\sum_{J=S_{d-1,p}+1}^{N_{d,p}}c_{J}\partial_{x_{1}}b_{J}(({{\mathbf{x}}}_{T})_{1},\cdot)=$
$\displaystyle\sum_{J=S_{d-1,p}+1}^{N_{d,p}}c_{J}\widetilde{f}_{J-S_{d-1,p}}=0.$

This implies that $c_{J}=0$ for $J=1,\dots,N_{d,p}$, since
$\left\\{\widehat{f}_{s}\right\\}_{s=1,\dots,S_{d-1,p}}$ and
$\left\\{\widetilde{f}_{s}\right\\}_{s=1,\dots,S_{d-1,p-1}}$ are linearly
independent. ∎

A consequence of Proposition 5.4.1 is that the conditions in the definition of
${\mathbb{Q\\!T}}^{p}(T)$ are linearly independent:

$\dim\big{(}\mathbb{P}^{p}_{d}(T)\big{)}-\\#\\{{\boldsymbol{i}}\in\mathbb{N}^{d}\mid|{\boldsymbol{i}}|\leq
p-2\\}=\begin{pmatrix}p+d\\\ d\end{pmatrix}-\begin{pmatrix}p+d-2\\\
d\end{pmatrix}=N_{d,p}=\dim\big{(}{\mathbb{Q\\!T}}^{p}(T)\big{)}.$

In particular, we have

$\dim\big{(}{\mathbb{Q\\!T}}^{p}(T)\big{)}=N_{d,p}=\begin{cases}2&d=1\\\
2p+1&d=2\\\
(p+1)^{2}&d=3\end{cases}\quad={\mathcal{O}}_{p\to\infty}(p^{d-1}).$

Notice that, in the one dimensional case, if we vary the polynomial degree
$p$, then the dimension of the quasi-Trefftz space remains the same but the
space changes.

By comparing the dimension of ${\mathbb{Q\\!T}}^{p}(T)$ with the dimension of
$\mathbb{P}^{p}_{d}(T)$, we obtain

$\displaystyle\dim\big{(}{\mathbb{Q\\!T}}^{p}(T)\big{)}=$
$\displaystyle\frac{(p+d-2)!(2p+d-1)}{(d-1)!p!}$ $\displaystyle=$
$\displaystyle{\mathcal{O}}_{p\to\infty}(p^{d-1})\ll\dim(\mathbb{P}^{p}(T))=\binom{p+d}{d}={\mathcal{O}}_{p\to\infty}(p^{d}).$

We highlight that for large polynomial degrees $p$, the dimension of the
quasi-Trefftz space is much smaller than the dimension of the full polynomial
space of the same degree.

We have shown in Theorem 5.2.1 that smooth solutions of the diffusion-
advection-reaction equation with variable coefficients are approximated in
${\mathbb{Q\\!T}}^{p}(T)$ and in $\mathbb{P}_{d}^{p}(T)$ with the same
convergence rates with respect to the meshsize $h$. This is the main advantage
of Trefftz and quasi-Trefftz schemes: offering the same level of accuracy for
much fewer degrees of freedom.

###  5.5 The construction of the basis functions

The aim of this section is the construction of the basis functions of
${\mathbb{Q\\!T}}^{p}(T)$. We describe in detail the iterative algorithm for
the case $d=1$, $d=2$ and for the general $d$-dimensional case. Every step of
the iterations is based on the recursive formula (5.4), hence depends on the
Taylor coefficients at ${{\mathbf{x}}}_{T}$ of the PDE coefficients
${\bm{K}}$, ${\boldsymbol{\beta}}$ and $\sigma$. We also stress that each step
uses values already computed.

####  5.5.1 Algorithm for the case $d=1$

Consider the case $d=1$. As noticed before, in this particular case the
dimension of the quasi-Trefftz space ${\mathbb{Q\\!T}}^{p}(T)$ is equal to
$N_{1,p}=2$, which does not depend on $p$, even though the space does. We
recall that, in the one dimensional case, the space of the analytical
solutions of a homogeneous linear differential equation of the second-order
with continuous coefficients is a two-dimensional space; this is reflected by
the quasi-Trefftz space. Since we are in the scalar case, we denote by x the
only variable, by $\text{x}_{T}$ the point ${{\mathbf{x}}}_{T}\in T$ and by
$i$ the multi-indices ${\bm{i}}$. The two quasi-Trefftz basis functions
$b_{1}$ and $b_{2}$ are degree-$p$ polynomials in one variable such that
satisfy

$\displaystyle b_{1}(\text{x}_{T})=$ $\displaystyle\widehat{f}\text{ and
}\partial_{x}b_{1}(\text{x}_{T})=0,$ $\displaystyle b_{2}(\text{x}_{T})=$
$\displaystyle 0\text{ and }\partial_{x}b_{2}(\text{x}_{T})=\widetilde{f},$

where $\widehat{f}$ and $\widetilde{f}$ are two non zero constants that have
been chosen. For $J=1,2$ we have

$b_{J}(\text{x})=\sum_{k\in\mathbb{N},k\leq
p}a_{k}\left(\frac{\text{x}-\text{x}_{T}}{h_{T}}\right)^{k},$

and we want to compute the coefficients $a_{k}$. The coefficients $a_{0}$ and
$a_{1}$ are determined by the choice of $\widehat{f}$ and $\widetilde{f}$.
Formula (5.4) allows to compute $a_{i+2}$ from the coefficients $a_{j}$ with
$j\leq i+1$. We can perform a single loop to compute $a_{i+2}$ for increasing
values of $i$. This procedure is described in Algorithm 1.

Algorithm

Data: $p$, $\text{x}_{T}$, $h_{T}$, $D^{\ell}{\bm{K}}(\text{x}_{T})$,
$D^{\ell}{\boldsymbol{\beta}}(\text{x}_{T})$ for $\ell\leq p-1$,
$D^{\ell}\sigma(\text{x}_{T})$ for $\ell\leq p-2$.

Fix coefficients $a_{0}$, $a_{1}$, choosing polynomial bases $\widehat{f}$,
$\widetilde{f}$.

For each $J=1,2$, we construct $b_{J}$ as follows:

for  _$i=0$ to $p-2$_ do

$\begin{split}a_{i+2}=&\frac{{h_{T}^{i+2}}}{{\bm{K}}(\text{x}_{T})(i+2)!}\left(-\sum_{0<\ell\leq
i+1}\frac{(i+1)!}{\ell!}D^{\ell}\bm{K}(\text{x}_{T})\frac{a_{i+2-\ell}(i+2-\ell)}{h_{T}^{i+2-\ell}}\right.\\\
&\left.+\sum_{0\leq\ell\leq
i+1}\frac{(i+1)!}{\ell!}D^{\ell}\bm{\beta}(\text{x}_{T})\frac{a_{i+1-\ell}}{h_{T}^{i+1-\ell}}\right.\\\
&\left.+\sum_{0\leq\ell\leq
i}\frac{i!}{\ell!}D^{\ell}\sigma(\text{x}_{T})\frac{a_{i-\ell}}{h_{T}^{i-\ell}}\right),\end{split}$

end for

$\displaystyle b_{J}(x)=\sum_{k\in\mathbb{N},k\leq
p}a_{k}\bigg{(}\frac{x-\text{x}_{T}}{h_{T}}\bigg{)}^{k}.$

Algorithm 1 The algorithm for the construction of $b_{J}$ in the case $d=1$.

In order to show that if we vary the polynomial degree $p$ the quasi-Trefftz
space changes, even though the dimension remains the same, we consider the
following example: a homogeneous diffusion-advection-reaction equation in
$T=(0,2)$ with coefficients

${\bm{K}}(x)=x^{4}+1,\qquad{\boldsymbol{\beta}}(x)=\text{sin}(x)+2,\qquad\sigma(x)=1.$
(5.18)

In Figure 5.1 we plot the two quasi-Trefftz basis functions $b_{1}$ and
$b_{2}$ for polynomial degree $p=1,2,3,4$ using Algorithm 1.

(a) $p=1$

(b) $p=2$

(c) $p=3$

(d) $p=4$

Figure 5.1: The two basis functions $b_{1}$ and $b_{2}$ of the Quasi-Trefftz
space ${\mathbb{Q\\!T}}^{p}((0,2))$ for polynomial degree $p=1,2,3,4$ for the
example (5.18).

####  5.5.2 Algorithm for the case $d=2$

We now focus on the two dimensional case and describe the iterative algorithm
to compute the $N_{2,p}=2p+1$ quasi-Trefftz basis functions. For
$J=1,\dots,2p+1$, we need to compute the coefficients of the monomial
expansion of the basis function $b_{J}$, i.e, $a_{k_{1},k_{2}}$ with
$k_{1},k_{2}\in\mathbb{N}$, $k_{1}+k_{2}\leq p$. In Figure 5.2, the
coefficients $a_{k_{1},k_{2}}$ are represented by the dots in the plane of
indices $(k_{1},k_{2})\in\mathbb{N}^{2}$. Under the constraint that
$k_{1}+k_{2}\leq p$, these dots form a triangular shape in the plane. To
initialize the algorithm we have to choose two polynomial bases:

$\left\\{\widehat{f}_{s}\right\\}_{s=1,\dots,p+1}\text{ basis for
}\mathbb{P}^{p}(\mathbb{R})\quad\text{ and
}\quad\left\\{\widetilde{f}_{s}\right\\}_{s=1,\dots,p}\text{ basis for
}\mathbb{P}^{p-1}(\mathbb{R}),$

which we call “Cauchy data”. This choice fixes the coefficients $a_{0,k_{2}}$
with $0\leq k_{2}\leq p$ and $a_{1,k_{2}}$ with $0\leq k_{2}\leq p-1$,
represented by the shaded yellow area in the figure. Once these are known, all
the other coefficients are uniquely defined and can be computed using the
relation (5.4). From this recurrence relation we can compute
$a_{i_{1}+2,i_{2}}$ from the coefficients $a_{j_{1},j_{2}}$ with $j_{1}\leq
i_{1}+2$, $j_{2}\leq i_{2}+2$ and $j_{1}+j_{2}\leq i_{1}+i_{2}+2$, different
from $a_{i_{1}+2,i_{2}}$ itself. We have to apply (5.4) following a precise
ordering of the multi-indices $(i_{1},i_{2})$, since at every step of the
iterations we need to use values already computed. This suggests to proceed
“diagonally”, computing the values $a_{k_{1},k_{2}}$ for $k_{1}+k_{2}=r$
increasingly from $r=2$ to $r=p$. On each of these diagonals we compute the
values of $a_{k_{1},r-k_{2}}$ for increasing $k_{1}$. This is done performing
a double loop: the external loop moves away from the origin across diagonals
($\nearrow$) and the inner loop moves from the $k_{2}$ axis to the $k_{1}$
axis along diagonals ($\searrow$). The ordering of the coefficients
$a_{i_{1}+2,i_{2}}$ computed is shown by the arrows in Figure 5.3. This
procedure is described in Algorithm 2.

Figure 5.2: A graphical representation of the coefficients of a quasi-Trefftz
basis function $b_{J}$ in the two dimensional case (here $p=6$). The black
dots $\bullet$ in the $(k_{1},k_{2})$ plane represent the coefficients
$a_{k_{1},k_{2}}$ of the function $b_{J}$. The coefficients with
$k_{1}\in\\{0,1\\}$, which correspond to the nodes in the shaded yellow area,
are determined from the “Cauchy data”, the second and third sets of equations
in (5.17), which we assign. The shape of the stencil formed by the nodes
surrounded by the blue diamond
$\mathop{\raisebox{-1.06558pt}{$\leavevmode\hbox to9.35pt{\vbox
to9.35pt{\pgfpicture\makeatletter\hbox{\hskip-3.07497pt\lower 10.82495pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{}{}{{}}{{}}\pgfsys@setlinewidth{0.3875pt}\pgfsys@invoke{
}\pgfsys@roundjoin\pgfsys@invoke{ }{}{{}}{} {}{} {}{} {}{} {}{}
{}\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0.5,1,1}\pgfsys@color@cmyk@stroke{0.5}{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{1.6pt}\pgfsys@invoke{
}{}\pgfsys@moveto{7.74997pt}{11.62495pt}\pgfsys@lineto{3.87497pt}{15.49992pt}\pgfsys@lineto{7.74997pt}{19.37492pt}\pgfsys@lineto{11.62494pt}{15.49992pt}\pgfsys@lineto{7.74997pt}{11.62495pt}\pgfsys@closepath\pgfsys@stroke\pgfsys@invoke{
} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$}}$ represents the
first equation in (5.17) for $(i_{1},i_{2})=(2,2)$ which, given
$a_{j_{1},j_{2}}$ with $j_{1}\leq i_{1}+2$, $j_{2}\leq i_{2}+2$ and
$j_{1}+j_{2}\leq i_{1}+i_{2}+2$ (except for $a_{i_{1}+2,i_{2}}$), allows to
compute the coefficient $a_{i_{1}+2,i_{2}}$ associated with the node
surrounded by the red circle . Each coefficient in the non-shaded region is
computed with formula (5.4) in this way, using its corresponding stencil of
blue
diamonds.$p$$0$$a_{i_{1}+2,i_{2}}$$a_{i_{1},i_{2}+2}$$k_{1}$$k_{2}$Cauchydata$i_{1}$$i_{2}$$p$

$p$$0$Cauchydata$k_{1}$$k_{2}$

Figure 5.3: An illustration of the ordering of the coefficients of a quasi-
Trefftz basis function $b_{J}$ computed through Algorithm 2. All the
coefficients in the non-shaded region are computed with formula (5.4) in a
double loop: first across diagonals $\nearrow$, and then along diagonals
$\searrow$. The ordering is shown by the arrows
$\mathop{\raisebox{-1.06558pt}{$\leavevmode\hbox to8.83pt{\vbox
to8.83pt{\pgfpicture\makeatletter\hbox{\hskip 4.27496pt\lower-0.68062pt\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}\definecolor{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@rgb@stroke{0}{0}{0}\pgfsys@invoke{
}\pgfsys@color@rgb@fill{0}{0}{0}\pgfsys@invoke{
}\pgfsys@setlinewidth{0.4pt}\pgfsys@invoke{ }\nullfont\hbox
to0.0pt{\pgfsys@beginscope\pgfsys@invoke{
}{}{}{{}}{{}}\pgfsys@setlinewidth{0.3875pt}\pgfsys@invoke{
}\pgfsys@roundjoin\pgfsys@invoke{ }{}{{}}{}
{}{}\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@setlinewidth{0.8pt}\pgfsys@invoke{
}\color[rgb]{1,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,1}\pgfsys@color@cmyk@stroke{0}{1}{0}{0}\pgfsys@invoke{
}\pgfsys@color@cmyk@fill{0}{1}{0}{0}\pgfsys@invoke{
}\definecolor[named]{pgffillcolor}{rgb}{1,0,1}{}{}{}{}{{}}\pgfsys@moveto{-3.87497pt}{7.74995pt}\pgfsys@lineto{3.3093pt}{0.56569pt}\pgfsys@stroke\pgfsys@invoke{
}{{}{{}}{}{}{{}}{{{}}{{{}}{\pgfsys@beginscope\pgfsys@invoke{
}\pgfsys@transformcm{0.7071}{-0.7071}{0.7071}{0.7071}{3.59212pt}{0.28287pt}\pgfsys@invoke{
}\pgfsys@invoke{ \lxSVG@closescope }\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope}}{{}}}} \pgfsys@invoke{\lxSVG@closescope }\pgfsys@endscope
\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope{}{}{}\hss}\pgfsys@discardpath\pgfsys@invoke{\lxSVG@closescope
}\pgfsys@endscope\hss}}\lxSVG@closescope\endpgfpicture}}$}}$.

Algorithm

Data: $p$, ${{\mathbf{x}}}_{T}$, $h_{T}$,
$D^{\bm{\ell}}{\bm{K}}({{\mathbf{x}}}_{T})$,
$D^{\bm{\ell}}{\boldsymbol{\beta}}({{\mathbf{x}}}_{T})$ for
$\left|\bm{\ell}\right|\leq p-1$, $D^{\bm{\ell}}\sigma({{\mathbf{x}}}_{T})$
for $\left|\bm{\ell}\right|\leq p-2$.

Fix coefficients $a_{0,k_{2}}$, $a_{1,k_{2}}$, choosing polynomial bases
$\left\\{\widehat{f}_{J}\right\\}$, $\left\\{\widetilde{f}_{J}\right\\}$.

For each $J=1,\dots,2p+1$, we construct $b_{J}$ as follows:

for  _$r=2$ to $p$ (loop across diagonals $\nearrow$) _ do

for  _$i_{1}=0$ to $r-2$ (loop along diagonals $\searrow$) _ do

set $i_{2}=r-i_{1}-2$ and compute
$\begin{split}a_{{\bm{i}}+2{{\mathbf{e}}}_{1}}=&\frac{{h_{T}^{\left|{\bm{i}}\right|+2}}}{{\bm{K}}_{11}({{\mathbf{x}}}_{T})({\bm{i}}+2{{\mathbf{e}}}_{1})!}\Biggl{(}-\sum_{\begin{subarray}{c}j,m\in\\{1,2\\}\\\
\bm{\ell}\leq\bm{i}+\bm{e}_{j}\\\
(j,m,\bm{\ell})\neq(1,1,\bm{0})\end{subarray}}\frac{(\bm{i}+\bm{e}_{j})!}{\bm{\ell}!}\times\\\
&D^{\bm{\ell}}\bm{K}_{jm}({{\mathbf{x}}}_{T})\frac{a_{\bm{i}+\bm{e}_{j}-\bm{\ell}+\bm{e}_{m}}(i_{m}+(\bm{e}_{j})_{m}-\ell_{m}+1)}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|+2}}\\\
&+\sum_{\begin{subarray}{c}j=1,2\\\
\bm{\ell}\leq\bm{i}+\bm{e}_{j}\end{subarray}}\frac{(\bm{i}+\bm{e}_{j})!}{\bm{\ell}!}D^{\bm{\ell}}\bm{\beta}_{j}({{\mathbf{x}}}_{T})\frac{a_{\bm{i}+\bm{e}_{j}-\bm{\ell}}}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|+1}}\\\
&+\sum_{\bm{\ell}\leq\bm{i}}\frac{{\bm{i}}!}{\bm{\ell}!}D^{\bm{\ell}}\sigma({{\mathbf{x}}}_{T})\frac{a_{\bm{i}-\bm{\ell}}}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|}}\Biggr{)},\end{split}$

end for

end for

$\displaystyle b_{J}(x_{1},x_{2})=\sum_{0\leq k_{1}+k_{2}\leq
p}a_{k_{1},k_{2}}\bigg{(}\frac{x_{1}-({{\mathbf{x}}}_{T})_{1}}{h_{T}}\bigg{)}^{k_{1}}\bigg{(}\frac{x_{2}-({{\mathbf{x}}}_{T})_{2}}{h_{T}}\bigg{)}^{k_{2}}.$

Algorithm 2 The algorithm for the construction of $b_{J}$ in the case $d=2$.

####  5.5.3 Algorithm for the case $d>2$

We describe the extension of Algorithm 2 to the general $d$-dimensional case.
The structure is the same, the only difference is that we need to perform
three loops instead of two. In fact, for each value of
$r=\left|(i_{2},\dots,i_{d})\right|+i_{1}+2$ and of $i_{1}$ we have not only
one coefficient to compute but one for each
$(i_{2},\dots,i_{d})\in\mathbb{N}^{d-1}$ with
$\left|(i_{2},\dots,i_{d})\right|=r-i_{1}-2$. For this reason we need a
further inner loop over $(i_{2},\dots,i_{d})$. We note that each coefficient
of this loop can be computed independently of the others and we preserve the
requirement that at every step of the iterations all the values needed have
already been computed. In Figure 5.4 the coefficients $a_{{\boldsymbol{k}}}$
are represented by the dots forming a simplex in the space of indices
${\boldsymbol{k}}\in\mathbb{N}^{d}$ for the case $d=3$. The coefficients in
the two yellow triangles are assigned by the “Cauchy data”. The general
coefficient, highlighted by the red diamond, is computed with the relation
(5.4) using the coefficients corresponding to the black dots inside the blue
stencil. This procedure is described in Algorithm 3.

$k_{1}$$k_{3}$$k_{2}$$p$$p$$p$

Figure 5.4:  A graphical representation of the coefficients
$a_{{\boldsymbol{k}}}$ of a quasi-Trefftz basis function $b_{J}$ in the three
dimensional case (here $p=6$). The black dots  in the $(k_{1},k_{2},k_{3})$
space represent the coefficients $a_{k_{1},k_{2},k_{3}}$ of the function
$b_{J}$. The two yellow triangles in the planes $k_{1}=0$ and $k_{1}=1$
contain the coefficients whose values are determined from the choice of
$\\{\widehat{f}_{s}\\}_{s=1,\dots,S_{d-1,p}}$ and
$\\{\widetilde{f}_{s}\\}_{s=1,\dots,S_{d-1,p-1}}$. The coefficient
$a_{i_{1}+2,i_{2},i_{3}}$, represented by the red diamond
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\blacklozenge}$,
is computed with formula (5.4) for $(i_{1},i_{2},i_{3})=(2,1,1)$ using all the
coefficients inside the blue stencil. The coefficients whose dots are outside
the two yellow triangles are computed with formula (5.4) from the coefficients
inside their corresponding blue stencil. Algorithm 3 performs three loops: the
first loop through triangles parallel to the violet one (corresponding to
$r=6$), the second loop through triangles parallel to the yellow ones, and the
third one along the segments determined by the intersection between the two
triangles.

Algorithm

Data: $p$, ${{\mathbf{x}}}_{T}$, $h_{T}$,
$D^{\bm{\ell}}{\bm{K}}({{\mathbf{x}}}_{T})$,
$D^{\bm{\ell}}{\boldsymbol{\beta}}({{\mathbf{x}}}_{T})$ for
$\left|\bm{\ell}\right|\leq p-1$, $D^{\bm{\ell}}\sigma({{\mathbf{x}}}_{T})$
for $\left|\bm{\ell}\right|\leq p-2$.

Fix coefficients $a_{0,k_{2},\dots,k_{d}}$, $a_{1,k_{2},\dots,k_{d}}$,
choosing polynomial bases $\left\\{\widehat{f}_{s}\right\\}$,
$\left\\{\widetilde{f}_{s}\right\\}$.

For each $J=1,\ldots,N_{d,p}$, we construct $b_{J}$ as follows:

for  _$r=2$ to $p$ (loop across $\\{i_{1}+i_{2}+\dots+i_{d}=r-2\\}$
hyperplanes $\nearrow$)_ do

for  _$i_{1}=0$ to $r-2$ (loop across constant-$i_{1}$ hyperplanes
$\rightarrow$)_ do

for  _ $(i_{2},\dots,i_{d})$ with $\left|(i_{2},\dots,i_{d})\right|=r-i_{1}-2$
_ do

$\begin{split}a_{{\bm{i}}+2{{\mathbf{e}}}_{1}}=&\frac{{h_{T}^{\left|{\bm{i}}\right|+2}}}{{\bm{K}}_{11}({{\mathbf{x}}}_{T})({\bm{i}}+2{{\mathbf{e}}}_{1})!}\Biggl{(}-\sum_{\begin{subarray}{c}j,m=1,\dots,d\\\
\bm{\ell}\leq\bm{i}+\bm{e}_{j}\\\
(j,m,\bm{\ell})\neq(1,1,\bm{0})\end{subarray}}\frac{(\bm{i}+\bm{e}_{j})!}{\bm{\ell}!}\times\\\
&D^{\bm{\ell}}\bm{K}_{jm}({{\mathbf{x}}}_{T})\frac{a_{\bm{i}+\bm{e}_{j}-\bm{\ell}+\bm{e}_{m}}(i_{m}+(\bm{e}_{j})_{m}-\ell_{m}+1)}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|+2}}\\\
&+\sum_{\begin{subarray}{c}j=1,\dots,d\\\
\bm{\ell}\leq\bm{i}+\bm{e}_{j}\end{subarray}}\frac{(\bm{i}+\bm{e}_{j})!}{\bm{\ell}!}D^{\bm{\ell}}\bm{\beta}_{j}({{\mathbf{x}}}_{T})\frac{a_{\bm{i}+\bm{e}_{j}-\bm{\ell}}}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|+1}}\\\
&+\sum_{\bm{\ell}\leq\bm{i}}\frac{{\bm{i}}!}{\bm{\ell}!}D^{\bm{\ell}}\sigma({{\mathbf{x}}}_{T})\frac{a_{\bm{i}-\bm{\ell}}}{h_{T}^{\left|{\bm{i}}-\bm{\ell}\right|}}\Biggr{)},\end{split}$

end for

end for

end for

$\displaystyle
b_{J}({{\mathbf{x}}})=\sum_{{\boldsymbol{k}}\in\mathbb{N}^{d},|{\boldsymbol{k}}|\leq
p}a_{\boldsymbol{k}}\left(\frac{\mathbf{x}-\mathbf{x}_{T}}{h_{T}}\right)^{{\boldsymbol{k}}}$.

Algorithm 3 The algorithm for the construction of $b_{J}$ in the
$d$-dimensional case.

##  Chapter 6 Numerical experiments

In this chapter we present some numerical tests in two dimensions. We begin
with a simple example with polynomial solution in Section 6.1. In Section 6.2
we consider two test cases with smooth solutions in order to validate the
$h$-convergence of the quasi-Trefftz DG method and compare the three different
formulations SIPG, IIPG and NIPG introduced in Section 3.1 that result for a
different choice of the symmetrization parameter $\epsilon$. Moreover, in
Section 6.3 we compare the quasi-Trefftz space against the full polynomial one
in terms of $p$-convergence and of number of degrees of freedom. Finally, in
Section 6.4 we consider advection-dominated and reaction-dominated examples to
show the behaviour of the quasi-Trefftz DG method also in these cases.

The DG scheme (3.8) with as discrete space the global quasi-Trefftz space
(5.9) for the diffusion-advection-reaction equation has been implemented in
MATLAB. We remark that we consider the homogeneous diffusion-advection-
reaction equation, i.e., with $f=0$. In all the experiments we use the SIPG
method, i.e., the symmetrization parameter $\epsilon=-1$, except in the
example of Section 6.2 where the three formulations are compared.

All the experiments are performed on the unit square $\Omega=(0,1)^{2}$, which
is partitioned using conforming meshes of triangular elements generated by
DistMesh, see [29], and satisfying the assumptions of Section 2.1. For all the
experiments we will say, with an abuse of notation, the meshsize $h$ used,
though this is just the input of the mesh generator and usually the meshsize
as defined in 2.1.4 is a little larger than this input. We calculate the
“correct” $h$ and we use this value in the plots of the $h$-convergence in
Section 6.1 and in Section 6.2.

We compute all the volume integrals in the assembly of the matrix and of the
right-hand side using the following quadrature rule. First, we use the
Gauss–Legendre quadrature rule on the interval $[0,1]$, then we consider the
tensor product of the nodes in $[0,1]^{2}$ and then with the Duffy
transformation [1, p. 21] we map the points into the reference triangle, whose
vertices are $(0,0)$, $(1,0)$ and $(0,1)$. Finally, the quadrature points are
mapped onto the generic triangle $T\in{\mathcal{T}}_{h}$ using the classical
affine map, the reference element onto $T$. For computing the integrals over
edges we simply use the Gauss–Legendre quadrature rule combined with this
affine map. Since the Gauss–Legendre quadrature rule with $n$ points is exact
for polynomials of degree less than or equal to $2n-1$, we use $p+1$ points in
$[0,1]$, hence $(p+1)^{2}$ on each triangles, where $p$ is the polynomial
degree used. The linear system is solved using the backslash command.

To construct the quasi-Trefftz basis functions, we use the symbolic
computations only once to calculate the derivative of the PDE coefficients
that we need for Algorithm 2 and to evaluate them in the barycentre of each
element. Their values on the barycentres are stored and passed in input to the
quasi-Trefftz algorithm. To initialize Algorithm 2, we choose scaled monomials
as polynomial bases:

$\displaystyle\left\\{\widehat{f}_{s}(x_{2})=\bigg{(}\frac{x_{2}-({{\mathbf{x}}}_{T})_{2}}{h_{T}}\bigg{)}^{s-1}\right\\}_{s=1,\dots,p+1}$
$\displaystyle\text{ basis for }\mathbb{P}^{p}(\mathbb{R}),$
$\displaystyle\left\\{\widetilde{f}_{s}(x_{2})=\frac{1}{h_{T}}\bigg{(}\frac{x_{2}-({{\mathbf{x}}}_{T})_{2}}{h_{T}}\bigg{)}^{s-1}\right\\}_{s=1,\dots,p}$
$\displaystyle\text{ basis for }\mathbb{P}^{p-1}(\mathbb{R}).$

###  6.1 Polynomial solution

The aim of this simple example is to test the correctness of the quasi-Trefftz
DG method implemented in the particular case when the exact solution belongs
to the discrete space. Consider the problem (1.6)-(1.8) with $f=0$ and with
coefficients

${\bm{K}}(x_{1},x_{2})=Id,\qquad{\boldsymbol{\beta}}(x_{1},x_{2})=\bm{0},\qquad\sigma(x_{1},x_{2})=\frac{4}{x_{1}^{2}+x_{2}^{2}+1},$

where $Id$ denotes the $2\times 2$ identity matrix.

We impose only Dirichlet boundary condition $g_{D}$ on $\partial\Omega$ such
that the exact solution is given by $u(x_{1},x_{2})=x_{1}^{2}+x_{2}^{2}+1$.
Since $u$ is a polynomial of degree $2$ which satisfy the partial differential
equation, we have that $u\in{\mathbb{Q\\!T}}^{2}({\mathcal{T}}_{h})$ for the
definition of the quasi-Trefftz space with $p=2$. This follows also from
Theorem 5.2.1 which states that the Taylor polynomial of degree $2$ of $u$,
which is $u$ itself, belongs to ${\mathbb{Q\\!T}}^{2}(T)$ for all
$T\in{\mathcal{T}}_{h}$.

We use ${\mathbb{Q\\!T}}^{2}({\mathcal{T}}_{h})$ as discrete space and choose
the penalty parameter $\gamma=32$ and the symmetrization parameter
$\epsilon=-1$. In Figure 6.1 we show the quasi-Trefftz DG error against the
meshsize $h$. We consider the values $h=2^{-j}$ for $j=1,\dots,5$ and we
measure the error in the $L^{2}$-norm, in the $H^{1}$-norm and in the
$L^{\infty}$-norm. We observe that the error is not null but reaches almost
the machine precision. Even though from a theoretical point of view the error
should be zero, this is what we expect because of round-off errors.

Figure 6.1: The quasi-Trefftz DG error for the problem of Section 6.1 using
$p=2$. The error in the $L^{2}$-norm (red line), in the $H^{1}$-norm (yellow
line) and in the $L^{\infty}$-norm (blue line) is plotted against the meshsize
$h$ for $h=2^{-j}$ for $j=1,\dots,5$. Since the exact solution
$u\in{\mathbb{Q\\!T}}^{2}({\mathcal{T}}_{h})$ the error almost reaches the
machine precision.

###  6.2 h-convergence

In this section we presents two test cases with smooth solutions in order to
validate the optimal convergence rates with respect to the mesh size $h$ of
the quasi-Trefftz DG method that we have proved in Theorem 5.3.1.

As first example, we consider the problem (1.6)-(1.8) with $f=0$ and with
coefficients

${\bm{K}}(x_{1},x_{2})=\left[\begin{matrix}e^{x_{1}-x_{2}}&0\\\
0&e^{x_{1}-x_{2}}\end{matrix}\right],\qquad{\boldsymbol{\beta}}(x_{1},x_{2})=\bm{0},\qquad\sigma(x_{1},x_{2})=0.$
(6.1)

We assign only Dirichlet boundary condition $g_{D}$ on $\partial\Omega$ such
that the exact solution is $u(x_{1},x_{2})=e^{-x_{1}+x_{2}}$.

If the solution is smooth, such as in this case, the approximation bounds of
Theorem 5.3.1 for the $\left\|\cdot\right\|_{dar}$-norm of the quasi-Trefftz
DG error hold. Here we compute the $L^{2}$-norm and the $H^{1}$-norm of the
quasi-Trefftz DG error, which are both bounded above by the
$\left\|\cdot\right\|_{dar}$-norm. In Figure 6.2 we plot the error against the
meshsize $h$ for polynomial degree $p=1,2,3,4$.

(a) $p=1$

(b) $p=2$

(c) $p=3$

(d) $p=4$

Figure 6.2: $h$-convergence for the test case with coefficients (6.1). The
quasi-Trefftz DG error is plotted against the meshsize $h$ for the polynomial
degree $p=1,2,3,4$. The red lines represent the error measured in the
$L^{2}$-norm, the yellow lines in the $H^{1}$-norm and the blue lines in the
$L^{\infty}$-norm. The dots correspond to $h=2^{-1},\dots,2^{-5}$.

For each polynomial degree we consider the values $h=2^{-j}$ for $j=1,\dots,5$
and we measure the error in the $L^{2}$-norm (red line) in the $H^{1}$-norm
(yellow line) and in the $L^{\infty}$-norm (blue line). We choose
$\gamma=8p^{2}$ as penalty parameter and $\epsilon=-1$ as symmetrization
parameter. We observe rates of convergence of order $\mathcal{O}(h^{p})$
(dashed violet line) in the $H^{1}$-norm and of order $\mathcal{O}(h^{p+1})$
(dashed green line) in the $L^{2}$-norm and in the $L^{\infty}$-norm. The
rates correspond to the theoretical approximation bounds: rates of convergence
of order $\mathcal{O}(h^{p})$ for the $\left\|\cdot\right\|_{dar}$-norm of the
error. Hence, we have the optimal approximation properties that one achieve
with the full polynomial space.

We now consider a second experiment with a smooth solution in order to test
the quasi-Trefftz DG method in a more general case with also Neumann boundary
condition and with all the terms of diffusion, advection and reaction. As
second example, we consider the problem (1.6)-(1.8) with $f=0$ and with
coefficients

${\bm{K}}(x_{1},x_{2})=\left[\begin{matrix}x_{1}+x_{2}+1&0\\\
0&x_{1}+x_{2}+1\end{matrix}\right],\quad{\boldsymbol{\beta}}(x_{1},x_{2})=\left[\begin{matrix}1\\\
0\end{matrix}\right],\quad\sigma(x_{1},x_{2})=\frac{3}{x_{1}+x_{2}+1}.$ (6.2)

We assign Dirichlet boundary condition $g_{D}$ on the inflow
$\Gamma_{-}=\\{{{\mathbf{x}}}\in\ (0,1)^{2}\mid x_{1}=0\\}$ and Neumann
boundary condition $g_{N}$ on the outflow $\Gamma_{+}=\\{{{\mathbf{x}}}\in\
(0,1)^{2}\mid x_{1}=1\text{ or }x_{2}=0\text{ or }x_{2}=1\\}$ such that the
exact solution is given by $u(x_{1},x_{2})=\frac{1}{x_{1}+x_{2}+1}$. Computing
the quasi-Trefftz DG error, we obtain results similar to the ones shown in
Figure 6.2 for the previous case.

We now use this problem to compare the three formulations that arise varying
the value of the symmetrization parameter $\epsilon$. We recall that for
$\epsilon=-1$ we have the SIPG method, for $\epsilon=0$ we have the IIPG
method and for $\epsilon=+1$ we have the NIPG method, as described in Section
3.1. In Table 6.1 and in Table 6.2 we show the quasi-Trefftz DG error measured
in the $L^{2}$-norm and in the $H^{1}$-norm, respectively. We consider the
polynomial degrees $p=2,3,4$ and the penalty parameter $\gamma=8p^{2}$. For
each polynomial degree we compute the error for the three methods for
$h=2^{-3},2^{-4},2^{-5}$.

quasi-Trefftz DG error in $L^{2}$-norm
---
| $h$ | SIPG | rate | IIPG | rate | NIPG | rate
$p=2$, $\gamma=32$ | $2^{-3}$ | $1.445\times 10^{-5}$ | - | $3.159\times 10^{-5}$ | - | $4.990\times 10^{-5}$ | -
$2^{-4}$ | $1.704\times 10^{-6}$ | 3.085 | $6.488\times 10^{-6}$ | 2.284 | $1.105\times 10^{-5}$ | 2.175
$2^{-5}$ | $1.991\times 10^{-7}$ | 3.098 | $1.491\times 10^{-6}$ | 2.122 | $2.616\times 10^{-6}$ | 2.079
$p=3$, $\gamma=72$ | $2^{-3}$ | $8.496\times 10^{-7}$ | - | $7.530\times 10^{-7}$ | - | $7.183\times 10^{-7}$ | -
$2^{-4}$ | $4.954\times 10^{-8}$ | 4.100 | $4.276\times 10^{-8}$ | 4.138 | $4.044\times 10^{-8}$ | 4.151
$2^{-5}$ | $2.947\times 10^{-9}$ | 4.071 | $2.479\times 10^{-9}$ | 4.108 | $2.289\times 10^{-9}$ | 4.143
$p=4$, $\gamma=128$ | $2^{-3}$ | $5.753\times 10^{-8}$ | - | $5.929\times 10^{-8}$ | - | $6.578\times 10^{-8}$ | -
$2^{-4}$ | $1.860\times 10^{-9}$ | 4.951 | $2.104\times 10^{-9}$ | 4.817 | $2.681\times 10^{-9}$ | 4.617
$2^{-5}$ | $5.991\times 10^{-11}$ | 4.956 | $8.486\times 10^{-11}$ | 4.632 | $1.296\times 10^{-10}$ | 4.371
Table 6.1: Comparison of the three formulation SIPG, IIPG and NIPG in terms of
$L^{2}$ error for polynomial degree $p=2,3,4$ for the example (6.2).

All the three methods present optimal convergence rates of order
$\mathcal{O}(h^{p})$ in the $H^{1}$-norm, as shown in Table 6.2. Observing the
rates in the $L^{2}$-norm in Table 6.1, we see that the SIPG method has
convergence rates of order $\mathcal{O}(h^{p+1})$ while the IIPG and the NIPG
methods present rates of order $\mathcal{O}(h^{p+1})$ if the polynomial degree
is odd and of order between $\mathcal{O}(h^{p})$ and $\mathcal{O}(h^{p+1})$ if
the polynomial degree is even. This is also observed numerically in [31, p.
57] for a problem with only the diffusion term and in [31, Theorem 2.14] is
proved that the $L^{2}$ optimal error estimate holds for the SIPG method while
the others two have suboptimal rates. For this reason and also since in the
problems with only diffusion it preserves the symmetry of the original
problem, the SIPG formulation is preferred.

quasi-Trefftz DG error in $H^{1}$-norm
---
| $h$ | SIPG | rate | IIPG | rate | NIPG | rate
$p=2$, $\gamma=32$ | $2^{-3}$ | $1.145\times 10^{-3}$ | - | $1.113\times 10^{-3}$ | - | $1.096\times 10^{-3}$ | -
$2^{-4}$ | $2.785\times 10^{-4}$ | 2.039 | $2.722\times 10^{-4}$ | 2.032 | $2.687\times 10^{-4}$ | 2.028
$2^{-5}$ | $6.640\times 10^{-5}$ | 2.069 | $6.528\times 10^{-5}$ | 2.060 | $6.467\times 10^{-5}$ | 2.055
$p=3$, $\gamma=72$ | $2^{-3}$ | $6.721\times 10^{-5}$ | - | $6.587\times 10^{-5}$ | - | $6.501\times 10^{-5}$ | -
$2^{-4}$ | $7.992\times 10^{-6}$ | 3.072 | $7.882\times 10^{-6}$ | 3.063 | $7.813\times 10^{-6}$ | 3.057
$2^{-5}$ | $9.229\times 10^{-7}$ | 3.114 | $9.157\times 10^{-7}$ | 3.106 | $9.111\times 10^{-7}$ | 3.100
$p=4$, $\gamma=128$ | $2^{-3}$ | $4.770\times 10^{-6}$ | - | $4.652\times 10^{-6}$ | - | $4.562\times 10^{-6}$ | -
$2^{-4}$ | $2.978\times 10^{-7}$ | 4.002 | $2.909\times 10^{-7}$ | 3.9995 | $2.854\times 10^{-7}$ | 3.999
$2^{-5}$ | $1.819\times 10^{-8}$ | 4.033 | $1.777\times 10^{-8}$ | 4.033 | $1.743\times 10^{-8}$ | 4.034
Table 6.2: Comparison of the three formulation SIPG, IIPG and NIPG in terms of
$H^{1}$ error for polynomial degree $p=2,3,4$ for the example (6.2).

###  6.3 p-convergence and number of DOFs

We now compare the quasi-Trefftz space
${\mathbb{Q\\!T}}^{p}({\mathcal{T}}_{h})$ with the full polynomial space
$\mathbb{P}_{2}^{p}({\mathcal{T}}_{h})$ in terms of $p$-convergence and of
number of degrees of freedom (DOFs). We consider the first example in Section
6.2 with coefficients (6.1), we fix the meshsize $h=2^{-4}$ and choose the
parameters $\gamma=8p^{2}$ and $\epsilon=-1$.

(a) Error in norm $L^{2}$

(b) Error in norm $H^{1}$

Figure 6.3: $p$-convergence for the problem with coefficients (6.1). The red
lines represent the DG error for the ${\mathbb{Q\\!T}}^{p}({\mathcal{T}}_{h})$
while the dashed blue lines the DG error for
$\mathbb{P}^{p}_{2}({\mathcal{T}}_{h})$. The error is measured in the
$L^{2}$-norm (a) and in the $H^{1}$-norm (b) and the meshsize is $h=2^{-4}$.

In Figure 6.3 we plot the error, in the $L^{2}$-norm (a) and in the
$H^{1}$-norm (b), obtained using the quasi-Trefftz space (red line) and the
full polynomial space (dashed blue line) for polynomial degree $p=1,2,3,4$.
These plots show that the quasi-Trefftz space, as the full polynomial space,
leads to exponential convergence in terms of polynomial degree $p$. However,
we observe that the full polynomial space presents a slightly smaller error.
The advantages of the quasi-Trefftz space are shown in Figure 6.4, where the
error is compared in terms of global number of degrees of freedom. We clearly
see that, for comparable number of degrees of freedom, with the quasi-Trefftz
space we can reach a higher accuracy, especially for high polynomial degrees
$p$.

(a) Error in norm $L^{2}$

(b) Error in norm $H^{1}$

Figure 6.4: Comparison between ${\mathbb{Q\\!T}}^{p}({\mathcal{T}}_{h})$ and
$\mathbb{P}_{2}^{p}({\mathcal{T}}_{h})$ in terms of degrees of freedom for the
problem with coefficients (6.1) using $h=2^{-4}$. The dots correspond to
$p=1,2,3,4$ and the errors are measured in the $L^{2}$-norm (a) and in the
$H^{1}$-norm (b). These plots show that with the quasi-Trefftz space we have
the same level of accuracy for much fewer degrees of freedom.

###  6.4 Advection-dominated and reaction-dominated problems

In this section we consider advection-dominated and reaction-dominated
problems in order to illustrate the approximate solutions obtained with the
quasi-Trefftz DG method in these particular cases. We recall that the exact
solution of this kind of problems can present boundary and internal layers,
which are narrow regions where the solution and its gradient change rapidly.
For this reason the numerical treatment of such problems can be difficult. It
is known that, if we use the standard Finite Element Method, the approximate
solution, in general, presents spurious oscillations and can be completely
different from the exact one.

First, we consider the problem (1.6)-(1.8) with $f=0$ and with coefficients

${\bm{K}}(x_{1},x_{2})=\nu
Id,\qquad{\boldsymbol{\beta}}(x_{1},x_{2})=\left[\begin{matrix}x_{2}+1\\\
-x_{1}+2\end{matrix}\right],\qquad\sigma(x_{1},x_{2})=0,$ (6.3)

for different values of the diffusion parameter $\nu=10^{-j}$ for
$j=1,\dots,4$. We impose the Dirichlet boundary condition on the inflow
$\Gamma_{-}=\\{{{\mathbf{x}}}\in\ \partial\Omega\mid x_{1}=0\text{ or
}x_{2}=0\\}$:

$g_{D}(x_{1},x_{2})=\begin{cases}1\quad\text{ on $\\{x_{1}=0,0\leq
x_{2}<1\\}$,}\\\ 1\quad\text{ on $\\{x_{2}=0,0\leq
x_{1}\leq\frac{1}{3}\\}$,}\\\ 0\quad\text{ on
$\\{x_{2}=0,\frac{1}{3}<x_{1}<1\\}$,}\\\ \end{cases}$ (6.4)

and the Neumann boundary condition $g_{N}=0$ on the outflow
$\Gamma_{+}=\\{{{\mathbf{x}}}\in\ \partial\Omega\mid x_{1}=1\text{ or
}x_{2}=1\\}$. We note that, as the value of the diffusion parameter $\nu$
decreases, the advection term ${\boldsymbol{\beta}}$ becomes dominant and for
$0<\nu\ll 1$ the exact solution presents an internal layer.

In Figure 6.5 we plot the quasi-Trefftz DG approximate solutions obtained
using the meshsize $h=2^{-4}$ and the polynomial degree $p=3$.

(a) $\nu=10^{-1}$, $\gamma=10$

(b) $\nu=10^{-2}$, $\gamma=1$

(c) $\nu=10^{-3}$, $\gamma=10^{-1}$

(d) $\nu=10^{-4}$, $\gamma=10^{-2}$

Figure 6.5: Approximate solutions of the problem with coefficients (6.3),
obtained using the meshsize $h=2^{-4}$ and the quasi-Trefftz space with
polynomial degree $p=3$.

We fix the symmetrization parameter $\epsilon=-1$ and for each value of
$\nu=10^{-1},\dots,10^{-4}$ we choose the penalty parameter
$\gamma=10,1,10^{-1},10^{-2}$, respectively. These choices of the penalty
parameter are made taking into consideration the two following aspects. On the
one hand, we cannot choose a too small $\gamma$ because from the theory of
Section 4.2 we know that, in order to have the discrete coercivity, the
penalty parameter needs to be greater than $\gamma_{\epsilon}$, defined in
(4.4), which in this case is equal to $3\nu C_{tr}$, with $C_{tr}$ defined in
(2.8). On the other hand, if the penalty parameter is too large we are
imposing the continuity across the mesh facets and we obtain similar results
to ones of the standard Galerkin method, which, as previously stated, it does
not perform well in the advection-dominated regime due to the spurious
oscillations. We notice that we do not have spurious oscillations spread over
all the domain but the flat part of the solution is well approximated and the
discontinuity at the boundary and the internal layer are well captured with
small oscillations.

We now consider the same example with coefficients (6.3) but with Dirichlet
boundary condition (6.4) on the inflow $\Gamma_{-}$ and homogeneous Dirichlet
condition on the outflow $\Gamma_{+}$. This is done in order to assess the
behaviour of the quasi-Trefftz DG method in a test case not covered by the
theory, since the assumption (1.14) is not valid. In this experiment in the
advection-dominated regime the exact solution presents a boundary layer on the
edge $\\{x_{2}=0\\}$ and also an internal layer.

(a) $\nu=10^{-1}$, $\gamma=10$

(b) $\nu=10^{-2}$, $\gamma=10^{-1}$

(c) $\nu=10^{-3}$, $\gamma=10^{-2}$

(d) $\nu=10^{-4}$, $\gamma=10^{-3}$

Figure 6.6: Approximate solutions of the problem with coefficients (6.3) with
only Dirichlet boundary condition, obtained using the meshsize $h=2^{-4}$ and
the quasi-Trefftz space with polynomial degree $p=3$.

In Figure 6.6 we plot the quasi-Trefftz DG approximate solutions obtained
using the meshsize $h=2^{-4}$ and the polynomial degree $p=3$. We fix the
symmetrization parameter $\epsilon=-1$ and for each value of
$\nu=10^{-1},\dots,10^{-4}$ we choose the penalty parameter
$\gamma=10,10^{-1},10^{-2},10^{-3}$, respectively. The internal layer is well
captured with only very small oscillations. The boundary layer, instead, is
completely missed in the advection-dominated regime but this is what we expect
when using DG methods due to the weak enforcement of boundary conditions [4].
In fact, we obtain very similar results if we use the DG scheme with the
classical full polynomial space.

We now consider the problem (1.6)-(1.8) with $f=0$ and with coefficients

${\bm{K}}(x_{1},x_{2})=\nu
Id,\qquad{\boldsymbol{\beta}}(x_{1},x_{2})=\bm{0},\qquad\sigma(x_{1},x_{2})=x_{1}+x_{2}+1,$
(6.5)

for different values of the diffusion parameter $\nu=10^{-j}$ for
$j=1,\dots,4$. We assign only Dirichlet boundary condition $g_{D}=1$ on
$\partial\Omega$. In this case, as the value of the diffusion parameter $\nu$
decreases, the reaction term $\sigma$ becomes dominant and for $0<\nu\ll 1$
the exact solution presents boundary layers on $\partial\Omega$. In Figure 6.6
we plot the quasi-Trefftz DG approximate solutions obtained using the meshsize
$h=2^{-4}$ and the polynomial degree $p=3$. We fix the symmetrization
parameter $\epsilon=-1$ and for each value of $\nu=10^{-1},\dots,10^{-4}$ we
choose the penalty parameter $\gamma=10,1,10^{-1},10^{-2}$, respectively. As
in the examples before, the flat part of the solution is well approximated and
the boundary layers are captured with small oscillations.

(a) $\nu=10^{-1}$, $\gamma=10$

(b) $\nu=10^{-2}$, $\gamma=1$

(c) $\nu=10^{-3}$, $\gamma=10^{-1}$

(d) $\nu=10^{-4}$, $\gamma=10^{-2}$

Figure 6.7: Approximate solutions of the problem with coefficients (6.5),
obtained using the meshsize $h=2^{-4}$ and the quasi-Trefftz space with
polynomial degree $p=3$.

##  Chapter 7 Conclusions and future developments

We have introduced a quasi-Trefftz discontinuous Galerkin method for
approximating the homogeneous diffusion-advection-reaction equation with
varying coefficients. The DG variational formulation, its well-posedness and
the quasi-optimality in a mesh-dependent norm have been analysed. We have
presented the polynomial quasi-Trefftz space for a general homogeneous linear
partial differential equation with sufficiently smooth coefficients and proved
that it reaches the same approximation orders with respect to the mesh size as
the full polynomial space but with fewer degrees of freedom. For the
homogeneous diffusion-advection-reaction equation, we have proved optimal
$h$-convergence of the quasi-Trefftz DG scheme and described an algorithm for
the construction of the quasi-Trefftz basis functions in dimension
$d\in\mathbb{N}$, $d\geq 1$. We have presented some numerical experiments in
two dimensions that validate the accuracy of the method.

We now consider some possible improvements of the method and future research
directions.

To better assess the computational performance of the method, the method
should be implemented in the three-dimensional case and in a high-performance
open-source software such as NGSolve, where the Trefftz method and quasi-
Trefftz method for the wave equation have already been implemented.

The algorithm for the construction of the quasi-Trefftz basis functions is
initialized assigning two polynomial bases. Those conditions could be
optimized, through a more accurate study of the properties of the basis
functions, in order to improve conditioning, computing time and accuracy.

The construction of non-polynomials quasi-Trefftz functions might also be
analysed; they could be useful, for example, to efficiently approximate
solutions that present boundary layers, such as in the case of advection-
dominated or reaction-dominated problems.

We have proved optimal convergence rate with respect to the mesh size
(h-convergence) and a challenging extension is the analysis of the
approximation properties for increasing polynomial degrees (p-convergence).
This is a difficult task since in the context of Trefftz schemes has been
achieved for time-harmonic equations (Helmholtz [17] and Maxwell [18]) but not
yet for the wave equation and neither in the context of the quasi-Trefftz
methods.

The analysis of the DG scheme considered includes the source term while the
quasi-Trefftz space is designed for the homogenous equation. A significant
improvement of the method is the extension to the inhomogeneous advection-
diffusion-reaction equation. We need to construct a particular approximate
local solution, considering a variant of the algorithm presented, and then
apply the quasi-Trefftz DG method in an adequate way to solve for the
difference.

Another interesting extension is the combination of the method presented and
the quasi-Trefftz DG method already implemented for the wave equation with
piecewise-smooth coefficients [23], in order to approximate PDEs whose nature
changes in the domain, such as the Euler-Tricomi equation
($\partial^{2}_{x}u+x\partial^{2}_{y}u=0$), useful in the study of transonic
flow. Other similar equations are used in plasma physics, for example,
$-\Delta u+\sigma u=0$, with $\sigma$ that changes sign, which is a Helmholtz
equation in a part of the domain and a diffusion-reaction equation in the
other [25].

## References

  * [1] P.. Antonietti, C. Marcati, I. Mazzieri and A. Quarteroni  “High order discontinuous Galerkin methods on simplicial elements for the elastodynamics equation”  In _Numerical Algorithms_ 71 Springer, 2016, pp. 181–206
  * [2] D.. Arnold  “An interior penalty finite element method with discontinuous elements”  In _SIAM journal on numerical analysis_ 19.4 SIAM, 1982, pp. 742–760
  * [3] D.. Arnold, F. Brezzi, B. Cockburn and L.. Marini  “Unified analysis of discontinuous Galerkin methods for elliptic problems”  In _SIAM journal on numerical analysis_ 39.5 SIAM, 2002, pp. 1749–1779
  * [4] B. Ayuso and L.. Marini  “Discontinuous Galerkin methods for advection-diffusion-reaction problems”  In _SIAM Journal on Numerical Analysis_ 47.2 SIAM, 2009, pp. 1391–1420
  * [5] S.. Brenner  “The mathematical theory of finite element methods”  Springer, 2008
  * [6] H. Brezis  “Functional analysis, Sobolev spaces and partial differential equations”  Springer, 2011
  * [7] F. Brezzi, L.. Marini and A. Russo  “On the choice of a stabilizing subgrid for convection–diffusion problems”  In _Computer methods in applied mechanics and engineering_ 194.2-5 Elsevier, 2005, pp. 127–148
  * [8] F. Brezzi, L.. Marini and E. Süli  “Discontinuous Galerkin methods for first-order hyperbolic problems”  In _Mathematical models and methods in applied sciences_ 14.12 World Scientific, 2004, pp. 1893–1903
  * [9] J.. Callahan  “Advanced calculus: a geometric view”  Springer, 2010
  * [10] B. Cockburn  “Discontinuous Galerkin methods”  In _ZAMM-Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik: Applied Mathematics and Mechanics_ 83.11 Wiley Online Library, 2003, pp. 731–754
  * [11] D.. Di Pietro and A. Ern  “Mathematical aspects of discontinuous Galerkin methods”  Springer Science & Business Media, 2011
  * [12] L.. Evans  “Partial differential equations (graduate studies in mathematics, vol. 19)”  In _Instructor_ 67 American Mathematics Society (1998) Text# 2 (Required):, 2009
  * [13] S. Gómez and A. Moiola  “A Space-Time Trefftz Discontinuous Galerkin Method for the Linear Schrödinger Equation”  In _SIAM Journal on Numerical Analysis_ 60.2 SIAM, 2022, pp. 688–714
  * [14] S. Gómez and A. Moiola  “A space-time DG method for the Schrödinger equation with variable potential”  In _arXiv preprint arXiv:2306.05780_ , 2023
  * [15] S. Gómez, A. Moiola, I. Perugia and P. Stocker  “On polynomial Trefftz spaces for the linear time-dependent Schrödinger equation”  In _Applied Mathematics Letters_ 146, 2023
  * [16] F. Heimann, C. Lehrenfeld, P. Stocker and H. Wahl  “Unfitted Trefftz discontinuous Galerkin methods for elliptic boundary value problems”  In _arXiv preprint arXiv:2212.12236_ , 2022
  * [17] R. Hiptmair, A. Moiola and I. Perugia  “Plane wave discontinuous Galerkin methods for the 2D Helmholtz equation: analysis of the p-version”  In _SIAM Journal on Numerical Analysis_ 49.1 SIAM, 2011, pp. 264–284
  * [18] R. Hiptmair, A. Moiola and I. Perugia  “Error analysis of Trefftz-discontinuous Galerkin methods for the time-harmonic Maxwell equations”  In _Mathematics of Computation_ 82.281, 2013, pp. 247–268
  * [19] R. Hiptmair, A. Moiola and I. Perugia  “A survey of Trefftz methods for the Helmholtz equation”  In _Building bridges: connections and challenges in modern approaches to numerical partial differential equations_ Springer, 2016, pp. 237–279
  * [20] P. Houston, C. Schwab and E. Süli  “Discontinuous hp-finite element methods for advection-diffusion-reaction problems”  In _SIAM Journal on Numerical Analysis_ 39.6 SIAM, 2002, pp. 2133–2163
  * [21] L.. Imbert-Gerard and G. Sylvand  “Three types of quasi-Trefftz functions for the 3D convected Helmholtz equation: construction and approximation properties”  In _arXiv preprint arXiv:2201.12993_ , 2022
  * [22] L.. Imbert-Gérard and B. Després  “A generalized plane-wave numerical method for smooth nonconstant coefficients”  In _IMA Journal of Numerical Analysis_ 34.3 Oxford University Press, 2014, pp. 1072–1103
  * [23] L.. Imbert-Gérard, A. Moiola and P. Stocker  “A space–time quasi-Trefftz DG method for the wave equation with piecewise-smooth coefficients”  In _Mathematics of Computation_ 92.341, 2023, pp. 1211–1249
  * [24] L.. Imbert-Gérard and G. Sylvand  “A roadmap for Generalized Plane Waves and their interpolation properties”  In _Numerische Mathematik_ 149.1 Springer, 2021, pp. 87–137
  * [25] L.M. Imbert-Gérard  “Analyse mathématique et numérique de problèmes d’ondes apparaissant dans les plasmas magnétiques”, 2013
  * [26] C. Lehrenfeld and P. Stocker  “Embedded Trefftz discontinuous Galerkin methods”  In _International Journal for Numerical Methods in Engineering_ Wiley Online Library, 2023
  * [27] A. Moiola  “Trefftz-discontinuous Galerkin methods for time-harmonic wave problems”, 2011
  * [28] A. Moiola and I. Perugia  “A space–time Trefftz discontinuous Galerkin method for the acoustic wave equation in first-order formulation”  In _Numerische Mathematik_ 138.2 Springer, 2018, pp. 389–435
  * [29] P.. Persson and G. Strang  “A simple mesh generator in MATLAB”  In _SIAM review_ 46.2 SIAM, 2004, pp. 329–345
  * [30] I. Perugia, J. Schöberl, P. Stocker and C. Wintersteiger  “Tent pitching and Trefftz-DG method for the acoustic wave equation”  In _Computers & Mathematics with Applications_ 79.10 Elsevier, 2020, pp. 2987–3000
  * [31] B. Rivière  “Discontinuous Galerkin methods for solving elliptic and parabolic equations: theory and implementation”  SIAM, 2008
  * [32] E. Trefftz  “Ein Gegenstück zum Ritzschen Verfahren. Verhandl. des 2. Intern. Kongresses f”  In _Technische Mechanik. Zürich_ , 1926