Embedded Trefftz discontinuous Galerkin methods

In Trefftz discontinuous Galerkin methods a partial differential equation is discretized using discontinuous shape functions that are chosen to be elementwise in the kernel of the corresponding differential operator. We propose a new variant, the embedded Trefftz discontinuous Galerkin method, which is the Galerkin projection of an underlying discontinuous Galerkin method onto a subspace of Trefftz-type. The subspace can be described in a very general way and to obtain it no Trefftz functions have to be calculated explicitly, instead the corresponding embedding operator is constructed. In the simplest cases the method recovers established Trefftz discontinuous Galerkin methods. But the approach allows to conveniently extend to general cases, including inhomogeneous sources and non-constant coefficient differential operators. We introduce the method, discuss implementational aspects and explore its potential on a set of standard PDE problems. Compared to standard discontinuous Galerkin methods we observe a severe reduction of the globally coupled unknowns in all considered cases, reducing the corresponding computing time significantly. Moreover, for the Helmholtz problem we even observe an improved accuracy similar to Trefftz discontinuous Galerkin methods based on plane waves.


Introduction
In this manuscript we propose a novel numerical method closely related to Trefftz discontinuous Galerkin methods. The main idea of Trefftz methods, originating from [55], is to choose optimal discretization spaces that provide the same approximation quality as comparable discrete spaces with a significant reduction of the number of degrees of freedom (ndofs). For an overview on Trefftz methods see [26,34,52]. Polynomial Trefftz functions have been obtained for several linear partial differential operators with constant coefficients such as Laplace equation [42,43], acoustic wave equation [37,38,46], heat equation [53], plate vibration and beam vibration equation [1,35], and time-dependent Maxwell's equation [27]. Efforts to generate Trefftz polynomials in a general case have been undertaken, see [21,25]. The idea of finding Trefftz polynomials via Taylor series, used in [36], was extended in [57] to construct Trefftz-like polynomials in the case of non-linear elliptic PDEs possibly with smooth coefficients. The method requires a Taylor expansion of the coefficients of the PDE to construct these polynomials via a recursive procedure for the coefficients of a Taylor polynomial.
Discontinuous Galerkin (DG) schemes can be easily combined with Trefftz functions as the basis construction of elements decouples. Trefftz-DG schemes for Laplace equations have been analyzed in [20,32,33]. Different time-dependent problems have been discussed recently, such as wave problems in one space dimension [28,29,51], the acoustic wave equation [3,46,50], elastoacoustics [5], a class of Friedrichs systems coming from linear transport [6,47], for time-dependent Maxwell's equation see [12,27], and the linear Schrödinger equation [16]. Very popular applications of Trefftz methods are wave propagation problems in frequency domain. There, no polynomial Trefftz space exists and plane wave functions are used instead. A DG scheme for Helmholtz equation has been presented and analyzed in [7,15,17,41,44], for more information on Trefftz methods for Helmholtz equation see the recent survey [19] and the references therein. Plane waves have also been used for linear elasticity [45] and time-harmonic Maxwell's equation [14,18,44].
Apart from combining Trefftz spaces with DG schemes, harmonic polynomials and plane waves have also been combined with virtual element methods, see [39,40,49].
In the case of a differential operator with smooth coefficients the local solutions in a Trefftz space are not sufficient to provide high-order approximation. This can be circumvented by weakening the requirements of the Trefftz space, providing again a sufficiently large basis. In [24] a quasi-Trefftz DG method is introduced and analyzed for the acoustic wave equation with smooth coefficients. Plane waves have been generalized to work with a DG method for the Helmholtz equation with smooth coefficient in [23].
Trefftz methods usually are applied to homogeneous problems (problems with no volume source term). To treat inhomogeneous problems with a Trefftz approach, one needs to construct a particular solution. Then it is possible to apply Trefftz methods to solve for the difference. In [56] the Poisson problem is treated in two separate ways, once a fundamental solutions are used as trial functions, while for the other approach radial basis functions are used to express the inhomogeneous part. Another approach for Helmholtz and time-harmonic Maxwell's equation is presented in [22], where a series of inhomogeneous local problems are solved using spectral elements to obtain local particular solutions.
1.1. Main contributions and outline of the paper. The embedded Trefftz method proposed here side-steps the explicit construction of a Trefftz space, instead we present a simple method of constructing a corresponding embedding. In its simplest form, the method proposed here can be seen as a convenient way to set up the linear system of a Trefftz DG discretization by means of a Galerkin projection of a standard DG method onto its Trefftz subspace. Instead of implementing Trefftz functions explicitly, the characterization of the Trefftz function space as the kernel of an associated differential operator is exploited to construct an embedding of a Trefftz subspace into the DG space in a very generic way. We denote this embedding as the Trefftz embedding. The construction requires only element-local operations of small matrices and is embarrassingly parallel. To compute the discrete kernel numerically several established methods exists, e.g. one can use QR factorizations, singular, or eigenvalue decompositions. We denote this approach as embedded Trefftz DG method.
When the embedded Trefftz DG approach is used to implement an existing Trefftz DG method one trades the convenient setup for slightly larger costs in the setup of the linear system. However, the strongest feature of the embedded Trefftz DG approach is the following: The generic way the linear systems are set up allows to apply the concept of Trefftz methods beyond their previous limitations. Two of these limitations that can be conveniently exceeded are: • In many interesting applications a standard Trefftz space is either not polynomial or unfeasible to set up at all limiting the application of Trefftz DG methods for these cases. • If a standard Trefftz DG subspace exists, still only homogeneous equations (equations with no volume source terms) are easily dealt with. Inhomogeneous equations are difficult to deal with for standard Trefftz DG methods. In cases where a polynomial Trefftz space is not embedded in the discretization space or when a Trefftz basis is unknown, the proposed embedded Trefftz DG method can still be applied. This is done by slightly weakening the properties of the Trefftz space we are looking to embed. With this we are able to treat for example the Helmholtz equation and PDEs with piecewise smooth coefficients, e.g. the acoustic wave equation.
The embedded Trefftz DG method offers a more generic way of obtaining local particular solutions exploiting the underlying DG discretization. Computing element-wise particular solutions in the DG space can be done with little computational effort and then be used to homogenize the DG system.
The paper continuous as follows: In Section 2 we present the method, recovering traditional Trefftz DG methods and its extensions to the case of weaker Trefftz spaces and inhomogeneous equations. In Section 3 we address implementational aspects and turn to applications of the approach to a set of PDE problems in Section 4. We present numerical examples for Laplace equation, Poisson equation, acoustic wave equation with piecewise constant and also with smooth coefficient, Helmholtz equation, and a linear transport equation. In Section 5 we compare the Trefftz DG method to another popular acceleration technique for DG methods, the Hybrid DG method. The method was implemented using NGSolve and NGSTrefftz 1 . Conclusion and outlook are presented in Section 6.

The method
We introduce the method in several steps. After some preliminaries, we first introduce the basic concepts of the method in Section 2.2. Here, we will assume that a suitable Trefftz space exists, and is embedded in the discrete DG space V p (T h ). Here, we understand as the Trefftz space the subspace of functions that locally fulfill the PDE pointwise. This case only applies to a small set of problems where (polynomial) Trefftz spaces exist, like the Laplace equation and the acoustic wave equation with piece-wise constant wavespeed. We then turn to cases where a Trefftz space is not embedded in the discrete space. In this case we consider a weak Trefftz space instead. This is presented in Section 2.3. How inhomogeneous equations can easily be dealt with is explained in Section 2.4 resulting in a method that allows to treat a large class of PDE problems in a uniform way.
2.1. Preliminaries. We consider a linear partial differential operator L and a corresponding boundary (and possibly initial) value problem supplemented by suitable boundary conditions and paired with a suitable discontinuous Galerkin formulation The DG formulation is set on a mesh T h of the domain Ω, and we assume that the space V p (T h ) is the product of local spaces V p (K) for K ∈ T h where the index p indicates a polynomial degree, e.g. V p (K) = P p (K) the space of polynomials up to degree p. We note that space-time DG formulations are allowed in this setting, so that d ∈ N is either the dimension of the space or space-time domain.
2.2. Embedded Trefftz method. At first, in this subsection, we make the following simplifying assumptions: Firstly, we assume f = 0. Second, we assume that the differential operator has the form L = d l=1 α l D β l l for α l ∈ R and β l ∈ N and that all mesh elements are straight. Valid for a non-constant field α. This assumption ensures that the subspace of V p (K) that we define next is sufficiently large to allow for reasonable approximations of the solution.
We define the Trefftz space The Trefftz version of (2) is its Galerkin projection to T p (T h ): In contrast to previous works on Trefftz methods, we do not construct the space T p from scratch based on accordingly defined basis function, but rather aim to build an embedding operator for the Trefftz space into the underlying piecewise polynomial space V p (K). We call this the embedded Trefftz DG method.
Construction of the embedding operator. Let {φ i } be the set of basis functions of V p (T h ), N = dim(V p (T h )) and G : ., N the canonical unit vectors in R n , so that G(e i ) = φ i , we define the following matrices and vector for i, j = 1, . . . , N where ·, · 0,h = K∈T h ·, · K is the element-wise L 2 -inner product. We are interested in the kernel of L (in an element-wise and pointwise sense) in V p (T h ) as this is the part where the Trefftz DG method operates. We note that and hence are looking for a basis of ker(W). As W is block diagonal, with blocks corresponding to the elements K ∈ T h we will construct the basis element by element which implies that the following calculations can be done in parallel. For a fixed K ∈ T h let W K ∈ R N K ×N K be the corresponding block with N K = dim(V p (K)). For the dimension of the kernel we have M K := dim(ker(L)) = N K − L K with L K = dim(range(L)) where (for the case of differential operator L in the assumed form) L K is typically known. Now, we can determine the kernel of W K and collect a set of orthogonal basis vectors in a matrix T K ∈ R N K ×M K so that ker(W K ) = T K · R M K . The kernel matrix T K can be determined numerically, e.g. by a QR decomposition or a singular value decomposition, we discuss this in more depth in Section 3.4.
Composing the individual element matrices T K into a block matrix T ∈ R N ×M with M = dim(ker(W)) = T ∈T h M T we easily obtain the characterization of the kernel as ker(W) = T·R M . Further note that there holds T T T = I M ×M .
We define the Galerkin isomorphism between R M and T p (T h ) as G T : R M → T p (T h ), x → G(Tx) and denote a discrete operator corresponding to T as the Trefftz embedding T : . To setup the linear systems corresponding to (4) we first assemble A, l and T and arrive at: The resulting linear system is well controlled in terms of its conditioning.
Lemma 1 (Conditioning of the embedded Trefftz method). The spectral relative condition number of the embedded Trefftz DG method is bounded by the corresponding condition of the DG method, Proof. By construction of T all its column vectors are orthogonal.
Remark 2 (Matrix representing the kernel of L). There are other possibilities to characterize the kernel of L in V p (K). Here, we chose the specific form as it most obviously displays that u ∈ ker W is equivalent to LG(u) = 0 pointwise as there has to hold LG(u), LG(u) 0,h = 0. In Section 2.3 we will replace W with a proper generalization.
Remark 3 (Degrees of freedom). One major advantage of the Trefftz method is the reduction of the number of (globally coupled) degrees of freedom (ndofs) without harming the approximation too much. Let us first discuss the reduction on an example. Assume L is a second order operator, say −∆, the unknown field is scalar and the mesh is triangular. Then N K = (p + 1)(p + 2)/2 and −∆V p (K) = V p−2 (K) and hence L K = (p − 1)p/2 yields M K = 2p + 1. This holds for every element and similar calculations can be made for different differential operators and vectorial problems. In general, the Trefftz method decreases ndofs from O(p d ) to O(p d−1 ) and is hence especially attractive for higher order methods. In the lowest order case p = 0 or also p = 1 (depending on the differential operator) one often has T p (T h ) = V p (T h ). A similar reduction of ndofs from O(p d ) to O(p d−1 ) in a discontinuous Galerkin setting is also achieved in methods that associate dofs primarily to facets and rely in some way on static condensation such as Hybrid DG methods [8,9] or Hybrid High order (HHO) methods [11]. We discuss the comparison between these approaches in more detail in Section 5.
Remark 4 (Uniqueness of the Trefftz embedding). We notice that while the embedding T : is unique the matrix T and hence, the basis for T p (T h ) resulting from the previous procedure is not.
Remark 5 (Volume integrals). One advantage of Trefftz DG methods that is sometimes advertised is that volume integrals that involve Lu (possibly after partial integration) can be removed during assembly as Lu = 0 on the Trefftz space. This is also possible when using the embedded Trefftz method, see also Remark 10.

Weak Trefftz embedding.
In the previous section we restricted the differential operator to have a special form, here we aim to lift some of the restrictions, specifically we want to treat non-constant coefficients and mixed differential orders also including id. In these situations the Trefftz space as defined in (3) will typically much too small, i.e. the requirement Lv = 0 on the polynomial spaces may restrict the discrete solution too much leading to a version of locking.
To circumvent this problem, we weaken our condition in the Trefftz space. We introduce a projection Π that is yet to be defined and define the weak Trefftz space and the embedded weak Trefftz DG method : This replaces the pointwise condition Lv = 0 with something that is potentially weaker. The projection Π will be designed to ensure that the resulting weak Trefftz space is reasonably large to allow for proper approximations of the PDE solution.
A projection that has proven fruitful is given by The multi-index q ∈ N d on V q denotes the polynomial space of maximal order q = (q 1 , . . . , q d ) in each variable (x 1 , . . . , x d ).
(If all entries of the multi-index are equal, we will continue to write V q .) Let us now consider the differential operator L = j∈N d α j D j with α j ∈ C(T h ). Then we choose q i equal to the largest appearing differential order, i.e. q i = p − max α j ≡0 j i , i = 1, . . . , d.
With these definitions, the analogue to W from Section 2.2, i.e. the matrix that can be used to numerically extract the weak Trefftz space, becomes for j = 1, ..., N and i = 1, ...,Ñ = dim V q (T h ) and {ψ i } a basis for V q (T h ) (W) ij = Lφ j , ψ i . (11) In Section 3.1 we present a slightly more convenient version (with the same kernel) that we also used in the numerical examples in this manuscript.
Remark 6 (Consistency with Section 2.2). We recover the Trefftz method from Section 2.2 if we choose V q = LV p . For instance for L = −∆ we can choose V q = ∆V p = V p−2 . From the constant coefficient case treated in Section 2.2, it becomes clear that smaller values in q would make the space larger than the Trefftz space, resulting in a sub-optimal reduction of the space.
Remark 7 (Quasi-Trefftz). In [24] weak Trefftz spaces have been used to treat the acoustic wave equation with smooth coefficients, we recall the quasi-Trefftz spaces used there in (30). We recover this quasi-Trefftz method if we replace Π in (9b) with a Taylor polynomial expansion of order p − 1 in the element center.
2.4. Inhomogeneous PDEs. The Trefftz methods discussed in Section 2.3 and Section 2.2 as well as those known in the literature are designed for the homogeneous problem Lu = 0 only and hence are not able to solve for problems with non-zero source term. A simple adjustment however allows to deal also with the inhomogeneous case for the embedded Trefftz DG method.
In this section we treat equations of the form in Ω for f ∈ L 2 (Ω), supplemented by suitable boundary conditions for which a suitable DG discretizations is given in the form where the linear form g(·) corresponds to the weak imposition of boundary conditions. On the continuous level of the PDE we can introduce an affine shift and decompose the solution to Lu = f with u = u 0 + u f . Here u f is a particular solution to Lu f = f and u 0 (uniquely) solves Lu 0 = 0 (with boundary data depending on u f ). As the embedded Trefftz (or weak Trefftz) method is based on an underlying DG discretization we are able to construct a reasonable particular solution in the DG space V p (T h ) and apply the same homogenization strategy.
We then write . u h,f is a (non-unique) particular solution to ΠLu h,f = Πf that we can compute in an element-local fashion, see Section 3.5. For the existence of a particular solution we require that ΠL : This essentially depends on the choice of L, p and q. It always holds true in the setting of Section 2.2.
In the remainder we assume that a particular solution exists.
For u h,f a particular solution, we are looking for a solution After homogenization this means that we are looking for This translates to the solution of the linear system Remark 8 (Degenerate cases). Let us briefly discuss two extreme cases of the Trefftz embedding. If L = id, then T p (T h ) = {0} and the discrete solution of Lu = u = f will effectively be computed when determining the -in this case unique -particular solution u h,f . The other extreme case is Assuming a coercive problem and a DG formulation that provides good discretization properties, on the embedded (weak) Trefftz DG subspace we can still apply all essential analysis tools to obtain a Céa-type result: Lemma 9. Let u ∈ V (Ω) be a weak solution to the PDE problem under consideration and u h ∈ T p (T h ) + u h,f be the solution to (13) and · h be a suitable norm on V p (T h ) ∪ V (Ω). Assume that problem (2) is well-posed, specifically a h and are continuous with respect to where we used coercivity and continuity and that Adding a triangle inequality and taking the infimum yields the result.
The crucial question is of course on the approximation quality of the (affinely shifted) weak Trefftz space. We leave that open for future research.

Implementational aspects
Let us consider several implementational aspects in this separate section. First, in Section 3.1 we discuss how to set up a version of W -the matrix that characterizes the (weak) Trefftz subspace -in a more convenient way than what directly follows from Section 2.3. Next, we discuss how the dimension of the (weak) Trefftz subspace can be obtained in an implementation. Afterwards, we present the algorithmic structure of the embedded Trefftz DG method in comparison to a standard DG method. We discuss how to numerically compute the kernel of the matrix W and the particular solutions in Sections 3.4 and 3.5. In Section 3.6 we discuss the computational complexity of the embedded Trefftz DG method.
3.1. Characterization and implementation of the (weak) Trefftz DG space. Involving a new space V q (T h ) in the implementation and assembling of the matrix W K as used in (11) may occur somewhat unneccessarily cumbersome. In the following, we propose a different but equivalent formulation that is more convenient and is also applied in the numerical examples later. Instead of using V q (T h ) we again use the finite element space V p (T h ), and find a differential operatorL withL : The operatorL is a modification of L which only keeps the highest order derivatives in each direction and has only constant coefficients. In the simplest cases, i.e. in the setting of Section 2.2 we simply chooseL = L. Otherwise, in the general setting of Section 2.3, for example for the Helmholtz equation with variable wavespeed, with the operator This way, we can re-define the matrix W from (11) for the general case as Note that this matrix W has dimension N × N as in the original definition (5b) in the setting of Section 2.2.

3.2.
Dimension of the (weak) Trefftz subspace. On each element K ∈ T h the dimension of the ker ΠL, M K , can be computed as N K − dim(range ΠL) also for the (weak) Trefftz spaces. In most cases this expression could easily be precomputed and be given to the algorithm for the numerical kernel extraction. Alternative -and this is a bit more convenient as it requires less user input -we can also try to read that information from the matrix directly. Next, we explain how this can be achieved and notice that this approach has been used in all numerical examples below. When numerically computing the kernel of a matrix, we use a QR, singular value or eigenvalue decomposition. In all these factorizations the diagonal of one of these factors (e.g. the diagonal matrix of the SVD or the triangular matrix of the QR decomposition) would have M K zeros assuming exact arithmetics. Due to inexact computer arithmetics this will not be the case exactly and hence we use a truncation parameter ε > 0 to determine which values are considered as (numerical) zeros. A numerical study on the choice of the threshold for the Laplace problem is done in Section 4.2.
3.3. The algorithm. In Algorithm 1 we show how the embedded Trefftz DG method is implemented. We display the method in pseudo-code and a python code example in NGSolve, using the operatorL as described in Section 3.1. Let us briefly comment on the algorithm.
1. First of all, we notice that the computation and the application of the Trefftz embedding is separated from the setup of the linear system of the DG method. On the one hand this obviously leaves room for (efficiency) improvements as an integration into the DG assembly could allow to reduce the need for storing the DG matrices (and vectors in the homogeneous case) in the first place. On the other hand it also facilitates the realization of the embedded Trefftz DG method in a DG software framework without the need to interfere with optimized routines such as the assembly. We focus on the second approach in the presentation of the implementation in Algorithm 1. Both approaches are implemented in NGSolve+NGSTrefftz and can be found in the software documentation, see [54].
2. Further, we observe that the setup of the matrix T is done element-by-element and can hence be carried out in an embarrassingly parallel manner. The resulting matrix is block diagonal such that the setup of T T AT and the r.h.s. vector can also be done efficiently.
3. Assuming a corresponding DG method already exists for a user, he only needs to additionally specify the operators L andL,the truncation parameter ε and the r.h.s. to obtain the corresponding embedded Trefftz DG method. All of these inputs are canonical in most cases rendering the approach and its implementation quite easily accessible.
Require: Basis functions {φi}i, DG formulation (a h , l), operators L,L, truncation parameter ε, r.h.s. f 1: function dg matrix 2: Algorithm 1: Embedded Trefftz DG algorithm in the version discussed in Section 3. Pseudo-code (left) and running python code (using NGSolve and NGSTrefftz, right). The parts used by a standard DG scheme are colored black, the parts using the embedded Trefftz DG method, as in Section 2.3, are colored in teal. Code needed to treat the inhomogeneous problem, presented in Section 2.4, is colored purple.

3.4.
Computing the kernel matrix T K . In this section we briefly discuss how to numerically extract the kernel matrix The elementwise matrix W K is of sizeÑ K × N K , whereÑ K depends on the approach: If the matrix W K stems from either (5b) or (16), then W K is a square matrix andÑ K = N K . If (11) is used, theñ N K = dim V q (K). We discuss two possibilities: One based on a QR decomposition and one a singular value decomposition.
QR decomposition. The most obvious option seems to be a QR decomposition of W T K , s.t. W T K = Q K · R K for an orthogonal matrix Q K , of size N K ×N K , and an upper triangular matrix R K , sized N K ×Ñ K , where we assume an ordering such that ( Singular value decomposition. Applying a singular value decomposition (SVD) to W K yields W K = U K Σ K V T K for two orthogonal matrices U K , V K and a diagonal matrix Σ K with (Σ K ) jj = 0 for j > L K . Then one easily sees that the last M K columns of V K span the kernel of W K and we choose this submatrix as 3.5. Computing particular DG solutions. To compute a (local) particular solution needed to solve inhomogeneous PDEs, see Section 2.4, on an element K ∈ T h . We assemble (w K ) i = (f, G(e i )) = (f,Lφ i ) and define (u f ) K = W † K w K . Here W † K denotes the pseudoinverse of the matrix W K , which can be obtained using the QR decomposition or SVD of the matrix W K which may have already been computed when numerically computing the kernel of W K , cf. Section 3.4.
3.6. Algorithmic complexity. The computation of the kernel of ΠL or a particular solution scales well in the mesh size and can easily be parallelized. However, the scaling in the local unknowns per element is severe. For a discretization of order p the unknowns on each element scale like O(p d ). Hence, the costs for element-wise QR or SV decompositions (including the setup of the adjusted matrices and vectors and neglecting parallelization effects) are O(h −d p 3d ). Let us compare these costs with the solver costs of standard DG and classical Trefftz DG methods. For a standard DG method we have O(h −d p d ) unknowns globally and assuming a solver complexity The proposed method has a bad asymptotic complexity if h = const and p → ∞. In this case the computational costs will soon be dominated by the QR or SVD operations in the Trefftz DG setup. However, in the opposite case, p = const and h → 0 asymptotically the costs of the embedded Trefftz setup are only relevant for an optimal solver complexity α = 1 and are negligible otherwise. In the applications below we observe that at least moderately high (fixed) polynomial degrees can be easily used with the embedded Trefftz DG method without dominating the computational costs by the computation of the embedding.

Applications
We consider a set of different PDE problems and aim to compare the embedded Trefftz DG method to the underlying DG method (with full polynomial basis) and to a standard Trefftz basis, whenever available. We provide the scripts used to obtain the numerical results presented in this section 2 and interactive code examples to try out the implementation online without any prerequisites 3 . In

Preliminaries and Notation.
Although there exist many DG schemes for each of the problems considered here, in order to focus on the comparison of the different bases, we will stick to one DG scheme per equation.
Further, we concentrated on simplicial meshes here, but the approach carries directly over to hexahedral and even polygonal meshes immediately. If not explicitly mentioned we do not exploit the possibility to remove volume integrals, cf. Remark 5.
In the following numerical experiments we will always use sparse direct solvers, i.e. umfpack [10] for non-symmetric and a sparsecholesky solver implemented in NGSolve for symmetric linear systems.
In all examples, for the embedded Trefftz DG method we used the local SVD to compute the numerical kernel and particular soluations.
For the purpose of brevity in the plots below we label the standard DG method with DGand the corresponding Trefftz method with T p .
For the description of the DG schemes in the following we introduce some standard DG notation. We denote by F h the set of facets and distinguish F int h , the set of interior facets, from F bnd h , the set of boundary facets. Let K and K be two neighboring elements sharing a facet F ∈ F int h . On F the functions u K and u K denote the two limits of a discrete function from the different sides of the element interfaces. n K and n K are the unit outer normals to K and K . We define and consider a symmetric IP-DG discretization, cf. [2], given by As interior penalty parameter we choose α = 4. We apply our method as described in Section 2.2 with L = −∆. The (embedded) Trefftz DG space, as in (3), is now the space of harmonic polynomials. For the numerical example we set the boundary condition g such that the exact solution is given by and consider nested simplicial unstructured meshes created by refining a coarse simplicial unstructured mesh of initial mesh size h ≈ 0.5.  In Figure 1 we observe exponential convergence in terms of p and ndofs. In [20] exponential convergence for harmonic polynomials in terms of ndofs is seen to be superior to standard polynomials, this can be seen in the center of Figure 1. We measure the error in the DG-norm given by In Figure 1 we also compare the condition number of the different system matrices. As expected from Lemma 1, the conditioning of the embedded Trefftz method is bounded by the conditioning of the full system, and even outperforms it. The very basic implementation of the harmonic polynomials that we used here appears very ill-conditioned. We show plots of the runtime in Figures 2 and 3 for the two and three dimensional example, respectively. In Figure 2 it can be seen that the embedded Trefftz DG method benefits greatly from parallelization. The runtime is broken into parts of assembling the linear system and solving the linear system in Figures 2 and 3. For the embedded Trefftz DG method we plot also the time spent finding the local kernels of the operator using SVD, as described in Section 3.4. Note that in the case of the embedded Trefftz DG method the solving step includes the matrix multiplications needed in (7), therefore the time spent is not equal to that of the standard Trefftz method. In 2d performing the SVD sequentially is time consuming, as shown in Figure 2, a QR decomposition could improve the performance. However, in 3d, see Figure 3, the time spent on the SVD is completely negligible compared to the solver costs.   It is possible to automatically find the dimension of the Trefftz space when computing the null space of the local operator. To account for numerical errors in the computations of the singular values we choose a threshold of ε = 10 −7 , identifying values below as zero singular values. In Figure 2 we plot the largest singular value that still needs to be identified as a zero singular value, as well as the smallest non-zero singular value. Note that prior knowledge of the dimension of the Trefftz space, which we do not exploit here in the implementation, would eliminate this problem completely, cf. Section 3.2.
In [32,33] convergence rates in h for a mixed formulation over harmonic polynomials are shown. The results in Figure 3 show that we recover the expected convergence rate of u − u h L 2 (Ω) = O(h min(m,p)+1 ) for u ∈ H m (Ω).

Poisson equation. Now, we consider the Poisson equation with Dirichlet boundary conditions
We can use the symmetric IP-DG discretization given in (17) with the new right hand side For the numerical example we set the right hand side and the boundary conditions such that the exact solution is given by To apply the embedded Trefftz method we use the approach described in Section 2.4 to find a particular solution and homogenize the system.  Results in terms of accuracy and computing time are shown in Figure 4. The embedded Trefftz DG method matches the convergence rates of the DG scheme using the standard polynomial space. We consider the runtime for fixed mesh sizes and varying polynomial degree from p = 1, . . . , 7. As discussed in Section 3.6 we expect worse performance for sparse meshes and large polyomial degree, due to the increased cost of computing the local kernel of the operator. We can observe this on the mesh of the unit cube with mesh size h = 0.5, where the embedded Trefftz method does not show significant improvement of the runtime. The runtime on the fine mesh with h = 0.125 is considerably improved even for rather large values of p.

Poisson equation with varying coefficient. We consider the Poisson equation with varying coefficients and Dirichlet boundary conditions, given by
We can use the symmetric IP-DG discretization given in (17) and right hand side (20), by inserting M at the appropriate locations. For the numerical example we fix the right hand side f and the boundary conditions such that the exact solution is given by u = sin(x) sin(y) with M = 1 + x 0 0 1 + y , and Ω = (0, 1) 2 .
To apply the embedded Trefftz method we apply again use the approach detailed in Section 2.4 to homogenize the system. Due to the varying coefficient we cannot find a standard Trefftz space, thus we follow the approach details in Section 2.3, embedding a weak Trefftz space. To construct the Trefftz embedding we apply Eqs. (10) and (11), resulting in We compare different choices for the polynomial degree of the space V q , with q = p−1, p−2, p−3.
Note that the considerations in Section 2.3 imply that the optimal order is given by q = p − 2 since ∆ : V p → V p−2 .  The results presented in Figure 5 show the expected behavior. Testing with an increased number of test functions, using q = p−1, leads to an embedding of a smaller subspace than the actual weak Trefftz space. This results in a loss of the good approximation properties, as we see in Figure 5 on the left. Using a smaller test space V q with q = p − 3, is less efficient since the approximation space holds more degrees of freedom, as shown in Figure 5 on the right. While the error for the choices q = p − 2 and q = p − 3 shows slight differences, most likely due to different conditioning, the behavior is the same, as predicted by Lemma 9. Specifically, the different computation of the particular solution does not affect the approximation. Recall that the particular solutio satisfies ΠLu h,f = Πf , i.e. Lu h,f , v = f, v ∀v ∈ V q , as described in Section 2.4. Hence, for the choice of q = p − 3, the approximation of the particular solution might be worse than for the choice q = p − 2, however, the loss in accuracy is made up for by the larger Trefftz space used after homogenization.

4.5.
Acoustic wave equation with piecewise constant wave speed. Next, we consider the first order wave equation in a space-time setting, given by where (σ, v) are the unknowns -acoustic speed and pressure, (σ 0 , v 0 ) are the initial conditions, g D is Dirichlet boundary data, and T is the final time. In this section we will assume that the wavespeed c is piecewise constant.
To apply our framework we take The local space-time Trefftz space was introduced and analyzed in [46], and is given by We consider the space-time DG-scheme used in [4,24,46,50]: with and For the numerical example we set boundary and initial conditions such that the exact solution with wavespeed c = 1 is given by The penalization parameters in (27) are chosen as α = β = 0.5, as in [50]. The considered spacetime meshes are made up from a simplicial unstructured mesh of the spatial domain and a tensor product mesh in time. The height of the time slabs is chosen approximately as the mesh size of the spatial domain. In Fig. 6 we observe that both Trefftz method maintain the accuracy of the underlying DG method while yielding almost an order of magnitude speed up in the computation time on the finest considered level. Also on the finest level the solver costs seem to dominate so that the embedded Trefftz DG method and the Trefftz method perform equally well.
4.6. Acoustic wave equation in inhomogeneous media. Instead of the piecewise constant case from the previous section, we now consider (24) with wavespeed smoothly varying in space c = c(x). The construction of a basis for (25) with smooth wavespeed is non-viable. A space with similar properties, for which a basis can be constructed, has been introduced in [24] -the quasi-Trefftz space for an element K ∈ T h is given by where we use the notation i = (i x , i t ) = (i x1 , . . . , i xn , i t ) ∈ N n+1 0 for integer non-negative multiindices and D i f := ∂ ix 1 x1 · · · ∂ ix n xn ∂ it t f for the derivatives. The polynomials in the space are constructed such that the Taylor polynomial of their image under the wave operator vanishes at the element center (x K , t K ) up to order p − 2.
The conditions of the quasi-Trefftz space are formulated in a way to allow for recursive construction of the basis functions. It is possible to write this space in the way of (9b), see Remark 7. However, this is not the most practical way to implement the embedding. Thus, we proceed with  the construction proposed in Eqs. (10) and (11). While an explicit construction of the weak Trefftz space is unfeasible, the embedded Trefftz method is still able to project on such a space. We construct the weak Trefftz embedding as described in Sections 2.3 and 3.3, using the following condition for the weak Trefftz space (27) can be dropped when using the embedded Trefftz DG method, which we also did in the numerical examples in this section, cf. Remark 10.
We consider the exact solution given by with κ = 2.5 and wavespeed c(x, y) = x+y+1 on the space-time domain Ω×(0, 1) with Ω = (0, 1) 2 .  In Fig. 7 we observe a similar performance as for the homogeneous case. When looking carefully we can now see a difference between Quasi-Trefftz and embedded Trefftz DG method. While in the homogeneous case the embedded Trefftz space coincides with the Trefftz space, in the inhomogeneous case the two spaces do not coincide, compare (9b) and (30). However, the asymptotics seem to be unaffected and all methods perform equally well in terms of accuracy on a given mesh. The timings are again in agreement with the experience from the homogeneous case.
Remark 10. The volume term in (27) can be omitted for the case of homogeneous media when using Trefftz or embedded Trefftz DG methods. In the inhomogeneous case, in [24] the quasi-Trefftz methods also require the volume terms for stability, however, even in the case of smooth coefficients, these terms can still be omitted when using embedded Trefftz DG method as described in Section 4.6. Note that in [4,24] there appear additional Galerkin-least squares terms needed for the analysis, however, already in [24] it was observed that they do not seem to play a significant role in the numerical examples, which is why they are neglected here. We consider the DG-scheme used in [7,15,17,41,44] with (bi)linear forms with the choices α = β = δ = 0.5 used in [7]. A Trefftz space for the Helmholtz equation in two dimensions is given by the (non-polynomial) space of plane wave functions The combination of DG (bi)linear forms and this trial and test space is known as plane wave DG method. In the numerical experiment we consider the exact solution given by where H (1) 0 is the zero-th order Hankel function of the first kind. Results are shown in Fig. 8 for ω = 1. Surprisingly, the embedded Trefftz DG solution is extremely close to the plane wave DG method although the one space is piecewise polynomial and the other one is not. While all methods exhibit the same convergence rate, both Trefftz methods yield the significantly smaller error compared to the DG method. This can certainly not be attributed to different approximation spaces as the embedded Trefftz DG space is a subspace of the DG space. Hence, the results suggest that the embedded Trefftz method has better stability properties than the standard DG method, probably comparable to that of the plane wave DG method, cf. also Remark 11. In Fig. 8 on the right, we consider p-convergence on two different meshes. On the finer mesh, with h = 2 −3 , the plane waves fail to converge for p = 5. We have implemented the simplest form of the plane wave basis functions, given by the functions in (33), which are notoriously haunted by ill-conditioning. While more stable constructions for the plane wave basis exist, see for example [44,Sec. 3.4.1], it is interesting to note that the embedded Trefftz space shows improved conditioning without any additional effort.
Remark 11 (Approximation of plane wave functions with the embedded Trefftz DG space). The remarkable proximity of both Trefftz DG method in the previous experiment may suggest that the embedded Trefftz space approximates the plane wave DG space. We investigate this in more detail for a one-dimensional example. On [0, 1] we consider the embedded Trefftz space for p = 1, ..., 5. The embedded Trefftz space is only two-dimensional as in 1D dim V p −dim V q = 2 (independent of p). On this space we approximate the two plane wave functions sin(ωx) and cos(ωx) for ω = 2π. The result is shown in Fig. 9 (top row). For comparison we approximate sin(ωx) and cos(ωx)  on the (much) larger space V p (T h ) (bottom row). We observe that starting from p = 2 the approximation quality of the embedded Trefftz space is close to the that of the DG space and for p = 5 the plane waves are extremely well resolved rendering the embedded Trefftz space in close proximity to the plane wave space.

Linear transport equation.
In this section we finally consider an example that is typically not related to Trefftz method: A scalar linear transport problem, the advection equation. It reads as for a given velocity field b which we assume to be divergence-free. As underlying DG discretization we choose the standard Upwind DG formulation which reads as Here, we used the upwind notationû(x) = lim t→0 + u(x − bt). We do not assume b to be piecewise constant so that there will not be a suitable polynomial Trefftz space in general. To L = b · ∇ we hence choose the weak Trefftz space with V q with q = p − 1.The local Trefftz space on an element K will hence have the dimension M = #P p (K) − #P p−1 (K) which in 1D is M = 1, in 2D is M = p + 1 and in 3D is (p + 1)(p + 2)/2. In comparison to scalar second order problems we hence only have approximately half the degrees of freedoms in the resulting weak Trefftz space. In Fig. 10 we display a set of weak Trefftz basis functions that is obtained in 2D for a non-constant flow field and p = 4. The basis functions are approximately -but not exactly -constant along the flow trajectories.
For the numerical study we choose Ω = (0, 1) 3 and b = (− sin(x 2 ), cos(x 1 ), x 1 ) T and choose the r.h.s. data f and u D so that u = sin(x 1 ) sin(x 2 ) sin(x 3 ) is the exact solution. Starting from a coarse simplicial and unstructured mesh with mesh size h ≈ 0.5 we apply successive uniform refinements to compare the standard DG and the weak embedded Trefftz DG method. In the left half of Fig. 11 we observe the convergence behavior of the embedded Trefftz DG method compared to the standard Upwind DG method for p = 3, 4, 5. Both methods converge with optimal rate and we observe that there is only a marginal difference between the results. Moreover, a careful look at the numbers reveals that the solution of the embedded Trefftz DG method has a slightly smaller error on few occasions (see for instance the zoomed region in the left plot of Fig. 11). This may appear surprising at first glance as the approximation space is (by construction) smaller (T p (T h ) ⊂ V p (T h ) ), i.e. the approximation quality has not been improved. This suggests that the stability has been slightly improved in these cases. On the right half of Fig. 11 we compare the runtime of the two approaches for p = 3, 4, 5. We observe that on sufficiently fine meshes the costs associated to the two methods separate so that on the finest mesh level the computational costs differ already by more than an order of magnitude. Moreover, we observe that the embedded Trefftz DG method for p = 5 is still much cheaper than the standard DG method for p = 4 and only slightly more expensive than the standard DG method for p = 3, i.e. we can see the step from DG to embedded Trefftz DG as a way to obtain also two orders of accuracy more for the same computation time, at least in the given example.
Although linear transport problems are not in the typical class of problems that are associated with Trefftz methods, the obtained results are very encouraging. Compared to the second order equations considered before the reduction of degrees of freedom is even higher and the gain in the computation time is remarkable.

Comparison to Hybrid DG
We saw in the previous examples that the embedded Trefftz DG approach has two applications: First, it can be seen as an acceleration technique to reduce the costs of DG methods when solving linear systems. Second, the restriction of the DG space to a suitable subset can also have an effect on the stability of the method, cf. the Helmholtz problem in Section 4.7. For the former aspect alternative approaches exists. One technique that became very popular in the last decade is the class of Hybrid DG methods, see [8,9,13]. In the remainder of this section we want to shed some light on the comparison between the two approaches. To this end, we will first recap the structure of HDG method and make a rough conceptual comparison especially in terms of asymptotic complexity in Section 5.1 and afterwards compare the methods on the scalar examples from Section 4.3 and Section 4.8 and close the section with a rough conclusion.

5.1.
Introduction. In Hybrid DG methods additional unknowns are introduced on the element interfaces in order to allow for a decoupling of element unknowns. Element unknowns can then be removed from global linear systems by static condensation reducing the globally coupled ndofs (asymptotically for increasing p) from O(p d ) to O(p d−1 ) and thereby reducing the computational costs dramatically.
Let us sketch the structure of an HDG method. To this end, we restrict to a scalar PDE here, but extensions to the vector case are obvious. Let F p (F h ) be the space of piecewise polynomials up to degree p on each facet of the mesh T h and F p D (F h ) and F p 0 (F h ) its subspaces with prescribed (inhomogeneous and homogeneous) values on Dirichlet-type boundaries. Then a typical discrete variational formulation of a Hybrid DG method is formulated in terms of the pair of volume and facet unknowns in the form: Here the bilinear form a h (·, ·) can be very similar to a corresponding DG formulation with respect to the integral terms. However, direct couplings between volume unknowns of adjacent elements are avoided by involving the facet unknowns for the inter-element communication, so that the set of volume unknowns of one element is completely determined by the facet unknowns on the element boundary (and the r.h.s. data). This enables one of the main features of Hybrid DG schemes: the possibility to do static condensation. This allows to eliminate the interior unknowns in V p (T h ) completely based on the facet unknowns in F p (F h ). This can be done on the level of the variational formulation where with u = u(u F ) one can formulate a discrete formulation solely based on u F of the form: After solving for u F element-local problems can be solved to re-obtain u(u F ) ∈ V p (T h ). This elimination is often done merely on the linear algebra level based on a Schur complement strategy. We note that in many Hybrid DG formulations for second order PDEs an auxiliary (flux) variable is also introduced locally which however can be eliminated alongside u so that the structure of the global linear system for u F is not affected by this.
The asymptotic complexity of the Hybrid DG approach is the same as the one of the Embedded Trefftz DG method. The total ndofs is determined by the element unknowns, i.e. O(h −d p d ), the globally coupled ndofs however only scales with O(h −d p d−1 ). Similarly the local operations, i.e. the SVD or QR decomposition for the Embedded Trefftz or the Schur complement strategy require O(h −d p 3d ) operations. Correspondingly, nnzesand arithmetic operations for general purpose linear algebra solvers have the same complexity, cf. Section 3.6.
We summarize that the Hybrid DG and the embedded Trefftz methods achieve the same asymptotic complexity, however the global unknowns are associated with element interfaces in the Hybrid DG case and associated with volume elements in the case of embedded Trefftz methods. In Fig. 12 Figure 12. Sketch of global/local dofs of Hybrid DG, standard DG and Embedded Trefftz DG for 2D simplex elements and polynomials degree p = 3. Here q diff is the maximum differentiation index of L.
In the first part of Section 3.3 we have discussed two possible implementations of the embedded Trefftz-DG method. For the timings in this section we use the faster method, which avoids the assembly of the full DG matrix.

Poisson equation.
To compare the Hybrid DG method with the Embedded Trefftz DG method on a concrete example, we reconsider the example from Section 4.3 and provide a simple hybrid DG discretization for it. A hybrid DG analogue to the interior penalty formulation takes the form of (36) with Note that there are no average and jump operators in (37) that directly involve neighboring element functions. Further, Dirichlet-type boundary conditions are directly imposed on the facet space which leads to the usual weak imposition of them on the element unknowns. As for the DG formulation chosen in Section 4.3 let us also mention that the Hybrid DG formulation is only one specific version (and a specifically simple one). For second order elliptic and diffusion dominated problems, based on the similarity to (hybrid) mixed methods, element-local postprocessing schemes can be devised to obtain an additional order of accuracy for a postprocessed solution field. We neglect this aspect here but refer to [8,30,31,48] for more details.  In Fig. 13 we compare the computational results for the three methods, DG, Hybrid DG and Embedded Trefftz DG for fixed polynomial degrees and successively refined meshes. We observe that the differences in the error is negligible whereas the computation time of the Hybrid DG method is still smaller than the one of the Embedded Trefftz DG method. Comparing the overall computation time with the time spend only to solve the global linear system we observe that this makes up for most of the time.

Linear transport equation.
Let us now come to the first order example, the linear transport equation as in Section 4.8. A Hybrid DG version of (35) takes the form of (36) with Here the upwind choiceû is as follows: On inflow boundaries (b · n x < 0) we setû = u F , i.e. the facet value is taken whereas on outflow boundaries (b · n x > 0) we setû = u. To ensure that the facet value taken on inflow boundaries is meaningful, it is glued together with the outflow trace of the corresponding adjacent element by the stabilization term that include v F . Note that testing with v = 0 yields that u F is exactly the upwind trace at a corresponding facet. Hence, the hybrid Upwind DG formulation is equivalent to the Upwind DG formulation used above, cf. [13]. Note that inflow boundary conditions are prescribed through the facet space here and do not appear in the linear form.  In Fig. 14 we again compare DG, Hybrid DG and Embedded Trefftz DG methods. Again, the errors are comparable. In fact, DG and Hybrid DG are equivalent and hence yield exactly the same result. For the computation time we now observe that the Embedded Trefftz DG method is the fastest solution method as soon as the overhead in the system setup becomes negligible. 5.4. Some number crunching. In the previous two sections we considered first and second order scalar PDEs. Next, we would like to compare the sparsity patterns of DG, HDG and Trefftz DG for a generic setting: We consider an unstructured simplicial mesh with 54 elements in 2D and 729 elements in 3D and compare different methods w.r.t. the number of degrees of freedom (ndofs) and the number of non-zero entries (nnzes) for p = 0, .., 5. We note that although ndofs is the simpler measure nnzes has a more direct implication on the computational costs that are associated to solving linear systems. In Table 1 A standard DG scheme, a corresponding HDG scheme and two embedded Trefftz DG schemes are compared. Here the embedded Trefftz DG schemes denoted as TDG(1) and TDG(2) are distinguished depending on the leading order Table 1. Comparison of different measures (ndofs, nnzes) for the computational overhead for the solution of the largest global linear system related to the methods DG, Hybrid DG (after static condenstaion) and (embedded) Trefftz DG.
of the differential operator involved. For first order problems a comparison of DG and HDG with TDG(1) is of interest, whereas TDG(2) is to be considered for comparison for second order formulations. For simplicity, we assume here that all element and facet unknowns of one element or facet couple with all unknowns of a corresponding facet or element (depending on the method) neighbor.
We observe that for first order problems, the embedded Trefftz DG method even outperforms the HDG method. For second order problems HDG and embedded Trefftz DG methods are similar, in 2D HDG is cheaper, in 3D the embedded Trefftz DG method. However, optimized HDG formulations are able to reduce the polynomial degree on the element interfaces by one order in many cases without effecting the order of accuracy. In these cases the HDG method will still be cheaper. However, compared to the plain DG method the improvement of the embedded Trefftz DG method is already remarkable.

Conclusion and possible extensions
We have presented a method to reduce the matrix of a DG scheme post assembly, by projecting the polynomial basis onto Trefftz spaces. Several numerical examples have been presented, showing that the method matches the convergence rates of polynomial Trefftz spaces and showcasing its flexibility and potential in case of smooth coefficients and inhomogeneous equations.
While the Trefftz spaces used in the polynomial case are well understood, analytical properties of the space used for the embedding in the case of smooth coefficients and plane wave Trefftz spaces will be addressed in a forthcoming paper. One remarkable finding are the convincing results that have been obtained for the linear transport problem, a problem that is typically not associated with Trefftz methods. Further, we saw that the method is more than an acceleration technique, as the restriction to a subspace can even improve stability as seen for the Helmholtz problem. The flexibility of the approach suggests to consider the approach in many more and possibly more complex cases.
We focused on using the local operator and its kernel to emulate Trefftz spaces. However, the technique could also be applied to impose other constraints on the solution as an alternative to complicated constructions of special basis functions, weak impositions through the DG formulation or Lagrange multiplier techniques. Furthermore, the extraction done here for the kernel could similarly be applied for the range of certain differential operators or orthogonal (w.r.t. to a localizable inner product) complements.
One obvious drawback is the confinement to DG schemes. It presents an alternative to other approaches that accelerate DG methods, such as HDG. However, using the method as a pure acceleration technique for DG methods may be less attractive, if sophisticated preconditioners and linear solvers for the DG discretization are available. Moreover, the scaling of the computing time of the embedding construction in p is not good, suggesting that really high polynomial degrees, e.g. beyond p = 10, may not be feasible. In these cases however, the method seems to at least offer a very convenient way to investigate Trefftz and Trefftz-type methods as a research tool.
Another difficulty is, that combining the approach with (partially) conforming methods (such as those based on H(div) or H(curl) spaces) is not directly possible, due to the support of basis functions spanning multiple elements and overlapping partially the construction of a local embedding unclear. Additionally, the application of the PDE operator in strong form is only feasible element-wise.