# A high-order discontinuous Galerkin approximation to ordinary differential equations with applications to elastodynamics

A high-order discontinuous Galerkin approximation to ordinary differential equations with... Abstract The aim of this work is to propose and analyse a new high-order discontinuous Galerkin finite element method for the time integration of a Cauchy problem for second-order ordinary differential equations. These equations typically arise after space semidiscretization of second-order hyperbolic-type differential problems, e.g., wave, elastodynamics and acoustics equations. After introducing the new method, we analyse its well-posedness and prove a priori error estimates in a suitable (mesh-dependent) norm. Numerical results are also presented to verify our theoretical estimates. 1. Introduction In this article, we develop a high-order discontinuous Galerkin finite element method for the numerical approximation of ordinary differential equations that arise after space semidiscretization of second-order hyperbolic problems. The applications we have in mind include, e.g., acoustic, elastic and electromagnetic wave propagation phenomena. Traditional approaches for the numerical integration of (second-order) ordinary differential systems rely on implicit and explicit finite difference, Runge–Kutta and Newmark schemes (see, e.g., LeVeque, 2007; Quarteroni et al., 2007; Butcher, 2008, for a comprehensive review). In many engineering applications, explicit methods are in general preferred to implicit ones. In fact, although the latter are unconditionally stable, the former are computationally cheaper. The main drawback of explicit methods is the time step limitation imposed by the Courant–Friedrichs–Lewy (CFL) condition. Such a constraint, which depends in general on the space discretization parameters and the media properties, can severely affect the computational efficiency. A possible way to alleviate this limitation is to introduce suitable local time-stepping (LTS) algorithms (Collino et al., 2003; Diaz & Grote, 2009; Grote & Mitkova, 2013), using a small time step only when needed. Another possibility is to adopt an explicit LTS method by extending the so-called arbitrary high-order derivatives discontinuous Galerkin (ADER-DG) approach (Dumbser et al., 2007; Taube et al., 2009). In this context, a proper time step can be tailored for each element of the time mesh. However, to correctly propagate the wave field from one element to the other, an additional (computationally demanding) synchronization process has to be operated. Differently than above, here, we derive an implicit arbitrarily high-order accurate time integration scheme based on a DG approach. DG methods (Reed & Hill, 1973; Lesaint & Raviart, 1974) have been firstly proposed to approximate (in space) hyperbolic problems (Reed & Hill, 1973) and then generalized to elliptic and parabolic equations (Wheeler, 1978; Arnold, 1982); see also Arnold et al. (2001–2002), Houston et al. (2000), Cockburn & Shu (1998), Rivière (2008), Hesthaven & Warburton (2008) and Di Pietro & Ern (2012). Relevant applications and analysis of DG schemes for the scalar wave equation can be found in Rivière & Wheeler (2003), Grote et al. (2006) and Baccouch (2012), whereas for elastodynamics problems we refer the reader to Dumbser et al. (2007), Wilcox et al. (2010), Antonietti et al. (2012, 2016a,b,c, 2017), Mazzieri et al. (2013), Delcourte & Glinsky (2015) and Ferroni et al. (2017). The DG approach has also been used to solve initial-value problems. In time-dependent problems, the information follows the positive direction of time and solutions are casual (they depend on the past but not on the future). In contrast with finite difference time integration schemes, for which the solution at the current step depends upon the previous steps, time-discontinuous Galerkin methods applied over time slabs $$[t_n, t_{n+1}]$$ lead to a casual system in which the solution at the current time slab depends only upon the solution at $$t_n^-$$. By coupling discontinuous Galerkin discretizations in both space and time leads to a fully space-time finite element formulation. Relevant works on this topic concern both parabolic and hyperbolic problems (see, e.g., Delfour et al., 1981; Werder et al., 2001; van der Vegt et al., 2006). For the latter, space-time finite elements are typically built upon reformulating the original problem as a system of first-order equations (see, e.g., Hughes & Hulbert, 1988; Johnson, 1993; Banks et al., 2014). The latter can be seen as the result of space semidiscretization of first-order hyperbolic problems or even second-order hyperbolic problems in which the problem is formulated in terms of the displacement (respectively velocity) and the stress (respectively strain) tensor fields. To the best of our knowledge, only a few recent results regarding finite element approximations of second-order differential systems are available in the literature (Hulbert & Hughes, 1990; Thompson & Pinsky, 1996; Adjerid & Temimi, 2011; Yang et al., 2012; Walkington, 2014). In Adjerid & Temimi (2011), a new DG approach based on the solution of the scalar wave equation (and higher-order differential equations) has been proposed and analysed. The stabilization terms appearing in Adjerid & Temimi (2011) introduced to penalize the jumps of the solution and its derivative across different time slabs are similar to those proposed in Hughes & Hulbert (1988), Hulbert & Hughes (1990) and Thompson & Pinsky (1996), where a Galerkin least square (GLS) approach is applied to stabilize the numerical scheme and prove its convergence. As an extension of the space-time formulation of Hughes & Stewart (1996), in Yang et al. (2012), the authors present an enriched version of the space-time finite element method to incorporate in the same model multiple temporal scale features. A combination of continuous and discontinuous Galerkin time-stepping approach is used in Walkington (2014) to develop arbitrary order approximation for second-order hyperbolic problems. Stability, convergence and accuracy are proved for scalar wave propagation with nonhomogeneous boundary conditions. In this work, a new DG method for the solution of systems of second-order ordinary differential equations is presented. The resulting weak formulation is obtained by imposing the continuity of tractions and velocities across time slabs weakly, without adding any extra GLS stabilization term. We show that the proposed formulation, in which the displacement field is the only unknown, is well posed and we prove a priori stability and error estimates in a suitable mesh-dependent norm. The obtained time-discontinuous scheme results in an implicit and unconditionally stable method. Moreover, allowing independent displacement interpolations between different time slabs, this method is naturally suited for an adaptive choice of the time discretization parameters, i.e., the use of high-order polynomials/small time steps only when the solution features a sharp temporal derivative. The article is organized as follows. In Section 2, we formulate the problem, discretize it and analyse its well-posedness. Then, we derive the corresponding algebraic system of equations. In Section 3, we carry out the convergence analysis providing suitable stability and error estimates. The application of the proposed method to the elastodynamics equations is described in Section 4. Here, the space-time finite element formulation is obtained combining the discontinuous Galerkin spectral element (DGSE) spatial discretization proposed in Antonietti et al. (2012) to the one presented here for the time integration. Numerical results are shown in Section 5. Throughout the article, we denote by $$\|\mathbf{a}\|$$ the Euclidean norm of a vector $$\mathbf{a} \in \mathbb{R}^d$$, $$d \geq 1$$ and by $$\|A\|_{\infty} = \max_{i=1, \dots, m} \sum_{j = 1}^{n}|a_{ij}|$$, the $$\ell^\infty$$-norm of a matrix $$A \in \mathbb{R}^{m\times n}$$, $$m,n \geq 1$$. Moreover, $$C$$ denotes a generic positive constant that may take different values in different places, but is always independent of the discretization parameters. The notation $$x \lesssim y$$ means $$x \leq C y$$ for a constant $$C$$ as before. For a given $$I \subset \mathbb{R}$$, for any $$v: I \rightarrow \mathbb{R}$$ we denote by $$L^p(I)$$ and $$H^p(I)$$, $$p \in \mathbb{N} \setminus 0$$ the usual Lebesgue and Hilbert spaces, respectively, and endow them with the usual norms (see Adams & Fournier, 2003). For $$p=0$$, we write $$L^2(I)$$ in place of $$H^0(I)$$. Finally, we use boldface type for vectorial functions. More precisely, the Lebesgue and Hilbert spaces of vector-valued functions are denoted by $${\bf L}^{p}(I)=[L^{p}(I)]^{d}$$ and $${\bf H}^{p}(I)=[H^{p}(I)]^{d}$$, respectively, $$d \geq 1$$. 2. A model problem and its discontinuous Galerkin approximation In this section, we introduce a high-order discontinuous Galerkin finite element method for second-order ordinary differential equations, prove its well-posedness and provide its algebraic form. 2.1 Problem statement and its DG discretization For $$T > 0$$, we consider the following model problem: find $$\mathbf{u}: (0,T] \rightarrow \mathbb{R}^d, d \geq 1$$, such that $$\ddot{\mathbf{u}}(t) + L\dot{\mathbf{u}}(t) + K{\mathbf{u}}(t) = \mathbf{f}(t) \qquad \forall t \in (0, T],$$ (2.1) where $$L, K \in \mathbb{R}^{d\times d}$$ are symmetric, positive definite matrices and $$\mathbf{f} \in {\bf{L}}^2(0,T]$$. We supplement problem (2.1) with the following initial conditions: \begin{align}\mathbf{u}(0) = \mathbf{u}_{0}, \qquad \dot{\mathbf{u}}(0) =\mathbf{u}_{1},\end{align} (2.2) where $$\mathbf{u}_0, \mathbf{u}_1 \in \mathbb{R}^d$$. Problem (2.1) is well posed and admits a unique solution $$\textbf{u} \in {\bf{H}}^2(0,T]$$ in the interval $$(0,T]$$ (see, e.g., Kroopnick, 1999). We partition the interval $$I = (0, T]$$ into $$N$$ time slabs $$I_n = (t_{n-1}, t_n]$$ having length $${\it{\Delta}} t_n = t_n - t_{n-1}$$, for $$n = 1,..,N$$, with $$t_0 = 0$$ and $$t_N=T$$, as it is shown in Fig. 1. Fig. 1. View largeDownload slide Example of time domain partition (top). Zoom of the time domain partition: values $$t_n^+$$ and $$t_n^-$$ are also reported (bottom). Fig. 1. View largeDownload slide Example of time domain partition (top). Zoom of the time domain partition: values $$t_n^+$$ and $$t_n^-$$ are also reported (bottom). In the following we will use the notation: \begin{equation*}({\bf{u}},{\bf{v}})_{I} = \int_{I} {\bf{u}}(s) \cdot {\bf{v}}(s)\,{\rm d}s, \quad \quad \langle {\bf{u}}, {\bf{v}} \rangle_t = {\bf{u}}(t) \cdot {\bf{v}}(t), \end{equation*} where $$\mathbf{a} \cdot \mathbf{b}$$ indicates the Euclidean scalar product between two vectors $$\mathbf{a}, \mathbf{b} \in \mathbb{R}^d$$. To deal with discontinuous functions, we also define, for (a regular enough) $$\mathbf{v}$$, the jump operator at $$t_n$$ as \begin{align*}[\mathbf{v}]_{n} &= \mathbf{v}(t_{n}^+) - \mathbf{v}(t_{n}^-) \quad {\rm for} \; n \geq 0, \end{align*} where \begin{align*}\mathbf{v}(t_{n}^\pm) &= \lim_{\epsilon \rightarrow 0^\pm} \mathbf{v}(t_{n} + \epsilon) \quad {\rm for} \; n \geq 0. \end{align*} Implicit with the above definition is that \begin{align*}[\mathbf{v}]_{0} = \mathbf{v}(0^+) - \mathbf{u}_0,\qquad [\dot{\mathbf{v}}]_{0} = \dot{\mathbf{v}}(0^+) - \mathbf{u}_1. \end{align*} Moreover, we use the symbols $$\mathbf{v}_{n}^+ = \mathbf{v}(t_{n}^+)$$ and $$\mathbf{v}_{n}^- = \mathbf{v}(t_{n}^-)$$ to represent the trace of (a regular enough) $$\mathbf{v}$$, taken within the interior of $$I_{n+1}$$ and $$I_{n}$$, respectively (cf. Fig. 1). Next, following a time integration approach, we incrementally build (on $$n$$) an approximation of the exact solution $$\mathbf{u}$$ on each time slab $$I_n$$. For this reason, we focus on the generic interval $$I_n$$, and assume the solution on $$I_{n-1}$$ to be known. Note that $$\textbf{u} \in {\bf{H}}^2(I_n)$$. If we multiply equation (2.1) by $$\dot{\textbf{v}}(t)$$, being $${\mathbf{v}}(t)$$ a regular enough function, we obtain $$(\ddot{\mathbf{u}},\dot{\mathbf{v}})_{I_n} + (L \dot{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n} + (K{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n} = ( \mathbf{f} , \dot{\mathbf{v}})_{I_n}.$$ (2.3) Finally, since $$\mathbf{u} \in {\bf{H}}^2(0,T)$$, we observe that $$[\mathbf{u}]_n = [\dot{\mathbf{u}}]_n = \textbf{0}$$, for $$n=1, \dots, N$$, we rewrite (2.3) adding suitable (strongly consistent) terms $$(\ddot{\mathbf{u}},\dot{\mathbf{v}})_{I_n} + (L \dot{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n} + (K{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n} + [\dot{\mathbf{u}}]_{n-1} \cdot \dot{\mathbf{v}}_{n-1}^+ + K [\mathbf{u}]_{n-1} \cdot \mathbf{v}_{n-1}^+ = ( \mathbf{f} , \dot{\mathbf{v}})_{I_n}.$$ (2.4) Next, we introduce the local finite dimensional space \begin{align*}{{\bf{V}}_n^r} = \left\{\mathbf{v}: I_n \rightarrow \mathbb{R}^d: \mathbf{v} \in [\mathbb{P}^r(I_n)]^d \right\}\!, \end{align*} where $$\mathbb{P}^r(I_n)$$ is the space of polynomials of degree not greater than or equal to $$r \geq 2$$ on $$I_n$$. Then, introducing $$\textbf{r}=(r_1,...,r_N) \in \mathbb{N}^N$$ the polynomial degree vector, we can define the DG finite element space as \begin{align*}{\mathcal{V}^\bf{r}} = \left\{\mathbf{v} \in {\bf{L}}^2(0,T): \mathbf{v}|_{I_n} \in {{\bf{V}}_n^{r_n}} \forall \, n = 1, \dots, N \right\}\!, \end{align*} whose dimension is $$\sum_{n=1}^N (r_n+1)d$$. Summing over all time slabs in (2.4), we are able to define the bilinear form $${\mathcal A} : {{\mathcal{V}^\bf{r}} \times {\mathcal{V}^\bf{r}}} \rightarrow \mathbb{R}$$ \begin{align}{\mathcal A}(\mathbf{u}, \mathbf{v}) & = \sum_{n=1}^N \Big( (\ddot{\mathbf{u}},\dot{\mathbf{v}})_{I_n} + (L \dot{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n} + (K{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n}\Big) + \sum_{n=1}^{N-1} \left( \dot{[\mathbf{u}]}_{n} \cdot \dot{\mathbf{v}}_{n}^+ + K [\mathbf{u}]_{n} \cdot \mathbf{v}_{n}^+ \right)\nonumber\\ &\quad{} + \dot{\mathbf{u}}_{0}^+ \cdot \dot{\mathbf{v}}_{0}^+ + K \mathbf{u}_{0}^+ \cdot \mathbf{v}_{0}^+,\end{align} (2.5) and the linear functional $$F : {{\mathcal{V}^\bf{r}}} \rightarrow \mathbb{R}$$ as \begin{align}F (\mathbf{v}) & = \sum_{n=1}^N ( \mathbf{f} , \dot{\mathbf{v}} )_{I_n} + {\mathbf{u}}_{1} \cdot \dot{\mathbf{v}}_{0}^+ + K \mathbf{u}_{0} \cdot \mathbf{v}_{0}^+,\end{align} (2.6) respectively. Notice that for $$n = 1$$, we implicitly adopted the convention that $$\mathbf{u}_0^- = \mathbf{u}_0$$ and $$\dot{\mathbf{u}}_0^- = \mathbf{u}_1$$. Moreover, we point out that the definition of the bilinear form $${\mathcal A}(\cdot, \cdot)$$ makes sense whenever its arguments are, at least, $$H^2(I_n)-$$ functions for any $$n=1,\ldots,N$$. The DG formulation of problem (2.1)–(2.2) reads as follows: find $$\mathbf{u}_{\textrm{DG}} \in {\mathcal{V}^\bf{r}}$$ such that \begin{align}{\mathcal A}(\mathbf{u}_{\small{DG}}, \mathbf{v}) = F(\mathbf{v}) \qquad \forall \, \mathbf{v}\in {\mathcal{V}^\bf{r}}.\end{align} (2.7) The existence and uniqueness of the discrete solution $$\textbf{u}_{\textrm{DG}} \in {\mathcal{V}^\bf{r}}$$ of problem (2.7) is a direct consequence of the following result. Proposition 2.1 The function $$\left\||{{\cdot}}\right\|| : {\mathcal{V}^\bf{r}} \rightarrow \mathbb{R}^+$$, defined as \begin{align}\left\|| {{\mathbf{v}}} \right\||^2 & = \sum_{n=1}^N \|L^\frac{1}{2} \dot{\mathbf{v}}\|_{{\bf{L}}^2(I_n)}^2 + \frac{1}{2} (\dot{\mathbf{v}}_{0}^+)^2 + \frac{1}{2}\sum_{n=1}^{N-1} \left(\dot{[\mathbf{v}]}_{n} \right)^2 + \frac{1}{2}(\dot{\mathbf{v}}_{T}^-)^2 \nonumber\\ &\quad{} + \frac{1}{2} (K^\frac{1}{2}{\mathbf{v}}_{0}^+)^2 + \frac{1}{2} \sum_{n=1}^{N-1} \left(K^\frac{1}{2} [\mathbf{v}]_{n} \right)^2 +\frac{1}{2} (K^\frac{1}{2}{\mathbf{v}}_{T}^-)^2,\end{align} (2.8) is a norm on $${\mathcal{V}^\bf{r}}$$. Proof Since the homogeneity and the subadditivity properties are satisfied, we just need to show that \begin{align*}\left\||{\mathbf{v}}\right\|| = 0 \Leftrightarrow \mathbf{v} = \mathbf{0}. \end{align*} The sufficient condition is trivial. We then show that $$\left\||{\mathbf{v}}\right\||= 0$$ implies $$\mathbf{v} = \textbf{0}$$. Clearly, if $$\left\||{\mathbf{v}}\right\|| = 0$$ then all the terms on the right-hand side of (2.8) are null. In particular, from the facts that $$K$$ and $$L$$ are positive definite and $$\|L^{\frac{1}{2}}\dot{\mathbf{v}} \|_{{\bf{L}}^2(I_0)} = 0$$ and $$( K^{\frac{1}{2}} \mathbf{v}_0^+ )^2 = 0$$, we conclude that $$\mathbf{v}$$ is such that \begin{align*}\begin{cases}\dot{\mathbf{v}}(t) = \mathbf{0} \qquad \forall \,t \in I_0 ,\\ \mathbf{v}_0^+ = \mathbf{0}. \end{cases}\end{align*} Therefore, we have that $$\mathbf{v}(t)|_{I_0} = \mathbf{0}$$. We now proceed by induction, and consider the interval $$I_n$$, supposing $$\mathbf{v}(t)|_{I_{n-1}} = \mathbf{0}$$. From $$\left\||{\mathbf{v}}\right\||= 0$$ we have $$(K^\frac{1}{2} [\mathbf{v}]_{n} )^2 = 0$$, which yields $$\mathbf{v}_{n}^+ = 0$$. Then, we infer that $$\mathbf{v}$$ is such that \begin{align*}\begin{cases}\dot{\mathbf{v}}(t) = \mathbf{0} \qquad \forall \,t \in I_n ,\\ \mathbf{v}_n^+ = \mathbf{0}. \end{cases}\end{align*} As a result, we conclude that $$\mathbf{v}$$ is the null function on any interval $$I_n$$, $$n = 0, \dots, N-1$$, and the proof is complete. □ The following result follows from Proposition 2.1. Remark 2.2 Taking $$\mathbf{u} = \mathbf{v}$$ in (2.5) and integrating by parts, we obtain \begin{equation*}{\mathcal A}(\mathbf{v}, \mathbf{v}) = \left\||{ \mathbf{v}}\right\||^2 \quad \forall \, \mathbf{v}\in {\mathcal{V}^\bf{r}}, \end{equation*} i.e., the bilinear form $${\mathcal A}(\cdot, \cdot)$$ defined in (2.5) is coercive with respect to the norm $$\left\||{\cdot}\right\||$$, with coercivity constant $$\alpha = 1$$. Therefore, the following result holds. Proposition 2.3 Problem (2.7) admits a unique solution $$\textbf{u}_{\textrm{DG}} \in {\mathcal{V}^\bf{r}}$$. Proof The thesis follows directly from Proposition 2.1, the bilinearity of $${\mathcal A}(\cdot, \cdot)$$ and the linearity of $$F(\cdot)$$. □ 2.2 Algebraic formulation We derive here the algebraic formulation corresponding to problem (2.4) for the time slab $$I_n$$, $$n\geq 1$$, where a local degree $${r_n}$$ is employed. We remark that the employment of discontinuous functions between a node $$t_n$$ allows to compute the solution of the problem separately for one time slab at a time. Indeed, the weak formulation (2.7) restricted to $${{\bf{V}}_n^{{r_n}}}$$ reads as: find $${\mathbf{u}_n} = {\textbf{u}_{\textrm{DG}}}_{|_{I_n}} \in {{\bf{V}}_n^{{r_n}}}$$ such that \begin{align*}& (\ddot{\bf u}_n,\dot{\mathbf{v}})_{I_n} + (L \dot{\bf u}_n, \dot{\mathbf{v}} )_{I_n} + (K{{{\mathbf{u}_n}}}, \dot{\mathbf{v}} )_{I_n} + \langle \dot{\bf u}_n, \dot{\mathbf{v}}\rangle_{t_{n-1}^+}+ \langle K {{\mathbf{u}_n}} , \mathbf{v}\rangle_{t_{n-1}^+}\\ &\quad{} = ( \mathbf{f} , \dot{\mathbf{v}})_{I_n} + \dot{\bf u}_n(t_{n-1}^-) \cdot \dot{\mathbf{v}}_{n-1}^+ + K {{\mathbf{u}_n}}(t_{n-1}^-) \cdot \mathbf{v}_{n-1}^+ \quad \forall \; \mathbf{v} \in {{\bf{V}}_n^{{r_n}}}, \end{align*} where on the right-hand side the values $$\dot{\bf u}_n(t_{n-1}^-)$$ and $${{\mathbf{u}_n}}(t_{n-1}^-)$$ computed for $$I_{n-1}$$ are used as initial conditions for the current slab. Notice that, for the slab $$I_1$$, we set $$\dot{\bf u}_1(t_{0}^-) = {\bf{u}}_1$$ and $${\bf u}_1(t_{0}^-) = {\bf{u}}_0$$. Focusing on the generic interval $$I_n$$, we introduce a basis $$\{\psi^\ell(t) \}_{\ell=1, \dots, {r_n}+1}$$ for the polynomial space $$\mathbb{P}^{r_n}_n(I_n)$$ and define $$D= d({r_n} + 1)$$ the dimension of the local finite dimensional space $${{\bf{V}}_n^{{r_n}}}$$. We also introduce the vectorial basis $$\{ {\bf {\it{\Psi}}_i^\ell}(t) \}_{i=1, \dots, d}^{\ell=1, \dots, {r_n}+1}$$, where $${\bf {\it{\Psi}}_i^\ell}(t)$$ is the $$d$$-dimensional vector whose $$i$$th component is $$\psi^\ell$$ and other components are zero. Using the notation just introduced, we can write the trial function $${\mathbf{u}_n}$$ as a linear combination of the basis functions, i.e., \begin{align*}{\mathbf{u}_n}(t) = \sum\limits_{j=1}^d \sum\limits_{m=1}^{{r_n}+1} \alpha_j^m {\bf {\it{\Psi}}_j^m}(t), \end{align*} where $$\alpha_j^m \in \mathbb{R}$$ for $$j = 1, \dots, d$$ and $$m=1, \dots, {r_n}+1$$. Next, we write equation (2.7) for any test function $${\bf {\it{\Psi}}_i^\ell}(t), i = 1, \dots, d, \ell = 1, \dots, {r_n}+1$$, obtaining the algebraic system \begin{align*}A \mathbf{U}_n = \mathbf{F}_n, \end{align*} where on the interval $$I_n$$, $$\mathbf{U}_n \in \mathbb{R}^D$$ is the solution vector, $$\mathbf{F}_n \in \mathbb{R}^D$$ corresponds to the data and is given componentwise as \begin{equation*}(\mathbf{F}_n)_i^\ell = ( \mathbf{f} , \dot{\bf {\it{\Psi}}}_i^\ell )_{I_n} + \dot{\bf u}_n(t_{n-1}^-) \cdot \dot{\bf {\it{\Psi}}}_i^\ell(t_{n-1}^+) + K {{\mathbf{u}_n}}(t_{n-1}^-) \cdot \bf{{\it{\Psi}}}_i^\ell(t_{n-1}^+) \quad {\rm for } \; i = 1, \dots, d, \ell = 1, \dots, {r_n}+1, \end{equation*} and $$A \in \mathbb{R}^{D\times D}$$ is the related local stiffness matrix. We next investigate the structure of the matrix $$A$$ by defining the following local matrices for $$\ell, m = 1, \dots, {r_n}+1$$, \begin{align*}& M^1_{\ell m} = ( \ddot{\psi}^m , \dot{\psi}^\ell )_{I_n}, \qquad M^{2}_{\ell m} = (\dot{\psi}^m ,\dot{\psi}^\ell )_{I_n}, \qquad M^3_{\ell m} = ( {\psi}^m, \dot\psi^\ell )_{I_n}, \\ & M^4_{\ell m} = \langle \dot{\psi}^m, \dot{\psi}^\ell \rangle_{t_{n-1}^+}, \quad M^{5}_{\ell m} = \langle \psi^m, \psi^\ell \rangle_{t_{n-1}^+}. \end{align*} Setting \begin{align*}M & = M^1 + M^4, \\ B_{ij} & = L_{ij} M^{2} + K_{ij} \left( M^3 + M^{5} \right)\!, \qquad i,j=1,...,d, \end{align*} with $$M, B_{ij} \in \mathbb{R}^{({r_n}+1)\times({r_n}+1)}$$ for any $$i,j = 1, \dots,d$$, we can rewrite the matrix $$A$$ as \begin{align}A = \begin{vmatrix}M & 0 & \cdots & 0\\ 0 & M & & \vdots \\ \vdots & & \ddots & 0 \\ 0 & \cdots & 0 & M \end{vmatrix}+ \begin{vmatrix}B_{11} & B_{12} & \cdots & B_{1d}\\ B_{21} & B_{22} & & \vdots \\ \vdots & & \ddots & \\ B_{d1} & \cdots & & B_{dd}\end{vmatrix}.\end{align} (2.9) Its sparsity clearly depends on the sparsity pattern of the matrices $$L$$ and $$K$$. 3. Convergence analysis In this section, we first analyse the stability of the DG method (2.7) and then we prove a priori error estimates. 3.1 Stability analysis Proposition 3.1 Let $$\mathbf{f} \in \textbf{L}^2(0,T]$$ and $$\mathbf{u}_0,\mathbf{u}_1 \in \mathbb{R}^d$$. Then, the solution of (2.7), $${\mathbf{u}_{\textrm{DG}}} \in {\mathcal{V}^\bf{r}}$$, satisfies $$\left\||{ {\mathbf{u}_{\textrm{DG}}} }\right\||\lesssim \left( \|L^{-\frac{1}{2}}\mathbf{f}\|_{{\bf{L}}^2(0,T)}^2 + (K^\frac{1}{2} \mathbf{u}_0 )^2 + (\mathbf{u}_1 )^2 \right)^{\frac{1}{2}}.$$ (3.1) Proof Using the definition (2.8) and arithmetic–geometric inequalities, we obtain \begin{align*}\left\||{{\mathbf{u}_{\textrm{DG}}}}\right\||^2 & = {\mathcal A}({\mathbf{u}_{\textrm{DG}}},{\mathbf{u}_{\textrm{DG}}}) = \sum_{n=1}^N ( \mathbf{f} , \dot{\mathbf{u}}_{\textrm{DG}} )_{I_n} + {\mathbf{u}}_{1} \cdot \dot{\mathbf{u}}_{\textrm{DG}}(0^+) + K \mathbf{u}_{0} \cdot {\mathbf{u}}_{\textrm{DG}}(0^+) \\ & \leq \frac{1}{2}\sum_{n=1}^N \left( \| L^{-\frac{1}{2}} \mathbf{f}\|_{{\bf{L}}^2(I_n)}^2 + \| L^{\frac{1}{2}} \dot{\mathbf{u}}_{\textrm{DG}} \|_{{\bf{L}}^2(I_n)}^2 \right) \\ & \quad + {\mathbf{u}}_{1}^2 + \frac{1}{4} \dot{\mathbf{u}}_{\textrm{DG}}^2(0^+) + (K^\frac{1}{2} \mathbf{u}_{0})^2 + \frac{1}{4} \left(K^\frac{1}{2}{\mathbf{u}}_{\textrm{DG}}(0^+)\right)^2, \end{align*} i.e., \begin{align*}\frac{1}{2}\left\||{ \textbf{u}_{\textrm{DG}}^2}\right\|| \leq \frac{1}{2}\sum_{n=1}^N \| L^{-\frac{1}{2}} \mathbf{f}\|_{{\bf{L}}^2(I_n)}^2 + {\mathbf{u}}_{1}^2 + (K^\frac{1}{2} \mathbf{u}_{0})^2. \end{align*} □ 3.2 Error estimates In this section, we derive a priori error estimates in the $$\left\||{{\cdot}}\right\||$$ norm defined in (2.8). We start by introducing some preliminary definitions and results. Definition 3.2 (Schötzau & Schwab, 2000, Definition 3.1) Let $$I = (-1, 1)$$. For a function $$w \in L^2(I)$$ which is continuous at $$t =1$$ we define the $$L^2$$-projection (with boundary preserving value) $${\mathcal{P}}^r w \in \mathbb{P}^r(I)$$, $$r \geq 0$$, by the $$r + 1$$ conditions \begin{align}({\mathcal{P}}^r w - w)(1) &= 0, \\ \end{align} (3.2a) \begin{align}\int_I (w -{\mathcal{P}}^r w) q \,{\rm d}t &= 0 \quad \forall q \in \mathbb{P}^{r-1}(I).\end{align} (3.2b) For $$r =0$$, only the first condition (3.2a) is necessary. $${\mathcal{P}}^rw$$ is well posed and satisfies optimal error estimates, see Schötzau & Schwab (2000, Theorem 3.10). Definition 3.3 Let $$I = (-1, 1)$$. For a function $$u \in H^1(I)$$ such that $$u'$$ is continuous at $$t=1$$, we define $${\it{\Pi}}^r u \in \mathbb{P}^r(I)$$, $$r \geq 1$$, by $${\it{\Pi}}^r u(t) = u(-1) + \int_{-1}^t {\mathcal{P}}^{r-1}u'(s)\,{\rm d}s \quad \forall t \in (-1,1),$$ (3.3) where $${\mathcal{P}}^{r-1}u'$$ is given in Definition 3.2. In the following, we denote by $$\{L_i\}_{i\geq 0}$$, $$L_i \in \mathbb{P}^i(I)$$, the Legendre polynomials on $$I = (-1, 1)$$, defined by $$L_0(t) = 1$$, $$L_1(t) = t$$ and $$\int_{-1}^t L_i(s)\,{\rm d}s =\frac{1}{2i+1}(L_{i+1} (t) - L_{i-1}(t))$$ for $$i \geq 1$$, see Canuto et al. (2006). Lemma 3.4 The projection $${\it{\Pi}}^r u$$ in Definition 3.3 is well defined. Proof We consider the Legendre series $$u' = \sum_{i=0}^\infty b_i L_i$$ with coefficients $$b_i \in \mathbb{R}$$. From Schötzau & Schwab (2000, Lemma 3.5), we find that $${\mathcal{P}}^{r-1}u'(t) = \sum_{i=0}^{r-2} b_i L_i(t) + \left( \sum_{i=r-1}^{\infty} b_i \right)L_{r-1}(t)$$ (3.4) and consequently, $${\it{\Pi}}^r u(t) = u(-1) + \sum_{i=0}^{r-2} b_i \int_{-1}^t L_i(s) \,{\rm d}s + \left( \sum_{i=r-1}^{\infty} b_i \right)\int_{-1}^t L_{r-1}(s) \,{\rm d}s.$$ (3.5) Then, \begin{align*}{\it{\Pi}}^r u(t) & = \left( u(-1) + b_0 -\frac{b_1}{3}\right)L_0(t) + \left(b_0 - \frac{b_2}{5} \right)L_1(t) + \sum_{i=2}^{r-3} \left(\frac{b_{i-1}}{2i-1} - \frac{b_{i+1}}{2i+3} \right)L_i(t) \\ &\quad{} + \left(\frac{b_{r-3}}{2r-5} - \sum_{i=r-1}^{\infty}\frac{b_{i}}{2r-1} \right)L_{r-2}(t) + \frac{b_{r-2}}{2r-3}L_{r-1}(t) + \sum_{i=r-1}^{\infty}\frac{b_i}{2r-1}L_r(t), \end{align*} or, equivalently, $${\it{\Pi}}^r u (t) = \sum_{i=0}^r u_i^* L_i(t),$$ (3.6) with \begin{eqnarray}u_0^* & = & u(-1) +b_0 -\frac{b_1}{3}, \\ \end{eqnarray} (3.7) \begin{eqnarray}u_1^* & = & b_0 - \frac{b_2}{5}, \\ \end{eqnarray} (3.8) \begin{eqnarray}u_i^* & = &\frac{b_{i-1}}{2i-1} - \frac{b_{i+1}}{2i+3} \quad {\rm for \;} i = 2,...,r-3,\\ \end{eqnarray} (3.9) \begin{eqnarray}u_{r-2}^* & = & \frac{b_{r-3}}{2r-5} - \sum_{i=r-1}^{\infty}\frac{b_{i}}{2r-1},\\ \end{eqnarray} (3.10) \begin{eqnarray}u_{r-1}^* & = & \frac{b_{r-2}}{2r-3},\\ \end{eqnarray} (3.11) \begin{eqnarray}u_r^* & = & \sum_{i=r-1}^{\infty}\frac{b_i}{2r-1} . \end{eqnarray} (3.12) The well-posedness of $${\it{\Pi}}^r u$$ follows from (3.6) and the orthogonality of the Legendre polynomials $$L_i$$. □ Now, we analyse the properties of the projection $${\it{\Pi}}^r u$$. Lemma 3.5 Let $$I = (-1, 1)$$. For any $$u \in H^1(I),$$ such that $$u'$$ is continuous at $$t=1$$ we have, for $$r \geq 1$$, \begin{equation*}\| {\it{\Pi}}^r u \|_{L^2(I)} \lesssim |u(1)| + |u'(1)| + \| u' \|_{L^2(I)}. \end{equation*} Proof We expand $$u'$$ as $$u' = \sum_{i=0}^\infty b_i L_i$$ with coefficients $$b_i \in \mathbb{R}$$; then, $$u$$ can be written as $$u(t) = u(-1) + \sum_{i=0}^\infty b_i \int_{-1}^t L_i(s) \,{\rm d}s$$. By orthogonality of the Legendre polynomials, it easily follows that $$u(1) = u(-1) + 2b_0.$$ (3.13) Next, from (3.7)–(3.12) and (3.13), we have that \begin{align}\| {\it{\Pi}}^r u \|_{L^2(I)}^2 & = \sum_{i=0}^r \frac{2}{2i+1} {u_i^*}^2 \nonumber \\ & = 2\left(u(-1)+b_0 -\frac{b_1}{3}\right)^2 + \frac{2}{3}\left(b_0 -\frac{b_2}{5}\right)^2 + \sum_{j=2}^{r-3}\frac{2}{2j+1} \left(\frac{b_{j-1}}{2j-1} - \frac{b_{j+1}}{2j+3} \right)^2 \nonumber\\ &\quad{} + \frac{2}{2r-3}\left(\frac{b_{r-3}}{2r-5}- \sum_{j=r-1}^{\infty} \frac{b_j}{2r-1}\right)^2 + \frac{2}{2r-1}\left(\frac{b_{r-2}}{2r-3}\right)^2 + \frac{2}{2r+1}\left(\sum_{j=r-1}^{\infty} \frac{b_{j}}{2r-1}\right)^2 \nonumber\\ &\leq 4u(1)^2 + 4b_0^2 + \frac{4}{3}b_1^2 + \frac{4}{3}b_0^2 + \frac{4}{75}b_2^2 + 2 \sum_{j=2}^{r-3}\frac{2}{2j+1}\frac{b_{j-1}^2}{(2j-1)^2} + 2\sum_{j=2}^{r-3}\frac{2}{2j+1}\frac{b_{j+1}^2}{(2j+3)^2} \nonumber \\ & \quad + \frac{4}{2r-3}\frac{b_{r-3}^2}{(2r-5)^2} + \frac{4}{2r-3}\left(\sum_{j=r-1}^{\infty} \frac{b_j}{2r-1}\right)^2 + \frac{2}{2r-1}\left(\frac{b_{r-2}}{2r-3}\right)^2\nonumber \\ &\quad + \frac{2}{2r+1}\left(\sum_{j=r-1}^{\infty} \frac{b_{j}}{2r-1}\right)^2.\end{align} (3.14) Then, we observe that \begin{gather*}\frac{4}{2r-3}\left(\sum_{j=r-1}^{\infty} \frac{b_j}{2r-1}\right)^2 \lesssim \frac{1}{r^3}\left(\sum_{j=r-1}^{\infty} b_j\right)^2 \lesssim |u'(1)|^2, \\ \frac{2}{2r+1}\left(\sum_{j=r-1}^{\infty} \frac{b_j}{2r-1}\right)^2 \lesssim \frac{1}{r^3}\left(\sum_{j=r-1}^{\infty} b_j\right)^2 \lesssim |u'(1)|^2,\\ 4b_0^2 + 2\sum_{j=2}^{r-3}\frac{2}{2j+1}\frac{b_{j-1}^2}{(2j-1)^2} + \frac{4}{2r-3}\frac{b_{r-3}^2}{(2r-5)^2} + \frac{2}{2r-1}\frac{b_{r-2}^2}{(2r-3)^2} \lesssim \sum_{j=0}^{r-2}\frac{2}{2j+1}b_j^2 \leq \| u' \|_{L^2(I)}^2, \end{gather*} and similarly \begin{equation*}\frac{4}{3}b_0^2 + \frac{4}{3}b_1^2 + \frac{4}{75}b_2^2 + 2\sum_{j=2}^{r-3}\frac{2}{2j+1}\frac{b_{j+1}^2}{(2j+3)^2} \lesssim \| u' \|_{L^2(I)}^2. \end{equation*} Using the above estimates into (3.14) yields $\| {\it{\Pi}}^r u \|_{L^2(I)}^2 \lesssim |u(1)|^2 + |u'(1)|^2 + \| u' \|_{L^2(I)}^2.$ □ Lemma 3.6 Under the same assumptions on $$u$$ made in Lemma 3.5 it holds \begin{align}({\it{\Pi}}^r u - u)(-1) &= 0,\\ \end{align} (3.15a) \begin{align}({\it{\Pi}}^r u - u)(1) &= 0, \\ \end{align} (3.15b) \begin{align}({\it{\Pi}}^r u - u)'(1) &= 0,\\ \end{align} (3.15c) \begin{align}\int_I (u -{\it{\Pi}}^r u)' q {\,\rm d}t &= 0 \quad \forall q \in \mathbb{P}^{r-2}(I).\end{align} (3.15d) Proof From Definition 3.3 we obtain (3.15a). Property (3.15b) can be proved as follows: \begin{equation*}{\it{\Pi}}^r u(1) = u(-1) + \sum_{i=0}^{r-2} b_i \int_{-1}^1 L_i(s) {\,\rm d}s + \left( \sum_{i=r-1}^{\infty} b_i \right)\int_{-1}^1 L_{r-1}(s) {\,\rm d}s = u(-1) + 2b_0 = u(1). \end{equation*} Equality (3.15c) follows by taking the derivative with respect to $$t$$ in (3.5), evaluating the latter at $$t=1$$ and using (3.4). Finally, (3.15d) comes from \begin{equation*}\int_I (u -{\it{\Pi}}^r u)' L_k{\,\rm d}t = \int_I (u' -{\mathcal{P}}^{r-1} u') L_k {\,\rm d}t = \int_I \left(\sum_{i=r}^{\infty} b_i L_i(t) - \left( \sum_{i=r}^{\infty} b_i \right)L_{r-1}\right)L_k{\,\rm d}t = 0, \end{equation*} for all $$k=0,...,r-2$$. □ On an arbitrary interval $$I = (a, b)$$ with $${\it{\Delta}} t = b - a > 0$$, we define $${\it{\Pi}}^r_{I}$$ via the linear map $${T} : (-1, 1) \rightarrow (a, b)$$, $$\xi \rightarrow y = \frac{1}{2} (a + b + \xi {\it{\Delta}} t )$$ as $${\it{\Pi}}^r_{I} u = [{\it{\Pi}}^r (u \circ {T})] \circ {T}^{-1},$$ (3.16) from which it follows that \begin{align}({\it{\Pi}}^r_I u - u)(a) &= 0,\\ \end{align} (3.17a) \begin{align}({\it{\Pi}}^r_I u - u)(b) &= 0, \\ \end{align} (3.17b) \begin{align}({\it{\Pi}}^r_I u - u)'(b) &= 0,\\ \end{align} (3.17c) \begin{align}\int_I (u -{\it{\Pi}}^r_I u)' q {\,\rm d}t &= 0 \quad \forall q \in \mathbb{P}^{r-2}(I).\end{align} (3.17d) Remark 3.7 Reasoning as in the proof of Lemma 3.4, on the generic interval $$I=(a,b)$$, $${\it{\Pi}}^r_I$$ can be written as \begin{equation*}{\it{\Pi}}^r_I u = \sum_{i=0}^{r} \widetilde{u}_i^* \widetilde{L}_i = u(a) + \sum_{i=0}^{r-2} \widetilde{b}_i \int_{a}^t \widetilde{L}_i(s) {\,\rm d}s + \left( \sum_{i=r-1}^{\infty} \widetilde{b}_i \right)\int_{a}^t \widetilde{L}_{r-1}(s) {\,\rm d}s, \end{equation*} where $$\widetilde{L}_i$$ represents the $$i$$th shifted Legendre polynomial (cf. Cohen & Tan, 2012), so that $$u = \sum_{i=0}^{\infty} \widetilde{u}_i \widetilde{L}_i$$, and the coefficients $$\widetilde{u}_{i}^*$$ are defined as in (3.7)–(3.12), provided that $$b_j$$ is replaced by $$\widetilde{b}_j$$, for any $$j$$. We recall the following approximation results from Babuska & Suri (1987), see also Schwab (1998) and Canuto et al. (2007). Proposition 3.8 Let $$I = (a, b)$$ with $${\it{\Delta}} t = b - a > 0$$. For every $$u \in H^s(I)$$ there exist a constant C, independent of $$u,r,{\it{\Delta}} t$$ and a sequence $$\{{Q}^r u\}_{r\geq 1}$$, with each $${Q}^r u \in \mathbb{P}^r (I)$$ such that, for any $$0 \leq q \leq s$$, \begin{align}\| u - {Q}^ru \|_{H^q(I)} & \lesssim \frac{{\it{\Delta}} t^{\mu-q}}{r^{s-q}} \| u \|_{H^s(I)} \quad s\geq 0, \\ \end{align} (3.18) \begin{align}\| u - {Q}^ru \|_{L^2(\partial I)} & \lesssim \frac{{\it{\Delta}} t^{\mu-\frac{1}{2}}}{r^{s-\frac{1}{2}}} \| u \|_{H^s(I)} \quad s> \frac{1}{2}, \\ \end{align} (3.19) \begin{align}\| u - {Q}^ru \|_{H^1(\partial I)} & \lesssim \frac{{\it{\Delta}} t^{\mu-\frac{3}{2}}}{r^{s-\frac{3}{2}}} \| u \|_{H^s(I)} \quad s> \frac{3}{2},\end{align} (3.20) where $$\mu = \min(r+1,s).$$ Next, we define the projector $$\mathcal{P}^{r}_Iu$$ on an arbitrary interval $$I = (a, b)$$ with $${\it{\Delta}} t = b - a > 0$$ via the linear map $${T} : (-1, 1) \rightarrow (a, b)$$, $$\xi \rightarrow y = \frac{1}{2} (a + b + \xi {\it{\Delta}} t )$$ as $$\mathcal{P}^r_{I} u = [\mathcal{P}^r (u \circ {T})] \circ {T}^{-1},$$ (3.21) for which the following approximation result holds, see Schötzau & Schwab (2000). Proposition 3.9 Let $$I = (a, b)$$ with $${\it{\Delta}} t = b - a > 0$$. For any $$u \in H^s(I)$$, $$s \geq 2$$, there holds \begin{align}\| u' - \mathcal{P}^{r-1}_Iu' \|_{L^2(I)} & \lesssim \frac{{\it{\Delta}} t^{\mu-1}}{r^{s-1}} \| u \|_{H^s(I)},\end{align} (3.22) where $$\mu = \min(r+1,s).$$ The following property holds. Lemma 3.10 Let $$I = (a, b)$$ with $${\it{\Delta}} t = b - a > 0$$. For any $$u \in H^s(I)$$, $$s\geq 2$$, let $${\it{\Pi}}^r_{I} u$$, $$r \geq 2$$ be the projector defined as in (3.16). Then, \begin{align}\| (u- {\it{\Pi}}_{I}^r u)' \|_{L^2(I)} & \lesssim \frac{{\it{\Delta}} t^{\mu-1}}{r^{s-1}} \| u \|_{H^s(I)} \\ \end{align} (3.23) \begin{align}\| u- {\it{\Pi}}_{I}^r u \|_{L^2(I)} & \lesssim \frac{{\it{\Delta}} t^{\mu}}{r^{s-1}} \| u \|_{H^s(I)},\end{align} (3.24) where $$\mu = \min(r+1,s).$$ Proof We start by proving (3.23), exploiting the definition (3.3) which directly yields $$({\it{\Pi}}_{I}^r u)' = \mathcal{P}_{I}^{r-1} (u')$$ \begin{align}\| (u- {\it{\Pi}}_{I}^r u)' \|_{L^2(I)} = \| u'- ({\it{\Pi}}_{I}^r u)' \|_{L^2(I)}= \| u'- \mathcal{P}_{I}^{r-1} (u') \|_{L^2(I)} \lesssim \frac{{\it{\Delta}} t^{\mu-1}}{r^{s-1}} \| u \|_{H^s(I)},\end{align} (3.25) where in the last inequality we have used (3.22). As concerns (3.24), employing the definition (3.3), we obtain \begin{align*}\| u- {\it{\Pi}}_{I}^r u \|_{L^2(I)}^2 = & \left\| u(a) + \int_a^tu'(s){\rm d}s - u(a) - \int_a^t \mathcal{P}_{I}^{r-1} (u')(s){\,\rm d}s \right\|_{L^2(I)}^2 \\ = & \left\| \int_a^t u'(s)\,{\rm d}s - \int_a^t \mathcal{P}_{I}^{r-1} (u')(s){\,\rm d}s \right\|_{L^2(I)}^2 = \int_a^b \left( \int_a^t(u'(s) - \mathcal{P}_{I}^{r-1} (u'))(s){\,\rm d}s \right)^2{\,\rm d}t \\ \leq & \int_a^b \left( \int_a^t\big|u'(s) - \mathcal{P}_{I}^{r-1} (u')(s)\big|{\,\rm d}s \right)^2{\,\rm d}t \\ \leq & \int_a^b \left( \left(\int_a^t {\,\rm d}s \right)^{\frac{1}{2}} \left(\int_a^t\big|u'(s) - \mathcal{P}_{I}^{r-1} (u')(s)\big|^2{\,\rm d}s\right)^{\frac{1}{2}} \right)^2{\,\rm d}t\\ = & \int_a^b \left( t-a \right) \left(\int_a^t\big|u'(s) - \mathcal{P}_{I}^{r-1} (u')(s)\big|^2{\,\rm d}s\right) {\,\rm d}t \\ \leq & \int_a^b {\it{\Delta}} t \left(\int_a^b\big|u'(s) - \mathcal{P}_{I}^{r-1} (u')(s)\big|^2{\,\rm d}s\right) {\,\rm d}t = {\it{\Delta}} t^2 \int_a^b \big|u'(s) - \mathcal{P}_{I}^{r-1} (u')(s)\big|^2{\,\rm d}s \\ = & {\it{\Delta}} t^2 \| u'- \mathcal{P}_{I}^{r-1} (u') \|_{L^2(I)}^2 \lesssim {\it{\Delta}} t^2 \frac{{\it{\Delta}} t^{2\mu-2}}{r^{2s-2}} \| u \|_{H^s(I)}^2, \end{align*} from which (3.24) directly follows. □ As a direct consequence of the above Lemma we have the following result. Corollary 3.11 Let $$I= (a,b)$$ with $${\it{\Delta}} t = b-a$$ and let $$u \in H^s(I)$$, with $$s \geq 2$$. There holds \begin{align}\| (u- {\it{\Pi}}^r_{I} u)'' \|_{L^2(I)} & \lesssim \frac{{\it{\Delta}} t^{ \mu - 2}}{r^{s-3}} \| u \|_{H^s(I)},\end{align} (3.26) where $$\mu = \min(r+1,s)$$. Proof From the triangle inequality, interpolation estimates (3.18) and (3.22), we get \begin{align*}\| (u- {\it{\Pi}}^r_{I} u)'' \|_{L^2(I)}&\leq \| (u- {Q}^r u)'' \|_{L^2(I)}+ \| ({Q}^r u- {\it{\Pi}}^r_{I} u)'' \|_{L^2(I)} \\ &\lesssim \frac{{\it{\Delta}} t^{\mu-2}}{r^{s-2}} \| u \|_{H^s(I)}+ \frac{r^2}{{\it{\Delta}} t}\| ({Q}^r u- {\it{\Pi}}^r_{I} u)' \|_{L^2(I)}\\ &= \frac{{\it{\Delta}} t^{\mu-2}}{r^{s-2}} \| u \|_{H^s(I)}+ \frac{r^2}{{\it{\Delta}} t} \| ({Q}^r u)'- {\mathcal{P}}^{r-1}_I u' \|_{L^2(I)}\\ &\lesssim \frac{{\it{\Delta}} t^{\mu-2}}{r^{s-2}} \| u \|_{H^s(I)}+ \frac{r^2}{{\it{\Delta}} t} \| (u - {Q}^r u)'\|_{L^2(I)}+ \frac{r^2}{{\it{\Delta}} t} \|u' - {\mathcal{P}}^{r-1}_I u' \|_{L^2(I)}\\ &\lesssim \frac{{\it{\Delta}} t^{\mu-2}}{r^{s-2}} \| u \|_{H^s(I)}+ \frac{r^2}{{\it{\Delta}} t} \frac{{\it{\Delta}} t^{\mu-1}}{r^{s-1}} \| u \|_{H^s(I)}+ \frac{r^2}{{\it{\Delta}} t} \frac{{\it{\Delta}} t^{\mu-1}}{r^{s-1}} \| u \|_{H^s(I)}\\ &\lesssim \frac{{\it{\Delta}} t^{\mu-2}}{r^{s-3}} \| u \|_{H^s(I)}, \end{align*} where, in the second step, we have also used the inverse inequality $$\| \phi \|_{H^1(I)} \lesssim {\it{\Delta}} t^{-1}r^{2} \| \phi \|_{L^2(I)}$$ valid for any $$\phi\in \mathbb{P}^r(I)$$ (cf. Schwab, 1998), e.g., □ before proving the main result of this section, we notice that if $$\mathbf{v} \in {\bf{H}}^s(I)$$, for $$s \geq 2$$, then its projection $${\it{\Pi}}^{r}_I {\mathbf{v}} \in [\mathbb{P}^r(I)]^d$$ is defined, componentwise, according to (3.16) and the same kind of estimates given in (3.24), (3.23) and (3.26) hold. Now, let $$I = (0,T)$$, and let $$\mathbf{u} \in {\bf{H}}^s(0,T]$$, $$s \geq 2$$, be the solution of (2.3) and let $${\it{\Pi}}^\mathbf{r}_{I} \mathbf{u} \in {\mathcal{V}^\bf{r}}$$ be its projection such that $${\it{\Pi}}^\mathbf{r}_{I} \mathbf{u}_{|_{I_n}} = {\it{\Pi}}^{r_n}_{I_n} \mathbf{u} \in {{\bf{V}}_n^{r_n}}$$ is defined according to (3.16), for $$n=1,...,N$$. Then, it is easy to see that for $$n = 1,...,N$$ it holds \begin{align}(\mathbf{u} - {\it{\Pi}}^\textbf{r}_{I})(t_n) = 0, & \\ \end{align} (3.27) \begin{align}(\mathbf{u} - {\it{\Pi}}^\textbf{r}_{I})'(t_n) = 0, & \\ \end{align} (3.28) \begin{align}((\mathbf{u} - {\it{\Pi}}^\textbf{r}_{I})', q)_{I_n} = 0 & \qquad \forall \; q \in \mathbb{P}^{r_n-2}(I_n).\end{align} (3.29) Finally, we observe that (2.7) is strongly consistent (cf. Quarteroni, 2014). Indeed, it holds \begin{align}\mathcal{A}( \mathbf{u} - \mathbf{u}_{\textrm{DG}}, \mathbf{v} ) = 0 \qquad \forall \; \mathbf{v} \in {\mathcal{V}^\bf{r}}.\end{align} (3.30) We are now ready to prove the following convergence result. Theorem 3.12 Let $$\mathbf{u}$$ be the solution of (2.1)–(2.2) and let $${\mathbf{u}_{\textrm{DG}}} \in {\mathcal{V}^\bf{r}}$$ be its DG finite element approximation based on employing formulation (2.7). If $$\mathbf{u}|_{I_n} \in {\bf{H}}^{s_n}(I_n)$$, for any $$n=1,...,N$$, with $$s_n \geq 2$$, then \begin{align}\left\||{\mathbf{u} - {\mathbf{u}_{\textrm{DG}}}}\right\|| \lesssim \sum_{n=1}^N \frac{{\it{\Delta}} t_n^{{\mu_n} - \frac{3}{2}}}{ r_n^{s_n-3}} \|\mathbf{u}\|_{{\bf{H}}^{s_n}(I_n)},\end{align} (3.31) where $$\mu_n = \min(r_n+1,s_n)$$ for any $$n=1,...,N$$, and where the hidden constant depends on the norm of the matrices $$L$$ and $$K$$. Proof We define the error on the interval $$I=(0,T)$$ as $$\mathbf{e} = \mathbf{u} - \mathbf{u}_{\textrm{DG}}$$, and we split it as $$\mathbf{e} = {\mathbf{e}^{\pi}} + {\mathbf{e}^h}$$, where $${\mathbf{e}^{\pi}}$$ is the projection error while $${\mathbf{e}^h} \in {\mathcal{V}^\bf{r}}$$ is the remainder, i.e., $${\mathbf{e}^{\pi}} = \mathbf{u} - {\it{\Pi}}^{\textbf{r}}_I \textbf{u}$$ and $${\mathbf{e}^h} = {\it{\Pi}}^{\textbf{r}}_I \textbf{u} - \mathbf{u}_{\textrm{DG}}$$. Clearly, we have $$\left\||{\mathbf{e}}\right\|| \leq \left\||{\textbf{e}^\pi}\right\|| + \left\||{\textbf{e}^h}\right\||.$$ By exploiting properties (3.27)–(3.29) and by employing estimates (3.23) and (3.26), we can bound the term $$\left\||{\textbf{e}^\pi}\right\||$$ as \begin{align*}\left\||{\textbf{e}^\pi}\right\||^2 & = \sum_{n=1}^N \|L^\frac{1}{2} {\dot{\mathbf{e}}^{\pi}} \|^2_{L^2(I_n)} + \frac{1}{2}\sum_{n=1}^N {\dot{\mathbf{e}}^{\pi}}(t_{n-1}^+)^2 = \sum_{n=1}^N \|L^\frac{1}{2} {\dot{\mathbf{e}}^{\pi}} \|^2_{L^2(I_n)} + \frac{1}{2}\sum_{n=1}^N \left( - \int_{t_{n-1}}^{t_n} {\ddot{\mathbf{e}}^{\pi}}(s) \,{\,\rm d} s \right)^2 \nonumber\\ & \lesssim \sum_{n=1}^N \left( \|{\dot{\mathbf{e}}^{\pi}} \|^2_{L^2(I_n)} + {{\it{\Delta}} t_n} \|{\ddot{\mathbf{e}}^{\pi}} \|^2_{L^2(I_n)}\right) \nonumber \lesssim \sum_{n=1}^N \frac{{\it{\Delta}} t_n^{2\mu_n-3}}{r^{2s_n-{6}}} \| {\bf{u}} \|_{H^{s_n}(I_n)}^2, \end{align*} with $$\mu_n = \min(r_n+1,s_n)$$ for any $$n=1,...,N$$. For the term $$\left\||{\textbf{e}^h}\right\||$$, the Galerkin orthogonality the Galerkin orthogonality (3.30) yields \begin{align*}\left\||{\textbf{e}^h}\right\||^2 & = {\mathcal A}( \textbf{e}^h, \textbf{e}^h) = - {\mathcal A}( \textbf{e}^\pi, \textbf{e}^h)\\ & = - \sum_{n=1}^N \Big( ( {\ddot{\mathbf{e}}^{\pi}} ,{\dot{\mathbf{e}}^h})_{I_n} + (L {\dot{\mathbf{e}}^{\pi}}, {\dot{\mathbf{e}}^h})_{I_n} +(K{\mathbf{e}^{\pi}}, {\dot{\mathbf{e}}^h} )_{I_n} + [ {\dot{\mathbf{e}}^{\pi}} ]_n \cdot {\dot{\mathbf{e}}^h} (t_n^+) + K[{\mathbf{e}^{\pi}}]_n \cdot {\mathbf{e}^h}(t_n^+) \Big)\\ & \quad- {\dot{\mathbf{e}}^{\pi}} (t_0^+) \cdot {\dot{\mathbf{e}}^h} (t_0^+) - K{\mathbf{e}^{\pi}}(t_0^+) \cdot {\mathbf{e}^h}(t_0^+). \end{align*} Integrating by parts the term $$( {\ddot{\mathbf{e}}^{\pi}} ,{\dot{\mathbf{e}}^h})_{I_n}$$ and rearranging the addends, we obtain \begin{align*}\left\||{\textbf{e}^h}\right\||^2 & = \sum_{n=1}^N \Big( ( {\dot{\mathbf{e}}^{\pi}} ,{\ddot{\mathbf{e}}^h})_{I_n} - (L {\dot{\mathbf{e}}^{\pi}}, {\dot{\mathbf{e}}^h})_{I_n} -(K{\mathbf{e}^{\pi}}, {\dot{\mathbf{e}}^h} )_{I_n} \Big) \\ & \qquad + \sum_{n=1}^{N-1} \Big( [ {\dot{\mathbf{e}}^h} ]_n \cdot {\dot{\mathbf{e}}^{\pi}} (t_n^-) - K [{\mathbf{e}^{\pi}}]_n \cdot {\mathbf{e}^h}(t_n^+) \Big) - {\dot{\mathbf{e}}^{\pi}} (T^-) \cdot {\dot{\mathbf{e}}^h} (T^-) - K{\mathbf{e}^{\pi}}(0^+) \cdot {\mathbf{e}^h}(0^+). \end{align*} Using now properties (3.27)–(3.29) into the above equation yields \begin{align*}\left\||{\textbf{e}^h}\right\||^2 & = \sum_{n=1}^N \Big( - (L {\dot{\mathbf{e}}^{\pi}}, {\dot{\mathbf{e}}^h})_{I_n} - (K{\mathbf{e}^{\pi}}, {\dot{\mathbf{e}}^h} )_{I_n} \Big). \end{align*} Applying Cauchy–Schwarz and arithmetic–geometric inequalities and using estimates (3.24) and (3.23), it is easy to obtain \begin{align*}\left\||{\textbf{e}^h}\right\||^2 & \lesssim \sum_{n=1}^N \Big( \|{\dot{\mathbf{e}}^{\pi}}\|_{L^2(In)}^2 + \|{\mathbf{e}^{\pi}}\|_{L^2(In)}^2 \Big) \lesssim { \sum_{n=1}^N \frac{{\it{\Delta}} t_n^{2\mu_n-2}}{{r_n}^{2s_n-{2}}}\| {\bf{u}} \|_{H^{s_n}(I_n)}^2, }\end{align*} where $$\mu_n = \min(r_n+1,s_n)$$ for any $$n=1,...,N$$, and where the hidden constant depends on the norms $$\|L^\frac{1}{2}\|_{\infty}$$ and $$\|L^{-\frac{1}{2}}K \|_{\infty}$$. The thesis follows by combining the above bound with the estimate on $$\left\||{\textbf{e}^\pi}\right\||$$. □ Remark 3.13 In the case when $${\it{\Delta}} t_n = {\it{\Delta}} t$$, $$r_n = r\geq 2$$ and $$s_n = s \geq 2$$, for $$n=1,...,N$$, the error bound of Theorem 3.12 becomes \begin{align*}\left\||{\mathbf{u} - {\mathbf{u}_{\textrm{DG}}}}\right\|| \lesssim \frac{{\it{\Delta}} t^{r-\frac{1}{2}}}{ r^{s-{3}} }\|\mathbf{u}\|_{{\bf{H}}^{s}(0,T)}, \end{align*} where the hidden constant depends on the norm of the matrices $$L$$ and $$K$$. 4. Application to the elastodynamics problem In this section, we apply the method presented in Section 2 to the simulation of elastic wave propagations in heterogeneous media. In particular, we will adopt the DG method previously developed to handle the time integration of the second-order differential system arising after space discretization obtained with the DGSE approach proposed in Antonietti et al. (2012). Since the focus of this article is on time integration, here, we simply report the mathematical model of linear viscoelasticity and the algebraic linear system resulting after the space discretization obtained employing the DGSE scheme of Antonietti et al. (2012). Finally, some numerical results are discussed. 4.1 Mathematical modelling of seismic wave propagations and its algebraic formulation For a given open-bounded domain $${\it{\Omega}} \subset \mathbb{R}^m$$, $$m=2,3$$, we consider the following problem: for $$T > 0$$ find $$\mathbf{u}: {\it{\Omega}} \times [0,T] \rightarrow \mathbb{R}^m$$ such that \begin{align}\displaystyle \rho \partial_{tt} {\bf{u}} + 2\rho \zeta \partial_{t}{\mathbf{u}} + \rho \zeta^2 \mathbf{u} -{\nabla \cdot{{\bf{\sigma}}}} & = {\bf{f}}, && \textrm{in} \; {\it{\Omega}} \times (0,T], \\ \end{align} (4.1a) \begin{align}\displaystyle {\bf{u}} & = {\bf{0}}, && \textrm{on} \; {\it{\Gamma}}_D\times (0,T], \\ \end{align} (4.1b) \begin{align}\displaystyle {\bf{\sigma}} {\bf{n}} & = {\bf{0}}, && \textrm{on} \; {\it{\Gamma}}_N\times (0,T], \\ \end{align} (4.1c) \begin{align}\displaystyle \partial_t {\bf{u}} & = {\bf{u}}_1, && \textrm{in} \; {\it{\Omega}} \times \{ 0\}, \\ \end{align} (4.1d) \begin{align}\displaystyle {\bf{u}} &= {\bf{u}}_0, &&\textrm{in} \; {\it{\Omega}} \times \{ 0\},\end{align} (4.1e) where $$\overline{\partial {\it{\Omega}}} = \overline{{\it{\Gamma}}_D} \cup \overline{{\it{\Gamma}}_N}$$ with $${\it{\Gamma}}_D \cap {\it{\Gamma}}_N =\emptyset$$, $$\mathbf{f} \in L^2 ((0,T]; {\bf{L}}^2({\it{\Omega}}))$$ is the source term, and $$\rho \in L^\infty({\it{\Omega}})$$ is such that $$\rho = \rho(\mathbf{x}) > 0$$ for almost any $$\mathbf{x} \in {\it{\Omega}}$$. The function $$\zeta = \zeta ({\bf{x}}) > 0$$ is a decay factor whose dimension is the inverse of time and is a piecewise constant function, then $$\zeta \in L^{\infty}({\it{\Omega}})$$. We suppose the stress tensor $${\bf{\sigma}}$$ to be related to the strain tensor $${\bf{\varepsilon}}(\mathbf{u}) = \frac{1}{2}\bigl( \nabla\mathbf{u} + \nabla\mathbf{u}^T \bigr)$$ through the Hooke’s law, i.e., \begin{align*}{\bf{\sigma}} = {2\mu}{\bf{\varepsilon}} + \lambda {\text{tr}}({\bf{\varepsilon}})\mathbb{I}, \end{align*} where $$\lambda = \lambda({\bf{x}})$$ and $$\mu=\mu({\bf{x}})$$ are the Lamé elastic coefficients of the medium, $${\text{tr}}(\cdot)$$ is the trace operator and $$\mathbb{I}\in \mathbb{R}^{m \times m}$$ is the identity tensor. Here and in the following, we will suppose $$\rho, \lambda$$ and $$\mu$$ to be uniformly bounded functions in $${\it{\Omega}}$$, i.e., $$\rho, \lambda, \mu \in L^{\infty}({\it{\Omega}})$$. The semidiscretization of problem (4.1a)–(4.1e) by the DGSE technique (Antonietti et al., 2012) results into the following second-order differential system for the nodal displacement $${\bf{U}}(t) \,{=} [{\bf{U}}^1(t),..., {\bf{U}}^m(t) ]^T$$ \begin{align}\begin{cases}M \ddot{{\bf{U}}}(t) + C \dot{{\bf{U}}}(t) + D{\bf{U}}(t) + A {\bf{U}}(t) = {\bf{F}}(t), & t \in (0,T], \\ \dot{\textbf{U}}(0) = {\bf{u}}_1, & \\ \textbf{U}(0) = {\bf{u}}_0, & \\ \end{cases}\end{align} (4.2) where $$\ddot{{\bf{U}}}(t)$$ (respectively, $$\dot{{\bf{U}}}(t)$$) represents the vector of nodal acceleration (respectively velocity) and $${\bf{F}}(t)$$ the vector of externally applied loads, i.e., $${\bf{F}}(t)= [{\bf{F}}^1(t),..., {\bf{F}}^m(t)]^T$$. Here, $$M$$ is the diagonal mass matrix, $$C$$ and $$D$$ stand for the structural damping (and have the same structure of $$M$$) while the action of the stiffness matrix $$A$$ to the displacement vector $${\bf{U}}$$ represents the internal (elastic) forces. We remark that the matrix $$A$$ in (4.2) is symmetric and positive definite (see Antonietti et al., 2012). The latter properties are also verified by the matrices $$M, C$$ and $$D$$. To rewrite (4.2) in the form (2.1), we multiply it by $$M^{-\frac{1}{2}}$$ and set $$\mathbf{Z}(t) = M^{\frac{1}{2}} \mathbf{U}(t)$$ and obtain \begin{align}\begin{cases}\ddot{\mathbf{Z}}(t) + L \dot{\mathbf{Z}}(t) + K \mathbf{Z}(t) = \mathbf{G}(t), & t \in (0,T], \\ \dot{\textbf{Z}}(0) = M^{\frac{1}{2}} {\bf{u}}_1, & \\ {\textbf{Z}}(0) = M^{\frac{1}{2}} {\bf{u}}_0, & \\ \end{cases}\end{align} (4.3) where $$L = M^{-\frac{1}{2}}C M^{-\frac{1}{2}}, \, K = M^{-\frac{1}{2}}(D + A) M^{-\frac{1}{2}}$$, are symmetric and positive definite and $$\mathbf{G}(t) = M^{-\frac{1}{2}}\mathbf{F}(t)$$. 5 Numerical results In this section, we present some numerical results to highlight the behavior of the DG method introduced in Section 2. In particular, we first focus our attention on a scalar benchmark, then we test it on the system of differential equations arisen by the space discretization of the elastodynamics problem. 5.1 Numerical results for a scalar problem We consider the interval $$I=(0,T]$$, $$T=10$$, subdivided into $$N$$ time slabs $$I_n$$, for $$n=1,..., N$$, of uniform length $${\it{\Delta}} t$$. Moreover, we suppose the polynomial approximation vector $$\textbf{r}$$ to be constant for each time slab, i.e., $$r_1=...=r_N=r\geq 2$$. We consider the following problem \begin{align}\begin{cases}\ddot{u}(t) + 5\dot{u}(t) + 6u(t) = 0 \qquad \forall t \in (0, 10], \\ u(0) = 2, \\ \dot{u}(0) = -5, \end{cases}\end{align} (5.1) whose exact solution is $$u(t) = e^{-3t} + e^{-2t}, \, t \in [0, 10].$$ We compute the error $$\left\||{u-u_{\textrm{DG}}}\right\||$$ versus $$1/{\it{\Delta}} t$$, with $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,2,3$$, for several polynomial degrees $$r = 2,3,4,5$$. The computed errors for the scalar test case are reported in Fig. 2, in log-log scale. As expected, a convergence rate of $${\it{\Delta}} t^{r - \frac{1}{2}}$$ is observed (cf. Theorem 3.12). This is also confirmed by the results shown in Table 1, where the computed convergence rates are reported. Although we do not have a theoretical proof, for completeness, we also report in Table 2 the computed error with respect to the $$\| \cdot \|_{L^2(0,T)}$$ and $$\| \cdot \|_{H^1(0,T)}$$ norms and the corresponding rates of convergence. The method achieves an optimal rate of convergence with respect to both norms, i.e., the error measured in the $$L^2$$-norm (respectively $$H^1$$-norm) decays as $${\it{\Delta}} t^r$$ (respectively $${\it{\Delta}} t^{r-1}$$) for $${\it{\Delta}} t$$ going to $$0$$. In Fig. 3, we display the computed errors versus the polynomial degree for different values of $${\it{\Delta}} t$$. As expected, an exponential convergence rate is observed. Fig. 2. View largeDownload slide Scalar test case. Computed errors $$\left\||{ u - u_{\textrm{DG}}}\right\||$$ versus $$1/{\it{\Delta}} t$$ for $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,2,3$$ (loglog scale) and $$r = 2,3,4,5$$. Fig. 2. View largeDownload slide Scalar test case. Computed errors $$\left\||{ u - u_{\textrm{DG}}}\right\||$$ versus $$1/{\it{\Delta}} t$$ for $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,2,3$$ (loglog scale) and $$r = 2,3,4,5$$. Fig. 3. View largeDownload slide Scalar test case. Computed errors $$\left\||{ u - u_{\textrm{DG}}}\right\||$$ versus $$r = 2,3,4,5$$ (semilogy scale) for $${\it{\Delta}} t = 0.025, 0.0125, 0.00625$$. Fig. 3. View largeDownload slide Scalar test case. Computed errors $$\left\||{ u - u_{\textrm{DG}}}\right\||$$ versus $$r = 2,3,4,5$$ (semilogy scale) for $${\it{\Delta}} t = 0.025, 0.0125, 0.00625$$. Table 1 Scalar test case. Computed convergence rates in the $$\left\||{\cdot}\right\||$$ norm with respect to the polynomial degree $$r=2,3,4,5$$ $${\it{\Delta}} t$$ $$r = 2$$ $$r = 3$$ $$r = 4$$ $$r = 5$$ 1.00e$$-$$0 — — — — 5.00e$$-$$1 1.4017 2.2646 3.1802 4.1246 2.50e$$-$$1 1.5224 2.4577 3.4305 4.4403 1.25e$$-$$1 1.5295 2.5101 3.6584 4.5516 $${\it{\Delta}} t$$ $$r = 2$$ $$r = 3$$ $$r = 4$$ $$r = 5$$ 1.00e$$-$$0 — — — — 5.00e$$-$$1 1.4017 2.2646 3.1802 4.1246 2.50e$$-$$1 1.5224 2.4577 3.4305 4.4403 1.25e$$-$$1 1.5295 2.5101 3.6584 4.5516 Table 1 Scalar test case. Computed convergence rates in the $$\left\||{\cdot}\right\||$$ norm with respect to the polynomial degree $$r=2,3,4,5$$ $${\it{\Delta}} t$$ $$r = 2$$ $$r = 3$$ $$r = 4$$ $$r = 5$$ 1.00e$$-$$0 — — — — 5.00e$$-$$1 1.4017 2.2646 3.1802 4.1246 2.50e$$-$$1 1.5224 2.4577 3.4305 4.4403 1.25e$$-$$1 1.5295 2.5101 3.6584 4.5516 $${\it{\Delta}} t$$ $$r = 2$$ $$r = 3$$ $$r = 4$$ $$r = 5$$ 1.00e$$-$$0 — — — — 5.00e$$-$$1 1.4017 2.2646 3.1802 4.1246 2.50e$$-$$1 1.5224 2.4577 3.4305 4.4403 1.25e$$-$$1 1.5295 2.5101 3.6584 4.5516 Table 2 Scalar test case. Computed error in the $$\| \cdot \|_{L^2(0,T)}$$ and $$\| \cdot \|_{H^1(0,T)}$$ norms and corresponding convergence rates with respect to the polynomial approximation degree $$r=2,3,4,5$$ $$r$$ $${\it{\Delta}} t$$ $$\| \cdot \|_{L^2(0,T)}$$ rate $$\| \cdot \|_{H^1(0,T)}$$ rate 1.00e$$-$$0 4.4902e$$-$$02 — 6.7961e$$-$$01 — 8.00e$$-$$1 2.9612e$$-$$02 1.8655 4.9201e$$-$$01 1.4476 2 4.00e$$-$$1 7.0288e$$-$$03 2.0749 1.5416e$$-$$01 1.6743 2.00.e$$-$$1 1.2044e$$-$$03 2.5450 4.2254e$$-$$02 1.8672 1.00.e$$-$$1 1.7331e$$-$$04 2.7968 1.0988e$$-$$02 1.9431 1.00e$$-$$0 6.9895e$$-$$03 — 1.4649e$$-$$01 — 8.00e$$-$$1 3.9782e$$-$$03 2.5257 8.6484e$$-$$02 2.3617 3 4.00e$$-$$1 4.5984e$$-$$04 3.1129 1.3821e$$-$$02 2.6455 2.00.e$$-$$1 3.6876e$$-$$05 3.6404 1.8863e$$-$$03 2.8732 1.00.e$$-$$1 2.5546e$$-$$06 3.8515 2.4343e$$-$$04 2.9540 1.00e$$-$$0 1.0209e$$-$$03 — 2.4885e$$-$$02 — 8.00e$$-$$1 4.4989e$$-$$04 3.6723 1.1885e$$-$$02 3.3119 4 4.00e$$-$$1 2.3662e$$-$$05 4.2489 9.6064e$$-$$04 3.6290 2.00.e$$-$$1 9.0456e$$-$$07 4.7092 6.5449e$$-$$05 3.8755 1.00.e$$-$$1 3.0658e$$-$$08 4.8829 4.2099e$$-$$06 3.9585 1.00e$$-$$0 1.2456e$$-$$04 — 3.4801e$$-$$03 — 8.00e$$-$$1 4.2903e$$-$$05 4.7766 1.3399e$$-$$03 4.2773 5 4.00e$$-$$1 1.0725e$$-$$06 5.3220 5.4621e$$-$$05 4.6165 2.00.e$$-$$1 1.9994e$$-$$08 5.7452 1.8603e$$-$$06 4.8759 1.00.e$$-$$1 3.3483e$$-$$10 5.9000 5.9740e$$-$$08 4.9607 $$r$$ $${\it{\Delta}} t$$ $$\| \cdot \|_{L^2(0,T)}$$ rate $$\| \cdot \|_{H^1(0,T)}$$ rate 1.00e$$-$$0 4.4902e$$-$$02 — 6.7961e$$-$$01 — 8.00e$$-$$1 2.9612e$$-$$02 1.8655 4.9201e$$-$$01 1.4476 2 4.00e$$-$$1 7.0288e$$-$$03 2.0749 1.5416e$$-$$01 1.6743 2.00.e$$-$$1 1.2044e$$-$$03 2.5450 4.2254e$$-$$02 1.8672 1.00.e$$-$$1 1.7331e$$-$$04 2.7968 1.0988e$$-$$02 1.9431 1.00e$$-$$0 6.9895e$$-$$03 — 1.4649e$$-$$01 — 8.00e$$-$$1 3.9782e$$-$$03 2.5257 8.6484e$$-$$02 2.3617 3 4.00e$$-$$1 4.5984e$$-$$04 3.1129 1.3821e$$-$$02 2.6455 2.00.e$$-$$1 3.6876e$$-$$05 3.6404 1.8863e$$-$$03 2.8732 1.00.e$$-$$1 2.5546e$$-$$06 3.8515 2.4343e$$-$$04 2.9540 1.00e$$-$$0 1.0209e$$-$$03 — 2.4885e$$-$$02 — 8.00e$$-$$1 4.4989e$$-$$04 3.6723 1.1885e$$-$$02 3.3119 4 4.00e$$-$$1 2.3662e$$-$$05 4.2489 9.6064e$$-$$04 3.6290 2.00.e$$-$$1 9.0456e$$-$$07 4.7092 6.5449e$$-$$05 3.8755 1.00.e$$-$$1 3.0658e$$-$$08 4.8829 4.2099e$$-$$06 3.9585 1.00e$$-$$0 1.2456e$$-$$04 — 3.4801e$$-$$03 — 8.00e$$-$$1 4.2903e$$-$$05 4.7766 1.3399e$$-$$03 4.2773 5 4.00e$$-$$1 1.0725e$$-$$06 5.3220 5.4621e$$-$$05 4.6165 2.00.e$$-$$1 1.9994e$$-$$08 5.7452 1.8603e$$-$$06 4.8759 1.00.e$$-$$1 3.3483e$$-$$10 5.9000 5.9740e$$-$$08 4.9607 Table 2 Scalar test case. Computed error in the $$\| \cdot \|_{L^2(0,T)}$$ and $$\| \cdot \|_{H^1(0,T)}$$ norms and corresponding convergence rates with respect to the polynomial approximation degree $$r=2,3,4,5$$ $$r$$ $${\it{\Delta}} t$$ $$\| \cdot \|_{L^2(0,T)}$$ rate $$\| \cdot \|_{H^1(0,T)}$$ rate 1.00e$$-$$0 4.4902e$$-$$02 — 6.7961e$$-$$01 — 8.00e$$-$$1 2.9612e$$-$$02 1.8655 4.9201e$$-$$01 1.4476 2 4.00e$$-$$1 7.0288e$$-$$03 2.0749 1.5416e$$-$$01 1.6743 2.00.e$$-$$1 1.2044e$$-$$03 2.5450 4.2254e$$-$$02 1.8672 1.00.e$$-$$1 1.7331e$$-$$04 2.7968 1.0988e$$-$$02 1.9431 1.00e$$-$$0 6.9895e$$-$$03 — 1.4649e$$-$$01 — 8.00e$$-$$1 3.9782e$$-$$03 2.5257 8.6484e$$-$$02 2.3617 3 4.00e$$-$$1 4.5984e$$-$$04 3.1129 1.3821e$$-$$02 2.6455 2.00.e$$-$$1 3.6876e$$-$$05 3.6404 1.8863e$$-$$03 2.8732 1.00.e$$-$$1 2.5546e$$-$$06 3.8515 2.4343e$$-$$04 2.9540 1.00e$$-$$0 1.0209e$$-$$03 — 2.4885e$$-$$02 — 8.00e$$-$$1 4.4989e$$-$$04 3.6723 1.1885e$$-$$02 3.3119 4 4.00e$$-$$1 2.3662e$$-$$05 4.2489 9.6064e$$-$$04 3.6290 2.00.e$$-$$1 9.0456e$$-$$07 4.7092 6.5449e$$-$$05 3.8755 1.00.e$$-$$1 3.0658e$$-$$08 4.8829 4.2099e$$-$$06 3.9585 1.00e$$-$$0 1.2456e$$-$$04 — 3.4801e$$-$$03 — 8.00e$$-$$1 4.2903e$$-$$05 4.7766 1.3399e$$-$$03 4.2773 5 4.00e$$-$$1 1.0725e$$-$$06 5.3220 5.4621e$$-$$05 4.6165 2.00.e$$-$$1 1.9994e$$-$$08 5.7452 1.8603e$$-$$06 4.8759 1.00.e$$-$$1 3.3483e$$-$$10 5.9000 5.9740e$$-$$08 4.9607 $$r$$ $${\it{\Delta}} t$$ $$\| \cdot \|_{L^2(0,T)}$$ rate $$\| \cdot \|_{H^1(0,T)}$$ rate 1.00e$$-$$0 4.4902e$$-$$02 — 6.7961e$$-$$01 — 8.00e$$-$$1 2.9612e$$-$$02 1.8655 4.9201e$$-$$01 1.4476 2 4.00e$$-$$1 7.0288e$$-$$03 2.0749 1.5416e$$-$$01 1.6743 2.00.e$$-$$1 1.2044e$$-$$03 2.5450 4.2254e$$-$$02 1.8672 1.00.e$$-$$1 1.7331e$$-$$04 2.7968 1.0988e$$-$$02 1.9431 1.00e$$-$$0 6.9895e$$-$$03 — 1.4649e$$-$$01 — 8.00e$$-$$1 3.9782e$$-$$03 2.5257 8.6484e$$-$$02 2.3617 3 4.00e$$-$$1 4.5984e$$-$$04 3.1129 1.3821e$$-$$02 2.6455 2.00.e$$-$$1 3.6876e$$-$$05 3.6404 1.8863e$$-$$03 2.8732 1.00.e$$-$$1 2.5546e$$-$$06 3.8515 2.4343e$$-$$04 2.9540 1.00e$$-$$0 1.0209e$$-$$03 — 2.4885e$$-$$02 — 8.00e$$-$$1 4.4989e$$-$$04 3.6723 1.1885e$$-$$02 3.3119 4 4.00e$$-$$1 2.3662e$$-$$05 4.2489 9.6064e$$-$$04 3.6290 2.00.e$$-$$1 9.0456e$$-$$07 4.7092 6.5449e$$-$$05 3.8755 1.00.e$$-$$1 3.0658e$$-$$08 4.8829 4.2099e$$-$$06 3.9585 1.00e$$-$$0 1.2456e$$-$$04 — 3.4801e$$-$$03 — 8.00e$$-$$1 4.2903e$$-$$05 4.7766 1.3399e$$-$$03 4.2773 5 4.00e$$-$$1 1.0725e$$-$$06 5.3220 5.4621e$$-$$05 4.6165 2.00.e$$-$$1 1.9994e$$-$$08 5.7452 1.8603e$$-$$06 4.8759 1.00.e$$-$$1 3.3483e$$-$$10 5.9000 5.9740e$$-$$08 4.9607 5.2 Numerical results for the elastodynamics problem The numerical results presented in the sequel have been obtained by using the open source software SPEED (http://speed.mox.polimi.it; accessed 11 October 2017) suitably adapted to apply the DG scheme presented in Section 2 to the system (4.3). For all the numerical simulations, we consider the interval $$I=(0,T]$$ subdivided into $$N$$ time slabs $$I_n$$, for $$n=1,..., N$$ having uniform length $${\it{\Delta}} t$$. 5.2.1 Test case 1 We consider $${\it{\Omega}} = (0,1)^2, {\it{\Gamma}}_D = \partial {\it{\Omega}}$$ and $$T = 50$$. We set the mass density $$\rho=1$$, the Lamé coefficients $$\lambda=\mu=1$$, $$\zeta=0.01$$ and choose the data $$\mathbf{f}, \mathbf{u}_0, \mathbf{u}_1$$ such that the exact solution of problem in (4.1a)–(4.1e) is given by \begin{align*}\mathbf{u}= \sin(\sqrt{2}\pi t) \begin{bmatrix}- \sin^2(\pi x) \sin(2\pi y) \\ \sin(2\pi x) \sin^2(\pi y) \end{bmatrix}\!. \end{align*} In the first example, we compute the errors $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ versus $$1/ {\it{\Delta}} t$$ for $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,2,3,4$$ and varying the polynomial degree $$r = 2,3,4$$. As of the space discretization of the domain $${\it{\Omega}}$$, we consider a Cartesian grid with characteristic size $$h = 0.125$$ and we set a polynomial approximation degree $$q=r+1$$. The computed errors are reported in Fig. 4. Fig. 4. View largeDownload slide Test case 1. Computed errors $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ as a function of $$1/ {\it{\Delta}} t$$ for $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,\dots, 4$$ (log-log scale), $$r=2,3,4$$, $$h=0.125$$ and $$q=r+1$$. Fig. 4. View largeDownload slide Test case 1. Computed errors $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ as a function of $$1/ {\it{\Delta}} t$$ for $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,\dots, 4$$ (log-log scale), $$r=2,3,4$$, $$h=0.125$$ and $$q=r+1$$. A convergence rate of $${\it{\Delta}} t^{r - \frac{1}{2}}$$ is observed, in accordance with the theoretical estimate (3.31). In Table 3, the computed convergence rates are reported for $$r = 2,3,4$$, respectively. Table 3 Test case 1. Computed errors $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ and computed convergence rates, $$r = 2,3,4$$, $$h=0.125$$ and $$q=r+1$$ $$r$$ $${\it{\Delta}} t$$ $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ rate 1.00e$$-$$0 4.1632e+02 — 5.00e$$-$$1 2.5100e+02 0.7351 2 2.50e$$-$$1 9.1706e+01 1.4474 1.25e$$-$$1 3.3960e+01 1.4331 6.25e$$-$$2 1.2156e+01 1.4821 1.00e$$-$$0 3.3909e+02 — 5.00e$$-$$1 9.2986e+01 1.8665 3 2.50e$$-$$1 1.4534e+01 2.6775 1.25e$$-$$1 2.5484e+00 2.5117 6.25e$$-$$2 4.4983e$$-$$01 2.5021 1.00e$$-$$0 1.7465e+02 — 5.00e$$-$$1 1.4155e+01 3.6250 4 2.50e$$-$$1 1.3873e+00 3.3509 1.25e$$-$$1 1.2308e$$-$$01 3.4946 6.25e$$-$$2 1.0550e$$-$$02 3.5441 $$r$$ $${\it{\Delta}} t$$ $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ rate 1.00e$$-$$0 4.1632e+02 — 5.00e$$-$$1 2.5100e+02 0.7351 2 2.50e$$-$$1 9.1706e+01 1.4474 1.25e$$-$$1 3.3960e+01 1.4331 6.25e$$-$$2 1.2156e+01 1.4821 1.00e$$-$$0 3.3909e+02 — 5.00e$$-$$1 9.2986e+01 1.8665 3 2.50e$$-$$1 1.4534e+01 2.6775 1.25e$$-$$1 2.5484e+00 2.5117 6.25e$$-$$2 4.4983e$$-$$01 2.5021 1.00e$$-$$0 1.7465e+02 — 5.00e$$-$$1 1.4155e+01 3.6250 4 2.50e$$-$$1 1.3873e+00 3.3509 1.25e$$-$$1 1.2308e$$-$$01 3.4946 6.25e$$-$$2 1.0550e$$-$$02 3.5441 Table 3 Test case 1. Computed errors $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ and computed convergence rates, $$r = 2,3,4$$, $$h=0.125$$ and $$q=r+1$$ $$r$$ $${\it{\Delta}} t$$ $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ rate 1.00e$$-$$0 4.1632e+02 — 5.00e$$-$$1 2.5100e+02 0.7351 2 2.50e$$-$$1 9.1706e+01 1.4474 1.25e$$-$$1 3.3960e+01 1.4331 6.25e$$-$$2 1.2156e+01 1.4821 1.00e$$-$$0 3.3909e+02 — 5.00e$$-$$1 9.2986e+01 1.8665 3 2.50e$$-$$1 1.4534e+01 2.6775 1.25e$$-$$1 2.5484e+00 2.5117 6.25e$$-$$2 4.4983e$$-$$01 2.5021 1.00e$$-$$0 1.7465e+02 — 5.00e$$-$$1 1.4155e+01 3.6250 4 2.50e$$-$$1 1.3873e+00 3.3509 1.25e$$-$$1 1.2308e$$-$$01 3.4946 6.25e$$-$$2 1.0550e$$-$$02 3.5441 $$r$$ $${\it{\Delta}} t$$ $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ rate 1.00e$$-$$0 4.1632e+02 — 5.00e$$-$$1 2.5100e+02 0.7351 2 2.50e$$-$$1 9.1706e+01 1.4474 1.25e$$-$$1 3.3960e+01 1.4331 6.25e$$-$$2 1.2156e+01 1.4821 1.00e$$-$$0 3.3909e+02 — 5.00e$$-$$1 9.2986e+01 1.8665 3 2.50e$$-$$1 1.4534e+01 2.6775 1.25e$$-$$1 2.5484e+00 2.5117 6.25e$$-$$2 4.4983e$$-$$01 2.5021 1.00e$$-$$0 1.7465e+02 — 5.00e$$-$$1 1.4155e+01 3.6250 4 2.50e$$-$$1 1.3873e+00 3.3509 1.25e$$-$$1 1.2308e$$-$$01 3.4946 6.25e$$-$$2 1.0550e$$-$$02 3.5441 5.2.2 Test case 2: a numerical test with nonreflective boundary conditions When simulating seismic wave propagations, ideally artificial boundaries should consent incident wave to be propagated without producing any reflection. A simple strategy consists in imposing nonreflective conditions on artificial boundaries $${\it{\Gamma}}_{NR} \subset \partial{\it{\Omega}}$$ such that $${\it{\Gamma}}_{NR} \cap {\it{\Gamma}}_D \cap {\it{\Gamma}}_N = \emptyset$$, $$\begin{cases}\dfrac{\partial ({\bf{u}} \cdot {\bf{n}})}{\partial n} = -\dfrac{1}{V_P} \dfrac{\partial ({\bf{u}} \cdot {\bf{n}})}{\partial t} + \dfrac{V_S - V_P}{V_P} \dfrac{\partial ({\bf{u}} \cdot {\bf{\tau}})}{\partial t}, \\ \dfrac{\partial ({\bf{u}} \cdot {\bf{\tau}})}{\partial n} = -\dfrac{1}{V_S} \dfrac{\partial ({\bf{u}} \cdot {\bf{\tau}})}{\partial t} + \dfrac{V_S - V_P}{V_S} \dfrac{\partial ({\bf{u}} \cdot {\bf{n}})}{\partial t}, \end{cases}$$ (5.2) where $$\mathbf{n} = (n_1, n_2)$$ (respectively $${\bf{\tau}} = (\tau_1, \tau_2)$$) is the unit normal (respectively tangential) vector to $$\partial {\it{\Omega}}$$. Here, $$V_{\rm P} = \sqrt{(\lambda+2\mu)/\rho}$$ and $$V_{\rm S} = \sqrt{\mu/\rho}$$ are the propagation velocities of compressional (P) and shear (S) waves, respectively. Equations (5.2) are first-order nonreflecting boundary conditions (see, e.g., Stacey, 1988). Their use breaks of symmetry in the matrix $$K$$ in (4.3). The test consists of propagating a pure shear plane wave through a viscoelastic, horizontally layered soil profile (cf. Fig. 5). The mechanical properties of the layers are summarized in Table 4. The incidence is orthogonal to the free surface and the excitation consists of a displacement Ricker wavelet $$\textbf{f} (\textbf{x}, t) = h(t) \delta_q(\textbf{x} - \textbf{x}_0)$$ with \begin{equation*}h(t) = h_0 \left(1-2 \left(\pi f_{\rm peak}\right)^2 \left(t-t_0\right)^2 \right)e^{-(\pi f_{\rm peak})^2 (t-t_0)^ 2}, \end{equation*} where $$f_{\rm peak} = 1$$ Hz, $$t_0 = 2$$ s and $$h_0=1$$ is the amplitude of the wave in time domain. Here, $$\delta_q$$ is the numerical delta function, i.e., a polynomial of degree $$q$$ that approximates the $$\delta$$ distribution. In this case, the forcing term is applied at points $$\textbf{x}_0$$ lying at the bottom of the domain (see Fig. 5). Table 4 Mechanical properties Material $$\rho$$ $$\lambda$$ $$\mu$$ $$\zeta$$ 1 2000 2.00e+07 2.00e+07 3.1416e$$-$$02 2 2000 5.00e+08 5.00e+08 3.1416e$$-$$03 Material $$\rho$$ $$\lambda$$ $$\mu$$ $$\zeta$$ 1 2000 2.00e+07 2.00e+07 3.1416e$$-$$02 2 2000 5.00e+08 5.00e+08 3.1416e$$-$$03 Table 4 Mechanical properties Material $$\rho$$ $$\lambda$$ $$\mu$$ $$\zeta$$ 1 2000 2.00e+07 2.00e+07 3.1416e$$-$$02 2 2000 5.00e+08 5.00e+08 3.1416e$$-$$03 Material $$\rho$$ $$\lambda$$ $$\mu$$ $$\zeta$$ 1 2000 2.00e+07 2.00e+07 3.1416e$$-$$02 2 2000 5.00e+08 5.00e+08 3.1416e$$-$$03 Fig. 5. View largeDownload slide Domain $${\it{\Omega}}$$ for the test case 2. Nonreflective boundary conditions are imposed at the bottom edge, $${\bf{\sigma}}{\bf{n}} = \textbf{0}$$ is imposed at the top edge (free-surface condition), while homogeneous Dirichlet boundary conditions (in y-direction) are imposed on the lateral sides. Fig. 5. View largeDownload slide Domain $${\it{\Omega}}$$ for the test case 2. Nonreflective boundary conditions are imposed at the bottom edge, $${\bf{\sigma}}{\bf{n}} = \textbf{0}$$ is imposed at the top edge (free-surface condition), while homogeneous Dirichlet boundary conditions (in y-direction) are imposed on the lateral sides. The plane wave rises from the bottom of $${\it{\Omega}}$$, reaches the top of the computational domain, increases its amplitude (because of the free surface condition) and then is propagated backward, completely absorbed from the bottom boundary, thanks to the absorbing condition (5.2). To prevent spurious oscillations inside the domain, homogeneous Dirichlet boundary conditions (in y-direction) are imposed on the lateral sides of $${\it{\Omega}}$$ (cf. Fig. 5). For the simulation, we fix $$T=20 s$$, $${\it{\Delta}} t = 0.01 s, h=10$$ and $$q=4$$ for the space discretization, while we use second-order polynomials ($$r=2$$) for the time integration. We compare the results obtained using the proposed DG method with the analogous ones obtained coupling the DGSE space discretization with the classical leap-frog method (cf. Quarteroni & Valli, 2008, e.g., with $${\it{\Delta}} t = 0.0001$$ s. Notice that this is the biggest time step allowed by the CFL condition (Courant et al., 1928). We plot the displacement time histories (in the x-direction) recorded by two monitors: $$R_1$$ set on the free surface and $$R_2$$ located across Material 1 and Material 2 (cf. Fig. 5). In Fig. 6, we show the results obtained with the leap-frog and the DG method along with the semianalytical solution $${\bf{u}}_{TH}$$ based on the Thomson–Haskell propagation matrix method (see, e.g., Haskell, 1953). We can clearly see that all the methods produce almost identical solutions. To quantify the distance between the curves, in Fig. 6, we plot the error $$| {\bf{u}}_{TH}(t)- {\bf{u}}_{*}(t)|$$ for $$t\in(0,T]$$, where $${\bf{u}}_{*}$$ is either the solution obtained with the leap-frog scheme or the one obtained with the DG method and $${\bf{u}}_{TH}$$ is the Thomson–Haskell semianalytical solution. Both methods achieve the same level of accuracy. Finally, we compare the computational time of our method with the analogous results obtained with the leap-frog time integration scheme. In Table 5, we show the computational costs for each time step (second column) as well as the overall computational time for the time loop (third column). As expected, each time step is more expensive in our approach, as the solution of a linear system is involved, but overall this can be balanced using larger times steps, as no stability constraints are required. We point out that for our DG method, the matrix is assembled and factorized once and for all before the time loop. A more thorough comparison is under investigation, which also includes suitable parallelization techniques and possible preconditioned iterative methods. Fig. 6. View largeDownload slide Displacement registered by $$R_1$$ (a) and $$R_2$$ (b). Fig. 6. View largeDownload slide Displacement registered by $$R_1$$ (a) and $$R_2$$ (b). Table 5 Computational times for the time loop (space polynomial degree $$q = 4$$) Single time step Time loop $${\it{\Delta}} t$$ $$\| {\bf{u}}_{TH}- {\bf{u}}_{*}\|_{L^2(0,T)}$$ leap-frog 4.57e$$-$$3 s 9.15 s 1.e$$-$$4 s 1.70e$$-$$3 ($$R_1$$), 6.94e$$-$$4 ($$R_2$$) DG 2.98e$$-$$2 s 7.45 s 8.e$$-$$2 s 1.99e$$-$$3 ($$R_1$$), 8.01e$$-$$4 ($$R_2$$) Single time step Time loop $${\it{\Delta}} t$$ $$\| {\bf{u}}_{TH}- {\bf{u}}_{*}\|_{L^2(0,T)}$$ leap-frog 4.57e$$-$$3 s 9.15 s 1.e$$-$$4 s 1.70e$$-$$3 ($$R_1$$), 6.94e$$-$$4 ($$R_2$$) DG 2.98e$$-$$2 s 7.45 s 8.e$$-$$2 s 1.99e$$-$$3 ($$R_1$$), 8.01e$$-$$4 ($$R_2$$) Table 5 Computational times for the time loop (space polynomial degree $$q = 4$$) Single time step Time loop $${\it{\Delta}} t$$ $$\| {\bf{u}}_{TH}- {\bf{u}}_{*}\|_{L^2(0,T)}$$ leap-frog 4.57e$$-$$3 s 9.15 s 1.e$$-$$4 s 1.70e$$-$$3 ($$R_1$$), 6.94e$$-$$4 ($$R_2$$) DG 2.98e$$-$$2 s 7.45 s 8.e$$-$$2 s 1.99e$$-$$3 ($$R_1$$), 8.01e$$-$$4 ($$R_2$$) Single time step Time loop $${\it{\Delta}} t$$ $$\| {\bf{u}}_{TH}- {\bf{u}}_{*}\|_{L^2(0,T)}$$ leap-frog 4.57e$$-$$3 s 9.15 s 1.e$$-$$4 s 1.70e$$-$$3 ($$R_1$$), 6.94e$$-$$4 ($$R_2$$) DG 2.98e$$-$$2 s 7.45 s 8.e$$-$$2 s 1.99e$$-$$3 ($$R_1$$), 8.01e$$-$$4 ($$R_2$$) 6. Conclusions In this work, we have developed a new high-order DG finite element method for the temporal discretization of Cauchy problems for second-order differential equations. To show the capabilities of our scheme, we have applied it for the time integration of second-order system of equations resulting after DG semidiscretization (in space) of the elastodynamics equation. Our formulation contains suitable stabilization terms, which allowed the construction of an appropriate energy norm that naturally arose by the variational formulation of the problem. We have studied the well-posedness of the resulting scheme and proved a priori error estimates that properly depend on the local polynomial approximation degree and the local regularity of the exact solution. Our theoretical results have been confirmed by the numerical experiments carried out for both simplified test cases and on examples of practical interest. The method presented is well suited for the adaptive choice the discretization parameters employed within each time slab, namely the local time step $${\it{\Delta}} t_n$$ and the local approximation degree $$r_n$$, $$n= 1, \dots, N$$. The development of suitable a posteriori error estimators, to be employed within (space)-time adaptive strategies, is under investigation and will make the subject of future research. Acknowledgements The authors are grateful to the referees for their thorough and constructive comments that have greatly contributed to the improvement of the article. Funding Scientific Independence of young Researchers (SIR; RBSI14VT0S in part to P.F. Antonietti and I. Mazzieri); ‘PolyPDEs: Non-conforming polyhedral finite element methods for the approximation of partial differential equations’ by Italian Ministry of Education, Universities and Research (MIUR). References Adams R. A. & Fournier J. J. F. ( 2003 ) Sobolev Spaces . Pure and Applied Mathematics , vol. 140, 2nd edn. Amsterdam : Elsevier , pp. xiv+305 . Adjerid S. & Temimi H. ( 2011 ) A discontinuous Galerkin method for the wave equation. Comput. Methods Appl. Mech. Engrg. , 200 , 837 – 849 . Google Scholar Crossref Search ADS Antonietti P. F. , Ayuso de Dios B. , Mazzieri I. & Quarteroni A. ( 2016a ) Stability analysis of discontinuous Galerkin approximations to the elastodynamics problem. J. Sci. Comput. , 68 , 143 – 170 . Google Scholar Crossref Search ADS Antonietti P. F. , Ferroni A. , Mazzieri I. , Paolucci R. , Quarteroni A. , Smerzini C. & Stupazzini M. ( 2016b ) A review on numerical modelling of earthquakes through Discontinuous Galerkin Spectral Element methods ( submitted to ESAIM: Proceedings and Surveys ). Antonietti P. F. , Marcati C. , Mazzieri I. & Quarteroni A. ( 2016c ) High order discontinuous Galerkin methods on simplicial elements for the elastodynamics equation. Numer. Algorithms , 71 , 181 – 206 . Google Scholar Crossref Search ADS Antonietti P. F. , Ferroni A. , Mazzieri I. & Quarteroni A. ( 2017 ) hp-version discontinuous Galerkin approximations of the elastodynamics equation. Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2016 ( Bittencourt M. et al. eds). ecture Notes in Computational Science and Engineering 119. https://doi.org/10.1007/978-3-319-65870-4_1 . Antonietti P. F. , Mazzieri I. , Quarteroni A. & Rapetti F. ( 2012 ) Non-conforming high order approximations of the elastodynamics equation. Comput. Methods Appl. Mech. Engrg. , 209–212 , 212 – 238 . Google Scholar Crossref Search ADS Arnold D. N. ( 1982 ) An interior penalty finite element method with discontinuous elements. SIAM J. Numer. Anal. , 19 , 742 – 760 . Google Scholar Crossref Search ADS Arnold D. N. , Brezzi F. , Cockburn B. & Marini L. D. ( 2001–2002 ) Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal. , 39 , 1749 – 1779 . Google Scholar Crossref Search ADS Babuska I. & Suri M. ( 1987 ) The h-p version of the finite element method with quasiuniform meshes. ESAIM Math. Model. Numer. Anal. , 21 , 199 – 238 . Google Scholar Crossref Search ADS Baccouch M. ( 2012 ) A local discontinuous Galerkin method for the second-order wave equation. Comput. Methods Appl. Mech. Engrg. , 209–212 , 129 – 143 . Google Scholar Crossref Search ADS Banks H. T. , Birch M. J. , Brewin M. P. , Greenwald S. E. , Hu S. , Kenz Z. R. , Kruse C. , Maischak M. , Shaw S. & Whiteman J. R. ( 2014 ) High-order space-time finite element schemes for acoustic and viscodynamic wave equations with temporal decoupling. Int. J. Numer. Methods Engrg ., 98 , 131 – 156 . Google Scholar Crossref Search ADS Butcher J. C. ( 2008 ) Numerical Methods for Ordinary Differential Equations . New York : John Wiley & Sons, Ltd, pp. 137 – 316 . Canuto C. , Hussaini M. , Quarteroni A. & Zang T. ( 2006 ) Spectral Methods (Fundamental in Single Domains ). Berlin : Springer. Canuto C. G. , Hussaini M. Y. , Quarteroni A. & Zang T. A. ( 2007 ) Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics . Berlin, Heidelberg-New York : Springer . Cockburn B. & Shu C.-W. ( 1998 ) The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J. Numer. Anal. , 35 , 2440 – 2463 . Google Scholar Crossref Search ADS Cohen M. A. & Tan C. O. ( 2012 ) A polynomial approximation for arbitrary functions. Appl. Math. Lett. , 25 , 1947 – 1952 . Google Scholar Crossref Search ADS Collino F. , Fouquet T. & Joly P. ( 2003 ) A conservative space-time mesh refinement method for the 1-D wave equation. Part I: construction. Numer. Math. , 95 , 197 – 221 . Google Scholar Crossref Search ADS Courant R. , Friedrichs K. & Lewy H. ( 1928 ) Über die partiellen Differenzengleichungen der mathematischen Physik. Math. Ann. , 100 , 32 – 74 . Google Scholar Crossref Search ADS Delcourte S. & Glinsky N. ( 2015 ) Analysis of a high-order space and time discontinuous Galerkin method for elastodynamic equations application to 3D wave propagation. ESAIM Math. Model. Numer. Anal. , 49 , 1085 – 1126 . Google Scholar Crossref Search ADS Delfour M. , Hager W. & Trochu F. ( 1981 ) Discontinuous Galerkin methods for ordinary differential equations. Math.Comput. , 36 , 455 – 473 . Google Scholar Crossref Search ADS Di Pietro D. A. & Ern A. ( 2012 ) Mathematical Aspects of Discontinuous Galerkin Methods. Mathématiques et Applications , vol. 69. Berlin : Springer . Diaz J. & Grote M. J. ( 2009 ) Energy conserving explicit local time stepping for second-order wave equations. SIAM J. Sci. Comput. , 31 , 1985 – 2014 . Google Scholar Crossref Search ADS Dumbser M. , Käser M. & Toro E. ( 2007 ) An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes - V. Local time stepping and p-adaptivity. Geophys. J. Int. , 171 , 695 – 717 . Google Scholar Crossref Search ADS Ferroni A. , Antonietti P. F. , Mazzieri I. & Quarteroni A. ( 2017 ) Dispersion-dissipation analysis of 3D continuous and discontinuous spectral element methods for the elastodynamics equation. Geophys. J. Int. Available at https://doi.org/10.1093/gji/ggx384. Grote M. J. & Mitkova T. ( 2013 ) High-order explicit local time-stepping methods for damped wave equations. J. Comput. Appl. Math. , 239 , 270 – 289 . Google Scholar Crossref Search ADS Grote M. J. , Schneebeli A. & Schötzau D. ( 2006 ) Discontinuous Galerkin finite element method for the wave equation. SIAM J. Numer. Anal. , 44 , 2408 – 2431 . Google Scholar Crossref Search ADS Haskell N. A. ( 1953 ) The dispersion of surface waves on multilayered media. Bull. Seismol. Soc. Am. , 43 , 17 – 43 . Hesthaven J. & Warburton T. ( 2008 ) Nodal Discontinuous Galerkin Methods . Texts in Applied Mathematics , vol. 54. Berlin : Springer , pp. xiv+501 . Houston P. , Schwab C. & Süli E. ( 2000 ) Stabilized hp-finite element methods for first-order hyperbolic problems. SIAM J. Numer. Anal. , 37 , 1618 – 1643 . Google Scholar Crossref Search ADS Hughes T. J. & Hulbert G. M. ( 1988 ) Space-time finite element methods for elastodynamics: formulations and error estimates. Comput. Methods Appl. Mech. Engrg. , 66 , 339 – 363 . Google Scholar Crossref Search ADS Hughes T. J. & Stewart J. R. ( 1996 ) A space-time formulation for multiscale phenomena. J. Comput. Appl. Math. , 74 , 217 – 229 . Google Scholar Crossref Search ADS Hulbert G. M. & Hughes T. J. ( 1990 ) Space-time finite element methods for second-order hyperbolic equations. Comput. Methods Appl. Mech. Engrg. , 84 , 327 – 348 . Google Scholar Crossref Search ADS Johnson C. ( 1993 ) Discontinuous Galerkin finite element methods for second order hyperbolic problems. Comput. Methods Appl. Mech. Engrg. , 107 , 117 – 129 . Google Scholar Crossref Search ADS Kroopnick A. ( 1999 ) Bounded and l2-solutions to a second order nonlinear differential equation with a square integrable forcing term. Int. J. Math. Sci. , 22 , 569 – 571 . Google Scholar Crossref Search ADS Lesaint P. & Raviart P. ( 1974 ) On a Finite Element Method for Solving the Neutron Transport Equation . Publications mathématiques et informatique de Rennes. S4 : 1 – 40 . LeVeque R. ( 2007 ) Finite Difference Methods for Ordinary and Partial Differential Equations . Philadelphia, PA, USA : SIAM - Society for Industrial and Applied Mathematics. Mazzieri I. , Stupazzini M. , Guidotti R. & Smerzini C. ( 2013 ) Speed: spectral elements in elastodynamics with discontinuous Galerkin: a non-conforming approach for 3D multi-scale problems. Int. J. Numer. Methods Engrg. , 95 , 991 – 1010 . Google Scholar Crossref Search ADS Quarteroni A. ( 2014 ) Numerical Models for Differential Problems . MS&A. Modeling, Simulation and Applications, vol. 8. Italia, Milan : Springer. Quarteroni A. , Sacco R. & Saleri F. ( 2007 ) Numerical Mathematics . Texts in Applied Mathematics , vol. 37 , 2nd edn. Berlin : Springer , pp. xviii+655 . Quarteroni A. & Valli A. ( 2008 ) Numerical Approximation of Partial Differential Equations . Springer Series in Computational Mathematics. Berlin, Heidelberg : Springer. Reed W. H. & Hill T. R. ( 1973 ) Triangular mesh methods for the neutron transport equation. Technical Report LA-UR-73-479 . Los Alamos, New Mexico : Los Alamos Scientific Laboratory. Rivière B. ( 2008 ) Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations . Philadelphia, PA, USA : Society for Industrial and Applied Mathematics. Rivière B. & Wheeler M. F. ( 2003 ) Discontinuous finite element methods for acoustic and elastic wave problems. Contemporary Mathematics , 329 , 271 – 282 . Google Scholar Crossref Search ADS Schötzau D. & Schwab C. ( 2000 ) Time discretization of parabolic problems by the hp-version of the discontinuous Galerkin finite element method. SIAM J. Numer. Anal. , 38 , 837 – 875 . Google Scholar Crossref Search ADS Schwab C. ( 1998 ) p- and hp- Finite Element Methods. Theory and Applications in Solid and Fluid Mechanics . New York : Oxford University Press. Stacey R. ( 1988 ) Improved transparent boundary formulations for the elastic-wave equation. Bull. Seismol. Soc. Am. , 78 , 2089 – 2097 . Taube A. , Dumbser M. , Munz C.-D. & Schneider R. ( 2009 ) A high-order discontinuous Galerkin method with time-accurate local time stepping for the Maxwell equations. Int. J. Numer. Model. El. , 22 , 77 – 103 . Google Scholar Crossref Search ADS Thompson L. L. & Pinsky P. M. ( 1996 ) A space-time finite element method for structural acoustics in infinite domains part 1: Formulation, stability and convergence. Comput. Methods Appl. Mech. Engrg. , 132 , 195 – 227 . Google Scholar Crossref Search ADS van der Vegt J. J. W. , Klaij C. M. , van der Bos F. & van der Ven H. ( 2006 ) Space-time discontinuous Galerkin method for the compressible Navier-Stokes equations on deforming meshes. In Wesseling P. Onate E. & Periaux J. (eds), Proceedings European Conference on Computational Fluid Dynamics ECCOMAS CFD 2006 . Delft : TU Delft. Walkington N. J. ( 2014 ) Combined dg–cg time stepping for wave equations. SIAM J. Numer. Anal. , 52 , 1398 – 1417 . Google Scholar Crossref Search ADS Werder T. , Gerdes K. , Schötzau D. & Schwab C. ( 2001 ) hp-discontinuous Galerkin time stepping for parabolic problems. Comput. Methods Appl. Mech. Engrg. , 190 , 6685 – 6708 . Google Scholar Crossref Search ADS Wheeler M. F. ( 1978 ) An elliptic collocation-finite element method with interior penalties. SIAM J. Numer. Anal. , 15 , 152 – 161 . Google Scholar Crossref Search ADS Wilcox L. C. , Stadler G. , Burstedde C. & Ghattas O. ( 2010 ) A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media. J. Comput. Phys. , 229 , 9373 – 9396 . Google Scholar Crossref Search ADS Yang Y. , Chirputkar S. , Alpert D. N. , Eason T. , Spottswood S. & Qian D. ( 2012 ) Enriched space time finite element method: a new paradigm for multiscaling from elastodynamics to molecular dynamics. Int. J. Numer. Methods Engrg. , 92 , 115 – 140 . Google Scholar Crossref Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png IMA Journal of Numerical Analysis Oxford University Press

# A high-order discontinuous Galerkin approximation to ordinary differential equations with applications to elastodynamics

, Volume 38 (4) – Oct 16, 2018
26 pages

/lp/ou_press/a-high-order-discontinuous-galerkin-approximation-to-ordinary-U7gY0v1GJr
Publisher
Oxford University Press
ISSN
0272-4979
eISSN
1464-3642
D.O.I.
10.1093/imanum/drx062
Publisher site
See Article on Publisher Site

### Abstract

Abstract The aim of this work is to propose and analyse a new high-order discontinuous Galerkin finite element method for the time integration of a Cauchy problem for second-order ordinary differential equations. These equations typically arise after space semidiscretization of second-order hyperbolic-type differential problems, e.g., wave, elastodynamics and acoustics equations. After introducing the new method, we analyse its well-posedness and prove a priori error estimates in a suitable (mesh-dependent) norm. Numerical results are also presented to verify our theoretical estimates. 1. Introduction In this article, we develop a high-order discontinuous Galerkin finite element method for the numerical approximation of ordinary differential equations that arise after space semidiscretization of second-order hyperbolic problems. The applications we have in mind include, e.g., acoustic, elastic and electromagnetic wave propagation phenomena. Traditional approaches for the numerical integration of (second-order) ordinary differential systems rely on implicit and explicit finite difference, Runge–Kutta and Newmark schemes (see, e.g., LeVeque, 2007; Quarteroni et al., 2007; Butcher, 2008, for a comprehensive review). In many engineering applications, explicit methods are in general preferred to implicit ones. In fact, although the latter are unconditionally stable, the former are computationally cheaper. The main drawback of explicit methods is the time step limitation imposed by the Courant–Friedrichs–Lewy (CFL) condition. Such a constraint, which depends in general on the space discretization parameters and the media properties, can severely affect the computational efficiency. A possible way to alleviate this limitation is to introduce suitable local time-stepping (LTS) algorithms (Collino et al., 2003; Diaz & Grote, 2009; Grote & Mitkova, 2013), using a small time step only when needed. Another possibility is to adopt an explicit LTS method by extending the so-called arbitrary high-order derivatives discontinuous Galerkin (ADER-DG) approach (Dumbser et al., 2007; Taube et al., 2009). In this context, a proper time step can be tailored for each element of the time mesh. However, to correctly propagate the wave field from one element to the other, an additional (computationally demanding) synchronization process has to be operated. Differently than above, here, we derive an implicit arbitrarily high-order accurate time integration scheme based on a DG approach. DG methods (Reed & Hill, 1973; Lesaint & Raviart, 1974) have been firstly proposed to approximate (in space) hyperbolic problems (Reed & Hill, 1973) and then generalized to elliptic and parabolic equations (Wheeler, 1978; Arnold, 1982); see also Arnold et al. (2001–2002), Houston et al. (2000), Cockburn & Shu (1998), Rivière (2008), Hesthaven & Warburton (2008) and Di Pietro & Ern (2012). Relevant applications and analysis of DG schemes for the scalar wave equation can be found in Rivière & Wheeler (2003), Grote et al. (2006) and Baccouch (2012), whereas for elastodynamics problems we refer the reader to Dumbser et al. (2007), Wilcox et al. (2010), Antonietti et al. (2012, 2016a,b,c, 2017), Mazzieri et al. (2013), Delcourte & Glinsky (2015) and Ferroni et al. (2017). The DG approach has also been used to solve initial-value problems. In time-dependent problems, the information follows the positive direction of time and solutions are casual (they depend on the past but not on the future). In contrast with finite difference time integration schemes, for which the solution at the current step depends upon the previous steps, time-discontinuous Galerkin methods applied over time slabs $$[t_n, t_{n+1}]$$ lead to a casual system in which the solution at the current time slab depends only upon the solution at $$t_n^-$$. By coupling discontinuous Galerkin discretizations in both space and time leads to a fully space-time finite element formulation. Relevant works on this topic concern both parabolic and hyperbolic problems (see, e.g., Delfour et al., 1981; Werder et al., 2001; van der Vegt et al., 2006). For the latter, space-time finite elements are typically built upon reformulating the original problem as a system of first-order equations (see, e.g., Hughes & Hulbert, 1988; Johnson, 1993; Banks et al., 2014). The latter can be seen as the result of space semidiscretization of first-order hyperbolic problems or even second-order hyperbolic problems in which the problem is formulated in terms of the displacement (respectively velocity) and the stress (respectively strain) tensor fields. To the best of our knowledge, only a few recent results regarding finite element approximations of second-order differential systems are available in the literature (Hulbert & Hughes, 1990; Thompson & Pinsky, 1996; Adjerid & Temimi, 2011; Yang et al., 2012; Walkington, 2014). In Adjerid & Temimi (2011), a new DG approach based on the solution of the scalar wave equation (and higher-order differential equations) has been proposed and analysed. The stabilization terms appearing in Adjerid & Temimi (2011) introduced to penalize the jumps of the solution and its derivative across different time slabs are similar to those proposed in Hughes & Hulbert (1988), Hulbert & Hughes (1990) and Thompson & Pinsky (1996), where a Galerkin least square (GLS) approach is applied to stabilize the numerical scheme and prove its convergence. As an extension of the space-time formulation of Hughes & Stewart (1996), in Yang et al. (2012), the authors present an enriched version of the space-time finite element method to incorporate in the same model multiple temporal scale features. A combination of continuous and discontinuous Galerkin time-stepping approach is used in Walkington (2014) to develop arbitrary order approximation for second-order hyperbolic problems. Stability, convergence and accuracy are proved for scalar wave propagation with nonhomogeneous boundary conditions. In this work, a new DG method for the solution of systems of second-order ordinary differential equations is presented. The resulting weak formulation is obtained by imposing the continuity of tractions and velocities across time slabs weakly, without adding any extra GLS stabilization term. We show that the proposed formulation, in which the displacement field is the only unknown, is well posed and we prove a priori stability and error estimates in a suitable mesh-dependent norm. The obtained time-discontinuous scheme results in an implicit and unconditionally stable method. Moreover, allowing independent displacement interpolations between different time slabs, this method is naturally suited for an adaptive choice of the time discretization parameters, i.e., the use of high-order polynomials/small time steps only when the solution features a sharp temporal derivative. The article is organized as follows. In Section 2, we formulate the problem, discretize it and analyse its well-posedness. Then, we derive the corresponding algebraic system of equations. In Section 3, we carry out the convergence analysis providing suitable stability and error estimates. The application of the proposed method to the elastodynamics equations is described in Section 4. Here, the space-time finite element formulation is obtained combining the discontinuous Galerkin spectral element (DGSE) spatial discretization proposed in Antonietti et al. (2012) to the one presented here for the time integration. Numerical results are shown in Section 5. Throughout the article, we denote by $$\|\mathbf{a}\|$$ the Euclidean norm of a vector $$\mathbf{a} \in \mathbb{R}^d$$, $$d \geq 1$$ and by $$\|A\|_{\infty} = \max_{i=1, \dots, m} \sum_{j = 1}^{n}|a_{ij}|$$, the $$\ell^\infty$$-norm of a matrix $$A \in \mathbb{R}^{m\times n}$$, $$m,n \geq 1$$. Moreover, $$C$$ denotes a generic positive constant that may take different values in different places, but is always independent of the discretization parameters. The notation $$x \lesssim y$$ means $$x \leq C y$$ for a constant $$C$$ as before. For a given $$I \subset \mathbb{R}$$, for any $$v: I \rightarrow \mathbb{R}$$ we denote by $$L^p(I)$$ and $$H^p(I)$$, $$p \in \mathbb{N} \setminus 0$$ the usual Lebesgue and Hilbert spaces, respectively, and endow them with the usual norms (see Adams & Fournier, 2003). For $$p=0$$, we write $$L^2(I)$$ in place of $$H^0(I)$$. Finally, we use boldface type for vectorial functions. More precisely, the Lebesgue and Hilbert spaces of vector-valued functions are denoted by $${\bf L}^{p}(I)=[L^{p}(I)]^{d}$$ and $${\bf H}^{p}(I)=[H^{p}(I)]^{d}$$, respectively, $$d \geq 1$$. 2. A model problem and its discontinuous Galerkin approximation In this section, we introduce a high-order discontinuous Galerkin finite element method for second-order ordinary differential equations, prove its well-posedness and provide its algebraic form. 2.1 Problem statement and its DG discretization For $$T > 0$$, we consider the following model problem: find $$\mathbf{u}: (0,T] \rightarrow \mathbb{R}^d, d \geq 1$$, such that $$\ddot{\mathbf{u}}(t) + L\dot{\mathbf{u}}(t) + K{\mathbf{u}}(t) = \mathbf{f}(t) \qquad \forall t \in (0, T],$$ (2.1) where $$L, K \in \mathbb{R}^{d\times d}$$ are symmetric, positive definite matrices and $$\mathbf{f} \in {\bf{L}}^2(0,T]$$. We supplement problem (2.1) with the following initial conditions: \begin{align}\mathbf{u}(0) = \mathbf{u}_{0}, \qquad \dot{\mathbf{u}}(0) =\mathbf{u}_{1},\end{align} (2.2) where $$\mathbf{u}_0, \mathbf{u}_1 \in \mathbb{R}^d$$. Problem (2.1) is well posed and admits a unique solution $$\textbf{u} \in {\bf{H}}^2(0,T]$$ in the interval $$(0,T]$$ (see, e.g., Kroopnick, 1999). We partition the interval $$I = (0, T]$$ into $$N$$ time slabs $$I_n = (t_{n-1}, t_n]$$ having length $${\it{\Delta}} t_n = t_n - t_{n-1}$$, for $$n = 1,..,N$$, with $$t_0 = 0$$ and $$t_N=T$$, as it is shown in Fig. 1. Fig. 1. View largeDownload slide Example of time domain partition (top). Zoom of the time domain partition: values $$t_n^+$$ and $$t_n^-$$ are also reported (bottom). Fig. 1. View largeDownload slide Example of time domain partition (top). Zoom of the time domain partition: values $$t_n^+$$ and $$t_n^-$$ are also reported (bottom). In the following we will use the notation: \begin{equation*}({\bf{u}},{\bf{v}})_{I} = \int_{I} {\bf{u}}(s) \cdot {\bf{v}}(s)\,{\rm d}s, \quad \quad \langle {\bf{u}}, {\bf{v}} \rangle_t = {\bf{u}}(t) \cdot {\bf{v}}(t), \end{equation*} where $$\mathbf{a} \cdot \mathbf{b}$$ indicates the Euclidean scalar product between two vectors $$\mathbf{a}, \mathbf{b} \in \mathbb{R}^d$$. To deal with discontinuous functions, we also define, for (a regular enough) $$\mathbf{v}$$, the jump operator at $$t_n$$ as \begin{align*}[\mathbf{v}]_{n} &= \mathbf{v}(t_{n}^+) - \mathbf{v}(t_{n}^-) \quad {\rm for} \; n \geq 0, \end{align*} where \begin{align*}\mathbf{v}(t_{n}^\pm) &= \lim_{\epsilon \rightarrow 0^\pm} \mathbf{v}(t_{n} + \epsilon) \quad {\rm for} \; n \geq 0. \end{align*} Implicit with the above definition is that \begin{align*}[\mathbf{v}]_{0} = \mathbf{v}(0^+) - \mathbf{u}_0,\qquad [\dot{\mathbf{v}}]_{0} = \dot{\mathbf{v}}(0^+) - \mathbf{u}_1. \end{align*} Moreover, we use the symbols $$\mathbf{v}_{n}^+ = \mathbf{v}(t_{n}^+)$$ and $$\mathbf{v}_{n}^- = \mathbf{v}(t_{n}^-)$$ to represent the trace of (a regular enough) $$\mathbf{v}$$, taken within the interior of $$I_{n+1}$$ and $$I_{n}$$, respectively (cf. Fig. 1). Next, following a time integration approach, we incrementally build (on $$n$$) an approximation of the exact solution $$\mathbf{u}$$ on each time slab $$I_n$$. For this reason, we focus on the generic interval $$I_n$$, and assume the solution on $$I_{n-1}$$ to be known. Note that $$\textbf{u} \in {\bf{H}}^2(I_n)$$. If we multiply equation (2.1) by $$\dot{\textbf{v}}(t)$$, being $${\mathbf{v}}(t)$$ a regular enough function, we obtain $$(\ddot{\mathbf{u}},\dot{\mathbf{v}})_{I_n} + (L \dot{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n} + (K{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n} = ( \mathbf{f} , \dot{\mathbf{v}})_{I_n}.$$ (2.3) Finally, since $$\mathbf{u} \in {\bf{H}}^2(0,T)$$, we observe that $$[\mathbf{u}]_n = [\dot{\mathbf{u}}]_n = \textbf{0}$$, for $$n=1, \dots, N$$, we rewrite (2.3) adding suitable (strongly consistent) terms $$(\ddot{\mathbf{u}},\dot{\mathbf{v}})_{I_n} + (L \dot{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n} + (K{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n} + [\dot{\mathbf{u}}]_{n-1} \cdot \dot{\mathbf{v}}_{n-1}^+ + K [\mathbf{u}]_{n-1} \cdot \mathbf{v}_{n-1}^+ = ( \mathbf{f} , \dot{\mathbf{v}})_{I_n}.$$ (2.4) Next, we introduce the local finite dimensional space \begin{align*}{{\bf{V}}_n^r} = \left\{\mathbf{v}: I_n \rightarrow \mathbb{R}^d: \mathbf{v} \in [\mathbb{P}^r(I_n)]^d \right\}\!, \end{align*} where $$\mathbb{P}^r(I_n)$$ is the space of polynomials of degree not greater than or equal to $$r \geq 2$$ on $$I_n$$. Then, introducing $$\textbf{r}=(r_1,...,r_N) \in \mathbb{N}^N$$ the polynomial degree vector, we can define the DG finite element space as \begin{align*}{\mathcal{V}^\bf{r}} = \left\{\mathbf{v} \in {\bf{L}}^2(0,T): \mathbf{v}|_{I_n} \in {{\bf{V}}_n^{r_n}} \forall \, n = 1, \dots, N \right\}\!, \end{align*} whose dimension is $$\sum_{n=1}^N (r_n+1)d$$. Summing over all time slabs in (2.4), we are able to define the bilinear form $${\mathcal A} : {{\mathcal{V}^\bf{r}} \times {\mathcal{V}^\bf{r}}} \rightarrow \mathbb{R}$$ \begin{align}{\mathcal A}(\mathbf{u}, \mathbf{v}) & = \sum_{n=1}^N \Big( (\ddot{\mathbf{u}},\dot{\mathbf{v}})_{I_n} + (L \dot{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n} + (K{\mathbf{u}}, \dot{\mathbf{v}} )_{I_n}\Big) + \sum_{n=1}^{N-1} \left( \dot{[\mathbf{u}]}_{n} \cdot \dot{\mathbf{v}}_{n}^+ + K [\mathbf{u}]_{n} \cdot \mathbf{v}_{n}^+ \right)\nonumber\\ &\quad{} + \dot{\mathbf{u}}_{0}^+ \cdot \dot{\mathbf{v}}_{0}^+ + K \mathbf{u}_{0}^+ \cdot \mathbf{v}_{0}^+,\end{align} (2.5) and the linear functional $$F : {{\mathcal{V}^\bf{r}}} \rightarrow \mathbb{R}$$ as \begin{align}F (\mathbf{v}) & = \sum_{n=1}^N ( \mathbf{f} , \dot{\mathbf{v}} )_{I_n} + {\mathbf{u}}_{1} \cdot \dot{\mathbf{v}}_{0}^+ + K \mathbf{u}_{0} \cdot \mathbf{v}_{0}^+,\end{align} (2.6) respectively. Notice that for $$n = 1$$, we implicitly adopted the convention that $$\mathbf{u}_0^- = \mathbf{u}_0$$ and $$\dot{\mathbf{u}}_0^- = \mathbf{u}_1$$. Moreover, we point out that the definition of the bilinear form $${\mathcal A}(\cdot, \cdot)$$ makes sense whenever its arguments are, at least, $$H^2(I_n)-$$ functions for any $$n=1,\ldots,N$$. The DG formulation of problem (2.1)–(2.2) reads as follows: find $$\mathbf{u}_{\textrm{DG}} \in {\mathcal{V}^\bf{r}}$$ such that \begin{align}{\mathcal A}(\mathbf{u}_{\small{DG}}, \mathbf{v}) = F(\mathbf{v}) \qquad \forall \, \mathbf{v}\in {\mathcal{V}^\bf{r}}.\end{align} (2.7) The existence and uniqueness of the discrete solution $$\textbf{u}_{\textrm{DG}} \in {\mathcal{V}^\bf{r}}$$ of problem (2.7) is a direct consequence of the following result. Proposition 2.1 The function $$\left\||{{\cdot}}\right\|| : {\mathcal{V}^\bf{r}} \rightarrow \mathbb{R}^+$$, defined as \begin{align}\left\|| {{\mathbf{v}}} \right\||^2 & = \sum_{n=1}^N \|L^\frac{1}{2} \dot{\mathbf{v}}\|_{{\bf{L}}^2(I_n)}^2 + \frac{1}{2} (\dot{\mathbf{v}}_{0}^+)^2 + \frac{1}{2}\sum_{n=1}^{N-1} \left(\dot{[\mathbf{v}]}_{n} \right)^2 + \frac{1}{2}(\dot{\mathbf{v}}_{T}^-)^2 \nonumber\\ &\quad{} + \frac{1}{2} (K^\frac{1}{2}{\mathbf{v}}_{0}^+)^2 + \frac{1}{2} \sum_{n=1}^{N-1} \left(K^\frac{1}{2} [\mathbf{v}]_{n} \right)^2 +\frac{1}{2} (K^\frac{1}{2}{\mathbf{v}}_{T}^-)^2,\end{align} (2.8) is a norm on $${\mathcal{V}^\bf{r}}$$. Proof Since the homogeneity and the subadditivity properties are satisfied, we just need to show that \begin{align*}\left\||{\mathbf{v}}\right\|| = 0 \Leftrightarrow \mathbf{v} = \mathbf{0}. \end{align*} The sufficient condition is trivial. We then show that $$\left\||{\mathbf{v}}\right\||= 0$$ implies $$\mathbf{v} = \textbf{0}$$. Clearly, if $$\left\||{\mathbf{v}}\right\|| = 0$$ then all the terms on the right-hand side of (2.8) are null. In particular, from the facts that $$K$$ and $$L$$ are positive definite and $$\|L^{\frac{1}{2}}\dot{\mathbf{v}} \|_{{\bf{L}}^2(I_0)} = 0$$ and $$( K^{\frac{1}{2}} \mathbf{v}_0^+ )^2 = 0$$, we conclude that $$\mathbf{v}$$ is such that \begin{align*}\begin{cases}\dot{\mathbf{v}}(t) = \mathbf{0} \qquad \forall \,t \in I_0 ,\\ \mathbf{v}_0^+ = \mathbf{0}. \end{cases}\end{align*} Therefore, we have that $$\mathbf{v}(t)|_{I_0} = \mathbf{0}$$. We now proceed by induction, and consider the interval $$I_n$$, supposing $$\mathbf{v}(t)|_{I_{n-1}} = \mathbf{0}$$. From $$\left\||{\mathbf{v}}\right\||= 0$$ we have $$(K^\frac{1}{2} [\mathbf{v}]_{n} )^2 = 0$$, which yields $$\mathbf{v}_{n}^+ = 0$$. Then, we infer that $$\mathbf{v}$$ is such that \begin{align*}\begin{cases}\dot{\mathbf{v}}(t) = \mathbf{0} \qquad \forall \,t \in I_n ,\\ \mathbf{v}_n^+ = \mathbf{0}. \end{cases}\end{align*} As a result, we conclude that $$\mathbf{v}$$ is the null function on any interval $$I_n$$, $$n = 0, \dots, N-1$$, and the proof is complete. □ The following result follows from Proposition 2.1. Remark 2.2 Taking $$\mathbf{u} = \mathbf{v}$$ in (2.5) and integrating by parts, we obtain \begin{equation*}{\mathcal A}(\mathbf{v}, \mathbf{v}) = \left\||{ \mathbf{v}}\right\||^2 \quad \forall \, \mathbf{v}\in {\mathcal{V}^\bf{r}}, \end{equation*} i.e., the bilinear form $${\mathcal A}(\cdot, \cdot)$$ defined in (2.5) is coercive with respect to the norm $$\left\||{\cdot}\right\||$$, with coercivity constant $$\alpha = 1$$. Therefore, the following result holds. Proposition 2.3 Problem (2.7) admits a unique solution $$\textbf{u}_{\textrm{DG}} \in {\mathcal{V}^\bf{r}}$$. Proof The thesis follows directly from Proposition 2.1, the bilinearity of $${\mathcal A}(\cdot, \cdot)$$ and the linearity of $$F(\cdot)$$. □ 2.2 Algebraic formulation We derive here the algebraic formulation corresponding to problem (2.4) for the time slab $$I_n$$, $$n\geq 1$$, where a local degree $${r_n}$$ is employed. We remark that the employment of discontinuous functions between a node $$t_n$$ allows to compute the solution of the problem separately for one time slab at a time. Indeed, the weak formulation (2.7) restricted to $${{\bf{V}}_n^{{r_n}}}$$ reads as: find $${\mathbf{u}_n} = {\textbf{u}_{\textrm{DG}}}_{|_{I_n}} \in {{\bf{V}}_n^{{r_n}}}$$ such that \begin{align*}& (\ddot{\bf u}_n,\dot{\mathbf{v}})_{I_n} + (L \dot{\bf u}_n, \dot{\mathbf{v}} )_{I_n} + (K{{{\mathbf{u}_n}}}, \dot{\mathbf{v}} )_{I_n} + \langle \dot{\bf u}_n, \dot{\mathbf{v}}\rangle_{t_{n-1}^+}+ \langle K {{\mathbf{u}_n}} , \mathbf{v}\rangle_{t_{n-1}^+}\\ &\quad{} = ( \mathbf{f} , \dot{\mathbf{v}})_{I_n} + \dot{\bf u}_n(t_{n-1}^-) \cdot \dot{\mathbf{v}}_{n-1}^+ + K {{\mathbf{u}_n}}(t_{n-1}^-) \cdot \mathbf{v}_{n-1}^+ \quad \forall \; \mathbf{v} \in {{\bf{V}}_n^{{r_n}}}, \end{align*} where on the right-hand side the values $$\dot{\bf u}_n(t_{n-1}^-)$$ and $${{\mathbf{u}_n}}(t_{n-1}^-)$$ computed for $$I_{n-1}$$ are used as initial conditions for the current slab. Notice that, for the slab $$I_1$$, we set $$\dot{\bf u}_1(t_{0}^-) = {\bf{u}}_1$$ and $${\bf u}_1(t_{0}^-) = {\bf{u}}_0$$. Focusing on the generic interval $$I_n$$, we introduce a basis $$\{\psi^\ell(t) \}_{\ell=1, \dots, {r_n}+1}$$ for the polynomial space $$\mathbb{P}^{r_n}_n(I_n)$$ and define $$D= d({r_n} + 1)$$ the dimension of the local finite dimensional space $${{\bf{V}}_n^{{r_n}}}$$. We also introduce the vectorial basis $$\{ {\bf {\it{\Psi}}_i^\ell}(t) \}_{i=1, \dots, d}^{\ell=1, \dots, {r_n}+1}$$, where $${\bf {\it{\Psi}}_i^\ell}(t)$$ is the $$d$$-dimensional vector whose $$i$$th component is $$\psi^\ell$$ and other components are zero. Using the notation just introduced, we can write the trial function $${\mathbf{u}_n}$$ as a linear combination of the basis functions, i.e., \begin{align*}{\mathbf{u}_n}(t) = \sum\limits_{j=1}^d \sum\limits_{m=1}^{{r_n}+1} \alpha_j^m {\bf {\it{\Psi}}_j^m}(t), \end{align*} where $$\alpha_j^m \in \mathbb{R}$$ for $$j = 1, \dots, d$$ and $$m=1, \dots, {r_n}+1$$. Next, we write equation (2.7) for any test function $${\bf {\it{\Psi}}_i^\ell}(t), i = 1, \dots, d, \ell = 1, \dots, {r_n}+1$$, obtaining the algebraic system \begin{align*}A \mathbf{U}_n = \mathbf{F}_n, \end{align*} where on the interval $$I_n$$, $$\mathbf{U}_n \in \mathbb{R}^D$$ is the solution vector, $$\mathbf{F}_n \in \mathbb{R}^D$$ corresponds to the data and is given componentwise as \begin{equation*}(\mathbf{F}_n)_i^\ell = ( \mathbf{f} , \dot{\bf {\it{\Psi}}}_i^\ell )_{I_n} + \dot{\bf u}_n(t_{n-1}^-) \cdot \dot{\bf {\it{\Psi}}}_i^\ell(t_{n-1}^+) + K {{\mathbf{u}_n}}(t_{n-1}^-) \cdot \bf{{\it{\Psi}}}_i^\ell(t_{n-1}^+) \quad {\rm for } \; i = 1, \dots, d, \ell = 1, \dots, {r_n}+1, \end{equation*} and $$A \in \mathbb{R}^{D\times D}$$ is the related local stiffness matrix. We next investigate the structure of the matrix $$A$$ by defining the following local matrices for $$\ell, m = 1, \dots, {r_n}+1$$, \begin{align*}& M^1_{\ell m} = ( \ddot{\psi}^m , \dot{\psi}^\ell )_{I_n}, \qquad M^{2}_{\ell m} = (\dot{\psi}^m ,\dot{\psi}^\ell )_{I_n}, \qquad M^3_{\ell m} = ( {\psi}^m, \dot\psi^\ell )_{I_n}, \\ & M^4_{\ell m} = \langle \dot{\psi}^m, \dot{\psi}^\ell \rangle_{t_{n-1}^+}, \quad M^{5}_{\ell m} = \langle \psi^m, \psi^\ell \rangle_{t_{n-1}^+}. \end{align*} Setting \begin{align*}M & = M^1 + M^4, \\ B_{ij} & = L_{ij} M^{2} + K_{ij} \left( M^3 + M^{5} \right)\!, \qquad i,j=1,...,d, \end{align*} with $$M, B_{ij} \in \mathbb{R}^{({r_n}+1)\times({r_n}+1)}$$ for any $$i,j = 1, \dots,d$$, we can rewrite the matrix $$A$$ as \begin{align}A = \begin{vmatrix}M & 0 & \cdots & 0\\ 0 & M & & \vdots \\ \vdots & & \ddots & 0 \\ 0 & \cdots & 0 & M \end{vmatrix}+ \begin{vmatrix}B_{11} & B_{12} & \cdots & B_{1d}\\ B_{21} & B_{22} & & \vdots \\ \vdots & & \ddots & \\ B_{d1} & \cdots & & B_{dd}\end{vmatrix}.\end{align} (2.9) Its sparsity clearly depends on the sparsity pattern of the matrices $$L$$ and $$K$$. 3. Convergence analysis In this section, we first analyse the stability of the DG method (2.7) and then we prove a priori error estimates. 3.1 Stability analysis Proposition 3.1 Let $$\mathbf{f} \in \textbf{L}^2(0,T]$$ and $$\mathbf{u}_0,\mathbf{u}_1 \in \mathbb{R}^d$$. Then, the solution of (2.7), $${\mathbf{u}_{\textrm{DG}}} \in {\mathcal{V}^\bf{r}}$$, satisfies $$\left\||{ {\mathbf{u}_{\textrm{DG}}} }\right\||\lesssim \left( \|L^{-\frac{1}{2}}\mathbf{f}\|_{{\bf{L}}^2(0,T)}^2 + (K^\frac{1}{2} \mathbf{u}_0 )^2 + (\mathbf{u}_1 )^2 \right)^{\frac{1}{2}}.$$ (3.1) Proof Using the definition (2.8) and arithmetic–geometric inequalities, we obtain \begin{align*}\left\||{{\mathbf{u}_{\textrm{DG}}}}\right\||^2 & = {\mathcal A}({\mathbf{u}_{\textrm{DG}}},{\mathbf{u}_{\textrm{DG}}}) = \sum_{n=1}^N ( \mathbf{f} , \dot{\mathbf{u}}_{\textrm{DG}} )_{I_n} + {\mathbf{u}}_{1} \cdot \dot{\mathbf{u}}_{\textrm{DG}}(0^+) + K \mathbf{u}_{0} \cdot {\mathbf{u}}_{\textrm{DG}}(0^+) \\ & \leq \frac{1}{2}\sum_{n=1}^N \left( \| L^{-\frac{1}{2}} \mathbf{f}\|_{{\bf{L}}^2(I_n)}^2 + \| L^{\frac{1}{2}} \dot{\mathbf{u}}_{\textrm{DG}} \|_{{\bf{L}}^2(I_n)}^2 \right) \\ & \quad + {\mathbf{u}}_{1}^2 + \frac{1}{4} \dot{\mathbf{u}}_{\textrm{DG}}^2(0^+) + (K^\frac{1}{2} \mathbf{u}_{0})^2 + \frac{1}{4} \left(K^\frac{1}{2}{\mathbf{u}}_{\textrm{DG}}(0^+)\right)^2, \end{align*} i.e., \begin{align*}\frac{1}{2}\left\||{ \textbf{u}_{\textrm{DG}}^2}\right\|| \leq \frac{1}{2}\sum_{n=1}^N \| L^{-\frac{1}{2}} \mathbf{f}\|_{{\bf{L}}^2(I_n)}^2 + {\mathbf{u}}_{1}^2 + (K^\frac{1}{2} \mathbf{u}_{0})^2. \end{align*} □ 3.2 Error estimates In this section, we derive a priori error estimates in the $$\left\||{{\cdot}}\right\||$$ norm defined in (2.8). We start by introducing some preliminary definitions and results. Definition 3.2 (Schötzau & Schwab, 2000, Definition 3.1) Let $$I = (-1, 1)$$. For a function $$w \in L^2(I)$$ which is continuous at $$t =1$$ we define the $$L^2$$-projection (with boundary preserving value) $${\mathcal{P}}^r w \in \mathbb{P}^r(I)$$, $$r \geq 0$$, by the $$r + 1$$ conditions \begin{align}({\mathcal{P}}^r w - w)(1) &= 0, \\ \end{align} (3.2a) \begin{align}\int_I (w -{\mathcal{P}}^r w) q \,{\rm d}t &= 0 \quad \forall q \in \mathbb{P}^{r-1}(I).\end{align} (3.2b) For $$r =0$$, only the first condition (3.2a) is necessary. $${\mathcal{P}}^rw$$ is well posed and satisfies optimal error estimates, see Schötzau & Schwab (2000, Theorem 3.10). Definition 3.3 Let $$I = (-1, 1)$$. For a function $$u \in H^1(I)$$ such that $$u'$$ is continuous at $$t=1$$, we define $${\it{\Pi}}^r u \in \mathbb{P}^r(I)$$, $$r \geq 1$$, by $${\it{\Pi}}^r u(t) = u(-1) + \int_{-1}^t {\mathcal{P}}^{r-1}u'(s)\,{\rm d}s \quad \forall t \in (-1,1),$$ (3.3) where $${\mathcal{P}}^{r-1}u'$$ is given in Definition 3.2. In the following, we denote by $$\{L_i\}_{i\geq 0}$$, $$L_i \in \mathbb{P}^i(I)$$, the Legendre polynomials on $$I = (-1, 1)$$, defined by $$L_0(t) = 1$$, $$L_1(t) = t$$ and $$\int_{-1}^t L_i(s)\,{\rm d}s =\frac{1}{2i+1}(L_{i+1} (t) - L_{i-1}(t))$$ for $$i \geq 1$$, see Canuto et al. (2006). Lemma 3.4 The projection $${\it{\Pi}}^r u$$ in Definition 3.3 is well defined. Proof We consider the Legendre series $$u' = \sum_{i=0}^\infty b_i L_i$$ with coefficients $$b_i \in \mathbb{R}$$. From Schötzau & Schwab (2000, Lemma 3.5), we find that $${\mathcal{P}}^{r-1}u'(t) = \sum_{i=0}^{r-2} b_i L_i(t) + \left( \sum_{i=r-1}^{\infty} b_i \right)L_{r-1}(t)$$ (3.4) and consequently, $${\it{\Pi}}^r u(t) = u(-1) + \sum_{i=0}^{r-2} b_i \int_{-1}^t L_i(s) \,{\rm d}s + \left( \sum_{i=r-1}^{\infty} b_i \right)\int_{-1}^t L_{r-1}(s) \,{\rm d}s.$$ (3.5) Then, \begin{align*}{\it{\Pi}}^r u(t) & = \left( u(-1) + b_0 -\frac{b_1}{3}\right)L_0(t) + \left(b_0 - \frac{b_2}{5} \right)L_1(t) + \sum_{i=2}^{r-3} \left(\frac{b_{i-1}}{2i-1} - \frac{b_{i+1}}{2i+3} \right)L_i(t) \\ &\quad{} + \left(\frac{b_{r-3}}{2r-5} - \sum_{i=r-1}^{\infty}\frac{b_{i}}{2r-1} \right)L_{r-2}(t) + \frac{b_{r-2}}{2r-3}L_{r-1}(t) + \sum_{i=r-1}^{\infty}\frac{b_i}{2r-1}L_r(t), \end{align*} or, equivalently, $${\it{\Pi}}^r u (t) = \sum_{i=0}^r u_i^* L_i(t),$$ (3.6) with \begin{eqnarray}u_0^* & = & u(-1) +b_0 -\frac{b_1}{3}, \\ \end{eqnarray} (3.7) \begin{eqnarray}u_1^* & = & b_0 - \frac{b_2}{5}, \\ \end{eqnarray} (3.8) \begin{eqnarray}u_i^* & = &\frac{b_{i-1}}{2i-1} - \frac{b_{i+1}}{2i+3} \quad {\rm for \;} i = 2,...,r-3,\\ \end{eqnarray} (3.9) \begin{eqnarray}u_{r-2}^* & = & \frac{b_{r-3}}{2r-5} - \sum_{i=r-1}^{\infty}\frac{b_{i}}{2r-1},\\ \end{eqnarray} (3.10) \begin{eqnarray}u_{r-1}^* & = & \frac{b_{r-2}}{2r-3},\\ \end{eqnarray} (3.11) \begin{eqnarray}u_r^* & = & \sum_{i=r-1}^{\infty}\frac{b_i}{2r-1} . \end{eqnarray} (3.12) The well-posedness of $${\it{\Pi}}^r u$$ follows from (3.6) and the orthogonality of the Legendre polynomials $$L_i$$. □ Now, we analyse the properties of the projection $${\it{\Pi}}^r u$$. Lemma 3.5 Let $$I = (-1, 1)$$. For any $$u \in H^1(I),$$ such that $$u'$$ is continuous at $$t=1$$ we have, for $$r \geq 1$$, \begin{equation*}\| {\it{\Pi}}^r u \|_{L^2(I)} \lesssim |u(1)| + |u'(1)| + \| u' \|_{L^2(I)}. \end{equation*} Proof We expand $$u'$$ as $$u' = \sum_{i=0}^\infty b_i L_i$$ with coefficients $$b_i \in \mathbb{R}$$; then, $$u$$ can be written as $$u(t) = u(-1) + \sum_{i=0}^\infty b_i \int_{-1}^t L_i(s) \,{\rm d}s$$. By orthogonality of the Legendre polynomials, it easily follows that $$u(1) = u(-1) + 2b_0.$$ (3.13) Next, from (3.7)–(3.12) and (3.13), we have that \begin{align}\| {\it{\Pi}}^r u \|_{L^2(I)}^2 & = \sum_{i=0}^r \frac{2}{2i+1} {u_i^*}^2 \nonumber \\ & = 2\left(u(-1)+b_0 -\frac{b_1}{3}\right)^2 + \frac{2}{3}\left(b_0 -\frac{b_2}{5}\right)^2 + \sum_{j=2}^{r-3}\frac{2}{2j+1} \left(\frac{b_{j-1}}{2j-1} - \frac{b_{j+1}}{2j+3} \right)^2 \nonumber\\ &\quad{} + \frac{2}{2r-3}\left(\frac{b_{r-3}}{2r-5}- \sum_{j=r-1}^{\infty} \frac{b_j}{2r-1}\right)^2 + \frac{2}{2r-1}\left(\frac{b_{r-2}}{2r-3}\right)^2 + \frac{2}{2r+1}\left(\sum_{j=r-1}^{\infty} \frac{b_{j}}{2r-1}\right)^2 \nonumber\\ &\leq 4u(1)^2 + 4b_0^2 + \frac{4}{3}b_1^2 + \frac{4}{3}b_0^2 + \frac{4}{75}b_2^2 + 2 \sum_{j=2}^{r-3}\frac{2}{2j+1}\frac{b_{j-1}^2}{(2j-1)^2} + 2\sum_{j=2}^{r-3}\frac{2}{2j+1}\frac{b_{j+1}^2}{(2j+3)^2} \nonumber \\ & \quad + \frac{4}{2r-3}\frac{b_{r-3}^2}{(2r-5)^2} + \frac{4}{2r-3}\left(\sum_{j=r-1}^{\infty} \frac{b_j}{2r-1}\right)^2 + \frac{2}{2r-1}\left(\frac{b_{r-2}}{2r-3}\right)^2\nonumber \\ &\quad + \frac{2}{2r+1}\left(\sum_{j=r-1}^{\infty} \frac{b_{j}}{2r-1}\right)^2.\end{align} (3.14) Then, we observe that \begin{gather*}\frac{4}{2r-3}\left(\sum_{j=r-1}^{\infty} \frac{b_j}{2r-1}\right)^2 \lesssim \frac{1}{r^3}\left(\sum_{j=r-1}^{\infty} b_j\right)^2 \lesssim |u'(1)|^2, \\ \frac{2}{2r+1}\left(\sum_{j=r-1}^{\infty} \frac{b_j}{2r-1}\right)^2 \lesssim \frac{1}{r^3}\left(\sum_{j=r-1}^{\infty} b_j\right)^2 \lesssim |u'(1)|^2,\\ 4b_0^2 + 2\sum_{j=2}^{r-3}\frac{2}{2j+1}\frac{b_{j-1}^2}{(2j-1)^2} + \frac{4}{2r-3}\frac{b_{r-3}^2}{(2r-5)^2} + \frac{2}{2r-1}\frac{b_{r-2}^2}{(2r-3)^2} \lesssim \sum_{j=0}^{r-2}\frac{2}{2j+1}b_j^2 \leq \| u' \|_{L^2(I)}^2, \end{gather*} and similarly \begin{equation*}\frac{4}{3}b_0^2 + \frac{4}{3}b_1^2 + \frac{4}{75}b_2^2 + 2\sum_{j=2}^{r-3}\frac{2}{2j+1}\frac{b_{j+1}^2}{(2j+3)^2} \lesssim \| u' \|_{L^2(I)}^2. \end{equation*} Using the above estimates into (3.14) yields $\| {\it{\Pi}}^r u \|_{L^2(I)}^2 \lesssim |u(1)|^2 + |u'(1)|^2 + \| u' \|_{L^2(I)}^2.$ □ Lemma 3.6 Under the same assumptions on $$u$$ made in Lemma 3.5 it holds \begin{align}({\it{\Pi}}^r u - u)(-1) &= 0,\\ \end{align} (3.15a) \begin{align}({\it{\Pi}}^r u - u)(1) &= 0, \\ \end{align} (3.15b) \begin{align}({\it{\Pi}}^r u - u)'(1) &= 0,\\ \end{align} (3.15c) \begin{align}\int_I (u -{\it{\Pi}}^r u)' q {\,\rm d}t &= 0 \quad \forall q \in \mathbb{P}^{r-2}(I).\end{align} (3.15d) Proof From Definition 3.3 we obtain (3.15a). Property (3.15b) can be proved as follows: \begin{equation*}{\it{\Pi}}^r u(1) = u(-1) + \sum_{i=0}^{r-2} b_i \int_{-1}^1 L_i(s) {\,\rm d}s + \left( \sum_{i=r-1}^{\infty} b_i \right)\int_{-1}^1 L_{r-1}(s) {\,\rm d}s = u(-1) + 2b_0 = u(1). \end{equation*} Equality (3.15c) follows by taking the derivative with respect to $$t$$ in (3.5), evaluating the latter at $$t=1$$ and using (3.4). Finally, (3.15d) comes from \begin{equation*}\int_I (u -{\it{\Pi}}^r u)' L_k{\,\rm d}t = \int_I (u' -{\mathcal{P}}^{r-1} u') L_k {\,\rm d}t = \int_I \left(\sum_{i=r}^{\infty} b_i L_i(t) - \left( \sum_{i=r}^{\infty} b_i \right)L_{r-1}\right)L_k{\,\rm d}t = 0, \end{equation*} for all $$k=0,...,r-2$$. □ On an arbitrary interval $$I = (a, b)$$ with $${\it{\Delta}} t = b - a > 0$$, we define $${\it{\Pi}}^r_{I}$$ via the linear map $${T} : (-1, 1) \rightarrow (a, b)$$, $$\xi \rightarrow y = \frac{1}{2} (a + b + \xi {\it{\Delta}} t )$$ as $${\it{\Pi}}^r_{I} u = [{\it{\Pi}}^r (u \circ {T})] \circ {T}^{-1},$$ (3.16) from which it follows that \begin{align}({\it{\Pi}}^r_I u - u)(a) &= 0,\\ \end{align} (3.17a) \begin{align}({\it{\Pi}}^r_I u - u)(b) &= 0, \\ \end{align} (3.17b) \begin{align}({\it{\Pi}}^r_I u - u)'(b) &= 0,\\ \end{align} (3.17c) \begin{align}\int_I (u -{\it{\Pi}}^r_I u)' q {\,\rm d}t &= 0 \quad \forall q \in \mathbb{P}^{r-2}(I).\end{align} (3.17d) Remark 3.7 Reasoning as in the proof of Lemma 3.4, on the generic interval $$I=(a,b)$$, $${\it{\Pi}}^r_I$$ can be written as \begin{equation*}{\it{\Pi}}^r_I u = \sum_{i=0}^{r} \widetilde{u}_i^* \widetilde{L}_i = u(a) + \sum_{i=0}^{r-2} \widetilde{b}_i \int_{a}^t \widetilde{L}_i(s) {\,\rm d}s + \left( \sum_{i=r-1}^{\infty} \widetilde{b}_i \right)\int_{a}^t \widetilde{L}_{r-1}(s) {\,\rm d}s, \end{equation*} where $$\widetilde{L}_i$$ represents the $$i$$th shifted Legendre polynomial (cf. Cohen & Tan, 2012), so that $$u = \sum_{i=0}^{\infty} \widetilde{u}_i \widetilde{L}_i$$, and the coefficients $$\widetilde{u}_{i}^*$$ are defined as in (3.7)–(3.12), provided that $$b_j$$ is replaced by $$\widetilde{b}_j$$, for any $$j$$. We recall the following approximation results from Babuska & Suri (1987), see also Schwab (1998) and Canuto et al. (2007). Proposition 3.8 Let $$I = (a, b)$$ with $${\it{\Delta}} t = b - a > 0$$. For every $$u \in H^s(I)$$ there exist a constant C, independent of $$u,r,{\it{\Delta}} t$$ and a sequence $$\{{Q}^r u\}_{r\geq 1}$$, with each $${Q}^r u \in \mathbb{P}^r (I)$$ such that, for any $$0 \leq q \leq s$$, \begin{align}\| u - {Q}^ru \|_{H^q(I)} & \lesssim \frac{{\it{\Delta}} t^{\mu-q}}{r^{s-q}} \| u \|_{H^s(I)} \quad s\geq 0, \\ \end{align} (3.18) \begin{align}\| u - {Q}^ru \|_{L^2(\partial I)} & \lesssim \frac{{\it{\Delta}} t^{\mu-\frac{1}{2}}}{r^{s-\frac{1}{2}}} \| u \|_{H^s(I)} \quad s> \frac{1}{2}, \\ \end{align} (3.19) \begin{align}\| u - {Q}^ru \|_{H^1(\partial I)} & \lesssim \frac{{\it{\Delta}} t^{\mu-\frac{3}{2}}}{r^{s-\frac{3}{2}}} \| u \|_{H^s(I)} \quad s> \frac{3}{2},\end{align} (3.20) where $$\mu = \min(r+1,s).$$ Next, we define the projector $$\mathcal{P}^{r}_Iu$$ on an arbitrary interval $$I = (a, b)$$ with $${\it{\Delta}} t = b - a > 0$$ via the linear map $${T} : (-1, 1) \rightarrow (a, b)$$, $$\xi \rightarrow y = \frac{1}{2} (a + b + \xi {\it{\Delta}} t )$$ as $$\mathcal{P}^r_{I} u = [\mathcal{P}^r (u \circ {T})] \circ {T}^{-1},$$ (3.21) for which the following approximation result holds, see Schötzau & Schwab (2000). Proposition 3.9 Let $$I = (a, b)$$ with $${\it{\Delta}} t = b - a > 0$$. For any $$u \in H^s(I)$$, $$s \geq 2$$, there holds \begin{align}\| u' - \mathcal{P}^{r-1}_Iu' \|_{L^2(I)} & \lesssim \frac{{\it{\Delta}} t^{\mu-1}}{r^{s-1}} \| u \|_{H^s(I)},\end{align} (3.22) where $$\mu = \min(r+1,s).$$ The following property holds. Lemma 3.10 Let $$I = (a, b)$$ with $${\it{\Delta}} t = b - a > 0$$. For any $$u \in H^s(I)$$, $$s\geq 2$$, let $${\it{\Pi}}^r_{I} u$$, $$r \geq 2$$ be the projector defined as in (3.16). Then, \begin{align}\| (u- {\it{\Pi}}_{I}^r u)' \|_{L^2(I)} & \lesssim \frac{{\it{\Delta}} t^{\mu-1}}{r^{s-1}} \| u \|_{H^s(I)} \\ \end{align} (3.23) \begin{align}\| u- {\it{\Pi}}_{I}^r u \|_{L^2(I)} & \lesssim \frac{{\it{\Delta}} t^{\mu}}{r^{s-1}} \| u \|_{H^s(I)},\end{align} (3.24) where $$\mu = \min(r+1,s).$$ Proof We start by proving (3.23), exploiting the definition (3.3) which directly yields $$({\it{\Pi}}_{I}^r u)' = \mathcal{P}_{I}^{r-1} (u')$$ \begin{align}\| (u- {\it{\Pi}}_{I}^r u)' \|_{L^2(I)} = \| u'- ({\it{\Pi}}_{I}^r u)' \|_{L^2(I)}= \| u'- \mathcal{P}_{I}^{r-1} (u') \|_{L^2(I)} \lesssim \frac{{\it{\Delta}} t^{\mu-1}}{r^{s-1}} \| u \|_{H^s(I)},\end{align} (3.25) where in the last inequality we have used (3.22). As concerns (3.24), employing the definition (3.3), we obtain \begin{align*}\| u- {\it{\Pi}}_{I}^r u \|_{L^2(I)}^2 = & \left\| u(a) + \int_a^tu'(s){\rm d}s - u(a) - \int_a^t \mathcal{P}_{I}^{r-1} (u')(s){\,\rm d}s \right\|_{L^2(I)}^2 \\ = & \left\| \int_a^t u'(s)\,{\rm d}s - \int_a^t \mathcal{P}_{I}^{r-1} (u')(s){\,\rm d}s \right\|_{L^2(I)}^2 = \int_a^b \left( \int_a^t(u'(s) - \mathcal{P}_{I}^{r-1} (u'))(s){\,\rm d}s \right)^2{\,\rm d}t \\ \leq & \int_a^b \left( \int_a^t\big|u'(s) - \mathcal{P}_{I}^{r-1} (u')(s)\big|{\,\rm d}s \right)^2{\,\rm d}t \\ \leq & \int_a^b \left( \left(\int_a^t {\,\rm d}s \right)^{\frac{1}{2}} \left(\int_a^t\big|u'(s) - \mathcal{P}_{I}^{r-1} (u')(s)\big|^2{\,\rm d}s\right)^{\frac{1}{2}} \right)^2{\,\rm d}t\\ = & \int_a^b \left( t-a \right) \left(\int_a^t\big|u'(s) - \mathcal{P}_{I}^{r-1} (u')(s)\big|^2{\,\rm d}s\right) {\,\rm d}t \\ \leq & \int_a^b {\it{\Delta}} t \left(\int_a^b\big|u'(s) - \mathcal{P}_{I}^{r-1} (u')(s)\big|^2{\,\rm d}s\right) {\,\rm d}t = {\it{\Delta}} t^2 \int_a^b \big|u'(s) - \mathcal{P}_{I}^{r-1} (u')(s)\big|^2{\,\rm d}s \\ = & {\it{\Delta}} t^2 \| u'- \mathcal{P}_{I}^{r-1} (u') \|_{L^2(I)}^2 \lesssim {\it{\Delta}} t^2 \frac{{\it{\Delta}} t^{2\mu-2}}{r^{2s-2}} \| u \|_{H^s(I)}^2, \end{align*} from which (3.24) directly follows. □ As a direct consequence of the above Lemma we have the following result. Corollary 3.11 Let $$I= (a,b)$$ with $${\it{\Delta}} t = b-a$$ and let $$u \in H^s(I)$$, with $$s \geq 2$$. There holds \begin{align}\| (u- {\it{\Pi}}^r_{I} u)'' \|_{L^2(I)} & \lesssim \frac{{\it{\Delta}} t^{ \mu - 2}}{r^{s-3}} \| u \|_{H^s(I)},\end{align} (3.26) where $$\mu = \min(r+1,s)$$. Proof From the triangle inequality, interpolation estimates (3.18) and (3.22), we get \begin{align*}\| (u- {\it{\Pi}}^r_{I} u)'' \|_{L^2(I)}&\leq \| (u- {Q}^r u)'' \|_{L^2(I)}+ \| ({Q}^r u- {\it{\Pi}}^r_{I} u)'' \|_{L^2(I)} \\ &\lesssim \frac{{\it{\Delta}} t^{\mu-2}}{r^{s-2}} \| u \|_{H^s(I)}+ \frac{r^2}{{\it{\Delta}} t}\| ({Q}^r u- {\it{\Pi}}^r_{I} u)' \|_{L^2(I)}\\ &= \frac{{\it{\Delta}} t^{\mu-2}}{r^{s-2}} \| u \|_{H^s(I)}+ \frac{r^2}{{\it{\Delta}} t} \| ({Q}^r u)'- {\mathcal{P}}^{r-1}_I u' \|_{L^2(I)}\\ &\lesssim \frac{{\it{\Delta}} t^{\mu-2}}{r^{s-2}} \| u \|_{H^s(I)}+ \frac{r^2}{{\it{\Delta}} t} \| (u - {Q}^r u)'\|_{L^2(I)}+ \frac{r^2}{{\it{\Delta}} t} \|u' - {\mathcal{P}}^{r-1}_I u' \|_{L^2(I)}\\ &\lesssim \frac{{\it{\Delta}} t^{\mu-2}}{r^{s-2}} \| u \|_{H^s(I)}+ \frac{r^2}{{\it{\Delta}} t} \frac{{\it{\Delta}} t^{\mu-1}}{r^{s-1}} \| u \|_{H^s(I)}+ \frac{r^2}{{\it{\Delta}} t} \frac{{\it{\Delta}} t^{\mu-1}}{r^{s-1}} \| u \|_{H^s(I)}\\ &\lesssim \frac{{\it{\Delta}} t^{\mu-2}}{r^{s-3}} \| u \|_{H^s(I)}, \end{align*} where, in the second step, we have also used the inverse inequality $$\| \phi \|_{H^1(I)} \lesssim {\it{\Delta}} t^{-1}r^{2} \| \phi \|_{L^2(I)}$$ valid for any $$\phi\in \mathbb{P}^r(I)$$ (cf. Schwab, 1998), e.g., □ before proving the main result of this section, we notice that if $$\mathbf{v} \in {\bf{H}}^s(I)$$, for $$s \geq 2$$, then its projection $${\it{\Pi}}^{r}_I {\mathbf{v}} \in [\mathbb{P}^r(I)]^d$$ is defined, componentwise, according to (3.16) and the same kind of estimates given in (3.24), (3.23) and (3.26) hold. Now, let $$I = (0,T)$$, and let $$\mathbf{u} \in {\bf{H}}^s(0,T]$$, $$s \geq 2$$, be the solution of (2.3) and let $${\it{\Pi}}^\mathbf{r}_{I} \mathbf{u} \in {\mathcal{V}^\bf{r}}$$ be its projection such that $${\it{\Pi}}^\mathbf{r}_{I} \mathbf{u}_{|_{I_n}} = {\it{\Pi}}^{r_n}_{I_n} \mathbf{u} \in {{\bf{V}}_n^{r_n}}$$ is defined according to (3.16), for $$n=1,...,N$$. Then, it is easy to see that for $$n = 1,...,N$$ it holds \begin{align}(\mathbf{u} - {\it{\Pi}}^\textbf{r}_{I})(t_n) = 0, & \\ \end{align} (3.27) \begin{align}(\mathbf{u} - {\it{\Pi}}^\textbf{r}_{I})'(t_n) = 0, & \\ \end{align} (3.28) \begin{align}((\mathbf{u} - {\it{\Pi}}^\textbf{r}_{I})', q)_{I_n} = 0 & \qquad \forall \; q \in \mathbb{P}^{r_n-2}(I_n).\end{align} (3.29) Finally, we observe that (2.7) is strongly consistent (cf. Quarteroni, 2014). Indeed, it holds \begin{align}\mathcal{A}( \mathbf{u} - \mathbf{u}_{\textrm{DG}}, \mathbf{v} ) = 0 \qquad \forall \; \mathbf{v} \in {\mathcal{V}^\bf{r}}.\end{align} (3.30) We are now ready to prove the following convergence result. Theorem 3.12 Let $$\mathbf{u}$$ be the solution of (2.1)–(2.2) and let $${\mathbf{u}_{\textrm{DG}}} \in {\mathcal{V}^\bf{r}}$$ be its DG finite element approximation based on employing formulation (2.7). If $$\mathbf{u}|_{I_n} \in {\bf{H}}^{s_n}(I_n)$$, for any $$n=1,...,N$$, with $$s_n \geq 2$$, then \begin{align}\left\||{\mathbf{u} - {\mathbf{u}_{\textrm{DG}}}}\right\|| \lesssim \sum_{n=1}^N \frac{{\it{\Delta}} t_n^{{\mu_n} - \frac{3}{2}}}{ r_n^{s_n-3}} \|\mathbf{u}\|_{{\bf{H}}^{s_n}(I_n)},\end{align} (3.31) where $$\mu_n = \min(r_n+1,s_n)$$ for any $$n=1,...,N$$, and where the hidden constant depends on the norm of the matrices $$L$$ and $$K$$. Proof We define the error on the interval $$I=(0,T)$$ as $$\mathbf{e} = \mathbf{u} - \mathbf{u}_{\textrm{DG}}$$, and we split it as $$\mathbf{e} = {\mathbf{e}^{\pi}} + {\mathbf{e}^h}$$, where $${\mathbf{e}^{\pi}}$$ is the projection error while $${\mathbf{e}^h} \in {\mathcal{V}^\bf{r}}$$ is the remainder, i.e., $${\mathbf{e}^{\pi}} = \mathbf{u} - {\it{\Pi}}^{\textbf{r}}_I \textbf{u}$$ and $${\mathbf{e}^h} = {\it{\Pi}}^{\textbf{r}}_I \textbf{u} - \mathbf{u}_{\textrm{DG}}$$. Clearly, we have $$\left\||{\mathbf{e}}\right\|| \leq \left\||{\textbf{e}^\pi}\right\|| + \left\||{\textbf{e}^h}\right\||.$$ By exploiting properties (3.27)–(3.29) and by employing estimates (3.23) and (3.26), we can bound the term $$\left\||{\textbf{e}^\pi}\right\||$$ as \begin{align*}\left\||{\textbf{e}^\pi}\right\||^2 & = \sum_{n=1}^N \|L^\frac{1}{2} {\dot{\mathbf{e}}^{\pi}} \|^2_{L^2(I_n)} + \frac{1}{2}\sum_{n=1}^N {\dot{\mathbf{e}}^{\pi}}(t_{n-1}^+)^2 = \sum_{n=1}^N \|L^\frac{1}{2} {\dot{\mathbf{e}}^{\pi}} \|^2_{L^2(I_n)} + \frac{1}{2}\sum_{n=1}^N \left( - \int_{t_{n-1}}^{t_n} {\ddot{\mathbf{e}}^{\pi}}(s) \,{\,\rm d} s \right)^2 \nonumber\\ & \lesssim \sum_{n=1}^N \left( \|{\dot{\mathbf{e}}^{\pi}} \|^2_{L^2(I_n)} + {{\it{\Delta}} t_n} \|{\ddot{\mathbf{e}}^{\pi}} \|^2_{L^2(I_n)}\right) \nonumber \lesssim \sum_{n=1}^N \frac{{\it{\Delta}} t_n^{2\mu_n-3}}{r^{2s_n-{6}}} \| {\bf{u}} \|_{H^{s_n}(I_n)}^2, \end{align*} with $$\mu_n = \min(r_n+1,s_n)$$ for any $$n=1,...,N$$. For the term $$\left\||{\textbf{e}^h}\right\||$$, the Galerkin orthogonality the Galerkin orthogonality (3.30) yields \begin{align*}\left\||{\textbf{e}^h}\right\||^2 & = {\mathcal A}( \textbf{e}^h, \textbf{e}^h) = - {\mathcal A}( \textbf{e}^\pi, \textbf{e}^h)\\ & = - \sum_{n=1}^N \Big( ( {\ddot{\mathbf{e}}^{\pi}} ,{\dot{\mathbf{e}}^h})_{I_n} + (L {\dot{\mathbf{e}}^{\pi}}, {\dot{\mathbf{e}}^h})_{I_n} +(K{\mathbf{e}^{\pi}}, {\dot{\mathbf{e}}^h} )_{I_n} + [ {\dot{\mathbf{e}}^{\pi}} ]_n \cdot {\dot{\mathbf{e}}^h} (t_n^+) + K[{\mathbf{e}^{\pi}}]_n \cdot {\mathbf{e}^h}(t_n^+) \Big)\\ & \quad- {\dot{\mathbf{e}}^{\pi}} (t_0^+) \cdot {\dot{\mathbf{e}}^h} (t_0^+) - K{\mathbf{e}^{\pi}}(t_0^+) \cdot {\mathbf{e}^h}(t_0^+). \end{align*} Integrating by parts the term $$( {\ddot{\mathbf{e}}^{\pi}} ,{\dot{\mathbf{e}}^h})_{I_n}$$ and rearranging the addends, we obtain \begin{align*}\left\||{\textbf{e}^h}\right\||^2 & = \sum_{n=1}^N \Big( ( {\dot{\mathbf{e}}^{\pi}} ,{\ddot{\mathbf{e}}^h})_{I_n} - (L {\dot{\mathbf{e}}^{\pi}}, {\dot{\mathbf{e}}^h})_{I_n} -(K{\mathbf{e}^{\pi}}, {\dot{\mathbf{e}}^h} )_{I_n} \Big) \\ & \qquad + \sum_{n=1}^{N-1} \Big( [ {\dot{\mathbf{e}}^h} ]_n \cdot {\dot{\mathbf{e}}^{\pi}} (t_n^-) - K [{\mathbf{e}^{\pi}}]_n \cdot {\mathbf{e}^h}(t_n^+) \Big) - {\dot{\mathbf{e}}^{\pi}} (T^-) \cdot {\dot{\mathbf{e}}^h} (T^-) - K{\mathbf{e}^{\pi}}(0^+) \cdot {\mathbf{e}^h}(0^+). \end{align*} Using now properties (3.27)–(3.29) into the above equation yields \begin{align*}\left\||{\textbf{e}^h}\right\||^2 & = \sum_{n=1}^N \Big( - (L {\dot{\mathbf{e}}^{\pi}}, {\dot{\mathbf{e}}^h})_{I_n} - (K{\mathbf{e}^{\pi}}, {\dot{\mathbf{e}}^h} )_{I_n} \Big). \end{align*} Applying Cauchy–Schwarz and arithmetic–geometric inequalities and using estimates (3.24) and (3.23), it is easy to obtain \begin{align*}\left\||{\textbf{e}^h}\right\||^2 & \lesssim \sum_{n=1}^N \Big( \|{\dot{\mathbf{e}}^{\pi}}\|_{L^2(In)}^2 + \|{\mathbf{e}^{\pi}}\|_{L^2(In)}^2 \Big) \lesssim { \sum_{n=1}^N \frac{{\it{\Delta}} t_n^{2\mu_n-2}}{{r_n}^{2s_n-{2}}}\| {\bf{u}} \|_{H^{s_n}(I_n)}^2, }\end{align*} where $$\mu_n = \min(r_n+1,s_n)$$ for any $$n=1,...,N$$, and where the hidden constant depends on the norms $$\|L^\frac{1}{2}\|_{\infty}$$ and $$\|L^{-\frac{1}{2}}K \|_{\infty}$$. The thesis follows by combining the above bound with the estimate on $$\left\||{\textbf{e}^\pi}\right\||$$. □ Remark 3.13 In the case when $${\it{\Delta}} t_n = {\it{\Delta}} t$$, $$r_n = r\geq 2$$ and $$s_n = s \geq 2$$, for $$n=1,...,N$$, the error bound of Theorem 3.12 becomes \begin{align*}\left\||{\mathbf{u} - {\mathbf{u}_{\textrm{DG}}}}\right\|| \lesssim \frac{{\it{\Delta}} t^{r-\frac{1}{2}}}{ r^{s-{3}} }\|\mathbf{u}\|_{{\bf{H}}^{s}(0,T)}, \end{align*} where the hidden constant depends on the norm of the matrices $$L$$ and $$K$$. 4. Application to the elastodynamics problem In this section, we apply the method presented in Section 2 to the simulation of elastic wave propagations in heterogeneous media. In particular, we will adopt the DG method previously developed to handle the time integration of the second-order differential system arising after space discretization obtained with the DGSE approach proposed in Antonietti et al. (2012). Since the focus of this article is on time integration, here, we simply report the mathematical model of linear viscoelasticity and the algebraic linear system resulting after the space discretization obtained employing the DGSE scheme of Antonietti et al. (2012). Finally, some numerical results are discussed. 4.1 Mathematical modelling of seismic wave propagations and its algebraic formulation For a given open-bounded domain $${\it{\Omega}} \subset \mathbb{R}^m$$, $$m=2,3$$, we consider the following problem: for $$T > 0$$ find $$\mathbf{u}: {\it{\Omega}} \times [0,T] \rightarrow \mathbb{R}^m$$ such that \begin{align}\displaystyle \rho \partial_{tt} {\bf{u}} + 2\rho \zeta \partial_{t}{\mathbf{u}} + \rho \zeta^2 \mathbf{u} -{\nabla \cdot{{\bf{\sigma}}}} & = {\bf{f}}, && \textrm{in} \; {\it{\Omega}} \times (0,T], \\ \end{align} (4.1a) \begin{align}\displaystyle {\bf{u}} & = {\bf{0}}, && \textrm{on} \; {\it{\Gamma}}_D\times (0,T], \\ \end{align} (4.1b) \begin{align}\displaystyle {\bf{\sigma}} {\bf{n}} & = {\bf{0}}, && \textrm{on} \; {\it{\Gamma}}_N\times (0,T], \\ \end{align} (4.1c) \begin{align}\displaystyle \partial_t {\bf{u}} & = {\bf{u}}_1, && \textrm{in} \; {\it{\Omega}} \times \{ 0\}, \\ \end{align} (4.1d) \begin{align}\displaystyle {\bf{u}} &= {\bf{u}}_0, &&\textrm{in} \; {\it{\Omega}} \times \{ 0\},\end{align} (4.1e) where $$\overline{\partial {\it{\Omega}}} = \overline{{\it{\Gamma}}_D} \cup \overline{{\it{\Gamma}}_N}$$ with $${\it{\Gamma}}_D \cap {\it{\Gamma}}_N =\emptyset$$, $$\mathbf{f} \in L^2 ((0,T]; {\bf{L}}^2({\it{\Omega}}))$$ is the source term, and $$\rho \in L^\infty({\it{\Omega}})$$ is such that $$\rho = \rho(\mathbf{x}) > 0$$ for almost any $$\mathbf{x} \in {\it{\Omega}}$$. The function $$\zeta = \zeta ({\bf{x}}) > 0$$ is a decay factor whose dimension is the inverse of time and is a piecewise constant function, then $$\zeta \in L^{\infty}({\it{\Omega}})$$. We suppose the stress tensor $${\bf{\sigma}}$$ to be related to the strain tensor $${\bf{\varepsilon}}(\mathbf{u}) = \frac{1}{2}\bigl( \nabla\mathbf{u} + \nabla\mathbf{u}^T \bigr)$$ through the Hooke’s law, i.e., \begin{align*}{\bf{\sigma}} = {2\mu}{\bf{\varepsilon}} + \lambda {\text{tr}}({\bf{\varepsilon}})\mathbb{I}, \end{align*} where $$\lambda = \lambda({\bf{x}})$$ and $$\mu=\mu({\bf{x}})$$ are the Lamé elastic coefficients of the medium, $${\text{tr}}(\cdot)$$ is the trace operator and $$\mathbb{I}\in \mathbb{R}^{m \times m}$$ is the identity tensor. Here and in the following, we will suppose $$\rho, \lambda$$ and $$\mu$$ to be uniformly bounded functions in $${\it{\Omega}}$$, i.e., $$\rho, \lambda, \mu \in L^{\infty}({\it{\Omega}})$$. The semidiscretization of problem (4.1a)–(4.1e) by the DGSE technique (Antonietti et al., 2012) results into the following second-order differential system for the nodal displacement $${\bf{U}}(t) \,{=} [{\bf{U}}^1(t),..., {\bf{U}}^m(t) ]^T$$ \begin{align}\begin{cases}M \ddot{{\bf{U}}}(t) + C \dot{{\bf{U}}}(t) + D{\bf{U}}(t) + A {\bf{U}}(t) = {\bf{F}}(t), & t \in (0,T], \\ \dot{\textbf{U}}(0) = {\bf{u}}_1, & \\ \textbf{U}(0) = {\bf{u}}_0, & \\ \end{cases}\end{align} (4.2) where $$\ddot{{\bf{U}}}(t)$$ (respectively, $$\dot{{\bf{U}}}(t)$$) represents the vector of nodal acceleration (respectively velocity) and $${\bf{F}}(t)$$ the vector of externally applied loads, i.e., $${\bf{F}}(t)= [{\bf{F}}^1(t),..., {\bf{F}}^m(t)]^T$$. Here, $$M$$ is the diagonal mass matrix, $$C$$ and $$D$$ stand for the structural damping (and have the same structure of $$M$$) while the action of the stiffness matrix $$A$$ to the displacement vector $${\bf{U}}$$ represents the internal (elastic) forces. We remark that the matrix $$A$$ in (4.2) is symmetric and positive definite (see Antonietti et al., 2012). The latter properties are also verified by the matrices $$M, C$$ and $$D$$. To rewrite (4.2) in the form (2.1), we multiply it by $$M^{-\frac{1}{2}}$$ and set $$\mathbf{Z}(t) = M^{\frac{1}{2}} \mathbf{U}(t)$$ and obtain \begin{align}\begin{cases}\ddot{\mathbf{Z}}(t) + L \dot{\mathbf{Z}}(t) + K \mathbf{Z}(t) = \mathbf{G}(t), & t \in (0,T], \\ \dot{\textbf{Z}}(0) = M^{\frac{1}{2}} {\bf{u}}_1, & \\ {\textbf{Z}}(0) = M^{\frac{1}{2}} {\bf{u}}_0, & \\ \end{cases}\end{align} (4.3) where $$L = M^{-\frac{1}{2}}C M^{-\frac{1}{2}}, \, K = M^{-\frac{1}{2}}(D + A) M^{-\frac{1}{2}}$$, are symmetric and positive definite and $$\mathbf{G}(t) = M^{-\frac{1}{2}}\mathbf{F}(t)$$. 5 Numerical results In this section, we present some numerical results to highlight the behavior of the DG method introduced in Section 2. In particular, we first focus our attention on a scalar benchmark, then we test it on the system of differential equations arisen by the space discretization of the elastodynamics problem. 5.1 Numerical results for a scalar problem We consider the interval $$I=(0,T]$$, $$T=10$$, subdivided into $$N$$ time slabs $$I_n$$, for $$n=1,..., N$$, of uniform length $${\it{\Delta}} t$$. Moreover, we suppose the polynomial approximation vector $$\textbf{r}$$ to be constant for each time slab, i.e., $$r_1=...=r_N=r\geq 2$$. We consider the following problem \begin{align}\begin{cases}\ddot{u}(t) + 5\dot{u}(t) + 6u(t) = 0 \qquad \forall t \in (0, 10], \\ u(0) = 2, \\ \dot{u}(0) = -5, \end{cases}\end{align} (5.1) whose exact solution is $$u(t) = e^{-3t} + e^{-2t}, \, t \in [0, 10].$$ We compute the error $$\left\||{u-u_{\textrm{DG}}}\right\||$$ versus $$1/{\it{\Delta}} t$$, with $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,2,3$$, for several polynomial degrees $$r = 2,3,4,5$$. The computed errors for the scalar test case are reported in Fig. 2, in log-log scale. As expected, a convergence rate of $${\it{\Delta}} t^{r - \frac{1}{2}}$$ is observed (cf. Theorem 3.12). This is also confirmed by the results shown in Table 1, where the computed convergence rates are reported. Although we do not have a theoretical proof, for completeness, we also report in Table 2 the computed error with respect to the $$\| \cdot \|_{L^2(0,T)}$$ and $$\| \cdot \|_{H^1(0,T)}$$ norms and the corresponding rates of convergence. The method achieves an optimal rate of convergence with respect to both norms, i.e., the error measured in the $$L^2$$-norm (respectively $$H^1$$-norm) decays as $${\it{\Delta}} t^r$$ (respectively $${\it{\Delta}} t^{r-1}$$) for $${\it{\Delta}} t$$ going to $$0$$. In Fig. 3, we display the computed errors versus the polynomial degree for different values of $${\it{\Delta}} t$$. As expected, an exponential convergence rate is observed. Fig. 2. View largeDownload slide Scalar test case. Computed errors $$\left\||{ u - u_{\textrm{DG}}}\right\||$$ versus $$1/{\it{\Delta}} t$$ for $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,2,3$$ (loglog scale) and $$r = 2,3,4,5$$. Fig. 2. View largeDownload slide Scalar test case. Computed errors $$\left\||{ u - u_{\textrm{DG}}}\right\||$$ versus $$1/{\it{\Delta}} t$$ for $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,2,3$$ (loglog scale) and $$r = 2,3,4,5$$. Fig. 3. View largeDownload slide Scalar test case. Computed errors $$\left\||{ u - u_{\textrm{DG}}}\right\||$$ versus $$r = 2,3,4,5$$ (semilogy scale) for $${\it{\Delta}} t = 0.025, 0.0125, 0.00625$$. Fig. 3. View largeDownload slide Scalar test case. Computed errors $$\left\||{ u - u_{\textrm{DG}}}\right\||$$ versus $$r = 2,3,4,5$$ (semilogy scale) for $${\it{\Delta}} t = 0.025, 0.0125, 0.00625$$. Table 1 Scalar test case. Computed convergence rates in the $$\left\||{\cdot}\right\||$$ norm with respect to the polynomial degree $$r=2,3,4,5$$ $${\it{\Delta}} t$$ $$r = 2$$ $$r = 3$$ $$r = 4$$ $$r = 5$$ 1.00e$$-$$0 — — — — 5.00e$$-$$1 1.4017 2.2646 3.1802 4.1246 2.50e$$-$$1 1.5224 2.4577 3.4305 4.4403 1.25e$$-$$1 1.5295 2.5101 3.6584 4.5516 $${\it{\Delta}} t$$ $$r = 2$$ $$r = 3$$ $$r = 4$$ $$r = 5$$ 1.00e$$-$$0 — — — — 5.00e$$-$$1 1.4017 2.2646 3.1802 4.1246 2.50e$$-$$1 1.5224 2.4577 3.4305 4.4403 1.25e$$-$$1 1.5295 2.5101 3.6584 4.5516 Table 1 Scalar test case. Computed convergence rates in the $$\left\||{\cdot}\right\||$$ norm with respect to the polynomial degree $$r=2,3,4,5$$ $${\it{\Delta}} t$$ $$r = 2$$ $$r = 3$$ $$r = 4$$ $$r = 5$$ 1.00e$$-$$0 — — — — 5.00e$$-$$1 1.4017 2.2646 3.1802 4.1246 2.50e$$-$$1 1.5224 2.4577 3.4305 4.4403 1.25e$$-$$1 1.5295 2.5101 3.6584 4.5516 $${\it{\Delta}} t$$ $$r = 2$$ $$r = 3$$ $$r = 4$$ $$r = 5$$ 1.00e$$-$$0 — — — — 5.00e$$-$$1 1.4017 2.2646 3.1802 4.1246 2.50e$$-$$1 1.5224 2.4577 3.4305 4.4403 1.25e$$-$$1 1.5295 2.5101 3.6584 4.5516 Table 2 Scalar test case. Computed error in the $$\| \cdot \|_{L^2(0,T)}$$ and $$\| \cdot \|_{H^1(0,T)}$$ norms and corresponding convergence rates with respect to the polynomial approximation degree $$r=2,3,4,5$$ $$r$$ $${\it{\Delta}} t$$ $$\| \cdot \|_{L^2(0,T)}$$ rate $$\| \cdot \|_{H^1(0,T)}$$ rate 1.00e$$-$$0 4.4902e$$-$$02 — 6.7961e$$-$$01 — 8.00e$$-$$1 2.9612e$$-$$02 1.8655 4.9201e$$-$$01 1.4476 2 4.00e$$-$$1 7.0288e$$-$$03 2.0749 1.5416e$$-$$01 1.6743 2.00.e$$-$$1 1.2044e$$-$$03 2.5450 4.2254e$$-$$02 1.8672 1.00.e$$-$$1 1.7331e$$-$$04 2.7968 1.0988e$$-$$02 1.9431 1.00e$$-$$0 6.9895e$$-$$03 — 1.4649e$$-$$01 — 8.00e$$-$$1 3.9782e$$-$$03 2.5257 8.6484e$$-$$02 2.3617 3 4.00e$$-$$1 4.5984e$$-$$04 3.1129 1.3821e$$-$$02 2.6455 2.00.e$$-$$1 3.6876e$$-$$05 3.6404 1.8863e$$-$$03 2.8732 1.00.e$$-$$1 2.5546e$$-$$06 3.8515 2.4343e$$-$$04 2.9540 1.00e$$-$$0 1.0209e$$-$$03 — 2.4885e$$-$$02 — 8.00e$$-$$1 4.4989e$$-$$04 3.6723 1.1885e$$-$$02 3.3119 4 4.00e$$-$$1 2.3662e$$-$$05 4.2489 9.6064e$$-$$04 3.6290 2.00.e$$-$$1 9.0456e$$-$$07 4.7092 6.5449e$$-$$05 3.8755 1.00.e$$-$$1 3.0658e$$-$$08 4.8829 4.2099e$$-$$06 3.9585 1.00e$$-$$0 1.2456e$$-$$04 — 3.4801e$$-$$03 — 8.00e$$-$$1 4.2903e$$-$$05 4.7766 1.3399e$$-$$03 4.2773 5 4.00e$$-$$1 1.0725e$$-$$06 5.3220 5.4621e$$-$$05 4.6165 2.00.e$$-$$1 1.9994e$$-$$08 5.7452 1.8603e$$-$$06 4.8759 1.00.e$$-$$1 3.3483e$$-$$10 5.9000 5.9740e$$-$$08 4.9607 $$r$$ $${\it{\Delta}} t$$ $$\| \cdot \|_{L^2(0,T)}$$ rate $$\| \cdot \|_{H^1(0,T)}$$ rate 1.00e$$-$$0 4.4902e$$-$$02 — 6.7961e$$-$$01 — 8.00e$$-$$1 2.9612e$$-$$02 1.8655 4.9201e$$-$$01 1.4476 2 4.00e$$-$$1 7.0288e$$-$$03 2.0749 1.5416e$$-$$01 1.6743 2.00.e$$-$$1 1.2044e$$-$$03 2.5450 4.2254e$$-$$02 1.8672 1.00.e$$-$$1 1.7331e$$-$$04 2.7968 1.0988e$$-$$02 1.9431 1.00e$$-$$0 6.9895e$$-$$03 — 1.4649e$$-$$01 — 8.00e$$-$$1 3.9782e$$-$$03 2.5257 8.6484e$$-$$02 2.3617 3 4.00e$$-$$1 4.5984e$$-$$04 3.1129 1.3821e$$-$$02 2.6455 2.00.e$$-$$1 3.6876e$$-$$05 3.6404 1.8863e$$-$$03 2.8732 1.00.e$$-$$1 2.5546e$$-$$06 3.8515 2.4343e$$-$$04 2.9540 1.00e$$-$$0 1.0209e$$-$$03 — 2.4885e$$-$$02 — 8.00e$$-$$1 4.4989e$$-$$04 3.6723 1.1885e$$-$$02 3.3119 4 4.00e$$-$$1 2.3662e$$-$$05 4.2489 9.6064e$$-$$04 3.6290 2.00.e$$-$$1 9.0456e$$-$$07 4.7092 6.5449e$$-$$05 3.8755 1.00.e$$-$$1 3.0658e$$-$$08 4.8829 4.2099e$$-$$06 3.9585 1.00e$$-$$0 1.2456e$$-$$04 — 3.4801e$$-$$03 — 8.00e$$-$$1 4.2903e$$-$$05 4.7766 1.3399e$$-$$03 4.2773 5 4.00e$$-$$1 1.0725e$$-$$06 5.3220 5.4621e$$-$$05 4.6165 2.00.e$$-$$1 1.9994e$$-$$08 5.7452 1.8603e$$-$$06 4.8759 1.00.e$$-$$1 3.3483e$$-$$10 5.9000 5.9740e$$-$$08 4.9607 Table 2 Scalar test case. Computed error in the $$\| \cdot \|_{L^2(0,T)}$$ and $$\| \cdot \|_{H^1(0,T)}$$ norms and corresponding convergence rates with respect to the polynomial approximation degree $$r=2,3,4,5$$ $$r$$ $${\it{\Delta}} t$$ $$\| \cdot \|_{L^2(0,T)}$$ rate $$\| \cdot \|_{H^1(0,T)}$$ rate 1.00e$$-$$0 4.4902e$$-$$02 — 6.7961e$$-$$01 — 8.00e$$-$$1 2.9612e$$-$$02 1.8655 4.9201e$$-$$01 1.4476 2 4.00e$$-$$1 7.0288e$$-$$03 2.0749 1.5416e$$-$$01 1.6743 2.00.e$$-$$1 1.2044e$$-$$03 2.5450 4.2254e$$-$$02 1.8672 1.00.e$$-$$1 1.7331e$$-$$04 2.7968 1.0988e$$-$$02 1.9431 1.00e$$-$$0 6.9895e$$-$$03 — 1.4649e$$-$$01 — 8.00e$$-$$1 3.9782e$$-$$03 2.5257 8.6484e$$-$$02 2.3617 3 4.00e$$-$$1 4.5984e$$-$$04 3.1129 1.3821e$$-$$02 2.6455 2.00.e$$-$$1 3.6876e$$-$$05 3.6404 1.8863e$$-$$03 2.8732 1.00.e$$-$$1 2.5546e$$-$$06 3.8515 2.4343e$$-$$04 2.9540 1.00e$$-$$0 1.0209e$$-$$03 — 2.4885e$$-$$02 — 8.00e$$-$$1 4.4989e$$-$$04 3.6723 1.1885e$$-$$02 3.3119 4 4.00e$$-$$1 2.3662e$$-$$05 4.2489 9.6064e$$-$$04 3.6290 2.00.e$$-$$1 9.0456e$$-$$07 4.7092 6.5449e$$-$$05 3.8755 1.00.e$$-$$1 3.0658e$$-$$08 4.8829 4.2099e$$-$$06 3.9585 1.00e$$-$$0 1.2456e$$-$$04 — 3.4801e$$-$$03 — 8.00e$$-$$1 4.2903e$$-$$05 4.7766 1.3399e$$-$$03 4.2773 5 4.00e$$-$$1 1.0725e$$-$$06 5.3220 5.4621e$$-$$05 4.6165 2.00.e$$-$$1 1.9994e$$-$$08 5.7452 1.8603e$$-$$06 4.8759 1.00.e$$-$$1 3.3483e$$-$$10 5.9000 5.9740e$$-$$08 4.9607 $$r$$ $${\it{\Delta}} t$$ $$\| \cdot \|_{L^2(0,T)}$$ rate $$\| \cdot \|_{H^1(0,T)}$$ rate 1.00e$$-$$0 4.4902e$$-$$02 — 6.7961e$$-$$01 — 8.00e$$-$$1 2.9612e$$-$$02 1.8655 4.9201e$$-$$01 1.4476 2 4.00e$$-$$1 7.0288e$$-$$03 2.0749 1.5416e$$-$$01 1.6743 2.00.e$$-$$1 1.2044e$$-$$03 2.5450 4.2254e$$-$$02 1.8672 1.00.e$$-$$1 1.7331e$$-$$04 2.7968 1.0988e$$-$$02 1.9431 1.00e$$-$$0 6.9895e$$-$$03 — 1.4649e$$-$$01 — 8.00e$$-$$1 3.9782e$$-$$03 2.5257 8.6484e$$-$$02 2.3617 3 4.00e$$-$$1 4.5984e$$-$$04 3.1129 1.3821e$$-$$02 2.6455 2.00.e$$-$$1 3.6876e$$-$$05 3.6404 1.8863e$$-$$03 2.8732 1.00.e$$-$$1 2.5546e$$-$$06 3.8515 2.4343e$$-$$04 2.9540 1.00e$$-$$0 1.0209e$$-$$03 — 2.4885e$$-$$02 — 8.00e$$-$$1 4.4989e$$-$$04 3.6723 1.1885e$$-$$02 3.3119 4 4.00e$$-$$1 2.3662e$$-$$05 4.2489 9.6064e$$-$$04 3.6290 2.00.e$$-$$1 9.0456e$$-$$07 4.7092 6.5449e$$-$$05 3.8755 1.00.e$$-$$1 3.0658e$$-$$08 4.8829 4.2099e$$-$$06 3.9585 1.00e$$-$$0 1.2456e$$-$$04 — 3.4801e$$-$$03 — 8.00e$$-$$1 4.2903e$$-$$05 4.7766 1.3399e$$-$$03 4.2773 5 4.00e$$-$$1 1.0725e$$-$$06 5.3220 5.4621e$$-$$05 4.6165 2.00.e$$-$$1 1.9994e$$-$$08 5.7452 1.8603e$$-$$06 4.8759 1.00.e$$-$$1 3.3483e$$-$$10 5.9000 5.9740e$$-$$08 4.9607 5.2 Numerical results for the elastodynamics problem The numerical results presented in the sequel have been obtained by using the open source software SPEED (http://speed.mox.polimi.it; accessed 11 October 2017) suitably adapted to apply the DG scheme presented in Section 2 to the system (4.3). For all the numerical simulations, we consider the interval $$I=(0,T]$$ subdivided into $$N$$ time slabs $$I_n$$, for $$n=1,..., N$$ having uniform length $${\it{\Delta}} t$$. 5.2.1 Test case 1 We consider $${\it{\Omega}} = (0,1)^2, {\it{\Gamma}}_D = \partial {\it{\Omega}}$$ and $$T = 50$$. We set the mass density $$\rho=1$$, the Lamé coefficients $$\lambda=\mu=1$$, $$\zeta=0.01$$ and choose the data $$\mathbf{f}, \mathbf{u}_0, \mathbf{u}_1$$ such that the exact solution of problem in (4.1a)–(4.1e) is given by \begin{align*}\mathbf{u}= \sin(\sqrt{2}\pi t) \begin{bmatrix}- \sin^2(\pi x) \sin(2\pi y) \\ \sin(2\pi x) \sin^2(\pi y) \end{bmatrix}\!. \end{align*} In the first example, we compute the errors $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ versus $$1/ {\it{\Delta}} t$$ for $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,2,3,4$$ and varying the polynomial degree $$r = 2,3,4$$. As of the space discretization of the domain $${\it{\Omega}}$$, we consider a Cartesian grid with characteristic size $$h = 0.125$$ and we set a polynomial approximation degree $$q=r+1$$. The computed errors are reported in Fig. 4. Fig. 4. View largeDownload slide Test case 1. Computed errors $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ as a function of $$1/ {\it{\Delta}} t$$ for $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,\dots, 4$$ (log-log scale), $$r=2,3,4$$, $$h=0.125$$ and $$q=r+1$$. Fig. 4. View largeDownload slide Test case 1. Computed errors $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ as a function of $$1/ {\it{\Delta}} t$$ for $${\it{\Delta}} t = 2^{-\ell}, \ell = 0,1,\dots, 4$$ (log-log scale), $$r=2,3,4$$, $$h=0.125$$ and $$q=r+1$$. A convergence rate of $${\it{\Delta}} t^{r - \frac{1}{2}}$$ is observed, in accordance with the theoretical estimate (3.31). In Table 3, the computed convergence rates are reported for $$r = 2,3,4$$, respectively. Table 3 Test case 1. Computed errors $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ and computed convergence rates, $$r = 2,3,4$$, $$h=0.125$$ and $$q=r+1$$ $$r$$ $${\it{\Delta}} t$$ $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ rate 1.00e$$-$$0 4.1632e+02 — 5.00e$$-$$1 2.5100e+02 0.7351 2 2.50e$$-$$1 9.1706e+01 1.4474 1.25e$$-$$1 3.3960e+01 1.4331 6.25e$$-$$2 1.2156e+01 1.4821 1.00e$$-$$0 3.3909e+02 — 5.00e$$-$$1 9.2986e+01 1.8665 3 2.50e$$-$$1 1.4534e+01 2.6775 1.25e$$-$$1 2.5484e+00 2.5117 6.25e$$-$$2 4.4983e$$-$$01 2.5021 1.00e$$-$$0 1.7465e+02 — 5.00e$$-$$1 1.4155e+01 3.6250 4 2.50e$$-$$1 1.3873e+00 3.3509 1.25e$$-$$1 1.2308e$$-$$01 3.4946 6.25e$$-$$2 1.0550e$$-$$02 3.5441 $$r$$ $${\it{\Delta}} t$$ $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ rate 1.00e$$-$$0 4.1632e+02 — 5.00e$$-$$1 2.5100e+02 0.7351 2 2.50e$$-$$1 9.1706e+01 1.4474 1.25e$$-$$1 3.3960e+01 1.4331 6.25e$$-$$2 1.2156e+01 1.4821 1.00e$$-$$0 3.3909e+02 — 5.00e$$-$$1 9.2986e+01 1.8665 3 2.50e$$-$$1 1.4534e+01 2.6775 1.25e$$-$$1 2.5484e+00 2.5117 6.25e$$-$$2 4.4983e$$-$$01 2.5021 1.00e$$-$$0 1.7465e+02 — 5.00e$$-$$1 1.4155e+01 3.6250 4 2.50e$$-$$1 1.3873e+00 3.3509 1.25e$$-$$1 1.2308e$$-$$01 3.4946 6.25e$$-$$2 1.0550e$$-$$02 3.5441 Table 3 Test case 1. Computed errors $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ and computed convergence rates, $$r = 2,3,4$$, $$h=0.125$$ and $$q=r+1$$ $$r$$ $${\it{\Delta}} t$$ $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ rate 1.00e$$-$$0 4.1632e+02 — 5.00e$$-$$1 2.5100e+02 0.7351 2 2.50e$$-$$1 9.1706e+01 1.4474 1.25e$$-$$1 3.3960e+01 1.4331 6.25e$$-$$2 1.2156e+01 1.4821 1.00e$$-$$0 3.3909e+02 — 5.00e$$-$$1 9.2986e+01 1.8665 3 2.50e$$-$$1 1.4534e+01 2.6775 1.25e$$-$$1 2.5484e+00 2.5117 6.25e$$-$$2 4.4983e$$-$$01 2.5021 1.00e$$-$$0 1.7465e+02 — 5.00e$$-$$1 1.4155e+01 3.6250 4 2.50e$$-$$1 1.3873e+00 3.3509 1.25e$$-$$1 1.2308e$$-$$01 3.4946 6.25e$$-$$2 1.0550e$$-$$02 3.5441 $$r$$ $${\it{\Delta}} t$$ $$\left\||{\mathbf{u}-{\mathbf{u}_{\textrm{DG}}}}\right\||$$ rate 1.00e$$-$$0 4.1632e+02 — 5.00e$$-$$1 2.5100e+02 0.7351 2 2.50e$$-$$1 9.1706e+01 1.4474 1.25e$$-$$1 3.3960e+01 1.4331 6.25e$$-$$2 1.2156e+01 1.4821 1.00e$$-$$0 3.3909e+02 — 5.00e$$-$$1 9.2986e+01 1.8665 3 2.50e$$-$$1 1.4534e+01 2.6775 1.25e$$-$$1 2.5484e+00 2.5117 6.25e$$-$$2 4.4983e$$-$$01 2.5021 1.00e$$-$$0 1.7465e+02 — 5.00e$$-$$1 1.4155e+01 3.6250 4 2.50e$$-$$1 1.3873e+00 3.3509 1.25e$$-$$1 1.2308e$$-$$01 3.4946 6.25e$$-$$2 1.0550e$$-$$02 3.5441 5.2.2 Test case 2: a numerical test with nonreflective boundary conditions When simulating seismic wave propagations, ideally artificial boundaries should consent incident wave to be propagated without producing any reflection. A simple strategy consists in imposing nonreflective conditions on artificial boundaries $${\it{\Gamma}}_{NR} \subset \partial{\it{\Omega}}$$ such that $${\it{\Gamma}}_{NR} \cap {\it{\Gamma}}_D \cap {\it{\Gamma}}_N = \emptyset$$, $$\begin{cases}\dfrac{\partial ({\bf{u}} \cdot {\bf{n}})}{\partial n} = -\dfrac{1}{V_P} \dfrac{\partial ({\bf{u}} \cdot {\bf{n}})}{\partial t} + \dfrac{V_S - V_P}{V_P} \dfrac{\partial ({\bf{u}} \cdot {\bf{\tau}})}{\partial t}, \\ \dfrac{\partial ({\bf{u}} \cdot {\bf{\tau}})}{\partial n} = -\dfrac{1}{V_S} \dfrac{\partial ({\bf{u}} \cdot {\bf{\tau}})}{\partial t} + \dfrac{V_S - V_P}{V_S} \dfrac{\partial ({\bf{u}} \cdot {\bf{n}})}{\partial t}, \end{cases}$$ (5.2) where $$\mathbf{n} = (n_1, n_2)$$ (respectively $${\bf{\tau}} = (\tau_1, \tau_2)$$) is the unit normal (respectively tangential) vector to $$\partial {\it{\Omega}}$$. Here, $$V_{\rm P} = \sqrt{(\lambda+2\mu)/\rho}$$ and $$V_{\rm S} = \sqrt{\mu/\rho}$$ are the propagation velocities of compressional (P) and shear (S) waves, respectively. Equations (5.2) are first-order nonreflecting boundary conditions (see, e.g., Stacey, 1988). Their use breaks of symmetry in the matrix $$K$$ in (4.3). The test consists of propagating a pure shear plane wave through a viscoelastic, horizontally layered soil profile (cf. Fig. 5). The mechanical properties of the layers are summarized in Table 4. The incidence is orthogonal to the free surface and the excitation consists of a displacement Ricker wavelet $$\textbf{f} (\textbf{x}, t) = h(t) \delta_q(\textbf{x} - \textbf{x}_0)$$ with \begin{equation*}h(t) = h_0 \left(1-2 \left(\pi f_{\rm peak}\right)^2 \left(t-t_0\right)^2 \right)e^{-(\pi f_{\rm peak})^2 (t-t_0)^ 2}, \end{equation*} where $$f_{\rm peak} = 1$$ Hz, $$t_0 = 2$$ s and $$h_0=1$$ is the amplitude of the wave in time domain. Here, $$\delta_q$$ is the numerical delta function, i.e., a polynomial of degree $$q$$ that approximates the $$\delta$$ distribution. In this case, the forcing term is applied at points $$\textbf{x}_0$$ lying at the bottom of the domain (see Fig. 5). Table 4 Mechanical properties Material $$\rho$$ $$\lambda$$ $$\mu$$ $$\zeta$$ 1 2000 2.00e+07 2.00e+07 3.1416e$$-$$02 2 2000 5.00e+08 5.00e+08 3.1416e$$-$$03 Material $$\rho$$ $$\lambda$$ $$\mu$$ $$\zeta$$ 1 2000 2.00e+07 2.00e+07 3.1416e$$-$$02 2 2000 5.00e+08 5.00e+08 3.1416e$$-$$03 Table 4 Mechanical properties Material $$\rho$$ $$\lambda$$ $$\mu$$ $$\zeta$$ 1 2000 2.00e+07 2.00e+07 3.1416e$$-$$02 2 2000 5.00e+08 5.00e+08 3.1416e$$-$$03 Material $$\rho$$ $$\lambda$$ $$\mu$$ $$\zeta$$ 1 2000 2.00e+07 2.00e+07 3.1416e$$-$$02 2 2000 5.00e+08 5.00e+08 3.1416e$$-$$03 Fig. 5. View largeDownload slide Domain $${\it{\Omega}}$$ for the test case 2. Nonreflective boundary conditions are imposed at the bottom edge, $${\bf{\sigma}}{\bf{n}} = \textbf{0}$$ is imposed at the top edge (free-surface condition), while homogeneous Dirichlet boundary conditions (in y-direction) are imposed on the lateral sides. Fig. 5. View largeDownload slide Domain $${\it{\Omega}}$$ for the test case 2. Nonreflective boundary conditions are imposed at the bottom edge, $${\bf{\sigma}}{\bf{n}} = \textbf{0}$$ is imposed at the top edge (free-surface condition), while homogeneous Dirichlet boundary conditions (in y-direction) are imposed on the lateral sides. The plane wave rises from the bottom of $${\it{\Omega}}$$, reaches the top of the computational domain, increases its amplitude (because of the free surface condition) and then is propagated backward, completely absorbed from the bottom boundary, thanks to the absorbing condition (5.2). To prevent spurious oscillations inside the domain, homogeneous Dirichlet boundary conditions (in y-direction) are imposed on the lateral sides of $${\it{\Omega}}$$ (cf. Fig. 5). For the simulation, we fix $$T=20 s$$, $${\it{\Delta}} t = 0.01 s, h=10$$ and $$q=4$$ for the space discretization, while we use second-order polynomials ($$r=2$$) for the time integration. We compare the results obtained using the proposed DG method with the analogous ones obtained coupling the DGSE space discretization with the classical leap-frog method (cf. Quarteroni & Valli, 2008, e.g., with $${\it{\Delta}} t = 0.0001$$ s. Notice that this is the biggest time step allowed by the CFL condition (Courant et al., 1928). We plot the displacement time histories (in the x-direction) recorded by two monitors: $$R_1$$ set on the free surface and $$R_2$$ located across Material 1 and Material 2 (cf. Fig. 5). In Fig. 6, we show the results obtained with the leap-frog and the DG method along with the semianalytical solution $${\bf{u}}_{TH}$$ based on the Thomson–Haskell propagation matrix method (see, e.g., Haskell, 1953). We can clearly see that all the methods produce almost identical solutions. To quantify the distance between the curves, in Fig. 6, we plot the error $$| {\bf{u}}_{TH}(t)- {\bf{u}}_{*}(t)|$$ for $$t\in(0,T]$$, where $${\bf{u}}_{*}$$ is either the solution obtained with the leap-frog scheme or the one obtained with the DG method and $${\bf{u}}_{TH}$$ is the Thomson–Haskell semianalytical solution. Both methods achieve the same level of accuracy. Finally, we compare the computational time of our method with the analogous results obtained with the leap-frog time integration scheme. In Table 5, we show the computational costs for each time step (second column) as well as the overall computational time for the time loop (third column). As expected, each time step is more expensive in our approach, as the solution of a linear system is involved, but overall this can be balanced using larger times steps, as no stability constraints are required. We point out that for our DG method, the matrix is assembled and factorized once and for all before the time loop. A more thorough comparison is under investigation, which also includes suitable parallelization techniques and possible preconditioned iterative methods. Fig. 6. View largeDownload slide Displacement registered by $$R_1$$ (a) and $$R_2$$ (b). Fig. 6. View largeDownload slide Displacement registered by $$R_1$$ (a) and $$R_2$$ (b). Table 5 Computational times for the time loop (space polynomial degree $$q = 4$$) Single time step Time loop $${\it{\Delta}} t$$ $$\| {\bf{u}}_{TH}- {\bf{u}}_{*}\|_{L^2(0,T)}$$ leap-frog 4.57e$$-$$3 s 9.15 s 1.e$$-$$4 s 1.70e$$-$$3 ($$R_1$$), 6.94e$$-$$4 ($$R_2$$) DG 2.98e$$-$$2 s 7.45 s 8.e$$-$$2 s 1.99e$$-$$3 ($$R_1$$), 8.01e$$-$$4 ($$R_2$$) Single time step Time loop $${\it{\Delta}} t$$ $$\| {\bf{u}}_{TH}- {\bf{u}}_{*}\|_{L^2(0,T)}$$ leap-frog 4.57e$$-$$3 s 9.15 s 1.e$$-$$4 s 1.70e$$-$$3 ($$R_1$$), 6.94e$$-$$4 ($$R_2$$) DG 2.98e$$-$$2 s 7.45 s 8.e$$-$$2 s 1.99e$$-$$3 ($$R_1$$), 8.01e$$-$$4 ($$R_2$$) Table 5 Computational times for the time loop (space polynomial degree $$q = 4$$) Single time step Time loop $${\it{\Delta}} t$$ $$\| {\bf{u}}_{TH}- {\bf{u}}_{*}\|_{L^2(0,T)}$$ leap-frog 4.57e$$-$$3 s 9.15 s 1.e$$-$$4 s 1.70e$$-$$3 ($$R_1$$), 6.94e$$-$$4 ($$R_2$$) DG 2.98e$$-$$2 s 7.45 s 8.e$$-$$2 s 1.99e$$-$$3 ($$R_1$$), 8.01e$$-$$4 ($$R_2$$) Single time step Time loop $${\it{\Delta}} t$$ $$\| {\bf{u}}_{TH}- {\bf{u}}_{*}\|_{L^2(0,T)}$$ leap-frog 4.57e$$-$$3 s 9.15 s 1.e$$-$$4 s 1.70e$$-$$3 ($$R_1$$), 6.94e$$-$$4 ($$R_2$$) DG 2.98e$$-$$2 s 7.45 s 8.e$$-$$2 s 1.99e$$-$$3 ($$R_1$$), 8.01e$$-$$4 ($$R_2$$) 6. Conclusions In this work, we have developed a new high-order DG finite element method for the temporal discretization of Cauchy problems for second-order differential equations. To show the capabilities of our scheme, we have applied it for the time integration of second-order system of equations resulting after DG semidiscretization (in space) of the elastodynamics equation. Our formulation contains suitable stabilization terms, which allowed the construction of an appropriate energy norm that naturally arose by the variational formulation of the problem. We have studied the well-posedness of the resulting scheme and proved a priori error estimates that properly depend on the local polynomial approximation degree and the local regularity of the exact solution. Our theoretical results have been confirmed by the numerical experiments carried out for both simplified test cases and on examples of practical interest. The method presented is well suited for the adaptive choice the discretization parameters employed within each time slab, namely the local time step $${\it{\Delta}} t_n$$ and the local approximation degree $$r_n$$, $$n= 1, \dots, N$$. The development of suitable a posteriori error estimators, to be employed within (space)-time adaptive strategies, is under investigation and will make the subject of future research. Acknowledgements The authors are grateful to the referees for their thorough and constructive comments that have greatly contributed to the improvement of the article. Funding Scientific Independence of young Researchers (SIR; RBSI14VT0S in part to P.F. Antonietti and I. Mazzieri); ‘PolyPDEs: Non-conforming polyhedral finite element methods for the approximation of partial differential equations’ by Italian Ministry of Education, Universities and Research (MIUR). References Adams R. A. & Fournier J. J. F. ( 2003 ) Sobolev Spaces . Pure and Applied Mathematics , vol. 140, 2nd edn. Amsterdam : Elsevier , pp. xiv+305 . Adjerid S. & Temimi H. ( 2011 ) A discontinuous Galerkin method for the wave equation. Comput. Methods Appl. Mech. Engrg. , 200 , 837 – 849 . Google Scholar Crossref Search ADS Antonietti P. F. , Ayuso de Dios B. , Mazzieri I. & Quarteroni A. ( 2016a ) Stability analysis of discontinuous Galerkin approximations to the elastodynamics problem. J. Sci. Comput. , 68 , 143 – 170 . Google Scholar Crossref Search ADS Antonietti P. F. , Ferroni A. , Mazzieri I. , Paolucci R. , Quarteroni A. , Smerzini C. & Stupazzini M. ( 2016b ) A review on numerical modelling of earthquakes through Discontinuous Galerkin Spectral Element methods ( submitted to ESAIM: Proceedings and Surveys ). Antonietti P. F. , Marcati C. , Mazzieri I. & Quarteroni A. ( 2016c ) High order discontinuous Galerkin methods on simplicial elements for the elastodynamics equation. Numer. Algorithms , 71 , 181 – 206 . Google Scholar Crossref Search ADS Antonietti P. F. , Ferroni A. , Mazzieri I. & Quarteroni A. ( 2017 ) hp-version discontinuous Galerkin approximations of the elastodynamics equation. Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2016 ( Bittencourt M. et al. eds). ecture Notes in Computational Science and Engineering 119. https://doi.org/10.1007/978-3-319-65870-4_1 . Antonietti P. F. , Mazzieri I. , Quarteroni A. & Rapetti F. ( 2012 ) Non-conforming high order approximations of the elastodynamics equation. Comput. Methods Appl. Mech. Engrg. , 209–212 , 212 – 238 . Google Scholar Crossref Search ADS Arnold D. N. ( 1982 ) An interior penalty finite element method with discontinuous elements. SIAM J. Numer. Anal. , 19 , 742 – 760 . Google Scholar Crossref Search ADS Arnold D. N. , Brezzi F. , Cockburn B. & Marini L. D. ( 2001–2002 ) Unified analysis of discontinuous Galerkin methods for elliptic problems. SIAM J. Numer. Anal. , 39 , 1749 – 1779 . Google Scholar Crossref Search ADS Babuska I. & Suri M. ( 1987 ) The h-p version of the finite element method with quasiuniform meshes. ESAIM Math. Model. Numer. Anal. , 21 , 199 – 238 . Google Scholar Crossref Search ADS Baccouch M. ( 2012 ) A local discontinuous Galerkin method for the second-order wave equation. Comput. Methods Appl. Mech. Engrg. , 209–212 , 129 – 143 . Google Scholar Crossref Search ADS Banks H. T. , Birch M. J. , Brewin M. P. , Greenwald S. E. , Hu S. , Kenz Z. R. , Kruse C. , Maischak M. , Shaw S. & Whiteman J. R. ( 2014 ) High-order space-time finite element schemes for acoustic and viscodynamic wave equations with temporal decoupling. Int. J. Numer. Methods Engrg ., 98 , 131 – 156 . Google Scholar Crossref Search ADS Butcher J. C. ( 2008 ) Numerical Methods for Ordinary Differential Equations . New York : John Wiley & Sons, Ltd, pp. 137 – 316 . Canuto C. , Hussaini M. , Quarteroni A. & Zang T. ( 2006 ) Spectral Methods (Fundamental in Single Domains ). Berlin : Springer. Canuto C. G. , Hussaini M. Y. , Quarteroni A. & Zang T. A. ( 2007 ) Spectral Methods: Evolution to Complex Geometries and Applications to Fluid Dynamics . Berlin, Heidelberg-New York : Springer . Cockburn B. & Shu C.-W. ( 1998 ) The local discontinuous Galerkin method for time-dependent convection-diffusion systems. SIAM J. Numer. Anal. , 35 , 2440 – 2463 . Google Scholar Crossref Search ADS Cohen M. A. & Tan C. O. ( 2012 ) A polynomial approximation for arbitrary functions. Appl. Math. Lett. , 25 , 1947 – 1952 . Google Scholar Crossref Search ADS Collino F. , Fouquet T. & Joly P. ( 2003 ) A conservative space-time mesh refinement method for the 1-D wave equation. Part I: construction. Numer. Math. , 95 , 197 – 221 . Google Scholar Crossref Search ADS Courant R. , Friedrichs K. & Lewy H. ( 1928 ) Über die partiellen Differenzengleichungen der mathematischen Physik. Math. Ann. , 100 , 32 – 74 . Google Scholar Crossref Search ADS Delcourte S. & Glinsky N. ( 2015 ) Analysis of a high-order space and time discontinuous Galerkin method for elastodynamic equations application to 3D wave propagation. ESAIM Math. Model. Numer. Anal. , 49 , 1085 – 1126 . Google Scholar Crossref Search ADS Delfour M. , Hager W. & Trochu F. ( 1981 ) Discontinuous Galerkin methods for ordinary differential equations. Math.Comput. , 36 , 455 – 473 . Google Scholar Crossref Search ADS Di Pietro D. A. & Ern A. ( 2012 ) Mathematical Aspects of Discontinuous Galerkin Methods. Mathématiques et Applications , vol. 69. Berlin : Springer . Diaz J. & Grote M. J. ( 2009 ) Energy conserving explicit local time stepping for second-order wave equations. SIAM J. Sci. Comput. , 31 , 1985 – 2014 . Google Scholar Crossref Search ADS Dumbser M. , Käser M. & Toro E. ( 2007 ) An arbitrary high-order Discontinuous Galerkin method for elastic waves on unstructured meshes - V. Local time stepping and p-adaptivity. Geophys. J. Int. , 171 , 695 – 717 . Google Scholar Crossref Search ADS Ferroni A. , Antonietti P. F. , Mazzieri I. & Quarteroni A. ( 2017 ) Dispersion-dissipation analysis of 3D continuous and discontinuous spectral element methods for the elastodynamics equation. Geophys. J. Int. Available at https://doi.org/10.1093/gji/ggx384. Grote M. J. & Mitkova T. ( 2013 ) High-order explicit local time-stepping methods for damped wave equations. J. Comput. Appl. Math. , 239 , 270 – 289 . Google Scholar Crossref Search ADS Grote M. J. , Schneebeli A. & Schötzau D. ( 2006 ) Discontinuous Galerkin finite element method for the wave equation. SIAM J. Numer. Anal. , 44 , 2408 – 2431 . Google Scholar Crossref Search ADS Haskell N. A. ( 1953 ) The dispersion of surface waves on multilayered media. Bull. Seismol. Soc. Am. , 43 , 17 – 43 . Hesthaven J. & Warburton T. ( 2008 ) Nodal Discontinuous Galerkin Methods . Texts in Applied Mathematics , vol. 54. Berlin : Springer , pp. xiv+501 . Houston P. , Schwab C. & Süli E. ( 2000 ) Stabilized hp-finite element methods for first-order hyperbolic problems. SIAM J. Numer. Anal. , 37 , 1618 – 1643 . Google Scholar Crossref Search ADS Hughes T. J. & Hulbert G. M. ( 1988 ) Space-time finite element methods for elastodynamics: formulations and error estimates. Comput. Methods Appl. Mech. Engrg. , 66 , 339 – 363 . Google Scholar Crossref Search ADS Hughes T. J. & Stewart J. R. ( 1996 ) A space-time formulation for multiscale phenomena. J. Comput. Appl. Math. , 74 , 217 – 229 . Google Scholar Crossref Search ADS Hulbert G. M. & Hughes T. J. ( 1990 ) Space-time finite element methods for second-order hyperbolic equations. Comput. Methods Appl. Mech. Engrg. , 84 , 327 – 348 . Google Scholar Crossref Search ADS Johnson C. ( 1993 ) Discontinuous Galerkin finite element methods for second order hyperbolic problems. Comput. Methods Appl. Mech. Engrg. , 107 , 117 – 129 . Google Scholar Crossref Search ADS Kroopnick A. ( 1999 ) Bounded and l2-solutions to a second order nonlinear differential equation with a square integrable forcing term. Int. J. Math. Sci. , 22 , 569 – 571 . Google Scholar Crossref Search ADS Lesaint P. & Raviart P. ( 1974 ) On a Finite Element Method for Solving the Neutron Transport Equation . Publications mathématiques et informatique de Rennes. S4 : 1 – 40 . LeVeque R. ( 2007 ) Finite Difference Methods for Ordinary and Partial Differential Equations . Philadelphia, PA, USA : SIAM - Society for Industrial and Applied Mathematics. Mazzieri I. , Stupazzini M. , Guidotti R. & Smerzini C. ( 2013 ) Speed: spectral elements in elastodynamics with discontinuous Galerkin: a non-conforming approach for 3D multi-scale problems. Int. J. Numer. Methods Engrg. , 95 , 991 – 1010 . Google Scholar Crossref Search ADS Quarteroni A. ( 2014 ) Numerical Models for Differential Problems . MS&A. Modeling, Simulation and Applications, vol. 8. Italia, Milan : Springer. Quarteroni A. , Sacco R. & Saleri F. ( 2007 ) Numerical Mathematics . Texts in Applied Mathematics , vol. 37 , 2nd edn. Berlin : Springer , pp. xviii+655 . Quarteroni A. & Valli A. ( 2008 ) Numerical Approximation of Partial Differential Equations . Springer Series in Computational Mathematics. Berlin, Heidelberg : Springer. Reed W. H. & Hill T. R. ( 1973 ) Triangular mesh methods for the neutron transport equation. Technical Report LA-UR-73-479 . Los Alamos, New Mexico : Los Alamos Scientific Laboratory. Rivière B. ( 2008 ) Discontinuous Galerkin Methods for Solving Elliptic and Parabolic Equations . Philadelphia, PA, USA : Society for Industrial and Applied Mathematics. Rivière B. & Wheeler M. F. ( 2003 ) Discontinuous finite element methods for acoustic and elastic wave problems. Contemporary Mathematics , 329 , 271 – 282 . Google Scholar Crossref Search ADS Schötzau D. & Schwab C. ( 2000 ) Time discretization of parabolic problems by the hp-version of the discontinuous Galerkin finite element method. SIAM J. Numer. Anal. , 38 , 837 – 875 . Google Scholar Crossref Search ADS Schwab C. ( 1998 ) p- and hp- Finite Element Methods. Theory and Applications in Solid and Fluid Mechanics . New York : Oxford University Press. Stacey R. ( 1988 ) Improved transparent boundary formulations for the elastic-wave equation. Bull. Seismol. Soc. Am. , 78 , 2089 – 2097 . Taube A. , Dumbser M. , Munz C.-D. & Schneider R. ( 2009 ) A high-order discontinuous Galerkin method with time-accurate local time stepping for the Maxwell equations. Int. J. Numer. Model. El. , 22 , 77 – 103 . Google Scholar Crossref Search ADS Thompson L. L. & Pinsky P. M. ( 1996 ) A space-time finite element method for structural acoustics in infinite domains part 1: Formulation, stability and convergence. Comput. Methods Appl. Mech. Engrg. , 132 , 195 – 227 . Google Scholar Crossref Search ADS van der Vegt J. J. W. , Klaij C. M. , van der Bos F. & van der Ven H. ( 2006 ) Space-time discontinuous Galerkin method for the compressible Navier-Stokes equations on deforming meshes. In Wesseling P. Onate E. & Periaux J. (eds), Proceedings European Conference on Computational Fluid Dynamics ECCOMAS CFD 2006 . Delft : TU Delft. Walkington N. J. ( 2014 ) Combined dg–cg time stepping for wave equations. SIAM J. Numer. Anal. , 52 , 1398 – 1417 . Google Scholar Crossref Search ADS Werder T. , Gerdes K. , Schötzau D. & Schwab C. ( 2001 ) hp-discontinuous Galerkin time stepping for parabolic problems. Comput. Methods Appl. Mech. Engrg. , 190 , 6685 – 6708 . Google Scholar Crossref Search ADS Wheeler M. F. ( 1978 ) An elliptic collocation-finite element method with interior penalties. SIAM J. Numer. Anal. , 15 , 152 – 161 . Google Scholar Crossref Search ADS Wilcox L. C. , Stadler G. , Burstedde C. & Ghattas O. ( 2010 ) A high-order discontinuous Galerkin method for wave propagation through coupled elastic-acoustic media. J. Comput. Phys. , 229 , 9373 – 9396 . Google Scholar Crossref Search ADS Yang Y. , Chirputkar S. , Alpert D. N. , Eason T. , Spottswood S. & Qian D. ( 2012 ) Enriched space time finite element method: a new paradigm for multiscaling from elastodynamics to molecular dynamics. Int. J. Numer. Methods Engrg. , 92 , 115 – 140 . Google Scholar Crossref Search ADS © The authors 2017. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

### Journal

IMA Journal of Numerical AnalysisOxford University Press

Published: Oct 16, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations