<figure><img src="image_1.png" /></figure>

HAL
open science

# A conservative Eulerian-Lagrangian Runge-Kutta discontinuous Galerkin method for linear hyperbolic system with large time stepping size.

Xue Hong, Jingmei Qiu

► To cite this version:

Xue Hong, Jingmei Qiu. A conservative Eulerian-Lagrangian Runge-Kutta discontinuous Galerkin method for linear hyperbolic system with large time stepping size.. 2023. hal-03920844

HAL Id: hal-03920844

https://hal.science/hal-03920844

Preprint submitted on 4 Jan 2023

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.
# A conservative Eulerian-Lagrangian Runge-Kutta discontinuous Galerkin method for linear hyperbolic system with large time stepping*

Xue Hong†

Jing-Mei Qiu‡

January 3, 2023

## Abstract

We propose an Eulerian-Lagrangian (EL) Runge-Kutta (RK) discontinuous Galerkin (DG) method for linear hyperbolic system. The method is designed based on the EL DG method for transport problems [J. Comput. Phy. 446: 110632, 2021.], which tracks solution along approximations to characteristics in the DG framework, allowing extra large time stepping sizes with stability with respect to the classical RK DG method. Considering each characteristic family, a straightforward application of EL DG for hyperbolic system will be to transform to the characteristic variables, evolve them on associated characteristic related space-time regions, and transform them back to the original variables. However, the conservation could not be guaranteed in a general setting. In this paper, we formulate a conservative semi-discrete EL DG method by decomposing each variable into two parts, each of them associated with a different characteristic family. As a result, four different quantities are evolved in EL fashion and recombinced to update the solution. The fully discrete scheme is formulated by using method-of-lines RK methods, with intermediate RK solutions updated on the background mesh. Numerical results for 1D and 2D wave equations are presented to demonstrate the performance of the proposed ELDG method. These include the high order spatial and temporal accuracy, stability with extra large time stepping size, and conservative property.

**Key words:** Eulerian-Lagrangian; discontinuous Galerkin; characteristic method; linear hyperbolic system; conservative property.

## 1 Introduction

In this paper, we propose an Eulerian-Lagrangian (EL) discontinuous Galerkin (DG) method for the first-order hyperbolic system in the form of

$$U_t + \sum_{j=1}^{d} (A_j(\mathbf{x}, t)U)_{x_j} = F(\mathbf{x}, t), \quad (\mathbf{x}, t) \in \mathbb{R}^d \times [0, T], \qquad (1.1)$$

<footnote>*Research of the first author is supported by the Centre Henri Lebesgue, program ANR-11-LABX-0020-0 and the Brittany council. Research of the second author is supported by NSF grant NSF-DMS-2111253, Air Force Office of Scientific Research FA9550-18-1-0257 and University of Delaware.</footnote>

<footnote>†Université de Rennes, Inria Rennes (Mingus team) and IRMAR UMR CNRS 6625, F-35042 Rennes, France. (xue.hong@inria.fr).</footnote>

<footnote>‡Department of Mathematical Sciences, University of Delaware, Newark, DE, 19716, USA. (jingqiu@udel.edu).</footnote>
where $d$ is the spatial dimension, $U : \mathbb{R}^d \times [0, T] \rightarrow \mathbb{R}^n$, and $A_j(\mathbf{x}, t) \in \mathbb{R}^{n \times d}$. Several systems of this form including wave equations, Maxwell's equations, linearized shallow water equations etc. which have important physical applications [16]. There are many efficient DG methods for solving the hyperbolic systems. A particularly powerful combination is applying the DG method in space and a Runge-Kutta scheme in time to obtain RK DG scheme. There are a series of works about DG [8] and space-time DG [10, 25] methods.

One effective numerical approach for solving the hyperbolic equation is the characteristic method, which updates time-dependent solution by tracking characteristics. For scalar hyperbolic equations, there have been many pioneering works. In [6], Celia etc. developed an Eulerian Lagrangian Localized Adjoint Methods (ELLAM) [6], which introduces an adjoint problem for each test function in the continuous finite element framework and has been applied to different problems [29, 28]. On the other hand, the EL DG [30, 5] and SL DG [3] are being developed in the discontinuous Galerkin finite element framework with a similar introduction of adjoint problems for test functions. More important, the EL DG method takes the advantage of linear approximation of the characteristics, which make it have stability with extra large time stepping size $\Delta t \sim \frac{\Delta x}{(1+k) \max |a-\alpha|}$ compared with $\Delta t \sim \frac{\Delta x}{(1+k) \max |a|}$ of the classical Eulerian explicit RK DG method for linear advection problems, where $a$ is the velocity of advection problem and $\alpha$ is a linear approximation of $a$. So the time step size of EL DG with stability can be very large when $\alpha$ approximate $a$ well. The EL DG method is also closely related to the Arbitrary Lagrangian Eulerian (ALE) DG method [15, 13] which is a moving mesh DG method.

In this paper, we propose a conservative EL RK DG method for hyperbolic systems (1.1) ($F(\mathbf{x}, t) = 0$) with large time stepping sizes compared with that for classical RKDG method. We start from 1D cases, for which we consider characteristic variables and the associated characteristic space-time regions. For each characteristic family, a straightforward application of EL DG will be to transform to characteristic variables, evolve them on associated space-time regions, and transform them back to original variables. However, the conservation could not be guaranteed in a general setting. We decompose each variable into two parts, each of them associated with different characteristic families; as a result, four different quantities are evolved in EL fashion and recombined to update the solution. The fully discrete scheme is formulated by using method-of-lines RK methods, with intermediate RK solutions updated on the background mesh. For 2D hyperbolic
systems, characteristic variables are no longer constant along the bicharacteristics, characteristic Galerkin [26] or evolution Galerkin [2, 18, 19] methods have been proposed by taking into account information propagated in all bicharacteristic directions. However, the algorithm implementation is very complex. In this paper, we use the dimensional splitting method for higher dimensional problem, which maintain the simplicity of EL DG method for 1D cases, as well as other great properties such as conservation, high order spatial and temporal accuracy, and allows for extra large time steps with stability.

This paper is organized as follows. In Section 2, we review the EL DG for one-dimensional (1D) linear transport problems. In Section 3, we develop the EL DG method for a first-order hyperbolic system by evolving each component associated with each characteristic families, recombining them to update the solution. We also develop it to 2D problems by dimensional splitting. Conservation of ELDG schemes are proved. In Section 4, performance of the proposed ELDG method is shown through extensive numerical tests. Finally, concluding remarks are made in Section 5.

## 2 Review of EL DG formulation for 1D linear transport problems

To illustrate the key idea of the EL DG scheme, we start from a 1D linear transport equation in the following form

$$u_t + (a(x,t)u)_x = 0, \quad x \in [x_a, x_b]. \qquad (2.1)$$

For simplicity, we assume periodic boundary conditions, and the velocity field $a(x,t)$ is a continuous function of space and time. $\Omega$ is the whole computational domain. We perform a partition of the computational domain $x_a = x_{1/2} < x_{3/2} < \cdots < x_{N+1/2} = x_b$ as the background mesh partition. Let $I_j = [x_{j-1/2}, x_{j+1/2}]$ denote an element of length $\Delta x_j = x_{j+1/2} - x_{j-1/2}$ and define $\Delta x = \max_j \Delta x_j$. We define the finite dimensional approximation space, $V_h^k = \{v_h : v_h|_{I_j} \in P^k(I_j)\}$, where $P^k(I_j)$ denotes the set of polynomials of degree at most $k$. We let $t^n$ be the $n$-th time level and $\Delta t = t^{n+1} - t^n$ to be the time-stepping size. There are a lot of notations we need to mention in this paper, now we use $A \doteq B$ to denote $A$ as $B$ for simplicity.

The key idea in the EL DG formulation is, design adjoint problems for test functions to take advantage of information propagation along characteristics. The EL DG method proposed in [5] is formulated on a space-time region $\Omega_j = \tilde{I}_j(t) \times [t^n, t^{n+1}]$ with

$$\tilde{I}_j(t) = [\tilde{x}_{j-1/2}(t), \tilde{x}_{j+1/2}(t)], \quad t \in [t^n, t^{n+1}]$$
being the dynamic interval with size $\Delta x_j(t) = \tilde{x}_{j+\frac{1}{2}}(t) - \tilde{x}_{j-\frac{1}{2}}(t)$, see Figure 2.1. Here $\tilde{x}_{j\pm\frac{1}{2}}(t) = x_{j\pm\frac{1}{2}} + (t - t^{n+1})\nu_{j\pm\frac{1}{2}}$ are straight lines emanating from cell boundaries $x_{j\pm\frac{1}{2}}$ with slopes $\nu_{j\pm\frac{1}{2}} = a(x_{j\pm\frac{1}{2}}, t^{n+1})$ and $I_j^* = \tilde{I}_j(t^n) = [x_{j-\frac{1}{2}}^*, x_{j+\frac{1}{2}}^*]$ is the upstream cell of $I_j$ at $t^n$. The dynamic interval of $\tilde{I}_j(t)$ can always be linearly mapped to a reference cell $I_j$ by a mapping. A local adjoint problem

<figure><img src="image_2.png" /><figcaption>Figure 2.1: Illustration for dynamic element $\tilde{I}_j(t)$ of ELDG.</figcaption></figure>

of (2.1) for all test function is defined as:

$$ \left\{ \begin{array}{l} \psi_t + \alpha(x,t)\psi_x = 0, \quad (x,t) \in \Omega_j, \\ \psi(t = t^{n+1}) = \Psi(x). \quad \forall \Psi \in P^k(I_j). \end{array} \right. \qquad (2.2) $$

Here $\alpha(x,t)$ is a bilinear function of $(x,t)$ with $\forall t \in [t^n, t^{n+1}]$,

$$ \alpha(x_{j-\frac{1}{2}}, t^{n+1}) = a(x_{j-\frac{1}{2}}, t^{n+1}) \doteq \nu_{j-\frac{1}{2}}, \quad \alpha(x_{j+\frac{1}{2}}, t^{n+1}) = a(x_{j+\frac{1}{2}}, t^{n+1}) \doteq \nu_{j+\frac{1}{2}}, \qquad (2.3) $$

and

$$ \alpha(x,t) = -\nu_{j-\frac{1}{2}} \frac{x - \tilde{x}_{j+\frac{1}{2}}(t)}{\Delta x_j(t)} + \nu_{j+\frac{1}{2}} \frac{x - \tilde{x}_{j-\frac{1}{2}}(t)}{\Delta x_j(t)} \in P^1(\tilde{I}_j(t)). \qquad (2.4) $$

Notice that the test function $\psi$ stays the same polynomial, if $\tilde{I}_j(t)$ is mapped to a reference interval $I_j$. The EL DG [5] scheme can be formulated by $\int_{\tilde{I}_j(t)} (2.2) \cdot u + (2.1) \cdot \psi$

$$ \frac{d}{dt} \int_{\tilde{I}_j(t)} (u\psi) dx = - (\hat{F}\psi) \Big|_{\tilde{x}_{j+\frac{1}{2}}(t)} + (\hat{F}\psi) \Big|_{\tilde{x}_{j-\frac{1}{2}}(t)} + \int_{\tilde{I}_j(t)} F\psi_x dx. \qquad (2.5) $$

where $F(u) \doteq (a - \alpha)u$ and $\hat{F}$ is the Lax-Friedrichs flux. A method-of-lines RK discretization can be used for high order temporal accuracy [5].

## 3 The EL DG algorithm for hyperbolic system

The design of the EL DG algorithm for hyperbolic system shares a similar spirit as the 1D scalar case. We start from a 1D hyperbolic system. We firstly formulate the EL DG scheme by tracking information along different characteristics families.
## 3.1 1D hyperbolic system

We consider the hyperbolic system

$$U_t + (A(x)U)_x = F(x, t), \qquad (3.1)$$

where $U = [u^1, u^2]^T$ is a column vector, $A$ is a 2 by 2 matrix, $F$ is also a 2 by 1 vector. We use the following notations for the eigen-decomposition of $A(x)$:

* eigenvalue: $\lambda^{(1)}(x)$, $\lambda^{(2)}(x)$.

* $A(x) = R(x)\Lambda R^{-1}(x)$, where $\Lambda(x) = diag(\lambda^{(1)}(x), \lambda^{(2)}(x))$,

$$R(x) \doteq [r^{(1)}(x) | r^{(2)}(x)] = \begin{bmatrix} r_{11}(x) & r_{12}(x) \\ r_{21}(x) & r_{22}(x) \end{bmatrix} \qquad (3.2)$$

contains the right column eigenvectors $r^{(1)}$, $r^{(2)}$, and

$$R^{-1}(x) \doteq \begin{bmatrix} l^{(1)T}(x) \\ l^{(2)T}(x) \end{bmatrix} = \begin{bmatrix} l_{11}(x) & l_{12}(x) \\ l_{21}(x) & l_{22}(x) \end{bmatrix} \qquad (3.3)$$

contains the left row eigenvectors $l^{(1)T}$, $l^{(2)T}$.

In the following, we propose a conservative EL DG scheme for the system (3.1) by the procedure below. The notion of its background meshes is the same as the 1D scalar case.

### (1) Two partitions of space-time regions $\Omega_j^{(1)}$ and $\Omega_j^{(2)}$
According to the first and second characteristic families, we partition the computational domain as two sets of space-time regions $\Omega_j^{(1)}$ and $\Omega_j^{(2)}$ respectively. Here $\Omega_j^{(1)} = \tilde{I}_j^{(1)}(t) \times [t^n, t^{n+1}]$ is related to the first characteristic family. $\tilde{I}_j^{(1)}(t) = [\tilde{x}_{j-\frac{1}{2}}^{(1)}(t), \tilde{x}_{j+\frac{1}{2}}^{(1)}(t)]$ is the dynamic interval emanating from cell boundaries $x_{j\pm\frac{1}{2}}$ with slopes $\nu_{j\pm\frac{1}{2}}^{(1)}$ approximating the first characteristic velocity, see Figure 3.2 (left). In general, we choose $\nu_{j\pm\frac{1}{2}}^{(1)} = \lambda^{(1)}(x_{j\pm\frac{1}{2}})$. $I_j^{*,(1)} \doteq \tilde{I}_j^{(1)}(t^n)$ is the upstream cell of $I_j$ from the first characteristic family at $t^n$. The dynamic interval $\tilde{I}_j^{(1)}(t)$ can be linearly mapped to a reference cell $\xi \in I_j$, (see Figure 2.1). Here, we let $\tilde{x}^{(1)}(\tau; (\xi, t^{n+1}))$ be the linear map from $\tilde{I}_j^{(1)}(t)$ to $I_j$. Similar definition can be made to $\Omega_j^{(2)}$, $\tilde{I}_j^{(2)}(t)$ and $I_j^{*,(2)}$ for the second characteristic family. See Figure 3.2 (right) for illustration of $\Omega_j^{(2)}$.

### (2) Adjoint Problems
We consider an adjoint problem for the first characteristic family on $\Omega_j^{(1)}$:

$$\begin{cases} (\psi^{(1)})_t + \alpha^{(1)}(\psi^{(1)})_x = 0, \quad t \in [t^n, t^{n+1}], \\ (\psi^{(1)})(t = t^{n+1}) = \Psi^{(1)}(x), \end{cases} \qquad (3.4)$$
<figure><img src="image_3.png" /><figcaption>Figure 3.2: Illustration for dynamic elements $\tilde{I}_j^{(1)}(t)$ (left) and $\tilde{I}_j^{(2)}(t)$ (right) of ELDG for the first and second characteristic families of the system.</figcaption></figure>

where

$$ \alpha^{(1)}(x,t) = -\nu_{j-\frac{1}{2}}^{(1)} \frac{x - \tilde{x}_{j+\frac{1}{2}}^{(1)}(t)}{\Delta x_j^{(1)}(t)} + \nu_{j+\frac{1}{2}}^{(1)} \frac{x - \tilde{x}_{j-\frac{1}{2}}^{(1)}(t)}{\Delta x_j^{(1)}(t)} \in P^1(\tilde{I}_j^{(1)}(t)). \quad (3.5) $$

Similarly on $\Omega_j^{(2)}$:

$$ \left\{ \begin{array}{l} (\psi^{(2)})_t + \alpha^{(2)}(\psi^{(2)})_x = 0, \quad t \in [t^n, t^{n+1}], \\ (\psi^{(2)})(t = t^{n+1}) = \Psi^{(2)}(x), \end{array} \right. \quad (3.6) $$

where

$$ \alpha^{(2)}(x,t) = -\nu_{j-\frac{1}{2}}^{(2)} \frac{x - \tilde{x}_{j+\frac{1}{2}}^{(2)}(t)}{\Delta x_j^{(2)}(t)} + \nu_{j+\frac{1}{2}}^{(2)} \frac{x - \tilde{x}_{j-\frac{1}{2}}^{(2)}(t)}{\Delta x_j^{(2)}(t)} \in P^1(\tilde{I}_j^{(2)}(t)). \quad (3.7) $$

The adjoint problems provide finite dimensional time-dependent test function space, please see more details in the Appendix.

### (3) Formulation of a conservative semi-discrete ELDG scheme.

For linear hyperbolic system, a straightforward generalization of ELDG is to transform the original variable to the characteristic variables by a localized cell dependent approximating eigen-decomposition. So we firstly formulate a EL DG scheme by a localized characteristic field, but that is proved not conservative. We show the details in Appendix. To obtain the conservation, a critical point is that eigen-decomposition has to be consistent among two characteristic families and independent of partitions $\Omega_j^{(1)}, \Omega_j^{(2)}$. We directly use the exact eigen-decomposition (3.2) and (3.3). The following equalities hold

$$ \begin{aligned} & (r_{11}(x)l^{(1)T}(x) + r_{12}(x)l^{(2)T})U(x) = u^1, \\ & (r_{21}(x)l^{(1)T}(x) + r_{22}(x)l^{(2)T})U(x) = u^2, \end{aligned} \quad (3.8) $$

from $R(x)R^{-1}(x) = R^{-1}(x)R(x) = I$. They are critical to design a conservative ELDG scheme.
Take the vector product of $r_{11}(x)l^{(1)T}(x)$ from left with (3.1), we have a scalar equation

$$r_{11}(x)l^{(1)T}(x)(U_t + (A(x)U)_x) = r_{11}(x)l^{(1)T}(x)F(x,t). \quad (3.9)$$

Multiply $\psi^{(1)}$ to the equation above,

$$r_{11}(x)l^{(1)T}(x)(U_t + (A(x)U)_x)\psi^{(1)} = r_{11}(x)l^{(1)T}(x)F(x,t)\psi^{(1)}. \quad (3.10)$$

Meanwhile, from $r_{11}(x)l^{(1)T}(x)U \cdot (3.4)$, we have

$$r_{11}(x)l^{(1)T}(x)U(\psi^{(1)})_t + r_{11}(x)l^{(1)T}(x)U\alpha^{(1)}(\psi^{(1)})_x = 0. \quad (3.11)$$

Next, sum equations (3.10) and (3.11), and integrate over the space-time interval $\Omega_j^{(1)}$, then

$$\int_{\Omega_j^{(1)}} \left( r_{11}(x)l^{(1)T}(x)U_t\psi^{(1)} + r_{11}(x)l^{(1)T}(x)U(\psi^{(1)})_t + r_{11}(x)l^{(1)T}(x)(A(x)U)_x\psi^{(1)} \right) dxdt \\ + \int_{\Omega_j^{(1)}} r_{11}(x)l^{(1)T}(x)U\alpha^{(1)}(\psi^{(1)})_x dxdt = \int_{\Omega_j^{(1)}} r_{11}(x)l^{(1)T}(x)F(x,t)\psi^{(1)}dxdt. \quad (3.12)$$

A further manipulation on the left hand side (L.H.S.) of (3.12) gives

$$\begin{aligned} & \int_{\Omega_j^{(1)}} \left( (r_{11}(x)l^{(1)T}(x)U\psi^{(1)})_t + r_{11}(x)l^{(1)T}(x)(A(x)U)_x\psi^{(1)} + r_{11}(x)l^{(1)T}(x)U\alpha^{(1)}(\psi^{(1)})_x \right) dxdt \\ = & \int_{\Omega_j^{(1)}} \left( (r_{11}(x)l^{(1)T}(x)U\psi^{(1)})_t + (r_{11}(x)l^{(1)T}(x)A(x)U\psi^{(1)})_x - (r_{11}(x)l^{(1)T}(x))_x A(x)U\psi^{(1)} \right) dxdt \\ - & \int_{\Omega_j^{(1)}} \left( r_{11}(x)l^{(1)T}(x)A(x)U(\psi^{(1)})_x - r_{11}(x)l^{(1)T}(x)U\alpha^{(1)}(\psi^{(1)})_x \right) dxdt \\ = & \int_{t^n}^{t^{n+1}} \left( \frac{d}{dt} \int_{\tilde{I}_j^{(1)}(t)} (r_{11}(x)l^{(1)T}(x)U\psi^{(1)})dx + [r_{11}(x)l^{(1)T}(x)A(x)U\psi^{(1)} - \nu^{(1)}r_{11}(x)l^{(1)T}(x)U\psi^{(1)}]_{j-\frac{1}{2}}^{j+\frac{1}{2}} \right) dt \\ - & \int_{t^n}^{t^{n+1}} \int_{\tilde{I}_j^{(1)}(t)} \left( (r_{11}(x)l^{(1)T}(x))_x A(x)U\psi^{(1)} + r_{11}(x)l^{(1)T}(x)(A(x)U - \alpha^{(1)}U)\psi_x^{(1)} \right) dxdt. \end{aligned} \quad (3.13)$$

Letting $f^{11}(U) = r_{11}l^{(1)T}(AU - \alpha^{(1)}U)$, the time differential form of (3.12) with (3.13) gives

$$\begin{aligned} & \frac{d}{dt} \int_{\tilde{I}_j^{(1)}(t)} (r_{11}(x)l^{(1)T}(x)U\psi^{(1)})dx + (f^{11}\psi^{(1)})\Big|_{\tilde{x}_{j+\frac{1}{2}}^{(1)}(t)} - (f^{11}\psi^{(1)})\Big|_{\tilde{x}_{j-\frac{1}{2}}^{(1)}(t)} - \int_{\tilde{I}_j^{(1)}(t)} f^{11}\psi_x^{(1)}dx \\ & - \int_{\tilde{I}_j^{(1)}(t)} (r_{11}(x)l^{(1)T}(x))_x A(x)U\psi^{(1)}dx = \int_{\tilde{I}_j^{(1)}(t)} r_{11}(x)l^{(1)T}(x)F(x,t)\psi^{(1)}dx. \end{aligned} \quad (3.14)$$
Similarly, we have an equation related to $\lambda^{(2)}$

$$ \frac{d}{dt} \int_{\tilde{I}_j^{(1)}(t)} (r_{12}(x)l^{(2)T}(x)U\psi^{(2)})dx + (f^{12}\psi^{(2)}) \left| \tilde{x}_{j+\frac{1}{2}}^{(2)}(t) - (f^{12}\psi^{(2)}) \right| \tilde{x}_{j-\frac{1}{2}}^{(2)}(t) - \int_{\tilde{I}_j^{(2)}(t)} f^{12}\psi_x^{(2)}dx \\ - \int_{\tilde{I}_j^{(2)}(t)} (r_{12}(x)l^{(2)T}(x))_x A(x)U\psi^{(2)}dx = \int_{\tilde{I}_j^{(2)}(t)} r_{12}(x)l^{(2)T}(x)F(x,t)\psi^{(2)}dx, \qquad (3.15) $$

where $f^{12}(U) = r_{12}l^{(2)T}(AU - \alpha^{(2)}U)$. Then, we can update $u^1$ by (3.8), (3.14), (3.15), taking $\Psi^{(1)}(x) = \Psi(x)$ in (3.4) and $\Psi^{(2)}(x) = \Psi(x)$ in (3.6):

$$ \int_{I_j} u^{1,n+1} \Psi(x) dx \stackrel{(3.8)}{=} \int_{I_j} r_{11} l^{(1)T} U^{n+1} \Psi(x) dx + \int_{I_j} r_{12} l^{(2)T} U^{n+1} \Psi(x) dx \\ = \int_{I_j^{*,(1)}} r_{11} l^{(1)T} U^n \psi^{(1)} dx - \int_{t^n}^{t^{n+1}} (f^{11}\psi^{(1)}) \left| \tilde{x}_{j+\frac{1}{2}}^{(1)}(t) + (f^{11}\psi^{(1)}) \right| \tilde{x}_{j-\frac{1}{2}}^{(1)}(t) dt \\ + \int_{t^n}^{t^{n+1}} \int_{\tilde{I}_j^{(1)}(t)} r_{11}(x) l^{(1)T}(x) F(x,t) \psi^{(1)} + (r_{11}(x)l^{(1)T}(x))_x A(x) U\psi^{(1)} + f^{11}\psi_x^{(1)} dxdt \qquad (3.16) \\ + \int_{I_j^{*,(2)}} r_{12} l^{(2)T} U^n \psi^{(2)} dx - \int_{t^n}^{t^{n+1}} (f^{12}\psi^{(2)}) \left| \tilde{x}_{j+\frac{1}{2}}^{(2)}(t) + (f^{12}\psi^{(2)}) \right| \tilde{x}_{j-\frac{1}{2}}^{(2)}(t) dt \\ + \int_{t^n}^{t^{n+1}} \int_{\tilde{I}_j^{(2)}(t)} r_{12}(x) l^{(2)T}(x) F(x,t) \psi^{(2)} + (r_{12}(x)l^{(2)T}(x))_x A(x) U\psi^{(2)} + f^{12}\psi_x^{(2)} dxdt. $$

As the EL DG method for scalar equations, the dynamic interval of $I_j^{(1)}(t)$ and $I_j^{(2)}(t)$ can always be linearly mapped to a reference cell $\xi \in I_j$, the EL DG discretization of eq. (3.16) is to find $u_h^1(x, t) \in P^k(I_j(t))$, so that

$$ \int_{I_j} u_h^1(x, t^{n+1}) \Psi(x) dx = \int_{I_j^{*,(1)}} r_{11} l^{(1)T} U_h^n \psi^{(1)} dx + \int_{t^n}^{t^{n+1}} L_{11}(U_h(t), t, \tilde{I}_j^{(1)}(t)) dt \\ + \int_{I_j^{*,(2)}} r_{12} l^{(2)T} U_h^n \psi^{(2)} dx + \int_{t^n}^{t^{n+1}} L_{12}(U_h(t), t, \tilde{I}_j^{(2)}(t)) dt, \qquad (3.17) $$

for $\psi^{(1)}(x, t)$ satisfying the adjoint problem (3.4) with $\forall \Psi(x) = \psi(x, t^{n+1}) \in P^k(I_j)$. Here

$$ \begin{aligned} L_{11}(U_h(t), t, \tilde{I}_j^{(1)}(t)) = & -\widehat{f_{j+\frac{1}{2}}^{11}} \psi_{j+\frac{1}{2}}^{(1),-} + \widehat{f_{j-\frac{1}{2}}^{11}} \psi_{j-\frac{1}{2}}^{(1),+} + \int_{\tilde{I}_j^{(1)}(t)} f^{11}\psi_x^{(1)}(x,t)dx \\ & + \int_{\tilde{I}_j^{(1)}(t)} (r_{11}l^{(1)T})_x A U_h \psi^{(1)}(x,t) + r_{11}l^{(1)T} F\psi^{(1)}(x,t)dx, \end{aligned} \qquad (3.18) \\ L_{12}(U_h(t), t, \tilde{I}_j^{(2)}(t)) = & -\widehat{f_{j+\frac{1}{2}}^{12}} \psi_{j+\frac{1}{2}}^{(2),-} + \widehat{f_{j-\frac{1}{2}}^{12}} \psi_{j-\frac{1}{2}}^{(2),+} + \int_{\tilde{I}_j^{(2)}(t)} f^{12}\psi_x^{(2)}(x,t)dx \\ & + \int_{\tilde{I}_j^{(2)}(t)} (r_{12}l^{(2)T})_x A U_h \psi^{(2)}(x,t) + r_{12}l^{(2)T} F\psi^{(2)}(x,t)dx, \end{aligned} $$
where $\psi_{j\pm\frac{1}{2}}^{(1),\pm} = \psi^{(1)}(x_{j\pm\frac{1}{2}}^{(1),\pm}(t), t)$ and $\widehat{f^{11}}$ at a cell boundary can be taken as a monotone flux, e.g. the Lax-Friedrichs flux

$$ \widehat{f_{j+\frac{1}{2}}^{11}}(U) = r_{11}(x_{j+\frac{1}{2}})l^{(1)T}(x_{j+\frac{1}{2}})(\widehat{A \cdot U - \nu_{j+\frac{1}{2}}^{(1)} U})_{j+\frac{1}{2}}. $$

$\psi_{j\pm\frac{1}{2}}^{(2),\pm} = \psi^{(2)}(x_{j\pm\frac{1}{2}}^{(2),\pm}(t), t)$ and $\widehat{f^{12}}$ can be similarly defined at a cell boundary. We can similarly obtain the EL DG scheme for $u_h^2$.

### (4) Fully discrete EL DG scheme with method-of-lines RK schemes.

To update (3.17) from $U_h^n$ to $U_h^{n+1}$, we first apply the forward Euler time discretization to get 1st order accuracy, then we generalize the scheme to general RK methods. There are two main steps involved here. In order to describe the implementation procedure of the fully discrete EL DG scheme, we define the $L^2$ projection.

**Definition 3.1.** ($L^2$ projection) Let $u \in L^2(\Omega)$, $M = \{I_j\}_{j=1}^N$ and $\tilde{M} = \{\tilde{I}_j\}_{j=1}^N$ be two different meshes of the whole computational domain $\Omega$. We have function spaces $V_h^k = \{u : u|_{I_j} \in P^k(I_j), \forall j\}$ and $\tilde{V}_h^k = \{\tilde{u} : \tilde{u}|_{\tilde{I}_j} \in P^k(\tilde{I}_j), \forall j\}$ corresponding to meshes $M$ and $\tilde{M}$. The $L^2$ projection of $u_M \in V_h^k$ onto space $\tilde{V}_h^k$ can be defined as, find $\tilde{u}_M \in \tilde{V}_h^k$, s.t.

$$ \int_{\tilde{I}_j} \tilde{u}_M(x)\varphi(x)dx = \int_{\tilde{I}_j} u_M(x)\varphi(x)dx, \forall \varphi \in \tilde{V}_h^k. \quad (3.19) $$

We denote $\tilde{u}_M(x) = \text{Proj}[u_M(x); M, \tilde{M}]$. The evaluation of the right hand side of (3.19) can be done in a subinterval-by-subinterval fashion. The implementation details can be found in [12].

Then, we propose a fully discrete EL RK DG scheme with procedure as described:

1. **Obtain the initial condition** on upstream meshes $\tilde{I}_j^{(1)}(t^n)$ and $\tilde{I}_j^{(2)}(t^n)$ of (3.17) by $U_h^{(1)}(t^n) = \text{Proj}[U_h^n; I_j, \tilde{I}_j^{(1)}(t^n)]$ and $U_h^{(2)}(t^n) = \text{Proj}[U_h^n; I_j, \tilde{I}_j^{(2)}(t^n)]$, which are the $L^2$ projections of solutions from the background mesh to the upstream mesh.

2. **Update** (3.17) from $U_h^n$ to $U_h^{n+1}$, component-by-component.

(a) Get the mesh information of the dynamic element $\tilde{I}_j^{(1)}(t^{(l)})$, $l = 0, ..., s$ on RK stages by $\tilde{x}_{j\pm\frac{1}{2}}^{(1)}(t) = x_{j\pm\frac{1}{2}} + (t - t^{n+1})\nu_{j\pm\frac{1}{2}}^{(1)}$. Here $s = 1$ for forward-Euler method and $s = 2$ for Heun's method (SSPRK2), see the blue domain in Figure 3.2 and for explicit midpoint RK2 with intermediate stage in Figure 3.3.
(b1) For forward-Euler method, compute

$$
\begin{align}
\int_{I_j} u_h^{1,n+1} \Psi(x) dx &= \int_{I_j^{*,(1)}} r_{11} l^{(1)T} U_h^n \psi^{(1),n} dx + \Delta t L_{11}(U_h^n, t^n, \tilde{I}_j^{(1)}(t^n)) \nonumber \\
&\quad + \int_{I_j^{*,(2)}} r_{12} l^{(2)T} U_h^n \psi^{(2),n} dx + \Delta t L_{12}(U_h^n, t^n, \tilde{I}_j^{(2)}(t^n)) \nonumber \\
&= \int_{I_j^{*,(1)}} r_{11} l^{(1)T} U_h^n \psi^{(1),n} dx + \Delta t \left( -\widehat{f_{j+\frac{1}{2}}^{11}}(U_h^{(1)}(t^n)) \psi_{j+\frac{1}{2}}^{(1),n,-} + \widehat{f_{j-\frac{1}{2}}^{11}}(U_h^{(1)}(t^n)) \psi_{j-\frac{1}{2}}^{(1),n,+} \right. \nonumber \\
&\quad \left. + \int_{I_j^{*,(1)}} f_{11}^{11}(U_h^{(1)}(t^n)) \psi_x^{(1)}(x, t^n) dx + \int_{I_j^{*,(1)}} (r_{11} l^{(1)T})_x A U_h^n \psi^{(1)}(x, t^n) dx \right. \nonumber \\
&\quad \left. + \int_{I_j^{*,(1)}} r_{11} l^{(1)T} F(x, t^n) \psi^{(1)}(x, t^n) dx \right) \nonumber \\
&\quad + \int_{I_j^{*,(2)}} r_{12} l^{(2)T} U_h^n \psi^{(2),n} dx + \Delta t \left( -\widehat{f_{j+\frac{1}{2}}^{12}}(U_h^{(2)}(t^n)) \psi_{j+\frac{1}{2}}^{(2),n,-} + \widehat{f_{j-\frac{1}{2}}^{12}}(U_h^{(2)}(t^n)) \psi_{j-\frac{1}{2}}^{(2),n,+} \right) \nonumber \\
&\quad + \int_{I_j^{*,(2)}} f_{12}^{12}(U_h^{(2)}(t^n)) \psi_x^{(2)}(x, t^n) dx + \int_{I_j^{*,(2)}} (r_{12} l^{(2)T})_x A U_h^n \psi^{(2)}(x, t^n) dx \nonumber \\
&\quad \left. + \int_{I_j^{*,(2)}} r_{12} l^{(2)T} F(x, t^n) \psi^{(2)}(x, t^n) dx \right) \tag{3.20}
\end{align}
$$

where $\Delta t^n = t^{n+1} - t^n$, $U_h^{(1)}(t^n) = \text{Proj}[U_h^n, \{I_j\}_{j=1}^N, \{I_j^{*,(1)}\}_{j=1}^N]$ and $U_h^{(2)}(t^n) = \text{Proj}[U_h^n, \{I_j\}_{j=1}^N, \{I_j^{*,(2)}\}_{j=1}^N]$. We compute the four integration terms of (3.20) $\int_{I_j^{*,(1)}} r_{11} l^{(1)T} U_h^n \psi^{(1),n} dx$, $\int_{I_j^{*,(2)}} r_{12} l^{(2)T} U_h^n \psi^{(2),n} dx$, $\int_{\tilde{I}_j^{(1)}} (r_{11} l^{(1)T})_x A U_h^n \psi^{(1)}(x, t^n) dx$ and $\int_{\tilde{I}_j^{(2)}} (r_{12} l^{(2)T})_x A U_h^n \psi^{(2)}(x, t^n) dx$ highlighted in blue with $U_h^n$ on background meshes in order to make the scheme satisfy conservation. Thus they have to be evaluated subinterval-by-subinterval since DG solution is discontinuous. In summary, we can get $u_h^{1,n+1}$ on $I_j$ by (3.20). Similarly, we can get $u_h^{2,n+1}$ on $I_j$.

$$
\begin{array}{c|cc}
0 & 1 & 1 \\
\hline
1 & \frac{1}{2} & \frac{1}{2}
\end{array}
\quad
\begin{array}{c|cc}
0 & 1 & 1 \\
\hline
\frac{1}{2} & 0 & \frac{1}{2}
\end{array}
$$

<caption>Table 3.1: Heun's method (SSPRK2) Butcher Table 3.2: Explicit midpoint RK2 Butcher Tableau Tableau</caption>

(b2) For SSPRK2 method with Butcher tableau: 3.1, we get $u_h^1(t^{(1)})$ from (3.20), then compute

$$
\begin{align}
\int_{I_j} u_h^{1,n+1} \Psi(x) dx &= \int_{I_j^{*,(1)}} r_{11} l^{(1)T} U_h^n \psi^{(1),n} dx + 0.5 \Delta t L_{11}(U_h(t^{(1)}), t^{(1)}, I_j) \nonumber \\
&\quad + \int_{I_j^{*,(2)}} r_{12} l^{(2)T} U_h^n \psi^{(2),n} dx + 0.5 \Delta t L_{12}(U_h(t^{(1)}), t^{(1)}, I_j), \tag{3.21}
\end{align}
$$

where $t^{(1)} = t^{n+1}$, $u_h^1(t^{(1)})$ and $u_h^2(t^{(1)})$ are defined on background mesh $I_j$.
(b3) For general RK methods with intermediate stages, we will update intermediate RK solutions on background mesh as in [9]. For example, for a 2nd order RK method with Butcher tableau 3.2. It has an intermediate stage at $t^{(1)} = t^n + \frac{\Delta t}{2}$, we propose the following steps as [9], see Figure 3.3.

i. We denote the dynamic domain tracking $I_j$ from $t^{(1)}$ to $t_n$ with speed $\nu_{j \pm \frac{1}{2}}^{(1)}$ at mesh point $x_{j \pm \frac{1}{2}}$ as $\tilde{I}_{j,(1)}^{(1)}(t)$, see the green domain in Figure 3.3. $\tilde{I}_{j,(1)}^{(2)}(t)$ related to second characteristic is defined similarly. Then we can update $U_h^{(1)}$ on $\tilde{I}_{j,(1)}^{(1)}(t)$ and $\tilde{I}_{j,(1)}^{(2)}(t)$ from $t^n$ to $t^{(1)}$ as in a forward-Euler method.

ii. We update $U_h^{n+1}$ on dynamic domain $\tilde{I}_j^{(1)}(t)$ and $\tilde{I}_j^{(2)}(t)$ from $U_h^n$ with projection onto $I_j^{*,(1)}$ and $U_h(t^{(1)})$ with projection onto $\tilde{I}_{j,(1)}^{(1)}(t^{(1)})$.

<figure><img src="image_4.png" /><figcaption>Figure 3.3: Update RK intermediate solution at background mesh (red line) from the first characteristic family of a hyperbolic system.</figcaption></figure>

**Theorem 3.2. (Conservation)** The proposed fully discrete ELDG scheme with strong stability preserving Runge-Kutta (SSPRK) time discretization for (3.1) with $F = 0$ is locally conservative. In particular, given a DG solution $u_h(x, t^n) \in V_h^k$ with a periodic boundary condition, we have

$$ \sum_{i=1}^{N} \int_{I_j} U_h(x, t^{n+1}) dx = \sum_{i=1}^{N} \int_{I_j} U_h(x, t^n) dx. $$

*Proof.* We firstly consider the forward Euler time discretization. Taking $\Psi = 1$ and $F = 0$ in the
scheme (3.20), we have

$$
\begin{align}
\int_{\Omega} u_h^{1,n+1} dx &= \sum_j \int_{I_j} u_h^{1,n+1} dx \nonumber \\
&= \sum_j \left[ \int_{I_j^{*,(1)}} r_{11} l^{(1)T} U_h^n dx + \Delta t \left( -\widehat{f_{j+\frac{1}{2}}^{11}}(U_h^{(1)}(t^n)) + \widehat{f_{j-\frac{1}{2}}^{11}}(U_h^{(1)}(t^n)) + \int_{I_j^{*,(1)}} (r_{11} l^{(1)T})_x A U_h^n dx \right) \right] \nonumber \\
&\quad + \sum_j \left[ \int_{I_j^{*,(2)}} r_{12} l^{(2)T} U_h^n dx + \Delta t \left( -\widehat{f_{j+\frac{1}{2}}^{12}}(U_h^{(2)}(t^n)) + \widehat{f_{j-\frac{1}{2}}^{12}}(U_h^{(2)}(t^n)) + \int_{I_j^{*,(2)}} (r_{12} l^{(2)T})_x A U_h^n dx \right) \right] \nonumber \\
&= \sum_j \left( \int_{I_j^{*,(1)}} r_{11} l^{(1)T} U_h^n dx + \int_{I_j^{*,(2)}} r_{12} l^{(2)T} U_h^n dx \right) \nonumber \\
&\quad + \Delta t \sum_j \left( \int_{I_j^{*,(1)}} (r_{11} l^{(1)T})_x A U_h^n dx + \int_{I_j^{*,(2)}} (r_{12} l^{(2)T})_x A U_h^n dx \right) \nonumber \\
&= \int_{\Omega} (r_{11} l^{(1)T} + r_{12} l^{(2)T}) U_h^n dx + \Delta t \int_{\Omega} (r_{11} l^{(1)T} + r_{12} l^{(2)T})_x A U_h^n dx \nonumber \\
&= \int_{\Omega} u_h^{1,n} dx, \tag{3.22}
\end{align}
$$

which follows from the cancellation of unique fluxes at cell boundaries, $r_{11}l^{(1)T} + r_{12}l^{(1)T} = [1,0]$ and (3.8) with integration in a subinterval-by-subinterval fashion. The conservation for the fully discrete ELDG scheme can be proved in a similar fashion.

**Remark 3.3.** To maintain the mass-conservative property, the choice of eigenvectors $R(x)$ is not necessarily exact for ELDG scheme, as long as $R(x)$ and $R^{-1}(x)$ are a consistent pair throughout the domain. We can also choose an approximation of the exact eigenvector if it is not easy to obtain.

## 3.2 2D linear hyperbolic system

The solution for high-dimensional hyperbolic systems is given by means of a characteristic cone, rather than individual characteristic lines [16]. Numerically, characteristic Galerkin [26] or evolution Galerkin [2, 18, 19] methods have been proposed and developed to solve high dimensional hyperbolic system. This method is constructed by taking into account information propagated in all bicharacteristic directions and involving integrals around the characteristic cone. However, the awkward integrals over the mantle, involving intermediate times, limit both the accuracy and the stability of the resulting schemes. Thus the finite volume evolution Galerkin (FVEG) schemes are introduced, which is in a predictor-corrector plus finite volume framework to get higher accuracy.
FVEG method has been developed in [24] and widely applied in [17, 24, 20], and the stability and accuracy have been investigated in [21, 24, 23, 22]. Even though the FVEG method can achieve high-order accuracy and stability with extra large step, the algorithm implementation is still very complex for high-dimensional problems. In this paper, we use the dimensional splitting method for higher dimensional problem.

Consider a first order 2D linear hyperbolic system

$$U_t + (A(x, y)U)_x + (B(x, y)U)_y = 0, \quad (x, y) \in \Omega. \qquad (3.23)$$

We assume that the computational domain $\Omega$ is rectangular, and it can be partitioned into rectangular meshes with each computational cell $\Omega_{ij} = [x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}] \times [y_{j-\frac{1}{2}}, y_{j+\frac{1}{2}}]$, where we use the piecewise $Q^k$ tensor-product polynomial spaces. Then we extend ELDG algorithm to 2D problems via dimensional splitting [27].

<figure><img src="image_5.png" /><figcaption>Figure 3.4: Illustration of the 2D ELDG scheme via Strang splitting. k = 3.</figcaption></figure>

1. We first locate $(k+1)^2$ tensor-product Gaussian nodes on cell $\Omega_{ij} : (x_{i,p}, y_{j,q})$, $p, q = 0, ..., k$. For example, see Figure 3.4 (left) for the case of $k = 3$.

2. Then, the equation (3.23) is split into two 1D hyperbolic problems based on the quadrature nodes in both $x-$ and $y-$ directions:

$$U_t + (A(x, y)U)_x = 0, \qquad (3.24)$$

$$U_t + (B(x, y)U)_y = 0. \qquad (3.25)$$

Based on a 1D ELDG formulation, the split equations (3.24) and (3.25) are evolved via Strang splitting over a time step $\Delta t$ as follows.
* Evolve 1D equation (3.24) at different $y_{j,q}$ points for a half time-step $\Delta t/2$, see Figure 3.4 (middle). For each $y_{j,q}$, the $(k+1)$ point values are mapped to a $P^k$ polynomial per cell, then the 1D system (3.24) is evolved by the proposed ELDG scheme. Finally, we can map the evolved $P^k$ polynomial back to the $(k+1)$ point values to update the solution.

* Evolve 1D system (3.25) at different $x_{i,p}$ points for a full time-step $\Delta t$ as above, see Figure 3.4 (right).

* Evolve 1D system (3.24) at different $y_{j,q}$ points for another half time-step $\Delta t/2$.

The splitting 2D ELDG formulation maintains high order accuracy in space, extra large time stepping size with stability and conservation; and has a second order splitting error.

# 4 Numerical results

In this section, we show numerical results of the proposed scheme for several linear strict hyperbolic systems including wave equation, Maxwell equation and linearized shallow water equation. We set the time stepping size as $\Delta t = \frac{CFL}{a}\Delta x$ for 1D and $\Delta t = \frac{CFL}{\frac{a}{\Delta x} + \frac{b}{\Delta y}}$ for 2D, where a and b are maximum eigenvalues of coefficient matrixes in x- and y-directions respectively. We use the classical Runge-Kutta fourth order (RK4) method [14] for time discretization. We study the following aspects: the spatial order of convergence by using small enough time stepping sizes, the spatial super-convergence of the post-processed solutions [1, 7] produced by convolving the ELDG solution with a suitable kernel consisting of B-splines for the purpose of improving spatial convergence rate, the temporal order of convergence and numerical stability under a large time stepping size by varying CFL for a fixed spatial mesh.

## 4.1 1D wave equations

In this part, we consider the 1D wave equation:

$$u_{tt} = (a^2(x)u_x)_x + f(x, t). \qquad (4.1)$$

For simplicity, we assume periodic boundary conditions, and the velocity field $a(x)$ is a continuous non-zero and periodic function of space.
Defining $u^1 = u_t$ and $u^2 = u_x$, we can rewrite (4.1) as a linear system (3.1) with

$$U = \begin{bmatrix} u^1 \\ u^2 \end{bmatrix}, \quad A(x) = \begin{bmatrix} 0 & -a^2(x) \\ -1 & 0 \end{bmatrix}, \quad F(x,t) = \begin{bmatrix} f(x,t) \\ 0 \end{bmatrix}.$$

We also have the corresponding notations for the eigen-decomposition of $A(x)$:

* eigenvalue: $\lambda^{(1)}(x)$, $\lambda^{(2)}(x)$, where $\lambda^{(1)}(x) = a(x)$, $\lambda^{(2)}(x) = -a(x)$ for the 1D wave equation.

* $A(x) = R(x)\Lambda R^{-1}(x)$, where $\Lambda(x) = diag(\lambda^{(1)}(x), \lambda^{(2)}(x))$,

$$R(x) = \begin{bmatrix} r_{11}(x) & r_{12}(x) \\ r_{21}(x) & r_{22}(x) \end{bmatrix} = \begin{bmatrix} -a(x) & a(x) \\ 1 & 1 \end{bmatrix} \quad (4.2)$$

contains the right eigenvectors, and

$$R^{-1}(x) = \begin{bmatrix} l_j^{(1)T} \\ l_j^{(2)T} \end{bmatrix} = \begin{bmatrix} \frac{-1}{2a(x)} & \frac{1}{2} \\ \frac{1}{2a(x)} & \frac{1}{2} \end{bmatrix}. \quad (4.3)$$

contains the left eigenvectors.

Corresponding to (6.1) and (6.2) related to non-conservative EL DG scheme, we also have the approximating constant eigen-decomposition

$$R_j = \begin{bmatrix} r_j^{11} & r_j^{12} \\ r_j^{21} & r_j^{22} \end{bmatrix} = \begin{bmatrix} -a_j & a_j \\ 1 & 1 \end{bmatrix} \quad (4.4)$$

and

$$R_j^{-1} = \begin{bmatrix} l_j^{(1)T} \\ l_j^{(2)T} \end{bmatrix} = \begin{bmatrix} \frac{-1}{2a_j} & \frac{1}{2} \\ \frac{1}{2a_j} & \frac{1}{2} \end{bmatrix}. \quad (4.5)$$

**Example 4.1.** (1D wave equation with constant coefficient.) We consider the 1D wave equation eq. (4.1) with constant coefficient $a(x) = 1$ and the source term $f(x,t) = 0$. The initial data is $u(x, 0) = \sin(x)$, $x \in [0, 2\pi]$ with periodic boundary condition. The exact solution is $u(x, t) = \sin(x + t)$. For the constant coefficient problem, if using exact characteristic velocity fields for space-time partition and exact eigenvectors, the proposed EL DG method is the same as SL DG, then it is unconditionally stable. Here we perturb the characteristic velocity $\nu_{j+\frac{1}{2}}^{(1)}$ in (3.5) at cell boundaries and/or $a(x)$ in (4.2) and (4.3) related to approximating eigenvectors to get ELDG, ELDG1, ELDG2 and ELDG3 schemes respectively. Similarly we implement the non-conservative
<table><caption>Table 4.3: The numerical parameters of ELDG, ELDG1, ELDG2, ELDG3 methods and NMC ELDG, NMC ELDG1, NMC ELDG2, NMC ELDG3 methods for $u_{tt} = u_{xx}$.</caption><thead><tr><th>setting</th><th>ELDG</th><th>ELDG1</th><th>ELDG2</th><th>ELDG3</th></tr></thead><tbody><tr><td>ν<sup>(1)</sup><sub>j+1/2</sub> in (3.5), ν<sup>(2)</sup><sub>j+1/2</sub> = -ν<sup>(1)</sup><sub>j+1/2</sub></td><td>1</td><td>1 + sin(x<sub>j</sub>)Δx</td><td>1</td><td>1 + sin(x<sub>j</sub>)Δx</td></tr><tr><td>a(x) in (4.2) and (4.3)</td><td>1</td><td>1</td><td>1 + sin(x)Δx</td><td>1 + sin(x)Δx</td></tr><tr><td>setting</td><td>NMC ELDG</td><td>NMC ELDG1</td><td>NMC ELDG2</td><td>NMC ELDG3</td></tr><tr><td>ν<sup>(1)</sup><sub>j+1/2</sub> in (3.5), ν<sup>(2)</sup><sub>j+1/2</sub> = -ν<sup>(1)</sup><sub>j+1/2</sub></td><td>1</td><td>1 + sin(x<sub>j</sub>)Δx</td><td>1</td><td>1 + sin(x<sub>j</sub>)Δx</td></tr><tr><td>a<sub>j</sub> in (4.4) and (4.5)</td><td>1</td><td>1</td><td>1 + sin(x<sub>j</sub>)Δx</td><td>1 + sin(x<sub>j</sub>)Δx</td></tr></tbody></table>

ELDG methods denoted as NMC ELDG, NMC ELDG1, NMC ELDG2 and NMC ELDG3. Related parameters of these ELDG methods are given in Table 4.3.

Table 4.4 and 4.5 report spatial accuracies of the ELDG, ELDG1, ELDG2 and ELDG3 methods for this example under the same time stepping size without and with post-processing technique [1, 7]. We can observe the optimal convergence rate $k+1$ and $2k+1$. We vary time stepping sizes, with fixed well-resolved spatial meshes, and plot error vs. CFL in Figure 4.5 and 4.6 for ELDG, ELDG1, ELDG2 and ELDG3 schemes without and with post-processed technique respectively, after a long time $T = 100$. The plots from post-processed ELDG schemes better show the fourth order temporal convergence. ELDG2 and ELDG3 perform comparably; they have a more restricted time step constraint than ELDG1. It indicates that, stability is affected by approximations of characteristic via the space-time partition and approximation of eigenvectors. We also note that, in both Figure 4.5 and 4.6, the CFL allowed with stability is much larger than that of the RK DG method which is $\frac{1}{2k+1}$. Further, we verify the conservative property of the ELDG schemes are around machine precision and the non-conservative property of the NMC ELDG schemes is presented in Figure 4.7.

<table><caption>Table 4.4: 1D wave equation with constant coefficient. $u_{tt} = u_{xx}$ with initial condition $u(x, 0) = \sin(x)$ at $T = 1$. We use CFL = 0.3 and CFL = 0.18 with RK4 time discretization for all $P^1$ and $P^2$ respectively. The error for only $u^1 = u_t$ was shown in this table.</caption><thead><tr><th>Mesh</th><th>L<sup>1</sup> error</th><th>Order</th><th>L<sup>1</sup> error</th><th>Order</th><th>L<sup>1</sup> error</th><th>Order</th><th>L<sup>1</sup> error</th><th>Order</th></tr></thead><tbody><tr><td></td><td colspan="2"><i>P</i><sup>1</sup> ELDG</td><td colspan="2"><i>P</i><sup>1</sup> ELDG1</td><td colspan="2"><i>P</i><sup>1</sup> ELDG2</td><td colspan="2"><i>P</i><sup>1</sup> ELDG3</td></tr><tr><td>20</td><td>2.54E-03</td><td>-</td><td>2.42E-03</td><td>-</td><td>2.55E-03</td><td>-</td><td>2.49E-03</td><td>-</td></tr><tr><td>40</td><td>6.18E-04</td><td>2.03</td><td>5.97E-04</td><td>2.02</td><td>6.18E-04</td><td>2.04</td><td>5.99E-04</td><td>2.06</td></tr><tr><td>80</td><td>1.58E-04</td><td>1.96</td><td>1.55E-04</td><td>1.94</td><td>1.58E-04</td><td>1.96</td><td>1.55E-04</td><td>1.95</td></tr><tr><td>160</td><td>3.66E-05</td><td>2.11</td><td>3.62E-05</td><td>2.10</td><td>3.66E-05</td><td>2.11</td><td>3.62E-05</td><td>2.10</td></tr><tr><td></td><td colspan="2"><i>P</i><sup>2</sup> ELDG</td><td colspan="2"><i>P</i><sup>2</sup> ELDG1</td><td colspan="2"><i>P</i><sup>2</sup> ELDG2</td><td colspan="2"><i>P</i><sup>2</sup> ELDG3</td></tr><tr><td>20</td><td>5.92E-05</td><td>-</td><td>6.91E-05</td><td>-</td><td>6.01E-05</td><td>-</td><td>7.02E-05</td><td>-</td></tr><tr><td>40</td><td>7.48E-06</td><td>-</td><td>7.83E-06</td><td>3.14</td><td>7.49E-06</td><td>3.00</td><td>7.81E-06</td><td>3.17</td></tr><tr><td>80</td><td>9.17E-07</td><td>3.03</td><td>9.29E-07</td><td>3.08</td><td>9.17E-07</td><td>3.03</td><td>9.29E-07</td><td>3.07</td></tr><tr><td>160</td><td>1.17E-07</td><td>2.97</td><td>1.18E-07</td><td>2.98</td><td>1.17E-07</td><td>2.97</td><td>1.18E-07</td><td>2.98</td></tr></tbody></table>
<table><caption>Table 4.5: 1D wave equation with constant coefficient. $u_{tt} = u_{xx}$ with initial condition $u(x, 0) = \sin(x)$ at $T = 1$. We use $CFL = 0.3$ and $CFL = 0.18$ with RK4 time discretization and post-processed technique for all $P^1$ and $P^2$ respectively. The error with post-processed technique for only $u^1 = u_t$ was shown in this table.</caption><thead><tr><th rowspan="2">Mesh</th><th colspan="2">L<sup>1</sup> error</th><th colspan="2">L<sup>1</sup> error</th><th colspan="2">L<sup>1</sup> error</th></tr><tr><th>ELDG</th><th>Order</th><th>ELDG1</th><th>Order</th><th>ELDG2</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>2.26E-04</td><td>-</td><td>2.49E-04</td><td>-</td><td>2.38E-04</td><td>-</td></tr><tr><td>40</td><td>2.36E-05</td><td>3.26</td><td>2.40E-05</td><td>3.38</td><td>2.40E-05</td><td>3.31</td></tr><tr><td>80</td><td>2.66E-06</td><td>3.15</td><td>2.64E-06</td><td>3.18</td><td>2.67E-06</td><td>3.16</td></tr><tr><td>160</td><td>3.15E-07</td><td>3.08</td><td>3.11E-07</td><td>3.08</td><td>3.15E-07</td><td>3.09</td></tr></tbody><thead><tr><th rowspan="2">Mesh</th><th colspan="2">P<sup>2</sup> error</th><th colspan="2">P<sup>2</sup> error</th><th colspan="2">P<sup>2</sup> error</th></tr><tr><th>ELDG</th><th>Order</th><th>ELDG1</th><th>Order</th><th>ELDG2</th><th>Order</th></tr></thead><tbody><tr><td>20</td><td>2.15E-06</td><td>-</td><td>2.27E-06</td><td>-</td><td>2.19E-06</td><td>-</td></tr><tr><td>40</td><td>3.63E-08</td><td>5.89</td><td>3.86E-08</td><td>5.87</td><td>3.67E-08</td><td>5.90</td></tr><tr><td>80</td><td>6.40E-10</td><td>5.83</td><td>6.79E-10</td><td>5.83</td><td>6.46E-10</td><td>5.83</td></tr><tr><td>160</td><td>1.27E-11</td><td>5.66</td><td>1.33E-11</td><td>5.68</td><td>1.28E-11</td><td>5.66</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 4.5: The $L^\infty$ error versus $CFL$ of ELDG methods, ELDG1, ELDG2 and ELDG3 methods for 1D wave equation with constant coefficient: $u_{tt} = u_{xx}$ with initial condition $u(x, 0) = \sin(x)$. A long time simulation is performed with $T = 100$ and mesh size $N = 160$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4.6: The $L^\infty$ error versus $CFL$ of ELDG methods, ELDG1, ELDG2 and ELDG3 methods with post-processed technique for 1D wave equation with constant coefficient: $u_{tt} = u_{xx}$ with initial condition $u(x, 0) = \sin(x)$. A long time simulation is performed with $T = 100$ and mesh size $N = 160$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 4.7: The error of mass versus time of ELDG2, ELDG3, NMC ELDG2 and NMC ELDG3 methods for 1D wave equation with constant coefficient: $u_{tt} = u_{xx}$ with initial condition $u(x, 0) = \sin(x)$. A long time simulation $T = 100$ is performed with meshes $N = 160$, $CFL = 0.1$ and RK4 time discretization.</figcaption></figure>

**Example 4.2.** Then, we test the ELDG schemes for eq. (4.1) with a Gaussian initial condition

$$u^1 = \exp\left(-\frac{x^2}{0.005}\right), u^2 = 0.$$

The computational domain is $[-1, 1]$ with the periodic boundary conditions. The exact solutions $u^1 = 0.5[\exp(-\frac{(x+t)^2}{0.005}) + \exp(-\frac{(x-t)^2}{0.005})]$ and $u^2 = 0.5[\exp(-\frac{(x+t)^2}{0.005}) - \exp(-\frac{(x-t)^2}{0.005})]$ are the superposition of two Gaussian functions with a periodic extensions. We plot $u_1$ ELDG3 with $P^1$ and $P^2$ numerical solutions at time $T = 50.5$ for this system in Figure 4.8. We can observed that there is no significant phase difference with a long time simulation, meanwhile the dissipation can be improved by the mesh refinement and higher order spatial approximation.

<figure><img src="image_9.png" /><figcaption>Figure 4.8: Plots of the exact and numerical solutions $u_1$ at time $T = 50.5$ of ELDG3 scheme for solving $u_{tt} = u_{xx}$ with Gaussian function initial condition. The mesh size of $N = 80$ and $N = 320$ are used. Left: $k = 1$ ELDG3 with $CFL = 1.5$. Right: $k = 2$ ELDG3 with $CFL = 0.9$.</figcaption></figure>

**Example 4.3.** Next, we test the ELDG schemes for eq. (4.1) on $[0, 2\pi]$ with the periodic boundary
conditions and the following discontinuous initial condition

$$ u_0^1(x) = \begin{cases} 1, & \text{if } 0.95\pi \le x \le 1.05\pi, \\ 0.5, & \text{otherwise,} \end{cases} \qquad (4.6) \\ u_0^2(x) = 1. $$

The exact solutions $u^1$ and $u^2$ are discontinuous piecewise constants with moving discontinuities. It is a challenging test for controlling oscillations around discontinuities. We adopt a simple TVD limiter on background mesh at each RK stages with $M = 0$ in [8] for all schemes. As shown in Figure 4.5, the CFL constraint with stability is slightly less than 1 for ELDG3 scheme. We plot the numerical solutions $u_1$ of ELDG3 scheme with $P^1$ and $P^2$, $CFL = 0.9$ in Figure 4.9. It is found that oscillations are well controlled with the TVD limiter and ELDG method performs well for large time stepping size. Moreover, we track the conservation of ELDG methods VS NMC ELDG methods for eq. (4.1) and present results in Figure 4.10. It shows that the ELDG schemes maintain the conservation at the level of machine error, while the NMC ELDG schemes do not.

<figure><img src="image_10.png" /><figcaption>Figure 4.9: Plots of the numerical solutions $u_1$ of ELDG3 scheme with TVD limiter for solving $u_{tt} = u_{xx}$ with step function initial condition. The final integration time $T$ is 2.85. The mesh of 160 is used. Left: $k = 1$ ELDG3+TVDlimiter with $CFL = 1.5$. Right: $k = 2$ ELDG3+TVDlimiter with $CFL = 0.9$.</figcaption></figure>

**Example 4.4.** (1D wave equation with variable coefficient and source term.) We consider the 1D wave equation eq. (4.1) with variable coefficient $a(x) = 2+\sin(x)$ and exact solution $u(x, t) = \sin(x-2t)$ is periodic on $[0, 2\pi]$. The source term is $f(x, t) = -4\sin(x-2t)+\sin(x-2t)(2+\sin(x))^2-2(2+\sin(x))\cos(x)\cos(x-2t)$. For computation, we choose mesh velocity $\nu_{j+\frac{1}{2}} = a(x_{j+\frac{1}{2}})$, $\nu_{j+\frac{1}{2}} = -\nu_{j+\frac{1}{2}}$ and exact eigenvectors with $a(x) = 2+\sin(x)$ in (4.2) and (4.3).

The expected optimal spatial accuracies of the ELDG methods without and with post-processing technique are shown in Table 4.6 and Table 4.7 respectively. In Figure 4.11 and 4.12, we plot the
<figure><img src="image_11.png" /><figcaption>Figure 4.10: The error of mass versus time of ELDG2, ELDG3, NMC ELDG2 and NMC ELDG3 methods with TVD limiter for 1D wave equation with constant coefficient: $u_{tt} = u_{xx}$ with initial condition step function. $T = 2.85$, $N = 160$, CFL = 0.9 and RK4 time discretization are performed for the simulation.</figcaption></figure>

<table><caption>Table 4.6: 1D wave equation with variable coefficient and source term. $u_{tt} = (a^2(x)u_x)_x + f(x, t)$ with initial condition $u(x, 0) = \cos(x)$. $T = 1$. We use CFL = 0.1 for $P^1$ and $P^2$ with RK4 time discretization. The error for only $u^1 = u_t$ was shown in this table.</caption><thead><tr><th>Mesh</th><th>L<sup>1</sup> error</th><th>Order</th><th>L<sup>2</sup> error</th><th>Order</th><th>L<sup>∞</sup> error</th><th>Order</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">P<sup>1</sup></td></tr><tr><td>20</td><td>5.20E-03</td><td>–</td><td>6.70E-03</td><td>–</td><td>2.39E-02</td><td>–</td></tr><tr><td>40</td><td>1.32E-03</td><td>1.98</td><td>1.74E-03</td><td>1.94</td><td>6.54E-03</td><td>1.87</td></tr><tr><td>80</td><td>3.28E-04</td><td>2.00</td><td>4.40E-04</td><td>1.99</td><td>1.68E-03</td><td>1.96</td></tr><tr><td>160</td><td>8.11E-05</td><td>2.02</td><td>1.09E-04</td><td>2.01</td><td>4.16E-04</td><td>2.01</td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">P<sup>2</sup></td></tr><tr><td>20</td><td>1.16E-04</td><td>–</td><td>1.61E-04</td><td>–</td><td>5.48E-04</td><td>–</td></tr><tr><td>40</td><td>1.48E-05</td><td>2.97</td><td>2.01E-05</td><td>3.00</td><td>6.70E-05</td><td>3.03</td></tr><tr><td>80</td><td>1.88E-06</td><td>2.98</td><td>2.46E-06</td><td>3.03</td><td>7.84E-06</td><td>3.10</td></tr><tr><td>160</td><td>2.31E-07</td><td>3.02</td><td>3.13E-07</td><td>2.98</td><td>1.04E-06</td><td>2.91</td></tr></tbody></table>

<table><caption>Table 4.7: 1D wave equation with variable coefficient and source term. $u_{tt} = (a^2(x)u_x)_x + f(x, t)$ with initial condition $u(x, 0) = \cos(x)$. $T = 1$. We use CFL = 0.1 for $P^1$ and $P^2$ with post-processed technique and RK4. The error for only $u^1 = u_t$ was shown in this table.</caption><thead><tr><th>Mesh</th><th>L<sup>1</sup> error</th><th>Order</th><th>L<sup>2</sup> error</th><th>Order</th><th>L<sup>∞</sup> error</th><th>Order</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">P<sup>1</sup></td></tr><tr><td>20</td><td>9.10E-04</td><td>–</td><td>1.07E-03</td><td>–</td><td>1.91E-03</td><td>–</td></tr><tr><td>40</td><td>1.07E-04</td><td>3.09</td><td>1.26E-04</td><td>3.09</td><td>2.26E-04</td><td>3.08</td></tr><tr><td>80</td><td>1.29E-05</td><td>3.05</td><td>1.53E-05</td><td>3.05</td><td>2.75E-05</td><td>3.04</td></tr><tr><td>160</td><td>1.58E-06</td><td>3.02</td><td>1.88E-06</td><td>3.02</td><td>3.39E-06</td><td>3.02</td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">P<sup>2</sup></td></tr><tr><td>20</td><td>5.34E-06</td><td>–</td><td>6.39E-06</td><td>–</td><td>1.56E-05</td><td>–</td></tr><tr><td>40</td><td>8.95E-08</td><td>5.90</td><td>1.03E-07</td><td>5.96</td><td>2.82E-07</td><td>5.79</td></tr><tr><td>80</td><td>1.73E-09</td><td>5.69</td><td>1.94E-09</td><td>5.73</td><td>3.32E-09</td><td>6.41</td></tr><tr><td>160</td><td>6.62E-11</td><td>4.71</td><td>7.59E-11</td><td>4.67</td><td>1.46E-10</td><td>4.51</td></tr></tbody></table>

$L^\infty$ error versus CFL of EL DG methods without and with post-processing technique respectively.

The following observations are made: (1) The high order accuracy of the RK method reduce the
<figure><img src="image_12.png" /><figcaption>Figure 4.11: The $L^\infty$ error versus CFL of ELDG method for $u_{tt} = (a^2(x)u_x)_x + f(x,t)$ with initial condition $u(x,0) = \cos(x)$. $T = 1$. $\Delta t = CFL \Delta x$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 4.12: The $L^\infty$ error versus CFL of ELDG method with post-processed technique for $u_{tt} = (a^2(x)u_x)_x + f(x,t)$ with initial condition $u(x,0) = \cos(x)$. $T = 1$. $\Delta t = CFL \Delta x$.</figcaption></figure>

error magnitude when large time stepping size is used; (2) The ELDG methods with RK4 time discretization perform well around and before $CFL = 1$, which is well above the stability constraint of the RK DG method $1/(2k+1)$ for $P^k$ approximations. (3) After $CFL = 1$ and before stability constraint of the method, the temporal convergence order is observed to be consistent with the order of RK discretization; (4) The EL DG methods with post-processing technique have smaller error magnitude than those without post-processing.

## 4.2 2D wave equations

**Example 4.5.** (Two-dimensional linear system with constant coefficient matrices.) The second order wave equation $u_{tt} = u_{xx} + u_{yy}$ can be written as the following first order linear hyperbolic
system:

$$
\left\{
\begin{array}{l}
\begin{pmatrix} u \\ v \end{pmatrix}_t + \begin{pmatrix} -1 & 0 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} u \\ v \end{pmatrix}_x + \begin{pmatrix} 0 & -1 \\ -1 & 0 \end{pmatrix} \begin{pmatrix} u \\ v \end{pmatrix}_y = \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \\
u(x, y, 0) = \frac{1}{2\sqrt{2}} \sin(x + y) - \frac{1}{2\sqrt{2}} \cos(x + y), \\
v(x, y, 0) = \frac{\sqrt{2}-1}{2\sqrt{2}} \sin(x + y) + \frac{\sqrt{2}+1}{2\sqrt{2}} \cos(x + y)
\end{array}
\right.
\qquad (4.7)
$$

with period boundary conditions in both *x* and *y* directions. The exact solution is

$$
\left\{
\begin{aligned}
u(x, y, t) &= \frac{1}{2\sqrt{2}} \sin(x + y + \sqrt{2}t) - \frac{1}{2\sqrt{2}} \cos(x + y - \sqrt{2}t), \\
v(x, y, t) &= \frac{\sqrt{2}-1}{2\sqrt{2}} \sin(x + y + \sqrt{2}t) + \frac{\sqrt{2}+1}{2\sqrt{2}} \cos(x + y - \sqrt{2}t).
\end{aligned}
\right.
\qquad (4.8)
$$

We notice that the two matrices in equation (4.7) don't commute, thus the linear system can not be reduced to 2-D scalar problems. We test accuracy for $Q^k$ ELDG methods with RK4 and 4th splitting method [4] at $T = 1$ for $k = 1, 2$ with $CFL = 0.1$ in Table 4.8. As expected, the $(k+1)$th order convergence is observed for these methods. We plot the $L^\infty$ error versus $CFL$ of ELDG methods with $Q^1$ (left) and $Q^2$ (right) polynomial spaces for this case with Strang splitting and 4th order splitting in Figure 4.13, which shows that second and fourth order splitting errors are dominant when time-stepping sizes are large enough. The $CFL$ constraint with stability for ELDG method is larger than that for general RK DG method when high order time discretization is applied.

<table><caption>Table 4.8: Two-dimensional linear system with constant coefficient matrices. $Q^k$ EL DG methods ($k = 1, 2$) with RK4 and 4th splitting time discretization methods for (4.7) with the smooth initial condition at $T = 1$. $CFL = 0.1$.</caption>
  <thead>
    <tr>
      <th>Mesh</th>
      <th>L<sup>1</sup> error</th>
      <th>Order</th>
      <th>L<sup>2</sup> error</th>
      <th>Order</th>
      <th>L<sup>∞</sup> error</th>
      <th>Order</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <td colspan="7" style="text-align:center; font-weight:bold;">Q<sup>1</sup></td>
    </tr>
    <tr>
      <td>20<sup>2</sup></td>
      <td>8.03E-04</td>
      <td>-</td>
      <td>9.47E-04</td>
      <td>-</td>
      <td>1.85E-03</td>
      <td>-</td>
    </tr>
    <tr>
      <td>40<sup>2</sup></td>
      <td>2.16E-04</td>
      <td>1.89</td>
      <td>2.50E-04</td>
      <td>1.92</td>
      <td>4.56E-04</td>
      <td>2.02</td>
    </tr>
    <tr>
      <td>80<sup>2</sup></td>
      <td>5.57E-05</td>
      <td>1.96</td>
      <td>6.40E-05</td>
      <td>1.97</td>
      <td>1.13E-04</td>
      <td>2.02</td>
    </tr>
    <tr>
      <td>160<sup>2</sup></td>
      <td>1.43E-05</td>
      <td>1.96</td>
      <td>1.64E-05</td>
      <td>1.97</td>
      <td>2.84E-05</td>
      <td>1.99</td>
    </tr>
    <tr>
      <td colspan="7" style="text-align:center; font-weight:bold;">Q<sup>2</sup></td>
    </tr>
    <tr>
      <td>20<sup>2</sup></td>
      <td>1.70E-04</td>
      <td>-</td>
      <td>1.90E-04</td>
      <td>-</td>
      <td>3.12E-04</td>
      <td>-</td>
    </tr>
    <tr>
      <td>40<sup>2</sup></td>
      <td>2.21E-05</td>
      <td>2.95</td>
      <td>2.47E-05</td>
      <td>2.94</td>
      <td>4.14E-05</td>
      <td>2.91</td>
    </tr>
    <tr>
      <td>80<sup>2</sup></td>
      <td>2.75E-06</td>
      <td>3.00</td>
      <td>3.08E-06</td>
      <td>3.00</td>
      <td>5.21E-06</td>
      <td>2.99</td>
    </tr>
    <tr>
      <td>160<sup>2</sup></td>
      <td>3.38E-07</td>
      <td>3.02</td>
      <td>3.80E-07</td>
      <td>3.02</td>
      <td>6.45E-07</td>
      <td>3.01</td>
    </tr>
  </tbody>
</table>

**Example 4.6.** (Two-dimensional linear system with variable coefficient matrices.) The second order wave equation $u_{tt} = (a^2(x, y)u_x)_x + (b^2(x, y)u_y)_y$ can be written as the following first order linear hyperbolic system by taking $u_1 = u_t$, $u_2 = u_x$, $u_3 = u_y$:

$$
\left\{
\begin{pmatrix} u_1 \\ u_2 \\ u_3 \end{pmatrix}_t + \left[ \begin{pmatrix} 0 & -(a(x,y))^2 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \end{pmatrix} \right]_x + \left[ \begin{pmatrix} 0 & 0 & -(b(x,y))^2 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \end{pmatrix} \begin{pmatrix} u_1 \\ u_2 \\ u_3 \end{pmatrix} \right]_y = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}. \tag{4.9}
$$
<figure><img src="image_14.png" /><figcaption>Figure 4.13: The $L^\infty$ error versus CFL of ELDG method with Strang splitting and 4th splitting, RK4 time discretization for (4.7) at $T = 1$.</figcaption></figure>

Next we consider the system (4.9) with the initial condition

$$ \begin{cases} u_1(x, y, 0) = 2 \cos(x + y), \\ u_2(x, y, 0) = \cos(x + y), \\ u_3(x, y, 0) = \cos(x + y), \end{cases} \qquad (4.10) $$

where $a(x, y) = 1+0.5 \sin(x+y)$, $b(x, y) = \sqrt{(4 - (1+0.5 \sin(x+y))^2)}$ and the boundary condition is periodic in both $x$ and $y$ directions. The exact solution is

$$ \begin{cases} u_1(x, y, t) = 2 \cos(x + y + 2t), \\ u_2(x, y, t) = \cos(x + y + 2t), \\ u_3(x, y, t) = \cos(x + y + 2t). \end{cases} \qquad (4.11) $$

We report the spatial accuracy of $Q^k$ ELDG methods in Table 4.9. The expected optimal convergence is observed. We plot the $L^\infty$ error versus CFL of ELDG methods in Figure 4.14. The ELDG methods perform as well as that for the linear system with constant coefficient matrices, and the CFL allowed with stability is much larger than that of the RK DG method.

## 4.3 2D Maxwell equations

**Example 4.7.** Consider the 2D Maxwell equations:

$$ \begin{cases} \frac{\partial H_x}{\partial t} + \frac{\partial E_z}{\partial y} = 0, \\ \frac{\partial H_y}{\partial t} - \frac{\partial E_z}{\partial x} = 0, \\ \frac{\partial E_z}{\partial t} - \frac{\partial H_y}{\partial x} + \frac{\partial H_x}{\partial y} = 0, \end{cases} \qquad (4.12) $$

which is a linear hyporbolic system and can be written as

$$ U_t + AU_x + BU_y = 0, \qquad (4.13) $$
<table><caption>Table 4.9: Two-dimensional linear system with variable matrices. $Q^k$ ELDG methods ($k = 1, 2$) with RK4 and 4th splitting time discretization methods for (4.9) with the smooth initial condition at $T = 0.1$, $CFL = 0.1$.</caption><thead><tr><th>Mesh</th><th>L<sup>1</sup> error</th><th>Order</th><th>L<sup>2</sup> error</th><th>Order</th><th>L<sup>∞</sup> error</th><th>Order</th></tr></thead><tbody><tr><td colspan="7" style="text-align:center; font-weight:bold;">Q<sup>1</sup></td></tr><tr><td>20<sup>2</sup></td><td>1.89E-03</td><td>-</td><td>2.46E-03</td><td>-</td><td>5.29E-03</td><td>-</td></tr><tr><td>40<sup>2</sup></td><td>4.58E-04</td><td>2.05</td><td>6.01E-04</td><td>2.03</td><td>1.26E-03</td><td>2.07</td></tr><tr><td>80<sup>2</sup></td><td>1.13E-04</td><td>2.02</td><td>1.50E-04</td><td>2.00</td><td>3.10E-04</td><td>2.02</td></tr><tr><td>160<sup>2</sup></td><td>2.81E-05</td><td>2.01</td><td>3.73E-05</td><td>2.01</td><td>7.67E-05</td><td>2.01</td></tr><tr><td colspan="7" style="text-align:center; font-weight:bold;">Q<sup>2</sup></td></tr><tr><td>20<sup>2</sup></td><td>2.95E-04</td><td>-</td><td>3.66E-04</td><td>-</td><td>8.43E-04</td><td>-</td></tr><tr><td>40<sup>2</sup></td><td>4.06E-05</td><td>2.86</td><td>4.91E-05</td><td>2.90</td><td>1.01E-04</td><td>3.07</td></tr><tr><td>80<sup>2</sup></td><td>5.15E-06</td><td>2.98</td><td>6.20E-06</td><td>2.99</td><td>1.30E-05</td><td>2.95</td></tr><tr><td>160<sup>2</sup></td><td>6.49E-07</td><td>2.99</td><td>7.80E-07</td><td>2.99</td><td>1.63E-06</td><td>2.99</td></tr></tbody></table>

<figure><img src="image_15.png" /><figcaption>Figure 4.14: The $L^\infty$ error versus CFL of ELDG method with Strang splitting and 4th splitting, RK4 time discretization for (4.9). $T = 1$, mesh size 40$^2$.</figcaption></figure>

where

$$ U = \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix} = \begin{bmatrix} E_z \\ H_x \\ H_y \end{bmatrix}, \quad A = \begin{bmatrix} 0 & 0 & -1 \\ 0 & 0 & 0 \\ -1 & 0 & 0 \end{bmatrix}, \quad B = \begin{bmatrix} 0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}. $$

We take the computational domain $[-1, 1] \times [-1, 1]$ with periodic boundary condition and the Gaussian function initial condition:

$$ \begin{cases} u_1(x, y, 0) = \exp(-\frac{x^2+y^2}{0.005}), \\ u_2(x, y, 0) = 0, \\ u_3(x, y, 0) = 0. \end{cases} \qquad (4.14) $$

For this example, we show the numerical ELDG $Q^2$ solution $u_1$ at times $T = 0.5, 1, 1.5, 2$ in Figure 4.15.
<figure><img src="image_16.png" /><figcaption>Figure 4.15: Plots of the ELDG numerical solutions $u_1 = E_z$ and their contour plots at $T = 0.5, 1, 1.5, 2$ for 2D Maxwell equation (4.13) with Gaussian function initial condition. The mesh of $80 \times 80$ is used with 4th splitting method and RK4 time discretization, $CFL = 1$.</figcaption></figure>
## 4.4 2D linearized shallow water equations

In this part, we consider the following linearized shallow water system which is a derivation from oceanic shallow water model [11]: Consider the 2D Maxwell equations:

$$ \frac{\partial}{\partial t} \begin{bmatrix} \phi \\ \Phi u \\ \Phi v \end{bmatrix} + \frac{\partial}{\partial x} \begin{bmatrix} \Phi u \\ \Phi \phi \\ 0 \end{bmatrix} + \frac{\partial}{\partial y} \begin{bmatrix} \Phi v \\ 0 \\ \Phi \phi \end{bmatrix} = \begin{bmatrix} 0 \\ f\Phi v - r\Phi u + \frac{\tau_x}{\rho} \\ -f\Phi u - r\Phi v + \frac{\tau_y}{\rho} \end{bmatrix}, \quad (4.15) $$

where $\phi$ is the geopotential height, $\Phi > 0$ is a constant mean flow geopotential height, $(u, v)$ is the perturbed velocity, $\gamma \ge 0$ is the bottom friction, $(\tau_x, \tau_y)$ is the wind stress, $\rho$ is the water density, and $f = f_0 + \beta(y - y_m)$ is the Coriolis parameter, where $f_0, \beta, y_m$ are constants. The linearized shallow water equations is a linear hyporbolic system

$$ U_t + AU_x + BU_y = F \quad \text{in } \Omega, \quad (4.16) $$

where

$$ U = \begin{bmatrix} \phi \\ \Phi u \\ \Phi v \end{bmatrix}, \quad A = \begin{bmatrix} 0 & 1 & 0 \\ \Phi & 0 & 0 \\ 0 & 0 & 0 \end{bmatrix}, \quad B = \begin{bmatrix} 0 & 0 & 1 \\ 0 & 0 & 0 \\ \Phi & 0 & 0 \end{bmatrix}, \quad F = \begin{bmatrix} 0 \\ f\Phi v - \gamma\Phi u + \frac{\tau_x}{\rho} \\ -f\Phi u - \gamma\Phi v + \frac{\tau_y}{\rho} \end{bmatrix}. $$

**Example 4.8.** We take $\Phi = 1$, $f = 0, \gamma = 0$, and $(\tau_x, \tau_y) = 0$, which implies $F = 0$. The computational domain $\Omega$ is taken as $[-1, 1] \times [-1, 1]$ with periodic boundary condition and the discontinuous initial condition:

$$ \left\{ \begin{array}{ll} \phi(x, y, 0) = \begin{cases} 1, & \text{if } y \ge x, \\ 0.5, & \text{otherwise,} \end{cases} & (4.17) \\ u(x, y, 0) = 1, \\ v(x, y, 0) = 1. \end{array} \right. $$

For this example, we also use TVD limiter and show the numerical ELDG $Q^2$ solution $(\phi, u, v)$ at times $T = 0.5$ in Figure 4.16. Here we use $CFL = 1$, and we can clearly observed the good results with this extra lager time step.

**Example 4.9.** We also consider the linear Kelvin wave by taking $\Phi = 1$, $f_0 = y_m = 0, \beta = 1, \gamma = 0$, and $(\tau_x, \tau_y) = 0$, which implies $F = [0, yv, -yu]^T$. The computational domain $\Omega$ is chosen as $[-10, 10] \times [-5, 5]$ with periodic boundary condition. We consider the following exact solution:

$$ \left\{ \begin{array}{l} \phi = 1 + \exp(-\frac{y^2}{2}) \exp(-\frac{(x+5-t)^2}{2}), \\ u = \exp(-\frac{y^2}{2}) \exp(-\frac{(x+5-t)^2}{2}), \\ v = 0. \end{array} \right. \quad (4.18) $$
<figure><img src="image_17.png" /><figcaption>Figure 4.16: Plots of the ELDG $Q^2$ numerical solutions $(\phi, u, v)$ (from top to bottom) and their contour plots at $T = 0.5$ for 2D linearized shallow water equation (4.15) with Gaussian function initial condition. The mesh of 80×80 is used with 4th splitting method and RK4 time discretization, $CFL = 1$.</figcaption></figure>

# 5 Conclusion

In this paper, we have developed a conservative Eulerian-Lagrangian discontinuous Galerkin (ELDG) method for linear hyperbolic systems. The new framework tracks the information of each characteristic family by the corresponding characteristic region, and combines in a conservative fashion. The method inherits advantages in stability under large time stepping sizes, and in conservation, compactness and high order accuracy. These advantages are numerically verified by extensive
numerical tests for 1D and 2D linear wave equations. Future works include further theoretic development and application of ELDG methods for nonlinear hyperbolic problems.

# 6 Appendix

## 6.1 Non-conservative EL DG scheme

In this part, we formulate the scheme by a localized characteristic field. In particular, a piecewise constant $a_j$ approximating $a(x)$ in (3.3) is defined on $I_j$, and the corresponding

$$R_j \doteq [r_j^{(1)} | r_j^{(2)}] = \begin{bmatrix} r_j^{11} & r_j^{12} \\ r_j^{21} & r_j^{22} \end{bmatrix} \qquad (6.1)$$

and

$$R_j^{-1} \doteq \begin{bmatrix} l_j^{(1)T} \\ l_j^{(2)T} \end{bmatrix} \qquad (6.2)$$

Define $l_j^{(1)T}$ is locally defined on $\Omega_j^{(1)}$ approximating $l_j^{(1)T}(x)$. For simplicity, we only present the first order EL DG scheme. Take the vector product of $l_j^{(1)T}$ from left with (3.1), we have a scalar equation

$$l_j^{(1)T}(U_t + (A(x)U)_x) = l_j^{(1)T} F(x, t). \qquad (6.3)$$

Next, integrating over the space-time interval $\tilde{I}_j^{(1)}(t)$, then we have

$$\begin{aligned} & \frac{d}{dt} \int_{\tilde{I}_j^{(1)}(t)} (l_j^{(1)T} U) dx + l_j^{(1)T} (A(x)U - \nu_{j+\frac{1}{2}}^{(1)} U)|_{\tilde{x}_{j+\frac{1}{2}}^{(1)}(t)} - l_j^{(1)T} (A(x)U - \nu_{j-\frac{1}{2}}^{(1)} U)|_{\tilde{x}_{j-\frac{1}{2}}^{(1)}(t)} \\ & = \int_{\tilde{I}_j^{(1)}(t)} l_j^{(1)T} F(x, t) dx. \end{aligned} \qquad (6.4)$$

The first order EL DG discretization of eq. (6.4) is to find $l_j^{(1)} U_h(x, t) \in P^0(\tilde{I}_j^{(1)}(t))$, so that

$$\begin{aligned} & \frac{d}{dt} \int_{\tilde{I}_j^{(1)}(t)} l_j^{(1)T} U_h dx = - \left[ l_j^{(1)T} \widehat{(A(x)U_h - \nu_{j+\frac{1}{2}}^{(1)} U_h)} \right]_{j+\frac{1}{2}} + \left[ l_j^{(1)T} \widehat{(A(x)U_h - \nu_{j-\frac{1}{2}}^{(1)} U_h)} \right]_{j-\frac{1}{2}} \\ & + \int_{\tilde{I}_j^{(1)}(t)} l_j^{(1)T} F(x, t) dx \doteq L_1(U_h(t), t, \tilde{I}_j^{(1)}(t)). \end{aligned} \qquad (6.5)$$

Here $(A(x)U_h - \nu_{j+\frac{1}{2}}^{(1)} U_h)$ at a cell boundary can be taken as a monotone flux, e.g. the Lax-Friedrichs flux

$$(AU - \nu_{j+\frac{1}{2}}^{(1)} U)_{j+\frac{1}{2}} = \frac{1}{2} \left( A(x_{j+\frac{1}{2}}^{+}) U_{j+\frac{1}{2}}^{+} - \nu_{j+\frac{1}{2}}^{(1)} U_{j+\frac{1}{2}}^{+} + A(x_{j+\frac{1}{2}}^{-}) U_{j+\frac{1}{2}}^{-} - \nu_{j+\frac{1}{2}}^{(1)} U_{j+\frac{1}{2}}^{-} - \alpha_{1,2} (U_{j+\frac{1}{2}}^{+} - U_{j+\frac{1}{2}}^{-}) \right)$$
where $\alpha_{1,2} = \max\{|\lambda^{(1)}(x_{j+\frac{1}{2}}) - \nu_{j+\frac{1}{2}}^{(1)}|, |\lambda^{(2)}(x_{j+\frac{1}{2}}) - \nu_{j+\frac{1}{2}}^{(2)}|\}$. Similarly, we can easily update $l_j^{(2)T} U_h$ related to $\lambda^{(2)}$ in the following:

$$ \frac{d}{dt} \int_{\tilde{I}_j^{(2)}(t)} l_j^{(2)T} U_h dx = - \left[ l_j^{(2)T} (A(x)\widehat{U_h - \nu_{j+\frac{1}{2}}^{(2)} U_h})|_{j+\frac{1}{2}} \right] + \left[ l_j^{(2)T} (A(x)\widehat{U_h - \nu_{j-\frac{1}{2}}^{(2)} U_h})|_{j-\frac{1}{2}} \right] \quad (6.6) \\ + \int_{\tilde{I}_j^{(2)}(t)} l_j^{(2)T} F(x,t)dx \doteq L_2(U_h(t), t, \tilde{I}_j^{(2)}(t)), $$

where $l_j^{(2)T}$ is a constant vector, locally defined on $\tilde{I}_j^{(2)}(t)$ approximating $l^{(2)T}(x)$. A simple first order EL DG scheme is composed by two evolution equations (6.5) and (6.6). That is, we can update $u_h^1$ by (6.5), (6.6) and $u_h^1 = (r_j^{11}l_j^{(1)T} + r_j^{12}l_j^{(2)T})U_h$:

$$ \int_{I_j} u_h^{1,n+1} dx = \int_{I_j} r_j^{11} l_j^{(1)T} U_h^{n+1} dx + \int_{I_j} r_j^{12} l_j^{(2)T} U_h^{n+1} dx \\ = r_j^{11} \int_{I_j^{*,(1)}} l_j^{(1)T} U_h^n dx + r_j^{11} \int_{t^n}^{t^{n+1}} L_1(U_h^{(1)}(t), t, \tilde{I}_j^{(1)}(t)) dt \quad (6.7) \\ + r_j^{12} \int_{I_j^{*,(2)}} l_j^{(2)T} U_h^n dx + r_j^{12} \int_{t^n}^{t^{n+1}} L_2(U_h^{(2)}(t), t, \tilde{I}_j^{(2)}(t)) dt, $$

where $U_h^n$ and $U_h^{n+1}$ are defined on the background mesh $I_j$, $U_h^{(1)}(t)$ and $U_h^{(2)}(t)$ are defined on the space-time dynamic meshes $\tilde{I}_j^{(1)}(t)$ and $\tilde{I}_j^{(2)}(t)$ respectively. Similarly, we can update $u_h^2$.

We apply forward-Euler method for time discretization with above EL DG scheme (6.7):

$$ \int_{I_j} u^{1,n+1} dx \\ = r_j^{11} \int_{I_j^{*,(1)}} l_j^{(1)T} U^n dx - \Delta t r_j^{11} l_j^{(1)T} \left[ (A(x)\widehat{U^n - \nu_{j+\frac{1}{2}}^{(1)} U^n})|_{x^{*,(1)}_{j+\frac{1}{2}}} - (A(x)\widehat{U - \nu_{j-\frac{1}{2}}^{(1)} U})|_{x^{*,(1)}_{j-\frac{1}{2}}} \right] \quad (6.8) \\ + r_j^{12} \int_{I_j^{*,(2)}} l_j^{(2)T} U^n dx - \Delta t r_j^{12} l_j^{(2)T} \left[ (A(x)\widehat{U^n - \nu_{j+\frac{1}{2}}^{(2)} U^n})|_{x^{*,(2)}_{j+\frac{1}{2}}} - (A(x)\widehat{U - \nu_{j-\frac{1}{2}}^{(2)} U})|_{x^{*,(2)}_{j-\frac{1}{2}}} \right]. $$

**Remark 6.1.** The above EL DG scheme is not conservative for two reasons:

(1) Flux terms can't cancel each other as $r_j^{11}l_j^{(1)T}$ and $r_j^{12}l_j^{(2)T}$ are discontinuous across cell boundary of $\Omega_j^{(1)}, \Omega_j^{(2)}$.

(2) $\sum_j r_j^{11} \int_{I_j^{*,(1)}} l_j^{(1)T} U^n dx + r_j^{12} \int_{I_j^{*,(2)}} l_j^{(2)T} U^n dx \neq \sum_j \int_{I_j} U^n dx$ because of the inconsistency in characteristic transformations between neighboring cells among $r_j^{11}l_j^{(1)T}$ and $r_j^{12}l_j^{(2)T}$.

## 6.2 The notation of test function

We give some notations which is used in the implementation of the fully discrete EL DG scheme with RK time discretization. For convenience, we only give the definitions related to $\Omega_j^{(1)}$ below
because we can similarly get the definitions related to $\Omega_j^{(2)}$. As scalar case, we take test function $\psi^{(1)}(x,t)$ as $\psi_{j,m}^{(1)}(x,t) = \Psi_{j,m}(x - \alpha^{(1)}(t - t^{n+1}))$ in the adjoint problem which is a set of basis of $P^k(\tilde{I}_j^{(1)}(t))$, where $\tilde{I}_j^{(1)}(t)$ is donated by a domain related $\lambda^{(1)}$ as $\tilde{I}_j(t)$ in scalar problem. Here, we also take $\Psi_{j,m}$ as orthogonal basis on $I_j$, and let

$$u_h^{1,(1)}(x,t) = \sum_{l=0}^{k} \hat{u}_j^{1,(1);(l)}(t) \psi_{j,l}^{(1)}(x,t), \quad \text{on} \quad \tilde{I}_j^{(1)}(t), \qquad (6.9)$$

where $\hat{u}^{1,(1);(l)}$ are coefficients for the basis. Let $\hat{U}_j^{1,(1)}(t) = (\hat{u}_j^{1,(1);(0)}(t), \dots, \hat{u}_j^{1,(1);(k)}(t))^T$ be the coefficient vector of size $(k+1) \times 1$. Then we have

$$\left[ \int_{\tilde{I}_j^{(1)}(t)} u_h^{1,(1)}(x,t) \psi_{j,0}^{(1)}(x,t) dx, \dots, \int_{\tilde{I}_j^{(1)}(t)} u_h^{1,(1)}(x,t) \psi_{j,k}^{(1)}(x,t) dx \right]^T = \hat{U}_j^{1,(1)}(t), \forall t \in [t^n, t^{n+1}].$$

$\hat{U}_j^{2,(1)}(t)$ can be similarly defined. Similar definition can be made to $\Omega_j^{(2)}$ and $\psi^{(2)}(x,t)$ for the second characteristics family.

## References

[1] J. H. Bramble and A. H. Schatz. Higher order local accuracy by averaging in the finite element method. *Mathematics of Computation*, 31(137):94–111, 1977.

[2] D. S. Butler. The numerical solution of hyperbolic systems of partial differential equations in three independent variables. *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences*, 255(1281):232–252, 1960.

[3] X. Cai, W. Guo, and J.-M. Qiu. A high order conservative semi-Lagrangian discontinuous Galerkin method for two-dimensional transport simulations. *Journal of Scientific Computing*, 73(2-3):514–542, 2017.

[4] X. Cai, J. Qiu, and J.-M. Qiu. A conservative semi-Lagrangian HWENO method for the Vlasov equation. *Journal of Computational Physics*, 323:95–114, 2016.

[5] X. Cai, J.-M. Qiu, and Y. Yang. An Eulerian-Lagrangian discontinuous Galerkin method for transport problems and its application to nonlinear dynamics. *Journal of Computational Physics*, 439:110392, 2021.
[6] M. Celia, T. Russell, I. Herrera, and R. Ewing. An Eulerian-Lagrangian localized adjoint method for the advection-diffusion equation. *Advances in water resources*, 13(4):187–206, 1990.

[7] B. Cockburn, M. Luskin, C.-W. Shu, and E. Suli. Enhanced accuracy by post-processing for finite element methods for hyperbolic equations. *Mathematics of Computation*, 72(242):577–606, 2003.

[8] B. Cockburn and C.-W. Shu. TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws. II. General framework. *Mathematics of computation*, 52(186):411–435, 1989.

[9] M. Ding, X. Cai, W. Guo, and J.-M. Qiu. A semi-Lagrangian discontinuous Galerkin (DG)-local DG method for solving convection-diffusion equations. *Journal of Computational Physics*, 409:109295, 2020.

[10] R. S. Falk and G. R. Richter. Explicit finite element methods for symmetric hyperbolic equations. *SIAM Journal on Numerical Analysis*, 36(3):935–952, 1999.

[11] F. X. Giraldo and T. Warburton. A high-order triangular discontinuous galerkin oceanic shallow water model. *International journal for numerical methods in fluids*, 56(7):899–925, 2008.

[12] W. Guo, R. Nair, and J.-M. Qiu. A conservative semi-Lagrangian discontinuous Galerkin scheme on the cubed-sphere. *Monthly Weather Review*, 142(1):457–475, 2013.

[13] X. Hong and Y. Xia. Arbitrary Lagrangian-Eulerian discontinuous Galerkin method for hyperbolic equations involving $\delta$-singularities. *SIAM Journal on Numerical Analysis*, 58(1):125–152, 2020.

[14] G.-S. Jiang and C.-W. Shu. Efficient implementation of weighted ENO schemes. *Journal of Computational Physics*, 126(1):202–228, 1996.

[15] C. Klingenberg, G. Schnücke, and Y. Xia. Arbitrary Lagrangian-Eulerian discontinuous Galerkin method for conservation laws: analysis and application in one dimension. *Mathematics of Computation*, 86(305):1203–1232, 2017.
[16] T. Kröger. *Multidimensional systems of hyperbolic conservation laws, numerical schemes, and characteristic theory: connections, differences, and numerical comparison*. PhD thesis, Aachen, Techn. Hochsch., Diss., 2004, 2004.

[17] T. Kröger and M. Lukáčová-Medvidová. An evolution Galerkin scheme for the shallow water magnetohydrodynamic equations in two space dimensions. *Journal of Computational Physics*, 206(1):122–149, 2005.

[18] P. Lax. Systems of conservation laws. Technical report, LOS ALAMOS NATIONAL LAB NM, 1959.

[19] M. Lukáčová-Medvid'ová, K. Morton, and G. Warnecke. Evolution Galerkin methods for hyperbolic systems in two space dimensions. *Mathematics of Computation*, 69(232):1355–1384, 2000.

[20] M. Lukáčová-Medvid'ová, S. Noelle, and M. Kraft. Well-balanced finite volume evolution Galerkin methods for the shallow water equations. *Journal of computational physics*, 221(1):122–147, 2007.

[21] M. Lukáčová-Medvid'ová, J. Saibertová, and G. Warnecke. Finite volume evolution Galerkin methods for nonlinear hyperbolic systems. *Journal of Computational Physics*, 183(2):533–562, 2002.

[22] M. Lukáčová-Medvid'ová and E. Tadmor. On the entropy stability of the Roe-type finite volume methods. In *Proceedings of the twelfth international conference on hyperbolic problems, American Mathematical Society, eds. J.-G. Liu et al*, volume 67, 2009.

[23] M. Lukáčová-Medviďová, G. Warnecke, and Y. Zahaykah. On the stability of evolution Galerkin schemes applied to a two-dimensional wave equation system. *SIAM journal on numerical analysis*, 44(4):1556–1583, 2006.

[24] M. Lukáčová-Medvid'ová, G. Warnecke, and Y. Zahaykah. Finite volume evolution Galerkin (FVEG) methods for three-dimensional wave equation system. *Applied numerical mathematics*, 57(9):1050–1064, 2007.
[25] P. Monk and G. R. Richter. A discontinuous Galerkin method for linear symmetric hyperbolic systems in inhomogeneous media. *Journal of Scientific Computing*, 22(1):443–477, 2005.

[26] S. Ostkamp. Multidimensional characteristic Galerkin methods for hyperbolic systems. *Mathematical methods in the applied sciences*, 20(13):1111–1125, 1997.

[27] J.-M. Qiu and C.-W. Shu. Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov-Poisson system. *Journal of Computational Physics*, 230(23):8386–8409, 2011.

[28] T. Russell and M. Celia. An overview of research on Eulerian-Lagrangian localized adjoint methods (ELLAM). *Advances in water resources*, 25(8):1215–1231, 2002.

[29] H. Wang, R. Ewing, G. Qin, S. Lyons, M. Al-Lawatia, and S. Man. A family of Eulerian-Lagrangian localized adjoint methods for multi-dimensional advection-reaction equations. *Journal of Computational Physics*, 152(1):120–163, 1999.

[30] K. Wang, H. Wang, and M. Al-Lawatia. An Eulerian-Lagrangian discontinuous Galerkin method for transient advection-diffusion equations. *Numerical Methods for Partial Differential Equations: An International Journal*, 23(6):1343–1367, 2007.