Explicit Local Time-Stepping for the Inhomogeneous Wave Equation with Optimal Convergence

Marcus J. Grote Department of Mathematics and Computer Science, University of Basel, Spiegelgasse 1, 4051 Basel, Switzerland, ([email protected]) Simon R. J. Michel Mathematisches Institut, Universität Bern, Sidlerstrasse 5, 3012 Bern, Switzerland, ([email protected]) Stefan A. Sauter Institute for Mathematics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland, ([email protected])

Abstract

Adaptivity and local mesh refinement are crucial for the efficient numerical simulation of wave phenomena in complex geometry. Local mesh refinement, however, can impose a tiny time-step across the entire computational domain when using explicit time integration. By taking smaller time-steps yet only inside locally refined regions, local time-stepping methods overcome the stringent CFL stability restriction imposed on the global time-step by a small fraction of the elements without sacrificing explicitness. In [21], a leapfrog based local time-stepping method was proposed for the inhomogeneous wave equation, which applies standard leapfrog time-marching with a smaller time-step inside the refined region. Here, to remove potential instability at certain time-steps, a stabilized version is proposed which leads to optimal $L^{2}$ -error estimates under a CFL condition independent of the coarse-to-fine mesh ratio. Moreover, a weighted transition is introduced to restore optimal $H^{1}$ -convergence when the source is nonzero across the coarse-to-fine mesh interface. Numerical experiments corroborate the theoretical error estimates and illustrate the usefulness of these improvements.

1 Introduction

For the simulation of wave phenomena, finite element methods, be it conforming or discontinuous, easily accommodate complex geometry or material interfaces. Local mesh refinement, however, seriously hampers the efficiency of explicit time integrators, such as the popular leapfrog (LF) method, when a small fraction of elements impose a tiny time-step for stability across the entire mesh. By taking smaller time-steps, but only inside those refined regions, while taking larger time-steps elsewhere, local time-stepping methods (LTS) permit to overcome that overly stringent CFL stability constraint on the global time-step without sacrificing the explicitness or parallel efficiency of explicit time integrators.

Following the classical IMEX (implicit-explicit) approach, first locally implicit methods that combine the (explicit) second-order LF (or Størmer-Verlet) method in the coarse part with the (implicit) Crank-Nicolson method in the fine part of the mesh were proposed and further analyzed in [30, 16, 33, 12, 23, 24]. In [5], fourth-order energy-preserving IMEX schemes were derived for the wave equation, when the computational domain is divided into a coarse and fine region using Lagrange multipliers along the interface.

Even earlier, an explicit local time-stepping (LTS) method was proposed for the wave equation, which also employs the explicit LF method in the fine part of the mesh, yet with a smaller step-size [9, 11]; error estimates were derived in [10, 26] in the one-dimensional case. Stability of this second-order method, however, is guaranteed by enforcing energy conservation which requires at every time-step the solution of a linear system for the shared unknowns at the coarse/fine intersection; hence, the overall scheme is not fully explicit.

Fully explicit second-order LTS methods based on Størmer-Verlet in time and DG discretizations in space were proposed for Maxwell’s equations in [30, 29]. Starting from the popular LF method, fully explicit local time-stepping methods of arbitrary order were proposed for the homogeneous wave equation in [13] and for the inhomogeneous wave equation with nonzero forcing in [21]. Inside the locally refined region, they apply standard leapfrog time-marching with a smaller time-step, which also implies evaluating the inhomogeneous right-hand side at all intermediate times. Hence, when combined with a mass-lumped conforming [8, 7] or discontinuous Galerkin FE discretization [22] in space, the resulting LF-LTS method is truly explicit, inherently parallel and also conserves a discrete energy in the absence of forcing; it was successfully applied to 3D seismic wave propagation [28]. A multilevel version was later proposed [14] and achieved high parallel efficiency [31].

In [19], optimal $L^{2}$ -error estimates for the LF-LTS method from [13] were derived for $H^{1}$ -conforming FEM discretizations (with mass-lumping), though under a CFL condition where the global time-step $\Delta t$ in fact depended on the smallest elements in the mesh. For the convergence analysis, the inner loop over the $p$ local LF steps of size $\Delta t/p$ was rewritten in terms of a single global time-step $\Delta t$ , which involves Chebyshev polynomials of the first kind. In general one cannot improve upon the stability constraint on $\Delta t$ , as the LF-LTS method may become unstable at certain values of $\Delta t$ .

By replacing those Chebyshev polynomials by “damped” Chebyshev polynomials with damping parameter $\nu>0$ , which slightly modifies the LF-LTS scheme inside the refined region, a stabilized LF-LTS( $\nu$ ) version was proposed in [20] for the homogeneous wave equation, and also independently in [2], which led to optimal $L^{2}$ -error estimates yet under a CFL condition independent of the coarse-to-fine aspect ratio. The same damped Chebyshev polynomials were also previously used to stabilize leapfrog-Chebyshev (LFC, discrete Gautschi-type) methods [18, 4], where the “non-stiff” (or slowly varying), possibly nonlinear source term is evaluated only at every global time-step throughout the entire computational domain. In fact, the error analysis for semilinear ODE’s in [2] and for DG discretizations of linear wave equations in [3] both apply to the inhomogeneous case, too, when the source is evaluated only once during each global time-step, akin to the approach in [4].

Main contributions. The leapfrog based local time-stepping (LF-LTS) method introduced in [21] for the inhomogeneous wave equation can become unstable for certain time-steps. To remove those potential instabilities, we propose a stabilized version which leads to optimal $L^{2}$ -error estimates under a CFL condition independent of the coarse-to-fine mesh ratio. Thus we effectively introduce stabilization into the original method from [21], which, in contrast to the “split-LFC” method from [3], corresponds to standard LF time-marching, yet with a smaller time-step, inside the refined region. Generalizing the proof of $L^{2}$ -convergence for the homogeneous wave equation [20] to the inhomogeneous case is not straightforward due to the source evaluations at all intermediate local time-steps. In fact it requires the development of new analytical techniques involving Chebyshev polynomials of the second kind to rewrite and analyze the effect of those local time-steps and thus finally prove the stability and optimal consistency order of the method.

In [6] it was shown that the expected optimal $H^{1}$ -convergence rate of the LF-LTS method [21] for the inhomogeneous wave equation is reduced by $h^{1/2}$ when the source is nonzero across the coarse-to-fine mesh interface. To overcome this difficulty and thereby restore optimal $H^{1}$ -convergence, we introduce a new weighted transition which relies on a discrete mesh-based distance function.

Outline. The rest of our article is structured as follows. In Section 2, we first recall the $H^{1}$ -conforming Galerkin formulation of the inhomogeneous wave equation together with needed notation as the mesh is split into a coarse and a locally refined part. In Section 2.3, this leads to the stabilized LF-LTS $(\nu)$ algorithm for the inhomogenous wave equation, which also evaluates the source term inside the locally refined region at all intermediate time-steps. In Section 3, we prove optimal $L^{2}$ -convergence rates under a CFL condition independent of the fine-to-coarse mesh ratio. To accomodate for the intermediate source evaluations, our error analysis requires both Chebyshev polynomials of the first and the second kind. Still our analysis also includes the somewhat simpler situation when the source is evaluated at global time-steps only as in [3]. To achieve optimal $H^{1}$ -convergence regardless of the source, we introduce in Section 4 the new weighted transition which relies on a discrete mesh-based distance function. In Section 5, we present numerical experiments which illustrate the optimal $L^{2}$ -convergence rates expected from theory. In doing so, we also compare the LF-LTS method with the alternative “split-LFC” approach from [3], where the source inside the locally refined region is only evaluated once at each global time-step. While both LTS methods, either with or without intermediate local source evaluations, are second-order convergent in $\Delta t$ , our approach is more accurate (but also slightly more expensive) than the LTS method from [3] when the source strongly varies locally. Finally, we revisit the original example from [6] and show that the new weighted transition restores the expected optimal $H^{1}$ -convergence rate even when the source is nonzero across the coarse-to-fine mesh interface.

2 Stabilized LF-LTS Method

In Sections 2.1 and 2.2, we adopt and recall the setting for finite element spaces as in [20].

2.1 The Wave Equation

Let $\Omega\subset\mathbb{R}^{d}$ be a bounded Lipschitz domain with boundary $\Gamma=\partial\Omega$ . For $1\leq q\leq\infty$ , let $L^{q}\left(\Omega\right)$ be the standard Lebesgue space with norm $\left\|\cdot\right\|_{L^{q}\left(\Omega\right)}$ . For $q=2$ , the scalar product in $L^{2}\left(\Omega\right)$ is denoted by $\left(\cdot,\cdot\right)$ and the norm by $\left\|\cdot\right\|:=\left\|\cdot\right\|_{L^{2}\left(\Omega\right)}$ . For $k\in\mathbb{N}_{0}$ and $1\leq q\leq\infty$ , let $W^{k,q}\left(\Omega\right)$ denote the classical Sobolev spaces equipped with the norm $\left\|\cdot\right\|_{W^{k,q}\left(\Omega\right)}$ . For $q=2$ , these spaces are Hilbert spaces and denoted by $H^{k}\left(\Omega\right):=W^{k,2}\left(\Omega\right)$ with scalar product $\left(\cdot,\cdot\right)_{H^{k}\left(\Omega\right)}$ and norm $\left\|\cdot\right\|_{H^{k}\left(\Omega\right)}:=\left(\cdot,\cdot\right)_{H^{% k}\left(\Omega\right)}^{1/2}=\left\|\cdot\right\|_{W^{k,2}\left(\Omega\right)}$ .

We assume that the boundary is split into a Dirichlet and a Neumann part: $\Gamma=\Gamma_{D}\cup\Gamma_{N}$ , $\Gamma_{D}\cap\Gamma_{N}=\emptyset$ with the convention that $\Gamma_{D}$ is (relatively) closed and has positive surface measure. On $\Gamma_{D}$ homogeneous Dirichlet boundary conditions are imposed and Neumann boundary conditions on $\Gamma_{N}$ . The space of Sobolev functions in $H^{1}\left(\Omega\right)$ which vanish on the Dirichlet part of the boundary $\partial\Omega$ is denoted by $H_{D}^{1}\left(\Omega\right):=\left\{w\in H^{1}\left(\Omega\right):\left.w% \right|_{\Gamma_{D}}=0\right\}$ . For $\Gamma_{D}=\Gamma$ , we use the shorthand $H_{0}^{1}\left(\Omega\right)$ . Throughout this paper we restrict ourselves to function spaces over the field of real numbers.

Let $V\subset H^{1}\left(\Omega\right)$ denote a closed subspace of $H^{1}\left(\Omega\right)$ , such as $V=H^{1}\left(\Omega\right)$ or $V=H_{0}^{1}\left(\Omega\right)$ , and $a:V\times V\rightarrow\mathbb{R}$ denote a bilinear form, which is symmetric, continuous, and coercive:

a\left(u,v\right)=a\left(v,u\right)\qquad\forall u,v\in V

(2.1a)

and

\left|a\left(u,v\right)\right|\leq C_{\operatorname*{cont}}\left\|u\right\|_{H% ^{1}\left(\Omega\right)}\left\|v\right\|_{H^{1}\left(\Omega\right)}\qquad% \forall u,v\in V

(2.1b)

and

a\left(u,u\right)\geq c_{\operatorname*{coer}}\left\|u\right\|_{H^{1}\left(% \Omega\right)}^{2}\qquad\forall u\in V.

(2.1c)

For given $u_{0}\in V,v_{0}\in L^{2}\left(\Omega\right)$ and $F:\left[0,T\right]\rightarrow V^{\prime}$ , we consider the wave equation: Find $u:\left[0,T\right]\rightarrow V$ such that

\left(\ddot{u},w\right)+a\left(u,w\right)=F\left(w\right)\quad\forall w\in V,t>0

(2.2)

with initial conditions

u\left(0\right)=u_{0}\quad\text{and\quad}\dot{u}\left(0\right)=v_{0}.

(2.3)

It is well known that (2.2)–(2.3) is well-posed for sufficiently regular $u_{0}$ , $v_{0}$ and $F$ [27]. In fact, the weak solution $u$ can be shown to be continuous in time, that is, $u\in C^{0}(0,T;V),\dot{u}\in C^{0}(0,T;L^{2}\left(\Omega\right))$ – see [27, Chapter III, Theorems 8.1 and 8.2] for details – which implies that the initial conditions (2.3) are well defined. Moreover, we always assume that there exists a function $f:\left[0,T\right]\rightarrow L^{2}\left(\Omega\right)$ such that

F\left(t\right)\left(w\right)=\left(f\left(t\right),w\right)\qquad\forall w\in V% \quad\forall t\in\left[0,T\right].

For the stability and convergence analysis we will impose further smoothness assumptions on $f$ .

Example 1

The classical second-order wave equation in strong form is given by

\begin{split}\partial_{t}^{2}u-\nabla\cdot(c^{2}\nabla u)&=f\qquad\;\,\mbox{in% }\Omega\times(0,T),\\ u&=0\qquad\;\,\mbox{on }\Gamma_{D}\times(0,T),\\ \frac{\partial u}{\partial\nu}&=0\qquad\;\,\mbox{on }\Gamma_{N}\times(0,T),\\ u|_{t=0}&=u_{0}\qquad\mbox{in }\Omega,\\ u_{t}|_{t=0}&=v_{0}\qquad\mbox{in }\Omega,\\ &\end{split}

(2.4)

where the velocity field $c(x)$ satisfies $0<c_{\min}\leq c(x)\leq c_{\max}$ . In this case, we have $V:=H_{D}^{1}\left(\Omega\right):=\left\{w\in H^{1}\left(\Omega\right):\left.w% \right|_{\Gamma_{D}}=0\right\}$ ; the bilinear form is given by $a\left(u,v\right):=\left(c^{2}\nabla u,\nabla u\right)$ and the right-hand side by $F\left(w\right)=\left(f,w\right)$ for all $w\in V$ .

2.2 Galerkin Finite Element Spatial Discretization

We begin by briefly introducing our notation for $hp$ finite element spaces. Let $\mathcal{T}$ be a simplicial finite element mesh for $\Omega$ which is conforming in the sense that there are no hanging nodes. The local and global mesh width is

h_{\tau}:=\operatorname*{diam}\tau,\quad h:=\max_{\tau\in\mathcal{T}}h_{\tau}.

As a convention the triangles in $\mathcal{T}$ are closed sets. Without loss of generality, we may assume that there is fixed $h_{0}>0$ such that

h\leq h_{0}.

(2.5)

The finite element space of continuous, piecewise polynomials of degree $m\geq 1$ is denoted by

S_{\mathcal{T}}^{m}:=\left\{u\in C^{0}\left(\Omega\right)\mid\forall\tau\in% \mathcal{T}\quad\left.u\right|_{\tau}\in\mathbb{P}_{m}\right\},

where $\mathbb{P}_{m}$ is the space of polynomials of maximal total degree $m$ . The standard local nodal basis is denoted by $b_{m,z}$ , $z\in\Sigma_{\mathcal{T}}^{m}$ , where $\Sigma_{\mathcal{T}}^{m}$ is the full set of nodal points of order $m$ . The finite element space with incorporated Dirichlet boundary conditions is given by $S:=S_{\mathcal{T}}^{m}\cap V$ . The set of nodal points for $S$ is the subset $\Sigma:=\left\{z\in\Sigma_{\mathcal{T}}^{m}\mid z\in\Gamma_{\operatorname*{D}}\right\}$ and the basis in $S$ is denoted by $b_{z}:=b_{m,z}$ , $z\in\Sigma$ . Furthermore, we use the notion of a mass-lumped scalar product

\left(u,v\right)_{\mathcal{T}}:=\sum_{z\in\Sigma}d_{z}u\left(z\right)v\left(z% \right),\quad\forall u,v\in S

for certain coefficients $d_{z}>0$ and induced norm $\|\cdot\|_{\mathcal{T}}:=\left(\cdot,\cdot\right)_{\mathcal{T}}^{1/2}$ which is assumed to be equivalent to the $L^{2}\left(\Omega\right)$ -norm: there exists a constant $C_{\operatorname*{eq}}>0$ such that

C_{\operatorname*{eq}}^{-1}\left\|u\right\|_{\mathcal{T}}\leq\left\|u\right\|% \leq C_{\operatorname*{eq}}\left\|u\right\|_{\mathcal{T}}\quad\forall u\in S% \text{.}

(2.6)

For more details we refer, e.g., to [20, Definition 2.3].

We also associate to the discrete bilinear form the linear operator $A^{S}:S\rightarrow S$ defined by

\left(A^{S}u,v\right)_{\mathcal{T}}=a\left(u,v\right)\qquad\forall u,v\in S

(2.7)

and introduce the discrete right-hand side function $f_{S}:[0,T]\rightarrow S$ defined by

\left(f_{S}(t),v\right)_{\mathcal{T}}=F(t)(v)\qquad\forall\,v\in S\quad\forall% \,t\in[0,T].

With the definitions and notation at hand, the semi-discrete wave equation (with mass-lumping) is given by: find $u_{S}:\left[0,T\right]\rightarrow S$ such that

\left(\ddot{u}_{S},v\right)_{\mathcal{T}}+a\left(u_{S},v\right)=\left(f_{S}(t)% ,v\right)_{\mathcal{T}}\quad\forall v\in S,t>0

(2.8a)

with initial conditions

u_{S}\left(0\right)=r_{S}u_{0}\quad\text{and\quad}\dot{u}_{S}\left(0\right)=r_% {S}v_{0},

(2.8b)

where $r_{S}u_{0}$ is a finite element approximation of $u_{0}$ (see, e.g., [20, (2.8)]).

As in [20], we consider finite element meshes $\mathcal{T}$ which consist of a quasi-uniform part $\mathcal{T}_{\operatorname*{c}}\subset\mathcal{T}$ and allow for local refinements so that $\mathcal{T}_{\operatorname*{f}}:=\mathcal{T}\backslash\mathcal{T}_{% \operatorname*{c}}$ might contain much smaller triangles than $\mathcal{T}_{\operatorname{c}}$ . In this way, we assume that the mesh size of $\mathcal{T}_{\operatorname{c}}$ satisfies $h_{\operatorname{c}}=h$ . Similarly, we define coarse and locally refined regions of the domain¹¹1By $\operatorname*{int}\left(M\right)$ we denote the open interior of a subset $M\subset\mathbb{R}^{d}$ . as

\Omega_{\operatorname*{c}}:=\operatorname*{int}\left({\textstyle\bigcup% \nolimits_{\tau\in\mathcal{T}_{\operatorname*{c}}}}\tau\right),\quad\Omega_{% \operatorname*{f}}:=\operatorname*{int}\left({\textstyle\bigcup\nolimits_{\tau% \in\mathcal{T}_{\operatorname*{f}}}}\tau\right),\quad\Omega_{\operatorname*{f}% }^{+1}:=\operatorname*{int}\left({\textstyle\bigcup\nolimits_{\begin{subarray}% {c}\tau\in\mathcal{T}\\ \tau\cap\overline{\Omega_{\operatorname*{f}}}\neq\emptyset\end{subarray}}}\,% \tau\right).

(2.9)

Hence, $\Omega_{\operatorname*{f}}^{+1}$ contains $\Omega_{\operatorname*{f}}$ and all elements directly adjacent to it. We note that $\Omega_{\operatorname*{f}}$ and $\Omega_{\operatorname*{c}}$ are disjoint, while their union covers all of $\Omega$ except for the coarse/fine interface, and that $\Omega_{\operatorname*{c}}\cup\Omega_{\operatorname*{f}}^{+1}=\Omega$ since the interface between $\Omega_{\operatorname*{f}}$ and $\Omega_{\operatorname*{c}}$ is contained in $\Omega_{\operatorname*{f}}^{+1}$ .

We also split the degrees of freedom associated with the fine or coarse parts of the mesh, respectively, as

\Sigma_{\operatorname*{f}}:=\Sigma\cap\overline{\Omega_{\operatorname*{f}}}% \quad\text{and\quad}\Sigma_{\operatorname*{c}}:=\Sigma\backslash\Sigma_{% \operatorname*{f}},

and introduce the corresponding FE subspaces

S_{\operatorname*{c}}:=\operatorname{span}\left\{b_{z}:z\in\Sigma_{% \operatorname*{c}}\right\}\quad\mbox{and}\quad S_{\operatorname*{f}}:=% \operatorname{span}\left\{b_{z}:z\in\Sigma_{\operatorname*{f}}\right\},

which gives rise to a direct sum decomposition: for every $u\in S$ there exist unique functions $u_{\operatorname*{c}}\in S_{\operatorname*{c}}$ and $u_{\operatorname*{f}}\in S_{\operatorname*{f}}$ such that

u=u_{\operatorname*{c}}+u_{\operatorname*{f}}.

(2.10)

Hence, we can define the projections $\Pi_{\operatorname*{c}}^{S}:S\rightarrow S_{\operatorname*{c}}$ and $\Pi_{\operatorname*{f}}^{S}:S\rightarrow S_{\operatorname*{f}}$ by

\Pi_{\operatorname*{c}}^{S}u:=u_{\operatorname*{c}}\quad\text{and}\quad\Pi_{% \operatorname*{f}}^{S}u:=u_{\operatorname*{f}}.

(2.11)

In fact, the decomposition (2.10) is orthogonal with respect to the $\left(\cdot,\cdot\right)_{\mathcal{T}}$ scalar product [20, Lemma 3.2], which is essential to prove sharp bounds for the eigenvalues of our discrete bilinearform [20, Theorem 3.7].

2.3 Stabilized LF-LTS( $\nu$ ) Algorithm

The original leapfrog (LF) based local time-stepping (LF-LTS $(0)$ ) methods for the numerical solution of the second-order wave equations (2.2) were proposed for homogeneous right-hand sides in [13] and for inhomogeneous right-hand sides in [21]. Optimal convergence rates for the LF-LTS $(0)$ method from [13] with $p$ local time steps were derived for a conforming FEM discretization, albeit under a CFL condition where $\Delta t$ in fact depends on the smallest elements in the mesh [19]. To prove optimal $L^{2}$ convergence rates under a CFL condition independent of $p$ , a stabilized algorithm LF-LTS $(\nu)$ was introduced recently in [20], and also independently in [2]. Here, $\nu\geq 0$ denotes a small stabilization parameter; typically, we set $\nu=0.01$ . Stability and convergence for homogeneous right-hand sides were derived in [20] for a CFL condition independent of the coarse/fine mesh ratio.

We shall now extend this stabilized LF-LTS $(\nu)$ method to inhomogeneous right-hand sides so that for $\nu=0$ the method coincides with the original LF-LTS $(0)$ method [21, Sect. 4.1] for the inhomogeneous wave equation. To do so, let the (global) time-step $\Delta t=T/N$ and $u_{S}^{\left(n\right)}$ denote the FE approximation at time $t_{n}=n\Delta t$ . Given the numerical solution at times $t_{n-1}$ and $t_{n}$ , the LF-LTS method then computes the numerical solution of (2.8) at $t_{n+1}$ by using a smaller time-step $\Delta\tau=\Delta t/p$ inside the regions of local refinement; here, $p\geq 2$ denotes the “coarse” to “fine” time-step ratio.

Inside the fine region, the source $f$ is also evaluated over the time interval $[t_{n}-\Delta t,t_{n}+\Delta t]$ at the intermediate times $t_{n+k/p}=\left(n+k/p\right)\Delta t$ and we let

f_{S,k}^{(n)}:=f_{S}\left(t_{n}+\frac{k}{p}\Delta t\right),\quad-p\leq k\leq p.

The cost of those extra evaluations inside the locally refined region is typically negligible compared to the overall solution process. Still, a locally strongly varying source typically calls for such a a higher sampling of the right-hand side – see Section 5.1.

Given the fixed stabilization parameter $0\leq\nu\leq 1/2$ , we now define the following constants determined by stabilized versions [25] of Chebyshev polynomials of the first and second kind [32], denoted by $T_{n}$ and $U_{n}$ , respectively:

\delta_{p,\nu}:=1+\frac{\nu}{p^{2}},\qquad\omega_{p,\nu}:=2\,\frac{T_{p}^{% \prime}\left(\delta_{p,\nu}\right)}{T_{p}\left(\delta_{p,\nu}\right)},

(2.12)

and

\beta_{p,\nu}^{(k,\ell)}:=\frac{T_{k+\ell}\left(\delta_{p,\nu}\right)}{T_{k+1}% \left(\delta_{p,\nu}\right)}\quad\text{and}\quad\gamma_{p,\nu}^{(k)}:=\frac{% \left(p-k\right)\beta_{p,\nu}^{\left(k,p-k\right)}}{U_{p-1-k}\left(\delta_{p,% \nu}\right)}

(2.13)

for $k=0,1,\ldots,p-1$ and $\ell=-1,0,\ldots,p-k$ ; formally we set $T_{-1}=0$ . Later the special case

\beta_{p,\nu}^{(0,0)}=\frac{\gamma_{p,\nu}^{\left(0\right)}}{2p^{2}}\omega_{p,% \nu}=\frac{1}{T_{1}\left(\delta_{p,\nu}\right)}=\frac{1}{\delta_{p,\nu}}

(2.14)

will be employed which follows from (2.13) and the formula $T_{n}^{\prime}=nU_{n-1}$ (see [15, 18.9.21]).

The stabilized LF-LTS algorithm for the inhomogeneous wave equation which computes the new solution $u_{S}^{(n+1)}$ from the current and previous values $u_{S}^{(n)}$ and $u_{S}^{(n-1)}$ , is now given by:

Algorithm 1 (LF-LTS $(\nu)$ Galerkin FE Algorithm)

Let $n\geq 1$ .

Given $u_{S}^{(n-1)}$ , $u_{S}^{(n)}$ , compute $w_{S}^{(n)}$ as

w_{S}^{(n)}=\Pi_{\operatorname*{c}}^{S}f_{S,0}^{(n)}-A^{S}\Pi_{\operatorname*{% c}}^{S}u_{S}^{(n)}.

(2.15)

Compute

z_{S,1}^{(n)}=u_{S}^{\left(n\right)}+\frac{1}{2}\left(\frac{\Delta t}{p}\right% )^{2}\left(\frac{2p^{2}}{\omega_{p,\nu}}\,\beta_{p,\nu}^{(0,0)}\left(w_{S}^{(n% )}-A^{S}\Pi_{\operatorname*{f}}^{S}u_{S}^{(n)}\right)+\gamma_{p,\nu}^{(0)}\Pi_% {\operatorname*{f}}^{S}f_{S,0}^{(n)}\right).

(2.16)

For $k=1,\ldots,p-1,$ compute

	$\displaystyle z_{S,k+1}^{(n)}$	$\displaystyle=\left(1+\beta_{p,\nu}^{(k,-1)}\right)z_{S,k}^{(n)}-\beta_{p,\nu}% ^{(k,-1)}z_{S,k-1}^{(n)}$
		$\displaystyle\quad+\left(\frac{\Delta t}{p}\right)^{2}\left(\frac{2p^{2}}{% \omega_{p,\nu}}\,\beta_{p,\nu}^{(k,0)}\left(w_{S}^{(n)}-A^{S}\Pi_{% \operatorname{f}}^{S}z_{S,k}^{(n)}\right)+\gamma_{p,\nu}^{(k)}\frac{1}{2}\Pi_% {\operatorname{f}}^{S}\left(f_{S,k}^{(n)}+f_{S,-k}^{(n)}\right)\right).$

Compute

u_{S}^{(n+1)}=-u_{S}^{(n-1)}+2z_{S,p}^{(n)}.

The $\operatorname*{LF-LTS}(\nu)$ method is completed by specifying the initial values

	$\displaystyle u_{S}^{\left(0\right)}$	$\displaystyle=r_{S}u_{0},$
	$\displaystyle u_{S}^{\left(1\right)}$	$\displaystyle=r_{S}u_{0}+\Delta t\,r_{S}v_{0}+\frac{\Delta t^{2}}{2}\left(f_{S% }^{(0)}-A^{S}u_{0}\right).$		(3.20b)

Here, the term $A^{S}u_{0}$ in the second equation could be replaced by $A^{S,p,\nu}u_{0}$ thereby allowing for local time-stepping already during the very first time-step. In that case, the analysis below also applies but requires a minor change as explained in [20, Rem. 2.5]. This modification neither affects the stability nor the convergence rate of the overall LF-LTS $(\nu)$ scheme.

Remark 1

The above LF-LTS algorithm evaluates the source $f$ inside the refined region at all intermediate times $t_{n+k/p}=t_{n}+(k/p)\Delta t$ , $k=1-p,\dots,p-1$ . Hence for $\nu=0$ , it coincides with the original LF-LTS method for the inhomogeneous wave equation [21, Sect. 4.1], since

\delta_{p,0}=1,\quad\omega_{p,0}=2p^{2},\quad\beta_{p,0}^{(k,\ell)}=1,\quad% \gamma_{p,0}^{(k)}=\frac{p-k}{U_{p-1-k}(1)}=1

holds for all $k=0,1,\ldots,p-1$ , $\ell=-1,0,\ldots,p-k$ . In particular, for $\nu=0$ it corresponds to a standard LF discretization of the inhomogeneous wave equation with the smaller time-step $\Delta\tau=\Delta t/p$ . In contrast, the “split-LFC” local time integration method in [3] omits those intermediate evaluations of $f$ inside the fine region, akin to the previous LFC method [4]. That second variant from [3] without intermediate source evaluations is easily obtained by replacing in the above algorithm (2.15) with

w_{S}^{(n)}=f_{S,0}^{(n)}-A^{S}\Pi_{\operatorname*{c}}^{S}u_{S}^{(n)}

and setting $\gamma_{p,\nu}^{(k)}=0$ for $k=0,1,\ldots,p-1.$ Clearly if $f=0$ , both algorithms coincide.

3 Convergence Theory in the $L^{2}$ -norm

3.1 Two-Step Formulation via Chebyshev Polynomials

For the analysis, we will first rewrite the above LF-LTS $(\nu)$ algorithm in a two-step leapfrog-like formulation as it was done in [20] for homogeneous right-hand sides. For this, we define the polynomials

$\displaystyle P_{p,\nu,k}\left(x\right)$	$\displaystyle:=2\left(1-\frac{T_{k}\left(\delta_{p,\nu}-\frac{x}{\omega_{p,\nu% }}\right)}{T_{k}\left(\delta_{p,\nu}\right)}\right),$	(3.1)
$\displaystyle P_{p,\nu,k}^{\Delta t}\left(x\right)$	$\displaystyle:=\frac{P_{p,\nu,k}\left(\Delta t^{2}x\right)}{\Delta t^{2}x},$	(3.2)
$\displaystyle Q_{p,\nu,r,k}^{\Delta t}\left(x\right)$	$\displaystyle:=\frac{p-k}{p^{2}}\frac{T_{p}\left(\delta_{p,\nu}\right)}{T_{k}% \left(\delta_{p,\nu}\right)}\frac{U_{k-1-r}\left(\delta_{p,\nu}-\frac{\Delta t% ^{2}x}{\omega_{p,\nu}}\right)}{U_{p-1-r}\left(\delta_{p,\nu}\right)},$	(3.3)

where formally we set $U_{-1}=0$ and write short $P_{p,\nu}:=P_{p,\nu,p}$ , $P_{p,\nu}^{\Delta t}:=P_{p,\nu,p}^{\Delta t}$ , and $Q_{p,\nu,r}^{\Delta t}:=Q_{p,\nu,r,p}^{\Delta t}$ in what follows.

Remark 2

The univariate polynomials $P_{p,\nu,k}^{\Delta t}$ , $Q_{p,\nu,r,k}^{\Delta t}$ will appear in an explicit representation formula for the functions $z_{S,k}^{(n)}$ at the local time steps. Some simple properties are listed below:

P_{p,\nu,k}^{\Delta t}\in\mathbb{P}_{k-1}\qquad Q_{p,\nu,r,k}^{\Delta t}\in% \mathbb{P}_{k-1-r}

and the first two terms in a Taylor series around $x=0$ are given by

	$\displaystyle P_{p,\nu,k}^{\Delta t}\left(x\right)$	$\displaystyle=\frac{2}{\omega_{p,\nu}}\frac{T_{k}^{\prime}\left(\delta_{p,\nu}% \right)}{T_{k}\left(\delta_{p,\nu}\right)}-x\frac{\Delta t^{2}}{\omega_{p,\nu}% ^{2}}\frac{T_{k}^{\prime\prime}\left(\delta_{p,\nu}\right)}{T_{k}\left(\delta_% {p,\nu}\right)}+\ldots$		(3.4)
	$\displaystyle Q_{p,\nu,r,k}^{\Delta t}\left(x\right)$	$\displaystyle=\frac{p-k}{p^{2}}\frac{T_{p}\left(\delta_{p,\nu}\right)}{T_{k}% \left(\delta_{p,\nu}\right)}\left(1-x\frac{\Delta t^{2}}{\omega_{p,\nu}}\frac{% U_{k-1-r}^{\prime}\left(\delta_{p,\nu}\right)}{U_{p-1-r}\left(\delta_{p,\nu}% \right)}+\ldots\right).$		(3.5)

Note that the expression for $Q_{p,\nu,r,k}^{\Delta t}\left(x\right)$ is valid only for $r\leq k-1$ while for $k=r$ we have $Q_{p,\nu,r,r+1}^{\Delta t}=0$ . Since $\delta_{p,\nu}$ lies outside the interval $\left]-1,1\right[$ the denominators in these expressions cannot be zero. For the special case $k=1$ it is straightforward to verify

P_{p,\nu,1}^{\Delta t}(x)=Q_{p,\nu,0,1}^{\Delta t}(x)=\frac{\gamma_{p,\nu}^{(0% )}}{p^{2}}.

(3.6)

The following recurrence relations for $P_{p,\nu,k}$ follow easily from the recurrence relations for Chebyshev polynomials as already stated in [4, Lem. 6.1] and [20, Lem. B.1].

Lemma 1

Let $r=0,1,\ldots,p-1$ . Then, the polynomials $P_{p,\nu,k}$ defined in (3.1) satisfy the recurrence relation

$\displaystyle P_{p,\nu,0}\left(x\right)$	$\displaystyle=0,$
$\displaystyle P_{p,\nu,1}\left(x\right)$	$\displaystyle=\frac{2x}{\delta_{p,\nu}\omega_{p,\nu}},$	(3.7)
$\displaystyle P_{p,\nu,k+1}\left(x\right)$	$\displaystyle=2\beta_{p,\nu}^{\left(k,0\right)}\left(\delta_{p,\nu}-\frac{x}{% \omega_{p,\nu}}\right)P_{p,\nu,k}\left(x\right)-\beta_{p,\nu}^{\left(k,-1% \right)}P_{p,\nu,k-1}\left(x\right)+\frac{4}{\omega_{p,\nu}}\beta_{p,\nu}^{% \left(k,0\right)}x,$

while $P_{p,\nu,k}^{\Delta t}$ satisfies

$\displaystyle P_{p,\nu,0}^{\Delta t}(x)$	$\displaystyle=0,$
$\displaystyle P_{p,\nu,1}^{\Delta t}(x)$	$\displaystyle=\frac{2}{\delta_{p,\nu}\omega_{p,\nu}},$	(3.8)
$\displaystyle P_{p,\nu,k+1}^{\Delta t}(x)$	$\displaystyle=2\beta_{p,\nu}^{\left(k,0\right)}\left(\delta_{p,\nu}-\frac{% \Delta t^{2}x}{\omega_{p,\nu}}\right)P_{p,\nu,k}^{\Delta t}\left(x\right)-% \beta_{p,\nu}^{\left(k,-1\right)}P_{p,\nu,k-1}^{\Delta t}\left(x\right)+\frac{% 4}{\omega_{p,\nu}}\beta_{p,\nu}^{\left(k,0\right)}.$	(3.9)

Lemma 2

Let $r=0,1,\ldots,p-1$ . Then, the polynomials $Q_{p,\nu,r,k}^{\Delta t}$ defined in (3.3) satisfy the recurrence relation

$\displaystyle Q_{p,\nu,r,r}^{\Delta t}(x)$	$\displaystyle=0,$		(3.10)
$\displaystyle Q_{p,\nu,r,r+1}^{\Delta t}(x)$	$\displaystyle=\frac{\gamma_{p,\nu}^{(r)}}{p^{2}},$		(3.11)
$\displaystyle Q_{p,\nu,r,k+1}^{\Delta t}(x)$	$\displaystyle=2\beta_{p,\nu}^{(k,0)}\,\left(\delta_{p,\nu}-\frac{\Delta t^{2}}% {\omega_{p,\nu}}x\right)Q_{p,\nu,r,k}^{\Delta t}(x)$
	$\displaystyle\quad-\beta_{p,\nu}^{(k,-1)}Q_{p,\nu,r,k-1}^{\Delta t}(x),$	$\displaystyle r+1\leq k\leq p-1.$	(3.12)

Proof. The first equation (3.10) holds by (3.3) and (3.11) follows directly since $U_{0}\equiv 1$ . Hence, let $r+1\leq k\leq p-1$ and insert definitions (2.13), (3.3) into the right-hand side of (3.12) to obtain

	$\displaystyle\left(1+\beta_{p,\nu}^{(k,-1)}-\left(\frac{\Delta t}{p}\right)^{2% }\frac{2p^{2}}{\omega_{p,\nu}}\beta_{p,\nu}^{(k,0)}\,x\right)Q_{p,\nu,r,k}^{% \Delta t}(x)-\beta_{p,\nu}^{(k,-1)}Q_{p,\nu,r,k-1}^{\Delta t}\left(x\right)$
	$\displaystyle=\left(\frac{T_{k+1}\left(\delta_{p,\nu}\right)+T_{k-1}\left(% \delta_{p,\nu}\right)-2\frac{\Delta t^{2}x}{\omega_{p,\nu}}T_{k}\left(\delta_{% p,\nu}\right)}{T_{k}\left(\delta_{p,\nu}\right)}U_{k-1-r}\left(\delta_{p,\nu}-% \frac{\Delta t^{2}x}{\omega_{p,\nu}}\right)\right.$
	$\displaystyle\qquad\left.-U_{k-2-r}\left(\delta_{p,\nu}-\frac{\Delta t^{2}x}{% \omega_{p,\nu}}\right)\right)\frac{\gamma_{p,\nu}^{(r)}}{p^{2}}\frac{T_{r+1}% \left(\delta_{p,\nu}\right)}{T_{k+1}\left(\delta_{p,\nu}\right)}.$

Using the standard recursions for Chebyshev polynomials of first and second kind [15, Table 18.9.1], respectively, we conclude

	$\displaystyle\left(1+\beta_{p,\nu}^{(k,-1)}-\left(\frac{\Delta t}{p}\right)^{2% }\frac{2p^{2}}{\omega_{p,\nu}}\beta_{p,\nu}^{(k,0)}\,x\right)Q_{p,\nu,r,k}^{% \Delta t}(x)-\beta_{p,\nu}^{(k,-1)}Q_{p,\nu,r,k-1}^{\Delta t}(x)$
	$\displaystyle=\frac{\gamma_{p,\nu}^{(r)}}{p^{2}}\frac{T_{r+1}\left(\delta_{p,% \nu}\right)}{T_{k+1}\left(\delta_{p,\nu}\right)}\left(2\left(\delta_{p,\nu}-% \frac{\Delta t^{2}x}{\omega_{p,\nu}}\right)U_{k-1-r}\left(\delta_{p,\nu}-\frac% {\Delta t^{2}x}{\omega_{p,\nu}}\right)-U_{k-2-r}\left(\delta_{p,\nu}-\frac{% \Delta t^{2}x}{\omega_{p,\nu}}\right)\right)$
	$\displaystyle=\frac{\gamma_{p,\nu}^{(r)}}{p^{2}}\frac{T_{r+1}\left(\delta_{p,% \nu}\right)}{T_{k+1}\left(\delta_{p,\nu}\right)}U_{k-r}\left(\delta_{p,\nu}-% \frac{\Delta t^{2}x}{\omega_{p,\nu}}\right)$
	$\displaystyle=Q_{p,\nu,r,k+1}^{\Delta t}(x).$

Lemma 3

For $0\leq k\leq p$ , the functions $z_{S,k}^{(n)}$ defined in Algorithm 1 can be written in the form

z_{S,k}^{(n)}=\mathfrak{L}_{S,k}^{(n)}u_{S}^{(n)}+\frac{\Delta t^{2}}{2}\left(% \frac{2}{\omega_{p,\nu}}\frac{T_{k}^{\prime}\left(\delta_{p,\nu}\right)}{T_{k}% \left(\delta_{p,\nu}\right)}\Pi_{\operatorname*{c}}^{S}f_{S,0}^{(n)}+\sum_{r=-% (k-1)}^{k-1}Q_{p,\nu,|r|,k}^{\Delta t}\left(A^{S}\Pi_{\operatorname*{f}}^{S}% \right)\Pi_{\operatorname*{f}}^{S}f_{S,r}^{(n)}\right)

(3.13)

for

\mathfrak{L}_{S,k}^{(n)}=\left(I^{S}-\frac{\Delta t^{2}}{2}P_{p,\nu,k}^{\Delta t% }\left(A^{S}\Pi_{\operatorname*{f}}^{S}\right)A^{S}\right).

Proof. In [20, Lemma B.2, (B.2)], relation (3.13) is shown for $f_{S}=0$ . In the inhomogeneous case here, it follows by the same arguments that

z_{S,k}^{(n)}=\mathfrak{L}_{S,k}^{(n)}u_{S}^{(n)}+\frac{\Delta t^{2}}{2}% \mathfrak{R}_{S,k}^{(n)},

(3.14)

where the term $\mathfrak{R}_{S,k}^{(n)}$ only depends on the right-hand side evaluations $\left(f_{S,r}^{(n)}\right)_{r=-\left(k-1\right)}^{k-1}$ .

Hence, it remains to prove that

\mathfrak{R}_{S,k}^{(n)}=\frac{2}{\omega_{p,\nu}}\frac{T_{k}^{\prime}\left(% \delta_{p,\nu}\right)}{T_{k}\left(\delta_{p,\nu}\right)}\Pi_{\operatorname*{c}% }^{S}f_{S,0}^{(n)}+\sum_{r=-(k-1)}^{k-1}Q_{p,\nu,|r|,k}^{\Delta t}\left(A^{S}% \Pi_{\operatorname*{f}}^{S}\right)\Pi_{\operatorname*{f}}^{S}f_{S,r}^{\left(n% \right)}.

(3.15)

From $\Pi_{\operatorname{f}}^{S}\Pi_{\operatorname{c}}^{S}f=0$ (see (2.11)) it follows:

P_{p,\nu,k}^{\Delta t}\left(A^{S}\Pi_{\operatorname*{f}}^{S}\right)\Pi_{% \operatorname*{c}}^{S}f_{S,0}^{(n)}=P_{p,\nu,k}^{\Delta t}\left(0\right)\Pi_{% \operatorname*{c}}^{S}f_{S,0}^{(n)}\overset{\text{(\ref{firsttermexp})}}{=}% \frac{2}{\omega_{p,\nu}}\frac{T_{k}^{\prime}\left(\delta_{p,\nu}\right)}{T_{k}% \left(\delta_{p,\nu}\right)}\Pi_{\operatorname*{c}}^{S}f_{S,0}^{(n)}

so that

\mathfrak{R}_{S,k}^{(n)}=P_{p,\nu,k}^{\Delta t}\left(A^{S}\Pi_{\operatorname*{% f}}^{S}\right)\Pi_{\operatorname*{c}}^{S}f_{S,0}^{(n)}+\sum_{r=-(k-1)}^{k-1}Q_% {p,\nu,|r|,k}^{\Delta t}\left(A^{S}\Pi_{\operatorname*{f}}^{S}\right)\Pi_{% \operatorname*{f}}^{S}f_{S,r}^{\left(n\right)}

(3.16)

implies (3.15).

We prove (3.16) by induction. For $k=0$ , this is trivial, and for $k=1$ , we use the definition of $z_{S,1}^{\left(n\right)}$ in Algorithm 1, insert (3.11) for $r=0$ (see also Remark 2), employ $w_{S}^{(n)}=\Pi_{\operatorname*{c}}^{S}f_{S,0}^{(n)}-A^{S}\Pi_{\operatorname*{% c}}^{S}u_{S}^{(n)}$ (see (2.15)) and (3.14) to obtain

	$\displaystyle z_{S,1}^{(n)}$	$\displaystyle=u_{S}^{\left(n\right)}+\frac{1}{2}\left(\frac{\Delta t}{p}\right% )^{2}\left(\frac{2p^{2}}{\omega_{p,\nu}}\,\beta_{p,\nu}^{(0,0)}\left(w_{S}^{(n% )}-A^{S}\Pi_{\operatorname{f}}^{S}u_{S}^{(n)}\right)+\gamma_{p,\nu}^{(0)}\Pi_% {\operatorname{f}}^{S}f_{S,0}^{(n)}\right)$
		$\displaystyle=\mathfrak{L}_{S,1}^{(n)}u_{S}^{(n)}+\frac{1}{2}\left(\frac{% \Delta t}{p}\right)^{2}\left(\frac{2p^{2}}{\omega_{p,\nu}}\,\beta_{p,\nu}^{(0,% 0)}\Pi_{\operatorname{c}}^{S}f_{S,0}^{(n)}+\gamma_{p,\nu}^{(0)}\Pi_{% \operatorname{f}}^{S}f_{S,0}^{(n)}\right)$
		$\displaystyle\overset{\text{(\ref{eq:recursion_Q_pnurk_1})}}{=}\mathfrak{L}_{S% ,1}^{(n)}u_{S}^{(n)}+\frac{\Delta t^{2}}{2}\left(\frac{2}{\omega_{p,\nu}}\beta% _{p,\nu}^{(0,0)}\Pi_{\operatorname{c}}^{S}f_{S,0}^{(n)}+Q_{p,\nu,0,1}^{\Delta t% }\left(A^{S}\Pi_{\operatorname{f}}^{S}\right)\Pi_{\operatorname*{f}}^{S}f_{S,% 0}\right).$

From $\beta_{p,\nu}^{\left(0,0\right)}=\frac{T_{1}^{\prime}\left(\delta_{p,\nu}% \right)}{T_{1}\left(\delta_{p,\nu}\right)}$ (cf. (2.14)) it follows that the second summand is $\frac{\Delta t^{2}}{2}\mathfrak{R}_{S,1}^{(n)}$ .

For $k>1$ , Step 3 in Algorithm 1 in combination with the induction hypothesis, the definition of $w_{S}^{\left(n\right)}$ , and the fact that the linear part with respect $u_{S}^{\left(n\right)}$ is known from (3.14) leads to

	$\displaystyle z_{S,k+1}^{(n)}$	$\displaystyle=\left(1+\beta_{p,\nu}^{(k,-1)}\right)z_{S,k}^{(n)}-\beta_{p,\nu}% ^{(k,-1)}z_{S,k-1}^{(n)}$
		$\displaystyle\quad+\left(\frac{\Delta t}{p}\right)^{2}\left(\frac{2p^{2}}{% \omega_{p,\nu}}\,\beta_{p,\nu}^{(k,0)}\left(w_{S}^{(n)}-A^{S}\Pi_{% \operatorname{f}}^{S}z_{S,k}^{(n)}\right)+\gamma_{p,\nu}^{(k)}\frac{1}{2}\Pi_% {\operatorname{f}}^{S}\left(f_{S,k}^{(n)}+f_{S,-k}^{(n)}\right)\right)$
		$\displaystyle=\mathfrak{L}_{S,k+1}^{(n)}u_{S}^{(n)}+\frac{\Delta t^{2}}{2}\rho$

for

	$\displaystyle\rho:=\,$	$\displaystyle 2\beta_{p,\nu}^{(k,0)}\left(\delta_{p,\nu}-\frac{\Delta t^{2}}{% \omega_{p,\nu}}A^{S}\Pi_{\operatorname*{f}}^{S}\right)\mathfrak{R}_{S,k}^{(n)}% -\beta_{p,\nu}^{(k,-1)}\mathfrak{R}_{S,k-1}^{(n)}$
		$\displaystyle\quad+4\frac{\,\beta_{p,\nu}^{(k,0)}}{\omega_{p,\nu}}\Pi_{% \operatorname{c}}^{S}f_{S,0}^{(n)}+\frac{\gamma_{p,\nu}^{(k)}}{p^{2}}\Pi_{% \operatorname{f}}^{S}\left(f_{S,k}^{(n)}+f_{S,-k}^{(n)}\right),$

where we used

\left(1+\beta_{p,\nu}^{(k,-1)}\right)/\left(2\beta_{p,\nu}^{(k,0)}\right)=% \frac{T_{k+1}\left(\delta_{p,\nu}\right)+T_{k-1}\left(\delta_{p,\nu}\right)}{2% T_{k}\left(\delta_{p,\nu}\right)}\overset{\text{\cite[cite]{[\@@bibref{}{NIST:% DLMF}{}{}, Table 18.9.1]}}}{=}\delta_{p,\nu}.

It remains to prove that $\rho=\mathfrak{R}_{S,k+1}^{(n)}$ . We use the induction hypotheses for $\mathfrak{R}_{S,k}^{(n)}$ as in (3.15), (3.11) and obtain

	$\displaystyle\rho$	$\displaystyle:=\left.\left(2\beta_{p,\nu}^{(k,0)}\left(\delta_{p,\nu}-\frac{% \Delta t^{2}}{\omega_{p,\nu}}x\right)P_{p,\nu,k}^{\Delta t}\left(x\right)-% \beta_{p,\nu}^{(k,-1)}P_{p,\nu,k-1}^{\Delta t}\left(x\right)+4\frac{\,\beta_{p% ,\nu}^{(k,0)}}{\omega_{p,\nu}}\right)\right\|_{x=A^{S}\Pi_{\operatorname{f}}^{% S}}\Pi_{\operatorname{c}}^{S}f_{S,0}^{(n)}$
		$\displaystyle+\sum_{r=-(k-1)}^{k-1}\left.\left(2\beta_{p,\nu}^{(k,0)}\left(% \delta_{p,\nu}-\frac{\Delta t^{2}}{\omega_{p,\nu}}x\right)Q_{p,\nu,\|r\|,k}^{% \Delta t}\left(x\right)-\beta_{p,\nu}^{(k,-1)}Q_{p,\nu,\|r\|,k-1}^{\Delta t}% \left(x\right)\right)\right\|_{x=A^{S}\Pi_{\operatorname{f}}^{S}}\Pi_{% \operatorname{f}}^{S}f_{S,r}^{\left(n\right)}$
		$\displaystyle+Q_{p,\nu,k,k+1}^{\Delta t}\left(A^{S}\Pi_{\operatorname{f}}^{S}% \right)\Pi_{\operatorname{f}}^{S}\left(f_{S,k}^{(n)}+f_{S,-k}^{(n)}\right)$

where we used $Q_{p,\nu,k-1,k-1}^{\Delta t}=0$ . The first row can be simplified by using the recurrence for $P_{p,\nu,k+1}^{\Delta t}$ (see (3.7)). For the second one we employ the recursion (3.12) for $Q_{p,\nu,|r|,k}^{\Delta t}$ so that the second row equals

	$\displaystyle\sum_{r=-(k-1)}^{k-1}\left.Q_{p,\nu,r,k+1}^{\Delta t}\left(x% \right)\right\|_{x=A^{S}\Pi_{\operatorname{f}}^{S}}\Pi_{\operatorname{f}}^{S}% f_{S,r}^{\left(n\right)}+Q_{p,\nu,k,k+1}^{\Delta t}\left(A^{S}\Pi_{% \operatorname{f}}^{S}\right)\Pi_{\operatorname{f}}^{S}\left(f_{S,k}^{(n)}+f_% {S,-k}^{(n)}\right)$
	$\displaystyle\qquad=\sum_{r=-k}^{k}\left.Q_{p,\nu,r,k+1}^{\Delta t}\left(x% \right)\right\|_{x=A^{S}\Pi_{\operatorname{f}}^{S}}\Pi_{\operatorname{f}}^{S}% f_{S,r}^{\left(n\right)}.$

In this way we have proved (3.16) and, in turn, (3.15).

Corollary 4

Let $a^{p,\nu}:S\times S\rightarrow\mathbb{R}$ denote the discrete bilinear form

a^{p,\nu}\left(u,v\right):=a\left(P_{p,\nu}^{\Delta t}\left(\Pi_{\operatorname% *{f}}^{S}A^{S}\right)u,v\right)\qquad\forall u,v\in S

(3.17)

with associated (linear) operator $A^{S,p,\nu}:S\rightarrow S$

	$\displaystyle A^{S,p,\nu}$	$\displaystyle:=\left(A^{S}\right)^{1/2}P_{p,\nu}^{\Delta t}\left(\left(A^{S}% \right)^{1/2}\Pi_{\operatorname*{f}}^{S}\left(A^{S}\right)^{1/2}\right)\left(A% ^{S}\right)^{1/2}$		(3.18)
		$\displaystyle=P_{p,\nu}^{\Delta t}\left(A^{S}\Pi_{\operatorname{f}}^{S}\right% )A^{S}=A^{S}P_{p,\nu}^{\Delta t}\left(\Pi_{\operatorname{f}}^{S}A^{S}\right)$

and let $R_{S}:[0,T]\rightarrow S$ by defined by

R_{S}\left(t\right):=P_{p,\nu}^{\Delta t}\left(A^{S}\Pi_{\operatorname*{f}}^{S% }\right)\Pi_{\operatorname*{c}}^{S}f_{S}(t)+\sum_{r=-(p-1)}^{p-1}Q_{p,\nu,|r|}% ^{\Delta t}\left(A^{S}\Pi_{\operatorname*{f}}^{S}\right)\Pi_{\operatorname*{f}% }^{S}f_{S}\left(t+\frac{r\Delta t}{p}\right).

(3.19)

Then, for $n\geq 1$ , the $\operatorname*{LF-LTS}(\nu)$ algorithm is equivalent to

\left(u_{S}^{(n+1)}-2\,u_{S}^{(n)}+u_{S}^{(n-1)},w\right)_{\mathcal{T}}=\Delta t% ^{2}\left[\left(R_{S}^{(n)},w\right)_{\mathcal{T}}-a^{p,\nu}\left(u_{S}^{(n)},% w\right)\right]\quad\forall\,w\in S,

(3.20a)

where $R_{S}^{(n)}:=R_{S}\left(t_{n}\right)$ .

Proof. The assertion follows directly from inserting (3.13) for $k=p$ into Step 4 of Algorithm 1.

3.2 Main Convergence Results

The time-steps $\Delta t$ and $\Delta\tau=\Delta t/p$ , $p\in\mathbb{N}$ , used in $\mathcal{T}_{\operatorname*{c}}$ and $\mathcal{T}_{\operatorname*{f}}$ , respectively, are each determined by the smallest triangle in either part of the mesh. Hence, the “coarse”-to-“fine” time-step ratio, $p$ , must satisfy

\dfrac{\min\left\{h_{\tau}:\tau\in\mathcal{T}_{\operatorname*{c}}\right\}}{h_{% \tau}}\leq p\qquad\forall\tau\in\mathcal{T}.

Given the quasi-uniformity constant for the coarse mesh $\mathcal{T}_{\operatorname*{c}}$ ,

C_{\operatorname*{qu},\operatorname*{c}}:=\frac{h_{\operatorname*{c}}}{\min% \left\{h_{\tau}:\tau\in\mathcal{T}_{\operatorname*{c}}\right\}},

the ratio between the largest and smallest elements of the mesh is thus bounded by

\dfrac{h_{\operatorname*{c}}}{h_{\tau}}\leq C_{\operatorname*{qu},% \operatorname*{c}}p\qquad\forall\tau\in\mathcal{T}.

(3.21)

In [20, Section 3.1], theoretical properties for the bilinear form $a^{p,\nu}$ are derived, which will be used in the error analysis below. We impose the same CFL-type stability condition as in [20, (3.8)] on the ratio $\Delta t/h_{\operatorname*{c}}$ :

\left(3+\frac{C_{\operatorname*{cont}}}{c_{\operatorname*{coer}}}\right)C_{% \operatorname*{cont}}C_{\operatorname*{inv}}^{2}\left(\frac{\Delta t}{h_{% \operatorname*{c}}}\right)^{2}\leq\frac{\nu}{\nu+1},

(3.22)

where $C_{\operatorname*{inv}}$ is a constant independent of $p$ such that the inverse inequalities holds (cf. [20, (3.3), (3.4)]):

\left.\begin{array}[c]{c}\left\|u\right\|_{H^{1}\left(\Omega_{\operatorname{c}% }\right)}\leq C_{\operatorname*{inv}}h_{\operatorname*{c}}^{-1}\left\|u\right% \|_{\mathcal{T}}\\ \left\|u\right\|_{H^{1}\left(\Omega\right)}\leq C_{\operatorname*{inv}}\dfrac{% p}{h_{\operatorname*{c}}}\left\|u\right\|_{\mathcal{T}}\end{array}\right\}% \quad\forall u\in S.

(3.23)

For a discussion and a possible simplification of the CFL condition (3.22) we refer to [20, (3.8)].

To derive a priori error estimates for the LF-LTS $(\nu)$ Galerkin FE solution of (3.20), we first introduce

v_{S}^{\left(n+1/2\right)}:=\frac{u_{S}^{\left(n+1\right)}-u_{S}^{\left(n% \right)}}{\Delta t},

and rewrite (3.20) as a one-step method

$\displaystyle\left(v_{S}^{\left(n+1/2\right)},q\right)_{\mathcal{T}}$	$\displaystyle=\left(v_{S}^{\left(n-1/2\right)},q\right)_{\mathcal{T}}-\Delta t% \,a^{p,\nu}\left(u_{S}^{\left(n\right)},q\right)+\Delta t\left(R_{S}^{(n)},q% \right)_{\mathcal{T}}\;$	$\displaystyle\forall q\in S,$	(3.24)
$\displaystyle\left(u_{S}^{\left(n+1\right)},r\right)_{\mathcal{T}}$	$\displaystyle=\left(u_{S}^{\left(n\right)},r\right)_{\mathcal{T}}+\Delta t% \left(v_{S}^{\left(n+1/2\right)},r\right)_{\mathcal{T}}\;$	$\displaystyle\forall r\in S,$
$\displaystyle\left(u_{S}^{\left(0\right)},w\right)_{\mathcal{T}}$	$\displaystyle=\left(u_{0},w\right)_{\mathcal{T}}$	$\displaystyle\forall w\in S,$
$\displaystyle\left(v_{S}^{\left(1/2\right)},w\right)_{\mathcal{T}}$	$\displaystyle=\left(v_{0},w\right)_{\mathcal{T}}+\frac{\Delta t}{2}\left(\left% (f_{S}^{(0)},w\right)_{\mathcal{T}}-a\left(u_{0},w\right)\right)\;$	$\displaystyle\forall w\in S.$

The first two equations in (3.24) correspond to the one-step iteration

\left[\begin{array}[c]{c}v_{S}^{(n+1/2)}\\ u_{S}^{(n+1)}\end{array}\right]=\mathfrak{B}\left[\begin{array}[c]{c}v_{S}^{(n% -1/2)}\\ u_{S}^{(n)}\end{array}\right]+\Delta t\begin{bmatrix}R_{S}^{(n)}\\ \Delta t\,R_{S}^{(n)}\end{bmatrix}\quad n\geq 1

with

\mathfrak{B}:=\left[\begin{array}[c]{cc}I_{S}&-\Delta tA^{S,p,\nu}\\ \Delta t\,I_{S}&I_{S}-\Delta t^{2}A^{S,p,\nu}\end{array}\right]

and $A^{S,p,\nu}$ as in (3.18).

Next, we denote the error by

\mathbf{e}^{\left(n+1\right)}:=\begin{bmatrix}e_{v}^{\left(n+1/2\right)}\\ e_{u}^{\left(n+1\right)}\end{bmatrix}:=\begin{bmatrix}v\left(t_{n+1/2}\right)-% v_{S}^{\left(n+1/2\right)}\\ u\left(t_{n+1}\right)-u_{S}^{\left(n+1\right)}\end{bmatrix},

where $u$ is the solution of (2.2)-(2.3) and $v$ the solution of the corresponding first-order formulation: Find $u,v:[0,T]\rightarrow V$ such that

	$\displaystyle\left(\dot{v},w\right)+a\left(u,w\right)$	$\displaystyle=\left(f\left(t\right),w\right)$	$\displaystyle\quad\forall\,w$	$\displaystyle\in V,\quad t>0,$
	$\displaystyle\left(v,w\right)$	$\displaystyle=\left(\dot{u},w\right)$	$\displaystyle\quad\forall\,w$	$\displaystyle\in V,\quad t>0$

with initial conditions $u(0)=u_{0}$ and $v(0)=v_{0}$ .

We shall split the error into a semi-discrete and a fully discrete contribution. To do so, we introduce the first-order formulation of the semi-discrete problem (2.8). Find $u_{S},v_{S}:\left[0,T\right]\rightarrow S$ such that

\begin{array}[c]{cl}\left.\begin{array}[c]{rcl}\left(\dot{v}_{S},w\right)_{% \mathcal{T}}+a\left(u_{S},w\right)&=&\left(f_{S}\left(t\right),w\right)_{% \mathcal{T}}\\ \left(v_{S},w\right)_{\mathcal{T}}&=&\left(\dot{u}_{S},w\right)_{\mathcal{T}}% \end{array}\right\}&\forall w\in S,\quad t>0,\\ \begin{array}[c]{rcl}\qquad u_{S}\left(0\right)&=&r_{S}u_{0}\\ \qquad v_{S}\left(0\right)&=&r_{S}v_{0}.\end{array}&\end{array}

Hence, we may write

\mathbf{e}^{\left(n+1\right)}=\mathbf{e}_{S}^{\left(n+1\right)}+\mathbf{e}_{S,% \Delta t}^{\left(n+1\right)}

(3.25)

with

	$\displaystyle\mathbf{e}_{S}^{\left(n+1\right)}$	$\displaystyle:=\left[\begin{array}[c]{c}e_{v,S}^{\left(n+1/2\right)}\\ e_{u,S}^{\left(n+1\right)}\end{array}\right]:=\left[\begin{array}[c]{c}v\left(% t_{n+1/2}\right)-v_{S}\left(t_{n+1/2}\right)\\ u\left(t_{n+1}\right)-u_{S}\left(t_{n+1}\right)\end{array}\right],$
	$\displaystyle\mathbf{e}_{S,\Delta t}^{\left(n+1\right)}$	$\displaystyle:=\left[\begin{array}[c]{c}e_{v,S,\Delta t}^{\left(n+1/2\right)}% \\ e_{u,S,\Delta t}^{\left(n+1\right)}\end{array}\right]:=\left[\begin{array}[c]{% c}v_{S}\left(t_{n+1/2}\right)-v_{S}^{\left(n+1/2\right)}\\ u_{S}\left(t_{n+1}\right)-u_{S}^{\left(n+1\right)}\end{array}\right].$

Following similar arguments as in [19, §3.2], we derive a recurrence relation for $\mathbf{e}_{S,\Delta t}^{(n+1)}$ , which can be solved and eventually yields the explicit error representation

	$\displaystyle\left[\begin{array}[c]{c}e_{v,S,\Delta t}^{\left(n+1/2\right)}\\ e_{u,S,\Delta t}^{\left(n+1\right)}\end{array}\right]$	$\displaystyle=\mathfrak{B}^{n}\left[\begin{array}[c]{c}e_{v,S,\Delta t}^{\left% (1/2\right)}\\ e_{u,S,\Delta t}^{\left(1\right)}\end{array}\right]+\Delta t\sum_{\ell=1}^{n-1% }\mathfrak{B}^{\ell}\left(\mbox{\boldmath$\sigma$}^{\left(n-\ell\right)}-\mbox% {\boldmath$\sigma$}^{\left(n-\ell+1\right)}\right)$		(3.30)
		$\displaystyle\quad+\Delta t\mbox{\boldmath$\sigma$}^{\left(n\right)}-\Delta t% \mathfrak{B}^{n}\mbox{\boldmath$\sigma$}^{\left(1\right)},$

where

\mbox{\boldmath$\sigma$}^{\left(n\right)}=\frac{1}{\Delta t}\begin{bmatrix}-% \dfrac{u_{S}\left(t_{n+1}\right)-u_{S}\left(t_{n}\right)}{\Delta t}+v_{S}\left% (t_{n+1/2}\right)\\ u_{S}\left(t_{n}\right)+\left(A^{S,p,\nu}\right)^{-1}\left(\dfrac{v_{S}\left(t% _{n+1/2}\right)-v_{S}\left(t_{n-1/2}\right)}{\Delta t}-R_{S}^{(n)}\right)\end{bmatrix}

(3.31)

with $R_{S}^{(n)}$ defined as in Corollary 4.

Before proving the main convergence theorem we provide some approximation estimates.

Lemma 5

Let $R_{S}$ be defined as in (3.19) and let $\theta\in\left[t_{n},t_{n+1}\right]$ , $\nu>0$ . Then, there exist constants $C_{1}$ , $C_{2}$ independent of $\Delta t$ , $h$ , $p$ , and $f_{S}$ such that

\left\|\left(A^{S,p,\nu}\right)^{-1}\partial_{t}^{3}R_{S}\left(\theta\right)% \right\|_{\mathcal{T}}\leq C_{1}\max_{t\in[t_{n-1},t_{n+2}]}{\left\|\partial_{% t}^{3}f_{S}\left(t\right)\right\|_{\mathcal{T}}}

(3.32)

and

\left\|\left(A^{S,p,\nu}\right)^{-1}\partial_{t}^{\ell}\left(f_{S}^{(n)}-R_{S}% ^{(n)}\right)\right\|_{\mathcal{T}}\leq C_{2}\Delta t^{2}\max_{k=0,2}\left\{% \max_{t\in[t_{n},t_{n+1}]}{\left\|\partial_{t}^{k+\ell}f_{S}\left(t\right)% \right\|_{\mathcal{T}}}\right\}

(3.33)

for $\ell=0,1$ .

Proof. For an operator $K:S\rightarrow S$ we introduce the operator norm

\left\|K\right\|_{S\leftarrow S}:=\sup_{w\in S\backslash\left\{0\right\}}\frac% {\left\|Kw\right\|_{\mathcal{T}}}{\left\|w\right\|_{\mathcal{T}}}.

We start from the definition

	$\displaystyle\left(A^{S,p,\nu}\right)^{-1}\partial_{t}^{\ell}R_{S}\left(t\right)$	$\displaystyle=\left(A^{S,p,\nu}\right)^{-1}P_{p,\nu}^{\Delta t}\left(A^{S}\Pi_% {\operatorname{f}}^{S}\right)\Pi_{\operatorname{c}}^{S}\partial_{t}^{\ell}f_% {S}\left(t\right)$		(3.34)
		$\displaystyle+\left(A^{S,p,\nu}\right)^{-1}\sum_{r=-(p-1)}^{p-1}Q_{p,\nu,\|r\|}^% {\Delta t}\left(A^{S}\Pi_{\operatorname{f}}^{S}\right)\Pi_{\operatorname{f}}% ^{S}\partial_{t}^{\ell}f_{S}\left(t+\frac{r\Delta t}{p}\right).$

Since $\Pi_{\operatorname{f}}^{S}\Pi_{\operatorname{c}}^{S}f_{S}=0$ we obtain from (3.4) and (2.12) for the first term in (3.34)

P_{p,\nu}^{\Delta t}\left(A^{S}\Pi_{\operatorname*{f}}^{S}\right)\Pi_{% \operatorname*{c}}^{S}\partial_{t}^{\ell}f_{S}\left(t\right)=\Pi_{% \operatorname*{c}}^{S}\partial_{t}^{\ell}f_{S}\left(t\right).

(3.35)

and hence

\left\|P_{p,\nu}^{\Delta t}\left(A^{S}\Pi_{\operatorname*{f}}^{S}\right)\Pi_{% \operatorname*{c}}^{S}\partial_{t}^{\ell}f_{S}\left(t\right)\right\|_{\mathcal% {T}}\leq\left\|\partial_{t}^{\ell}f_{S}\left(t\right)\right\|_{\mathcal{T}}

From [20, Theorem 3.7] it follows

\left\|\left(A^{S,p,\nu}\right)^{-1}v\right\|_{\mathcal{T}}\leq\frac{1}{c_{\nu% }}\left\|v\right\|_{\mathcal{T}},

with

c_{\nu}=\frac{1}{2}\min\left\{\frac{c_{\operatorname*{coer}}}{C_{\operatorname% *{eq}}^{2}},\frac{2\nu}{\left(\nu+1\right)\Delta t^{2}}\right\}

(3.36)

and, in turn,

\left\|\left(A^{S,p,\nu}\right)^{-1}\partial_{t}^{3}P_{p,\nu}^{\Delta t}\left(% A^{S}\Pi_{\operatorname*{f}}^{S}\right)\Pi_{\operatorname*{c}}^{S}\partial_{t}% ^{\ell}f_{S}\left(t\right)\right\|_{\mathcal{T}}\leq\frac{1}{c_{\nu}}\left\|% \partial_{t}^{\ell}f_{S}\left(t\right)\right\|_{\mathcal{T}}.

(3.37)

Next we consider the second summand in (3.34). The decomposition operator $E^{S}:S\rightarrow S_{\operatorname{c}}\times S_{\operatorname{f}}$ is defined by $E^{S}u=\left(\Pi_{\operatorname{c}}^{S}u,\Pi_{\operatorname{f}}^{S}u\right)$ (cf. (2.11)) so that its left-inverse $C^{S}:S_{\operatorname{c}}\times S_{\operatorname{f}}\rightarrow S$ is given by $C^{S}\left(u_{\operatorname{c}},u_{\operatorname{f}}\right):=u_{\operatorname{% c}}+u_{\operatorname{f}}$ . For an operator $K:S\rightarrow S$ and $s,t\in\left\{\operatorname{c},\operatorname{f}\right\}$ we set $K_{st}:=\Pi_{s}^{S}K\Pi_{t}^{S}$ .

Then, it follows from [20, (3.17), Lem. 3.6]

	$\displaystyle\left(A^{S,p,\nu}\right)^{-1}$	$\displaystyle=\left(A^{S}\right)^{-1}+\left(\left(A^{S,p,\nu}\right)^{-1}-% \left(A^{S}\right)^{-1}\right)$
		$\displaystyle=\left(A^{S}\right)^{-1}+C^{S}\left[\begin{array}[c]{cc}0&0\\ 0&\left(A_{\operatorname{f}\operatorname{f}}^{S,p,\nu}\right)^{-1}-\left(A_{% \operatorname{f}\operatorname{f}}^{S}\right)^{-1}\end{array}\right]E^{S}$
		$\displaystyle=\left(A^{S}\right)^{-1}+\left(\left(A_{\operatorname{f}% \operatorname{f}}^{S,p,\nu}\right)^{-1}-\left(A_{\operatorname{f}\operatorname% {f}}^{S}\right)^{-1}\right)\Pi_{\operatorname{f}}^{S}$

and

	$\displaystyle\left(A^{S,p,\nu}\right)^{-1}Q_{p,\nu,\|r\|}^{\Delta t}\left(A^{S}% \Pi_{\operatorname{f}}^{S}\right)\Pi_{\operatorname{f}}^{S}$	$\displaystyle=\left(A^{S}\right)^{-1}Q_{p,\nu,\|r\|}^{\Delta t}\left(A^{S}\Pi_{% \operatorname{f}}^{S}\right)\Pi_{\operatorname{f}}^{S}$		(3.38)
		$\displaystyle+\left(\left(A_{\operatorname{f}\operatorname{f}}^{S,p,\nu}\right% )^{-1}-\left(A_{\operatorname{f}\operatorname{f}}^{S}\right)^{-1}\right)Q_{p,% \nu,\|r\|}^{\Delta t}\left(A_{\operatorname{f}\operatorname{f}}^{S}\right)\Pi_{% \operatorname*{f}}^{S},$

where we used that $\Pi_{\operatorname{f}}^{S}$ is a projection. Clearly $A_{\operatorname{f}\operatorname{f}}^{S}$ and $Q_{p,\nu,|r|}^{\Delta t}\left(A_{\operatorname{f}\operatorname{f}}^{S}\right)$ are self-adjoint with respect to $\left(\cdot,\cdot\right)_{\mathcal{T}}$ and positive semidefinite. Furthermore, [20, Lem. 3.6] implies

\left\|\left(A_{\operatorname{f}\operatorname{f}}^{S,p,\nu}\right)^{-1}-\left(% A_{\operatorname{f}\operatorname{f}}^{S}\right)^{-1}\right\|_{S\leftarrow S}% \leq\frac{\nu+1}{2\nu}\Delta t^{2}

and by applying well-known spectral theory we get

\left\|\left(\left(A_{\operatorname{f}\operatorname{f}}^{S,p,\nu}\right)^{-1}-% \left(A_{\operatorname{f}\operatorname{f}}^{S}\right)^{-1}\right)Q_{p,\nu,|r|}% ^{\Delta t}\left(A_{\operatorname{f}\operatorname{f}}^{S}\right)\Pi_{% \operatorname*{f}}^{S}\right\|_{\mathcal{T\leftarrow T}}\leq\frac{\nu+1}{2\nu}% \Delta t^{2}\max_{x\in\mathcal{\sigma}\left(A_{\operatorname{f}\operatorname{f% }}^{S}\right)}\left|Q_{p,\nu,|r|}^{\Delta t}\left(x\right)\right|.

(3.39)

Next, we estimate the maximum on the right-hand side. The combination of (3.3) and (2.13) with $Q_{p,\nu,r}^{\Delta t}:=Q_{p,\nu,r,p}^{\Delta t}$ leads to

Q_{p,\nu,\left|r\right|}^{\Delta t}(x)=\frac{p-\left|r\right|}{p^{2}}\frac{\,U% _{p-1-\left|r\right|}\left(\delta_{p,\nu}-\frac{\Delta t^{2}x}{\omega_{p,\nu}}% \right)}{U_{p-1-\left|r\right|}\left(\delta_{p,\nu}\right)}.

(3.40)

From [20, (3.21)] it follows that the spectrum of $A_{\operatorname{f}\operatorname{f}}^{S}$ is contained in $\left[0,\frac{2\delta_{p,\nu}\omega_{p,\nu}}{\Delta t^{2}}\right]$ so that the argument $\delta_{p,\nu}-\frac{\Delta t^{2}x}{\omega_{p,\nu}}$ in (3.40) is contained in $\left[-\delta_{p,\nu},\delta_{p,\nu}\right]$ . From [15, 18.7.3, 18.14.1] we conclude that

\left\|U_{p-1-\left|r\right|}\right\|_{L^{\infty}\left(\left[-1,1\right]\right% )}=U_{p-1-\left|r\right|}\left(1\right)

and since $U_{p-1-\left|r\right|}$ behaves monotonically outside $\left[-1,1\right]$ we obtain by symmetry

\,\left|U_{p-1-\left|r\right|}\left(\delta_{p,\nu}-\frac{\Delta t^{2}x}{\omega% _{p,\nu}}\right)\right|\leq U_{p-1-\left|r\right|}\left(\delta_{p,\nu}\right)% \quad\forall x\in\mathcal{\sigma}\left(A_{\operatorname{f}\operatorname{f}}^{S% }\right).

(3.41)

This leads to

\left|Q_{p,\nu,\left|r\right|}^{\Delta t}\left(x\right)\right|\leq\frac{p-% \left|r\right|}{p^{2}}\quad\forall x\in\mathcal{\sigma}\left(A_{\operatorname{% f}\operatorname{f}}^{S}\right).

By inserting this into (3.39) we get

\left\|\left(\left(A_{\operatorname{f}\operatorname{f}}^{S,p,\nu}\right)^{-1}-% \left(A_{\operatorname{f}\operatorname{f}}^{S}\right)^{-1}\right)Q_{p,\nu,|r|}% ^{\Delta t}\left(A_{\operatorname{f}\operatorname{f}}^{S}\right)\Pi_{% \operatorname*{f}}^{S}\right\|_{\mathcal{T\leftarrow T}}\leq\frac{\nu+1}{2\nu}% \Delta t^{2}\frac{p-\left|r\right|}{p^{2}}.

Next we estimate the first summand in (3.38) and write

\left(A^{S}\right)^{-1}Q_{p,\nu,|r|}^{\Delta t}\left(A^{S}\Pi_{\operatorname*{% f}}^{S}\right)\Pi_{\operatorname*{f}}^{S}=Q_{p,\nu,|r|}^{\Delta t}\left(0% \right)\left(A^{S}\right)^{-1}\Pi_{\operatorname*{f}}^{S}+\tilde{Q}_{p,\nu,|r|% }^{\Delta t}\left(A_{\operatorname{f}\operatorname{f}}^{S}\right)\Pi_{% \operatorname*{f}}^{S},

(3.42)

where

\tilde{Q}_{p,\nu,|r|}^{\Delta t}\left(x\right):=\frac{Q_{p,\nu,|r|}^{\Delta t}% \left(x\right)-Q_{p,\nu,|r|}^{\Delta t}\left(0\right)}{x}.

The expansion (3.5) leads to

	$\displaystyle\left\\|Q_{p,\nu,\|r\|}^{\Delta t}\left(0\right)\left(A^{S}\right)^{% -1}\Pi_{\operatorname*{f}}^{S}\right\\|_{S\leftarrow S}$	$\displaystyle=\frac{p-r}{p^{2}}\left\\|\left(A^{S}\right)^{-1}\Pi_{% \operatorname*{f}}^{S}\right\\|_{S\leftarrow S}$		(3.43)
		$\displaystyle\overset{\text{\cite[cite]{[\@@bibref{}{grote2021stabilized}{}{},% (3.7)]}}}{\leq}\frac{C_{\operatorname{eq}}^{2}}{c_{\operatorname{coer}}}% \frac{p-r}{p^{2}}\left\\|\Pi_{\operatorname{f}}^{S}\right\\|_{S\leftarrow S}% \leq\frac{C_{\operatorname*{eq}}^{2}}{c_{\operatorname{coer}}}\frac{p-r}{p^{2}}.$

For the estimate of $\tilde{Q}_{p,\nu,|r|}^{\Delta t}$ we employ (3.40) for

\tilde{Q}_{p,\nu,|r|}^{\Delta t}\left(x\right)=\frac{p-\left|r\right|}{p^{2}}% \frac{\,U_{p-1-\left|r\right|}\left(\delta_{p,\nu}-\frac{\Delta t^{2}x}{\omega% _{p,\nu}}\right)-U_{p-1-\left|r\right|}\left(\delta_{p,\nu}\right)}{U_{p-1-% \left|r\right|}\left(\delta_{p,\nu}\right)x}.

The combination of a Taylor argument and the monotonicity of orthogonal polynomials outside $\left[-1,1\right]$ results in the estimate

\left|\tilde{Q}_{p,\nu,|r|}^{\Delta t}\left(x\right)\right|\leq\frac{p-\left|r% \right|}{p^{2}}\frac{U_{p-1-\left|r\right|}^{\prime}\left(\delta_{p,\nu}\right% )}{U_{p-1-\left|r\right|}\left(\delta_{p,\nu}\right)}\frac{\Delta t^{2}}{% \omega_{p,\nu}}\quad\forall x\in\mathcal{\sigma}\left(A_{\operatorname{f}% \operatorname{f}}^{S}\right).

We employ Lemma 8 to estimate the ratio of the Chebyshev polynomials and obtain

\left|\tilde{Q}_{p,\nu,|r|}^{\Delta t}\left(x\right)\right|\leq\frac{p-\left|r% \right|}{p^{2}}\frac{\left(p-\left|r\right|\right)^{2}\operatorname*{e}^{\nu/2% }}{3}\frac{\Delta t^{2}}{\omega_{p,\nu}}\quad\forall x\in\mathcal{\sigma}\left% (A_{\operatorname{f}\operatorname{f}}^{S}\right).

(3.44)

The estimate [20, Lem. A.2] of $\omega_{p,\nu}$ implies

\left|\tilde{Q}_{p,\nu,|r|}^{\Delta t}\left(x\right)\right|\leq\frac{p-\left|r% \right|}{p^{2}}\frac{\left(p-\left|r\right|\right)^{2}\operatorname*{e}^{\nu}}% {6p^{2}}\Delta t^{2}\leq\frac{1}{6}\frac{p-\left|r\right|}{p^{2}}\operatorname% *{e}\nolimits^{\nu}\Delta t^{2}.

(3.45)

This finishes the estimate of the left-hand side in (3.42) by combining (3.43) and (3.45):

\left\|\left(A^{S}\right)^{-1}Q_{p,\nu,|r|}^{\Delta t}\left(A^{S}\Pi_{% \operatorname*{f}}^{S}\right)\Pi_{\operatorname*{f}}^{S}\right\|_{S\leftarrow S% }\leq\frac{p-r}{p^{2}}\left(\frac{C_{\operatorname*{eq}}^{2}}{c_{\operatorname% {coer}}}+\frac{\operatorname*{e}\nolimits^{\nu}\Delta t^{2}}{6}\right).

With this estimate at hand, the second term in (3.34) can be estimated by

	$\displaystyle\left\\|\left(A^{S,p,\nu}\right)^{-1}\sum_{r=-(p-1)}^{p-1}Q_{p,\nu% ,\|r\|}^{\Delta t}\left(A^{S}\Pi_{\operatorname{f}}^{S}\right)\Pi_{% \operatorname{f}}^{S}\partial_{t}^{\ell}f_{S}\left(t+\frac{r\Delta t}{p}% \right)\right\\|_{\mathcal{T}}$		(3.46)
	$\displaystyle\qquad\leq\left(\frac{C_{\operatorname{eq}}^{2}}{c_{% \operatorname{coer}}}+\frac{\operatorname{e}\nolimits^{\nu}\Delta t^{2}}{6}% \right)\sum_{r=-(p-1)}^{p-1}\frac{p-\left\|r\right\|}{p^{2}}\max_{t\in[t_{m-1},t% _{m+2}]}\left\\|\partial_{t}^{\ell}f_{S}\left(t\right)\right\\|_{\mathcal{T}}$
	$\displaystyle\qquad=\left(\frac{C_{\operatorname{eq}}^{2}}{c_{\operatorname{% coer}}}+\frac{\operatorname{e}\nolimits^{\nu}\Delta t^{2}}{6}\right)\max_{t% \in[t_{m-1},t_{m+2}]}\left\\|\partial_{t}^{\ell}f_{S}\left(t\right)\right\\|_{% \mathcal{T}}.$

For the second estimate (3.33) we use

	$\displaystyle\partial_{t}^{\ell}R_{S}\left(t\right)$	$\displaystyle=P_{p,\nu}^{\Delta t}\left(A^{S}\Pi_{\operatorname{f}}^{S}\right% )\Pi_{\operatorname{c}}^{S}f_{S}^{\left(\ell\right)}(t)+\left(\sum_{r=-(p-1)}% ^{p-1}Q_{p,\nu,\|r\|}^{\Delta t}\left(A^{S}\Pi_{\operatorname{f}}^{S}\right)% \right)\Pi_{\operatorname{f}}^{S}\partial_{t}^{\ell}f_{S}\left(t\right)$		(3.47)
		$\displaystyle+\sum_{r=1}^{p-1}Q_{p,\nu,r}^{\Delta t}\left(A^{S}\Pi_{% \operatorname{f}}^{S}\right)\Pi_{\operatorname{f}}^{S}\left(\partial_{t}^{% \ell}f_{S}\left(t+\frac{r\Delta t}{p}\right)+\partial_{t}^{\ell}f_{S}\left(t-% \frac{r\Delta t}{p}\right)-2\partial_{t}^{\ell}f_{S}\left(t\right)\right).$

A Taylor argument for the last difference yields

\left\|\partial_{t}^{\ell}f_{S}\left(t+\frac{r\Delta t}{p}\right)+\partial_{t}% ^{\ell}f_{S}\left(t-\frac{r\Delta t}{p}\right)-2\partial_{t}^{\ell}f_{S}\left(% t\right)\right\|_{\mathcal{T}}\leq\sup_{t\in\left[t-\Delta t,t+\Delta t\right]% }\left\|\partial_{t}^{\ell+2}f_{S}\left(t\right)\right\|_{\mathcal{T}}\left(% \frac{r\Delta t}{p}\right)^{2}

and as in (3.46) we get

	$\displaystyle\left\\|\left(A^{S,p,\nu}\right)^{-1}\sum_{r=1}^{p-1}Q_{p,\nu,\|r\|}% ^{\Delta t}\left(A^{S}\Pi_{\operatorname{f}}^{S}\right)\Pi_{\operatorname{f}% }^{S}\left(\partial_{t}^{\ell}f_{S}\left(t+\frac{r\Delta t}{p}\right)+\partial% _{t}^{\ell}f_{S}\left(t-\frac{r\Delta t}{p}\right)-2\partial_{t}^{\ell}f_{S}% \left(t\right)\right)\right\\|_{\mathcal{T}}$
	$\displaystyle\qquad\leq\Delta t^{2}\left(\frac{C_{\operatorname{eq}}^{2}}{c_{% \operatorname{coer}}}+\frac{\operatorname{e}\nolimits^{\nu}\Delta t^{2}}{6}% \right)\sup_{t\in\left[t-\Delta t,t+\Delta t\right]}\left\\|\partial_{t}^{\ell+% 2}f_{S}\left(t\right)\right\\|_{\mathcal{T}}.$

For the first term in (3.47) we use (3.35)

P_{p,\nu}^{\Delta t}\left(A^{S}\Pi_{\operatorname*{f}}^{S}\right)\Pi_{% \operatorname*{c}}^{S}\partial_{t}^{\ell}f_{S}(t)=\Pi_{\operatorname*{c}}^{S}% \partial_{t}^{\ell}f_{S}(t).

For the second one we write

\sum_{r=-(p-1)}^{p-1}Q_{p,\nu,|r|}^{\Delta t}\left(x\right)=1+\rho\left(x\right)

for

\rho\left(x\right):=\sum_{r=-(p-1)}^{p-1}\frac{p-\left|r\right|}{p^{2}}\frac{% \,U_{p-1-\left|r\right|}\left(\delta_{p,\nu}-\frac{\Delta t^{2}x}{\omega_{p,% \nu}}\right)-U_{p-1-\left|r\right|}\left(\delta_{p,\nu}\right)}{U_{p-1-\left|r% \right|}\left(\delta_{p,\nu}\right)}.

Then,

$\displaystyle\left\\|\left(A^{S,p,\nu}\right)^{-1}\partial_{t}^{\ell}\left(f_{S% }-R_{S}\right)\right\\|_{\mathcal{T}}$	$\displaystyle\leq\left\\|\left(A^{S,p,\nu}\right)^{-1}\partial_{t}^{\ell}\left(% f_{S}-\Pi_{\operatorname{c}}^{S}f_{S}-\Pi_{\operatorname{f}}^{S}f_{S}\right)% \right\\|_{\mathcal{T}}$	(3.48)
	$\displaystyle+\left\\|\left(A^{S,p,\nu}\right)^{-1}\rho\left(A^{S}\Pi_{% \operatorname{f}}^{S}\right)\Pi_{\operatorname{f}}^{S}\partial_{t}^{\ell}f_{% S}\left(t\right)\right\\|_{\mathcal{T}}$
	$\displaystyle+\Delta t^{2}\left(\frac{C_{\operatorname{eq}}^{2}}{c_{% \operatorname{coer}}}+\frac{\operatorname{e}\nolimits^{\nu}\Delta t^{2}}{6}% \right)\sup_{t\in\left[t-\Delta t,t+\Delta t\right]}\left\\|\partial_{t}^{\ell+% 2}f_{S}\left(t\right)\right\\|_{\mathcal{T}}.$

Clearly, the first term vanishes and it remains to estimate the second one. Analogous arguments as for (3.38) lead to

\left(A^{S,p,\nu}\right)^{-1}\rho\left(A^{S}\Pi_{\operatorname*{f}}^{S}\right)% \Pi_{\operatorname*{f}}^{S}=\left(A^{S}\right)^{-1}\rho\left(A^{S}\Pi_{% \operatorname*{f}}^{S}\right)\Pi_{\operatorname*{f}}^{S}+\left(\left(A_{% \operatorname{f}\operatorname{f}}^{S,p,\nu}\right)^{-1}-\left(A_{\operatorname% {f}\operatorname{f}}^{S}\right)^{-1}\right)\rho\left(A_{\operatorname{f}% \operatorname{f}}^{S}\right)\Pi_{\operatorname*{f}}^{S}

(3.49)

and, in turn, as in (3.39):

\left\|\left(\left(A_{\operatorname{f}\operatorname{f}}^{S,p,\nu}\right)^{-1}-% \left(A_{\operatorname{f}\operatorname{f}}^{S}\right)^{-1}\right)\rho\left(A_{% \operatorname{f}\operatorname{f}}^{S}\right)\Pi_{\operatorname*{f}}^{S}\right% \|_{\mathcal{T\leftarrow T}}\leq\frac{\nu+1}{2\nu}\Delta t^{2}\max_{x\in% \mathcal{\sigma}\left(A_{\operatorname{f}\operatorname{f}}^{S}\right)}\left|% \rho\left(x\right)\right|.

It holds

\left|\rho\left(x\right)\right|\leq\sum_{r=-(p-1)}^{p-1}\frac{p-\left|r\right|% }{p^{2}}\left|\frac{\,U_{p-1-\left|r\right|}\left(\delta_{p,\nu}-\frac{\Delta t% ^{2}x}{\omega_{p,\nu}}\right)}{U_{p-1-\left|r\right|}\left(\delta_{p,\nu}% \right)}-1\right|

and as in (3.41) we obtain

\left|\frac{\,U_{p-1-\left|r\right|}\left(\delta_{p,\nu}-\frac{\Delta t^{2}x}{% \omega_{p,\nu}}\right)}{U_{p-1-\left|r\right|}\left(\delta_{p,\nu}\right)}% \right|\leq 1\qquad\forall x\in\mathcal{\sigma}\left(A_{\operatorname{f}% \operatorname{f}}^{S}\right).

The combination of these estimates implies

\left\|\left(\left(A_{\operatorname{f}\operatorname{f}}^{S,p,\nu}\right)^{-1}-% \left(A_{\operatorname{f}\operatorname{f}}^{S}\right)^{-1}\right)\rho\left(A_{% \operatorname{f}\operatorname{f}}^{S}\right)\Pi_{\operatorname*{f}}^{S}\right% \|_{\mathcal{T\leftarrow T}}\leq\frac{\nu+1}{\nu}\Delta t^{2}.

It remains to estimate the first summand in the right-hand side of (3.49) and we proceed as in (3.42):

\left(A^{S}\right)^{-1}\rho\left(A^{S}\Pi_{\operatorname*{f}}^{S}\right)\Pi_{% \operatorname*{f}}^{S}=\rho\left(0\right)\left(A^{S}\right)^{-1}\Pi_{% \operatorname*{f}}^{S}+\tilde{\rho}\left(A_{\operatorname{f}\operatorname{f}}^% {S}\right)\Pi_{\operatorname*{f}}^{S}

for

\tilde{\rho}\left(x\right):=\frac{\rho\left(x\right)-\rho\left(0\right)}{x}.

Since it holds by definition that $\rho\left(0\right)=0$ , only the term containing $\tilde{\rho}$ has to be estimated. It is easy to see that

\tilde{\rho}\left(x\right)=\frac{\rho\left(x\right)}{x}=\sum_{r=-(p-1)}^{p-1}% \tilde{Q}_{p,\nu,|r|}^{\Delta t}\left(x\right).

The estimate (3.45) then leads to

\left\|\tilde{\rho}\left(A_{\operatorname{f}\operatorname{f}}^{S}\right)\Pi_{% \operatorname*{f}}^{S}\right\|_{\mathcal{T\leftarrow T}}\leq\sum_{r=-(p-1)}^{p% -1}\frac{1}{6}\frac{p-\left|r\right|}{p^{2}}\operatorname*{e}\nolimits^{\nu}% \Delta t^{2}=\frac{\operatorname*{e}\nolimits^{\nu}}{6}\Delta t^{2}.

This finishes the proof of an estimate for the left-hand side in (3.49):

\left\|\left(A^{S,p,\nu}\right)^{-1}\rho\left(A^{S}\Pi_{\operatorname*{f}}^{S}% \right)\Pi_{\operatorname*{f}}^{S}\right\|_{\mathcal{T}\leftarrow\mathcal{T}}% \leq\left(\frac{\operatorname*{e}\nolimits^{\nu}}{6}+\frac{\nu+1}{\nu}\right)% \Delta t^{2}.

By inserting this into (3.48) proves the second estimate (3.33):

\left\|\left(A^{S,p,\nu}\right)^{-1}\partial_{t}^{\ell}\left(f_{S}-R_{S}\right% )\right\|_{\mathcal{T}}\leq C_{3}\Delta t^{2}\max_{k\in\left\{0,2\right\}}\sup% _{t\in\left[t-\Delta t,t+\Delta t\right]}\left\|\partial_{t}^{\ell+k}f_{S}% \left(t\right)\right\|_{\mathcal{T}},

for $C_{3}=\tfrac{C_{\operatorname*{eq}}^{2}}{c_{\operatorname{coer}}}+\tfrac{\nu+1% }{\nu}+\tfrac{\operatorname*{e}\nolimits^{\nu}}{6}\left(1+\Delta t^{2}\right)$ .

Remark 3

The “split-LFC” method proposed in [3] for the inhomogeneous wave equation evaluates the right-hand side only at each global time-step. The convergence analysis for such a coarser sampling of $f$ easily follows as a special case from of our theory: In Lemma 3, the representation for that variant is obtained by setting $\Pi_{\operatorname{c}}^{S}$ to the identity operator and the last summation $\sum_{-k-1}^{k+1}\ldots$ to zero. Lemma 5 then holds verbatim while its proof already ends with estimate (3.37), since the analysis of the Chebyshev polynomials of second kind is only needed for a fully resolved right-hand side. As a consequence, the following two Theorems 6 and 7 also hold for that variant, that is, when $f$ is evaluated only once per global time-step, since the proofs remain identical.

Theorem 6

Assume that (2.1), (2.5), (2.6), (3.21) and (3.22) hold and let the solution and right-hand side of the semi-discrete equation (2.8) satisfy $u_{S}\in W^{5,\infty}\left(\left[0,T\right];L^{2}\left(\Omega\right)\right)$ and $f_{S}\in W^{3,\infty}\left(\left[0,T\right];L^{2}\left(\Omega\right)\right)$ , respectively. Then the fully discrete solution $u_{S}^{\left(n+1\right)}$ of (3.20) satisfies the error estimate

\left\|e_{u,S,\Delta t}^{\left(n+1\right)}\right\|\leq C\Delta t^{2}\left(1+T% \right)\mathcal{M}\left(u_{S},f_{S}\right)

with

\mathcal{M}\left(u_{S},f_{S}\right):=\max\left\{\max_{2\leq\ell\leq 5}\left\|% \partial_{t}^{\ell}u_{S}\right\|_{L^{\infty}\left(\left[0,T\right];L^{2}(% \Omega)\right)},\max_{0\leq\ell\leq 3}\left\|\partial_{t}^{\ell}f_{S}\right\|_% {L^{\infty}\left(\left[0,T\right];L^{2}(\Omega)\right)}\right\}

and a constant $C$ , which is independent of $n$ , $\Delta t$ , $T$ , $h$ , $p$ , $u_{S}$ , and $f_{S}$ .

Proof. The proof follows the lines of the proofs of [19, Theorem 16] and [20, Theorem 3.10]. The LF-LTS $(\nu)$ scheme (3.24) is stable [20, Theorem 3.9], i.e., for $f_{S}\equiv 0$ ,

\left\|v_{S}^{\left(n{+1/2}\right)}\right\|_{\mathcal{T}}+\left\|u_{S}^{\left(% n\right)}\right\|_{\mathcal{T}}\leq C_{0}\left(\left\|v_{S}^{\left(1/2\right)}% \right\|_{\mathcal{T}}+\left\|u_{S}^{\left(1\right)}\right\|_{\mathcal{T}}% \right),

where $C_{0}$ is independent of $n$ , $\Delta t$ , $h$ , and $T$ . We apply the stability estimate to the second component of the error representation (3.30) and obtain²²2For a pair of functions $\mathbf{v}=\left(v_{1},v_{2}\right)^{\intercal}\in S\times S$ we use the notation $\left\|\mathbf{v}\right\|_{\ell^{1}}:=\left\|v_{1}\right\|_{\mathcal{T}}+\left% \|v_{2}\right\|_{\mathcal{T}}$ .

	$\displaystyle\left\\|e_{u,S,\Delta t}^{\left(n+1\right)}\right\\|_{\mathcal{T}}$	$\displaystyle\leq C_{0}\left\\|e_{S,\Delta t}^{\left(1\right)}\right\\|_{% \mathcal{T}}+C_{0}\Delta t\sum_{\ell=1}^{n-1}\left\\|\mbox{\boldmath$\sigma$}^{% \left(n-\ell\right)}-\mbox{\boldmath$\sigma$}^{\left(n-\ell+1\right)}\right\\|_% {\ell^{1}}$
		$\displaystyle+C_{0}\Delta t\left\\|\mbox{\boldmath$\sigma$}^{\left(1\right)}% \right\\|_{\ell^{1}}+\Delta t\left\\|\mbox{\boldmath$\sigma$}^{\left(n\right)}% \right\\|_{\ell^{1}}.$		(3.50)

First, we consider the last term in the right-hand side of (3.50). By similar arguments as in the proof of [20, Theorem 3.10], we obtain

\Delta t\left\|\mbox{\boldmath$\sigma$}^{\left(n\right)}\right\|_{\ell^{1}}% \leq\dfrac{\Delta t^{2}}{24}\left(24\dfrac{\nu+1}{\nu}+\dfrac{1}{c_{\nu}}+1% \right)\mathcal{M}_{n}\left(u_{S},f_{S}\right)+\left\|\left(A^{S,p,\nu}\right)% ^{-1}\left(f_{S}^{(n)}-R_{S}^{(n)}\right)\right\|_{\mathcal{T}}

(3.51)

with $c_{\nu}$ as in (3.36) and

\mathcal{M}_{n}\left(u_{S},f_{S}\right):=\max\left\{\max_{2\leq\ell\leq 5}% \left\{\max_{t\in\left[t_{n-1/2},t_{n+2}\right]}\left\|\partial_{t}^{\ell}u_{S% }(t)\right\|_{\mathcal{T}}\right\},\max_{0\leq\ell\leq 3}\left\{\max_{t\in% \left[t_{n-1},t_{n+2}\right]}\left\|\partial_{t}^{\ell}f_{S}(t)\right\|_{% \mathcal{T}}\right\}\right\}.

We apply Lemma 5 to (3.51) to obtain

\Delta t\left\|\mbox{\boldmath$\sigma$}^{\left(n\right)}\right\|_{\ell^{1}}% \leq\dfrac{\Delta t^{2}}{24}\left(24\dfrac{\nu+1}{\nu}+\dfrac{1+C_{2}}{c_{\nu}% }+1\right)\mathcal{M}_{n}\left(u_{S},f_{S}\right).

(3.52)

Next, we investigate the summands in the second term of the right-hand side of (3.50),

	$\displaystyle\mbox{\boldmath$\sigma$}^{\left(m\right)}-\mbox{\boldmath$\sigma$% }^{\left(m+1\right)}=$
	$\displaystyle\!\!\left[\!\!\!\begin{array}[c]{c}\dfrac{u_{S}\left(t_{m+2}% \right)-2u_{S}\left(t_{m+1}\right)+u_{S}\left(t_{m}\right)}{\Delta t^{2}}+% \dfrac{v_{S}\left(t_{m+1/2}\right)-v_{S}\left(t_{m+3/2}\right)}{\Delta t}\\ \dfrac{u_{S}\left(t_{m}\right)-u_{S}\left(t_{m+1}\right)}{\Delta t}-\left(A^{S% ,p,\nu}\right)^{-1}\!\left(\dfrac{v_{S}\left(t_{m+3/2}\right)-2v_{S}\left(t_{m% +1/2}\right)+v_{S}\left(t_{m-1/2}\right)}{\Delta t^{2}}+\dfrac{R_{S}^{(m)}-R_{% S}^{(m+1)}}{\Delta t}\right)\end{array}\!\!\!\right],$

where $R_{S}^{(m)}$ is defined as in Corollary 4. By similar arguments as in the proof of [20, Theorem 3.10] and by Lemma 5, we obtain

	$\displaystyle\left\\|\mbox{\boldmath$\sigma$}^{\left(m\right)}-\mbox{\boldmath$% \sigma$}^{\left(m+1\right)}\right\\|_{\ell^{1}}$	$\displaystyle\leq\dfrac{\Delta t^{2}}{24}\left(6\dfrac{\nu+1}{\nu}+\dfrac{3}{c% _{\nu}}+4\right)\mathcal{M}_{m}\left(u_{S},f_{S}\right)$
		$\displaystyle\qquad+\left\\|\left(A^{S,p,\nu}\right)^{-1}\left(\dot{f}_{S}^{(m+% 1/2)}-\dot{R}_{S}^{(m+1/2)}\right)\right\\|_{\mathcal{T}}.$		(3.53)

We apply Lemma5 to (3.53) to obtain

\left\|\mbox{\boldmath$\sigma$}^{\left(m\right)}-\mbox{\boldmath$\sigma$}^{% \left(m+1\right)}\right\|_{\ell^{1}}\leq\dfrac{\Delta t^{2}}{24}\left(6\dfrac{% \nu+1}{\nu}+\dfrac{3+C_{2}}{c_{\nu}}+4\right)\mathcal{M}_{m}\left(u_{S},f_{S}\right)

(3.54)

The rest of the proof follows by the same arguments as in the proof of [20, Theorem 3.10].

Theorem 7

Assume that (2.1), (2.5), (2.6), (3.21) and (3.22) hold and let the solution and right-hand side of (2.2) satisfy $u\in W^{8,\infty}\left(\left[0,T\right];H^{m+1}\left(\Omega\right)\right)$ and $f\in W^{6,\infty}\left(\left[0,T\right];H^{m+1}\left(\Omega\right)\right)$ , respectively. Then, the corresponding fully discrete Galerkin FE formulation with local time-stepping (3.20) has a unique solution which satisfies the error estimate

\left\|u\left(t_{n+1}\right)-u_{S}^{(n+1)}\right\|\leq C(1+T)\left(h^{m+1}+% \Delta t^{2}\right)\mathcal{Q}(u,f)

with

Q(u,f):=\left(1+Ch^{m+1}\left(1+T\right)\right)\max\left\{\left\|u\right\|_{W^% {8,\infty}\left(\left[0,T\right];H^{m+1}\left(\Omega\right)\right)},\left\|f% \right\|_{W^{6,\infty}\left(\left[0,T\right];H^{m+1}\left(\Omega\right)\right)% }\right\}

and constants $C$ independent of $n$ , $\Delta t$ , $h$ , $p$ , and the final time $T$ .

Proof. We adapt the proof in [20, Thm .3.11] to the case of an inhomogeneous right-hand side. We again use [1, Thm 4.1] to estimate the error between the exact and semi-discrete solution. All assumptions in that theorem are verified in the proof of [20, Thm. 3.11] except for the condition for the right-hand side. As usual, we denote below by $C$ a generic constant which may differ from one instance to the next.

From [7, Lem. 5.2] it follows that

\left|\left(\chi,\psi\right)-\left(\chi,\psi\right)_{\mathcal{T}}\right|\leq Ch% ^{r+\mu}\left(\sum_{\tau\in\mathcal{T}}\left\|\chi\right\|_{H^{r}\left(\tau% \right)}^{2}\right)^{1/2}\left(\sum_{\tau\in\mathcal{T}}\left\|\psi\right\|_{H% ^{\mu}\left(\tau\right)}^{2}\right)^{1/2}\quad\forall\chi,\psi\in S

(3.55)

for all $1\leq r,\mu\leq m-1$ . To deduce condition (3.5) in [1] we transform the local quantities $\left\|\chi\right\|_{H^{r}\left(\tau\right)}^{2}$ in (3.55) to the unit reference simplex $\hat{\tau}$ . By $F_{\tau}:\hat{\tau}\rightarrow\tau$ , we denote an affine bijection from the reference element to the simplex $\tau$ . For a function $g$ with domain $\tau,$ its pullback is denoted by $\hat{g}_{\tau}=g\circ F_{\tau}$ . Let $\left|\tau\right|$ denote the volume of the simplex $\tau\in\mathcal{T}$ . A standard scaling argument implies:

\left\|\chi\right\|_{H^{r}\left(\tau\right)}^{2}\leq C\left|\tau\right|\sum_{% \left|\boldsymbol{\alpha}\right|\leq r}h_{\tau}^{-2\left|\boldsymbol{\alpha}% \right|}\left\|\partial^{\boldsymbol{\alpha}}\hat{\chi}\right\|_{L^{2}\left(% \hat{\tau}\right)}^{2}.

For $\mu\in\mathbb{N}$ and $q\geq 2$ , the embedding $W^{\mu,q}\left(\hat{\tau}\right)\hookrightarrow L^{2}\left(\hat{\tau}\right)$ is continuous so that

	$\displaystyle\left\\|\chi\right\\|_{H^{r}\left(\tau\right)}^{2}$	$\displaystyle\leq C\left\|\tau\right\|\sum_{\left\|\boldsymbol{\alpha}\right\|\leq r% }h_{\tau}^{-2\left\|\boldsymbol{\alpha}\right\|}\left\\|\partial^{\boldsymbol{% \alpha}}\hat{\chi}\right\\|_{W^{\mu,q}\left(\hat{\tau}\right)}^{2}$
		$\displaystyle\leq C\left\|\tau\right\|\sum_{\left\|\boldsymbol{\alpha}\right\|\leq r% }h_{\tau}^{-2\left\|\boldsymbol{\alpha}\right\|}\sum_{\left\|\boldsymbol{\beta}% \right\|\leq\mu}\left\\|\partial^{\boldsymbol{\alpha}+\boldsymbol{\beta}}\hat{% \chi}\right\\|_{L^{q}\left(\hat{\tau}\right)}^{2}.$

Transforming back to the physical simplex $\tau$ yields

\left\|\chi\right\|_{H^{r}\left(\tau\right)}^{2}\leq C\left|\tau\right|\sum_{% \left|\boldsymbol{\alpha}\right|\leq r}h_{\tau}^{-2\left|\boldsymbol{\alpha}% \right|}\sum_{\left|\boldsymbol{\beta}\right|\leq\mu}h^{2\left|\boldsymbol{% \alpha}+\boldsymbol{\beta}\right|}\left|\tau\right|^{-1/q}\left\|\partial^{% \boldsymbol{\alpha}+\boldsymbol{\beta}}\chi\right\|_{L^{q}\left(\tau\right)}^{% 2}.

We use

h_{\tau}^{-2\left|\boldsymbol{\alpha}\right|}h_{\tau}^{2\left|\boldsymbol{% \alpha}+\boldsymbol{\beta}\right|}\leq 1\quad\text{and\quad}\left|\tau\right|% \left|\tau\right|^{-2/q}\leq 1

to obtain

\left\|\chi\right\|_{H^{r}\left(\tau\right)}^{2}\leq C\sum_{\left|\boldsymbol{% \alpha}\right|\leq r}\sum_{\left|\boldsymbol{\beta}\right|\leq\mu}\left\|% \partial^{\boldsymbol{\alpha}+\boldsymbol{\beta}}\chi\right\|_{L^{q}\left(\tau% \right)}^{2}\leq C\left\|\chi\right\|_{W^{r+\mu,q}\left(\tau\right)}^{2}.

(3.56)

Since the $\left\|\cdot\right\|_{\ell^{2}}-$ and the $\left\|\cdot\right\|_{\ell^{q}}-$ norm are equivalent in finite-dimensional Euclidean (index) spaces, the combination of (3.55) with (3.56) results in

\left|\left(\chi,\psi\right)-\left(\chi,\psi\right)_{\mathcal{T}}\right|\leq Ch% ^{r+\mu}\left(\sum_{\tau\in\mathcal{T}}\left\|\chi\right\|_{W^{r+\mu,q}\left(% \tau\right)}^{q}\right)^{1/q}\left(\sum_{\tau\in\mathcal{T}}\left\|\psi\right% \|_{H^{\mu}\left(\tau\right)}^{2}\right)^{1/2}\quad\forall\chi,\psi\in S

(3.57)

and this implies condition (3.5) in [1, Thm 4.1].

To estimate the effect of mass lumping in the right-hand side, let $f:=f\left(\cdot,t\right)\in W^{r+\mu,q}\left(\Omega\right)$ for $0<t\leq T$ and $\mu\in\left\{1,2\right\}$ . We employ the quasi-interpolation operator $\mathcal{I}_{h}^{\operatorname*{av}}:L^{1}\left(\Omega\right)\rightarrow S$ as in [17, (5.1)] and the splitting:

\left(f,\psi\right)-\left(f,\psi\right)_{\mathcal{T}}=\left(f-\mathcal{I}_{h}^% {\operatorname*{av}}f,\psi\right)-\left(f-\mathcal{I}_{h}^{\operatorname*{av}}% f,\psi\right)_{\mathcal{T}}+\left(\mathcal{I}_{h}^{\operatorname*{av}}f,\psi% \right)-\left(\mathcal{I}_{h}^{\operatorname*{av}}f,\psi\right)_{\mathcal{T}}.

(3.58)

For the last difference we use (3.57) to obtain

\left|\left(\mathcal{I}_{h}^{\operatorname*{av}}f,\psi\right)-\left(\mathcal{I% }_{h}^{\operatorname*{av}}f,\psi\right)_{\mathcal{T}}\right|\leq Ch^{r+\mu}% \left(\sum_{\tau\in\mathcal{T}}\left\|\mathcal{I}_{h}^{\operatorname*{av}}f% \right\|_{W^{r+\mu,q}\left(\tau\right)}^{q}\right)^{1/q}\left(\sum_{\tau\in% \mathcal{T}}\left\|\psi\right\|_{H^{\mu}\left(\tau\right)}^{2}\right)^{1/2}.

(3.59)

For the other terms in (3.58), the continuity of $\left(\cdot,\cdot\right)$ and $\left(\cdot,\cdot\right)_{\mathcal{T}}$ yields

\left|\left(f-\mathcal{I}_{h}^{\operatorname*{av}}f,\psi\right)\right|+\left|% \left(f-\mathcal{I}_{h}^{\operatorname*{av}}f,\psi\right)_{\mathcal{T}}\right|% \overset{\eqref{ceqCeq}}{\leq}\left(1+C_{\operatorname*{eq}}^{2}\right)\left\|% f-\mathcal{I}_{h}^{\operatorname*{av}}f\right\|\left\|\psi\right\|.

(3.60)

The approximation and stability properties of $\mathcal{I}_{h}^{\operatorname*{av}}$ are studied in [17, Lem. 5.1, Thm. 5.2] and lead to

\left\|f-\mathcal{I}_{h}^{\operatorname*{av}}f\right\|\leq Ch^{r+\mu}\left\|f% \right\|_{H^{r+\mu}\left(\Omega\right)}

(3.61)

since $r+\mu\leq m+1$ and

\left\|\mathcal{I}_{h}^{\operatorname*{av}}f\right\|_{W^{r+\mu,q}\left(\tau% \right)}\leq C\left\|f\right\|_{W^{r+\mu,q}\left(\tau\right)}.

(3.62)

We combine (3.58)–(3.62) to deduce the estimate

\left|\left(f,\psi\right)-\left(f,\psi\right)_{\mathcal{T}}\right|\leq Ch^{r+% \mu}\left(\left\|f\right\|_{H^{r+\mu}\left(\Omega\right)}\left\|\psi\right\|+% Ch^{r+\mu}\left\|f\right\|_{W^{r+\mu,q}\left(\Omega\right)}\left(\sum_{\tau\in% \mathcal{T}}\left\|\psi\right\|_{H^{\mu}\left(\tau\right)}^{2}\right)^{1/2}% \right).

Sobolev’s embedding theorem implies $W^{r+\mu,q}\left(\Omega\right)\hookrightarrow H^{r+\mu}\left(\Omega\right)$ and $\left\|\psi\right\|\leq\sqrt{\sum_{\tau\in\mathcal{T}}\left\|\psi\right\|_{H^{% \mu}\left(\tau\right)}^{2}}$ for $\mu\in\left\{1,2\right\}$ and, finally, (3.8) in [1, Thm 4.1] follows. The application of this theorem yields the error estimate between the exact and semi-discrete solution

	$\displaystyle\left\\|\partial_{t}^{\ell}\left(u-u_{S}\right)\right\\|_{L^{\infty% }\left(\left[0,T\right];L^{2}\left(\Omega\right)\right)}$	$\displaystyle\leq C_{\ell}h^{m+1}\left(\sum_{k=\ell}^{\ell+1}\left\\|\partial_{% t}^{k}f\right\\|_{L^{2}\left(\left[0,T\right];H^{m+1}\left(\Omega\right)\right)% }\right.$
		$\displaystyle\qquad\left.+\sum_{k=\ell}^{\ell+3}\left\\|\partial_{t}^{k}u\right% \\|_{L^{2}\left(\left[0,T\right];H^{m+1}\left(\Omega\right)\right)}\right)$		(3.63)

for $\ell=0$ . Inspection of the proof in [1, Thm 4.1] implies that (3.63) also holds for any $\ell\in\mathbb{N}$ provided the right-hand side in (3.63) exists. Proceeding as in the proof of [20, Thm. 3.11] we arrive at

	$\displaystyle\left\\|\partial_{t}^{\ell}\left(u-u_{S}\right)\right\\|_{L^{\infty% }\left(\left[0,T\right];L^{2}\left(\Omega\right)\right)}$	$\displaystyle\leq C_{\ell}h^{m+1}\left(1+T\right)\left(\left\\|f\right\\|_{W^{% \ell+1,\infty}\left(\left[0,T\right];H^{m+1}\left(\Omega\right)\right)}\right.$
		$\displaystyle\qquad\left.+\left\\|u\right\\|_{W^{\ell+3,\infty}\left(\left[0,T% \right];H^{m+1}\left(\Omega\right)\right)}\right).$		(3.64)

As in the proof of Theorem [20, Thm. 3.11], the combination with (3.64) leads to

	$\displaystyle\left\\|u\left(t_{n+1}\right)-u_{S}^{\left(n+1\right)}\right\\|$	$\displaystyle\leq C\left(1+T\right)\left(h^{m+1}+\Delta t^{2}\right)\max\left% \{\mathcal{M}\left(u_{S},f_{S}\right),\right.$
		$\displaystyle\quad\left.\left\\|u\right\\|_{W^{3,\infty}\left(\left[0,T\right];H% ^{m+1}\left(\Omega\right)\right)},\left\\|f\right\\|_{W^{1,\infty}\left(\left[0,% T\right];H^{m+1}\left(\Omega\right)\right)}\right\}.$		(3.65)

The final estimate is obtained by recalling the definition of $\mathcal{M}$ as in Thm. 6:

	$\displaystyle\left\\|\partial_{t}^{\ell}u_{S}\right\\|_{L^{\infty}\left(\left[0,% T\right];L^{2}\left(\Omega\right)\right)}\leq$	$\displaystyle\left\\|\partial_{t}^{\ell}u\right\\|_{L^{\infty}\left(\left[0,T% \right];L^{2}\left(\Omega\right)\right)}+\left\\|\partial_{t}^{\ell}\left(u_{S}% -u\right)\right\\|_{L^{\infty}\left(\left[0,T\right];H^{m+1}\left(\Omega\right)% \right)}$
	$\displaystyle\overset{\text{(\ref{deltumus})}}{\leq}$	$\displaystyle\left\\|\partial_{t}^{\ell}u\right\\|_{L^{\infty}\left(\left[0,T% \right];L^{2}\left(\Omega\right)\right)}+$
		$\displaystyle+C_{\ell}h^{m+1}\left(1+T\right)\max\left\{\left\\|u\right\\|_{W^{% \ell+3,\infty}\left(\left[0,T\right];H^{m+1}\left(\Omega\right)\right)},\right.$
		$\displaystyle\quad\left.\left\\|f\right\\|_{W^{\ell+1,\infty}\left(\left[0,T% \right];H^{m+1}\left(\Omega\right)\right)}\right\}.$

Since $f\in W^{3,\infty}\left(\left[0,T\right],L^{2}\left(\Omega\right)\right)$ , we obtain from the relation

\left(\partial_{t}^{\ell}f_{S}\left(t\right),v\right)=\left(\partial_{t}^{\ell% }f\left(t\right),v\right)\quad\forall v\in S

and upon setting $v=\partial_{t}^{\ell}f_{S}\left(t\right)\in S$ :

\left\|\partial_{t}^{\ell}f_{S}\right\|_{L^{\infty}\left(\left[0,T\right];L^{2% }\left(\Omega\right)\right)}\leq\left\|\partial_{t}^{\ell}f\right\|_{L^{\infty% }\left(\left[0,T\right];L^{2}\left(\Omega\right)\right)}.

Hence,

\mathcal{M}\left(u_{S},f_{S}\right)\leq\left(1+Ch^{m+1}\left(1+T\right)\right)% \max\left\{\left\|u\right\|_{W^{8,\infty}\left(\left[0,T\right];H^{m+1}\left(% \Omega\right)\right)},\left\|f\right\|_{W^{6,\infty}\left(\left[0,T\right];H^{% m+1}\left(\Omega\right)\right)}\right\}

and the combination with (3.65) leads to the assertion.

4 Towards Optimal $H^{1}$ Convergence

The LF-LTS( $\nu$ ) Algorithm 1 achieves the optimal convergence rate ${\cal O}(\Delta t^{2}+h^{m+1})$ with respect to the $L^{2}$ -norm with piecewise polynomial FE of degree $m$ also in the presence of a source $f$ – see Theorem 7. In [6], however, it was shown that the original LF-LTS method from [21, Algorithm 1], that is, the LF-LTS( $\nu$ ) Algorithm 1 with $\nu=0$ , may not achieve the expected optimal convergence rate in the $H^{1}$ -norm whenever the source $f$ is nonzero across the coarse-to-fine interface. To restore optimal $H^{1}$ -convergence regardless of $f$ , we shall replace the somewhat abrupt projections from Section 2.2 by a new more gradual weighted transition.

Let the set of simplex vertices be denoted by $\mathcal{V}\left(\mathcal{T}\right)$ and let $\mathcal{E}\left(\mathcal{T}\right)$ denote the set of edges if $d\geq 2$ while for $d=1$ we set $\mathcal{E}\left(\mathcal{T}\right)=\mathcal{T}$ . The definition of $S_{\mathcal{T}}^{m}$ gives rise to the interpolation operator $I_{\mathcal{T}}^{m}:C^{0}\left(\overline{\Omega}\right)\rightarrow S_{\mathcal% {T}}^{m}$ . For example, for $m=1$ the space $S_{\mathcal{T}}^{1}$ consists of the continuous, piecewise linear polynomials and $I_{\mathcal{T}}^{1}$ is the nodal interpolant

I_{\mathcal{T}}^{1}u=\sum_{z\in\mathcal{V}\left(\mathcal{T}\right)}u\left(z% \right)b_{1,z}\text{.}

For two points $y,z$ in Euclidean space, we denote by $\left[y,z\right]$ the straight line connecting $y$ and $z$ .

Definition 1 (connecting sequence)

For $y,z\in\mathcal{V}\left(\mathcal{T}\right)$ , a sequence $\mathbf{x}=\left(x_{i}\right)_{i=0}^{m}\subset\mathcal{V}\left(\mathcal{T}\right)$ is a connecting sequence of $y$ and $z$ if

x_{0}=y,\quad x_{m}=z,\quad\text{and for }1\leq j\leq m:\text{ }\left[x_{j-1},% x_{j}\right]\in\mathcal{E}\left(\mathcal{T}\right)\text{.}

The length of $\mathbf{x}$ is $\operatorname{length}\left(\mathbf{x}\right):=m$ .

Definition 2 (graph distance)

The graph distance between $y,z\in\mathcal{V}\left(\mathcal{T}\right)$ is given by

\operatorname{dist}_{\mathcal{T}}\left(y,z\right):=\min\left\{\operatorname{% length}\left(\mathbf{x}\right):\mathbf{x}\text{ is a connecting sequence of }y\text{ and }z\right\}.

The graph distance between $\mathcal{T}_{\operatorname{f}}$ and some $z\in\mathcal{V}\left(\mathcal{T}\right)$ is given by

\operatorname{dist}\left(\mathcal{T}_{\operatorname{f}},z\right):=\min_{y\in% \mathcal{T}_{\operatorname{f}}}\operatorname{dist}_{\mathcal{T}}\left(y,z% \right).

Definition 3 (discrete distance function)

Let $\mathcal{T}$ be a conforming simplicial finite element mesh and let the subset $\mathcal{T}_{\operatorname{f}}\subset\mathcal{T}$ be fixed. The discrete distance function for $s$ simplex layers about $\mathcal{T}_{\operatorname{f}}$ , $s\in\mathbb{N}_{\geq 1}$ , is given as follows. For $z\in\mathcal{V}$ , its value is set to

\eta_{s}\left(z\right):=\left(1-\frac{\operatorname{dist}\left(\mathcal{T}_{% \operatorname{f}},z\right)}{s}\right)_{+}

with $\left(x\right)_{+}:=\max\left\{x,0\right\}$ . On $\Omega$ , the function is then given by

\eta_{s}=\sum_{z\in\mathcal{V}\left(\mathcal{T}\right)}\eta_{s}\left(z\right)b% _{1,z}.

Definition 4 (fine/coarse mesh mapping)

The fine mesh mapping $\Pi_{\operatorname{f}}^{S}:S\rightarrow S$ and coarse mesh mapping $\Pi_{\operatorname{c}}^{S}:S\rightarrow S$ are given by

\Pi_{\operatorname{f}}^{S}u:=\sum_{z\in\Sigma}u_{z}\eta_{s}\left(z\right)b_{z}% \quad\text{and\quad}\Pi_{\operatorname{c}}^{S}:=I-\Pi_{\operatorname{f}}^{S}.

(4.1)

Remark 4

Note that the above definition of $\Pi_{\operatorname{f}}^{S}$ and $\Pi_{\operatorname{c}}^{S}$ generalizes the respective operators as defined in [20] and (2.11): For $s=1$ , both definitions coincide whereas for $s>1$ the coarse-to-fine transition is now spread across several layers.

Algorithm 1 formally remains unchanged; however, the operators $\Pi_{\operatorname{f}}^{S}$ and $\Pi_{\operatorname{c}}^{S}$ have to replaced by the new definition (4.1). This comes at an extra cost since the local time stepping is applied in $s$ additional simplex layers around the original fine region $\Omega_{\operatorname*{f}}$ as in (2.9) which are recursively defined by $\Omega_{\operatorname*{f}}^{+1}$ as in (2.9) and for $r=2,\ldots s$ by

\Omega_{\operatorname*{f}}^{+r}:=\operatorname*{int}\left(\bigcup{\left\{\tau% \in\mathcal{T}:\tau\cap\overline{\Omega_{\operatorname*{f}}^{+\left(r-1\right)% }}\neq\emptyset\right\}}\right).

In Section 5.2 we will report on numerical experiments for this modification.

5 Numerical Experiments

First, we present numerical experiments which corroborate the convergence theory from Section 3 and hence demonstrate that the LF-LTS( $\nu$ ) method for the inhomogeneous wave equation defined in Algorithm 1 achieves the expected optimal rates of convergence with respect to the $L^{2}$ norm. In doing so, we also compare the LF-LTS method with the alternative “split-LFC” approach from [3], where the source inside the locally refined region is only evaluated once at each global time-step – see Remark 1.

Next, in Section 5.2, we show via two numerical experiments that the weighted transition, introduced in Section 4, indeed restores optimal $H^{1}$ -convergence regardless of the source and its position with respect to the refined region in the mesh.

5.1 Optimal $L^{2}$ Convergence

We consider the one-dimensional wave equation (2.4) in $\Omega=(0,4)$ with homogeneous Dirichlet boundary conditions, i.e. $\Gamma=\Gamma_{D}$ , $c\equiv 1$ , and inhomogeneous source term

f(x,t)=250\,e^{-400\left(\left(x-2\right)^{2}+\left(t-0.1\right)^{2}\right)},

(5.1)

a narrow space-time Gaussian pulse centered about $x=2$ and $t=0.1$ . The initial conditions are set to zero, i.e., $u_{0}\equiv 0$ , $v_{0}\equiv 0$ .

For the numerical solution, we use piecewise linear $H^{1}$ -conforming finite elements with mass-lumping in space and the LF-LTS( $\nu$ ) method with global time-step $\Delta t=\operatorname*{e}^{-\nu}h_{\operatorname*{c}}$ and $\nu=0.01$ in time. Next, we split the computational domain $\Omega$ into a coarse part, $\Omega_{\operatorname*{c}}=(0,1.6)\cup(2.4,4)$ , and a locally refined part, $\Omega_{\operatorname*{f}}=[1.6,2.4]$ , inside each of which we use an equidistant mesh with mesh sizes $h_{\operatorname*{c}}$ or $h_{\operatorname*{f}}=h_{\operatorname*{c}}/2$ , respectively. Hence inside $\Omega_{\operatorname*{f}}$ , the LF-LTS( $\nu$ ) method takes $p=2$ local time-steps of size $\Delta t/p$ for each global time-step of size $\Delta t$ .

In Fig. 1, we display the LF-LTS( $\nu$ ) numerical solution for $h_{\operatorname*{c}}=0.01$ at times $t\approx 0.05$ , $t\approx 0.1$ , and $t=0.15$ , together with the relative $L^{2}$ -error at $t=0.15$ for a sequence of locally refined meshes. The $L^{2}$ -error indeed displays the expected optimal second-order convergence rate $\mathcal{O}\left(h_{\operatorname*{c}}^{2}\right)$ , as proved in Theorem 7.

Refer to caption — Figure 1: Snapshots of the LF-LTS( $\nu$ ) solution of (2.4) with $\nu=0.01$ and $f$ as in (5.1) (top left, top right and bottom left). Relative $L^{2}$ -error vs. $h=h_{\operatorname*{c}}$ for $\mathbb{P}_{1}$ finite elements (bottom right).

In Fig. 2, that same relative $L^{2}$ -error is compared to that of the alternative approach from [3], denoted here by LFC-LTS, which omits the intermediate source evaluations inside the refined region. Although both methods converge with the optimal rate $\mathcal{O}\left(h_{\operatorname*{c}}^{2}\right)$ , as expected, the LF-LTS( $\nu$ ) algorithm from Section 2.3 here yields a threefold reduction in the error.

5.2 Optimal $H^{1}$ -Convergence and weighted transition

In [6], it was shown for a particular (spatially constant) source term $f$ that the original LF-LTS method from [21, Algorithm 1], that is, the LF-LTS( $\nu$ ) Algorithm 1 with $\nu=0$ , may not always achieve the expected optimal convergence rate in the $H^{1}$ -norm. Here, we carefully study this behavior for two different source terms and show that the weighted transition from Section 4 restores the optimal $H^{1}$ -convergence regardless of $f$ .

5.2.1 Space-time local forcing

Again, we consider the wave equation (2.4) with zero initial conditions and source $f$ as in (5.1) from the previous section. For spatial discretization, however, we use piecewise quadratic $H^{1}$ -conforming finite elements with mass-lumping ( $\mathbb{P}_{2}$ -FE). The computational domain separates into a coarse part $\Omega_{\operatorname*{c}}$ and a locally refined part $\Omega_{\operatorname*{f}}$ with corresponding mesh sizes $h_{\operatorname*{c}}$ or $h_{\operatorname*{f}}=h_{\operatorname*{c}}/5$ , respectively.

Now, while keeping the source fixed, we shall progressively shift the location of the refined region rightward and thus mimic the following three distinct situations: (i) source $f$ located entirely inside $\Omega_{\operatorname*{f}}$ , (ii) source $f$ nonzero across coarse-to-fine mesh interface, (iii) source $f$ located entirely inside $\Omega_{\operatorname*{c}}$ . In each case we monitor the $H^{1}$ -convergence rates of the LF-LTS( $\nu$ ) method for $p=5$ and $\nu=0.01$ either with, or without, weighted transition. Hence the fine mesh mapping $\Pi^{S}_{\operatorname*{f}}$ is either defined by (4.1) with weighted transition across $s\sim 1/h_{\operatorname*{c}}$ elements, or without weighted transition, i.e., as in (4.1) with $s=1$ – see also Remark 4.

First, we set $\Omega_{\operatorname*{f}}=[1.6,2.4]$ so that the support of $f$ completely lies inside $\Omega_{\operatorname*{f}}$ (up to machine precision), as illustrated in Fig. 3. Here, the LF-LTS( $\nu$ ) method achieves the optimal second-order rate $\mathcal{O}(h_{\operatorname*{c}}^{2})$ both with, or without, weighted transition, as shown in Fig. 3. Recall that overall second-order convergence vs. $h_{\operatorname*{c}}$ (or $\Delta t$ ) in the $H^{1}$ -norm is expected here since we use $\mathbb{P}_{2}$ -FE in space.

Next, we shift the locally refined part rightwards and set $\Omega_{\operatorname*{f}}=[2,2.4]$ so that $f$ is non-zero across the interface between $\Omega_{\operatorname*{f}}$ and $\Omega_{\operatorname*{c}}$ , as shown in Fig. 4; in fact, the coarse-to-fine mesh interface now coincides with the maximum of $f$ at $x=2$ . In Fig. 4, we observe that the LF-LTS( $\nu$ ) method without weighted transition only converges as $\mathcal{O}(h_{\operatorname*{c}}^{1.5})$ , while using the weighted transition across $s\sim 1/h_{\operatorname*{c}}$ layers of elements restores the expected optimal second-order convergence.

Finally, we shift the locally refined part even farther to the right and set $\Omega_{\operatorname*{f}}:=[2.2,2.4]$ so that the support of the source term essentially lies inside the coarse part of the mesh, as shown in Fig. 5. Here, as in the first case, we again observe optimal second-order convergence in the $H^{1}$ -norm, be it with or without weighted transition.

In summary, the LF-LTS( $\nu$ ) method without weighted transition achieves in the considered examples the expected optimal convergence rates in the $H^{1}$ -norm when the source vanishes at the coarse-to-fine mesh transition; otherwise, when the source is non-zero across the mesh interface, the $H^{1}$ -convergence rate is reduced by one half. The weighted transition introduced in Section 4 restores the expected optimal $H^{1}$ -convergence even when $f$ is nonzero across the coarse-to-fine mesh interface.

5.2.2 Spatially constant solution

Finally, we revisit the example from [6, Section 7.1.3] and consider the wave equation (2.4) in $\Omega=(0,1)$ with homogeneous Neumann boundary conditions and zero initial conditions. Hence, we set the right-hand side $f$ as in [6] such that the solution is given by the spatially constant function

u(x,t)=\left\{\begin{aligned} &0,&&t\leq 0.1,\\ &\frac{1}{1+\operatorname*{e}^{0.8\left(\frac{1}{t-0.1}+\frac{1}{t-0.9}\right)% }},&&0.1<t<0.9,\\ &1,&&t\geq 0.9.\end{aligned}\right.

(5.2)

Now, we discretize (2.4) with $\mathbb{P}_{2}$ -FE in space and the LF-LTS( $\nu$ ) Algorithm 1 with $\nu=0.01$ in time. The refined region is (arbitrarily) set to the right half of the computational domain as $\Omega_{\operatorname*{f}}=[0.5,1]$ with a coarse-to-fine mesh size ratio of $p=2$ . Clearly, here the spatially constant source is nonzero across the coarse-to-fine mesh transition.

In Figure 6, we compare the $L^{2}$ - and $H^{1}$ -errors on a sequence of meshes of size $h=h_{\operatorname*{c}}$ either with, or without, weighted transition. Although the $L^{2}$ -convergence theory from Section 3 does not apply here due to the Neumann boundary conditions, we still observe optimal second-order convergence both with or without weighted transition. Note that the overall second-order rate of convergence is dictated here by the time accuracy of the LF-LTS scheme. In contrast for the relative $H^{1}$ -error, we observe an order reduction by $h^{1/2}$ without weighted transition. By using the weighted transition across $s\sim 1/h$ element layers from Section 4, however, optimal second-order convergence is again restored – recall that $\mathbb{P}_{2}$ -FE are used here for the spatial discretization.

To exhibit the mechanism that leads to this subtle order reduction by $h^{1/2}$ , which only affects the $H^{1}$ and not the $L^{2}$ -convergence, we now strongly magnify in Fig. 7 the two (nearly constant) numerical solutions across the coarse-to-fine mesh transition for two different mesh sizes. Both variants of the LF-LTS( $\nu$ ) method, with or without weighted transition, introduce a small perturbation near $x=0.5$ , whose amplitude tends to zero as $h\rightarrow 0$ . Without weighted transition, however, this perturbation increasingly steepens as $h\to 0$ , which results in increased gradients at the intersection of $\Omega_{\operatorname*{c}}$ and $\Omega_{\operatorname*{f}}$ . Hence the weighted transition prevents the formation of increasingly steep gradients as $h\rightarrow 0$ at the coarse-to-fine mesh meeting point and thus restores the expected optimal $H^{1}$ -convergence.

Appendix A Appendix

Here we prove some technical estimates about Chebyshev polynomials.

Lemma 8

Let $\delta_{p,\nu}$ be as in (2.12) and $0\leq m\leq n\leq p-1$ . Then

\left|U_{n}^{\left(m\right)}\left(\delta_{p,\nu}\right)\right|\leq\frac{\left(% n+1\right)^{2m+1}\operatorname*{e}^{\nu/2}}{\left(2m+1\right)!!}.

(A.1)

An estimate from below is given by

U_{n}^{\left(m\right)}\left(\delta_{p,\nu}\right)\geq 2^{m}m!\binom{n+m+1}{n-m% }\left(\frac{2\nu}{p^{2}}\right)^{m}

with the particular case

U_{n}\left(\delta_{p,\nu}\right)\geq n+1.

Proof. From [15, 18.9.21], we know that $U_{n}^{\left(m\right)}=\frac{1}{n+1}T_{n+1}^{\left(m+1\right)}$ and (A.1) is equivalent to

T_{n+1}^{\left(m+1\right)}\left(\delta_{p,\nu}\right)\leq\frac{\left(n+1\right% )^{2m+2}\operatorname*{e}^{\nu/2}}{\left(2m+1\right)!!}.

The definition $\delta_{p,\nu}=1+\nu/p^{2}$ along the side constraints for $m,n,p-1$ imply

\delta_{p,\nu}\leq\delta_{n+1,\nu}.

Since $T_{n+1}^{\left(m+1\right)}\left(x\right)$ is increasing for $x\geq 1$ we obtain from [20, Lem. A.1]

T_{n+1}^{\left(m+1\right)}\left(\delta_{p,\nu}\right)\leq T_{n+1}^{\left(m+1% \right)}\left(\delta_{n+1,\nu}\right)\leq\frac{\left(n+1\right)^{2\left(m+1% \right)}\operatorname*{e}^{\nu/2}}{\left(2m+1\right)!!}.

For the estimate from below we use the representation in [15, 18.7.4 and 18.5.7]

	$\displaystyle U_{n}^{\left(m\right)}\left(x\right)$	$\displaystyle=\frac{m!}{2^{m}}\frac{n!}{\left(\frac{3}{2}\right)_{n}}\sum_{% \ell=0}^{n}\frac{{\left(n+2\right)_{\ell}}{\left(\ell+3/2\right)_{n-\ell}}}{% \ell!(n-\ell)!}\binom{\ell}{m}\left(\frac{x-1}{2}\right)^{\ell-m}$
		$\displaystyle=2^{m}m!\sum_{\ell=0}^{n}\binom{n+\ell+1}{n-\ell}\binom{\ell}{m}% \left(2\left(x-1\right)\right)^{\ell}.$

For $x=\delta_{p,\nu}$ , this is a sum of positive terms and hence a lower bound is obtained by truncating it after the first non-zero term, i.e., for $\ell=m:$

U_{n}^{\left(m\right)}\left(\delta_{p,\nu}\right)\geq 2^{m}m!\binom{n+m+1}{n-m% }\left(\frac{2\nu}{p^{2}}\right)^{m}.

References

[1] G. A. Baker and V. A. Dougalis. The effect of quadrature errors on finite element approximations for second order hyperbolic equations. SIAM Journal on Numerical Analysis, 13(4):577–598, 1976.
[2] C. Carle and M. Hochbruck. Error analysis of multirate leapfrog-type methods for second-order semilinear ODEs. SIAM J. Numer. Anal., 60(5):2897–2924, 2022.
[3] C. Carle and M. Hochbruck. Error analysis of second-order local time integration methods for discontinuous Galerkin discretizations of linear wave equations. Math. Comp., 93(350):2611–2641, 2024.
[4] C. Carle, M. Hochbruck, and A. Sturm. On leapfrog-Chebyshev schemes. SIAM J. Numer. Anal., 58(4):2404–2433, 2020.
[5] J. Chabassier and S. Imperiale. Fourth-order energy-preserving locally implicit time discretization for linear wave equations. Internat. J. Numer. Methods Engrg., 106(8):593–622, 2016.
[6] Chabassier, Juliette and Imperiale, Sébastien. Construction and convergence analysis of conservative second order local time discretisation for linear wave equations. ESAIM: M2AN, 55(4):1507–1543, 2021.
[7] G. Cohen, P. Joly, J. E. Roberts, and N. Tordjman. Higher order triangular finite elements with mass lumping for the wave equation. SIAM J. Numer. Anal., 38(6):2047–2078, 2001.
[8] G. C. Cohen. Higher Order Numerical Methods for Transient Wave Equations. Springer Verlag, 2002.
[9] F. Collino, T. Fouquet, and P. Joly. A conservative space-time mesh refinement method for the 1-D wave equation. I. Construction. Numer. Math., 95(2):197–221, 2003.
[10] F. Collino, T. Fouquet, and P. Joly. A conservative space-time mesh refinement method for the 1-D wave equation. II. Analysis. Numer. Math., 95(2):223–251, 2003.
[11] F. Collino, T. Fouquet, and P. Joly. Conservative space-time mesh refinement methods for the FDTD solution of Maxwell’s equations. J. Comput. Phys., 211(1):9–35, 2006.
[12] S. Descombes, S. Lantéri, and L. Moya. Locally implicit discontinuous Galerkin method for time domain electromagnetics. J. Sci. Comp., 56:190–218, 2013.
[13] J. Diaz and M. J. Grote. Energy conserving explicit local time stepping for second-order wave equations. SIAM J. Sci. Comput., 31(3):1985–2014, 2009.
[14] J. Diaz and M. J. Grote. Multi-level explicit local time-stepping methods for second-order wave equations. Comput. Methods Appl. Mech. Engrg., 291:240–265, 2015.
[15] NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov/, Release 1.0.13 of 2016-09-16. F. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller and B. V. Saunders, eds.
[16] V. Dolean, H. Fahs, L. Fezoui, and S. Lanteri. Locally implicit discontinuous Galerkin method for time domain electromagnetics. J. Comput. Phys., 229:512–526, 2010.
[17] Ern, Alexandre and Guermond, Jean-Luc. Finite element quasi-interpolation and best approximation. ESAIM: M2AN, 51:1367–1385, 2017.
[18] J. C. Gilbert and P. Joly. Higher order time stepping for second order hyperbolic problems and optimal CFL conditions. In R. Glowinski and P. Neittaanmäki, editors, Partial Differential Equations, volume 16 of Comput. Methods Appl. Sci., pages 67–93. Springer-Verlag, Dordrecht, The Netherlands, 2008.
[19] M. J. Grote, M. Mehlin, and S. A. Sauter. Convergence analysis of energy conserving explicit local time-stepping methods for the wave equation. SIAM J. Numer. Anal., 56(2):994–1021, 2018.
[20] M. J. Grote, S. Michel, and S. A. Sauter. Stabilized leapfrog based local time-stepping method for the wave equation. Math. Comp., 90(332):2603–2643, 2021.
[21] M. J. Grote and T. Mitkova. Explicit local time-stepping methods for Maxwell’s equations. J. Comput. Appl. Math., 234(12):3283–3302, 2010.
[22] M. J. Grote, A. Schneebeli, and D. Schötzau. Discontinuous Galerkin finite element method for the wave equation. SIAM J. Numer. Anal., 44(6):2408–2431, 2006.
[23] M. Hochbruck and A. Sturm. Error analysis of a second-order locally implicit method for linear Maxwell’s equations. SIAM Journal on Numerical Analysis, 54(5):3167–3191, 2016.
[24] M. Hochbruck and A. Sturm. Upwind discontinuous Galerkin space discretization and locally implicit time integration for linear Maxwell’s equations. Math. Comp., 88(317):1121–1153, 2019.
[25] W. Hundsdorfer and J. Verwer. Numerical solution of time-dependent advection-diffusion-reaction equations, volume 33 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 2003.
[26] P. Joly and J. Rodríguez. An error analysis of conservative space-time mesh refinement methods for the one-dimensional wave equation. SIAM J. Numer. Anal., 43(2):825–859, 2005.
[27] J. Lions and E. Magenes. Non-Homogeneous Boundary Value Problems and Applications. Springer-Verlag, Berlin, 1972.
[28] S. Minisini, E. Zhebel, A. Kononov, and W. A. Mulder. Local time stepping with the discontinuous Galerkin method for wave propagation in 3D heterogeneous media. Geophysics, 78:T67–T77, 2013.
[29] E. Montseny, S. Pernet, X. Ferriéres, and G. Cohen. Dissipative terms and local time-stepping improvements in a spatial high order Discontinuous Galerkin scheme for the time-domain Maxwell’s equations. J. Comput. Phys., 227:6795–6820, 2008.
[30] S. Piperno. Symplectic local time-stepping in non-dissipative DGTD methods applied to wave propagation problems. M2AN Math. Model. Numer. Anal., 40(5):815–841, 2006.
[31] M. Rietmann, M. J. Grote, D. Peter, and O. Schenk. Newmark local time stepping on high-performance computing architectures. J. Comput. Phys., 334:308–326, 2017.
[32] T. J. Rivlin. The Chebyshev Polynomials. Wiley, New York, 1974.
[33] J. G. Verwer. Component splitting for semi-discrete Maxwell equations. BIT, 51:427–445, 2011.

	$\displaystyle z_{S,1}^{(n)}$	$\displaystyle=u_{S}^{\left(n\right)}+\frac{1}{2}\left(\frac{\Delta t}{p}\right% )^{2}\left(\frac{2p^{2}}{\omega_{p,\nu}}\,\beta_{p,\nu}^{(0,0)}\left(w_{S}^{(n% )}-A^{S}\Pi_{\operatorname{f}}^{S}u_{S}^{(n)}\right)+\gamma_{p,\nu}^{(0)}\Pi_% {\operatorname{f}}^{S}f_{S,0}^{(n)}\right)$
		$\displaystyle=\mathfrak{L}_{S,1}^{(n)}u_{S}^{(n)}+\frac{1}{2}\left(\frac{% \Delta t}{p}\right)^{2}\left(\frac{2p^{2}}{\omega_{p,\nu}}\,\beta_{p,\nu}^{(0,% 0)}\Pi_{\operatorname{c}}^{S}f_{S,0}^{(n)}+\gamma_{p,\nu}^{(0)}\Pi_{% \operatorname{f}}^{S}f_{S,0}^{(n)}\right)$
		$\displaystyle\overset{\text{(\ref{eq:recursion_Q_pnurk_1})}}{=}\mathfrak{L}_{S% ,1}^{(n)}u_{S}^{(n)}+\frac{\Delta t^{2}}{2}\left(\frac{2}{\omega_{p,\nu}}\beta% _{p,\nu}^{(0,0)}\Pi_{\operatorname{c}}^{S}f_{S,0}^{(n)}+Q_{p,\nu,0,1}^{\Delta t% }\left(A^{S}\Pi_{\operatorname{f}}^{S}\right)\Pi_{\operatorname*{f}}^{S}f_{S,% 0}\right).$

	$\displaystyle\rho$	$\displaystyle:=\left.\left(2\beta_{p,\nu}^{(k,0)}\left(\delta_{p,\nu}-\frac{% \Delta t^{2}}{\omega_{p,\nu}}x\right)P_{p,\nu,k}^{\Delta t}\left(x\right)-% \beta_{p,\nu}^{(k,-1)}P_{p,\nu,k-1}^{\Delta t}\left(x\right)+4\frac{\,\beta_{p% ,\nu}^{(k,0)}}{\omega_{p,\nu}}\right)\right\|_{x=A^{S}\Pi_{\operatorname{f}}^{% S}}\Pi_{\operatorname{c}}^{S}f_{S,0}^{(n)}$
		$\displaystyle+\sum_{r=-(k-1)}^{k-1}\left.\left(2\beta_{p,\nu}^{(k,0)}\left(% \delta_{p,\nu}-\frac{\Delta t^{2}}{\omega_{p,\nu}}x\right)Q_{p,\nu,\|r\|,k}^{% \Delta t}\left(x\right)-\beta_{p,\nu}^{(k,-1)}Q_{p,\nu,\|r\|,k-1}^{\Delta t}% \left(x\right)\right)\right\|_{x=A^{S}\Pi_{\operatorname{f}}^{S}}\Pi_{% \operatorname{f}}^{S}f_{S,r}^{\left(n\right)}$
		$\displaystyle+Q_{p,\nu,k,k+1}^{\Delta t}\left(A^{S}\Pi_{\operatorname{f}}^{S}% \right)\Pi_{\operatorname{f}}^{S}\left(f_{S,k}^{(n)}+f_{S,-k}^{(n)}\right)$

	$\displaystyle\left(A^{S,p,\nu}\right)^{-1}Q_{p,\nu,\|r\|}^{\Delta t}\left(A^{S}% \Pi_{\operatorname{f}}^{S}\right)\Pi_{\operatorname{f}}^{S}$	$\displaystyle=\left(A^{S}\right)^{-1}Q_{p,\nu,\|r\|}^{\Delta t}\left(A^{S}\Pi_{% \operatorname{f}}^{S}\right)\Pi_{\operatorname{f}}^{S}$		(3.38)
		$\displaystyle+\left(\left(A_{\operatorname{f}\operatorname{f}}^{S,p,\nu}\right% )^{-1}-\left(A_{\operatorname{f}\operatorname{f}}^{S}\right)^{-1}\right)Q_{p,% \nu,\|r\|}^{\Delta t}\left(A_{\operatorname{f}\operatorname{f}}^{S}\right)\Pi_{% \operatorname*{f}}^{S},$

	$\displaystyle\left\\|Q_{p,\nu,\|r\|}^{\Delta t}\left(0\right)\left(A^{S}\right)^{% -1}\Pi_{\operatorname*{f}}^{S}\right\\|_{S\leftarrow S}$	$\displaystyle=\frac{p-r}{p^{2}}\left\\|\left(A^{S}\right)^{-1}\Pi_{% \operatorname*{f}}^{S}\right\\|_{S\leftarrow S}$		(3.43)
		$\displaystyle\overset{\text{\cite[cite]{[\@@bibref{}{grote2021stabilized}{}{},% (3.7)]}}}{\leq}\frac{C_{\operatorname{eq}}^{2}}{c_{\operatorname{coer}}}% \frac{p-r}{p^{2}}\left\\|\Pi_{\operatorname{f}}^{S}\right\\|_{S\leftarrow S}% \leq\frac{C_{\operatorname*{eq}}^{2}}{c_{\operatorname{coer}}}\frac{p-r}{p^{2}}.$

	$\displaystyle\left\\|\left(A^{S,p,\nu}\right)^{-1}\sum_{r=-(p-1)}^{p-1}Q_{p,\nu% ,\|r\|}^{\Delta t}\left(A^{S}\Pi_{\operatorname{f}}^{S}\right)\Pi_{% \operatorname{f}}^{S}\partial_{t}^{\ell}f_{S}\left(t+\frac{r\Delta t}{p}% \right)\right\\|_{\mathcal{T}}$		(3.46)
	$\displaystyle\qquad\leq\left(\frac{C_{\operatorname{eq}}^{2}}{c_{% \operatorname{coer}}}+\frac{\operatorname{e}\nolimits^{\nu}\Delta t^{2}}{6}% \right)\sum_{r=-(p-1)}^{p-1}\frac{p-\left\|r\right\|}{p^{2}}\max_{t\in[t_{m-1},t% _{m+2}]}\left\\|\partial_{t}^{\ell}f_{S}\left(t\right)\right\\|_{\mathcal{T}}$
	$\displaystyle\qquad=\left(\frac{C_{\operatorname{eq}}^{2}}{c_{\operatorname{% coer}}}+\frac{\operatorname{e}\nolimits^{\nu}\Delta t^{2}}{6}\right)\max_{t% \in[t_{m-1},t_{m+2}]}\left\\|\partial_{t}^{\ell}f_{S}\left(t\right)\right\\|_{% \mathcal{T}}.$

Explicit Local Time-Stepping for the Inhomogeneous Wave Equation with Optimal Convergence

Abstract

1 Introduction

2 Stabilized LF-LTS Method

2.1 The Wave Equation

Example 1

2.2 Galerkin Finite Element Spatial Discretization

2.3 Stabilized LF-LTS(ν𝜈\nuitalic_ν) Algorithm

Algorithm 1 (LF-LTS(ν)𝜈(\nu)( italic_ν ) Galerkin FE Algorithm)

Remark 1

3 Convergence Theory in the L2superscript𝐿2L^{2}italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT-norm

3.1 Two-Step Formulation via Chebyshev Polynomials

Remark 2

Lemma 1

Lemma 2

Lemma 3

Corollary 4

3.2 Main Convergence Results

Lemma 5

Remark 3

Theorem 6

Theorem 7

4 Towards Optimal H1superscript𝐻1H^{1}italic_H start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT Convergence

Definition 1 (connecting sequence)

Definition 2 (graph distance)

Definition 3 (discrete distance function)

Definition 4 (fine/coarse mesh mapping)

Remark 4

5 Numerical Experiments

5.1 Optimal L2superscript𝐿2L^{2}italic_L start_POSTSUPERSCRIPT 2 end_POSTSUPERSCRIPT Convergence

5.2 Optimal H1superscript𝐻1H^{1}italic_H start_POSTSUPERSCRIPT 1 end_POSTSUPERSCRIPT-Convergence and weighted transition

5.2.1 Space-time local forcing

5.2.2 Spatially constant solution

Appendix A Appendix

Lemma 8

References

2.3 Stabilized LF-LTS( $\nu$ ) Algorithm

Algorithm 1 (LF-LTS $(\nu)$ Galerkin FE Algorithm)

3 Convergence Theory in the $L^{2}$ -norm

4 Towards Optimal $H^{1}$ Convergence

5.1 Optimal $L^{2}$ Convergence

5.2 Optimal $H^{1}$ -Convergence and weighted transition