# Efficient and Parallel Solution of High-order Continuous Time Galerkin for
Dissipative and Wave Propagation Problems111This work is supported in part by
China National Key Technologies R&D Program under the grant 2019YFA0709600 and
China NSF under the grant 118311061, 12288201, 12201621.

Zhiming Chen222LSEC, Institute of Computational Mathematics, Academy of
Mathematics and System Sciences and School of Mathematical Science, University
of Chinese Academy of Sciences, Chinese Academy of Sciences, Beijing 100190,
China. E-mail<EMAIL_ADDRESS>Yong Liu333LSEC, Institute of
Computational Mathematics, Academy of Mathematics and Systems Science, Chinese
Academy of Sciences, Beijing 100190, P.R. China. E-mail<EMAIL_ADDRESS>
Abstract. We propose efficient and parallel algorithms for the implementation
of the high-order continuous time Galerkin method for dissipative and wave
propagation problems. By using Legendre polynomials as shape functions, we
obtain a special structure of the stiffness matrix which allows us to extend
the diagonal Padé approximation to solve ordinary differential equations with
source terms. The unconditional stability, $hp$ error estimates, and $hp$
superconvergence at the nodes of the continuous time Galerkin method are
proved. Numerical examples confirm our theoretical results.

Key words. Implicit time discretization; Padé approximation; Parallel
implementation

AMS classification. 65M60

##  1 Introduction

In this paper, we study the following system of ordinary differential
equations (ODEs)

$\displaystyle{\mathbf{Y}}^{\prime}(t)=\mathbb{D}\mathbf{Y}(t)+\mathbf{R}(t)\text{
in }(0,T),\ \ \ \ \mathbf{Y}(0)=\mathbf{Y}_{0},$ (1.1)

which is obtained from the method-of-lines approach for linear partial
differential equations (PDEs) after space discretization. Here $T>0$ is the
length of the time interval, $\mathbf{Y},\mathbf{R}\in\mathbb{R}^{M}$, and
$\mathbb{D}$ is an $M\times M$ real constant matrix, where $M$ is the number
of degrees of freedom of the spatial discretization. Without loss of
generality, we assume

$\displaystyle\mathbb{D}+\mathbb{D}^{T}\leq 0,$ (1.2)

that is, $\mathbb{D}+\mathbb{D}^{T}$ is a semi-negative definite matrix. This
condition is satisfied by a large class of linear PDEs including the
dissipative problems such as the parabolic equations and the wave propagation
problems such as the wave equation and Maxwell equations.

Let $0=t_{0}<t_{1}<\cdots<t_{N}=T$ be a partition of $(0,T)$. If the source
$\mathbf{R}=\bf{0}$ in (1.1), the exact solution in each time interval
$(t_{n},t_{n+1})$ is $\mathbf{Y}(t)=e^{\mathbb{D}(t-t_{n})}\mathbf{Y}(t_{n})$
for which Padé approximation to the exponential function can be used to
construct and analyze numerical schemes to solve (1.1). In [12], by using the
partial fraction formula for the Padé approximation, the $[r/r]$, $r\geq 1$,
Padé approximation leads to the following method

$\mathbf{Y}(t_{n+1})\approx\frac{P_{r}(\tau_{n}\mathbb{D})}{P_{r}(-\tau_{n}\mathbb{D})}\mathbf{Y}(t_{n})=\left[(-1)^{r}\mathbb{I}+\sum^{r}_{j=1}-\frac{P_{r}(-\zeta_{j})}{P^{\prime}_{r}(\zeta_{j})}(\zeta_{j}\mathbb{I}+\tau_{n}\mathbb{D})^{-1}\right]\mathbf{Y}(t_{n}),$
(1.3)

where $\tau_{n}=t_{n+1}-t_{n}$, $\mathbb{I}\in\mathbb{R}^{M\times M}$ is the
identity matrix, $P_{r}(z)$ is the numerator of the $[r/r]$ Padé approximation
to the exponential function $e^{z}$, and $\\{\zeta_{1},\cdots,\zeta_{r}\\}$
are zeros of $P_{r}(z)$ which are known to be simple and lie in the left-half
plane. (1.3) indicates that one can compute the approximation of the solution
$\mathbf{Y}(t_{n+1})$ in each time step by solving $k$ complex matrix problems
and $r-2k$ real matrix problems of the form
$\zeta\mathbb{I}+\tau_{n}\mathbb{D}$ in parallel, where $k$, $0\leq k\leq
r/2$, is the number of complex zeros of $P_{r}(z)$ (see Remark 3.1 below). The
purpose of this paper is to construct algorithms sharing this very desirable
property for solving (1.1) with general nonzero sources $\mathbf{R}(t)$.

There exists a large literature on implicit single-step time-stepping methods
for solving (1.1) (see, e.g.,[16] and the references therein). The following
continuous time Galerkin method proposed in [17] is probably the simplest

$\mathbf{Y}^{\prime}_{r}=\mathbb{D}\mathcal{P}_{r-1}\mathbf{Y}_{r}+\mathcal{P}_{r-1}\mathbf{R}\
\ \ \mbox{in }(t_{n},t_{n+1}),\ \ \ 0\leq n\leq N-1,$ (1.4)

where $\mathbf{Y}_{r}$ is a piecewise polynomial of degree $r\geq 1$ in each
interval $(t_{n},t_{n+1})$, continuous at the nodes $t=t_{n}$, and
$\mathcal{P}_{r-1}$ is the local $L^{2}$ projection to the space of
polynomials of degree $(r-1)$ in each interval. It is shown in [17] that (1.4)
is equivalent to the $r$-stage Gauss collocation method at the nodes when
$\mathbf{R}=\bf{0}$ and has the highest classical order $2r$ among all
$r$-stage Runge-Kutta methods [16, Table 5.12]. The continuous time Galerkin
method, together with finite element discretization in space, is used in [2]
for the heat equation and in [11], [15] for the wave equation. We refer to [1]
for a unified framework and the comparison of the most popular implicit
single-step time-stepping methods including also the discontinuous time
Galerkin method and various Runge-Kutta methods.

The difficulty in using the high-order continuous time Galerkin method or any
implicit time Runge-Kutta methods is that a straightforward implementation
requires to solve a system of linear equations of the size $rM\times rM$,
which is not feasible in most time for PDE problems. In a recent work [23],
efficient iterative algorithms are developed based on optimal preconditioning
of the stage matrix for finding the stage vectors of the implicit Runge-Kutta
methods for solving (1.1). For an $r$-stage implicit Runge-Kutta method, the
stage matrix is an $r\times r$ block matrix with each block being a $M\times
M$ matrix. One can find further references in [23] for developing efficient
algorithms implementing the high order implicit time discretization methods in
the literature. We also refer to [22], [19] for the implementation of the
discontinuous time Galerkin method based on the block diagnalization of the
stiffness matrix.

In this paper we propose an efficient realization of the method (1.4) which
uses Legendre polynomials as shape functions to obtain a new stiffness matrix
which is different from the stage matrix in [23] applying to the Gauss
collocation method. By exploiting the special structure of the stiffness
matrix, we construct an algorithm which computes the solution
$\mathbf{Y}_{r}(t_{n})$ at each node by solving $k$ complex matrix problems
and $r-2k$ real matrix problems in parallel, where $k$, $0\leq k\leq r/2$, is
the number of complex zeros of the $[r/r]$ Padé numerator $P_{r}(z)$.
Moreover, a parallel-in-time algorithm is proposed to compute the other
coefficients of the solution $\mathbf{Y}_{r}$ in each time interval
$(t_{n},t_{n+1})$ which solves in parallel $kr$ complex matrix problems and
$(r-2k)r$ real matrix problems. For the dissipative system, in which
$\mathbb{D}+\mathbb{D}^{T}$ is negative definite, in the parallel-in-time
algorithm, only $k$ complex matrix problems and $(r-2k)$ real matrix problems
need to be solved. We remark that our parallel-in-time algorithm is different
from the other parallel-in-time algorithms based on domain decomposition or
space-time multigrid techniques in the literature (see, e.g., [13]).

As a by-product of our analysis, we obtain the following formula (Theorem 3.3)
to compute the nodal values $\mathbf{Y}_{r}(t_{n+1})$, $0\leq n\leq N-1$, of
the solution of (1.4)

$\mathbf{Y}_{r}(t_{n+1})=\frac{P_{r}(\tau_{n}\mathbb{D})}{P_{r}(-\tau_{n}\mathbb{D})}\mathbf{Y}_{r}(t_{n})+\sum^{r}_{k=1}(-1)^{k+1}\frac{\phi_{k1}(\tau_{n}\mathbb{D})}{P_{r}(-\tau_{n}\mathbb{D})}\,\mathbf{b}_{k-1}+\tau_{n}\mathbf{R}_{0},$
(1.5)

where for $k=1,\cdots,r$, $\phi_{k1}(\lambda)$ is a polynomial of degree $r$
satisfying some recurrence relations, and $\mathbf{b}_{k}$, $\mathbf{R}_{0}$
are vectors depending on the source $\mathbf{R}$. (1.5) can be viewed as a
generalization of the $[r/r]$ Padé approximation (1.3) for solving the ODE
system without sources.

The layout of the paper is as follows. In section 2 we introduce the
continuous time Galerkin method for (1.1) and prove the strong stability and
derive a $hp$ error estimate. In section 3 we propose our parallel algorithms
to implement the continuous time Galerkin method. In section 4 we consider an
alternative implementation for the dissipative system. In section 5 we prove
the optimal stability and error estimates in terms of $r$ when $\mathbb{D}$ is
a symmetric or skew-symmetric matrix. In section 6 we consider the application
of the algorithms in this paper to solve the linear convection-diffusion
equation by using the local discontinuous Galerkin method and the wave
equation with discontinuous coefficients by using the unfitted finite element
spatial discretization.

##  2 Implicit time discretization

In this section, we introduce the continuous time Galerkin method for solving
(1.1). Let $0=t_{0}<t_{1}<\ldots<t_{N}=T$ be a partition of the time interval
$(0,T)$ with time steps $\tau_{n}=t_{n+1}-t_{n}$, $0\leq n\leq N-1$. We set
$I_{n}=(t_{n},t_{n+1})$ and $\tau=\displaystyle\max_{0\leq n\leq
N-1}\\{\tau_{n}\\}$. For any integer $m\geq 1$, we define the finite element
space

$\displaystyle\mathbf{V}_{\tau}^{m}:=\\{\mathbf{v}\in[C(0,T)]^{M}:\mathbf{v}|_{I_{n}}\in[P^{m}]^{M},0\leq
n\leq N-1\\},$

where $P^{m}$ is the set of polynomials whose degree is at most $m$. Define
the local projection $\mathcal{P}_{m}$, $m\geq 0$, such that in each time
interval $I_{n}$, $\mathcal{P}_{m}:[L^{2}(I_{n})]^{M}\rightarrow[P^{m}]^{M}$
satisfies

$\displaystyle\int_{I_{n}}(\mathcal{P}_{m}\mathbf{v},\mathbf{w})\,dt=\int_{I_{n}}(\mathbf{v},\mathbf{w})\,dt\quad\forall\mathbf{w}\in[P^{m}]^{M},$

where we denote $(\cdot,\cdot)$ the inner product of $\mathbb{R}^{M}$. It is
well-known (see, e.g., Schwab [21]) that for any $s\geq 0$, $m\geq 0$,

$\displaystyle\|\mathbf{v}-\mathcal{P}_{m}\mathbf{v}\|_{L^{2}(I_{n})}\leq
C\frac{\,\tau_{n}^{\min(m+1,s)}}{(m+1)^{s}}\|\mathbf{v}\|_{H^{s}(I_{n})}\ \
\forall\mathbf{v}\in[H^{s}(I_{n})]^{M},$ (2.1)

where the constant $C$ is independent of $m,\tau_{n}$ but may depend on $s$.
In this paper, for any integer $d\geq 1$ and Banach space $X$, we denote
$\|\cdot\|_{X}$ both the norm of $X$ and $[X]^{d}$.

For any integer $r\geq 1$, the continuous time Galerkin method for solving
(1.1) is to find the function $\mathbf{Y}_{r}\in\mathbf{V}_{\tau}^{r}$ such
that $\mathbf{Y}_{r}(0)=\mathbf{Y}_{0}$ and

$\displaystyle\mathbf{Y}_{r}^{\prime}=\mathbb{D}\mathcal{P}_{r-1}\mathbf{Y}_{r}+\mathcal{P}_{r-1}\mathbf{R}\
\ \mbox{in }I_{n},\ \ 0\leq n\leq N-1.$ (2.2)

The following stability lemma extends an idea in Griesmaier and Monk [15]
where the continuous time Galerkin discretization in time and hybridizable
discontinuous Galerkin method in space for the wave equation are considered.

######  Lemma 2.1.

The problem (2.2) has a unique solution
$\mathbf{Y}_{r}\in\mathbf{V}^{r}_{\tau}$ which satisfies

$\max_{1\leq n\leq
N}\|\mathbf{Y}_{r}(t_{n})\|_{\mathbb{R}^{M}}\leq\|\mathbf{Y}_{0}\|_{\mathbb{R}^{M}}+CT^{1/2}\|\mathbf{R}\|_{L^{2}(0,T)},$
(2.3) $\max_{0\leq t\leq T}\|\mathbf{Y}_{r}\|_{\mathbb{R}^{M}}\leq
Cr^{2}(\|\mathbf{Y}_{0}\|_{\mathbb{R}^{M}}+T^{1/2}\|\mathbf{R}\|_{L^{2}(0,T)}).$
(2.4)

where the constant $C$ is independent of $r,\tau,\mathbb{D}$ and $\mathbf{R}$.

###### Proof.

At each time step, (2.2) is equivalent to a linear system of equations whose
existence and uniqueness of the solution follow from the stability estimate
(2.4). To prove the stability estimates (2.3)-(2.4), we denote by
$\\{L_{j}\\}^{\infty}_{j=0}$ the Legendre polynomials on $(-1,1)$ and define
$\widetilde{L}_{j}=L_{j}\circ\psi^{-1}$, where $\psi:(-1,1)\to(t_{n},t_{n+1})$
is the mapping $\psi(\xi)=\frac{t_{n}+t_{n+1}}{2}+\frac{t_{n+1}-t_{n}}{2}\xi$
for $\xi\in(-1,1)$. Then $\\{\widetilde{L}_{j}\\}^{\infty}_{j=0}$ are
orthogonal in $L^{2}(I_{n})$,
$\widetilde{L}_{r}(t_{n})=(-1)^{r},\widetilde{L}_{r}(t_{n+1})=1$, and

$\int_{I_{n}}|\widetilde{L}_{r}|^{2}dt=\frac{\tau_{n}}{2r+1},\ \
\int_{I_{n}}|\widetilde{L}_{r}^{\prime}|^{2}dt=\frac{2r(r+1)}{\tau_{n}}.$
(2.5)

For $n=0,\cdots,N-1$, let $\mathbf{Y}_{r}^{n}=\mathbf{Y}_{r}(t_{n})$ and
$\hat{\mathbf{Y}}_{r}\in[P^{r}]^{M}$ satisfy

$\hat{\mathbf{Y}}^{\prime}_{r}=\mathbb{D}\mathcal{P}_{r-1}\hat{\mathbf{Y}}_{r}\
\ \mbox{in }I_{n},\ \ \hat{\mathbf{Y}}_{r}(t_{n})=\mathbf{Y}^{n}_{r}.$ (2.6)

By multiplying (2.6) by $\hat{\mathbf{Y}}_{r}$ and integrating over $I_{n}$,
we obtain easily by (1.2) that

$\displaystyle\frac{1}{2}\|\hat{\mathbf{Y}}_{r}(t_{n+1})\|_{\mathbb{R}^{M}}^{2}-\frac{1}{2}\|\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}}^{2}=\int_{I_{n}}(\mathbb{D}\mathcal{P}_{r-1}\hat{\mathbf{Y}}_{r},\mathcal{P}_{r-1}\hat{\mathbf{Y}}_{r})dt\leq
0.$

This implies

$\|\hat{\mathbf{Y}}_{r}(t_{n+1})\|_{\mathbb{R}^{M}}\leq\|\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}}.$
(2.7)

Since $\hat{\mathbf{Y}}_{r}\in[P^{r}]^{M}$ in $I_{n}$, we have the following
decomposition introduced in [15]

$\hat{\mathbf{Y}}_{r}=(-1)^{r}\mathbf{Y}^{n}_{r}\widetilde{L}_{r}+(t-t_{n})\widetilde{\mathbf{Y}}_{r},\
\ \widetilde{\mathbf{Y}}_{r}\in[P^{r-1}]^{M}.$ (2.8)

Notice that
$\mathcal{P}_{r-1}\hat{\mathbf{Y}}_{r}=\mathcal{P}_{r-1}[(t-t_{n})\widetilde{\mathbf{Y}}_{r}]$,
substituting this decomposition to (2.6), we have

$\displaystyle(-1)^{r}\mathbf{Y}^{n}_{r}\widetilde{L}_{r}^{\prime}+\widetilde{\mathbf{Y}}_{r}+(t-t_{n})\widetilde{\mathbf{Y}}_{r}^{\prime}=\mathbb{D}\mathcal{P}_{r-1}[(t-t_{n})\widetilde{\mathbf{Y}}_{r}]\
\ \mbox{in }I_{n}.$

Multiply the equation by $\widetilde{\mathbf{Y}}_{r}\in[P^{r-1}]^{M}$ and
integrate over $I_{n}$, we have by (2.5) that

$\displaystyle\frac{1}{2}\|\widetilde{\mathbf{Y}}_{r}\|_{L^{2}(I_{n})}^{2}+\frac{1}{2}\tau_{n}\|\widetilde{\mathbf{Y}}_{r}(t_{n+1})\|_{\mathbb{R}^{M}}^{2}\leq
C\tau_{n}^{-1/2}r\|\mathbf{Y}_{r}^{n}\|_{\mathbb{R}^{M}}\|\widetilde{\mathbf{Y}}_{r}\|_{L^{2}(I_{n})},$

where we have used the fact that by (1.2)

$\int_{I_{n}}(\mathbb{D}\mathcal{P}_{r-1}[(t-t_{n})\widetilde{\mathbf{Y}}_{r}],\widetilde{\mathbf{Y}}_{r})dt=\int_{I_{n}}(t-t_{n})(\mathbb{D}\widetilde{\mathbf{Y}}_{r},\widetilde{\mathbf{Y}}_{r})dt\leq
0.$ (2.9)

This yields $\|\widetilde{\mathbf{Y}}_{r}\|_{L^{2}(I_{n})}\leq
C\tau_{n}^{-1/2}r\|\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}}$ and thus by using
(2.5)

$\|\hat{\mathbf{Y}}_{r}\|_{L^{2}(I_{n})}\leq
C\tau_{n}^{1/2}r\|\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}}.$ (2.10)

On the other hand, it follows from (2.2) and (2.6) that

$(\mathbf{Y}_{r}-\hat{\mathbf{Y}}_{r})^{\prime}=\mathbb{D}\mathcal{P}_{r-1}(\mathbf{Y}_{r}-\hat{\mathbf{Y}}_{r})+\mathcal{P}_{r-1}\mathbf{R}\
\ \mbox{in }I_{n},\ \ (\mathbf{Y}_{r}-\hat{\mathbf{Y}}_{r})(t_{n})=0.$ (2.11)

Then $\mathbf{Y}_{r}-\hat{\mathbf{Y}}_{r}=(t-t_{n})\mathbf{W}_{r}$ for some
$\mathbf{W}_{r}\in[P^{r-1}]^{M}$. By substituting this relation into the
equation (2.11) we have

$\displaystyle\mathbf{W}_{r}+(t-t_{n})\mathbf{W}_{r}^{\prime}=\mathbb{D}\mathcal{P}_{r-1}[(t-t_{n})\mathbf{W}_{r}]+\mathcal{P}_{r-1}\mathbf{R}\
\ \mbox{in }I_{n}.$

By multiplying the equation by $\mathbf{W}_{r}$ and integrating over $I_{n}$,
we obtain by a similar bound as in (2.9) that

$\displaystyle\frac{1}{2}\|\mathbf{W}_{r}\|_{L^{2}(I_{n})}^{2}+\frac{1}{2}\tau_{n}\|\mathbf{W}_{r}(t_{n+1})\|_{\mathbb{R}^{M}}^{2}\leq\|\mathbf{R}\|_{L^{2}(I_{n})}\|\mathbf{W}_{r}\|_{L^{2}(I_{n})}.$

This yields $\|\mathbf{W}_{r}\|_{L^{2}(I_{n})}\leq
2\|\mathbf{R}\|_{L^{2}(I_{n})}$ and thus

$\|\mathbf{Y}_{r}-\hat{\mathbf{Y}}_{r}\|_{L^{2}(I_{n})}\leq
2\tau_{n}\|\mathbf{R}\|_{L^{2}(I_{n})}.$ (2.12)

Now by multiplying (2.11) by $\mathbf{Y}_{r}-\hat{\mathbf{Y}}_{r}$ and
integrating over $I_{n}$ we obtain by (1.2) and (2.12) that

$\displaystyle\frac{1}{2}\|(\mathbf{Y}_{r}-\hat{\mathbf{Y}}_{r})(t_{n+1})\|_{\mathbb{R}^{M}}^{2}\leq\int_{I_{n}}(\mathcal{P}_{r-1}\mathbf{R},\mathbf{Y}_{r}-\hat{\mathbf{Y}}_{r})dt\leq
2\tau_{n}\|\mathbf{R}\|_{L^{2}(I_{n})}^{2},$

which implies by the triangle inequality and (2.7) that

$\displaystyle\|\mathbf{Y}_{r}(t_{n+1})\|_{\mathbb{R}^{M}}\leq\|\mathbf{Y}_{r}^{n}\|_{\mathbb{R}^{M}}+2\tau_{n}^{1/2}\|\mathbf{R}\|_{L^{2}(I_{n})}.$

This yields (2.3). Next by using the triangle inequality, (2.10), and (2.12),
we have

$\|\mathbf{Y}_{r}\|_{L^{2}(I_{n})}\leq
C\tau^{1/2}_{n}r\|\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}}+2\tau_{n}\|\mathbf{R}\|_{L^{2}(I_{n})},$
(2.13)

which implies by the $hp$ inverse estimate that

$\displaystyle\max_{t_{n}\leq t\leq
t_{n+1}}\|\mathbf{Y}_{r}\|_{\mathbb{R}^{M}}\leq
C\tau_{n}^{-1/2}r\|\mathbf{Y}_{r}\|_{L^{2}(I_{n})}\leq
Cr^{2}\|\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}}+C\tau_{n}^{1/2}r\|\mathbf{R}\|_{L^{2}(I_{n})}.$

This shows (2.4) and completes the proof of the lemma. ∎

To derive an $hp$ a priori error estimate for the continuous time Galerkin
method (2.2), we first recall an interpolation operator in the literature
(see, e.g., [21, Theorem 3.17]).

######  Lemma 2.2.

There exists an interpolation operator
$\Pi_{r}:[H^{1}(0,T)]^{M}\to\mathbf{V}^{r}_{\tau}$ such that for any
$\mathbf{v}\in[W^{1+s,\infty}(0,T)]^{M}$, $s\geq 1$, and $n=0,1,\cdots,N-1$,

$\displaystyle(\Pi_{r}\mathbf{v})(t_{n})=\mathbf{v}(t_{n}),\ \
(\Pi_{r}\mathbf{v})(t_{n+1})=\mathbf{v}(t_{n+1}),\ \
(\Pi_{r}\mathbf{v})^{\prime}=\mathcal{P}_{r-1}\mathbf{v}^{\prime}\ \ \mbox{in
}I_{n},$ (2.14)
$\displaystyle\|\mathbf{v}-\Pi_{r}\mathbf{v}\|_{L^{2}(I_{n})}\leq
C\frac{\tau^{\min(r+1,s)}}{r^{s}}\|\mathbf{v}\|_{H^{s}(I_{n})},$ (2.15)
$\displaystyle\max_{t_{n}\leq t\leq
t_{n+1}}\|\mathbf{v}-\Pi_{r}\mathbf{v}\|_{\mathbb{R}^{M}}\leq
C\frac{\tau^{\min(r,s)+1}}{r^{s}}\|\mathbf{v}^{\prime}\|_{W^{s,\infty}(I_{n})},$
(2.16)

where the constant $C$ is independent of $\tau,r$ but may depend on $s$.

###### Proof.

The interpolation operator is defined as

$\displaystyle\Pi_{r}\mathbf{v}=\mathbf{v}(t_{n})+\int_{t_{n}}^{t}(\mathcal{P}_{r-1}\mathbf{v}^{\prime})dt\
\ \ \ \forall t\in I_{n}.$

(2.14) follows easily from this definition. Next by using (2.1), we have

$\displaystyle\max_{t_{n}\leq t\leq
t_{n+1}}\|\mathbf{v}-\Pi_{r}\mathbf{v}\|_{\mathbb{R}^{M}}\leq\tau_{n}^{1/2}\|\mathbf{v}^{\prime}-\mathcal{P}_{r-1}\mathbf{v}^{\prime}\|_{L^{2}(I_{n})}\leq
C\frac{\tau_{n}^{\min(r,s)+1}}{r^{s}}\|\mathbf{v}^{\prime}\|_{W^{s,\infty}(I_{n})}.$

This shows (2.16).

The estimate (2.15) is proved for $s\geq 2$ in [21]. Here we use the duality
argument to show (2.15) also from $s\geq 1$. Let $\mathbf{w}\in
H^{1}_{0}(I_{n})$ be the solution of the problem

$\displaystyle-\mathbf{w}^{\prime\prime}=\mathbf{v}-\Pi_{r}\mathbf{v}\ \
\mbox{in }I_{n}.$

It is easy to see that $\|\mathbf{w}\|_{H^{2}(I_{n})}\leq
C\|\mathbf{v}-\Pi_{r}\mathbf{v}\|_{L^{2}(I_{n})}$. Since
$(\mathbf{v}-\Pi_{r}\mathbf{v})(t_{n})=\mathbf{0},(\mathbf{v}-\Pi_{r}\mathbf{v})(t_{n+1})=\mathbf{0}$,
we multiply the equation by $\mathbf{v}-\Pi_{r}\mathbf{v}$, integrate over
$I_{n}$, and use (2.14) to obtain

$\displaystyle\|\mathbf{v}-\Pi_{r}\mathbf{v}\|_{L^{2}(I_{n})}^{2}=\int_{I_{n}}(\mathbf{w}^{\prime},\mathbf{v}^{\prime}-(\Pi_{r}\mathbf{v})^{\prime})dt$
$\displaystyle=$
$\displaystyle\int_{I_{n}}(\mathbf{w}^{\prime}-\mathcal{P}_{r-1}\mathbf{w}^{\prime},\mathbf{v}^{\prime}-\mathcal{P}_{r-1}\mathbf{v}^{\prime})dt$
$\displaystyle\leq$ $\displaystyle
C\frac{\tau^{\min(r+1,s)}}{r^{s}}\|\mathbf{w}\|_{H^{2}(I_{n})}\|\mathbf{v}\|_{H^{s}(I_{n})}.$

This completes the proof by using $\|\mathbf{w}\|_{H^{2}(I_{n})}\leq
C\|\mathbf{v}-\Pi_{r}\mathbf{v}\|_{L^{2}(I_{n})}$. ∎

The following theorem on the $hp$ error estimate is the main result of this
section.

######  Theorem 2.1.

Let $s\geq 1$. Assume that $\mathbf{R}\in[H^{s}(0,T)]^{M}$,
$\mathbf{Y}\in[W^{1+s,\infty}(0,T)]^{M}$ and
$\mathbf{Y}_{r}\in\mathbf{V}^{r}_{\tau}$ is the solution of the problem (2.2),
we have

$\displaystyle\max_{1\leq n\leq
N}\|(\mathbf{Y}-\mathbf{Y}_{r})(t_{n})\|_{\mathbb{R}^{M}}\leq
CT^{1/2}\frac{\tau^{\min(r+1,s)}}{r^{s}}\|\mathbb{D}\mathbf{Y}\|_{H^{s}(0,T)},$
$\displaystyle\max_{0\leq t\leq
T}\|\mathbf{Y}-\mathbf{Y}_{r}\|_{\mathbb{R}^{M}}\leq
C(1+T^{1/2})\frac{\tau^{\min(r+1,s)}}{r^{s-2}}(T^{1/2}\|\mathbf{Y}\|_{W^{s+1,\infty}(0,T)}+\|\mathbf{R}\|_{H^{s}(0,T)}),$

where the constant $C$ is independent of $\tau,r,\mathbb{D}$ but may depend on
$s$.

###### Proof.

Let $\Pi_{r}\mathbf{Y}\in\mathbf{V}^{r}_{\tau}$ be the interpolation of
$\mathbf{Y}$ defined in Lemma 2.2. Since
$(\Pi_{r}\mathbf{Y})^{\prime}=\mathcal{P}_{r-1}\mathbf{Y}^{\prime}$ in
$I_{n}$, we have

$\displaystyle(\Pi_{r}\mathbf{Y})^{\prime}=\mathcal{P}_{r-1}(\mathbb{D}\mathbf{Y}+\mathbf{R})=\mathbb{D}\mathcal{P}_{r-1}(\mathbf{Y}-\Pi_{r}\mathbf{Y})+\mathbb{D}\mathcal{P}_{r-1}(\Pi_{r}\mathbf{Y})+\mathcal{P}_{r-1}\mathbf{R}\
\ \mbox{in }I_{n}.$

Thus by (2.2) we have

$\mathbf{Y}_{r}^{\prime}-(\Pi_{r}\mathbf{Y})^{\prime}=\mathbb{D}\mathcal{P}_{r-1}(\mathbf{Y}_{r}-\Pi_{r}\mathbf{Y})-\mathbb{D}\mathcal{P}_{r-1}(\mathbf{Y}-\Pi_{r}\mathbf{Y}).$
(2.17)

As $(\mathbf{Y}_{r}-\Pi_{r}\mathbf{Y})(0)=\mathbf{0}$, we use (2.3) and (2.15)
to obtain

$\displaystyle\max_{1\leq n\leq
N}\|(\mathbf{Y}_{r}-\Pi_{r}\mathbf{Y})(t_{n})\|_{\mathbb{R}^{M}}$
$\displaystyle\leq$ $\displaystyle
CT^{1/2}\|\mathbb{D}\mathcal{P}_{r-1}(\mathbf{Y}-\Pi_{r}\mathbf{Y})\|_{L^{2}(0,T)}$
$\displaystyle\leq$ $\displaystyle
CT^{1/2}\frac{\tau^{\min(r+1,s)}}{r^{s}}\|\mathbb{D}\mathbf{Y}\|_{H^{s}(0,T)}.$

This shows the first estimate as $\mathbf{Y}(t_{n})=\Pi_{r}\mathbf{Y}(t_{n})$.
The second estimate can be proved similarly by using (2.4), (2.16), and
$\mathbf{Y}^{\prime}=\mathbb{D}\mathbf{Y}+\mathbf{R}$. ∎

We remark that the first estimate in Theorem 2.1 is optimal in $\tau$ and $r$
and the second estimate in the theorem is optimal in $\tau$ but suboptimal in
$r$ which is due to the stability estimate (2.4) in Lemma 2.1. In section 5 we
will show that the stability in the $L^{2}$ norm can be improved to remove the
dependence on $r$ in (2.13) when $\mathbb{D}$ is symmetric or skew-symmetric
by using the explicit formulas of $\mathbf{Y}_{r}(t)$ in section 3. We remark
that many spatial discretization matrices of the wave-like equations satisfy
the property that $\mathbb{D}$ is skew-symmetric, such as the energy
conserving mixed finite element methods for solving the Hodge wave equation in
Wu and Bai [25] and the unfitted finite element method of the acoustic wave
equation in Chen et. al. [7].

The classical order of Runge-Kutta methods is the convergence order at the
nodes $t=t_{n}$. For the continuous time Galerkin method, it is proved to be
${2r}$ when $r\geq 2$ in Hulme [17] for nonlinear ODEs and in Aziz and Monk
[2] for parabolic equations. The following theorem shows the $hp$
superconvergence of the continuous time Galerkin method at the nodes by using
the idea of quasi-projection in [2, §4].

######  Theorem 2.2.

Let $s\geq 1$. Assume that $\mathbb{D}^{r}\mathbf{Y}\in[H^{s}(0,T)]^{M}$ and
$\mathbf{Y}_{r}\in\mathbf{V}^{r}_{\tau}$ is the solution of the problem (2.2),
we have

$\displaystyle\max_{1\leq n\leq
N}\|(\mathbf{Y}-\mathbf{Y}_{r})(t_{n})\|_{\mathbb{R}^{M}}\leq
CT^{1/2}\frac{\tau^{\min(2r,s+r-1)}}{r^{s}}\|\mathbb{D}^{r}\mathbf{Y}\|_{H^{s}(0,T)},$

where the constant $C$ is independent of $\tau,r$ but may depend on $s$.

###### Proof.

If $r=1$, the theorem follows from the first estimate of Theorem 2.1. Now we
assume $r\geq 2$. Let $\Pi_{r}\mathbf{Y}\in\mathbf{V}^{r}_{\tau}$ be the
interpolation of $\mathbf{Y}$ defined in Lemma 2.2. Denote
$\bm{\omega}_{0}=\mathbf{Y}-\Pi_{r}\mathbf{Y}$. For $1\leq i\leq r-1$, we
define correction functions $\bm{\omega}_{i}$ such that

$\displaystyle\bm{\omega}_{i}(t_{n})=\mathbf{0},\,\bm{\omega}_{i}^{\prime}=\mathbb{D}\mathcal{P}_{r-1}\bm{\omega}_{i-1}\
\ \text{ in }I_{n},\ \ n=0,1,\cdots,N-1.$ (2.18)

We claim that $\bm{\omega}_{i}(t_{n+1})=\mathbf{0}$ so that
$\bm{\omega}_{i}\in\mathbf{V}^{r}_{\tau}$. In fact, by (2.14), we have
$(\bm{\omega}_{0}^{\prime},\mathbf{v})_{I_{n}}=0$ for any
$\mathbf{v}\in[P^{r-1}]^{M}$, where $(\cdot,\cdot)_{I_{n}}$ is the inner
product of $[L^{2}(I_{n})]^{M}$. Since
$\bm{\omega}_{0}(t_{n})=\bm{\omega}_{0}(t_{n+1})=\bm{0}$, we obtain by
integration by parts that
$(\bm{\omega}_{1}^{\prime},\mathbf{v}^{\prime})_{I_{n}}=(\mathbb{D}\bm{\omega}_{0},\mathbf{v}^{\prime})_{I_{n}}=0$.
Therefore, $(\bm{\omega}_{1}^{\prime},\mathbf{v})_{I_{n}}=0$ for any
$\mathbf{v}\in[P^{r-2}]^{M}$ and consequently,
$\bm{\omega}_{1}(t_{n+1})=\int_{I_{n}}\bm{\omega}_{1}^{\prime}dt=\mathbf{0}$.
By mathematical induction, we know easily by the same argument that
$(\bm{\omega}_{i}^{\prime},\mathbf{v})_{I_{n}}=0$ for any
$\mathbf{v}\in[P^{r-i-1}]^{M}$ and $\bm{\omega}_{i}(t_{n+1})=\mathbf{0}$. This
shows the claim.

Let
$\displaystyle\bm{\omega}=\sum_{i=1}^{r-1}\bm{\omega}_{i}\in\mathbf{V}^{r}_{\tau}$.
By (2.17) and (2.18), we have

$\displaystyle\mathbf{Y}_{r}^{\prime}-(\Pi_{r}\mathbf{Y})^{\prime}+\bm{\omega}^{\prime}=\mathbb{D}\mathcal{P}_{r-1}(\mathbf{Y}_{r}-\Pi_{r}\mathbf{Y}+\bm{\omega})-\mathbb{D}\mathcal{P}_{r-1}(\bm{\omega}_{r-1}).$

As $(\mathbf{Y}_{r}-\Pi_{r}\mathbf{Y}+\bm{\omega})(0)=\mathbf{0}$, we use
(2.3) to obtain

$\displaystyle\|(\mathbf{Y}_{r}-\Pi_{r}\mathbf{Y}+\bm{\omega})(t_{n})\|_{\mathbb{R}^{M}}$
$\displaystyle\leq$ $\displaystyle
CT^{1/2}\|\mathbb{D}\mathcal{P}_{r-1}(\bm{\omega}_{r-1})\|_{L^{2}(0,T)}\leq
CT^{1/2}\|\mathbb{D}\bm{\omega}_{r-1}\|_{L^{2}(0,T)}.$

Now it follows from (2.18) that

$\displaystyle\|\mathbb{D}\bm{\omega}_{i}\|_{L^{2}(I_{n})}\leq\tau\|\mathbb{D}^{2}\bm{\omega}_{i-1}\|_{L^{2}(I_{n})},\quad
1\leq i\leq r-1.$

By using (2.15) we have then

$\displaystyle\|\mathbb{D}\bm{\omega}_{r-1}\|_{L^{2}(I_{n})}\leq\tau^{r-1}\|\mathbb{D}^{r}\bm{\omega}_{0}\|_{L^{2}(I_{n})}\leq
C\frac{\tau^{\min(2r,s+r-1)}}{r^{s}}\|\mathbb{D}^{r}\mathbf{Y}\|_{H^{s}(I_{n})}.$

This completes the proof since by (2.14) and (2.18),
$(\mathbf{Y}_{r}-\Pi_{r}\mathbf{Y}+\bm{\omega})(t_{n})=(\mathbf{Y}_{r}-\mathbf{Y})(t_{n})$,
$1\leq n\leq N$. ∎

The correction function $\bm{\omega}=\sum^{r-1}_{i=1}\bm{\omega}_{i}$ is
introduced in [2] which is related to the idea of quasi-projection in Douglas
Jr. et al [10]. Our new observation is that $\bm{\omega}_{i}=\bf{0}$ at the
nodes for $1\leq i\leq r-1$, which simplifies the proof.

To conclude this section, we recall some facts about Padé approximation to the
exponential function which can be found in Saff and Varga [20] and the
references therein. For any integers $m,n\geq 0$, the $[m/n]$ Padé
approximation to $e^{z}$ is defined as the polynomials $P_{m}(z)\in P^{m}$,
$Q_{n}(z)\in P^{n}$, $Q_{n}(0)=1$, for which

$\displaystyle e^{z}-\frac{P_{m}(z)}{Q_{n}(z)}=O(|z|^{m+n+1})\ \ \ \ \mbox{as
}|z|\to 0.$

It is known that

$P_{m}(z)=\sum^{m}_{j=0}\frac{(m+n-j)!m!z^{j}}{(m+n)!j!(m-j)!},\ \
Q_{n}(z)=\sum^{n}_{j=0}\frac{(m+n-j)!n!(-z)^{j}}{(m+n)!j!(n-j)!}.$ (2.19)

Obviously, $Q_{n}(z)=P_{n}(-z)$. When $m=n$, $P_{m}(z),Q_{m}(z)$ are called
diagonal Padé numerator and denominator of type $[m/m]$ for $e^{z}$. The
following lemma follows easily from (2.19)

######  Lemma 2.3.

The diagonal Padé numerator of type $[m/m]$ for $e^{z}$ satisfies
$P_{0}(z)=1,P_{1}(z)=1+\frac{1}{2}z$,
$P_{2}(z)=1+\frac{1}{2}z+\frac{1}{12}{z^{2}}$, and

$\displaystyle P_{m}(z)=P_{m-1}(z)+\frac{z^{2}}{4(2m-1)(2m-3)}P_{m-2}(z),\ \ \
\ m\geq 2.$

The following lemma is proved in Hairer and Wanner [16, Theorem 4.12], [20,
Theorem 2.4]. It is essential in proving the A-stability of numerical methods
for ODEs based on the Padé approximation of the exponential function.

######  Lemma 2.4.

All zeros of the diagonal Padé numerator of type $[m/m]$, $m\geq 1$, for
$e^{z}$ are simple and lie in the half-plane $\\{z\in\mathbb{C}:{\rm
Re}(z)\leq-2\\}$.

For $m\geq 1$, denote $\zeta_{1},\cdots,\zeta_{m}\in\mathbb{C}$ the zeros of
$P_{m}(z)$, the diagonal Padé numerator of type $[m/m]$ for $e^{z}$, then by
(2.19)

$\displaystyle P_{m}(z)=\frac{m!}{(2m)!}(z-\zeta_{1})\cdots(z-\zeta_{m}),\ \
Q_{m}(z)=(-1)^{m}\frac{m!}{(2m)!}(z+\zeta_{1})\cdots(z+\zeta_{m}).$

Recall that any polynomial $F\in P^{m-1}$ can be expanded as the Lagrange
interpolation function at $m$ distinct zeros of $P_{m}(-z)$. This yields the
following partial fraction formula (see, e.g., Szegö [24, Theorem 3.3.5])

$\frac{F(z)}{P_{m}(-z)}=\sum^{m}_{j=1}-\frac{F(-\zeta_{j})}{P_{m}^{\prime}(\zeta_{j})}\,\frac{1}{z+\zeta_{j}}.$
(2.20)

Since $P_{m}(z)-(-1)^{m}P_{m}(-z)\in P^{m-1}$, we obtain the partial fraction
formula for the $[m/m]$ Padé approximation of $e^{z}$ (see Gallopoulos and
Saad [12])

$R_{m,m}(z)=\frac{P_{m}(z)}{P_{m}(-z)}=(-1)^{m}+\sum^{m}_{j=1}-\frac{P_{m}(-\zeta_{j})}{P_{m}^{\prime}(\zeta_{j})}\,\frac{1}{z+\zeta_{j}}.$
(2.21)

Recall that if $p(z)=\sum^{m}_{i=0}a_{i}z^{i}$ is a polynomial of degree $m$,
then $p(\mathbb{X}):=\sum^{m}_{i=0}a_{i}\mathbb{X}^{i}$ for any matrix
$\mathbb{X}\in\mathbb{R}^{d\times d}$, $d\geq 1$. Obviously, if
$p(z)=p_{1}(z)+p_{2}(z)$ or $q(z)=p_{1}(z)\cdot p_{2}(z)$, where $p_{1},p_{2}$
are polynomials, then $p(\mathbb{X})=p_{1}(\mathbb{X})+p_{2}(\mathbb{X})$,
$q(\mathbb{X})=p_{1}(\mathbb{X})\cdot p_{2}(\mathbb{X})$. It follows now from
(2.20)-(2.21) that for any $F\in P^{m-1}$ and any matrix
$\mathbb{X}\in\mathbb{R}^{d\times d}$ such that $P_{m}(-\mathbb{X})$ is
invertible,

$\displaystyle\frac{F(\mathbb{X})}{P_{m}(-\mathbb{X})}=\sum^{m}_{j=1}-\frac{F(-\zeta_{j})}{P_{m}^{\prime}(\zeta_{j})}(\zeta_{j}\mathbb{I}+\mathbb{X})^{-1},$
(2.22)
$\displaystyle\frac{P_{m}(\mathbb{X})}{P_{m}(-\mathbb{X})}=(-1)^{m}\mathbb{I}+\sum^{m}_{j=1}-\frac{P_{m}(-\zeta_{j})}{P_{m}^{\prime}(\zeta_{j})}\,(\zeta_{j}\mathbb{I}+\mathbb{X})^{-1}.$
(2.23)

The identity (2.23) is the basis of the method (1.3) in the introduction.

##  3 Parallel implementation

In this section, we propose parallel algorithms to implement the problem (2.2)
based on finding analytic formulas of the determinant and all factors of the
stiffness matrix of the continuous time Galerkin method at each time step. To
form the stiffness matrix, we use the Legendre polynomials
$\\{\widetilde{L}_{j}\\}^{r}_{j=0}$ in $I_{n}$ as the basis functions. We
assume

$\displaystyle\mathbf{Y}_{r}(t)=\sum_{j=0}^{r}\mathbf{a}_{j}\widetilde{L}_{j}(t),\quad\mathcal{P}_{r-1}\mathbf{R}=\sum_{j=0}^{r-1}\mathbf{R}_{j}\widetilde{L}_{j}(t),$

where $\mathbf{a}_{j},\mathbf{R}_{j}\in\mathbb{R}^{M}$. Since
$\mathcal{P}_{r-1}\widetilde{L}_{r}=0$ in $I_{n}$, from (2.2), we have

$\displaystyle\sum_{j=0}^{r}\mathbf{a}_{j}\widetilde{L}_{j}^{\prime}(t)=\mathbb{D}\sum_{j=0}^{r-1}\mathbf{a}_{j}\widetilde{L}_{j}(t)+\sum_{j=0}^{r-1}\mathbf{R}_{j}\widetilde{L}_{j}(t).$

For any $k\geq 1$, multiplying the equation by
$(t-t_{n})(t_{n+1}-t)\widetilde{L}^{\prime}_{k}(t)$ and integrating over
$I_{n}$, we obtain

$\displaystyle\mathbf{a}_{k}\frac{\tau_{n}}{2}\frac{k(k+1)}{k+\frac{1}{2}}=$
$\displaystyle\mathbb{D}\sum_{j=0}^{r-1}\int_{t_{n}}^{t_{n+1}}\mathbf{a}_{j}\widetilde{L}_{j}(t)(t-t_{n})(t_{n+1}-t)\widetilde{L}_{k}^{\prime}(t)\,dt$
$\displaystyle+\sum_{j=0}^{r-1}\int_{t_{n}}^{t_{n+1}}\mathbf{R}_{j}\widetilde{L}_{j}(t)(t-t_{n})(t_{n+1}-t)\widetilde{L}_{k}^{\prime}(t)\,dt,$
(3.1)

where we have used the fact that

$\displaystyle\int_{t_{n}}^{t_{n+1}}\widetilde{L}_{j}^{\prime}\widetilde{L}_{k}^{\prime}(t-t_{n})(t_{n+1}-t)\,dt=\frac{\tau_{n}}{2}\int^{1}_{-1}L_{j}^{\prime}L_{k}^{\prime}(1-t^{2})dt=\frac{\tau_{n}}{2}\frac{k(k+1)}{k+\frac{1}{2}}\delta_{j,k}.$

Here $\delta_{j,k}$ is the Kronecker delta function. By the recursion relation
$(2k+1)\widetilde{L}_{k}(t)=\frac{\tau_{n}}{2}(\widetilde{L}^{\prime}_{k+1}(t)-\widetilde{L}^{\prime}_{k-1}(t))$,

$\displaystyle\int_{t_{n}}^{t_{n+1}}\widetilde{L}_{j}(t)(t-t_{n})(t_{n+1}-t)\widetilde{L}^{\prime}_{k}(t)\,dt$
$\displaystyle=$
$\displaystyle\,\frac{1}{2j+1}\int_{t_{n}}^{t_{n+1}}\frac{\tau_{n}}{2}(\widetilde{L}^{\prime}_{j+1}(t)-\widetilde{L}^{\prime}_{j-1}(t))\widetilde{L}^{\prime}_{k}(t)(t-t_{n})(t_{n+1}-t)\,dt$
$\displaystyle=$
$\displaystyle\,\frac{\tau_{n}^{2}}{4}\frac{k(k+1)}{k+\frac{1}{2}}\left(\frac{1}{2k-1}\delta_{j+1,k}-\frac{1}{2k+3}\delta_{j-1,k}\right).$

Substituting the identity into (3.1), we have

$\displaystyle\mathbf{a}_{k}=\frac{\tau_{n}}{2}\mathbb{D}\left(\frac{\mathbf{a}_{k-1}}{2k-1}-\frac{\mathbf{a}_{k+1}}{2k+3}\right)+\frac{\tau_{n}}{2}\left(\frac{\mathbf{R}_{k-1}}{2k-1}-\frac{\mathbf{R}_{k+1}}{2k+3}\right),\quad
1\leq k\leq r-2,$ (3.2)
$\displaystyle\mathbf{a}_{k}=\frac{\tau_{n}}{2}\mathbb{D}\frac{\mathbf{a}_{k-1}}{2k-1}+\frac{\tau_{n}}{2}\frac{\mathbf{R}_{k-1}}{2k-1},\quad
k=r-1,r.$ (3.3)

By the condition $\mathbf{Y}_{r}(t_{n})=\mathbf{Y}_{r}^{n}$, we also have

$\displaystyle\sum_{j=0}^{r}(-1)^{j}\mathbf{a}_{j}=\mathbf{Y}_{r}^{n}.$ (3.4)

(3.2)-(3.4) can be written as a system of linear equations

$\displaystyle\mathbb{A}\mathbf{X}=\mathbf{B},$ (3.5)

where
$\mathbf{X}=(\mathbf{a}_{0}^{T},\mathbf{a}_{1}^{T},\cdots,\mathbf{a}_{r}^{T})^{T}$,
$\mathbf{B}=(\mathbf{b}_{0}^{T},\mathbf{b}_{1}^{T},\cdots,\mathbf{b}_{r}^{T})^{T}$
with

$\displaystyle\mathbf{b}_{k-1}=\left\\{\begin{array}[]{ll}-\frac{\tau_{n}}{2}\left(\frac{\mathbf{R}_{k-1}}{2k-1}-\frac{\mathbf{R}_{k+1}}{2k+3}\right)&\mbox{if
}1\leq k\leq r-2,\\\ -\frac{\tau_{n}}{2}\frac{\mathbf{R}_{k-1}}{2k-1}&\mbox{if
}k=r-1,r,\\\ \mathbf{Y}^{n}_{r}&\mbox{if }k=r+1,\end{array}\right.$

and

$\displaystyle\mathbb{A}=\left(\begin{array}[]{ccccccc}\frac{\tau_{n}}{2}\mathbb{D}&-\mathbb{I}&-\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{5}&\cdots&\cdots&\cdots&0\\\
0&\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{3}&-\mathbb{I}&-\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{7}&\cdots&\cdots&0\\\
\vdots&\vdots&\ddots&\ddots&\ddots&\vdots&\vdots\\\
0&\cdots&\cdots&\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{2r-5}&-\mathbb{I}&-\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{2r-1}&0\\\
0&\cdots&\cdots&\cdots&\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{2r-3}&-\mathbb{I}&0\\\
0&\cdots&\cdots&\cdots&\cdots&\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{2r-1}&-\mathbb{I}\\\
\mathbb{I}&-\mathbb{I}&\mathbb{I}&-\mathbb{I}&\cdots&\cdots&(-1)^{r+2}\mathbb{I}\end{array}\right).$

Here $\mathbb{I}\in\mathbb{R}^{M\times M}$ is the identity matrix. By Lemma
2.1, (3.5) has a unique solution. Since $\mathbb{A}\in\mathbb{R}^{M(r+1)\times
M(r+1)}$, it is expensive to solve (3.5) directly when $M\gg 1$. Here we
propose efficient and parallel algorithms to solve (3.5).

Notice that $\mathbb{A}$ is a $(r+1)\times(r+1)$ block matrix. For
$\lambda\in\mathbb{R}$, we define
$\mathbb{E}_{r+1}(\lambda)\in\mathbb{R}^{(r+1)\times(r+1)}$ by

$\displaystyle\mathbb{E}_{r+1}(\lambda):=\left(\begin{array}[]{ccccccc}a_{1}&-1&b_{1}&\cdots&\cdots&\cdots&0\\\
0&a_{2}&-1&b_{2}&\cdots&\cdots&0\\\
\vdots&\vdots&\ddots&\ddots&\ddots&\vdots&\vdots\\\
0&\cdots&\cdots&a_{r-2}&-1&b_{r-2}&0\\\ 0&\cdots&\cdots&\cdots&a_{r-1}&-1&0\\\
0&\cdots&\cdots&\cdots&\cdots&a_{r}&-1\\\
c_{1}&c_{2}&\cdots&\cdots&\cdots&c_{r}&c_{r+1}\\\ \end{array}\right)$ (3.14)

with

$a_{k}=\frac{\lambda}{2}\frac{1}{2k-1},\ \
b_{k}=-\frac{\lambda}{2}\frac{1}{2k+3},\ \ c_{k}=(-1)^{k+1},\ \
k=1,\cdots,r+1.$ (3.15)

Then $\mathbb{A}=\mathbb{E}_{r+1}(\tau_{n}\mathbb{D})$ by replacing each
element $e_{ij}(\lambda)$ of $\mathbb{E}_{r+1}(\lambda)$ by the $M\times M$
matrix $e_{ij}(\tau_{n}\mathbb{D})\mathbb{I}$, $i,j=1,\cdots,r+1$.

We are going to solve (3.5) by extending the Cramer rule for the block
matrices. We first introduce some notation. For any matrix
$\mathbb{X}=(X_{ij})^{d}_{i,j=1}$, $d\geq 1$, we denote $\mathbb{X}_{i,j}$ the
matrix obtained by removing the $i$-th row and $j$-th column, $i,j=1,\dots,d$.
We also denote $\mathbb{X}_{(i_{1},\cdots,i_{k}),(j_{1},\cdots,j_{l})}$ the
matrix obtained by removing the $i_{1},\cdots,i_{k}$-th rows and the
$j_{1},\cdots,j_{l}$-th columns, where $1\leq i_{1}<\cdots<i_{k}\leq d,1\leq
j_{1}<\cdots<j_{l}\leq d$. The following property about the adjugate matrix is
well known

$(\det\mathbb{X})\delta_{i,j}=\sum^{d}_{k=1}(-1)^{i+k}(\det\mathbb{X}_{k,i})X_{kj}.$
(3.16)

Denote by
$\mathbb{H}=\mathbb{E}_{r+1}(\lambda)_{(r-1,r,r+1),(r-1,r,r+1)}\in\mathbb{R}^{(r-2)\times(r-2)}$.
Then the matrix $\mathbb{E}_{r+1}(\lambda)$ can be partitioned as

$\displaystyle\mathbb{E}_{r+1}(\lambda)=\left(\begin{array}[]{ccc|ccc}&&&&&\\\
&\mathbb{H}&&b_{r-3}&0&0\\\ &&&-1&b_{r-2}&0\\\
\cline{1-3}\cr\lx@intercol\hfil\hfil\lx@intercol&a_{r-1}&-1&0\\\
\lx@intercol\hfil\hfil\lx@intercol&0&a_{r}&-1\\\
c_{1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r-1}&c_{r}&c_{r+1}\end{array}\right).$ (3.23)

Similarly, we have

$\displaystyle\mathbb{E}_{r}(\lambda)=\left(\begin{array}[]{ccc|cc}&&&&\\\
&\mathbb{H}&&b_{r-3}&0\\\ &&&-1&0\\\
\cline{1-3}\cr\lx@intercol\hfil\hfil\lx@intercol&a_{r-1}&-1\\\
c_{1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r-1}&c_{r}\end{array}\right),\ \
\mathbb{E}_{r-1}(\lambda)=\left(\begin{array}[]{ccc|c}&&&\\\ &\mathbb{H}&&0\\\
&&&-1\\\ \cline{1-3}\cr c_{1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r-1}\end{array}\right).$ (3.33)

The following simple identities will play an important role in our analysis

$\mathbb{E}_{r+1}(\lambda)_{r,(r,r+1)}=\mathbb{E}_{r}(\lambda)_{*,r},$ (3.34)

where for any $\mathbb{X}\in\mathbb{R}^{d\times d}$, we denote
$\mathbb{X}_{*,j}\in\mathbb{R}^{d\times(d-1)}$ the matrix by removing the
$j-$th column of $\mathbb{X}$. Similarly, we denote
$\mathbb{X}_{i,*}\in\mathbb{R}^{(d-1)\times d}$ the matrix by removing the
$i$-th row of $\mathbb{X}$.

The following elementary lemma is useful in our analysis.

######  Lemma 3.1.

For any $r\geq 3$ and $1\leq j\leq r-2$, we have

$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{(r-2,r,r+1),(j,r,r+1)}]$
$\displaystyle=-a_{r-1}\det[\mathbb{E}_{r-1}(\lambda)_{r-1,j}]$
$\displaystyle=a_{r-1}c_{r-1}\det[\mathbb{E}_{r-1}(\lambda)_{r-2,j}],$ (3.35)
$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{(r-1,r,r+1),(j,r,r+1)}]$
$\displaystyle=-\det[\mathbb{E}_{r}(\lambda)_{r,j}]$
$\displaystyle=c_{r}\det[\mathbb{E}_{r}(\lambda)_{r-1,j}].$ (3.36)

###### Proof.

We only prove (3.35). The identity (3.36) can be proved similarly. By (3.34),

$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{(r-2,r,r+1),(j,r,r+1)}]$
$\displaystyle=$ $\displaystyle\det[\mathbb{E}_{r}(\lambda)_{(r-2,r),(j,r)}]$
$\displaystyle=$ $\displaystyle\det\left(\begin{array}[]{ccc|c}&&&\\\
&\mathbb{H}_{r-2,j}&&\\\ &&&b_{r-3}\\\
\cline{1-3}\cr&&\lx@intercol\hfil\hfil\lx@intercol&a_{r-1}\end{array}\right)$
$\displaystyle=$ $\displaystyle a_{r-1}\det[\mathbb{H}_{r-2,j}].$

On the other hand, it is easy to see from (3.34) that

$\displaystyle\det[\mathbb{E}_{r-1}(\lambda)_{r-1,j}]=-\det\mathbb{H}_{r-2,j},\
\ \det[\mathbb{E}_{r-1}(\lambda)_{r-2,j}]=c_{r-1}\det[\mathbb{H}_{r-2,j}].$

This completes the proof. ∎

To proceed, we note that

$\mathbb{E}_{2}(\lambda)=\left(\begin{array}[]{cc}a_{1}&-1\\\
1&-1\end{array}\right),\ \ \ \
\mathbb{E}_{3}(\lambda)=\left(\begin{array}[]{ccc}a_{1}&-1&0\\\ 0&a_{2}&-1\\\
1&-1&1\end{array}\right).$ (3.38)

######  Lemma 3.2.

Let $\varphi_{r}(\lambda)=\det\mathbb{E}_{r+1}(\lambda)$ be the determinant of
$\mathbb{E}_{r+1}(\lambda)$. Then $\varphi_{1}(\lambda)=1-a_{1}$,
$\varphi_{2}(\lambda)=1-a_{1}+a_{1}a_{2}$, and

$\displaystyle\varphi_{r}(\lambda)=\varphi_{r-1}(\lambda)+a_{r}a_{r-1}\varphi_{r-2}(\lambda),\quad
r\geq 3.$ (3.39)

Moreover, $\varphi_{r}(\lambda)=P_{r}(-\lambda)$, where $P_{r}(\lambda)$ is
the numerator of $[r/r]$ Padé approximation of $e^{\lambda}$.

###### Proof.

The determinants of $\mathbb{E}_{r+1}(\lambda)$ for $r=1,2$ follow easily from
(3.38). Since by (3.15), $a_{r}=-b_{r-2}$, we obtain by adding the $r$-th row
to the $(r-2)$-th row of $\mathbb{E}_{r+1}(\lambda)$ in (3.23) and then
expanding the determinant by the $r$-th row that

$\displaystyle\varphi_{r}(\lambda)$ $\displaystyle=$
$\displaystyle\det\left(\begin{array}[]{ccc|ccc}&&&&&\\\
&\mathbb{H}&&b_{r-3}&0&0\\\ &&&-1&0&-1\\\
\cline{1-3}\cr\lx@intercol\hfil\hfil\lx@intercol&a_{r-1}&-1&0\\\
\lx@intercol\hfil\hfil\lx@intercol&0&a_{r}&-1\\\
c_{1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r-1}&c_{r}&c_{r+1}\end{array}\right)$
$\displaystyle=$
$\displaystyle\det\mathbb{E}_{r}(\lambda)+a_{r}\det\left(\begin{array}[]{ccc|cc}&&&&\\\
&\mathbb{H}&&b_{r-3}&0\\\ &&&-1&-1\\\
\cline{1-3}\cr\lx@intercol\hfil\hfil\lx@intercol&a_{r-1}&0\\\
c_{1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r-1}&c_{r}\end{array}\right)$ $\displaystyle=$
$\displaystyle\det\mathbb{E}_{r}(\lambda)+a_{r}a_{r-1}\det\mathbb{E}_{r-1}(\lambda),$

where in the last equality we have expanded the determinant by the $(r-1)$-th
row. Since $a_{1}=\lambda/2$, $a_{2}=\lambda/6$, we use Lemma 2.3 to conclude
$\varphi_{r}(\lambda)=P_{r}(-\lambda)$, where $P_{r}(\lambda)$ is the
numerator of $[r/r]$ Padé approximation of $e^{\lambda}$. ∎

######  Theorem 3.1.

The matrix $\varphi_{r}(\tau_{n}\mathbb{D})$ is invertible. The solution of
(3.5) $\mathbf{X}=(\mathbf{a}_{0}^{T},\cdots,\mathbf{a}_{r}^{T})^{T}$
satisfies that for $i=1,\cdots,r+1$,

$\displaystyle\mathbf{a}_{i-1}=(-1)^{r}\delta_{i,r+1}\mathbf{b}_{r}+\sum^{r}_{j=1}(\zeta_{j}\mathbb{I}+\tau_{n}\mathbb{D})^{-1}\left(\sum^{r+1}_{k=1}(-1)^{i+k+1}\frac{\phi_{ki}(-\zeta_{j})}{P^{\prime}_{r}(\zeta_{j})}\,\mathbf{b}_{k-1}\right),$

where $\phi_{ki}(\lambda)=\det[\mathbb{E}_{r+1}(\lambda)_{k,i}]$,
$k,i=1,\cdots,r+1$, are the minors of $\mathbb{E}_{r+1}(\lambda)$.

###### Proof.

By Lemma 3.2,
$\varphi_{r}(\lambda)=P_{r}(-\lambda)=(-1)^{r}\frac{r!}{(2r)!}(\lambda+\zeta_{1})\cdots(\lambda+\zeta_{r})$,
where $\zeta_{1},\cdots,\zeta_{r}\in\mathbb{C}$ are zeros of the diagonal Padé
numerator of type $[r/r]$ for $e^{\lambda}$. By Lemma 2.4, ${\rm
Re}(\zeta_{k})\leq-2$, $k=1,\cdots,r$. On the other hand, (1.2) implies that
the eigenvalues of $\mathbb{D}$ lie in the left half-plane. Thus the
eigenvalues of $\zeta_{k}\mathbb{I}+\tau_{n}\mathbb{D}$, $1\leq k\leq r$, lie
in the half-plane $\\{z\in\mathbb{C}:{\rm Re}(z)\leq-2\\}$. This shows
$\varphi_{r}(\tau_{n}\mathbb{D})$ is invertible.

Now by (3.16) we have

$\displaystyle[\det\mathbb{E}_{r+1}(\lambda)]\delta_{i,j}=\sum^{r+1}_{k=1}(-1)^{i+k}[\det\mathbb{E}_{r+1}(\lambda)_{k,i}]e_{kj}(\lambda),$

where $e_{kj}(\lambda)$ is the $(k,j)$ element of $\mathbb{E}_{r+1}(\lambda)$.
By replacing $\lambda$ by $\tau_{n}\mathbb{D}$ in above equality, we have

$\varphi_{r}(\tau_{n}\mathbb{D})\delta_{i,j}=\sum^{r+1}_{k=1}(-1)^{i+k}\phi_{ki}(\tau_{n}\mathbb{D})e_{kj}(\tau_{n}\mathbb{D}).$
(3.42)

From (3.5) we have

$\displaystyle\sum^{r+1}_{j=1}e_{ij}(\tau_{n}\mathbb{D})\mathbf{a}_{j-1}=\mathbf{b}_{i-1},\
\ \ i=1,\cdots,r+1.$

Thus multiplying (3.42) by $\mathbf{a}_{j-1}$ and summing over $j$ from $1$ to
$r+1$, we obtain

$\displaystyle\varphi_{r}(\tau_{n}\mathbb{D})\mathbf{a}_{i-1}$
$\displaystyle=$
$\displaystyle\sum^{r+1}_{k=1}(-1)^{i+k}\phi_{ki}(\tau_{n}\mathbb{D})\cdot\sum^{r+1}_{j=1}e_{kj}(\tau_{n}\mathbb{D})\mathbf{a}_{j-1}$
(3.43) $\displaystyle=$
$\displaystyle\sum^{r+1}_{k=1}(-1)^{i+k}\phi_{ki}(\tau_{n}\mathbb{D})\mathbf{b}_{k-1}.$

Note that for $(k,i)\not=(r+1,r+1)$, $\phi_{ki}(\lambda)\in P^{m}$, $m\leq
r-1$, by (2.22) we have

$\displaystyle\frac{\phi_{ki}(\tau_{n}\mathbb{D})}{\varphi_{r}(\tau_{n}\mathbb{D})}=\sum^{r}_{j=1}-\frac{\phi_{ki}(-\zeta_{j})}{P_{r}^{\prime}(\zeta_{j})}(\zeta_{j}\mathbb{I}+\tau_{n}\mathbb{D})^{-1}.$

For $(k,i)=(r+1,r+1)$, we have from (3.14) that
$\det[E_{r+1}(\lambda)_{r+1,r+1}]=\prod_{j=1}^{r}a_{j}=\frac{r!}{(2r)!}\lambda^{r}$.
Thus $\phi_{r+1,r+1}(\lambda)-(-1)^{r}\varphi_{r}(\lambda)\in P^{r-1}$, by
using (2.22) again we obtain

$\displaystyle\frac{\phi_{r+1,r+1}(\tau_{n}\mathbb{D})}{\varphi_{r}(\tau_{n}\mathbb{D})}=(-1)^{r}\mathbb{I}+\sum^{r}_{j=1}-\frac{\phi_{r+1,r+1}(-\zeta_{j})}{P_{r}^{\prime}(\zeta_{j})}(\zeta_{j}\mathbb{I}+\tau_{n}\mathbb{D})^{-1}.$

This completes the proof of the theorem. ∎

We remark that (3.43) can also be proved by using an abstract result in Brown
[4, Theorem 2.19 and Corollary 2.21] where linear algebra when matrix elements
are defined over a space of commuting matrices are studied.

From this theorem we know that the discrete problem (2.2) can be solved by
solving $r(r+1)$ linear systems of equations of order $M\times M$ in parallel
once all minors of $\mathbb{E}_{r+1}(\lambda)$ at
$\lambda=-\zeta_{j},j=1,\cdots,r$, are known. In the following, we will find
recursive formulas to computing these minors.

Let $\mathbb{G}_{r}(\lambda)=\mathbb{E}_{r+1}(\lambda)_{r,r+1}$, $r\geq 1$.
The determinants of $\mathbb{G}_{r}(\lambda)$ for $r=1,2$ can be calculated by
(3.38). We have the following lemma for $\det\mathbb{G}_{r}(\lambda)$ for
$r\geq 3$.

######  Lemma 3.3.

For $r\geq 3$, we have

$\displaystyle\det\mathbb{G}_{r}(\lambda)=\det\mathbb{G}_{r-1}(\lambda)-a_{r-1}b_{r-2}\det\mathbb{G}_{r-2}(\lambda)+c_{r}\prod^{r-1}_{k=1}a_{k}.$

###### Proof.

By definition and the partition in (3.23), we know that

$\displaystyle\mathbb{G}_{r}(\lambda)=\left(\begin{array}[]{ccc|cc}&&&&\\\
&\mathbb{H}&&b_{r-3}&0\\\ &&&-1&b_{r-2}\\\
\cline{1-3}\cr\lx@intercol\hfil\hfil\lx@intercol&a_{r-1}&-1\\\
c_{1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r-1}&c_{r}\end{array}\right).$

By expanding the determinant by the $(r-1)$-th row and use (3.34), we obtain

$\displaystyle\det\mathbb{G}_{r}(\lambda)$ $\displaystyle=$
$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{(r-1,r),(r,r+1)}]+a_{r-1}\det\left(\begin{array}[]{ccc|c}&&&\\\
&\mathbb{H}&&\\\ &&&b_{r-2}\\\ \cline{1-3}\cr c_{1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r}\end{array}\right)$ $\displaystyle=$
$\displaystyle\det[\mathbb{E}_{r}(\lambda)_{r-1,r}]+a_{r-1}(c_{r}\det\mathbb{H}-b_{r-2}\det[\mathbb{E}_{r+1}(\lambda)_{(r-2,r-1.r),(r-1,r,r+1)}])$
$\displaystyle=$
$\displaystyle\det\mathbb{G}_{r-1}(\lambda)+c_{r}\prod^{r-1}_{k=1}a_{k}-a_{r-1}b_{r-2}\det[\mathbb{E}_{r}(\lambda)_{(r-2,r-1),(r-1,r)}]$
$\displaystyle=$
$\displaystyle\det\mathbb{G}_{r-1}(\lambda)+c_{r}\prod^{r-1}_{k=1}a_{k}-a_{r-1}b_{r-2}\det[\mathbb{E}_{r-1}(\lambda)_{r-2,r-1}]$
$\displaystyle=$
$\displaystyle\det\mathbb{G}_{r-1}(\lambda)-a_{r-1}b_{r-2}\det\mathbb{G}_{r-2}(\lambda)+c_{r}\prod^{r-1}_{k=1}a_{k},$

where we have used the fact that $\det\mathbb{H}=\Pi^{r-2}_{k=1}a_{k}$ and
expanded the determinant by the last column in the second equality. This
completes the proof. ∎

The minors of $\mathbb{G}_{r}(\lambda)$ for $r=1,2$ can be computed directly
by (3.38). The following lemma gives the recursive formulas for some of the
minors of $\mathbb{G}_{r}(\lambda)$ which will be used to compute the minors
of $\mathbb{E}_{r+1}(\lambda)$.

######  Lemma 3.4.

For $r\geq 3$, we have

$\displaystyle\det[\mathbb{G}_{r}(\lambda)_{i,r}]$ $\displaystyle=$
$\displaystyle(-1)^{r-i-1}\frac{a_{1}\cdots a_{r-1}}{a_{1}\cdots
a_{i}}\det\mathbb{G}_{i}(\lambda),\ \ 1\leq i\leq r-1,$
$\displaystyle\det[\mathbb{G}_{r}(\lambda)_{r-1,j}]$ $\displaystyle=$
$\displaystyle\det[\mathbb{E}_{r}(\lambda)_{r-1,j}]-b_{r-2}\det[\mathbb{G}_{r-1}(\lambda)_{r-2,j}],\
\ 1\leq j\leq r-1,$

###### Proof.

For $i=1,\cdots,r-2$, by definition, we have

$\displaystyle\det[\mathbb{G}_{r}(\lambda)_{i,r}]$ $\displaystyle=$
$\displaystyle\det\left(\begin{array}[]{ccccccccc}a_{1}&-1&b_{1}&&&&&&\\\
&\ddots&\ddots&\ddots&&&&\\\ &&a_{i-1}&-1&b_{i-1}&&&&\\\
&&&0&a_{i+1}&-1&b_{i+1}&&\\\ &&&&\ddots&\ddots&\ddots&\ddots&\\\
&&&&&0&a_{r-3}&-1&b_{r-3}\\\ &&&&&&0&a_{r-2}&-1\\\ &&&&&&&0&a_{r-1}\\\
c_{1}&\cdots&\cdots&\cdots&\cdots&\cdots&c_{r-3}&c_{r-2}&c_{r-1}\\\
\end{array}\right)$ $\displaystyle=$
$\displaystyle(-a_{r-1})\cdots(-a_{i+1})\det\left(\begin{array}[]{ccccc}a_{1}&-1&b_{1}&&\\\
&\ddots&\ddots&\ddots&\\\ &&a_{i-2}&-1&b_{i-2}\\\ &&&a_{i-1}&-1\\\
c_{1}&\cdots&\cdots&c_{i-1}&c_{i}\\\ \end{array}\right)$ $\displaystyle=$
$\displaystyle(-1)^{r-i-1}(\prod^{r-1}_{k=i+1}a_{k})\det\mathbb{G}_{i}(\lambda).$

Finally, for $\mathbb{G}(\lambda)_{i,r}$, $i=r-1$, we have by the definition
and using the first identity in (3.34)

$\displaystyle\mathbb{G}_{r}(\lambda)_{r-1,r}=\mathbb{E}_{r+1}(\lambda)_{(r-1,r),(r,r+1)}=\mathbb{E}_{r}(\lambda)_{r-1,r}=\mathbb{G}_{r-1}(\lambda).$

This shows the first equality of the lemma. To show the second equality, for
any $1\leq j\leq r-2$, we have by using the partition (3.23) that

$\displaystyle\mathbb{G}_{r}(\lambda)_{r-1,j}=\left(\begin{array}[]{cccccc|cc}&&&&&&&\\\
&&\mathbb{H}_{*,j}&&&&b_{r-3}&0\\\ &&&&&&-1&b_{r-2}\\\
\cline{1-6}\cr\lx@intercol\hfil\hfil\lx@intercol&&\\\
c_{1}&\cdots&c_{j-1}&c_{j+1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r-1}&c_{r}\end{array}\right).$

Thus by expanding the determinant by the last column, we have by using (3.34),
(3.36) and the definition of $\mathbb{G}_{r-1}(\lambda)$ that

$\displaystyle\det[\mathbb{G}_{r}(\lambda)_{r-1,j}]$ $\displaystyle=$
$\displaystyle
c_{r}\det[\mathbb{E}_{r+1}(\lambda)_{(r-1,r,r+1),(j,r,r+1)}]-b_{r-2}\det[\mathbb{E}_{r+1}(\lambda)_{(r-2,r-1,r),(j,r,r+1)}]$
$\displaystyle=$
$\displaystyle\det[\mathbb{E}_{r}(\lambda)_{r-1,j}]-b_{r-2}\det[\mathbb{E}_{r}(\lambda)_{(r-2,r-1),(j,r)}]$
$\displaystyle=$
$\displaystyle\det[\mathbb{E}_{r}(\lambda)_{r-1,j}]-b_{r-2}\det[\mathbb{G}_{r-1}(\lambda)_{r-2,j}].$

This completes the proof. ∎

The minors of $\mathbb{E}_{r+1}(\lambda)$ for $r=1,2$ can be easily computed
from (3.38). The following theorem gives the recursive formulas for computing
all minors of $\mathbb{E}_{r+1}(\lambda)$ for $r\geq 3$.

######  Theorem 3.2.

Let $r\geq 3$, we have
$1^{\circ}$ For $i=1,\cdots,r-2$,

$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{i,j}]=\det[\mathbb{E}_{r}(\lambda)_{i,j}]+a_{r}a_{r-1}\det[\mathbb{E}_{r-1}(\lambda)_{i,j}],\
\ 1\leq j\leq r-2,$
$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{i,j}]=(-1)^{j-i-1}\frac{a_{1}\cdots
a_{j-1}}{a_{1}\cdots a_{i}}\det\mathbb{G}_{i}(\lambda),\ \ j=r-1,r,r+1.$

$2^{\circ}$ For $i=r-1$, we have

$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{r-1,j}]=-c_{r+1}a_{r}\det[\mathbb{E}_{r}(\lambda)_{r,j}]+\det[\mathbb{G}_{r}(\lambda)_{r-1,j}],\
\ 1\leq j\leq r-1,$
$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{r-1,j}]=(-1)^{j-i-1}\frac{a_{1}\cdots
a_{j-1}}{a_{1}\cdots a_{i}}\det\mathbb{G}_{i}(\lambda),\ \ j=r,r+1.$

$3^{\circ}$ For $i=r$, we have

$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{r,j}]=\det[\mathbb{E}_{r}(\lambda)_{r-1,j}]-a_{r-1}b_{r-2}\det[\mathbb{E}_{r-1}(\lambda)_{r-2,j}],\
\ 1\leq j\leq r-2,$
$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{r,j}]=c_{r+1}(-1)^{i-j}\prod^{j-1}_{k=1}a_{k},\
\ j=r-1,r,$
$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{r,r+1}]=\det\mathbb{G}_{r-1}(\lambda)-a_{r-1}b_{r-2}\det\mathbb{G}_{r-2}(\lambda)+c_{r}\prod^{r-1}_{k=1}a_{k}.$

$4^{\circ}$ For $i=r+1$, we have

$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{r+1,j}]=-\det[\mathbb{E}_{r}(\lambda)_{r,j}]-a_{r-1}b_{r-2}\det[\mathbb{E}_{r-1}(\lambda)_{r-1,j}],\
\ 1\leq j\leq r-2,$
$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{r+1,j}]=(-1)^{i-j}\prod^{j-1}_{k=1}a_{k},\
\ j=r-1,r,r+1.$

###### Proof.

The proof is divided into 4 steps.
Step 1. The first equality in $1^{\circ}$ can be proved by the same argument
as that in Lemma 3.2. Here we omit the details. We only prove the second
equality in $1^{\circ}$ when $j=r+1$. The other cases can be proved similarly.
By the partition in (3.23), we know that for $1\leq i\leq r-2$,

$\displaystyle\mathbb{E}_{r+1}(\lambda)_{i,r+1}=\left(\begin{array}[]{ccc|cc}&&&&\\\
&\mathbb{H}_{i,*}&&\ \ \ \ \ \mathbb{F}_{i,*}&\hskip-5.69054pt\\\ &&&&\\\
\cline{1-3}\cr\cline{4-5}\cr\lx@intercol\hfil\hfil\lx@intercol&a_{r-1}&-1\\\
\lx@intercol\hfil\hfil\lx@intercol&0&a_{r}\\\ c_{1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r-1}&c_{r}\end{array}\right),\ \ \ \
\mathbb{F}=\left(\begin{array}[]{ccc}0&0&0\\\ \vdots&\ \vdots&\vdots\\\
0&0&0\\\ b_{r-3}&0&0\\\ -1&b_{r-2}&0\\\ \end{array}\right),$

where $\mathbb{F}\in\mathbb{R}^{(r-2)\times 3}$. By expanding the determinant
first by the $r$-th and then by the $(r-1)$-th row, we know by (3.34) that

$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{i,r+1}]$ $\displaystyle=$
$\displaystyle
a_{r}a_{r-1}\det[\mathbb{E}_{r+1}(\lambda)_{(i,r-1,r),(r-1,r,r+1)}]$
$\displaystyle=$ $\displaystyle
a_{r}a_{r-1}\det[\mathbb{E}_{r}(\lambda)_{(i,r-1),(r-1,r)}]$ $\displaystyle=$
$\displaystyle a_{r}a_{r-1}\det[\mathbb{G}_{r-1}(\lambda)_{i,r-1}].$

This shows the second equality in $1^{\circ}$ by the first identity in Lemma
3.4.

Step 2. We only prove the first equality in $2^{\circ}$ when $1\leq j\leq
r-2$. The other cases can be proved similarly. By the partition in (3.23), we
have for $1\leq j\leq r-2$,

$\displaystyle\mathbb{E}_{r+1}(\lambda)_{r-1,j}=\left(\begin{array}[]{ccc|ccc}&&&&&\\\
&\mathbb{H}_{*,j}&&b_{r-3}&0&0\\\ &&&-1&b_{r-2}&0\\\
\cline{1-3}\cr\lx@intercol\hfil\hfil\lx@intercol&0&a_{r}&-1\\\
c_{1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r-1}&c_{r}&c_{r+1}\end{array}\right)$

Expanding the determinant by the last column, we obtain

$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{r-1,j}]$ $\displaystyle=$
$\displaystyle
c_{r+1}a_{r}\det[\mathbb{E}_{r+1}(\lambda)_{(r-1,r,r+1),(j,r,r+1)}]+\det[\mathbb{E}_{r+1}(\lambda)_{(r-1,r),(j,r+1)}]$
$\displaystyle=$ $\displaystyle-
c_{r+1}a_{r}\det[\mathbb{E}_{r}(\lambda)_{r,j}]+\det[\mathbb{G}_{r}(\lambda)_{r-1,j}].$

This shows the first equality in $2^{\circ}$ for $1\leq j\leq r-2$.

Step 3. The second equality in $3^{\circ}$ can be easily proved. The last
equality is shown in Lemma 3.3 since
$\mathbb{E}_{r+1}(\lambda)_{r,r+1}=\mathbb{G}_{r}(\lambda)$ by definition. To
show the first equality in $3^{\circ}$, we again use the partition in (3.23)
to obtain for $1\leq j\leq r-2$,

$\displaystyle\mathbb{E}_{r+1}(\lambda)_{r,j}=\left(\begin{array}[]{cccccc|ccc}&&&&&&&&\\\
&&\mathbb{H}_{*,j}&&&&b_{r-3}&0&0\\\ &&&&&&-1&b_{r-2}&0\\\
\cline{1-6}\cr\lx@intercol\hfil\hfil\lx@intercol&a_{r-1}&-1&0\\\
c_{1}&\cdots&c_{j-1}&c_{j+1}&\cdots&\lx@intercol\hfil
c_{r-2}\hfil\lx@intercol&c_{r-1}&c_{r}&c_{r+1}\end{array}\right).$

By expanding the determinant successively by the last columns, we have

$\displaystyle\det[\mathbb{E}_{r+1}(\lambda)_{r,j}]$ $\displaystyle=$
$\displaystyle-
c_{r+1}\det[\mathbb{E}_{r+1}(\lambda)_{(r-1,r,r+1),(j,r,r+1}]-c_{r+1}b_{r-2}\det[\mathbb{E}_{r+1}(\lambda)_{(r-2,r,r+1),(j,r,r+1)}].$

This shows the first equality in $3^{\circ}$ by Lemma 3.1.

Step 4. The first equality in $4^{\circ}$ can be proved by the same argument
as that Step 3. The second equality can be easily proved. Here we omit the
details. ∎

This theorem indicates that all minors of $\mathbb{E}_{r+1}(\lambda)$ can be
computed once one knows all minors of $\mathbb{E}_{m}(\lambda)$, $1\leq m\leq
r$, and the minors $\det[\mathbb{G}_{r}(\lambda)_{r-1,j}]$, $1\leq j\leq r-1$,
which can be computed by Lemma 3.4 recursively based on the information of the
minors $E_{m}(\lambda)$, $1\leq m\leq r$.

The following lemma indicates that the nodal values of the solution to (2.2)
depends only on the coefficient $\mathbf{a}_{0}$.

######  Lemma 3.5.

Let $\mathbf{Y}_{r}(t)$, $r\geq 1$, be the solution of the problem (2.2). Then

$\displaystyle\mathbf{Y}_{r}(t_{n+1})=\mathbf{Y}_{r}(t_{n})+\tau_{n}\mathbb{D}\mathbf{a}_{0}+\tau_{n}\mathbf{R}_{0},\
\ \ n=1,\cdots,N-1.$

###### Proof.

We integrate (2.2) over $I_{n}$ and use the orthogonality of Legendre
polynomials to obtain

$\displaystyle\mathbf{Y}_{r}(t_{n+1})$
$\displaystyle=\mathbf{Y}_{r}(t_{n})+\int_{I_{n}}\mathbb{D}\,\left(\sum_{j=0}^{r-1}\mathbf{a}_{j}\widetilde{L}_{j}(t)\right)\,dt+\int_{I_{n}}\sum_{j=0}^{r-1}\mathbf{R}_{j}\widetilde{L}_{j}(t)\,dt,$
$\displaystyle=\mathbf{Y}_{r}(t_{n})+\tau_{n}\mathbb{D}\mathbf{a}_{0}+\tau_{n}\mathbf{R}_{0}.$

This completes proof. ∎

By Theorem 3.1, we have then

$\displaystyle\mathbf{Y}_{r}(t_{n+1})=\mathbf{Y}_{r}(t_{n})+\sum^{r}_{j=1}(\zeta_{j}\mathbb{I}+\tau_{n}\mathbb{D})^{-1}\left(\sum^{r+1}_{k=1}(-1)^{k}\frac{\phi_{k1}(-\zeta_{j})}{P_{r}^{\prime}(\zeta_{j})}(\tau_{n}\mathbb{D})\,\mathbf{b}_{k-1}\right)+\tau_{n}\mathbf{R}_{0}.$

This leads to the following parallel algorithm to compute the nodal values of
the solution $\mathbf{Y}_{r}$ to the problem (2.2).

######  Algorithm 3.1.

Given $\mathbf{Y}_{r}(t_{0})=\mathbf{Y}_{0}$. For $n=1,\cdots,N-1$, do the
following.
$1^{\circ}$ Compute $\mathbf{v}_{j}\in\mathbb{R}^{M}$, $j=1,\cdots,r$, in
parallel, where

$\displaystyle\mathbf{v}_{j}=\sum^{r+1}_{k=1}(-1)^{k}\frac{\phi_{k1}(-\zeta_{j})}{P_{r}^{\prime}(\zeta_{j})}(\tau_{n}\mathbb{D})\,\mathbf{b}_{k-1}.$

$2^{\circ}$ Solve
$(\tau_{n}\mathbb{D}+\zeta_{j}\mathbb{I})\mathbf{w}_{j}=\mathbf{v}_{j}$,
$j=1,\cdots,r$, in parallel.
$3^{\circ}$ Compute

$\displaystyle\mathbf{Y}_{r}(t_{n+1})=\mathbf{Y}_{r}(t_{n})+\sum_{j=0}^{r}\mathbf{w}_{j}+\tau_{n}\mathbf{R}_{0}.$

The following parallel-in-time algorithm computes the solution of the problem
(2.2) inside each time interval.

######  Algorithm 3.2.

Given $\mathbf{Y}_{r}(t_{0})=\mathbf{Y}_{0}$.
$1^{\circ}$ Call Algorithm 3.1 to obtain $\mathbf{Y}_{r}(t_{n})$,
$n=1,\ldots,N$.
$2^{\circ}$ Compute the coefficients $\mathbf{a}_{1},\cdots,\mathbf{a}_{r}$ of
$\mathbf{Y}_{r}$ in each time interval $I_{n}$, $n=1,\cdots,N-1$, in parallel
as follows.
(i) Compute $\mathbf{v}_{ij}\in\mathbb{R}^{M}$, $i=2,\cdots,r+1$, $j=1,\cdots
r$, in parallel, where

$\displaystyle\mathbf{v}_{ij}=\sum^{r+1}_{k=1}(-1)^{i+k+1}\frac{\phi_{ki}(-\zeta_{j})}{P_{r}^{\prime}(\zeta_{j})}\,\mathbf{b}_{k-1}.$

(ii) Solve
$(\tau_{n}\mathbb{D}+\zeta_{j}\mathbb{I})\mathbf{w}_{ij}=\mathbf{v}_{ij}$,
$i=2,\cdots,r+1,j=1,\cdots,r$, in parallel.
(iii) Compute
$\mathbf{a}_{i-1}=(-1)^{r}\delta_{i,r+1}\mathbf{b}_{r}+\sum_{j=1}^{r}\mathbf{w}_{ij}$,
$i=2,\cdots,r+1$, in parallel.

######  Remark 3.1.

In [12], it is observed that the zeros of $P_{r}(z)$ come in complex conjugate
pairs if they are complex. If $\zeta_{j^{\prime}}=\bar{\zeta}_{j}$,
$j,j^{\prime}=1,\cdots,r$, then
$\mathbf{v}_{j}=\bar{\mathbf{v}}_{j^{\prime}}$, and

$\displaystyle\mathbf{w}_{j}+\mathbf{w}_{j^{\prime}}=\frac{\mathbf{v}_{j}}{\tau_{n}\mathbb{D}+\zeta_{j}\mathbb{I}}+\frac{\bar{\mathbf{v}}_{j}}{\tau_{n}\mathbb{D}+\bar{\zeta}_{j}\mathbb{I}}=2\,{\rm
Re}\,\left[\frac{\mathbf{v}_{j}}{\tau_{n}\mathbb{D}+\zeta_{j}\mathbb{I}}\right].$

Thus one need only to solve $k$ complex matrix problems instead of $2k$ in
Algorithm 3.1 ($2^{\circ}$) and $kr$ complex matrix problems instead of $2kr$
in Algorithm 3.2 ($2^{\circ}$), where $2k$, $0\leq k\leq r/2$, are the number
of complex zeros of $P_{r}(z)$.

######  Remark 3.2.

If $\zeta=a+{\bf i}b$ is a complex zero of $P_{r}(z)$, then $a\leq-2$ by Lemma
2.3. By Remark 3.1, without loss of generality, we can choose one of the zeros
such that $b<0$. Let $\mathbf{w}=\mathbf{w}_{1}+{\bf i}\mathbf{w}_{2}$,
$\mathbf{v}=\mathbf{v}_{1}+{\bf i}\mathbf{v}_{2}$, where
$\mathbf{v}_{i},\mathbf{w}_{i}\in\mathbb{R}^{M}$, $i=1,2$, satisfy
$(\tau_{n}\mathbb{D}+\zeta\mathbb{I})\mathbf{w}=\mathbf{v}$. Then

$\displaystyle\widetilde{\mathbb{D}}\left(\begin{array}[]{c}\mathbf{w}_{1}\\\
\mathbf{w}_{2}\\\
\end{array}\right):=\left(\begin{array}[]{cc}\tau_{n}\mathbb{D}+a\mathbb{I}&-b\,\mathbb{I}\\\
-b\,\mathbb{I}&-(\tau_{n}\mathbb{D}+a\mathbb{I})\\\
\end{array}\right)\left(\begin{array}[]{c}\mathbf{w}_{1}\\\ \mathbf{w}_{2}\\\
\end{array}\right)=\left(\begin{array}[]{c}\mathbf{v}_{1}\\\
-\mathbf{v}_{2}\\\ \end{array}\right).$

Let $\mathbb{F}={\rm
diag}(\tau_{n}\mathbb{D}+(a+b)\,\mathbb{I},\tau_{n}\mathbb{D}+(a+b)\,\mathbb{I})\in\mathbb{R}^{2M\times
2M}$ be the diagonal matrix. It is shown in Chen et al [5, Lemma 4.1] that the
condition number $\kappa(\mathbb{F}^{-1}\widetilde{\mathbb{D}})\leq\sqrt{2}$.
Therefore, the complex system
$(\tau_{n}\mathbb{D}+\zeta\mathbb{I})\mathbf{w}=\mathbf{v}$ can be efficiently
solved if one has the efficient solver for the real matrix
$\tau_{n}\mathbb{D}+(a+b)\,\mathbb{I}$, where $a+b\leq-2$. Notice that the
eigenvalues of $\mathbb{D}$ lie in the left half-plane due to the assumption
$\mathbb{D}+\mathbb{D}^{T}\leq 0$.

######  Remark 3.3.

For non-standard ODE systems of the form,

$\displaystyle\mathbb{M}\mathbf{Y}^{\prime}=\mathbb{D}\mathbf{Y}+\mathbf{R}\ \
\text{ in }(0,T),\quad\mathbf{Y}(0)=\mathbf{Y}_{0},$ (3.53)

one can use the transformations
$\widetilde{\mathbf{Y}}=\mathbb{M}^{\frac{1}{2}}\mathbf{Y}$,
$\widetilde{\mathbb{D}}=\mathbb{M}^{-\frac{1}{2}}\mathbb{D}\mathbb{M}^{-\frac{1}{2}}$,
$\widetilde{\mathbf{R}}=\mathbb{M}^{-\frac{1}{2}}\mathbf{R}$ to transform the
problem (3.53) to (1.1) and use above algorithms to solve the transformed
problem. This leads to the following algorithm which is similar to Algorithm
3.1 to find the nodal values of the solution of the continuous time Galerkin
method for solving (3.53). A similar algorithm to Algorithm 3.2 can also be
formulated.

######  Algorithm 3.3.

Given $\mathbf{Y}_{r}(t_{0})=\mathbf{Y}_{0}$. For $n=1,\cdots,N-1$, do the
following.
$1^{\circ}$ Compute $\mathbf{v}_{j}\in\mathbb{R}^{M}$, $j=1,\cdots,r$, in
parallel, where

$\displaystyle\mathbf{v}_{j}=\sum^{r+1}_{k=1}(-1)^{k}\frac{\phi_{k1}(-\zeta_{j})}{P_{r}^{\prime}(\zeta_{j})}(\tau_{n}\mathbb{D})\mathbb{M}^{-1}\,\mathbf{b}_{k-1}.$

$2^{\circ}$ Solve
$(\tau_{n}\mathbb{D}+\zeta_{j}\mathbb{M})\mathbf{w}_{j}=\mathbf{v}_{j}$,
$j=1,\cdots,r$, in parallel.
$3^{\circ}$ Compute

$\displaystyle\mathbf{Y}_{r}(t_{n+1})=\mathbf{Y}_{r}(t_{n})+\sum_{j=1}^{r}\mathbf{w}_{j}+\tau_{n}\mathbb{M}^{-1}\mathbf{R}_{0}.$

To conclude this section, we prove the following theorem for finding the nodal
values of the solution (2.2) which extends (1.3) for solving the ODE system
(1.1) when $\mathbf{R}=\bf{0}$.

######  Theorem 3.3.

Let $\mathbf{Y}_{r}\in\mathbf{V}^{r}_{\tau}$, $r\geq 1$, be the solution of
the problem (2.2). Then for $n=1,\cdots,N-1$,

$\displaystyle\mathbf{Y}_{r}(t_{n+1})=\frac{P_{r}(\tau_{n}\mathbb{D})}{P_{r}(-\tau_{n}\mathbb{D})}\mathbf{Y}_{r}(t_{n})+\sum^{r}_{k=1}(-1)^{k+1}\frac{\phi_{k1}(\tau_{n}\mathbb{D})}{P_{r}(-\tau_{n}\mathbb{D})}\,\mathbf{b}_{k-1}+\tau_{n}\mathbf{R}_{0},$

where $\phi_{k1}(\lambda)=\det[\mathbb{E}_{r+1}(\lambda)_{k,1}]$.

###### Proof.

By Lemma 3.2 and (3.43) we have

$\displaystyle
P_{r}(-\tau_{n}\mathbb{D})\mathbf{a}_{0}=\sum^{r+1}_{k=1}(-1)^{k+1}\phi_{k1}(\tau_{n}\mathbb{D})\mathbf{b}_{k-1}.$

Since $\mathbf{b}_{r+1}=\mathbf{Y}_{r}(t_{n})$, by Lemma 3.5,

$\displaystyle\mathbf{Y}_{r}(t_{n+1})$ $\displaystyle=$
$\displaystyle\left[\mathbb{I}+(-1)^{r}\frac{\tau_{n}\mathbb{D}\phi_{r+1,1}(\tau_{n}\mathbb{D})}{P_{r}(-\tau_{n}\mathbb{D})}\right]\mathbf{Y}_{r}(t_{n})$
$\displaystyle+\sum^{r}_{k=1}(-1)^{k+1}\phi_{k1}(\tau_{n}\mathbb{D})\mathbf{b}_{k-1}+\tau_{n}\mathbf{R}_{0}.$

Denote by
$\psi_{r}(\lambda)=(-1)^{r}\lambda\phi_{r+1,1}(\lambda)=(-1)^{r}\lambda\det[\mathbb{E}_{r+1}(\lambda)_{r+1,1}]$.
By Theorem 3.2, $4^{\circ}$, we know that $\psi_{r}(\lambda)$ satsifies

$\displaystyle\psi_{r}(\lambda)=\psi_{r-1}(\lambda)+\frac{\lambda^{2}}{4}\frac{1}{(2r-1)(2r-3)}\psi_{r-2}(\lambda),\
\ r\geq 3.$

On the other hand, by (3.38), we have
$\psi_{1}(\lambda)=\lambda,\psi_{2}(\lambda)=\lambda$. This implies by Lemma
2.3 that $\psi_{r}(\lambda)=P_{r}(\lambda)-P_{r}(-\lambda)$. Thus

$\displaystyle\mathbb{I}+(-1)^{r}\frac{\tau_{n}\mathbb{D}\phi_{r+1,1}(\tau_{n}\mathbb{D})}{P_{r}(-\tau_{n}\mathbb{D})}=\mathbb{I}+\frac{P_{r}(\tau_{n}\mathbb{D})-P_{r}(-\tau_{n}\mathbb{D})}{P_{r}(-\tau_{n}\mathbb{D})}=P_{r}(\tau_{n}\mathbb{D}).$

This completes the proof. ∎

##  4 The dissipative system

In this section, we propose an alternative way to compute the coefficients
$\mathbf{a}_{1},\cdots,\mathbf{a}_{r}$ of the solution $\mathbf{Y}_{r}$ of the
problem (2.2) when the ODE system (1.1) is dissipative
$\mathbb{D}+\mathbb{D}^{T}<0$. The algorithm is based on the block tridiagonal
structure of the matrix and is less expensive than the step $2^{\circ}$ in
Algorithm 3.2.

Let
$\widetilde{\mathbf{X}}=(\mathbf{a}_{1}^{T},\cdots,\mathbf{a}_{r}^{T})^{T}\in\mathbb{R}^{rM}$
and
$\widetilde{\mathbf{B}}=(\widetilde{\mathbf{b}}^{T}_{1},\mathbf{b}_{2}^{T},\cdots,\mathbf{b}_{r}^{T})^{T}$
with
$\widetilde{\mathbf{b}}_{1}=\mathbf{b}_{1}-\frac{1}{2}(\tau_{n}\mathbb{D})\mathbf{a}_{0}$.
It follows from (3.5) that $\widetilde{\mathbf{X}}\in\mathbb{R}^{rM}$
satisfies

$\widetilde{\mathbb{A}}\widetilde{\mathbf{X}}=\widetilde{\mathbf{B}},$ (4.1)

where

$\displaystyle\widetilde{\mathbb{A}}=\left(\begin{array}[]{cccccc}-\mathbb{I}&-\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{5}&\cdots&\cdots&\cdots&0\\\
\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{3}&-\mathbb{I}&-\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{7}&\cdots&\cdots&0\\\
\vdots&\ddots&\ddots&\ddots&\vdots&\vdots\\\
\cdots&\cdots&\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{2r-5}&-\mathbb{I}&-\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{2r-1}&0\\\
\cdots&\cdots&\cdots&\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{2r-3}&-\mathbb{I}&0\\\
\cdots&\cdots&\cdots&\cdots&\frac{\tau_{n}}{2}\mathbb{D}\frac{1}{2r-1}&-\mathbb{I}\\\
\end{array}\right).$

It is easy to see that
$\widetilde{\mathbb{A}}=\mathbb{E}_{r+1}(\tau_{n}\mathbb{D})_{r+1,1}$, let
$\phi_{r}(\lambda)=\det[\mathbb{E}_{r+1}(\lambda)_{r+1,1}]$. The goal of this
section is to show that $\phi_{r}(\tau_{n}\mathbb{D})\in\mathbb{R}^{M\times
M}$ is invertible so that the standard chasing algorithm for block tridiagonal
matrices (cf., e.g., Golub and Van Load [14, §4.5]) can be used to solve
(4.1).

We start with the following lemma.

######  Lemma 4.1.

Let $u_{n}\in P^{n}$, $n\geq 0$, such that (i) $u_{0}=0,u_{1}=1$ or (ii)
$u_{0}=1,u_{1}=1+A_{1}t$, $A_{1}\in\mathbb{R}$, and for $n\geq 2$,

$u_{n}(t)=(1+A_{n}t)u_{n-1}(t)-C_{n}t^{2}u_{n-2},$ (4.3)

where $A_{n}\in\mathbb{R},C_{n}>0$ for $n\geq 2$. Then we have, for $n\geq 2$,

$\displaystyle u_{n}(u_{n-1}+tu^{\prime}_{n-1})-tu_{n-1}u^{\prime}_{n}>0\ \ \
\mbox{in }\mathbb{R}\backslash\\{0\\}.$

###### Proof.

We set, for $n\geq 1$,

$\displaystyle
G_{n}(t,t^{\prime})=\frac{t^{\prime}u_{n}(t)u_{n-1}(t^{\prime})-tu_{n-1}(t)u_{n}(t^{\prime})}{t^{\prime}-t}.$

It is easy to see by (4.3) that $G_{1}(t,t^{\prime})=u_{0}$, and for $n\geq
2$,

$\displaystyle
G_{n}(t,t^{\prime})=u_{n-1}(t)u_{n-1}(t^{\prime})+C_{n}tt^{\prime}G_{n-1}(t,t^{\prime}).$

This implies easily

$\displaystyle
G_{n}(t,t^{\prime})=u_{n-1}(t)u_{n-1}(t^{\prime})+\sum^{n-2}_{k=0}C_{n}\cdots
C_{k+2}(tt^{\prime})^{n-1-k}u_{k}(t)u_{k}(t^{\prime}),$

where we have used $G_{1}(t,t^{\prime})=u_{0}(t)u_{0}(t^{\prime})$ in both
cases (i) and (ii). By letting $t^{\prime}\to t$, we obtain for $n\geq 2$,

$\displaystyle\lim_{t^{\prime}\to
t}G_{n}(t,t^{\prime})=u_{n}(t)[tu_{n-1}(t)]^{\prime}-tu_{n-1}(t)u_{n}^{\prime}(t)>0\
\ \mbox{in }\mathbb{R}\backslash\\{0\\},$

where we have used the condition $u_{1}=1$ in the case (i) and $u_{0}=1$ in
the case (ii). This completes the proof. ∎

The following theorem is the main result of this section.

######  Theorem 4.1.

Let $\mathbb{D}\in\mathbb{R}^{M\times M}$ satisfy
$\mathbb{D}+\mathbb{D}^{T}<0$, that is, $\mathbb{D}+\mathbb{D}^{T}$ is
negative definite. Then the matrix
$\phi_{r}(\tau_{n}\mathbb{D})\in\mathbb{R}^{M\times M}$, $r\geq 1$, $1\leq
n\leq N-1$, is invertible.

###### Proof.

We are going to show that all zeros of $\phi_{r}(\lambda)$ locate at the
imaginary axis $\\{z\in\mathbb{C}:{\rm Re}(z)=0\\}$. This implies easily
$\phi_{r}(\tau_{n}\mathbb{D})$ is invertible since the eigenvalues of
$\mathbb{D}$ lie in the left-half plane due to the dissipative property
$\mathbb{D}+\mathbb{D}^{T}<0$.

To study the property of the zeros of $\phi_{r}(\lambda)$, we denote
$\psi_{r}(\lambda)=(-1)^{r}\phi_{r}(\lambda)$. By (3.38) and Theorem 3.2
($4^{\circ}$) we know that $\psi_{1}(\lambda)=1,\psi_{2}(\lambda)=1$, and

$\psi_{r}(\lambda)=\psi_{r-1}(\lambda)+d_{r}\lambda^{2}\psi_{r-2}(\lambda),\ \
\ d_{r}=\frac{1}{4(2r-1)(2r-3)},\ \ \ \mbox{for }r\geq 3.$ (4.4)

We note that (4.4) is also valid for $r=2$ if we define $\psi_{0}(\lambda)=0$.
Set $t=\lambda^{2}$ and for $m\geq 0$, define
$f_{m}(t)=\psi_{2m}(\lambda),g_{m}=\psi_{2m+1}(\lambda)$. Then

$\displaystyle f_{m}(t)=g_{m-1}(t)+d_{2m}tf_{m-1}(t),\ \
f_{0}(t)=0,f_{1}(t)=1,$ $\displaystyle g_{m}(t)=f_{m}(t)+d_{2m+1}tg_{m-1}(t),\
\ g_{0}(t)=1,g_{1}(t)=1+d_{3}t.$

This implies that, for $m\geq 2$,

$\displaystyle f_{m}(t)=(1+A_{m}t)f_{m-1}(t)-C_{m}t^{2}f_{m-2}(t),$ (4.5)
$\displaystyle
g_{m}(t)=(1+\tilde{A}_{m}t)g_{m-1}(t)-\tilde{C}_{m}t^{2}f_{m-2}(t),$ (4.6)

where
$A_{m}=d_{2m}+d_{2m-1},C_{m}=d_{2m-1}d_{2m-2},\tilde{A}_{m}=d_{2m+1}+d_{2m},\tilde{C}_{m}=d_{2m}d_{2m-1}$.
We observe that $f_{m},g_{m}$ satisfy the same recurrence relation but with
different coefficients and initial values. In the following we will only prove
$f_{m}$, $m\geq 2$, has $m$ real zeros in $(-\infty,0)$ which then implies
that $\phi_{2m}(\lambda)=f_{m}(\lambda^{1/2})$ has all zeros on the imaginary
axis. The proof for $\phi_{2m+1}(\lambda)$, $m\geq 1$, is similar and we omit
the details.

We extend the argument in [24, §3.3 (4)] for orthogonal polynomials to show
that $f_{m},m\geq 2$, has $m$ zeros in $(-\infty,0)$ by using Sturm theorem
(cf., e.g., Perron [18, pp.7-9]) based on the recurrence formula (4.5). We
first note that if $f_{m}(t)=\sum^{m}_{k=0}\theta_{k}t^{k}$, then
$\theta_{m}>0$ since by (4.4) the leading coefficients of the polynomial
$\psi_{r}(\lambda)$ are positive, and $\theta_{0}=f_{m}(0)=1$ by (4.5). Now we
claim that

$f_{m}(t),f_{m-1}(t),\cdots,f_{1}(t)$ (4.7)

form a Sturmian sequence in $[-M,-\delta]$ for sufficiently large $M>0$ and
sufficiently small $\delta>0$ in the following sense. (i) $f_{1}(t)=1$ has no
zeros in $[-M,-\delta]$. (ii) $f_{m}(-M)f_{m}(-\delta)\not=0$ for $M\gg 1$ and
$\delta\ll 1$ since $\theta_{m}>0$ and $\theta_{0}=1$. (iii) If
$c\in[-M,-\delta]$ is a zero of $f_{k}(t)$, $1\leq k\leq m-1$, then
$f_{k+1}(c)f_{k-1}(c)<0$. In fact, By Lemma 4.1, $f_{k-1}(c)\not=0$. By (4.5),
$f_{k+1}(c)=-C_{m}c^{2}f_{k-1}(c)$. (iv) If $c\in[-M,-\delta]$ such that
$f_{m}(c)=0$, then $f_{m}^{\prime}(c)f_{m-1}(c)>0$, which is a direct
consequence of Lemma 4.1. Now the number of variations of sign in (4.7) at
$t=-M$ is $m$ for sufficiently large $M$ since $\theta_{m}>0$; it is zero at
$t=-\delta$ for sufficiently small $\delta>0$ since $\theta_{0}=1$. Thus by
Sturm theorem we conclude that $f_{m}(t)$ has exactly $m$ zeros in
$[-M,-\delta]$. This completes the proof. ∎

Based on this theorem, we can use the following parallel-in-time algorithm to
compute the solution $\mathbf{Y}_{r}$ to the problem (2.2).

######  Algorithm 4.1.

Given $\mathbf{Y}_{r}(t_{0})=\mathbf{Y}_{0}$.
$1^{\circ}$ Call Algorithm 3.1 to obtain $\mathbf{Y}_{r}(t_{n})$ for
$n=1,\ldots,N$.
$2^{\circ}$ Compute $\mathbf{Y}_{r}$ in each time interval $I_{n}$,
$n=1,\cdots,N-1$, in parallel by solving (4.1) to obtain
$\mathbf{a}_{1},\cdots,\mathbf{a}_{r}$ using the LU decomposition for block
tridiagonal matrices.

We remark that the algorithm of $LU$ decomposition for block tridiagonal
matrices for solving (4.1) requires to solve $r$ systems of linear equations
of size $M$ in sequential instead of to solve $r^{2}$ systems of linear
equations of size $M$ in parallel in Step $2^{\circ}$ of Algorithm 3.2.

##  5 Optimal stability and error esstimates

In this section, we show optimal stability and error estimates of the
continuous time Galerkin method (2.2) in terms of $r$ when $\mathbb{D}$ is a
symmetric or skew-symmetric matrix. This will be achieved by using the
explicit formulas in Theorem 3.1. We start by studying further properties of
the minors of the stiffness matrix of the continuous time Galerkin method
$\mathbb{A}=\mathbb{E}_{r+1}(\tau_{n}\mathbb{D})$.

Let $\chi_{r+1,j}(\lambda)=(-1)^{r+1}\det[\mathbb{E}_{r+1}(\lambda)_{r+1,j}]$,
then by (3.38) we have

$\chi_{2,1}=-1,\chi_{2,2}=a_{1},\chi_{3,1}=-1,\chi_{3,2}=a_{1},\chi_{3,3}=-a_{1}a_{2}.$
(5.1)

For $r\geq 3$, by $4^{\circ}$ in Theorem 3.2, we have the following recursive
formulas

$\displaystyle\chi_{r+1,j}(\lambda)=\chi_{r,j}(\lambda)+a_{r}a_{r-1}\chi_{r-1,j}(\lambda),\
\ 1\leq j\leq r-2,$ (5.2)
$\displaystyle\chi_{r+1,j}(\lambda)=(-1)^{j}\prod^{j-1}_{k=1}a_{k},\ \
j=r-1,r,r+1.$ (5.3)

Let $\varphi_{0}=1$ and $\varphi_{r}(\lambda)=\det\mathbb{E}_{r+1}(\lambda)$,
$r\geq 1$. Then by Lemma 3.2, $\varphi_{1}=1-a_{1}$,

$\displaystyle\varphi_{r+1}=\varphi_{r}+a_{r}a_{r+1}\varphi_{r-1},\ \ r\geq
1.$ (5.4)

######  Lemma 5.1.

For any $r\geq 1$ and $\lambda\in\mathbb{R}$, we have

$\displaystyle\sum_{j=1}^{r-1}[\chi_{r,j}(-\lambda)\chi_{r+1,j}(\lambda)+\chi_{r,j}(\lambda)\chi_{r+1,j}(-\lambda)]\frac{1}{2j-1}$
$\displaystyle=\varphi_{r-1}(-\lambda)\varphi_{r}(\lambda)+\varphi_{r-1}(\lambda)\varphi_{r}(-\lambda)+\frac{(-1)^{r}2r}{2r-1}(a_{1}\cdots
a_{r-1})^{2},$ (5.5)
$\displaystyle\sum_{j=1}^{r+1}\chi_{r+1,j}(-\lambda)\chi_{r+1,j}(\lambda)\frac{1}{2j-1}=\varphi_{r}(-\lambda)\varphi_{r}(\lambda)+\frac{(-1)^{r+1}2r}{2r+1}(a_{1}\cdots
a_{r})^{2}.$ (5.6)

###### Proof.

We denote

$\displaystyle A_{r}:$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{r-1}\left[\chi_{r,j}(-\lambda)\chi_{r+1,j}(\lambda)+\chi_{r,j}(\lambda)\chi_{r+1,j}(-\lambda)\right]\frac{1}{2j-1}$
$\displaystyle-\left[\varphi_{r-1}(-\lambda)\varphi_{r}(\lambda)+\varphi_{r-1}(\lambda)\varphi_{r}(-\lambda)\right],$
$\displaystyle B_{r+1}:$ $\displaystyle=$
$\displaystyle\sum_{j=1}^{r+1}\chi_{r+1,j}(-\lambda)\chi_{r+1,j}(\lambda)\frac{1}{2j-1}-\varphi_{r}(-\lambda)\varphi_{r}(\lambda).$

We will argue by induction. First (5.5)-(5.6) are obvious for $r=1,2$ by
(5.1). Now we assume (5.5)-(5.6) are valid for all $r\leq n$, $n\geq 2$. Since
by (5.3), $\chi_{n+2,n}(\lambda)=\chi_{n+1,n}(\lambda)$, we have by (5.2) and
(5.4) that

$\displaystyle
A_{n+1}=2B_{n+1}-2\chi_{n+1,n+1}(-\lambda)\chi_{n+1,n+1}(\lambda)\frac{1}{2n+1}+a_{n}a_{n+1}A_{n}.$

Now by (5.3) and the induction assumption that (5.5)-(5.6) are valid for
$r=n$, we obtain

$\displaystyle A_{n+1}=\frac{(-1)^{n+1}2(n+1)}{2n+1}(a_{1}\cdots a_{n})^{2}.$

This shows (5.5) for $r=n+1$. Similarly, we can prove by (5.2)-(5.4) that

$\displaystyle B_{n+2}$ $\displaystyle=$ $\displaystyle
B_{n+1}+(a_{n}a_{n+1})^{2}B_{n}+a_{n}a_{n+1}A_{n}$
$\displaystyle+\chi_{n+2,n+2}(-\lambda)\chi_{n+2,n+2}(\lambda)\frac{1}{2n+3}-(a_{n}a_{n+1})^{2}\chi_{n,n}(-\lambda)\chi_{n,n}(\lambda)\frac{1}{2n-1}.$

Now by the induction assumption (5.6) for $r=n,n+1$ and (5.5) for $r=n$, we
obtain by using (5.3) that

$\displaystyle B_{n+2}=(-1)^{n+2}\frac{2(n+1)}{2n+3}(a_{1}\cdots
a_{n+1})^{2}.$

This completes the proof. ∎

######  Lemma 5.2.

Let $r\geq 1$. For any $\lambda\leq 0$, we have

$\sum^{r-1}_{j=1}\chi_{r,j}(\lambda)\chi_{r+1,j}(\lambda)\frac{1}{2j-1}\leq\varphi_{r-1}(\lambda)\varphi_{r}(\lambda),\
\
\sum^{r+1}_{j=1}\chi_{r+1,j}(\lambda)^{2}\frac{1}{2j-1}\leq\varphi_{r}(\lambda)^{2}.$
(5.7)

###### Proof.

For any $r\geq 1$, we denote

$\displaystyle
C_{r}:=\sum^{r-1}_{j=1}\chi_{r,j}\chi_{r+1,j}\frac{1}{2j-1}-\varphi_{r-1}\varphi_{r},\
\ D_{r+1}:=\sum^{r+1}_{j=1}\chi^{2}_{r+1,j}\frac{1}{2j-1}-\varphi^{2}_{r}.$

We again argue by induction. First (5.7) is obvious for $r=1,2$ by (5.1) since
$\lambda\leq 0$. Now we assume (5.7) is valid for all $r\leq n$, $n\geq 2$. By
(5.2)-(5.4), it is easy to see that

$\displaystyle
C_{n+1}=-\chi^{2}_{n+1,n+1}\frac{1}{2n+1}+D_{n+1}+a_{n}a_{n+1}C_{n},$

where we have used $\chi_{n+2,n}(\lambda)=\chi_{n+1,n}(\lambda)$. Thus if
$C_{n}\leq 0,D_{n+1}\leq 0$, then $C_{n+1}\leq 0$.

On the other hand, by (5.2)-(5.4), we have

$\displaystyle
D_{n+2}=D_{n+1}+a_{n}a_{n+1}D_{n}+2a_{n}a_{n+1}C_{n}+(a_{1}\cdots
a_{n+1})^{2}(\frac{1}{2n+3}-\frac{1}{2n-1}).$

Thus $D_{n+2}\leq 0$ if $D_{n}\leq 0,D_{n+1}\leq 0$, and $C_{n}\leq 0$. This
completes the proof. ∎

The following theorem is the main result of this section.

######  Theorem 5.1.

Let $\mathbb{D}$ be a symmetric or skew-symmetric matrix and
$\mathbf{Y}_{r}\in\mathbf{V}^{r}_{\tau}$ is the solution of the problem (2.2).
Then we have

$\displaystyle\|\mathbf{Y}_{r}\|_{L^{2}(0,T)}\leq
T^{1/2}\|\mathbf{Y}_{0}\|_{\mathbb{R}^{M}}+CT\|\mathbf{R}\|_{L^{2}(0,T)},$
(5.8) $\displaystyle\max_{0\leq t\leq
T}\|\mathbf{Y}_{r}\|_{\mathbb{R}^{M}}\leq
Cr(\|\mathbf{Y}_{0}\|_{\mathbb{R}^{M}}+T^{1/2}\|\mathbf{R}\|_{L^{2}(0,T)}),$
(5.9)

where the constant $C$ is independent of $\tau,r,\mathbb{D}$, and
$\mathbf{R}$.

###### Proof.

Let $\hat{\mathbf{Y}}_{r}\in[P^{r}]^{M}$ be defined in (2.6) of Lemma 2.1, we
claim that

$\displaystyle\|\hat{\mathbf{Y}}_{r}\|_{L^{2}(I_{n})}\leq\tau_{n}^{1/2}\|\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}},$
(5.10)

which improves the bound (2.10) in the proof of Lemma 2.1. To show (2.6), by
Theorem 3.1, we have

$\displaystyle\hat{\mathbf{Y}}_{r}(t)$
$\displaystyle=\sum_{j=1}^{r+1}\mathbf{a}_{j-1}\tilde{L}_{j-1}(t)=\sum_{j=1}^{r+1}(-1)^{j}\frac{\chi_{r+1,j}(\tau_{n}\mathbb{D})}{\varphi_{r}(\tau_{n}\mathbb{D})}\mathbf{Y}_{r}^{n}\tilde{L}_{j-1}(t)\
\ \forall t\in I_{n}.$

Thus by (2.5),

$\displaystyle\|\hat{\mathbf{Y}}_{r}\|_{L^{2}(I_{n})}^{2}=\sum_{j=1}^{r+1}\|\chi_{r+1,j}(\tau\mathbb{D})\varphi_{r}(\tau_{n}\mathbb{D})^{-1}\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}}^{2}\frac{\tau_{n}}{2j-1}.$
(5.11)

Denote
$\mathbf{Z}_{r}^{n}=\varphi_{r}(\tau_{n}\mathbb{D})^{-1}\mathbf{Y}_{r}^{n}$.
If $\mathbb{D}$ is skew-symmetric $\mathbb{D}^{T}=-\mathbb{D}^{T}$, by (2.5),
(5.6), we have

$\displaystyle\|\hat{\mathbf{Y}}_{r}\|_{L^{2}(I_{n})}^{2}$
$\displaystyle=\sum_{j=1}^{r+1}\|\chi_{r+1,j}(\tau_{n}\mathbb{D})\mathbf{Z}_{r}^{n}\|_{\mathbb{R}^{M}}^{2}\frac{\tau_{n}}{2j-1}$
$\displaystyle=\tau_{n}\|\varphi_{r}(\tau_{n}\mathbb{D})\mathbf{Z}_{r}^{n}\|_{\mathbb{R}^{M}}^{2}-\tau_{n}\left\|\sqrt{\frac{2r}{2r+1}}\frac{r!}{(2r)!}(\tau_{n}\mathbb{D})^{r}\mathbf{Z}_{r}^{n}\right\|_{\mathbb{R}^{M}}^{2}$
$\displaystyle\leq\tau_{n}\|\mathbf{Y}_{r}^{n}\|_{\mathbb{R}^{M}}^{2}.$

This shows the claim (5.10) when $\mathbb{D}$ is skew-symmetric.

If $\mathbb{D}$ is symmetric, the eigenvalues of $\mathbb{D}$ are non-positive
since $\mathbb{D}+\mathbb{D}^{T}\leq 0$. By Lemma 5.2, it is easy to show that

$\displaystyle\sum^{r+1}\|\chi_{r+1,j}(\tau_{n}\mathbb{D})\mathbf{Z}^{n}_{r}\|_{\mathbb{R}^{M}}^{2}\frac{\tau_{n}}{2j-1}-\tau_{n}\|\varphi_{r}(\tau_{n}\mathbb{D})\mathbf{Z}^{n}_{r}\|_{\mathbb{R}^{M}}^{2}\leq
0,$

which yields

$\displaystyle\sum^{r+1}_{j=1}\|\chi_{r+1,j}(\tau_{n}\mathbb{D})\varphi_{r}(\tau_{n}\mathbb{D})^{-1}\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}}^{2}\frac{\tau_{n}}{2j-1}\leq\tau_{n}\|\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}}^{2}.$

Now it follows from (5.11) that
$\|\hat{\mathbf{Y}}_{r}\|_{L^{2}(I_{n})}\leq\tau_{n}^{1/2}\|\mathbf{Y}^{n}_{r}\|_{\mathbb{R}^{M}}$.
This shows the claim (5.10) when $\mathbb{D}$ is a symmetric matrix.

It follows from (5.10), (2.12) and (2.3) that

$\displaystyle\|{\mathbf{Y}}_{r}\|_{L^{2}(I_{n})}$
$\displaystyle\leq\tau_{n}^{1/2}\|\mathbf{Y}_{r}^{n}\|_{\mathbb{R}^{M}}+2\tau_{n}\|\mathbf{R}\|_{L^{2}(I_{n})}$
$\displaystyle\leq\tau_{n}^{1/2}\|\mathbf{Y}_{0}\|_{\mathbb{R}^{M}}+C\tau_{n}^{1/2}T^{1/2}\|\mathbf{R}\|_{L^{2}(0,T)}$

This implies (5.8) easily. Now by the $hp$ inverse estimate,

$\displaystyle\max_{t_{n}\leq t\leq
t_{n+1}}\|\mathbf{Y}_{r}\|_{\mathbb{R}^{M}}\leq
C\tau_{n}^{-1/2}r\|\mathbf{Y}_{r}\|_{L^{2}(I_{n})}\leq
Cr(\|\mathbf{Y}_{0}\|_{\mathbb{R}^{M}}+T^{1/2}\|\mathbf{R}\|_{L^{2}(0,T)}).$

This shows (5.9). This completes the proof. ∎

The following theorem which improves the error estimates in Theorem 2.1 can be
proved by the argument in Theorem 2.1 by using Theorem 5.1 instead of Lemma
2.1. Here we omit the details.

######  Theorem 5.2.

Let $\mathbb{D}$ be a symmetric or skew-symmetric matrix. Assume that
$\mathbf{R}\in[H^{s}(0,T)]^{M}$, $\mathbf{Y}\in[W^{1+s,\infty}(0,T)]^{M}$,
$s\geq 1$, and $\mathbf{Y}_{r}\in\mathbf{V}^{r}_{\tau}$ is the solution of the
problem (2.2), we have

$\displaystyle\|\mathbf{Y}-\mathbf{Y}_{r}\|_{L^{2}(0,T)}\leq
C(1+T)\frac{\tau^{\min(r+1,s)}}{r^{s}}(\|\mathbf{Y}\|_{H^{s}(0,T)}+\|\mathbb{D}\mathbf{Y}\|_{H^{s}(0,T)}),$
$\displaystyle\max_{0\leq t\leq
T}\|\mathbf{Y}-\mathbf{Y}_{r}\|_{\mathbb{R}^{M}}\leq
C(1+T^{1/2})\frac{\tau^{\min(r+1,s)}}{r^{s-1}}(T^{1/2}\|\mathbf{Y}\|_{W^{s+1,\infty}(0,T)}+\|\mathbf{R}\|_{H^{s}(0,T)}),$

where the constant $C$ is independent of $\tau,r,\mathbb{D}$, and $\mathbf{R}$
but may depend on $s$.

We remark that the first estimate in Theorem 5.2 is optimal both in $\tau$ and
$r$.

##  6 Numerical examples

In this section, we provide some numerical examples to confirm the theoretical
results in this paper.

######  Example 1.

(Dissipative problem) Let $\Omega=(0,1)\times(0,1)$ and $T=1$. We consider the
following constant coefficient convection-diffusion problem

$\displaystyle\left\\{\begin{aligned}
&u_{t}+\nabla\cdot(\bm{\beta}u-\epsilon\nabla u)=f\ \ \mbox{in
}\Omega\times(0,T),\\\ &u(\mathbf{x},0)=u_{0}(\mathbf{x})\ \ \mbox{in
}\Omega.\end{aligned}\right.$ (6.1)

The boundary condition is set to be periodic. The source term $f$ is chosen
such that the exact solution is
$u(\mathbf{x},t)=\exp(-t)\sin(4\pi(x_{1}-t))\cos(4\pi(x_{2}-t))$.

We choose $\bm{\beta}=(1,1)^{T}$ and $\epsilon=1$ in (6.1). For spatial
discretizations, we apply the local discontinuous Galerkin (LDG) method in
Cockburn and Shu [9] by using purely upwind fluxes for convection terms and
alternating fluxes for diffusion terms. For the sake of completeness, we
recall the method for solving (6.1) here.

Let $\mathcal{M}$ denote a uniform Cartesian mesh of $\Omega$ with $h$ the
length of the sides of the elements. $\mathcal{E}=\mathcal{E}^{\rm
side}\cup\mathcal{E}^{\rm bdy}$, where $\mathcal{E}^{\rm side}:=\\{e=\partial
K\cap\partial K^{\prime}:K,K^{\prime}\in\mathcal{M}\\}$, $\mathcal{E}^{\rm
bdy}:=\\{e=\partial K\cap\partial\Omega:K\in\mathcal{M}\\}$. For any subset
$\widehat{\mathcal{M}}\subset\mathcal{M}$ and
$\widehat{\mathcal{E}}\subset\mathcal{E}$, we use the notation

$\displaystyle(u,v)_{\widehat{\mathcal{M}}}=\sum_{K\in\widehat{\mathcal{M}}}(u,v)_{K},\
\ \langle
u,v\rangle_{\widehat{\mathcal{E}}}=\sum_{e\in\widehat{\mathcal{E}}}\langle
u,v\rangle_{e},$

where $(\cdot,\cdot)_{K}$ and $\langle\cdot,\cdot\rangle_{e}$ denote the inner
product of $L^{2}(K)$ and $L^{2}(e)$, respectively.

For any $e\in\mathcal{E}$, we fix a unit normal vector $\mathbf{n}_{e}$ of $e$
with the convention that $\mathbf{n}_{e}$ is the unit outer normal to
$\partial\Omega$ if $e\in\mathcal{E}^{\rm bdy}$. For any $v\in
H^{1}(\mathcal{M}):=\\{v:v\in H^{1}(K),K\in\mathcal{M}\\}$, we define the jump
operator of $v$ across $e$:

$\displaystyle[\\![v]\\!]_{e}:=v^{-}-v^{+}\ \ \forall e\in\mathcal{E}^{\rm
side},\ \ \ \ [\\![v]\\!]_{e}:=v^{-}\ \ \forall e\in\mathcal{E}^{\rm bdy},$

where $v^{\pm}(\mathbf{x}):=\lim_{\varepsilon\rightarrow
0^{+}}v(\mathbf{x}\pm\varepsilon\mathbf{n}_{e})\ \forall\mathbf{x}\in e$. For
any integer $p\geq 0$, we define the finite element space

$\displaystyle V_{h}^{p}:=\\{v\in L^{2}(\Omega):v|_{K}\in
Q^{p}(K),K\in\mathcal{M}\\},$

where $Q^{p}(K)$ denotes the space of polynomials of degree at most $p$ in
each variable in $K$.

The semi-discrete problem is to find
$(u_{h},\mathbf{q}_{h})\in[V_{h}^{p}]^{3}$ such that, for all test functions
$(v_{h},\mathbf{r}_{h})\in[V_{h}^{p}]^{3}$,

$\displaystyle(\partial_{t}u_{h},v_{h})_{\mathcal{M}}+\mathcal{G}(\bm{\beta}u_{h},v_{h})=\sqrt{\epsilon}\,\left[-(\mathbf{q}_{h},\nabla
v_{h})_{\mathcal{M}}+\langle{\mathbf{q}}^{-}_{h}\cdot\mathbf{n},[\\![v_{h}]\\!]\rangle_{\mathcal{E}}\right]+(f,v_{h})_{\mathcal{M}},$
$\displaystyle(\mathbf{q}_{h},\mathbf{r}_{h})_{\mathcal{M}}=\sqrt{\epsilon}\,\left[-(u_{h},\text{div}\mathbf{r}_{h})_{\mathcal{M}}+\langle{u}^{+}_{h},[\\![\mathbf{r}_{h}]\\!]\cdot\mathbf{n}\rangle_{\mathcal{E}}\right],$
$\displaystyle u_{h}(\mathbf{x},0)=(\mathcal{P}_{h}u_{0})(\mathbf{x})\ \
\mbox{in }\Omega.$

Here $\mathcal{P}_{h}:L^{2}(\Omega)\to V_{h}^{p}$ is the standard $L^{2}$
projection operator, and

$\displaystyle\mathcal{G}(\bm{\beta}u_{h},v_{h})=-(\bm{\beta}u_{h},\nabla
v_{h})_{\mathcal{M}}+\langle\check{u}_{h}\bm{\beta}\cdot\mathbf{n},[\\![v_{h}]\\!]\rangle_{\mathcal{E}},$

where $\check{u}_{h}$ is chosen as the upwind flux: $\check{u}_{h}=u_{h}^{-}$
if $\bm{\beta}\cdot\mathbf{n}>0$, $\check{u}=u^{+}_{h}$ if
$\bm{\beta}\cdot\mathbf{n}<0$. For $e\in\mathcal{E}^{\rm bdy}$, we use the
periodic boundary condition to define ${u}_{h}^{+}$.

The optimal $L^{2}$-norm error estimate of order $p+1$ of the semi-discrete
scheme for quasi-uniform Cartesian meshes can be found in Cheng et al [8,
Theorem 2.4], where it is shown that $\max_{0\leq t\leq
T}\|u-u_{h}\|_{L^{2}(\Omega)}\leq C(1+T)h^{p+1}$. Therefore, combined with the
continuous time Galerkin scheme, we know that the fully discrete scheme has
$O(h^{p+1}+\tau^{r+1})$ accuracy in the norm
$\|\cdot\|_{L^{\infty}(0,T;L^{2}(\Omega))}$ and $O(h^{p+1}+\tau^{2r})$ in the
$L^{2}$ norm at nodes $t=t_{n}$, $n=1,\cdots,N$.

To test the accuracy at the nodes, we set $\tau=h^{\frac{p+1}{2r}}$ and thus
$N=T/\tau=T\beta^{\frac{1}{r}}$, where $\beta=h^{-\frac{p+1}{2}}$. The
numerator of the $[r/r]$ Padé approximation $P_{r}(z)$ has $2k$ complex zeros
and $1$ real root if $r=2k+1,k\geq 1$, and $2k$ complex zeros if $r=2k$,
$k\geq 1$. Denote by $C(2M)$ the costs of solving the matrix problem
$\tau_{n}\mathbb{D}+\zeta_{j}\mathbb{I}$ with $\zeta_{j}$ being complex and
$C(M)$ the costs of solving the matrix problem
$\tau_{n}\mathbb{D}+\zeta_{j}\mathbb{I}$ with $\zeta_{j}$ real, where
$\zeta_{j}$, $j=1,\cdots,r$, are zeros of $P_{r}(z)$. Then the computational
time in each time step of Algorithm 3.1 is proportion to $C(2M)$ for the
parallel computation and proportion to $kC(2M)+(r-2k)C(M)$ for the sequential
computation. The wall time of using Algorithm 3.1 using parallel machines is
then proportion to $N=T\beta^{\frac{1}{r}}$ which is decreasing in $r$. Thus
high order time discretization is preferred for parallel computations. On the
other hand, for the sequential computation, the wall time of using Algorithm
3.1 is proportion to $rN=Tr\beta^{\frac{1}{r}}$ which minimizes at
$r=\ln\beta$ for $r>0$. This implies that the optimal choice of the order for
the sequential computation is $r=\lfloor\ln\beta\rfloor+1$, where $\lfloor
a\rfloor$ is the maximum integer strictly less than $a>0$. Table 1 shows the
error $\|(u-u_{h})(\cdot,T)\|_{L^{2}(\Omega)}$ at the terminal time when
$r=\lfloor\ln\beta\rfloor+1$. The optimal $(p+1)$-th order is observed which
confirms our theoretical results. We observe that the errors of high order
schemes are significant smaller than the low order schemes.

To test the accuracy in the $\|\cdot\|_{L^{\infty}(0,T;L^{2}(\Omega))}$ norm,
we set $\tau=h^{\frac{p+1}{r+1}}$ and thus $N=T/\tau=T\gamma^{\frac{1}{r+1}}$,
where $\gamma=h^{-(p+1)}$. The wall time of using Algorithm 3.1 and Algorithm
3.2 for the parallel computation is proportion to $N=T\gamma^{\frac{1}{r+1}}$
which decreases in $r$. On the other hand, for the sequential computation, the
wall time of using Algorithm 3.1 and Algorithm 4.1 is
$rN=Tr\gamma^{\frac{1}{r+1}}$ which is increasing in $r$ if $\ln\gamma\leq 4$
and minimizes at $r^{*}=[-(2-\ln\gamma)+\sqrt{(2-\ln\gamma)^{2}-4}\,]/2$ if
$\ln\gamma\geq 4$. Since $r^{*}\geq 1$ is equivalent to $\ln\gamma\geq 4$, the
optimal choice of the order for minimizing the computation wall time is
$r=\max(1,\lfloor r^{*}\rfloor+1)$. We note that for the sequential
computation, Algorithm 4.1 is cheaper than Algorithm 3.2. Table 2 shows the
error

$\displaystyle\max_{0\leq n\leq N-1,1\leq k\leq
10}\|(u-u_{h})(\cdot,t_{n}+0.1k\tau_{n})\|_{L^{2}(\Omega)}$

as the approximation of $\|u-u_{h}\|_{L^{\infty}(0,T;L^{2}(\Omega))}$ when
$r=\max(1,\lfloor r^{*}\rfloor+1)$. We again observe the optimal $(p+1)$-th
order convergence and that high order methods perform much better than low
order methods.

Table 1: Example 1: numerical errors of $\|(u-u_{h})(\cdot,T)\|_{L^{2}(\Omega)}$ and orders. | $p=3$ | $p=4$ | $p=5$
---|---|---|---
$h$ | error | order | error | order | error | order
$1/4$ | 8.06E-03 | – | 1.29E-03 | – | 1.71E-04 | –
$1/8$ | 5.64E-04 | 3.84 | 4.40E-05 | 4.88 | 2.84E-06 | 5.91
$1/16$ | 3.55E-05 | 3.99 | 1.37E-06 | 5.01 | 4.33E-08 | 6.04
$1/32$ | 2.03E-06 | 4.13 | 4.26E-08 | 5.00 | 6.94E-10 | 5.96
Table 2: Example 1: numerical errors in $\|\cdot\|_{L^{\infty}(0,T;L^{2}(\Omega))}$ norm and orders. | $p=3$ | $p=4$ | $p=5$
---|---|---|---
$h$ | error | order | error | order | error | order
$1/4$ | 2.77E-02 | – | 3.99E-03 | – | 5.81E-04 | –
$1/8$ | 2.00E-03 | 3.79 | 1.36E-04 | 4.87 | 1.02E-05 | 5.83
$1/16$ | 1.14E-04 | 4.14 | 4.40E-06 | 4.95 | 1.69E-07 | 5.91
$1/32$ | 7.18E-06 | 3.98 | 1.33E-07 | 5.05 | 2.73E-09 | 5.95

######  Example 2.

(Wave propagation problem) Let $\Omega=(-2,2)\times(-2,2)$ and $T=1$. We
consider the following wave equation with discontinuous coefficients

$\displaystyle\left\\{\begin{aligned} &\frac{1}{\rho
c^{2}}\partial_{t}u=\mbox{\rm div}\,\mathbf{q}+f,\ \
\rho\partial_{t}\mathbf{q}=\nabla u\ \ \ \ \mbox{in }\Omega\times(0,T),\\\
&[\\![u]\\!]=0,\quad[\\![\mathbf{q}\cdot\mathbf{n}]\\!]=0\ \ \ \ \text{ on
}\Gamma\times(0,T),\\\ &u=0\quad\text{ on }\partial\Omega\times(0,T),\\\
&u(\mathbf{x},0)=u_{0}(\mathbf{x}),\ \
\mathbf{q}(\mathbf{x},0)=\mathbf{q}_{0}(\mathbf{x})\ \ \ \ \mbox{in
}\Omega.\end{aligned}\right.$ (6.2)

We assume the interface $\Gamma$ is the union of two closely located ellipses.
We take
$\Omega_{1}=\\{\mathbf{x}\in\Omega:\frac{(x_{1}-d_{1})^{2}}{a^{2}}+\frac{x_{2}^{2}}{b^{2}}<1\text{
or }\frac{(x_{1}-d_{2})^{2}}{a^{2}}+\frac{x_{2}^{2}}{b^{2}}<1\\}$, which is
the union of two disks, and $\Omega_{2}=\Omega\backslash\bar{\Omega}_{1}$.
Here $d_{1}=-0.82$, $d_{2}=0.82$, $a=0.81$, and $b=0.51$. The distance between
two ellipses is $0.02$. We consider the wave equation (6.2) with
$\rho_{1}=1/2,\,\rho_{2}=1$, $c_{1}=c_{2}=1$, and the source $f$ is chosen
such that the exact solution is

$\displaystyle
u(\mathbf{x},t)=\left\\{\begin{array}[]{cc}\cos(3t)\sin(r_{1}-1)\sin(r_{2}-1)\sin(3\pi
x_{1})\sin(3\pi x_{2})&\text{ in }\Omega_{1},\\\ \\\
2\cos(3t)\sin(r_{1}-1)\sin(r_{2}-1)\sin(3\pi x_{1})\sin(3\pi x_{2})&\text{ in
}\Omega_{2},\end{array}\right.$

where $r_{1}=\frac{(x_{1}-d_{1})^{2}}{a^{2}}+\frac{x_{2}^{2}}{b^{2}}\quad
r_{2}=\frac{(x_{1}-d_{2})^{2}}{a^{2}}+\frac{x_{2}^{2}}{b^{2}}$. The exact
solution $\mathbf{q}(\mathbf{x},t)$ is computed by (6.2) with the initial
condition $\mathbf{q}_{0}=0$.

We use the unfitted finite element method in Chen et al [7] to discretize the
problem in space. Let $\mathcal{M}$ be an induced mesh which is constructed
from a Cartesian partition $\mathcal{T}$ of the domain $\Omega$ with possible
local refinements and hanging nodes so that the elements are large with
respect to both domains $\Omega_{1},\Omega_{2}$. Let
$\mathcal{M}^{\Gamma}:=\\{K\in\mathcal{M}:K\cap\Gamma\not=\emptyset\\}$ and
$\mathcal{E}=\mathcal{E}^{\rm
side}\cup\mathcal{E}^{\Gamma}\cup\mathcal{E}^{\rm bdy}$, where
$\mathcal{E}^{\Gamma}:=\\{\Gamma_{\\!\\!K}=\Gamma\cap K:K\in\mathcal{M}\\}$.

For any $K\in\mathcal{M}^{\Gamma}$, $i=1,2$, let $K_{i}=K\cap\Omega_{i}$ and
$K_{i}^{h}$ the polygonal approximation of $K_{i}$ bounded by the sides of $K$
and $\Gamma_{K}^{h}$ which is the line segment connecting two intersection
points of $\Gamma_{K}\cap\partial K$. $K_{i}^{h}$ is the union of shape
regular triangles $K_{ij}^{h}$, $1\leq J_{i}^{K}\leq 3$, whose sides are the
sides of $K_{i}^{h}$ and $\Gamma_{K}^{h}$. We always set $K_{i1}^{h}$ the
element having $\Gamma_{K}^{h}$ as one of its sides. From $K_{ij}^{h}$ we
define the curved element $\widetilde{K}_{ij}^{h}$ by

$\displaystyle\widetilde{K}_{i1}^{h}=(K_{i}\cap
K_{i1}^{h})\cup(K_{i}\backslash\bar{K}_{i1}^{h}),\ \
\widetilde{K}_{ij}^{h}=K_{i}\cap K_{ij}^{h},\ \ j=2,\cdots,J_{i}^{K}.$

Then we know that $K$ is the union of curved triangles
$\widetilde{K}_{ij}^{h}$, $i=1,2,j=1,\cdots,J_{i}^{K}$.

For any integers $p,q\geq 1$, the space $P^{p}(K)$ denotes the space of
polynomials of degree at most $p$ in $K$ and $Q^{p,q}(K)$ denotes the space of
polynomials of degree at most $p$ for the first variable and $q$ for the
second variable in $K$. For any $K\in\mathcal{M}^{\Gamma}$, we define the
interface finite element spaces

$\displaystyle W_{p}(K)=\\{\varphi:\varphi|_{\widetilde{K}_{ij}^{h}}\in
P^{p}(\widetilde{K}_{ij}^{h}),\ i=1,2,\ j=1,\cdots,J_{i}^{K}\\},$

and $X_{p}(K)=W_{p}(K)\cap H^{1}(K_{1}\cup K_{2})$. Notice that the functions
in $X_{p}(K)$ are conforming in each $K_{i},i=1,2$. Now we define the
following unfitted finite element spaces

$\displaystyle X_{p}(\mathcal{M}):=\\{v\in
H^{1}(\Omega_{1}\cup\Omega_{2}):v|_{K}\in X_{p}(K)\ \ \forall
K\in\mathcal{M}^{\Gamma},$ $\displaystyle\hskip 136.5733ptv|_{K}\in Q^{p}(K)\
\ \forall K\in\mathcal{M}\backslash\mathcal{M}^{\Gamma}\\},$
$\displaystyle\mathbf{W}_{p}(\mathcal{M}):=\\{\bm{\psi}:\bm{\psi}|_{K}\in[W_{p}(K)]^{2}\
\ \forall K\in\mathcal{M}^{\Gamma},$ $\displaystyle\hskip
122.34692pt\bm{\psi}|_{K}\in Q^{p-1,p}(K)\times Q^{p,p-1}(K)\ \ \forall
K\in\mathcal{M}\backslash\mathcal{M}^{\Gamma}\\}.$

Let $X_{p}^{0}(\mathcal{M})=X_{p}(\mathcal{M})\cap
H^{1}_{0}(\Omega_{1}\cup\Omega_{2})$, where
$H^{1}_{0}(\Omega_{1}\cup\Omega_{2})=\\{v\in
H^{1}(\Omega_{1}\cup\Omega_{2}):v=0\mbox{ on }\partial\Omega\\}$.

The semi-discrete unfitted finite element method for solving (6.2) is then to
find $(u_{h},\mathbf{q}_{h})\in
X_{p}^{0}(\mathcal{M})\times\mathbf{W}_{p}(\mathcal{M})$ such that for all
$(\varphi_{h},\bm{\psi}_{h})\in
X_{p}^{0}(\mathcal{M})\times\mathbf{W}_{p}(\mathcal{M})$,

$\displaystyle\Big{(}\frac{1}{\rho
c^{2}}\partial_{t}u_{h},\varphi_{h}\Big{)}_{\mathcal{M}}=-(\mathbf{q}_{h},\nabla\varphi_{h})_{\mathcal{M}}+\langle{\mathbf{q}}_{h}^{-}\cdot\mathbf{n},[\\![\varphi_{h}]\\!]\rangle_{\mathcal{E}^{\Gamma}}+(f,\varphi_{h})_{\mathcal{M}},$
(6.3)
$\displaystyle(\rho\partial_{t}\mathbf{q}_{h},\bm{\psi}_{h})_{\mathcal{M}}=-(u_{h},\text{div}\bm{\psi}_{h})_{\mathcal{M}}+\langle{u}_{h}^{+},[{\hskip-1.5pt}[\bm{\psi}_{h}]{\hskip-1.5pt}]\cdot\mathbf{n}\rangle_{\mathcal{E}},$
(6.4) $\displaystyle u_{h}(\mathbf{x},0)=(\mathcal{P}_{h}u_{0})(\mathbf{x}),\
\ \mathbf{q}_{h}(\mathbf{x},0)=(\bm{P}_{h}\mathbf{q}_{0})(\mathbf{x})\ \
\mbox{in }\Omega,$ (6.5)

where $\mathcal{P}_{h}:L^{2}(\Omega)\to X_{p}^{0}(\mathcal{M})$ and
$\bm{P}_{h}:[L^{2}(\Omega)]^{2}\to\mathbf{W}_{p}(\mathcal{M})$ are the
standard $L^{2}$ projection operators.

It is shown in [7, Theorem 2.2] that the following energy error of the semi-
discrete scheme

$\displaystyle
E_{en}(t):=(\|(u-{u}_{h})(\cdot,t)\|_{L^{2}(\Omega)}^{2}+\|(\mathbf{q}-{\mathbf{q}}_{h})(\cdot,t)\|_{L^{2}(\Omega)}^{2})^{1/2}.$

has $p$-th order convergence. The semi-discrete problem (6.3)-(6.5) is an ODE
system which is solved by the continuous time Galerkin method in this paper.
By Theorem 2.1 and Theorem 2.2, we know that the energy error has
$O(h^{p}+\tau^{r+1})$ convergence rate in the norm $\max_{0\leq t\leq
T}E_{en}(t)$ and $O(h^{p}+\tau^{2r})$ in the $E_{en}(t)$ at nodes $t=t_{n}$,
$n=1,\ldots,N$.

Note that these two ellipses are close but not tangent. To resolve the
interface $\Gamma$ well, we locally refine the mesh near the interface such
that the interface deviation $\eta_{K}\leq\eta_{0}=0.05$ for all
$K\in\mathcal{M}^{\Gamma}$. For the concept of the interface deviation we
refer to [7, Definition 2.2], see also Chen et al [6, Definition 2.2]. As an
illustration, we show the computational mesh for $h=1/4$ in Figure 1.

In this example, we test the accuracy of the error in the
$\|\cdot\|_{L^{\infty}(0,T;L^{2}(\Omega))}$ norm and $E_{en}(T)$. As in
Example 1, to test the accuracy at nodes, we set $\tau=\frac{p}{2r}$, thus the
wall time of using Algorithm 3.1 for the sequential computation is proportion
to $Tr\nu^{\frac{1}{r}}$ which minimizes at $r=\ln\nu$ for $r>0$, where
$\nu=h^{-\frac{p}{2}}$. Table 3 shows the error $E_{en}(T)$ at the terminal
time when $r=\lfloor\ln\nu\rfloor+1$.

To test the accuracy in the $\|\cdot\|_{L^{\infty}(0,T;L^{2}(\Omega))}$ norm,
we set $\tau=h^{\frac{p}{r+1}}$ and thus $N=T/\tau=T\mu^{\frac{1}{r+1}}$,
where $\mu=h^{-p}$. As in Example 1, the wall time of using Algorithm 3.1 and
Algorithm 3.2 for the parallel computation is proportion to
$N=T\mu^{\frac{1}{r+1}}$ which decreases in $r$. However, since the wave
equation is not dissipative, we cannot use Algorithm 4.1 for the sequential
computation. In this case, the wall time of using Algorithm 3.1 and Algorithm
3.2 for the sequential computation is proportion to
$rN+r^{2}N=Tr(r+1)\mu^{\frac{1}{r+1}}$ which is increasing if $\ln\mu\leq 6$
and minimizes at $r^{**}=[-(3-\ln\mu)+\sqrt{(3-\ln\mu)^{2}-8}\,]/4$ if
$\ln\mu\geq 6$. Since $r^{**}\geq 1$ is equivalent to $\ln\mu\geq 6$, the
optimal choice of the order for minimizing the computation wall time is
$r=\max(1,\lfloor r^{**}\rfloor+1)$. Table 3 shows the error

$\displaystyle\max_{0\leq n\leq N-1,1\leq k\leq 10}E_{en}(t_{n}+0.1k\tau_{n})$

as the approximation of $\max_{0\leq t\leq T}E_{en}(t)$ when $r=\max(1,\lfloor
r^{**}\rfloor+1)$. We clearly observe the optimal $p$-th order convergence and
the superior performance of high order methods from Tables 3-4.

Figure 1:  Illustration of the computational domain and the mesh (left) and the corresponding zoomed local mesh (right) with $h=1/4$ in Example 2. Table 3: Example 2: numerical errors of $E_{en}(T)$ and orders. | $p=3$ | $p=4$ | $p=5$
---|---|---|---
$h$ | error | order | error | order | error | order
$1/4$ | 4.22E-01 | – | 2.21E-01 | – | 1.25E-01 | –
$1/8$ | 8.97E-02 | 2.24 | 2.64E-02 | 3.06 | 6.20E-03 | 4.33
$1/16$ | 1.32E-02 | 2.76 | 1.81E-03 | 3.87 | 3.82E-05 | 4.96
$1/32$ | 1.66E-03 | 2.99 | 1.13E-04 | 4.00 | 1.19E-06 | 5.00
Table 4: Example 2: numerical errors of $\max_{0\leq t\leq T}E_{en}(t)$ and orders. | $p=3$ | $p=4$ | $p=5$
---|---|---|---
$h$ | error | order | error | order | error | order
$1/4$ | 5.23E-01 | – | 2.95E-01 | – | 1.69E-01 | –
$1/8$ | 1.26E-01 | 2.05 | 3.70E-02 | 3.00 | 8.87E-03 | 4.25
$1/16$ | 1.90E-02 | 2.72 | 2.58E-03 | 3.84 | 2.93E-04 | 4.92
$1/32$ | 2.45E-03 | 2.96 | 1.66E-04 | 3.95 | 9.21E-06 | 4.99

## References

  * [1] G. Akrivis, C. Makridakis, and R. H. Nochetto,  Galerkin and Runge-Kutta methods: unified formulation, a posteriori error estimates and nodal superconvergence,  Numer. Math., 118: 429–456, 2011.
  * [2] A.K. Aziz and P. Monk,  Continuous finite elements in space and time for the heat equation,  Math. Comp., 52: 255-274, 1989.
  * [3] C. Bernardi and Y. Maday,  Spectral methods,  Handbook of Numerical Analysis, 5: 209–485, 1997.
  * [4] W.C. Brown,  Matrices over Commutative Rings, Marcel Dekker, New York, 1993.
  * [5] J. Chen, Z. Chen, T. Cui, and L. Zhang,  An adaptive finite element method for the eddy current model with circuit/field couplings,  SIAM J. Sci. Comput., 32: 1020–1042, 2010.
  * [6] Z. Chen, K. Li, and X. Xiang,  An adaptive high-order unfitted finite element method for elliptic interface problems, Numer. Math., 149: 507-548, 2021.
  * [7] Z. Chen, Y. Liu, and X. Xiang,  A high order explicit time finite element method for the acoustic wave equation with discontinuous coefficients,  arXiv:2112.02867v2.
  * [8] Y. Cheng, X. Meng, and Q. Zhang,  Application of generalized Gauss-Radau projections for the local discontinuous Galerkin method for linear convection-diffusion equations,  Math. Comp., 86: 1233–1267, 2017.
  * [9] B. Cockburn and C.-W. Shu,  The local discontinuous Galerkin method for time-dependent convection-diffusion systems,  SIAM J. Numer. Anal., 35: 2440–2463, 1998.
  * [10] J. Douglas Jr., T. Dupont, and M.F. Wheeler,  A quasi projection analysis of Galerkin methods for parabolic and hyperbolic equations,  Math. Comp., 32: 345-362, 1978.
  * [11] D. A. French and T. E. Peterson,  A continuous space-time finite element method for the wave equation,  Math. Comp., 65: 491–506, 1996.
  * [12] E. Gallopoulos, and Y. Saad,  On the parallel solution of parabolic equations,  Proceedings of the 3rd international conference on Supercomputing, 17–28, 1989.
  * [13] M.J. Gander,  50 years of time parallel time integration,  in Multiple shooting and time domain decomposition, T. Carraro, M. Geiger, S. Körkel, and R. Rannacher, eds., Springer-Verlag, Berlin, 2015, pp. 69–114.
  * [14] G.H. Golub and C.F. Van Loan,  Matrix Computations, third edition,  The John Hopkins University Press, Baltimore, Maryland, 1996
  * [15] R. Griesmaier, and P. Monk,  Discretization of the Wave Equation Using Continuous Elements in Time and a Hybridizable Discontinuous Galerkin Method in Space,  J. Sci. Comput., 58: 472–498, 2014.
  * [16] E. Hairer and G. Wanner,  Solving Ordinary Differential Equations II, Springer-Verlag, Berlin, Heidelberg, 1996.
  * [17] B.L. Hulme,  One-step piecewise polynomial Galerkin methods for initial value problems,  Math. Comp., 26: 415–426, 1972.
  * [18] O. Perron,  Algebra, Volume II,  Guschens Bücherei, Berlin, 1933.
  * [19] T. Richter, A. Springer, and Vexler, Efficient numerical realization of discontinuous Galerkin methods for temporal discretization of parabolic problems, Numer. Math., 124: 151–182, 2013.
  * [20] E. B. Saff and R. S. Varga,  On the zeros and Poles of Padé Approximation to $e^{z}$,  Numer. Math., 25: 1–14, 1975.
  * [21] D. Schötzau and Ch. Schwab,  p- and hp- Finite Element Methods, Oxford Science Publications, New York, 1998.
  * [22] D. Schötzau and Ch. Schwab,  The time discretization of parabolic problems by the hp-version of the discontinuous Galerkin finite element method,  SIAM J. Numer. Anal., 38: 837–875, 2000.
  * [23] B.N. Southworth, O. Krzysik, W. Pazner, and H. De Sterck,  Fast solution of fully implicit Runge-Kutta and discontinuous Galerkin in time for numerical PDEs, Part I: The linear setting,  SIAM J. Sci. Comput., 44: A416–A443, 2022.
  * [24] G. Szegö,  Orthogonal Polynomials,  American Mathematical Society, New York, 1939.
  * [25] Y. Wu and Y. Bai,  Error analysis of energy-preserving mixed finite element methods for the Hodge wave equation,  SIAM J. Numer. Anal., 59: 1433–1454, 2021.