# hp-version discontinuous Galerkin time-stepping schemes for diffusive-viscous wave equation

Min Zhang, Zhaonan Dong, Wenjing Yan

► To cite this version:

Min Zhang, Zhaonan Dong, Wenjing Yan. hp-version discontinuous Galerkin time-stepping schemes for diffusive-viscous wave equation. Computers & Mathematics with Applications, 2025, 200, pp.145–166. 10.1016/j.camwa.2025.09.021. hal-05007733v2

HAL Id: hal-05007733

https://inria.hal.science/hal-05007733v2

Submitted on 3 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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.

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

Distributed under a Creative Commons Attribution 4.0 International License
hp-version discontinuous Galerkin time-stepping schemes for
diffusive-viscous wave equation

Min Zhang*

Zhaonan Dong†

Wenjing Yan‡

September 3, 2025

Abstract

This work introduces a fully discrete scheme for the diffusion-viscous wave equation (DVWe)
in the second-order formulation, combining *hp*-DG time-stepping schemes with conforming finite
element methods (FEM). Two major theoretical contributions are presented: (1) *hp*-version a priori
error estimates in both the energy-norm and DG-norm, which are optimal in the spatial mesh size
$h$, temporal step size $\tau$, and temporal polynomial order $q$, yet suboptimal by one order in the spatial
polynomial order $p$. Furthermore, for solutions exhibiting weak singularities in time, exponential
convergence in terms of the total number of temporal degrees of freedom is proven using the *hp*-
refinement strategy. (2) An energy decay estimate that offers explicit bounds involving the model
parameters, discretization parameters, and the Poincaré inequality constant. A series of numerical
experiments are presented to validate the practical performance of the proposed approach.

**Key words:** diffusive-viscous wave equation, hp-discontinuous Galerkin time-stepping schemes, space-time finite element

# 1 Introduction

The diffusion viscous wave equation (DVWe) has demonstrated a large capacity to characterize frequency-
dependent seismic reflections in hydrocarbon reservoirs, and it has consequently found wide application
in seismic exploration (see [5, 22, 16]). As modern seismic acquisition experiments encounter increasingly
complex geological terrains, simulating seismic wave propagation through subsurface media via numerical
techniques has become a critical challenge in practical exploration.

Despite the practical importance of the diffusion-viscous wave equation (DVWe), relatively few the-
oretical investigations have explored its well-posedness. Under certain positivity assumptions, [14] es-
tablished theoretical results on well-posedness and stability. From a numerical standpoint, most existing
studies rely on finite volume methods [25, 32] and high-order spatial discretizations [23, 15, 36], typi-
cally paired with low-order time discretizations. In particular, the diffusive-viscous term (a third-order
derivative) makes the construction of stable explicit time-stepping schemes particularly difficult, as the
resulting CFL-like condition resemble those in parabolic problems and significantly impede practical
implementation.

These challenges motivate the use of Galerkin-based time discretizations for DVWe. Such approaches
are unconditionally stable and naturally allow for arbitrary-order approximations. A well-known example
is the discontinuous Galerkin (DG) time-stepping scheme, which has been extensively employed—often in
conjunction with finite element (FEM) or DG spatial discretizations—to solve parabolic problems in the
$h$-version framework [19, 30, 7, 24, 12]. More recently, *hp*-version DG methods have been proposed and
analyzed in [27, 26, 10], demonstrating the capability to achieve exponential convergence rates for locally
analytic solutions by appropriately increasing the polynomial order *p* and/or refining the time partition

<footnote>*Department of Mathematics and Statistics, Xi'an Jiaotong University, Xianning West Road 28th, P.R.China, minzcn@qq.com.</footnote>

<footnote>†Inria, 48 rue Barrault, 75013 Paris, France, and CERMICS, Ecole des Ponts, 77455 Marne-la-Vallée 2, France, zhaonan.dong@inria.fr.</footnote>

<footnote>‡Department of Mathematics and Statistics, Xi'an Jiaotong University, Xianning West Road 28th, P.R.China, wenjingyan@mail.xjtu.edu.cn.</footnote>
near singularities. Moreover, DG time-stepping schemes have also been applied to solve second-order wave equations in combination with finite element or DG spatial discretizations [17, 18, 20, 2, 3, 1, 13].

The aim of this work is to propose a fully discrete scheme for the diffusion-viscous wave equation in its second-order formulation, combining hp-DG time-stepping scheme with conforming finite element methods (FEM) in space. The main theoretical contributions of this paper can be summarized as follows: (1). We derive hp-version a priori error bounds in both energy-norm and DG-norm. The error estimates are optimal with respect to the spatial mesh size h, the temporal mesh size τ, and the temporal polynomial order q. They are, however, suboptimal by one order in the spatial polynomial order p. Moreover, for solutions exhibiting weak singularities in time, we establish exponential convergence in terms of the total number of temporal degrees of freedom by using hp-refinement strategy. (2). Building on the approach in [7], we prove the exponential energy decay estimates for the DG time-stepping scheme applied to the DVWe in the diffusive dominant regime. Notably, our results yield explicit energy decay bounds that depend on the diffusive parameter, temporal discretization parameters, and the Poincaré inequality constant. To the best of our knowledge, this appears to be the first time such explicit energy decay estimates have been presented in the literature.

The remainder of this work is structured as follows. Section 2 introduces the diffusion-viscous wave equation model, the discrete framework for space-time discretization, and the analysis tools. In Section 3.4, we analyze the stability and derive a priori error bounds for the semi-discrete scheme in time. We also present the exponential energy decay estimates of the hp-DG time-stepping scheme. Section 4 extends the analysis to the fully discrete scheme and establishes corresponding error estimates for DVWe. Finally, in Section 5, we test the practical performance of the proposed DG time-stepping scheme in combination with FEM spatial discretization through a series of numerical experiments, exploring the convergence rates in space and time and examining energy conservation and dissipation under various model parameters.

## 2 Model problem, discrete setting and analysis tools

We use standard notation throughout for Sobolev and Bochner spaces. Let $v^{(l)}$ denote the $l$-th time partial derivative, with the special notations $v' := v^{(1)}$ and $v'' := v^{(2)}$ for the first- and second-order time derivatives, respectively. Let $\Omega$ be a Lipschitz domain in $\mathbb{R}^d$, $d = 1, 2, 3$, with boundary $\partial\Omega$. For an integer $m \ge 0$, $W^{m,p}(\Omega)$ is the Sobolev space equipped with the responding usual norm $\|\cdot\|_{W^{m,p}(\Omega)}$ and semi-norm $|\cdot|_{W^{m,p}(\Omega)}$. In particular, when $p = 2$, we denote $W^{m,2}(\Omega) = H^m(\Omega)$ and the associated norm $\|\cdot\|_{H^m(\Omega)}$ and semi-norm $|\cdot|_{H^m(\Omega)}$. Moreover, for the Lebesgue space $L^2(\Omega) = H^0(\Omega)$, we denote the norm and the inner product as $\|\cdot\|$ and $(\cdot, \cdot)$ respectively. We here need the subspace $H_0^1(\Omega)$, which is Hilbert space equipped the $H_0^1(\Omega)$ norm and its dual space $H^{-1}(\Omega)$, denote $\langle \cdot, \cdot \rangle$ for the duality pairing in $H^{-1}(\Omega) \times H_0^1(\Omega)$ and dual norm $\|f\|_{H^{-1}(\Omega)} := \sup_{v \in H_0^1(\Omega), \|v\|_{H_0^1(\Omega)}=1} \langle f, v \rangle$. For a Banach space $X$, e.g., $X$ is $H_0^1(\Omega)$, $L^2(\Omega)$ or their subspaces, we introduce some Bochner spaces of $X$-valued functions $L^p(I; X)$ and for integer $m \ge 0$

$$W^{m,p}(I; X) = \{v \in L^p(I; X) \mid \|v^{(l)}\|_{L^p(I; X)} < \infty, \forall l \le m\},$$

with the norm

$$\|v\|_{W^{m,p}(I;X)} = \begin{cases} \left( \int_0^T \sum_{0 \le l \le m} \|v^{(l)}\|_X^p dt \right)^{1/p} & \text{if } 1 \le p < \infty, \\ \max_{0 \le l \le m} \text{ess sup}_{0 \le t \le T} \|v^{(l)}(t)\|_X & \text{if } p = \infty, \end{cases}$$

and semi-norm

$$|v|_{W^{m,p}(I;X)} = \begin{cases} \left( \int_0^T \|v^{(m)}\|_X^p dt \right)^{1/p} & \text{if } 1 \le p < \infty, \\ \text{ess sup}_{0 \le t \le T} \|v^{(m)}(t)\|_X & \text{if } p = \infty. \end{cases}$$

We also give a notation $H^m(I; X) := W^{m,2}(I; X)$ and use its extension for non-integer $m$, which is well-defined by K-method of interpolation. For $a \le cb$ ($a \ge cb$), with constant $c$ is dependent of discrete parameters, we abbreviate it as $a \lesssim b$ ($a \gtrsim b$).
## 2.1 Diffusive-viscous wave equation

Let $\Omega \subset \mathbb{R}^d$ be a bounded polytopic Lipschitz domain in $\mathbb{R}^d$, with $d = 1, 2, 3; T$ a positive final time; we denote $Q = \Omega \times I$ with $I = (0, T]$ as the space-time domain. We recall the classical formulation of the initial boundary value problem of the diffusive-viscous wave equation

$$u'' + \gamma u' - \nabla \cdot (\eta \nabla u') - \nabla \cdot (\zeta^2 \nabla u) = f \quad \text{in } Q, \qquad (2.1)$$

$$u = 0 \quad \text{on } \partial\Omega \times I, \qquad (2.2)$$

$$u = u_0 \quad \text{on } \Omega \times \{t = 0\}, \qquad (2.3)$$

$$u' = u_1 \quad \text{on } \Omega \times \{t = 0\}, \qquad (2.4)$$

where $\gamma = \gamma(\boldsymbol{x})$, $\eta = \eta(\boldsymbol{x})$ and $\zeta = \zeta(\boldsymbol{x})$ are the diffusive attenuation parameter, the viscous attenuation parameter, and the wave propagation speed in the non-dispersive medium, respectively. We assume parameters $\gamma$, $\zeta$ and $\eta$ are non-negative constants and the following regularity assumption holds

$$f \in L^2(I; L^2(\Omega)), \quad u_0 \in H_0^1(\Omega), \quad u_1 \in L^2(\Omega).$$

We consider the following weak formulation of problem: Find a solution $u$ defined on $Q$, such that $u \in L^\infty(I; H_0^1(\Omega))$, $u' \in L^\infty(I; L^2(\Omega)) \cap L^2(I; H_0^1(\Omega))$, $u'' \in L^2(I; H^{-1}(\Omega))$, and

$$\int_0^T \left( \langle u''(t), v \rangle + (\gamma u'(t), v) + (\eta \nabla u'(t), \nabla v) + (\zeta^2 \nabla u(t), \nabla v) \right) dt = \int_0^T (f(t), v) dt, \quad \forall v \in \mathcal{V}, \qquad (2.5)$$

with $u(0) = u_0$, $u'(0) = u_1$ and the space $\mathcal{V} := L^\infty(I; H_0^1(\Omega)) \cap W^{1,\infty}(I; L^2(\Omega)) \cap H^1(I; H_0^1(\Omega))$. The problem (2.5) is well posed: see, e.g., [14].

**Remark 2.1** In [14], the well-posedness is proven for mixed boundary conditions with strictly positive model parameters. To simplify the presentation, we consider only homogeneous Dirichlet boundary conditions and constant model parameters. However, we will numerically explore cases where some of the model parameters are equal to zero.

## 2.2 Discrete setting

We consider a partition $\mathcal{T}_\tau$ of the time interval $\bar{I}$ made by $N_T$ time-intervals such that $0 = t_0 < t_1 < \dots < t_{N_T} = T$ and define $I_n = (t_{n-1}, t_n]$ with time step size $\tau_n = t_n - t_{n-1}$ for all $n = 1, \dots, N_T$. For each time interval $I_n$, we associate a local polynomial degree $q_n$. We collect such time steps and polynomial orders in vectors $\boldsymbol{\tau}$ with $\boldsymbol{\tau}|_{I_n} := \tau_n$ and $\boldsymbol{q}$ with $\boldsymbol{q}|_{I_n} := q_n$, respectively. For any positive integer $k$, $\boldsymbol{q} + k$ is the vector of entries $q_n + k$.

We define the time semi-discrete space $V_{\boldsymbol{\tau}}^{\boldsymbol{q}}(I; X) := \bigoplus_{n=1}^{N_T} V_{\tau_n}^{q_n}(I_n; X)$, with

$$V_{\tau_n}^{q_n}(I_n; X) := \{w \in L^2(I_n; X) : w|_{I_n} \in \mathbb{P}^{q_n}(I_n; X)\} \quad \forall n = 1, \dots, N_T,$$

where $X$ is a Hilbert space. Note that these functions are discontinuous at temporal nodal points. For a function $u \in V_{\boldsymbol{\tau}}^{\boldsymbol{q}}(I; X)$, we use the symbols $u(t_n^\pm) = \lim_{\epsilon \to 0^\pm} u(t_n + \epsilon)$, for $n \ge 1$, and define the jumps across the nodes $t_n$ as $[u]_n = u(t_n^+) - u(t_n^-)$. Moreover, we use the symbol $u_n^+ = u(t_n^+)$ and $u_n^- = u(t_n^-)$ to present the trace of $w$, taken within the interior of $I_{n+1}$ and $I_n$, respectively. For $n = 0$, the jump across the time node $t_0$ is defined as $[u]_0 = u(t_0^+) - u_0$.

Furthermore, We consider a simplicial or tensor-product conforming mesh $\mathcal{T}_h$ of the domain $\Omega$, with mesh size $h_K$ for each $K \in \mathcal{T}_h$ and $h := \max_{K \in \mathcal{T}_h} \{h_K\}$. We assume the existence of a constant $\rho \in (0, 1)$ such that each element $K$ of $\mathcal{T}_h$ is star-shaped with respect to a ball of radius at least $\rho h_K$. Furthermore, we assume that the mesh is globally quasi-uniform, i.e., for the diameters $h_K$ and $h_{K'}$ of any two elements $K$ and $K'$ with $h_K \le h_{K'}$, one has $h_{K'} \le \sigma h_K$ for some strictly positive constant $\sigma$ independent of the mesh. We define the $H^1$-conforming finite element space $V_h^p$ of uniform order $p$, i.e., $V_h^p := \{w \in H_0^1(\Omega) : w|_K \in \mathbb{R}^p(K), \forall K \in \mathcal{T}_h\}$, where $\mathbb{R}^p(K)$ is the polynomial space $\mathbb{P}^p(K)$ for simplicial meshes and $\mathbb{R}^p(K)$ is the tensor-product polynomial space $\mathbb{Q}^p(K)$ for parallelogram and parallelepiped meshes.

In addition, we also define the space-time tensor product space $V_{\boldsymbol{\tau}}^{\boldsymbol{q}}(I; V_h^p) := \bigoplus_{n=1}^{N_T} V_{\tau_n}^{q_n}(I_n; V_h^p)$, with

$$V_{\tau_n}^{q_n}(I_n; V_h^p) := \{w \in L^2(I_n; V_h^p) : w|_{I_n} \in \mathbb{P}^{q_n}(I_n; V_h^p)\} \quad \forall n = 1, \dots, N_T.$$
## 2.3 Analysis tools

Let us briefly review the main hp-analysis tools used in this work. We first introduce the analysis tools related to the time discretization.

**Lemma 2.2** For any polynomial $v \in \mathbb{P}^{q_n}(I_n; X)$, $q_n \ge 0$, it holds

$$ \|v(t_n^-)\|_X^2 \le \frac{(q_n + 1)^2}{\tau_n} \|v\|_{L^2(I_n; X)}^2. \qquad (2.6) $$

**Proof.** The proof of (2.6) can be found in [33]. □

Next, we recall a commonly used interpolation operator for the numerical approximation for wave problem. The construction of the interpolation can be found in [31].

**Definition 2.3** For any $u \in H^2(I_n, X)$, with $n = 1, \dots, N_T$, there exists an interpolation operator $I_{\tau_n}^{q_n} u \in \mathbb{P}^{q_n}(I_n, X)$ such that the following properties hold,

$$ (u - I_{\tau_n}^{q_n} u)(t_n^{\pm}) = 0, \qquad (2.7) $$

$$ (u - I_{\tau_n}^{q_n} u)'(t_n^-) = 0, \qquad (2.8) $$

$$ \int_{I_n} ((u - I_{\tau_n}^{q_n} u)', w) \, dt = 0, \quad \forall w \in \mathbb{P}^{q_n-2}(I_n; X) \quad q_n \ge 2. \qquad (2.9) $$

where $I_{\tau_n}^{q_n} u(t_0^-) := u(t_0^-) = u_0$ and $(I_{\tau_n}^{q_n} u)'(t_0^-) := u'(t_0^-) = v_0$. For $q_n = 1$, the orthogonality property (2.9) is not taken into account. Now let $I = (0, T)$, $u|_{I_n} \in H^2(I_n; X)$ and $I_{\tau_n}^q u \in V_{\tau}^q(I; X)$ be its global projection such that $(I_{\tau_n}^q u)|_{I_n} := I_{\tau_n}^{q_n}(u|_{I_n})$.

**Lemma 2.4** For any $u \in H^{s_{0,n}+1}(I_n; X)$, with $n = 1, \dots, N_T$, any integer $1 \le s_n \le \min\{q_n, s_{0,n}\}$, with $q_n \ge 2$, we have

$$ \begin{aligned} \|u - I^{q_n} u\|_{L^2(I_n; X)}^2 &\lesssim \left(\frac{\tau_n}{2}\right)^{2s_n+2} \frac{(q_n - s_n)!}{q_n^2(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; X)}^2, \\ \|(u - I^{q_n} u)'\|_{L^2(I_n; X)}^2 &\lesssim \left(\frac{\tau_n}{2}\right)^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; X)}^2, \qquad (2.10) \\ \|(u - I^{q_n} u)'(t_{n-1}^+)\|_X^2 &\lesssim \left(\frac{\tau_n}{2}\right)^{2s_n-1} \frac{q_n(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; X)}^2. \end{aligned} $$

**Proof.** The above results can be found in [11, Proposition 2.5.], [28, Theorem 3.10], [6, Lemma 3.3], respectively. See also [34, Lemma 3.5]. □

Next, we introduce the approximation results for the spatial discretisation.

**Definition 2.5** For any $u \in H_0^1(\Omega)$, there exists the L$^2$-projection operator $P_h^p: L^2(\Omega) \to V_h^p$ such that

$$ (P_h^p u, v) = (u, v), \quad \forall v \in V_h^p, \qquad (2.11) $$

Similarly, there exists the elliptic-projection operator $\Pi_h^p: H_0^1(\Omega) \to V_h^p$ such that

$$ (\nabla \Pi_h^p u, \nabla v) = (\nabla u, \nabla v), \quad \forall v \in V_h^p. \qquad (2.12) $$

The following approximation and stability results hold.

**Lemma 2.6** For any $u \in H_0^1(\Omega)$, the following stability bound holds for the L$^2$-projection operator,

$$ \|\nabla P_h^p u\| \lesssim p \|\nabla u\|. \qquad (2.13) $$

In addition, assuming that $u \in H^{r_0+1}(\Omega)$, with $r_0 \ge 1$, following bounds hold:

$$ \|u - P_h^p u\| = \inf_{v \in V_h^p} \|u - v\| \lesssim \frac{h^{r+1}}{p^{r_0+1}} \|u\|_{H^{r_0+1}(\Omega)}, \qquad (2.14) $$
$$ \|\nabla(u - P_h^p u)\| \lesssim \frac{h^r}{p^{r_0-1}} \|u\|_{H^{r_0}(\Omega)}, \qquad (2.15) $$

$$ \|\nabla(u - \Pi_h^p u)\| = \inf_{v \in V_h^p} \|\nabla(u - v)\| \lesssim \frac{h^r}{p^{r_0}} \|u\|_{H^{r_0+1}(\Omega)}, \qquad (2.16) $$

with $r = \min\{p, r_0\}$.

**Proof.** First, we start with (2.13), we split the error bound:

$$ \|\nabla P_h^p u\| \le \|\nabla(P_h^p u - \mathcal{I}_{\mathrm{KM}} u)\| + \|\nabla \mathcal{I}_{\mathrm{KM}} u\|, $$

where $\mathcal{I}_{\mathrm{KM}}$ denotes the $H^1$-conforming Karkulik-Melenk operator. We refer to the original work [21] and the modified version in [9, Corollary 2.5], with the following bounds:

$$ \frac{p}{h} \|\nabla(u - \mathcal{I}_{\mathrm{KM}} u)\| + \|\nabla \mathcal{I}_{\mathrm{KM}} u\| \lesssim \|\nabla u\|. \qquad (2.17) $$

Using the above relation (2.17), $H^1$-norm to $L^2$-norm inverse inequality, $P_h^p \circ \mathcal{I}_{\mathrm{KM}} u = \mathcal{I}_{\mathrm{KM}} u$, and the stability of $P_h^p$, we infer

$$ \begin{aligned} \|\nabla P_h^p u\| &\le \|\nabla(P_h^p u - \mathcal{I}_{\mathrm{KM}} u)\| + \|\nabla \mathcal{I}_{\mathrm{KM}} u\| \\ &\lesssim \frac{p^2}{h} \|P_h^p(u - \mathcal{I}_{\mathrm{KM}} u)\| + \|\nabla u\| \lesssim \frac{p^2}{h} \|u - \mathcal{I}_{\mathrm{KM}} u\| + \|\nabla u\| \lesssim p \|\nabla u\|. \end{aligned} $$

Next, the error bounds (2.14) and (2.16) can be derived by the classical Babuška-Suri operator [4]. The error bound (2.15) can be derived by the stability result (2.13) and the classical Babuška-Suri operator [4]. □

**Remark 2.7 (Quasi-uniform mesh assumption)** We point out that the global quasi-uniform mesh assumption can be relaxed to the meshes with bounded mesh grading. The $H^1$-stability of $L^2$-projection onto FEM space with shape-regular mesh under bounded mesh grading assumption can be found in [8]. Therefore, the above hp-approximation results (2.14), (2.15) and (2.16) also hold on the shape regular meshes with bounded mesh grading.

# 3 Semi-discrete DG time-stepping scheme

In this section we introduce the semi-discrete DG time-stepping scheme for DVWe. The stability and error estimates will be given in DG-norm and energy-norm. Moreover, we present an exponential decay of the initial data energy in time for DVWe under the assumption wave propagation speed parameter $\zeta = 0$.

## 3.1 DG time-stepping scheme

By testing the PDE (2.1) by any function $v'$, $v \in V_{\tau}^q(I; H_0^1(\Omega))$, and integrate it in space and time over the cylinder $Q_n := I_n \times \Omega$,

$$ \int_{I_n} \langle u''(t), v' \rangle + (\gamma u'(t), v') + (\eta \nabla u'(t), \nabla v') + (\zeta^2 \nabla u(t), \nabla v') \, dt = \int_{I_n} (f, v') \, dt. $$

Next, we impose the initial conditions by adding the upwind jump terms ($[u']_{n-1}, v'(t_{n-1}^+)$) and also $(\zeta^2 \nabla [u]_{n-1}, \nabla v(t_{n-1}^+))$. So the discrete formulation on $Q_n$, $n = 1, 2, \dots, N_T$, is defined as: Find $u_\tau \in V_{\tau_n}^{q_n}(I; H_0^1(\Omega))$ such that

$$ B_n(u_\tau, v) = \ell_n(v) \quad v \in V_{\tau_n}^{q_n}(I_n; H_0^1(\Omega)), \qquad (3.1) $$

where the local bilinear form $B_n(u_\tau, v)$ is defined as

$$ \begin{aligned} B_n(u_\tau, v) := & \int_{I_n} (u_\tau'', v') + (\gamma u_\tau', v') + (\eta \nabla u_\tau', \nabla v') \\ & + (\zeta^2 \nabla u_\tau, \nabla v') \, dt + ([u_\tau']_{n-1}, v'(t_{n-1}^+)) + (\zeta^2 \nabla [u_\tau]_{n-1}, \nabla v(t_{n-1}^+)), \end{aligned} \qquad (3.2) $$
and the linear form is defined as $\ell_n(v) := \int_{I_n} (f, v') dt$. For $n = 1$, we define $u'_\tau(t_0^-) := u_1$ and $u_\tau(t_0^-) := u_0$.

Summing up all time intervals, then $u_\tau$ is a solution of the following discrete formulation: Find $u_\tau \in V_\tau^q(I; H_0^1(\Omega))$ such that

$$B(u_\tau, v) = \ell(v) \quad v \in V_\tau^q(I; H_0^1(\Omega)), \qquad (3.3)$$

where

$$B(u, v) := \sum_{n=1}^{N_T} \int_{I_n} (u''_{\tau}, v') + (\gamma u'_{\tau}, v') + (\eta \nabla u'_{\tau}, \nabla v') + (\zeta^2 \nabla u_{\tau}, \nabla v') dt \\ + \sum_{n=2}^{N_T} (([u'_{\tau}]_{n-1}, v'(t_{n-1}^{+})) + (\zeta^2 \nabla [u_{\tau}]_{n-1}, \nabla v(t_{n-1}^{+}))) + (u'_{\tau}(t_0^{+}), v'(t_0^{+})) + (\zeta^2 \nabla u_{\tau}(t_0^{+}), \nabla v(t_0^{+})).$$

and

$$\ell(v) := \sum_{n=1}^{N_T} \int_{I_n} (f, v') dt + (u_1, v'(t_0^{+})) + (\zeta^2 \nabla u_0, \nabla v(t_0^{+})).$$

## 3.2 Stability analysis

For any function $v \in V_{\tau_n}^{q_n}(I_n; H_0^1(\Omega))$, $n = 1, \dots, N_T$, the local stability bound for the local bilinear form (3.2), obtained via integration by parts in time, is defined as:

$$B_n(v, v) = \int_{I_n} \|\gamma^{\frac{1}{2}} v'\|^2 + \|\eta^{\frac{1}{2}} \nabla v'\|^2 dt \\ + \frac{1}{2} (\| [v']_{n-1} \|^2 + \|\zeta \nabla [v]_{n-1} \|^2 + \|v'(t_n^-)\|^2 + \|\zeta \nabla v(t_n^-)\|^2 - \|v'(t_{n-1}^-)\|^2 - \|\zeta \nabla v(t_{n-1}^-)\|^2). \qquad (3.4)$$

Similar, for any function $v \in V_\tau^q(I; H_0^1(\Omega))$, we have the following stability bound,

$$B(v, v) := |||v|||^2, \qquad (3.5)$$

with the DG-norm $||| \cdot |||$ defined as,

$$|||v|||^2 := \|v\|_E^2 + \frac{1}{2}\|v'(t_N^-)\|^2 + \frac{1}{2}\|\zeta \nabla v(t_N^-)\|^2 \\ + \frac{1}{2}\sum_{n=2}^{N_T} (\| [v']_{n-1} \|^2 + \|\zeta \nabla [v]_{n-1} \|^2) + \frac{1}{2}\|v'(t_0^+)\|^2 + \frac{1}{2}\|\zeta \nabla v(t_0^+)\|^2, \qquad (3.6)$$

and we define the energy-norm

$$\|v\|_E^2 := \int_0^T \|\gamma^{\frac{1}{2}} v'\|^2 + \|\eta^{\frac{1}{2}} \nabla v'\|^2 dt. \qquad (3.7)$$

It is easy to check that $||| \cdot |||$ is a norm on $V_\tau^q(I; H_0^1(\Omega))$.

Next, we present the stability bound of the DG solution in the DG-norm.

**Proposition 3.1** Let $\gamma\eta \neq 0$ and $u_\tau$ be the solution of semi-discrete scheme (3.3). It holds

$$|||u_\tau||| \lesssim \|u_1\| + \|\zeta \nabla u_0\| + \left( \min\left\{\frac{1}{\gamma}, \frac{C_p^2}{\eta}\right\} \int_0^T \|f\|^2 ds \right)^{\frac{1}{2}}.$$
**Proof.** Recall (3.3), choosing $v = u_τ$, using the Cauchy-Schwartz inequality, we infer

$$ |||u_{\tau}|||^2 = B(u_{\tau}, u_{\tau}) = \ell(u_{\tau}) \le \int_{0}^{T} \|f\| \|v'\| dt + \|u_1\| \|v'(t_0^+)\| + \|\zeta \nabla u_0\| \|\zeta \nabla v(t_0^+)\|. $$

From

$$ \begin{aligned} \int_{0}^{T} \|f\| \|v'\| \, dt &\le \frac{1}{2\gamma} \int_{0}^{T} \|f\|^2 \, dt + \frac{\gamma}{2} \int_{0}^{T} \|v'\|^2 \, dt, \\ \int_{0}^{T} \|f\| \|v'\| \, dt &\le \frac{C_p^2}{2\eta} \int_{0}^{T} \|f\|^2 \, dt + \frac{\eta}{2C_p^2} \int_{0}^{T} \|v'\|^2 \, dt \le \frac{C_p^2}{2\eta} \int_{0}^{T} \|f\|^2 \, dt + \frac{\eta}{2} \int_{0}^{T} \|\nabla v'\|^2 \, dt, \end{aligned} $$

then the bound of Proposition 3.1 is derived by using the Young's inequality. □

The well-posedness of the proposed DG time-stepping scheme is proven by using the above stability results.

### 3.3 Energy estimate for the DG time-stepping scheme with $\zeta = 0$

In this section, we consider the energy estimate for the diffusive wave equation in the diffusive dominant regime with $\zeta = 0$ and $f = 0$.

**Lemma 3.2** Let $u_τ$ be the solution of semi-discrete scheme (3.3). If $\zeta = 0$, and $f = 0$, we have

$$ \|u'_{\tau}(t_n^-)\|^2 \le \prod_{n=1}^{N_T} (1 + \Lambda_i)^{-1} \|u_1\|^2, \quad \text{with} \quad \Lambda_n = \frac{2\tau_n}{(q_n+1)^2} (\gamma + C_p^{-2}\eta), \quad (3.8) $$

where $C_p$ is the Poincaré inequality constant $\|v\| \le C_p \|\nabla v\|$ for all $v \in H_0^1(\Omega)$.

**Proof.** Recalling the proposed methods (3.1), we have $B_n(u_τ, u_τ) = 0$, for $n = 1, \dots, N_T$. Then, invoking the local discrete stability (3.4) and applying the inverse trace inequality (2.6), we have

$$ \begin{aligned} \|u'_{\tau}(t_{n-1}^-)\|^2 &= \int_{I_n} 2\|\gamma^{\frac{1}{2}} u'_{\tau}\|^2 + 2\|\eta^{\frac{1}{2}} \nabla u'_{\tau}\|^2 dt + \| [u'_{\tau}]_{n-1} \|^2 + \|u'_{\tau}(t_n^-)\|^2 \\ &\ge \int_{I_n} 2(C_p^{-2}\eta + \gamma) \|u'_{\tau}\|^2 dt + \|u'_{\tau}(t_n^-)\|^2 \ge (1 + \Lambda_n) \|u'_{\tau}(t_n^-)\|^2. \end{aligned} $$

Therefore, using the reconstruction relation, we arrive the energy estimates. □

**Remark 3.3** From (3.8), assuming constant time step $\tau_n = \tau$, constant polynomial degree $q_n = q$, for all $n = 1, \dots, N_T$, we deduce

$$ \|u'_{\tau}(t_n^-)\| \le e^{-\frac{\gamma+C_p^{-2}\eta}{(q+1)^2}t_n} \|u_1\|, \quad (3.9) $$

as $\tau \to 0$.

### 3.4 Error estimates for the DG time-stepping scheme

In this subsection we form the abstract error bound first, and then use the interpolation approximation results in the previous subsection to derive hp error bound.

First, we split the error as follows:

$$ e := u - u_{\tau} = (u - I_{\tau}^{q}u) + (I_{\tau}^{q}u - u_{\tau}) = \rho + \theta, \quad (3.10) $$

where the interpolation operator $I_{\tau}^q$ is defined in Definition 2.3.
**Lemma 3.4** Let $u_τ$ and $u$ be the solution of the semi-discrete scheme (3.3) and the problem (2.5), respectively. Then it holds

$$ \|e\|_{\mathrm{E}}^{2} \lesssim \int_{0}^{T} \|\gamma^{\frac{1}{2}} \rho'\|^{2} + \|\eta^{\frac{1}{2}} \nabla \rho'\|^{2} + \frac{\zeta^{4}}{\eta} \|\nabla \rho\|^{2} \, dt, \quad (3.11) $$

and

$$ |||e|||^2 \lesssim \int_0^T \|\gamma^{\frac{1}{2}} \rho'\|^2 + \|\eta^{\frac{1}{2}} \nabla \rho'\|^2 + \frac{\zeta^4}{\eta} \|\nabla \rho\|^2 dt + \sum_{n=0}^{N_T-1} \|\rho'(t_n^+)\|^2. \quad (3.12) $$

**Proof.** Using the consistency of the DG time-stepping scheme, we have $B(e, v) = 0$, for any $v \in V_{\tau}^q(I; H_0^1(\Omega))$. Therefore, choosing $v = \theta$, we obtain the following error equation

$$ |||\theta|||^2 = B(\theta, \theta) = -B(\rho, \theta). \quad (3.13) $$

Next, using integration by parts for the first term in $I_n$, we obtain

$$ \begin{aligned} B(\rho, \theta) &= \sum_{n=1}^{N_T} \int_{I_n} (\rho'', \theta') + (\gamma\rho', \theta') + (\eta\nabla\rho', \nabla\theta') + (\zeta^2\nabla\rho, \nabla\theta') \, dt \\ &\quad + \sum_{n=1}^{N_T-1} ([\rho']_n, v'(t_n^+)) + (\zeta^2\nabla[\rho]_n, \nabla\theta(t_n^+)) + (\rho'(t_0^+), \theta'(t_0^+)) + (\zeta^2\nabla\rho(t_0^+), \nabla\theta(t_0^+)) \\ &= \sum_{n=1}^{N_T} \int_{I_n} -(\rho', \theta'') + (\gamma\rho', \theta') + (\eta\nabla\rho', \nabla\theta') + (\zeta^2\nabla\rho, \nabla\theta') \, dt \\ &\quad - \sum_{n=1}^{N_T-1} (\rho'(t_n^-), [v']_n) + (\rho'(t_{N_T}^-), \theta'(t_{N_T}^-)) + \sum_{n=1}^{N_T-1} (\zeta^2\nabla[\rho]_n, \nabla\theta(t_n^-)) + (\zeta^2\nabla\rho(t_0^+), \nabla\theta(t_0^+)). \end{aligned} \quad (3.14) $$

Using the interpolation properties (2.7)-(2.9), it gives that

$$ B(\rho, \theta) = \sum_{n=1}^{N_T} \int_{I_n} (\gamma\rho', \theta') + (\eta\nabla\rho', \nabla\theta') + (\zeta^2\nabla\rho, \nabla\theta') \, dt. \quad (3.15) $$

Using the Cauchy-Schwartz inequality, we infer

$$ \begin{aligned} |B(\rho, \theta)| &\le \sum_{n=1}^{N_T} \int_{I_n} \|\gamma^{\frac{1}{2}} \rho'\| \|\gamma^{\frac{1}{2}} \theta'\| + \|\eta^{\frac{1}{2}} \nabla \rho'\| \|\eta^{\frac{1}{2}} \nabla \theta'\| + \zeta^2 \|\nabla \rho\| \|\nabla \theta'\| \, dt, \\ &= \sum_{n=1}^{N_T} \int_{I_n} \|\gamma^{\frac{1}{2}} \rho'\| \|\gamma^{\frac{1}{2}} \theta'\| + \|\eta^{\frac{1}{2}} \nabla \rho'\| \|\eta^{\frac{1}{2}} \nabla \theta'\| + \frac{\zeta^2}{\eta^{\frac{1}{2}}} \|\nabla \rho\| \|\eta^{\frac{1}{2}} \nabla \theta'\| \, dt \end{aligned} $$

and

$$ |||\theta|||^2 \le |B(\rho, \theta)| \lesssim \left( \int_0^{t_n} \|\gamma^{\frac{1}{2}} \rho'\|^2 + \|\eta^{\frac{1}{2}} \nabla \rho'\|^2 + \frac{\zeta^4}{\eta} \|\nabla \rho\|^2 \, dt \right)^{\frac{1}{2}} |||\theta|||. $$

By using triangle inequality, we obtain the abstract error bounds. □

**Theorem 3.5** Let $u_τ$ and $u$ be the solution of the semi-discrete scheme (3.3) and the problem (2.1), respectively. Assume that $u|_{I_n} \in H^{k_n+1}(I_n; H_0^1(\Omega))$, for all $n = 1, \dots, N_T$, the following results hold

$$ \|e\|_{\mathrm{E}}^{2} \lesssim \sum_{n=1}^{N_T} \left( \gamma \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2 + \left( \eta + \frac{\tau_n^2 \zeta^4}{q_n^2 \eta} \right) \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; H_0^1(\Omega))}^2 \right), \quad (3.16) $$
and

$$ |||e|||^2 \lesssim \sum_{n=1}^{N_T} \left( \gamma \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2 + \left( \eta + \frac{\tau_n^2 \zeta^4}{q_n^2 \eta} \right) \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; H_0^1(\Omega))}^2 \right. \\ \left. + \tau_n^{2s_n-1} q_n \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2 \right), \quad (3.17) $$

with $1 \le s_n \le \min\{q_n, k_n\}$.

**Proof.** Recalling the Lemma 3.4, the DG-norm error (3.12) is bounded by four terms $T_i$, $i = 1, 2, 3, 4$. We will bound them by using the approximation property (2.10) separately.

For $T_1$ and $T_2$ terms, using second bound in (2.10), we have

$$ T_1 + T_2 \lesssim \sum_{n=1}^{N_T} \left( \gamma \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2 + \eta \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; H_0^1(\Omega))}^2 \right). $$

For the $T_3$ term, using first bound in (2.10), we have

$$ T_3 \lesssim \sum_{n=1}^{N_T} \left( \frac{\zeta^4}{\eta} \frac{\tau_n^2}{q_n^2} \right) \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; H_0^1(\Omega))}^2. $$

Similarly, the $T_4$ term can be bounded by using the third bound in (2.10)

$$ T_4 \lesssim \sum_{n=1}^{N_T} \tau_n^{2s_n-1} q_n \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2. $$

Combining the above bounds, the error estimates (3.16) and (3.17) are derived. $\square$

**Corollary 3.6** If we assume time steps $\tau_n = \tau$, polynomial degrees $q_n = q$, regularity indices $k_n = k$, $s_n = s$, for $n = 1, \dots, N_T$ and using the Stirling formula. The above error bounds read

$$ \|e\|_{E} \approx \mathcal{O}\left(\frac{\tau^{s}}{q^{s}}\right) $$

and

$$ |||e||| \approx \mathcal{O}\left(\frac{\tau^{s-\frac{1}{2}}}{q^{s-\frac{1}{2}}}\right), $$

both of the above two error bounds are optimal time step size $\tau$ and polynomial degree $q$.

## 3.5 Exponential convergence rate for the singular solutions

In this subsection, we consider the hp-DG time-stepping scheme for solutions contains the weak start-up singularities at time $t = 0$, but analytic for $t > 0$. We proof that combining the geometrically graded time meshes and the linearly increasing polynomial degrees on different time interval, yields exponential convergence rates in terms of number of temporal degrees of freedom $N_{DoFs,t}$ for the piecewise analytic solutions exhibiting weak singularities in time.

We assume that the following regularity holds for the solution of (2.5), there are constants $C, d > 0$ such that for any real number $t \in I$, with $I := (0, 1)$ and $X \in \{L^2(\Omega), H_0^1(\Omega)\}$, it holds

$$ \|u^{(s)}(t)\|_{X} \le C d^{s} \Gamma(s+1) t^{\theta-s}, \theta > 1, \quad (3.18) $$

where $\Gamma(\cdot)$ is the Gamma function. Note that real number $\theta \ge 1$ suffices to ensure $u' \in L^\infty(I; L^2(\Omega))$ and $\nabla u \in L^\infty(I; L^2(\Omega))$.

**Definition 3.7** The geometric partition $\mathcal{M}_{N_T, \sigma} = \{I_n\}_{n=1}^{N_T}$ of $I$ with grading factor $\sigma \in (0, 1)$ and $N_T$ time intervals $I_n = (t_{n-1}, t_n)$ is given by the nodes $t_0 = 0, t_n = \sigma^{N_T-n}, \forall n \in \mathbb{N}_+$.
The time step $\tau_1 = \sigma^{N_T-1}$ and for $2 \le n \le N_T$, $\tau_n = t_n - t_{n-1} = \lambda t_{n-1}$ with $\lambda = (1 - \sigma)/\sigma$. From the assumption (3.18) and above mesh setting, we have the following estimates for the solution $u$.

**Lemma 3.8** By assuming (3.18), let $M_{N_T, \sigma} = \{I_n\}_{n=1}^{N_T}$ be a geometric partition of I. Then we have

$$ \begin{aligned} |u|_{H^{s_1+1}(I_1;X)}^2 &\le C d^{2s_1} \Gamma(2s_1+1) \sigma^{(N_T-1)(2\theta-2s_1-1)}, \\ |u|_{H^{s_n+1}(I_n;X)}^2 &\le C d^{2s_n} \Gamma(2s_n+1) \sigma^{(N_T-n+1)(2\theta-2s_n-1)}, \quad s_n \ge 1, \forall n=2,3,\dots,N_T. \end{aligned} $$

The constant $C, d$ are independent of $N_T, n$ and $s_n$.

**Proof.** This is a consequence of assumption (3.18), geometric partition Definition 3.7 and properties of the $\Gamma$-function. $\square$

**Definition 3.9** An approximation degree vector $\mathbf{q} = \{q_n\}_{n=1}^{N_T}$ is called linear with slope $\mu > 0$ on the geometric partition $M_{N_T, \sigma}$, if $q_n = \lfloor \mu n \rfloor$ for $1 \le n \le N_T$.

**Theorem 3.10** Assume (3.18). Let $u_\tau \in V_\tau^q(I; H_0^1(\Omega))$ be the solution of (3.3) for a geometric partition $M_{N_T, \sigma}$ and a linearly increasing degree vector $\mathbf{q}$. Then there exists $\mu_0 > 0$, such that for any $\mathbf{q}$ with slope $\mu \ge \mu_0$, the following error estimates hold.

$$ |||e||| \lesssim \exp(-BN_{DoFs,t}^{\frac{1}{2}}), $$

with constants $B$ independent of temporal degrees of freedom $N_{DoFs,t} = \sum_{n=1}^{N_T}(q_n + 1)$.

**Proof.** From the Theorem 3.5, we need to bound

$$ A_{a,b} := \sum_{m=1}^{N_T} \left(\frac{\tau_m}{2}\right)^{2s_m-b} q_m^a \frac{(q_m - s_m)!}{(q_m + s_m)!} |u|_{H^{s_m+1}(I_m;X)}^2 $$

with $0 \le s_m \le \min\{s_{0,m}, q_m\}$, $a \in \mathbb{R}$, $b \le 2\theta - 1$. Lemma 3.8 gives that

$$ \begin{aligned} A_{a,b} &\le C \left(\frac{\lambda \sigma^{N_T-1}}{2}\right)^{2s_1-b} q_1^a \frac{(q_1 - s_1)!}{(q_1 + s_1)!} d^{2s_1} \Gamma(2s_1 + 1) \sigma^{(N_T-1)(2\theta-2s_1-1)} \\ &\quad + C \sum_{m=2}^{N_T} \left(\frac{\lambda \sigma^{N_T-m+1}}{2}\right)^{2s_m-b} q_m^a \frac{(q_m - s_m)!}{(q_m + s_m)!} d^{2s_m} \Gamma(2s_m + 1) \sigma^{(N_T-m+1)(2\theta-2s_m-1)} \\ &\quad = C \sum_{m=1}^{N_T} \sigma^{(N_T-m+1)(2\theta-1-b)} q_m^a \left(\frac{d\lambda}{2}\right)^{2s_m} \frac{(q_m - s_m)!}{(q_m + s_m)!} \Gamma(2s_m + 1), \end{aligned} $$

where $\lambda^{-b} \le 1$ is hidden in the constant $C$. Now setting $s_m = \alpha_m q_m$ with $\alpha_m \in (0, 1)$, we obtain with Stirling's formula

$$ A_{a,b} \le C \sum_{m=1}^{N_T} \sigma^{(N_T-m+1)(2\theta-1-b)} q_m^a \left[ \left( \frac{\alpha d \lambda}{2} \right)^{2s_m} \frac{(1-\alpha_m)^{1-\alpha_m}}{(1+\alpha_m)^{1+\alpha_m}} \right]^{q_m} $$

We denote the function $f_{\lambda,d}(\alpha) = (\alpha d\lambda/2)^{2\alpha} \frac{(1-\alpha)^{1-\alpha}}{(1+\alpha)^{1+\alpha}}$, which satisfies $0 < \inf_{0<\alpha<1} f_{\lambda,d}(\alpha) \le f_{\lambda,d}(\alpha_{\min}) < 1$ with $\alpha_{\min} = \frac{1}{\sqrt{1+\lambda^2 d^2/4}}$. For the completeness, the proof is postponed to the Appendix ??.

Next, by setting $f_{\min} := f_{\lambda,d}(\alpha_{\min})$ and select $\alpha_m = \alpha_{\min}$ for $1 \le m \le N_T$. Hence, for $q_m = \lfloor \mu m \rfloor$ we have

$$ \begin{aligned} A_{a,b} &\le C \sigma^{(N_T-1)(2\theta-1-b)} q_1^a f_{\min}^{q_1} + C \sum_{m=2}^{N_T} \sigma^{(N_T-m+1)(2\theta-1-b)} q_m^a f_{\min}^{q_m} \\ &\le C \sigma^{(N_T-1)(2\theta-1-b)} (\mu)^a f_{\min}^{\mu-1} + C \sum_{m=2}^{N_T} \sigma^{(N_T-m+1)(2\theta-1-b)} (\mu m)^a f_{\min}^{q_m} \end{aligned} $$
$$ \leq C \sigma^{(N_T-1)(2\theta-1-b)} (\mu N_T)^a \left( f_{\min}^{\mu-1} + \sum_{m=2}^{N_T} \sigma^{-m(2\theta-1-b)} f_{\min}^{\mu m} \right). $$

Let

$$ \mu > \frac{(2\theta - 1 - b) \ln(\sigma)}{\ln(f_{\min})}. \qquad (3.19) $$

Then it holds $l := \frac{f_{\min}^{\mu}}{\sigma^{2\theta-1-b}} < 1$ and, consequently,

$$ A_{a,b} \leq C \sigma^{(N_T-1)(2\theta-1-b)} (\mu N_T)^a \left( 1 + \sum_{m=2}^{N_T} l^m \right) \leq C \sigma^{(N_T-1)(2\theta-1-b)} (\mu N_T)^a \frac{l^2}{1-l}. \qquad (3.20) $$

Observing that $N_{DoFs,t} \approx \mathcal{O}(N_T^2)$, from Theorem 3.5, we bound errors by letting $a = 0, b = 0$ for energy norm error $\|e\|_{E}$ and $a = 1, b = 1$ for DG norm error $|||e|||$.

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

**Remark 3.11** For a piecewise analytical solution, we can show that both the energy norm error and the DG norm error decay exponentially under p-refinement in time for a fixed temporal mesh size $\tau_n := \tau$ and temporal polynomial degree $q_n := q$. Specifically, we obtain

$$ \|e\|_{E} \lesssim \exp(-\tilde{B}N_{DoFs,t}), $$

and

$$ |||e||| \lesssim \exp(-\tilde{B}N_{DoFs,t}), $$

where $\tilde{B}$ depends only on $\tau$, and $N_{DoFs,t} = N_T(q+1)$.

# 4 Fully-discrete scheme and its error estimates

In this section we analyze a space-time fully discretization that combine the DG time-stepping scheme in Section 3.1 and the conforming finite element method in space.

As the FEM space $V_h^p$ is a subspace of the $H_0^1(\Omega)$ space, the local fully discrete formulation on $Q_n := I_n \times \Omega$, for all $n = 1, \dots, N_T$, is defined as: Find $u_{\tau h} \in V_{\tau_n}^{q_n}(I; V_h^p)$ such that

$$ \int_{I_n} (u''_{\tau h}, v') + (\gamma u'_{\tau h}, v') + (\eta \nabla u'_{\tau h}, \nabla v') + (\zeta^2 \nabla u_{\tau h}, \nabla v') dt \\ + ([u'_{\tau h}]_{n-1}, v'(t_{n-1}^{+})) + (\zeta^2 \nabla [u_{\tau h}]_{n-1}, \nabla v(t_{n-1}^{+})) = \int_{I_n} (f, v') dt, \qquad (4.1) $$

for all $v \in V_{\tau_n}^{q_n}(I; V_h^p)$, with the initial condition $u'_{\tau h}(t_0^-) := P_h^p u_1$ and $u_{\tau h}(t_0^-) := \Pi_h^p u_0$, where $P_h^p$ is the L$^2$-projection operator defined in (2.11) and $\Pi_h^p$ is the elliptic projection operator defined in (2.12).

Similarly, following the semi-discrete scheme (3.3), the fully-discrete formulation for all the time steps: Find $u_{\tau h} \in V_{\tau}^q(I; V_h^p)$ such that

$$ B(u_{\tau h}, v) = \ell(v) \quad \forall v \in V_{\tau}^{q}(I; V_{h}^{p}). \qquad (4.2) $$

We mention that the stability estimates in the DG-norm (3.5) and energy estimates (3.8) still hold for the fully discrete scheme.

## 4.1 Abstract error bounds

We start with deriving the abstract error bound for in the DG-norm estimates. First, we split the error as follow:

$$ e := u - u_{\tau h} = (u - I_{\tau}^{q} P_{h}^{p} u) + (I_{\tau}^{q} P_{h}^{p} u - u_{\tau h}) = \rho + \theta, \qquad (4.3) $$

where the interpolation operator $I_{\tau}^q$ is defined in Definition 2.3.
By using the regularity result of (2.5), we have $u \in H^1(I; H_0^1(\Omega)) \hookrightarrow u \in C^0(I; H_0^1(\Omega))$. Therefore, using the properties (2.7), we have

$$ \nabla[\rho]_n = -\nabla[I_{\tau}^q P_h^p u]_n = \nabla(P_h^p u(t_n^+) - P_h^p u(t_n^-)) = 0, \quad (4.4) $$

for any $n = 1, \dots, N_T - 1$. In addition, using the properties (2.7)-(2.9) of the operator $I_{\tau}^q$ and the $L^2$-orthogonal property of the $P_h^p$ (2.11), the following relations hold for any $n = 1, \dots, N_T$, such that:

$$ (\rho'(t_n^-), \xi'(t_n^-)) = ((u - P_h^p u)'(t_n^-), \xi'(t_n^-)) = 0, \quad \text{and} \quad \int_{I_n} (\rho', \xi'') \, dt = 0 \quad \forall \xi \in V_{\tau_n}^{q_n}(I_n; V_h^p). \quad (4.5) $$

**Lemma 4.1** Let $u_{\tau h}$ and $u$ be the solution of (4.2) and (2.1), respectively. It holds that

$$ \|e\|_{E}^{2} \lesssim \int_{0}^{T} \|\gamma^{\frac{1}{2}} \rho'\|^{2} + \|\eta^{\frac{1}{2}} \nabla \rho'\|^{2} + \frac{\zeta^{4}}{\eta} \|\nabla \rho\|^{2} \, dt + \|\zeta \nabla \rho(t_{0}^{+})\|^{2} \quad (4.6) $$

and

$$ |||e|||^2 \lesssim \int_0^T \|\gamma^{\frac{1}{2}} \rho'\|^2 + \|\eta^{\frac{1}{2}} \nabla \rho'\|^2 + \frac{\zeta^4}{\eta} \|\nabla \rho\|^2 \, dt + \|\zeta \nabla \rho(t_0^+)\|^2 + \|\rho'(t_0^+)\|^2 + \|\rho'(t_{N_T}^-)\|^2 + \sum_{n=2}^{N_T} \|[\rho']_{n-1}\|^2. \quad (4.7) $$

**Proof.** Using the consistency of the DG in time methods, we have $B(e, v_{\tau h}) = 0$, for any $v_{\tau h} \in V_{\tau}^{q}(I; V_h^p)$. Therefore, choosing $v_{\tau h} = \theta$, we obtain the following error equation

$$ |||\theta|||^2 = B(\theta, \theta) = -B(\rho, \theta). \quad (4.8) $$

Next, using integration by parts for the first term in $I_n$, we obtain

$$ \begin{aligned} B(\rho, \theta) &= \sum_{n=1}^{N_T} \int_{I_n} (\rho'', \theta') + (\gamma\rho', \theta') + (\eta\nabla\rho', \nabla\theta') + (\zeta^2\nabla\rho, \nabla\theta') \, dt \\ &\quad + \sum_{n=1}^{N_T-1} ([\rho']_n, v'(t_n^+)) + (\zeta^2\nabla[\rho]_n, \nabla\theta(t_n^+)) + (\rho'(t_0^+), \theta'(t_0^+)) + (\zeta^2\nabla\rho(t_0^+), \nabla\theta(t_0^+)) \\ &= \sum_{n=1}^{N_T} \int_{I_n} -(\rho', \theta'') + (\gamma\rho', \theta') + (\eta\nabla\rho', \nabla\theta') + (\zeta^2\nabla\rho, \nabla\theta') \, dt \\ &\quad - \sum_{n=1}^{N_T-1} (\rho'(t_n^-), [v']_n) + (\rho'(t_{N_T}^-), \theta'(t_{N_T}^-)) + \sum_{n=1}^{N_T-1} (\zeta^2\nabla[\rho]_n, \nabla\theta(t_n^-)) + (\zeta^2\nabla\rho(t_0^+), \nabla\theta(t_0^+)). \end{aligned} \quad (4.9) $$

Using the interpolation properties (4.4) and (4.5), it gives the local estimate

$$ B(\rho, \theta) = \sum_{n=1}^{N_T} \int_{I_n} (\gamma\rho', \theta') + (\eta\nabla\rho', \nabla\theta') + (\zeta^2\nabla\rho, \nabla\theta') \, dt + (\zeta^2\nabla\rho(t_0^+), \nabla\theta(t_0^+)). \quad (4.10) $$

Using the Cauchy-Schwartz inequality, we have

$$ \begin{aligned} |B(\rho, \theta)| &\le \sum_{n=1}^{N_T} \int_{I_n} \|\gamma^{\frac{1}{2}} \rho'\| \|\gamma^{\frac{1}{2}} \theta'\| + \|\eta^{\frac{1}{2}} \nabla \rho'\| \|\eta^{\frac{1}{2}} \nabla \theta'\| + \zeta^2 \|\nabla \rho\| \|\nabla \theta'\| \, dt + \|\zeta \nabla \rho(t_0^+)\| \|\zeta \nabla \theta(t_0^+)\|, \\ &\le \sum_{n=1}^{N_T} \int_{I_n} \|\gamma^{\frac{1}{2}} \rho'\| \|\gamma^{\frac{1}{2}} \theta'\| + \|\eta^{\frac{1}{2}} \nabla \rho'\| \|\eta^{\frac{1}{2}} \nabla \theta'\| + \frac{\zeta^2}{\eta^{\frac{1}{2}}} \|\nabla \rho\| \|\eta^{\frac{1}{2}} \nabla \theta'\| \, dt + \|\zeta \nabla \rho(t_0^+)\| \|\zeta \nabla \theta(t_0^+)\| \end{aligned} $$

and

$$ |||\theta|||^2 \le |B(\rho, \theta)| \lesssim \left( \sum_{n=1}^{N_T} \int_{I_n} \|\gamma^{\frac{1}{2}} \rho'\|^2 + \|\eta^{\frac{1}{2}} \nabla \rho'\|^2 + \frac{\zeta^4}{\eta} \|\nabla \rho\|^2 \, dt + \|\zeta \nabla \rho(t_0^+)\|^2 \right)^{\frac{1}{2}} |||\theta|||. \quad (4.11) $$

By using triangle inequality, $\|e\|_E \le \|\theta\|_E + \|\rho\|_E \le |||\theta||| + \|\rho\|_E$ together with the above relation (4.11), we derive the error bound (4.6). Similarly, we derive the error bound (4.7) by using the triangle inequality $|||e||| \le |||\theta||| + |||\rho|||$. $\square$
## 4.2 Convergence rate for the fully discrete scheme

We derive the convergence estimates for the fully-discretization corresponding to (4.1), which are explicit in the spatial mesh size, time step size, and the polynomial degrees in space and time.

**Theorem 4.2** Let $u_{h\tau}$ and $u$ be the solution of (4.1) and (2.5), respectively. Assume that $u|_{I_n} \in H^{k_n+1}(I_n; H_0^1(\Omega)) \cap H^2(I_n; H^{r_0+1}(\Omega))$, $k_n, r_0 \ge 1$. There it holds that

$$ \|e\|_{E}^{2} \lesssim \sum_{n=1}^{N_T} \left( \gamma \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2 + p^2 \left( \eta + \frac{\tau_n^2 \zeta^4}{q_n^2 \eta} \right) \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; H_0^1(\Omega))}^2 \right. \\ \left. + \left( \gamma \frac{h^2}{p^4} + \eta \right) \frac{h^{2r}}{p^{2r_0-2}} \|u'\|_{L^2(I_n; H^{r_0+1}(\Omega))}^2 + \frac{\zeta^4}{\eta} \frac{h^{2r}}{p^{2r_0-2}} \|u\|_{L^2(I_n; H^{r_0+1}(\Omega))}^2 \right) \\ + \zeta^2 \frac{h^{2r}}{p^{2r_0-2}} \|u(t_0^+)\|_{H^{r_0+1}(\Omega)}^2 \tag{4.12} $$

and

$$ \|e\|^2 \lesssim \sum_{n=1}^{N_T} \left( (1 + \gamma \tau_n) \tau_n^{2s_n-1} q_n \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2 \right. \\ + p^2 \left( \eta + \frac{\tau_n^2 \zeta^4}{q_n^2 \eta} \right) \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; H_0^1(\Omega))}^2 + \left( \gamma \frac{h^2}{p^4} + \eta \right) \frac{h^{2r}}{p^{2r_0-2}} \|u'\|_{L^2(I_n; H^{r_0+1}(\Omega))}^2 \\ + \frac{\zeta^4}{\eta} \frac{h^{2r}}{p^{2r_0-2}} \|u\|_{L^2(I_n; H^{r_0+1}(\Omega))}^2 \right) + \zeta^2 \frac{h^{2r}}{p^{2r_0-2}} \|u(t_0^+)\|_{H^{r_0+1}(\Omega)}^2 \\ + \frac{h^{2r+2}}{p^{2r_0+2}} \|u'(t_0^+)\|_{H^{r_0+1}(\Omega)}^2 + \frac{h^{2r+2}}{p^{2r_0+2}} \|u'(t_{N_T}^-)\|_{H^{r_0+1}(\Omega)}^2 + \sum_{n=2}^{N_T} \left( \frac{h^{2r+2}}{p^{2r_0+2}} \|u'(t_{n-1})\|_{H^{r_0+1}(\Omega)}^2 \right) \tag{4.13} $$

for any integer $1 \le s_n \le \min\{q_n, k_n\}$, $1 \le r \le \min\{p, r_0\}$.

**Proof.** Recalling the Lemma 4.1, the DG-norm error (4.7) is bounded by seven terms $T_i$, $i = 1, \dots, 7$. We will bound them by using the approximation property (2.10), (2.15) and (2.14) separately.

For the $T_1$ term, using the triangle inequality, orthogonality and stability of spatial $L^2$-orthogonal projection $P_h^p$ in $L^2$-norm, second bound in (2.10) and (2.14), we have

$$ \begin{aligned} T_1 &\le \sum_{n=1}^{N_T} \gamma \left( \|u' - P_h^p(u')\|_{L^2(I_n; L^2(\Omega))}^2 + \|P_h^p(u - I^{q_n}(u))'\|_{L^2(I_n; L^2(\Omega))}^2 \right) \\ &\le \sum_{n=1}^{N_T} \gamma \left( \|u' - P_h^p(u')\|_{L^2(I_n; L^2(\Omega))}^2 + \|(u - I^{q_n}(u))'\|_{L^2(I_n; L^2(\Omega))}^2 \right) \tag{4.14} \\ &\lesssim \sum_{n=1}^{N_T} \gamma \left( \frac{h^{2r+2}}{p^{2r_0+2}} \|u'\|_{L^2(I_n; H^{r_0+1}(\Omega))}^2 + \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2 \right). \end{aligned} $$

Next, using the triangle inequality, stability of spatial $L^2$-orthogonal projection $P_h^p$ in $H^1$-norm (2.13), second bound in (2.10) and (2.15), we have

$$ \begin{aligned} T_2 &\le \sum_{n=1}^{N_T} 2\eta \left( \|\nabla(u' - P_h^p(u'))\|_{L^2(I_n; L^2(\Omega))}^2 + \|\nabla(P_h^p(u - I^{q_n}(u))')\|_{L^2(I_n; L^2(\Omega))}^2 \right) \\ &\lesssim \sum_{n=1}^{N_T} \eta \left( \|\nabla(u' - P_h^p(u'))\|_{L^2(I_n; L^2(\Omega))}^2 + p^2 \|\nabla(u - I^{q_n}(u))'\|_{L^2(I_n; L^2(\Omega))}^2 \right) \tag{4.15} \\ &\lesssim \sum_{n=1}^{N_T} \eta \left( \frac{h^{2r}}{p^{2r_0-2}} \|u'\|_{L^2(I_n; H^{r_0+1}(\Omega))}^2 + p^2 \tau_n^{2s_n} \frac{(q_n - s_n)!}{(q_n + s_n)!} |\nabla u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2 \right). \end{aligned} $$
Similarly, using the triangle inequality, stability of spatial L²-orthogonal projection $P_h^p$ in H¹-norm (2.13), first bound in (2.10) and (2.15), we have

$$
\begin{align}
T_3 &\le \sum_{n=1}^{N_T} 2 \frac{\zeta^4}{\eta} \left( \|\nabla(u - P_h^p(u))\|_{L^2(I_n; L^2(\Omega))}^2 + \|\nabla(P_h^p(u - I^{q_n}(u)))\|_{L^2(I_n; L^2(\Omega))}^2 \right) \nonumber \\
&\lesssim \sum_{n=1}^{N_T} \frac{\zeta^4}{\eta} \left( \|\nabla(u - P_h^p(u))\|_{L^2(I_n; L^2(\Omega))}^2 + p^2 \|\nabla(u - I^{q_n}(u))\|_{L^2(I_n; L^2(\Omega))}^2 \right) \tag{4.16} \\
&\lesssim \sum_{n=1}^{N_T} \frac{\zeta^4}{\eta} \left( \frac{h^{2r}}{p^{2r_0-2}} \|u\|_{L^2(I_n; H^{r_0+1}(\Omega))}^2 + p^2 \frac{\tau_n^{2s_n+2}}{q_n^2} \frac{(q_n - s_n)!}{(q_n + s_n)!} |\nabla u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2 \right). \nonumber
\end{align}
$$

For the $T_4$, using the interpolation property (2.7) and approximation result (2.15), we have

$$ T_4 = \zeta^2 \left( \|\nabla(u - P_h^p(u))(t_0^+)\|_{L^2(\Omega)}^2 \right) \lesssim \zeta^2 \frac{h^{2r}}{p^{2r_0-2}} \|u(t_0^+)\|_{H^{r_0+1}(\Omega)}^2. \quad (4.17) $$

For the $T_5$ term, using the triangle inequality, orthogonality and stability of spatial L²-orthogonal projection $P_h^p$ in L²-norm, third bound in (2.10) and (2.14), we have

$$
\begin{align}
T_5 &\le \left( \| (u' - P_h^p u')(t_0^+) \|_{L^2(\Omega)}^2 + \| P_h^p (u - I_{\tau_n}^{q_n} u)'(t_0^+) \|_{L^2(\Omega)}^2 \right) \nonumber \\
&\le \left( \| (u' - P_h^p u')(t_0^+) \|_{L^2(\Omega)}^2 + \| (u - I_{\tau_n}^{q_n} u)'(t_0^+) \|_{L^2(\Omega)}^2 \right) \tag{4.18} \\
&\lesssim \left( \frac{h^{2r+2}}{p^{2r_0+2}} \| u'(t_0^+) \|_{L^2(I_1; H^{r_0+1}(\Omega))}^2 + \tau_1^{2s_1-1} q_1 \frac{(q_1-s_1)!}{(q_1+s_1)!} |u|_{H^{s_1+1}(I_1; L^2(\Omega))}^2 \right). \nonumber
\end{align}
$$

For the $T_6$ term, using property (2.8) for temporal projection $I_{\tau_{NT}}^{q_{NT}}$ and spatial L²-orthogonal projection $P_h^p$ in L²-norm, and (2.14), we have

$$ T_6 \le \|\rho'(t_N^-)\|^2 = \|(u - P_h^p u)'(t_N^-)\|^2 \lesssim \left( \frac{h^{2r+2}}{p^{2r_0+2}} \|u'(t_N^-)\|_{L^2(I_{NT}; H^{r_0+1}(\Omega))}^2 \right). \quad (4.19) $$

For the $T_7$ term, using the triangle inequality, the property (2.8) for temporal projection $I_{N_T}^{q_{N_T}}$, orthogonality and stability of spatial L²-orthogonal projection $P_h^p$ in L²-norm, third bound in (2.10), (2.14), and the regularity $u'(t_{n-1}^+) = u'(t_{n-1}^-)$, we have

$$
\begin{align}
T_7 &\lesssim \sum_{n=2}^{N_T} (\|\rho'(t_{n-1}^+)\|^2 + \|\rho'(t_{n-1}^-)\|^2) \nonumber \\
&\lesssim \sum_{n=2}^{N_T} \left( \| (u' - P_h^p u')(t_{n-1}^+) \|_{L^2(\Omega)}^2 + \| P_h^p (u - I_{\tau_n}^{q_n} u)'(t_{n-1}^+) \|_{L^2(\Omega)}^2 + \| (u' - P_h^p u')(t_{n-1}^-) \|_{L^2(\Omega)}^2 \right) \tag{4.20} \\
&\le \sum_{n=2}^{N_T} \left( \| (u' - P_h^p u')(t_{n-1}^+) \|_{L^2(\Omega)}^2 + \| (u - I_{\tau_n}^{q_n} u)'(t_{n-1}^+) \|_{L^2(\Omega)}^2 + \| (u' - P_h^p u')(t_{n-1}^-) \|_{L^2(\Omega)}^2 \right) \nonumber \\
&\lesssim \sum_{n=2}^{N_T} \left( \frac{h^{2r+2}}{p^{2r_0+2}} \| u'(t_{n-1}) \|_{L^2(I_n; H^{r_0+1}(\Omega))}^2 + \tau_n^{2s_n-1} q_n \frac{(q_n-s_n)!}{(q_n+s_n)!} |u|_{H^{s_n+1}(I_n; L^2(\Omega))}^2 \right). \nonumber
\end{align}
$$

Finally, combining about above bounds (4.14)–(4.20), we arrive the error estimates (4.12) and (4.13). $\square$

**Corollary 4.3** Let $u_{h\tau}$ and $u$ be the solution of (4.2) and (2.1) respectively. If we assume time steps $\tau_n = \tau$, polynomial degrees $q_n = q$, regularity indices $k_n = k$, $s_n = s$, for $n = 1, \dots, N_T$ and quasi-uniform mesh assumption, the above error bounds read

$$ \left( \int_0^T \| \gamma^{\frac{1}{2}} e' \|^2 + \| \eta^{\frac{1}{2}} \nabla e' \|^2 dt \right)^{\frac{1}{2}} \approx \mathcal{O} \left( (1+p) \frac{\tau^s}{q^s} \right) + \mathcal{O} \left( \frac{h^r}{p^{r_0-1}} \right), $$
and

$$ |||e||| \approx \mathcal{O} \left( \frac{\tau^{s-\frac{1}{2}}}{q^{s-\frac{1}{2}}} + (1+p) \frac{\tau^s}{q^s} \right) + \mathcal{O} \left( \frac{h^r}{p^{r_0-1}} \right), $$

For a fixed $p$, both of the error bounds mentioned above are optimal with respect to the time step size $\tau$, mesh size $h$, and polynomial degree $q$. However, if we assume $q = p$, then while the error bounds in the energy-norm and DG-norm remain optimal in $h$ and $\tau$, they become suboptimal by one order in $p$ and by one-half order in $q$, respectively.

**Remark 4.4 (Elliptic projection)** The above error bound can be derived by using the elliptic projection rather than the $L^2$-orthogonal projection. Under the assumption that the computation domain is convex, it is possible to derive the DG-norm error bound (4.13) with an hp-optimal bound.

# 5 Numerical results

In this section, we numerically investigate the a priori error estimates proven in Theorem 4.2 and Corollary 4.3. In addition, we also numerically investigate the norms which are not covered by the theory in this work.

Let the spatial domain $\Omega = (0, 1)^2$ be partitioned into uniform tensor-product meshes and the time interval be $I = (0, 1]$. The choice of time steps will be different in subsection 5.1 and 5.2 with the parameters $\gamma = \eta = \zeta = 1$. In subsection 5.3, 5.4 and 5.5, the parameters $\gamma, \eta, \zeta$ will be different.

## 5.1 Semi-discretization in time

In this subsection, we investigate the semi-discretization of the DG time-stepping scheme by solving a manufactured solution that is a quadratic polynomial function in the spatial variables. In the practical computation, we will use a tensor-product $Q_2$ polynomial basis with 4 spatial elements.

**Problem (I).** We choose the analytical solution $u$ to be a smooth function defined as follows:

$$ u(x, y, t) = \sin(2\pi t)x(1-x)y(1-y). $$

First, we explore the performance of the DG time-stepping scheme with a fixed polynomial order $q \in \{1, 2, 3, 4\}$ for all equidistant time steps $N_T \in \{2, 4, 8, 16, 32, 64, 128\}$. The total number of time DoFs is $N_{DoFs} = N_T(q+1)$ with $N_T = 1/\tau$.

In Figure 1, we present the convergence plots for the error $e = u - u_\tau$ in four norms, namely $|||e|||$, $||e||_E$, $||e||_{W_\infty^1(I;L^2(\Omega))}$, and $||\nabla e||_{L^\infty(I;L^2(\Omega))}$. We observe that the convergence rates in the DG-norm $|||e|||$ and the energy-norm $||e||_E$ are $O(\tau^{q-\frac{1}{2}})$ and $O(\tau^q)$, respectively. These rates are optimal and consistent with Corollary 4.3. Moreover, the convergence rate in $||e||_{W_\infty^1(I;L^2(\Omega))}$ is optimal at $O(\tau^q)$ for $q \in \{1, 2, 3, 4\}$. In addition, the rate in $||\nabla e||_{L^\infty(I;L^2(\Omega))}$ is $O(\tau^{q+1})$ for $q \in \{2, 3, 4\}$, which is also optimal. For $q=1$, the convergence rate is $O(\tau)$, which is suboptimal by one order in $\tau$.

Next, we will explore the p-version of the DG time-stepping scheme. So we fix the equidistant partition in time with $N_T = 5, 10$ steps for all uniform time polynomial orders $q = 1, 2, 3, 4, 5, 6, 7$. In Figure 2, we observe exponential convergence in terms of the number of temporal DoFs in four norms, respectively, as mentioned in Remark 3.11.

**Problem (II).** Next, we choose the singular analytical solution $u$ defined as follow:

$$ u(x, y, t) = t^\theta x(1-x)y(1-y), $$

with parameter $\theta > 1$. For any $\theta > 0$, the solution $u \in H^{\theta+0.5-\epsilon}(I; V)$ for any arbitrarily small $\epsilon > 0$.

In this example, we choose $\theta = 3.51$ in Problem (II). In Figure 3, the numerical results under uniform $\tau$-refinement with time steps $N_T \in \{2, 4, 8, 16, 32, 64, 128, 256, 512, 1024\}$ show that the convergence rates are $O(\tau^{\min\{q,3.01-\epsilon\}-0.5})$ and $O(\tau^{\min\{q,3.01-\epsilon\}})$ for the DG-norm $|||e|||$ and the energy-norm $||e||_E$, respectively. These numerical results match the theory in Corollary 4.3, and the convergence rates are optimal. Moreover, we observe that the convergence rate in $||e||_{W_\infty^1(I;L^2(\Omega))}$ is optimal, with $O(\tau^{\min\{q,2.51\}})$ for $q \in \{1, 2, 3, 4\}$. In addition, the convergence rate in $||\nabla e||_{L^\infty(I;L^2(\Omega))}$ is optimal, at $O(\tau^{\min\{q+1,3.51\}})$ for $q \in \{2, 3, 4\}$. For $q=1$, the convergence rate is $O(\tau)$, which is suboptimal by one order in $\tau$.
<figure><img src="image_3.png" /></figure>

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

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

<figure><img src="image_6.png" /><figcaption>Fig. 1. Problem (I). Uniform $\tau$-refinement for the DG time-stepping scheme for temporal polynomial order $q = 1, 2, 3, 4$.</figcaption></figure>

Next, we consider the *p*-version of the DG time-stepping scheme. Due to the limited regularity in time, we observe algebraic convergence rates with respect to the time polynomial order $q$ in Figure 4. The convergence rates in the DG-norm $|||e|||$ and the energy-norm $||e||_E$ are about $O(q^{6.5})$ and $O(q^{6.5})$, respectively. These rates are almost twice as large as the theoretical expectation. This doubling order convergence phenomenon is well documented in *p*-version FEM for solutions with singularities located at the mesh endpoints; see [29, 35]. Moreover, similar doubling order convergence rates are also observed in $||e||_{W^1(\Omega)}$ and $||\nabla e||_{L^\infty(\Omega)}$, with rates of $O(q^{6.0})$ and $O(q^{4.5})$, respectively.

Finally, we explore the *hp*-version DG time-stepping scheme for the aforementioned singular solution. In accordance with Theorem 3.10, we employ geometrically refined time intervals $M_{N_T, \sigma} = \{I_n\}_{n=1}^{N_T}$ toward the singularity, using a grading factor $\sigma$, combined with linearly increasing approximation orders $q_n = \lceil \mu n \rceil$. In Figure 5, for a fixed mesh grading factor $\sigma = 0.25$, we observe exponential convergence in terms of DoFs $\frac{1}{2}$ across four norms for different values of $\mu \in \{1, 1.5, 1.75, 2, 2.25, 2.5\}$. The numerical results indicate that the exponent of the exponential convergence depends on $\mu$, and the optimal value appears to be $\mu = 2$.

In addition, we explore the exponential convergence for a fixed value $\mu = 2$, with a varying value of $\sigma \in \{0.15, 0.17, 0.2, 0.25, 0.3, 0.35, 0.5\}$. In Figure 6, we observe exponential convergence in terms of DoFs $\frac{1}{2}$ in four norms, for various values of $\sigma$. The optimal mesh grading factor appears to be $\sigma = 0.25$. It is worth mentioning that *hp*-refinement with $\sigma = 0.5$, which commonly generated by the adaptive mesh refinement algorithm, has a much slower convergence rate for all norms compare to the other choice of $\sigma$.
<figure><img src="image_7.png" /><figcaption>Fig. 2. Problem (I). Uniform q-refinement for the DG time-stepping scheme for each equidistant time partition $N_T = 5, 10$.</figcaption></figure>

## 5.2 Space-time uniform mesh refinement

In this subsection, we consider the fully discrete error bound for the DG time-stepping scheme and FEM in space. We choose the manufactured solution as follows:

**Problem (III).** The analytic solution

$$u(x, y, t) = \sin(\pi x) \sin(\pi y) \sin(\pi t).$$

Following the theory in Corollary 4.3, we set $q_n = q = p$ and $\tau_n = \tau = h$, with $h \in \{\frac{1}{2}, \frac{1}{4}, \frac{1}{8}, \frac{1}{16}, \frac{1}{32}, \frac{1}{64}\}$ and $p \in \{1, 2, 3, 4\}$. In Figure 7, we observe that the optimal rates for DG-norm $|||e|||$ and the energy-norm $||e||_E$ are $O(h^{p-\frac{1}{2}})$ and $O(h^p)$, respectively.

## 5.3 Spatial convergence rate with different value of $\eta$

In this subsection, we investigate the dependency for the spatial convergence rate on the parameter $\eta$. As in practical relevant models, $\eta$ is always chosen to be much smaller than $\gamma$ and $\zeta$.

**Problem (IV).** We here consider the exact solution

$$u(x, y, t) = \sin(2\pi x) \sin(2\pi y) t^2$$

and set $\gamma = \zeta = 1$ and variable $\eta = 10^{-i}(i = 0, 3, 4, 5, 6, \infty)$. In order to control the discretization error in time, we use 4 time steps with $q = 2$. For spatial discretization, we use a sequences square meshes consist of $\{4, 16, 64, 256, 1024, 4096\}$ elements with polynomial order $p \in \{2, 3\}$ to test the spatial convergence rates. In Fig. 8, we observe the energy-norm $||e||_E$ convergence rate is $O(h^{p+1})$ for $\eta = 0$ which is optimal (recall that $||e||_E$ only contains the $L^2$-norm in space.) For other choice of $\eta > 0$, the
<figure><img src="image_8.png" /><figcaption>Fig. 3. Problem (II). Uniform $\tau$-refinement for the DG time-stepping scheme for temporal polynomial order $q = 1, 2, 3, 4$.</figcaption></figure>

convergence rate changes from $O(h^{p+1})$ to $O(h^p)$ under the mesh refinement. Numerical experiments show that this change of convergence rate happens when the mesh size $h^2 \approx c\eta$, which is the optimal rate in the energy-norm $\|e\|_{E}$.

## 5.4 Energy conservation with $\gamma = \eta = 0$

In this subsection, we will explore the energy conservation for the case with $\gamma = \eta = 0$, $f = 0$ and $\zeta = 1$. By setting the initial conditions are the same as problem (I) and the final time $t_{N_T} = 100$. We plot the space energy of approximated solution (3.3) as a function of time $t_n$. In Figure 9, we observe that proposed the DG time-stepping scheme and FEM in space is dissipative. Moreover, this dissipation can be reduced by increasing the number of time steps for a fixed time polynomial order $q \in \{2, 3\}$. Additionally, the dissipation for the higher-order scheme is less compared to that of the lower-order scheme in time. This numerical observation confirms the theoretical result reported in [2].

## 5.5 Energy decay for the diffusive wave with $\zeta = 0$

In this subsection, we explore the energy decay for the case with $\zeta = 0$, $f = 0$ and variable $\gamma$ and $\eta$. By setting the initial condition $u_1 = 2\pi x(1-x)y(1-y)$, $N_T = 50$ with $p = 2$ and 4 spatial elements in the current problem, we plot the value of $||u'(t_n^-)||$ in Figure 10. Recalling the result in Lemma 3.2, the estimate $\ln(||u'(t_n^-)||) \le \ln(|v_0|) - \frac{\gamma+\frac{\eta}{C_p}}{4(q+1)^2}t_n$ implies that $||u'(t_n^-)||$ decays exponentially with respect to $t_n$, and the exponent depends linearly on $(\gamma + \frac{\eta}{C_p})$.

We test two cases for the exponential decay of the $||u'(t_n^-)||$. The first case is to set $N_T = 50$, $q \in \{2, 4, 6\}$, $\gamma = 0$, and a variable value of $\eta \in \{0, 0.01, 0.1, 0.5, 1\}$. We observe that the exponent of
<figure><img src="image_9.png" /></figure>

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

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

<figure><img src="image_12.png" /><figcaption>Fig. 4. Problem (II). Uniform q-refinement for the DG time-stepping scheme.</figcaption></figure>

the exponential decay is greater for larger $\eta$. In the second case, we set $N_T = 50$, $q \in \{2, 4, 6\}$, $\eta = 0$, and a variable value of $\gamma \in \{0, 0.01, 0.1, 0.5, 1\}$. We observe that the exponent of the exponential decay is greater for larger $\gamma$.

## Acknowledgements

This work was partially supported by National Key R & D Program of the Ministry of Science and Technology of China (Grant NOs. 2020YFA0713400, 2020YFA0713401), National Natural Science Foundation of China (Grant NO. 12471388).

## References

[1] P. F. ANTONIETTI, M. BOTTI, AND I. MAZZIERI, *Discontinuous Galerkin discretization of coupled poroelasticity-elasticity problems*, IMA Journal of Numerical Analysis, (2024), p. drae093.

[2] P. F. ANTONIETTI, I. MAZZIERI, N. DAL SANTO, AND A. QUARTERONI, *A high-order discontinuous Galerkin approximation to ordinary differential equations with applications to elastodynamics*, IMA J. Numer. Anal., 38 (2017), pp. 1709–1734.

[3] P. F. ANTONIETTI, I. MAZZIERI, AND F. MIGLIORINI, *A space-time discontinuous Galerkin method for the elastic wave equation*, J. Comput. Phys., 419 (2020), p. 109685.

[4] I. BABUŠKA AND M. SURI, *The h-p version of the finite element method with quasi-uniform meshes*, RAIRO Modél. Math. Anal. Numér., 21 (1987), pp. 199–238.
<figure><img src="image_13.png" /><figcaption>Fig. 5. Problem (II). $\tau q$-refinement for the DG time-stepping scheme with $\sigma = 0.25$ and variable $\mu$.</figcaption></figure>

[5] J. M. CARCIONE, G. C. HERMAN, AND A. P. E. TEN KROODE, *Seismic modeling*, Geophysics, 67 (2002), pp. 1304–1325.

[6] P. CASTILLO, B. COCKBURN, D. SCHÖTZAU, AND C. SCHWAB, *Optimal a priori error estimates for the hp-version of the local discontinuous Galerkin method for convection-diffusion problems*, Math. Comp., 71 (2002), pp. 455–478.

[7] K. CHRYSAFINOS AND L. S. HOU, *Error estimates for semidiscrete finite element approximations of linear and semilinear parabolic equations under minimal regularity assumptions*, SIAM J. Numer. Anal., 40 (2002), pp. 282–306.

[8] L. DIENING, J. STORN, AND T. TSCHERPEL, *On the Sobolev and $L^p$-stability of the $L^2$-projection*, SIAM J. Numer. Anal., 59 (2021), pp. 2571–2607.

[9] Z. DONG AND A. ERN, *hp-error analysis of mixed-order hybrid high-order methods for elliptic problems on simplicial meshes*. working paper or preprint, Oct. 2024.

[10] Z. DONG, E. H. GEORGOULIS, AND P. J. HERBERT, *A hypocoercivity-exploiting stabilised finite element method for Kolmogorov equation*, Accepted in SIAM J. Numer. Anal., (2025).

[11] Z. DONG, L. MASCOTTO, AND Z. WANG, *A priori and a posteriori error estimates of a DG-CG method for the wave equation in second order formulation*, 2024.

[12] A. ERN AND J.-L. GUERMOND, *Finite Elements III: first-order and time-dependent PDEs*, Springer Cham, 2021.
<figure><img src="image_14.png" /></figure>

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

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

<figure><img src="image_17.png" /><figcaption>Fig. 6. Problem (II). Under $\tau q$-refinement for the DG time-stepping scheme with $\mu = 2$ and variable $\sigma$.</figcaption></figure>

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

<figure><img src="image_19.png" /><figcaption>Fig. 7. Problem (III). Space-time uniform mesh refinement with $p = q$ and $h = \tau$.</figcaption></figure>

[13] C. GABRIELE, J. G. MARTIN, AND M. ILARIO, *Discontinuous Galerkin time integration for second-order differential problems: formulations, analysis, and analogies*, 2025.

[14] W. HAN, J. GAO, Y. ZHANG, AND W. XU, *Well-posedness of the diffusive-viscous wave equation arising in geophysics*, J. Math. Anal. Appl., 486 (2020), p. 123914.

[15] W. HAN, C. SONG, F. WANG, AND J. GAO, *Numerical analysis of the diffusive-viscous wave equation*, Comput. Math. Appl., 102 (2021), pp. 54–64.
<figure><img src="image_20.png" /><figcaption>Fig. 8. Problem (IV). The convergence rates' varying with $\eta$ for the Cartesian mesh sequences 4, 16, 64, 256, 1024, 4096 elements and $p \in \{2, 3\}$ in energy-norm.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 9. Computed energy $\|u'_{\tau}\|^2 + \|\zeta \nabla u_{\tau}\|^2$ for $q = 2, 3$.</figcaption></figure>

[16] Z. HE, X. XIONG, AND L. BIAN, *Numerical simulation of seismic low-frequency shadows and its application*, Appl. Geophys., 5 (2008), pp. 301–306.

[17] T. J. HUGHES AND G. M. HULBERT, *Space-time finite element methods for elastodynamics: Formulations and error estimates*, Comput. Methods Appl. Mech. Engrg., 66 (1988), pp. 339–363.

[18] G. M. HULBERT AND T. J. HUGHES, *Space-time finite element methods for second-order hyperbolic equations*, Comput. Methods Appl. Mech. Engrg., 84 (1990), pp. 327–348.

[19] P. JAMET, *Galerkin-type approximations which are discontinuous in time for parabolic equations in a variable domain*, SIAM J. Numer. Anal., 15 (1978), pp. 912–928.

[20] C. JOHNSON, *Discontinuous Galerkin finite element methods for second order hyperbolic problems*, Comput. Methods Appl. Mech. Engrg., 107 (1993), pp. 117–129.

[21] M. KARKULIK AND J. M. MELENK, *Local high-order regularization and applications to hp-methods*, Comput. Math. Appl., 70 (2015), pp. 1606–1639.

[22] V. A. KORNEEV, G. M. GOLOSHUBIN, T. M. DALEY, AND D. B. SILIN, *Seismic low-frequency effects in monitoring fluid-saturated reservoirs*, Geophysics, 69 (2004), pp. 522–532.

[23] D. LING, C.-W. SHU, AND W. YAN, *Local discontinuous Galerkin methods for diffusive-viscous wave equations*, J. Comp. Appl. Math., 419 (2023), p. 114690.
<figure><img src="image_22.png" /></figure>

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

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

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

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

<figure><img src="image_27.png" /><figcaption>Fig. 10. The exponential decay rates of ||$u'(t_n^-)$|| for varying parameters $\gamma$ and $\eta$ with $q=2, 4, 6$.</figcaption></figure>

[24] C. MAKRIDAKIS AND R. H. NOCHETTO, *A posteriori error analysis for higher order dissipative methods for evolution problems*, Num. Math., 104 (2006), pp. 489–514.

[25] V. MENSAH, A. HIDALGO, AND R. M. FERRO, *Numerical modelling of the propagation of diffusive-viscous waves in a fluid-saturated reservoir using finite volume method*, Geophys. J. Int., 218 (2019), pp. 33–44.

[26] K. MUSTAPHA, H. BRUNNER, H. MUSTAPHA, AND D. SCHÖTZAU, *An hp-version discontinuous Galerkin method for integro-differential equations of parabolic type*, SIAM J. Numer. Anal., 49 (2011), pp. 1369–1396.
[27] D. SCHÖTZAU AND C. SCHWAB, *An hp a priori error analysis of the DG time-stepping method for initial value problems*, Calcolo, 37 (2000), pp. 207–232.

[28] ——, *Time discretization of parabolic problems by the hp-version of the discontinuous Galerkin finite element method*, SIAM J. Numer. Anal., 38 (2000), pp. 837–875.

[29] C. SCHWAB, *p-and hp-finite element methods. Theory and applications in solid and fluid mechanics*, The Clarendon Press, Oxford University Press, 1998.

[30] V. THOMÉE, *Galerkin finite element methods for parabolic problems. Second edition*, Springer-Verlag, Berlin, 2006.

[31] N. J. WALKINGTON, *Combined DG-CG time stepping for wave equations*, SIAM J. Numer. Anal., 52 (2014), pp. 1398–1417.

[32] W. WANG, W. YAN, AND D. YANG, *A cell-centered finite volume scheme for the diffusive-viscous wave equation on general polygonal meshes*, Appl. Math. Lett., 133 (2022), p. 108274.

[33] T. WARBURTON AND J. S. HESTHAVEN, *On the constants in hp-finite element trace inverse inequalities*, Comput. Methods Appl. Mech. Engrg., 192 (2003), pp. 2765–2773.

[34] L. WEI, YICHEN; YI, *An hp-version of the $C^0$-continuous petrov-Galerkin time stepping method for nonlinear second-order initial value problems*, Adv. Comp. Math., 46 (2020), pp. No. 56, 25 pp.

[35] T. WERDER, K. GERDES, D. SCHÖTZAU, AND C. SCHWAB, *hp-discontinuous Galerkin time stepping for parabolic problems*, Comput. Methods Appl. Mech. Engrg., 190 (2001), pp. 6685–6708.

[36] M. ZHANG, W. YAN, F. JING, AND H. ZHAO, *Discontinuous Galerkin method for the diffusive-viscous wave equation*, Appl. Numer. Math., 183 (2023), pp. 118–139.