# Stability of conforming space–time
isogeometric methods for the wave equation

Matteo Ferrari
Faculty of Mathematics
University of Vienna
1090 Vienna, Austria
<EMAIL_ADDRESS>
& Sara Fraschini
Faculty of Mathematics
University of Vienna
1090 Vienna, Austria
<EMAIL_ADDRESS>

###### Abstract

We consider a family of conforming space–time finite element discretizations
for the wave equation based on splines of maximal regularity in time.
Traditional techniques may require a CFL condition to guarantee stability.
Recent works by O. Steinbach and M. Zank (2018), and S. Fraschini, G. Loli, A.
Moiola, and G. Sangalli (2023), have introduced unconditionally stable schemes
by adding non-consistent penalty terms to the underlying bilinear form.
Stability and error analysis have been carried out for lowest order discrete
spaces. While higher order methods have shown promising properties through
numerical testing, their rigorous analysis was still missing. In this paper,
we address this stability analysis by studying the properties of the condition
number of a family of matrices associated with the time discretization. For
each spline order, we derive explicit estimates of both the CFL condition
required in the unstabilized case and the penalty term that minimises the
consistency error in the stabilized case. Numerical tests confirm the
sharpness of our results.

##  1 Introduction

The foundations of space–time Finite Element Methods (FEMs) for wave
propagation problems have been laid in the seminal paper [15]. Unlike standard
time-stepping techniques, space–time methods involve the simultaneous
discretization of space and time variables, yielding numerical solutions
available at all time instants. The popularity of these space–time schemes has
been promoted in the last decades by recent advancements in computer
technology, along with the remarkably properties of space–time
discretizations. These include effective handling of moving boundaries, high
order accuracy both in space and in time, multilevel preconditioning, local
mesh refinement and adaptivity, and space–time parallelization.

The majority of space–time methods for the wave equation belong to the
discontinuous Galerkin family (see, e.g. [17, 21, 22] and references therein).
In contrast, there has been much less work dedicated to space–time conforming
Galerkin methods, for which only a few works are available. Among these, we
mention the recent works [2, 9, 19, 25, 26]. It has been shown in [8, 25] that
stability of conforming space–time methods for the second order acoustic wave
equation may be subject to a CFL condition. This limitation has been overcome
in [25, 29], where an unconditionally stable space–time piecewise continuous
FEM has been designed by stabilizing the bilinear form coming from integration
by parts. In [8, 9], generalizing the idea from [25], a high order
unconditionally stable space–time isogeometric method has been devised.

The isogeometric method, also known as Isogeometric Analysis (IgA), has been
introduced in [14] as an advancement from classical FEMs, aiming to streamline
the integration of computer-aided design with numerical simulations.
Essentially, the isoparametric concept is invoked, since IgA employs spline
functions, or their generalizations, to both parametrize the computational
domain and approximate the differential equation of interest. Leveraging the
smoothness of splines leads to increased accuracy per degree of freedom when
compared to classical $C^{0}$ piecewise polynomial discretizations (see, e.g.
[3, 7]), and superior spectral approximation properties [16].

In this paper, we consider the space–time isogeometric method introduced in
[9]. We provide the first theoretical justification of its unconditional
stability. As in [25], we carry out our analysis focusing on an ordinary
differential equation strongly related to the time part of the acoustic wave
equation. Our theoretical results are derived by exploiting the algebraic
structure of the matrices involved, and are aimed at studying the properties
of their condition numbers. The analysis we perform combines two main tools:
properties of the so-called symbols [12] of the matrices associated with IgA
discretizations, and the behaviour of the condition number of banded Toeplitz
matrices characterized in [1]. All at once we are able to obtain, for each
spline degree, sharp estimates on the CFL condition required in the
unstabilized case, and derive what is the penalty term that guarantees
unconditional stability, while minimising the consistency error. Moreover, we
apply our analysis technique to different choices for the penalty term, and
show that, among those, the one proposed in [9] is the best one.

The paper is structured as follows. In Section 2, we revisit the framework
considered in [8, 9, 25], detailing the variational setting and the proposed
stabilized bilinear forms. In Section 3, we provide an overview of splines
with maximal regularity and of the properties of the associated Galerkin
matrices. In Section 4, we recall the results on the conditioning behaviour of
banded Toeplitz matrices, with emphasis on those exhibiting particular
symmetries in their coefficients. In Section 5, we present the proof of our
main result on the CFL condition and the optimal penalty term. In Section 6,
we consider and analyze other possible stabilizations terms. Finally, in
Section 7, we present numerical tests that demonstrate that our estimates are
sharp.

##  2 Model problem and stabilized variational formulations

We consider the linear acoustic wave equation defined on the space–time
cylinder $Q=\Omega\times(0,T)$,

$\begin{cases}\partial_{t}^{2}U(\boldsymbol{x},t)-\Delta_{\boldsymbol{x}}U(\boldsymbol{x},t)=F(\boldsymbol{x},t),&(\boldsymbol{x},t)\in
Q,\\\
U(\boldsymbol{x},t)=0,&(\boldsymbol{x},t)\in\partial\Omega\times(0,T),\\\
U(\boldsymbol{x},0)=0,\quad\partial_{t}U(\boldsymbol{x},t)_{|_{t=0}}=0,&\boldsymbol{x}\in\Omega,\end{cases}$
(1)

where $\Omega\subset\mathbb{R}^{d}$ ($d\in\mathbb{N}$) is a bounded Lipschitz
domain, $T>0$ is a finite time, and $F\in L^{2}(Q)$ is a given source term.

To recall the space–time variational formulations considered in [9, 25], we
introduce the Hilbert spaces

$\displaystyle H^{1,1}_{0;0,\bullet}(Q)$
$\displaystyle=L^{2}\big{(}0,T;H^{1}_{0}(\Omega)\big{)}\cap
H^{1}_{0,\bullet}\big{(}0,T;L^{2}(\Omega)\big{)},\quad
H^{1,1}_{0;\bullet,0}(Q)$
$\displaystyle=L^{2}\big{(}0,T;H^{1}_{0}(\Omega)\big{)}\cap
H^{1}_{\bullet,0}\big{(}0,T;L^{2}(\Omega)\big{)},$

where we denote $H^{1}_{0}(\Omega)=\left\\{w\in H^{1}(\Omega):\
w_{|_{\partial\Omega}}=0\right\\}$, and the Sobolev-Bochner spaces

$\displaystyle H^{1}_{0,\bullet}\big{(}0,T;L^{2}(\Omega)\big{)}$
$\displaystyle=\\{W\in H^{1}\big{(}0,T;L^{2}(\Omega)\big{)}:\
W(\cdot,0)=0\\},$ $\displaystyle
H^{1}_{\bullet,0}\big{(}0,T;L^{2}(\Omega)\big{)}$ $\displaystyle=\\{W\in
H^{1}\big{(}0,T;L^{2}(\Omega)\big{)}:\ W(\cdot,T)=0\\}.$

After integration by parts, the space–time variational formulation for (1)
reads as: find $U\in H^{1,1}_{0;0,\bullet}(Q)$ such that

$-(\partial_{t}U,\partial_{t}V)_{L^{2}(Q)}+(\nabla_{\boldsymbol{x}}U,\nabla_{\boldsymbol{x}}V)_{L^{2}(Q)}=(F,V)_{L^{2}(Q)},\quad\text{
for all }V\in H^{1,1}_{0;\bullet,0}(Q),$ (2)

where $(\cdot,\cdot)_{L^{2}(Q)}$ denotes the inner product in $L^{2}(Q)$. Note
that the initial condition $U(\boldsymbol{x},0)=0$ and the homogeneous
Dirichlet condition are considered in a strong sense, while
$\partial_{t}U(\boldsymbol{x},t)_{|_{t=0}}=0$ is incorporated in the
variational formulation.

Exploiting the Fourier expansion of the trial and tests functions, it has been
shown (see, e.g. [28]) that there exists a unique solution $U\in
H^{1,1}_{0;0,\bullet}(Q)$ of the variational formulation (2), stable with
respect to the source term $F$. In greater details, let
$\\{\psi_{\ell}\\}_{\ell}$ be an orthogonal basis in
$H^{1}_{0,\bullet}(0,T)=\\{w\in H^{1}(0,T):\ w(0)=0\\}$ composed by the
eigenfunctions of the second time derivative with zero initial condition, and
let $\\{\lambda_{j}\\}_{j}$ be the eigenfunctions of the space Laplacian with
Dirichlet boundary conditions. Any $U\in H^{1,1}_{0;0,\bullet}(Q)$ admits the
representation

$U(\boldsymbol{x},t)=\sum_{j}\sum_{\ell}U_{j,\ell}\psi_{\ell}(t)\lambda_{j}(\boldsymbol{x})=\sum_{j}u_{j}(t)\lambda_{j}(\boldsymbol{x}),\quad
u_{j}(t)=\sum_{\ell}U_{j,\ell}\psi_{\ell}(t),$

with suitable coefficients $\\{U_{j,\ell}\\}_{j,\ell}\subset\mathbb{R}$.
Choosing in (2) a test function $v(t)\lambda_{j}(\boldsymbol{x})$ with $v\in
H^{1}_{\bullet,0}(0,T)=\\{w\in H^{1}(0,T):\ w(T)=0\\}$, it follows that the
variational problem (2) is equivalent to finding, for each $j$, the
coefficient functions $u_{j}\in H^{1}_{0,\bullet}(0,T)$ such that

$-(\partial_{t}u_{j},\partial_{t}v)_{L^{2}(0,T)}+\mu_{j}(u_{j},v)_{L^{2}(0,T)}=(f_{j},v)_{L^{2}(0,T)},\quad\text{
for all }v\in H^{1}_{\bullet,0}(0,T),$

where $\\{\mu_{j}\\}_{j}$ is the non-decreasing, positive and divergent
sequence of the eigenvalues of the space Dirichlet Laplacian, and
$f_{j}(t)=(F(\cdot,t),\lambda_{j})_{L^{2}(\Omega)}$.

Our main focus lies in providing a finite element discretization for the
initial value problem: find $u\in H_{0,\bullet}^{1}(0,T)$ such that

$a_{\mu}(u,v)=(f,v)_{L^{2}(0,T)}\quad\text{for all~{}}\,v\in
H^{1}_{\bullet,0}(0,T),$ (3)

with $f\in L^{2}(0,T)$, and the bilinear form
$a_{\mu}:H^{1}_{0,\bullet}(0,T)\times H^{1}_{\bullet,0}(0,T)\to\mathbb{R}$
given by

$a_{\mu}(u,v)=-(\partial_{t}u,\partial_{t}v)_{L^{2}(0,T)}+\mu(u,v)_{L^{2}(0,T)},\quad\text{for~{}}\mu>0.$
(4)

If this method is stable independently of the mesh parameter that
characterizes the finite element subspace and $\mu$, then we expect
unconditional stability for the associated space–time discretization of the
wave problem (2).

In [26, Theorem 4.1] it has been shown unique solvability of (3) by a compact
perturbation argument. However, this result is not enough to guarantee
unconditional stability (with respect to $\mu$) for each numerical
discretization of (3). It has been observed in [26, Remark 4.9] that a
piecewise-linear continuous Galerkin-Petrov FEM discretization of (3), with a
uniform mesh size $h$, is stable if and only if $h^{2}\mu<12$. This mesh
condition turns out to be a CFL of the type $h_{t}<C_{\Omega}h_{x}$ in the
associated space–time variational formulation for the wave equation, with
$C_{\Omega}>0$ depending on the domain $\Omega$, and $h_{t}$ and $h_{x}$ the
temporal and spatial mesh parameters, respectively. To overcome this mesh
condition, the bilinear form (4) is stabilized in [25]. Specifically, an
unconditionally stable Galerkin-Petrov discretization is proposed by
perturbing the bilinear form (4) with a suitable penalty term [25, Corollary
17.1], i.e., considering

$a_{\mu,h}(u_{h},v_{h})=-(\partial_{t}u_{h},\partial_{t}v_{h})_{L^{2}(0,T)}+\mu(u_{h},v_{h})_{L^{2}(0,T)}-\frac{\mu
h^{2}}{12}(\partial_{t}u_{h},\partial_{t}v_{h})_{L^{2}(0,T)}.$

For all $h$ and $\mu$, the perturbed bilinear form $a_{\mu,h}$ satisfies a
discrete inf–sup condition [25, Lemma 17.6], and ensures optimal error
estimates in $H^{1}$ and $L^{2}$ norms. An extension of this stabilization to
high order Galerkin-Petrov $C^{0}$ finite elements has been proposed and
numerically verified in [29].

In [8, 9], in the same spirit, a stabilized numerical scheme for isogeometric
discretizations with maximal regularity splines has been proposed adding a
suitable non-consistent penalty term. For a spline space of degree $p$ with
maximal regularity $C^{p-1}$, the proposed stabilized bilinear form reads as

$a^{p}_{\mu,h}(u_{h},v_{h})=-(\partial_{t}u_{h},\partial_{t}v_{h})_{L^{2}(0,T)}+\mu(u_{h},v_{h})_{L^{2}(0,T)}+\mu\delta
h^{2p}(\partial_{t}^{p}u_{h},\partial_{t}^{p}v_{h})_{L^{2}(0,T)}\,,$ (5)

for a suitable parameter $\delta\in\mathbb{R}$. The stabilized bilinear form
(5) is then employed in the space–time finite element discretization for the
wave equation (1) defining

$a^{p}_{h}(U_{h},V_{h})=-(\partial_{t}U_{h},\partial_{t}V_{h})_{L^{2}(Q)}+(\nabla_{\boldsymbol{x}}U_{h},\nabla_{\boldsymbol{x}}V_{h})_{L^{2}(Q)}+\delta
h_{t}^{2p}(\nabla_{\boldsymbol{x}}\partial_{t}^{p}U_{h},\nabla_{\boldsymbol{x}}\partial_{t}^{p}V_{h})_{L^{2}(Q)}.$
(6)

Even if the proposed stabilized isogeometric formulation of [8, 9, 10] is
supported by a wide range of numerical experiments highlighting its stability
and approximation properties, rigorous theoretical results are still missing.
In the next sections, we address this theoretical analysis.

Let us define the following two numbers, depending on the spline degree $p$:

$\delta_{p}=-\frac{(2^{2(p+1)}-1)}{2^{2p-1}\pi^{2(p+1)}}\zeta\left(2(p+1)\right),\quad\rho_{p}=4\pi^{2}\frac{(2^{2p}-1)}{(2^{2(p+1)}-1)}\frac{\zeta(2p)}{\zeta(2(p+1))},$
(7)

where $\zeta$ is the Riemann zeta function. In Theorem 5.9, we show that the
Galerkin method associated with the bilinear form (5) is stable, in a sense
that will be specified later,

  * •

for all $h$ and $\mu$, if and only if $\delta\leq\delta_{p}$,

  * •

if and only if $\mu h^{2}\leq\rho_{p}$, when $\delta=0$.

In [9, Figure 1], it is reported the consistency error of the space–time
method related to (5) with respect to the perturbation coefficients $\delta$.
As one can observe, the optimal parameters that guarantee stability scale like
$\delta\approx-10^{-p}$. This is coherent with (7), as we show that
$\delta_{p}\sim-8/\pi^{2}\cdot\pi^{-2p}\approx-0.81\cdot 9.87^{-p}$ as
$p\to\infty$. In Section 7, we numerically validate that (7) are both sharp.

##  3 Maximal regularity splines and associated matrices

In this section, we describe the maximal regularity spline spaces and leverage
certain properties of the matrices representing the bilinear form (5).

For $N\in\mathbb{N}$, we introduce the uniform mesh parameter
$h=\nicefrac{{T}}{{N}}$ and the break points $t_{\ell}=\ell h$, for
$\ell=0,\ldots,N$. For a positive integer $p\geq 1$, let
$\\{\xi_{j}\\}_{j=0}^{2p+N}\subset[0,T]$ be an open knot vector. According to
the Cox-De Boor recursion formulas, spline basis functions are piecewise
polynomials defined recursively as

$\displaystyle\varphi^{p}_{j}(t)=\begin{cases}\dfrac{t-\xi_{j}}{\xi_{j+p}-\xi_{j}}\varphi_{j}^{p-1}(t)+\dfrac{\xi_{j+p+1}-t}{\xi_{j+p+1}-\xi_{j+1}}\varphi_{j+1}^{p-1}(t)&{\text{if
}}t\in[\xi_{j},\xi_{j+p+1}),\vspace{0.3cm}\\\
0&\text{otherwise},\end{cases}\quad\quad j=0,\ldots,N+p-1,$

beginning with $\varphi^{0}_{j}(t)=1$, if $t\in[\xi_{j},\xi_{j+1})$, and
$\varphi^{0}_{j}(t)=0$ otherwise. The knot vector that guarantees maximal
$C^{p-1}$ regularity is built as

$\xi_{0}=\xi_{1}=\ldots=\xi_{p}=t_{0}=0,\quad\xi_{p+1}=t_{1},\quad\ldots\quad\xi_{p+N-1}=t_{N-1},\quad\xi_{p+N}=\ldots=\xi_{2p+N}=t_{N}=T.$

The space of splines generated by $\\{\varphi_{j}^{p}\\}_{j=0}^{N+p-1}$ is
denoted by $S_{h}^{p}(0,T)$, while the trial and test spaces for the discrete
variational formulation (3) are, respectively,

$S^{p}_{h,0,\bullet}(0,T)=S_{h}^{p}(0,T)\cap
H_{0,\bullet}^{1}(0,T)\quad\text{and}\quad
S^{p}_{h,\bullet,0}(0,T)=S_{h}^{p}(0,T)\cap H_{\bullet,0}^{1}(0,T).$

From the $p^{\operatorname{th}}$ basis function to the
$(N-1)^{\operatorname{th}}$, each spline maintains the same shape but is
supported in different intervals. We say that $\varphi_{j}^{p}$ is central if
it satisfies $p\leq j\leq N-1$. The support of the central spline
$\varphi_{j}^{p}$ is $[t_{j-p},t_{j+1}]$.

(a) $p=1$

(b) $p=2$

(c) $p=3$

Figure 1: The spline basis of $S_{h}^{p}(0,1)$ with $N=5$. The break points
are indicated by vertical dashed lines, and central splines are highlighted
with thicker lines.

We now examine the matrix representation of the stabilized bilinear form
$a_{\mu,h}^{p}$ (5) with respect to the discrete spaces
$S^{p}_{h,0,\bullet}(0,T)\times S^{p}_{h,\bullet,0}(0,T)$. Let us define the
three matrices

$\displaystyle\mathbf{M}^{p}_{h}[\ell,j]=\int_{0}^{T}$
$\displaystyle\varphi_{j}^{p}(t)\varphi^{p}_{\ell-1}(t)\differential
t,\quad\mathbf{B}^{p}_{h}[\ell,j]=\int_{0}^{T}\partial_{t}\varphi_{j}^{p}(t)\partial_{t}\varphi^{p}_{\ell-1}(t)\differential
t,$ (8)
$\displaystyle\mathbf{D}^{p}_{h}[\ell,j]=\int_{0}^{T}\partial_{t}^{p}\varphi_{j}^{p}(t)\partial_{t}^{p}\varphi^{p}_{\ell-1}(t)\differential
t$

with $\ell,j=1,\ldots N+p-1$.

######  Example 3.1.

Here we explicitly write the entries of the matrices
$\mathbf{M}_{h}^{2},\mathbf{B}_{h}^{2},\mathbf{D}_{h}^{2}\in\mathbb{R}^{(N+1)\times(N+1)}$
to show their structure. With the exception of $5$ entries at the top-left and
$5$ at the bottom-right corners, they are banded Toeplitz matrices:

$\displaystyle\mathbf{M}_{h}^{2}$
$\displaystyle=\frac{h}{120}\begin{pmatrix}14&2\\\ 40&25&1\\\ 25&66&26&1\\\
1&26&66&26&1\\\ &\ddots&\ddots&\ddots&\ddots&\ddots\\\ &&1&26&66&26&1\\\
&&&1&26&66&25&2\\\ &&&&1&25&40&14\\\ \end{pmatrix}_{(N+1)\times(N+1)},$
$\displaystyle\mathbf{B}_{h}^{2}$
$\displaystyle=\frac{1}{6h}\begin{pmatrix}-6&-2\\\ 8&-1&-1\\\ -1&6&-2&-1\\\
-1&-2&6&-2&-1\\\ &\ddots&\ddots&\ddots&\ddots&\ddots\\\ &&-1&-2&6&-2&-1\\\
&&&-1&-2&6&-1&-2\\\ &&&&-1&-1&8&-6\\\ \end{pmatrix}_{(N+1)\times(N+1)}$

and

$\mathbf{D}_{h}^{2}=\frac{1}{h^{3}}\begin{pmatrix}-6&2\\\ 10&-5&1\\\
-5&6&-4&1\\\ 1&-4&6&-4&1\\\ &\ddots&\ddots&\ddots&\ddots&\ddots\\\
&&1&-4&6&-4&1\\\ &&&1&-4&6&-5&2\\\ &&&&1&-5&10&-6\\\
\end{pmatrix}_{(N+1)\times(N+1)}.$

Let $\Phi_{p}$ denote the cardinal B-spline of degree $p$, namely the spline
function defined over the uniform knot sequence $\\{0,\ldots,p+1\\}$ as

$\Phi_{p}(s)=\frac{s}{p}\Phi_{p-1}(s)+\frac{p+1-s}{p}\Phi_{p-1}(s-1),\quad
p\geq 1$

beginning with $\Phi_{0}(t)=1$ if $t\in[0,1]$ and $\Phi_{0}(t)=0$ otherwise.
Observe that the central splines are essentially cardinal splines that takes
into account the scaling parameter $h$ and a translation. Specifically, for
$j=p,\ldots,N-1$, each central spline $\varphi_{j}^{p}(t)$ can be expressed as
$\Phi_{p}(t/h-j+p)$. Consequently, we can represent the matrices in (8) as

$\mathbf{M}^{p}_{h}[\ell,j]=\int_{0}^{T}\Phi_{p}(t/h-j+p)\Phi_{p}(t/h-\ell+1+p)\differential
t=h\int_{\mathbb{R}}\Phi_{p}(s-j)\Phi_{p}(s-\ell+1)\differential s,$ (9)

for $j=p,\ldots,N-1$ and $\ell=p+1\ldots,N$, and similarly

$\mathbf{B}^{p}_{h}[\ell,j]=\frac{1}{h}\int_{\mathbb{R}}\partial_{s}\Phi_{p}(s-j)\partial_{s}\Phi_{p}(s-\ell+1)\differential
s,\quad\mathbf{D}^{p}_{h}[\ell,j]=\frac{1}{h^{2p-1}}\int_{\mathbb{R}}\partial_{s}^{p}\Phi_{p}(s-j)\partial_{s}^{p}\Phi_{p}(s-\ell+1)\differential
s.$ (10)

In light of the previous expressions, we outline some considerations on the
scale and structures of these matrices that are consequences of well-known
results on cardinal B-splines (see [4, Section 4.2] and [11]).

######  Proposition 3.2.

For all $p\geq 1$, let $\mathbf{M}_{h}^{p},\mathbf{B}_{h}^{p}$ and
$\mathbf{D}_{h}^{p}$ defined in (8). Then, the following properties hold true.

  * •

The entries within matrices $\frac{1}{h}\mathbf{M}_{h}^{p}$,
$h\mathbf{B}_{h}^{p}$ and $h^{2p-1}\mathbf{D}_{h}^{p}$ are independent of the
mesh parameter $h$.

  * •

Except for a fixed number of entries (determined by $p$ rather than $h$ and
$T$) located at the top-left and bottom-right corners, the matrices
$\mathbf{M}_{h}^{p}$, $\mathbf{B}_{h}^{p}$, and $\mathbf{D}_{h}^{p}$ exhibit a
banded Toeplitz structure.

  * •

In their purely banded Toeplitz parts, the matrices exhibit symmetry in the
entries with respect to the lower codiagonal. Specifically, the central non-
vanishing elements can be expressed as

$\displaystyle\frac{1}{h}\mathbf{M}^{p}_{h}[\ell,\ell-1\pm
j]=\Phi_{2p+1}(p+1-j),$ $\displaystyle\quad h\mathbf{B}^{p}_{h}[\ell,\ell-1\pm
j]=-\partial_{s}^{2}\Phi_{2p+1}(p+1-j),$ (11) $\displaystyle
h^{2p-1}\mathbf{D}^{p}_{h}[\ell,\ell-1\pm j]$
$\displaystyle=(-1)^{p}\partial_{s}^{2p}\Phi_{2p+1}(p+1-j),$

for $j=0,\ldots,p$ and $\ell=2p+1,\ldots,N-p$.

###### Proof.

The first two properties readily follow from the definition of the matrices.
On the other hand, representations (11) are consequences of (9)-(10), along
with the characterization for inner products of derivatives of cardinal
B-splines [11, Lemma 4]

$\int_{\mathbb{R}}\partial^{k_{1}}_{s}\Phi_{p}(s)\partial^{k_{2}}_{s}\Phi_{p}(s+j)\differential
s=(-1)^{k_{2}}\partial^{(k_{1}+k_{2})}_{s}\Phi_{2p+1}(p+1-j)$ (12)

and their symmetry property [11, Lemma 3]

$\partial_{s}^{k}\Phi_{2p+1}(p+1+s)=(-1)^{k}\partial_{s}^{k}\Phi_{2p+1}(p+1-s).$
(13)

∎

######  Remark 3.3.

The indices of the first and last rows representing integrals involving only
central splines are the $(2p+1)^{\operatorname{th}}$ and
$(N-p)^{\operatorname{th}}$, respectively. There exists at least one row of
that type if and only if $N-p\geq 2p+1$, i.e., $N\geq 3p+1$.

The system matrix representing the bilinear form $a_{\mu,h}^{p}$, as defined
in (5), with respect to the spaces $S^{p}_{h,0,\bullet}(0,T)\times
S^{p}_{h,\bullet,0}(0,T)$ is given by

$\mathbf{K}_{\mu,h}^{p}(\delta)=-\mathbf{B}_{h}^{p}+\mu\mathbf{M}_{h}^{p}+\mu\delta
h^{2p}\mathbf{D}_{h}^{p}.$

Upon looking at the scaled matrix $h\mathbf{K}_{\mu,h}^{p}(\delta)$, in
accordance with Proposition 3.2, its dependence on $\mu$ and $h$ is solely
through the product $\mu h^{2}$. Consequently, we introduce the following
definitions

$\rho=\mu
h^{2}\,,\quad\mathbf{K}_{n}^{p}(\rho,\delta)=h\mathbf{K}_{\mu,h}^{p}(\delta),\quad
n=N+p-1$ (14)

where $\mathbf{K}_{n}^{p}(\rho,\delta)\in\mathbb{R}^{n\times n}$. By virtue of
(11), the non-zero entries corresponding to the symmetric banded Toeplitz
structure of $\mathbf{K}_{n}^{p}(\rho,\delta)$ can be expressed as

$\mathbf{K}_{n}^{p}(\rho,\delta)[\ell,\ell-1\pm
j]=\partial_{s}^{2}\Phi_{2p+1}(p+1-j)+\rho\left(\Phi_{2p+1}(p+1-j)+\delta(-1)^{p}\partial_{s}^{2p}\Phi_{2p+1}(p+1-j)\right)$
(15)

for $j=0,\ldots,p$ and $\ell=2p+1,\ldots,n-2p+1$.

In the next sections, for fixed $p\geq 1$, we investigate the spectral
properties of the family of matrices
$\\{\mathbf{K}^{p}_{n}(\rho,\delta)\\}_{n}$ with respect to the size $n$ and
to the parameters $\rho$ and $\delta$. We derive necessary and sufficient
conditions that ensure the well-conditioning (in a sense that will be
specified later) of these matrix families in two distinct cases:

  * •

we find $\delta_{p}\in\mathbb{R}$ such that
$\\{\mathbf{K}_{n}^{p}(\rho,\delta)\\}_{n}$ is well-conditioned for all
$\rho>0$ if and only if $\delta\leq\delta_{p}$;

  * •

we find $\rho_{p}\in\mathbb{R}$ such that
$\\{\mathbf{K}_{n}^{p}(\rho,0)\\}_{n}$ is well-conditioned if and only if
$\rho\leq\rho_{p}$.

##  4 Conditioning of banded Toeplitz matrices

In this section, after reviewing some results related to the conditioning of
banded Toeplitz matrices (Section 4.1), we derive an easily verifiable
necessary condition for the well-conditioning of this family of matrices when
they exhibit symmetry with respect to each lower codiagonal (Section 4.2).

###  4.1 Conditioning of general banded Toeplitz matrices

A family of non-singular matrices $\\{\mathbf{A}_{n}\\}_{n}$, where each
$\mathbf{A}_{n}\in\mathbb{R}^{n\times n}$, is said to be well-conditioned if
the condition numbers ${\kappa(\mathbf{A}_{n})}$ are uniformly bounded with
respect to $n$. It is said to be weakly well-conditioned if
${\kappa(\mathbf{A}_{n})}$ grows as a power of $n$. Here and in the subsequent
sections, the condition numbers are computed using either the $\|\cdot\|_{1}$
or $\|\cdot\|_{\infty}$ matrix norms.

Given the coefficients $\\{a_{-m},\ldots,a_{0},\ldots,a_{k}\\}$, consider the
family of banded Toeplitz matrices $\\{\mathbf{A}_{n}\\}_{n}$,

$\mathbf{A}_{n}=\begin{pmatrix}a_{0}&\ldots&a_{k}&&\\\
\vdots&\ddots&&\ddots\\\ a_{-m}&&\ddots&&a_{k}\\\ &\ddots&&\ddots&\vdots\\\
&&a_{-m}&\ldots&a_{0}\\\ \end{pmatrix}_{n\times n}.$ (16)

We associate to this family the polynomial
$p^{\mathbf{A}}\in\mathbb{P}_{k+m}(\mathbb{R})$

$p^{\mathbf{A}}(z)=\sum_{j=-m}^{k}a_{j}z^{m+j}.$ (17)

Without loss of generality, we assume that $a_{-m}a_{k}\neq 0$. In [1], the
authors show that the well-conditioning (or weakly well-conditioning) of the
family $\\{\mathbf{A}_{n}\\}_{n}$ depends on the location of zeros of
$p^{\mathbf{A}}$.

We define a polynomial as being of type $(s,u,\ell)$ if it has $s$ zeros with
modulus smaller than $1$, $u$ zeros with unit modulus, and $\ell$ zeros with
modulus larger than $1$.

######  Theorem 4.1.

[1, Theorem 3] Let $\\{\mathbf{A}_{n}\\}_{n}$ represent a family of non-
singular banded Toeplitz matrices, and let $p^{\mathbf{A}}$ in (17) be the
associated polynomial. Then, the family $\\{\mathbf{A}_{n}\\}_{n}$ exhibits
the following behaviour:

  * •

it is well-conditioned if $p^{\mathbf{A}}$ is of type $(m,0,k)$;

  * •

it is weakly well-conditioned if $p^{\mathbf{A}}$ is of type $(m_{1},m_{2},k)$
or $(m,k_{1},k_{2})$, where $m_{1}+m_{2}=m$ and $k_{1}+k_{2}=k$. In this
scenario, $\kappa(\mathbf{A}_{n})$ grows at most as $\mathcal{O}(n^{\eta})$,
with $\eta$ representing the highest multiplicity among the zeros of unit
modulus of $p^{\mathbf{A}}$.

Additionally, if $p^{\mathbf{A}}$ has more than $k$ roots outside the unit
circle or more than $m$ roots inside it, then $\kappa(\mathbf{A}_{n})$ grows
exponentially as $n$ increases.

The proof of Theorem 4.1 relies on an auxiliary result. Let
$\\{z_{1},\ldots,z_{m+k}\\}$ denote the zeros of $p^{\mathbf{A}}$, arranged
such that $0<|z_{1}|\leq\ldots\leq|z_{m}|\leq\ldots\leq|z_{m+k}|$. We define
the polynomials $p_{\ell}\in\mathbb{P}_{m}(\mathbb{R})$ and
$p_{u}\in\mathbb{P}_{k}(\mathbb{R})$ as follows

$p_{\ell}(z)=\prod_{j=1}^{m}(z-z_{j})=\sum_{j=-m}^{0}\ell_{j}z^{m+j}\quad\quad\text{and}\quad\quad
p_{u}(z)=a_{k}\prod_{j=1}^{k}(z-z_{m+j})=\sum_{j=0}^{k}u_{j}z^{j},$

along with the Toeplitz triangular matrices

$\displaystyle\mathbf{L}_{n}=\begin{pmatrix}\ell_{0}\\\ \vdots&\ddots\\\
\ell_{-m}&&\ddots\\\ &\ddots&&\ddots\\\ &&\ell_{-m}&\dots&\ell_{0}\\\
\end{pmatrix}_{n\times
n}\quad\text{and}\quad\mathbf{U}_{n}=\begin{pmatrix}u_{0}&\dots&u_{k}\\\
&\ddots&&\ddots\\\ &&\ddots&&u_{k}\\\ &&&\ddots&\vdots\\\ &&&&u_{0}\\\
\end{pmatrix}_{n\times n}.$

When $p^{\mathbf{A}}$ is of type $(m_{1},m_{2},k)$ or $(m,k_{1},k_{2})$, where
$m_{1}+m_{2}=m$ or $k_{1}+k_{2}=k$, respectively, it can be shown (see [1,
Lemma 2]), that
$\mathbf{A}_{n}^{-1}=\mathbf{U}_{n}^{-1}\mathbf{L}_{n}^{-1}+\mathbf{H}_{n}$,
where for sufficiently large $n$, the $(j_{1},j_{2})^{\text{th}}$ entry of
$\mathbf{H}_{n}$ behaves like

$\mathcal{O}\left(j_{1}^{\xi-1}\left|z_{m}\right|^{j_{1}}j_{2}^{\nu-1}\left|z_{m+1}\right|^{-j_{2}}\right)$

where $\xi=\max\\{\xi_{j}:\xi_{j}\text{~{}multiplicity
of~{}}z_{j},|z_{j}|=|z_{m}|\leq 1\\}$ and
$\nu=\max\\{\nu_{j}:\nu_{j}\text{~{}multiplicity
of~{}}z_{j},|z_{j}|=|z_{m+1}|\geq 1\\}.$ This result shows that the well-
conditioning of $\mathbf{L}_{n}$ and $\mathbf{U}_{n}$ guarantees the well-
conditioning of $\mathbf{A}_{n}$.

Let us now consider another scenario. Suppose that two blocks,
$\mathbf{C}_{1}$ and $\mathbf{C}_{2}$, possibly different, are present in
every matrix within the family $\\{\widetilde{\mathbf{A}}_{n}\\}_{n}$:
$\mathbf{C}_{1}$ in the top-left corner and $\mathbf{C}_{2}$ the bottom-right
corner

$\widetilde{\mathbf{A}}_{n}=\left(\begin{array}[]{c|c|c}\makebox[0.0pt]{$\mathbf{C}_{1}$}&\leavevmode\hbox{\set@color$\begin{matrix}&&\\\
a_{k}&&\\\ &\ddots&\end{matrix}$}&\makebox[0.0pt]{$\phantom{A_{1}}$}\\\
\hline\cr\leavevmode\hbox{\set@color$\begin{matrix}&a_{-m}&\\\ &&\ddots\\\
&&\end{matrix}$}&\leavevmode\hbox{\set@color$\begin{matrix}a_{0}&&\ddots\\\
&\ddots&\\\
\ddots&&a_{0}\end{matrix}$}&\leavevmode\hbox{\set@color$\begin{matrix}&&\\\
\ddots&&\\\ \vspace{-0.3cm}&a_{k}&\end{matrix}$}\\\
\hline\cr\makebox[0.0pt]{$\phantom{A_{1}}$}&\leavevmode\hbox{\set@color$\begin{matrix}&\,\,\,\,\ddots&\\\
&&a_{-m}\\\ &&\end{matrix}$}&\makebox[0.0pt]{$\mathbf{C}_{2}$}\\\
\end{array}\right)_{n\times n}.$ (18)

Furthermore, we assume that these blocks $\mathbf{C}_{1}$ and $\mathbf{C}_{2}$
are independent of $n$. In this scenario, it can be proved that
$\widetilde{\mathbf{A}}_{n}^{-1}=\mathbf{U}_{n}^{-1}\mathbf{L}_{n}^{-1}+\widetilde{\mathbf{H}}_{n}$,
with $\mathbf{L}_{n}$ and $\mathbf{U}_{n}$ as before, and the
$(j_{1},j_{2})^{\operatorname{th}}$ entry of $\widetilde{\mathbf{H}}_{n}$ that
exhibits the behaviour

$\mathcal{O}\left(j_{1}^{\xi-1}\left|z_{m}\right|^{j_{1}}j_{2}^{\nu-1}\left|z_{m+1}\right|^{-j_{2}}\right)+\mathcal{O}\left(j_{2}^{\xi-1}\left|z_{m}\right|^{j_{2}}{j_{1}}^{\nu-1}\left|z_{m+1}\right|^{-j_{1}}\right)$
(19)

with $\xi$ and $\nu$ as defined previously. In fact, based on the proof of [1,
Lemma 2], it directly follows that the lemma remains valid when perturbing
each $\mathbf{A}_{n}$ only with $\mathbf{C}_{1}$ in the top-left corner.
Furthermore, adding $\mathbf{C}_{2}$ to the bottom-right corner is equivalent
to adding a perturbation only to the top-left corner and then transposing it,
and this generates the second term of (19). Finally, Theorem 4.1 also holds
true for the family of perturbed matrices
$\\{\widetilde{\mathbf{A}}_{n}\\}_{n}$, yielding the following corollary.

######  Corollary 4.2.

Let $\\{\widetilde{\mathbf{A}}_{n}\\}_{n}$ represent a family of matrices like
in (18). Then, the family $\\{\widetilde{\mathbf{A}}_{n}\\}_{n}$ exhibits the
same behaviour as the family $\\{\mathbf{A}_{n}\\}_{n}$ defined in (16).

Figure 2: The behaviour of $|{\widetilde{\mathbf{H}}}|$ (left plot) and
$|\mathbf{H}|$ (right plot) corresponding to the matrix
$\mathbf{K}^{2}_{\mu,h}(0)$ (refer to Example 3.1) and its purely banded
Toeplitz version, respectively, with $N=200$, $T=1$ and $\mu=180000$.

######  Remark 4.3.

Note that $\mathbf{A}_{n}-\widetilde{\mathbf{A}}_{n}$ is a low-rank correction
term. Specifically, if $\mathbf{C}_{1}\in\mathbb{R}^{r_{1}\times
g_{1}},\mathbf{C}_{2}\in\mathbb{R}^{r_{2}\times g_{2}}$, then
$\rank(\mathbf{A}_{n}-\widetilde{\mathbf{A}}_{n})\leq\min\\{r_{1}+r_{2},g_{1}+g_{2}\\}$
independent on $n$. Moreover, the family of matrices
$\\{\mathbf{K}_{n}^{p}(\rho,\delta)\\}_{n}$ defined in (14) follows the
structure outlined in (18) with
$\mathbf{C}_{1}=\mathbf{C}_{2}^{\top}\in\mathbb{R}^{2p\times(p-1)}$.

###  4.2 Conditioning of banded Toeplitz matrices with symmetries

For a given $p\geq 1$ and a set of $p+1$ real numbers $\\{k_{j}\\}_{j=0}^{p}$,
we consider the family of banded Toeplitz matrices denoted by
$\\{\mathbf{K}_{n}^{p}\\}_{n}$, with

$\mathbf{K}_{n}^{p}=\begin{pmatrix}k_{p-1}&k_{p-2}&\dots&k_{1}&k_{0}&&&\\\
k_{p}&\ &\ &\ &\ &\ddots\\\ k_{p-1}&\ &\ &\ &\ &\ &\ddots\\\ \vdots&\ &\ &\ &\
&\ &&k_{0}\\\ k_{1}&\ &\ &\ &\ &\ &&k_{1}\\\ k_{0}&\ &\ &\ &\ &\ &&\vdots\\\ \
&\ddots&\ &\ &\ &\ &&k_{p-2}\\\ \ &\
&k_{0}&k_{1}&\dots&k_{p-1}&k_{p}&k_{p-1}\end{pmatrix}_{n\times n}.$ (20)

These matrices align with the notations established in Theorem 4.1, with
$m=p+1$ and $k=p-1$. Consequently, investigating the conditioning behaviour of
such matrices entails analyzing the associated polynomial of degree $2p$

$q^{\mathbf{K}^{p}}(z)=k_{0}+k_{1}z+\ldots+k_{p}z^{p}+\ldots+k_{1}z^{2p-1}+k_{0}z^{2p}.$
(21)

Leveraging the specific structure of the polynomial $q^{\mathbf{K}^{p}}$, we
deduce the following result.

######  Lemma 4.4.

The family of matrices $\\{\mathbf{K}_{n}^{p}\\}_{n}$ in (20) is weakly well-
conditioned if and only if the associated polynomial $q^{\mathbf{K}^{p}}$, as
defined in (21), has exactly two zeros of unit modulus.

###### Proof.

Let us assume $k_{0}\neq 0$. According to Theorem 4.1, the matrix family
$\\{\mathbf{K}_{n}^{p}\\}_{n}$ exhibits weakly well-conditioned, or well-
conditioned, behaviour if and only if $q^{\mathbf{K}^{p}}$ falls into one of
the following types: $(r_{1},p+1-r_{1},p-1)$ for some
$r_{1}\in\\{0,\ldots,p+1\\}$ or $(p+1,p-1-r_{2},r_{2})$ for some
$r_{2}\in\\{0,\ldots,p-1\\}$. It is worth noting that if $\xi$ is a zero of
$q^{\mathbf{K}^{p}}$, then $\xi^{-1}$ is also a zero. Hence,
$q^{\mathbf{K}^{p}}$ necessarily is a polynomial of type $(r,2p-2r,r)$ for
some $r\in\\{0,\ldots,p\\}$. Consequently, the only admissible choice for
weakly well-conditioning is $q^{\mathbf{K}^{p}}$ of type $(p-1,2,p-1)$. If
$k_{0}=0$ and $k_{1}\neq 0$, then $q^{\mathbf{K}^{p}}(z)/z$ must be of type
$(p-2,2,p-2)$, and so forth. ∎

We aim to investigate the distribution of zeros of the polynomial
$q^{\mathbf{K}^{p}}$. This polynomial belongs to the category of self-
reciprocal polynomials (refer to [18]). In other words, $q^{\mathbf{K}^{p}}$
can be expressed as $q^{\mathbf{K}^{p}}(z)=\sum_{j=0}^{2p}k_{j}z^{j}$, where
$k_{j}\in\mathbb{R}$ and it exhibits symmetry $k_{j}=k_{2p-j}$ for
$j=0,\dots,p$. These polynomials can be represented as

$q^{\mathbf{K}^{p}}(z)=z^{p}\left[k_{0}\left(\frac{1}{z^{p}}+z^{p}\right)+k_{1}\left(\frac{1}{z^{p-1}}+z^{p-1}\right)+\ldots+k_{p-1}\left(\frac{1}{z}+z\right)+k_{p}\right].$

Moreover, each term $z^{-j}+z^{j}$ can be expressed as a polynomial in the new
variable $x(z)=z^{-1}+z$. Specifically, let $T_{j}$ denote the
$j^{\operatorname{th}}$ Chebyshev polynomial defined via the recurrence
relation [20, Section 1.2.1]

$T_{0}(x)=1,\quad T_{1}(x)=x,\quad
T_{j+1}(x)=2xT_{j}(x)-T_{j-1}(x),\quad\text{for }j\geq 1.$

This recursive formula has the closed form
$T_{j}(x)=\bigl{(}(x+\sqrt{x^{2}-1})^{j}+(x-\sqrt{x^{2}-1})^{j}\bigr{)}/2$.
Therefore, by substituting $x(z)=z^{-1}+z$, we obtain
$2T_{j}(\nicefrac{{x(z)}}{{2}})=z^{-j}+z^{j}$. Exploiting this property, we
introduce the so-called Chebyshev transform
$\mathcal{T}q^{\mathbf{K}^{p}}\in\mathbb{P}_{p}(\mathbb{R})$ of the self-
reciprocal polynomial $q^{\mathbf{K}^{p}}$ in (21), defined as

$\mathcal{T}q^{\mathbf{K}^{p}}(x)=2k_{0}T_{p}\left(\frac{x}{2}\right)+2k_{1}T_{p-1}\left(\frac{x}{2}\right)+\ldots+2k_{p-1}T_{1}\left(\frac{x}{2}\right)+k_{p}T_{0}\left(\frac{x}{2}\right).$
(22)

By definition, we have
$\mathcal{T}q^{\mathbf{K}^{p}}(x(z))=\nicefrac{{q^{\mathbf{K}^{p}}(z)}}{{z^{p}}}$
where $x(z)=z^{-1}+z$. Additionally, observe that the leading coefficient of
$\mathcal{T}q^{\mathbf{K}^{p}}$ is $k_{0}$. This follows from the well-known
fact that the leading coefficient of $T_{p}(x)$ is $2^{p-1}$ (see [24,
Equation 1.14])

As noted in [27, Theorem 3], for every pair of roots $(\xi,\xi^{-1})$ of
$q^{\mathbf{K}^{p}}(z)$ located on the unit disk, there exists a corresponding
zero of the polynomial $\mathcal{T}q^{\mathbf{K}^{p}}(x)$ within the interval
$[-2,2]$ on the real line. This leads us to the following conclusion.

######  Lemma 4.5.

The family of matrices $\\{\mathbf{K}_{n}^{p}\\}_{n}$ as defined in (20)
exhibits weakly well-conditioned behaviour if and only if the Chebyshev
transform $\mathcal{T}q^{\mathbf{K}^{p}}$ (22) of the polynomial
$q^{\mathbf{K}^{p}}$ in (21) has exactly one real root in the interval
$[-2,2]$.

Let $\\{x_{1},\ldots,x_{p}\\}$ be the roots of
$\mathcal{T}q^{\mathbf{K}^{p}}$. A necessary condition such that only one root
lies in the interval $[-2,2]$ reads as

$\prod_{j=1}^{p}(x_{j}^{2}-4)\leq 0.$ (23)

This condition becomes evident when all the roots are real, and even in the
case of a pair of complex conjugate roots $x$ and $\bar{x}$, (23) holds true
since

$(x^{2}-4)(\bar{x}^{2}-4)=(\text{Re}(x^{2})-4)^{2}+\text{Im}(x^{2})^{2}\geq
0.$

In the following, we obtain a simple criterion for evaluating the product in
(23).

######  Lemma 4.6.

Consider
$q^{\mathbf{K}^{p}}(z)=k_{0}z^{2p}+k_{1}z^{2p-1}+\ldots+k_{1}z+k_{0}$, a self-
reciprocal polynomial of degree $2p$. Let $\\{x_{1},\ldots,x_{p}\\}$ denote
the roots of its Chebyshev transform $\mathcal{T}q^{\mathbf{K}^{p}}$. Assuming
$k_{0}\neq 0$, the following formulas are valid

$\displaystyle\prod_{j=1}^{p}(x_{j}-2)=(-1)^{p}\,\frac{q^{\mathbf{K}^{p}}(1)}{k_{0}},\quad\quad\quad\prod_{j=1}^{p}(x_{j}+2)=\frac{q^{\mathbf{K}^{p}}(-1)}{k_{0}}.$
(24)

###### Proof.

To prove the first formula, let us consider
$g(x)=\mathcal{T}q^{\mathbf{K}^{p}}(x+2)$. For all roots $x_{j}$ of
$\mathcal{T}q^{\mathbf{K}^{p}}$, we have
$g(x_{j}-2)=\mathcal{T}q^{\mathbf{K}^{p}}(x_{j})=0$. Therefore, recalling
Vieta’s formulas the product $(x_{1}-2)\dots(x_{p}-2)$ is equal to
$(-1)^{p}g(0)/k_{0}$. We can compute
$g(0)=\mathcal{T}q^{\mathbf{K}^{p}}(2)=q^{\mathbf{K}^{p}}(1)$, the latter
because $x(z)=z^{-1}+z=2$ if and only if $z=1$. For the second formula, the
proof follows a similar reasoning. ∎

We conclude collecting the results from the previous lemmas.

######  Proposition 4.7.

If the family of matrices $\\{\mathbf{K}_{n}^{p}\\}_{n}$ in (20) is weakly
well-conditioned, then the associated polynomial $q^{\mathbf{K}^{p}}$, as
defined in (21), satisfies the following condition

$(-1)^{p}q^{\mathbf{K}^{p}}(-1)q^{\mathbf{K}^{p}}(1)\leq 0.$ (25)

###### Proof.

In case $k_{0}\neq 0$, the result follows from Lemmas 4.5 and 4.6. Conversely,
if $k_{0}=0$, analogous formulas to (24) can be derived for the roots of the
polynomial $\mathcal{T}q^{\mathbf{K}^{p}}(z)/z$ of degree $p-1$. The claim is
then established in a similar way. ∎

In general, (25) is clearly not sufficient. Indeed, if
$q^{\mathbf{K}^{p}}(1)=0$, then this condition is always satisfied. In the
next section, we will explicitly calculate the quantities
$q^{\mathbf{K}^{p}}(\pm 1)$ associated to the matrices (14), and we will show
that (25) is also sufficient for the well-conditioning in this particular
case.

##  5 Conditioning of matrices associated with maximal regularity splines

In this section, we present and prove the key results of this paper. We
analyze the polynomial $q^{\mathbf{K}^{p}(\rho,\delta)}$ associated with the
family of matrices $\\{\mathbf{K}_{n}^{p}(\rho,\delta)\\}_{n}$ defined in
(14). We investigate the conditions on $\rho$ and $\delta$ necessary for (25)
to hold true. Furthermore, we establish that, within this framework, this
condition not only proves necessary for the well-conditioning of the family
$\\{\mathbf{K}_{n}^{p}(\rho,\delta)\\}_{n}$, but is also sufficient.

Recalling the explicit formula (15), the polynomial
$q^{\mathbf{K}^{p}(\rho,\delta)}$ is given by

$\displaystyle q^{\mathbf{K}^{p}(\rho,\delta)}(z)$
$\displaystyle=k_{0}^{p}(\rho,\delta)+k_{1}^{p}(\rho,\delta)z+\ldots+k_{p}^{p}(\rho,\delta)z^{p}+\ldots+k_{1}^{p}(\rho,\delta)z^{2p-1}+k_{0}^{p}(\rho,\delta)z^{2p}$
$\displaystyle=z^{p}\left(k_{p}^{p}(\rho,\delta)+\sum_{j=1}^{p}k_{p-j}^{p}(\rho,\delta)(z^{-j}+z^{j})\right)$

where the coefficients are determined in terms of the cardinal spline
$\Phi_{2p+1}$ as

$k_{p-j}^{p}(\rho,\delta)=\partial_{s}^{2}\Phi_{2p+1}(p+1-j)+\rho\bigl{(}\Phi_{2p+1}(p+1-j)+\delta(-1)^{p}\partial_{s}^{2p}\Phi_{2p+1}(p+1-j)\bigr{)}$

for $j=0,\ldots,p$.

To investigate the number of zeros of $q^{\mathbf{K}^{p}(\rho,\delta)}$ on the
unit disk, we express $q^{\mathbf{K}^{p}(\rho,\delta)}(z)$ explicitly for
$z=e^{\mathrm{i}\theta}$ and $\theta\in[-\pi,\pi]$.

######  Proposition 5.1.

For $p\geq 1$ and $\theta\in[-\pi,\pi]$, let the functions
$M_{p},B_{p},D_{p}:[-\pi,\pi]\to\mathbb{R}$ be defined as

$\displaystyle M_{p}(\theta)$
$\displaystyle=(2-2\cos\theta)^{p+1}\sum_{j\in\mathbb{Z}}\frac{1}{(\theta+2j\pi)^{2p+2}},$
$\displaystyle B_{p}(\theta)$
$\displaystyle=-(2-2\cos\theta)^{p+1}\sum_{j\in\mathbb{Z}}\frac{1}{(\theta+2j\pi)^{2p}},$
$\displaystyle D_{p}(\theta)$ $\displaystyle=(-1)^{p}(2-2\cos\theta)^{p}.$

Then, the following relations hold true

$\displaystyle\Phi_{2p+1}(p+1)+\sum_{j=1}^{p}\Phi_{2p+1}(p+1-j)(e^{-\mathrm{i}j\theta}+e^{\mathrm{i}j\theta})$
$\displaystyle=M_{p}(\theta),$ (26)
$\displaystyle\partial_{s}^{2}\Phi_{2p+1}(p+1)+\sum_{j=1}^{p}\partial_{s}^{2}\Phi_{2p+1}(p+1-j)(e^{-\mathrm{i}j\theta}+e^{\mathrm{i}j\theta})$
$\displaystyle=B_{p}(\theta),$ (27)
$\displaystyle\partial_{s}^{2p}\Phi_{2p+1}(p+1)+\sum_{j=1}^{p}\partial_{s}^{2p}\Phi_{2p+1}(p+1-j)(e^{-\mathrm{i}j\theta}+e^{\mathrm{i}j\theta})$
$\displaystyle=D_{p}(\theta).$ (28)

Therefore, we can express
$q^{\mathbf{K}^{p}(\rho,\delta)}(e^{\mathrm{i}\theta})$ as

$q^{\mathbf{K}^{p}(\rho,\delta)}(e^{\mathrm{i}\theta})=e^{\mathrm{i}p\theta}\left(B_{p}(\theta)+\rho\left(M_{p}(\theta)+\delta(-1)^{p}D_{p}(\theta)\right)\right).$
(29)

###### Proof.

There is already a proof for the explicit expressions of $M_{p}(\theta)$ and
$B_{p}(\theta)$ in [11].For the sake of completeness, we report here the
details for $M_{p}(\theta)$, and we present a different proof for
$B_{p}(\theta)$, which is useful for calculating $D_{p}(\theta)$.

We recall that the Fourier transform of a cardinal spline $\Phi_{p}$ is,

$\widehat{\Phi}_{p}(\xi)=\left(\frac{1-e^{-\mathrm{i}\xi}}{\mathrm{i}\xi}\right)^{p+1},\quad\left|\widehat{\Phi}_{p}(\xi)\right|^{2}=\left(\frac{2-2\cos\xi}{\xi^{2}}\right)^{p+1},\quad\quad\xi\in\mathbb{R}$

(see [4, Example 3.4]), and satisfies the hypothesis of a Poisson summation
formula [4, Theorem 2.28], i.e., for all $\theta\in[-\pi,\pi]$

$\sum_{j\in\mathbb{Z}}\left(\int_{\mathbb{R}}\Phi_{p}(s)\Phi_{p}(s+j)\differential
s\right)e^{\mathrm{i}j\theta}=\sum_{j\in\mathbb{Z}}\left|\widehat{\Phi}_{p}(\theta+2j\pi)\right|^{2}.$

To show (26), recalling (12) and (13), we compute for $\theta\in[-\pi,\pi]$,

$\displaystyle\Phi_{2p+1}(p+1)$
$\displaystyle+\sum_{j=1}^{p}\Phi_{2p+1}(p+1-j)(e^{-\mathrm{i}j\theta}+e^{\mathrm{i}j\theta})=\sum_{j\in\mathbb{Z}}\Phi_{2p+1}(p+1-j)e^{\mathrm{i}j\theta}$
$\displaystyle=\sum_{j\in\mathbb{Z}}\left(\int_{\mathbb{R}}\Phi_{p}(s)\Phi_{p}(s+j)\differential
s\right)e^{\mathrm{i}j\theta}=\sum_{j\in\mathbb{Z}}\left|\widehat{\Phi}_{p}\left(\theta+2j\pi\right)\right|^{2}$
$\displaystyle=(2-2\cos\theta)^{p+1}\sum_{j\in\mathbb{Z}}\frac{1}{(\theta+2j\pi)^{2(p+1)}}.$

Using the recurrence relation for derivatives of cardinal splines [4, Theorem
4.3]

$\partial_{s}^{k}\Phi_{p}(s)=\partial_{s}^{k-1}\Phi_{p-1}(s)-\partial_{s}^{k-1}\Phi_{p-1}(s-1),$

and considering that
$e^{\mathrm{i}j\theta}-2e^{\mathrm{i}(j+1)\theta}+e^{\mathrm{i}(j+2)\theta}=e^{\mathrm{i}j\theta}(1-e^{\mathrm{i}\theta})^{2}=-e^{\mathrm{i}(j+1)\theta}(2-2\cos\theta)$,
we can compute

$\displaystyle\partial_{s}^{2}\Phi_{2p+1}(p+1)$
$\displaystyle+\sum_{j=1}^{p}\partial_{s}^{2}\Phi_{2p+1}(p+1-j)(e^{-\mathrm{i}j\theta}+e^{\mathrm{i}j\theta})$
$\displaystyle=\sum_{j\in\mathbb{Z}}\partial_{s}^{2}\Phi_{2p+1}(p+1-j)e^{\mathrm{i}j\theta}=\sum_{j\in\mathbb{Z}}\left(\partial_{s}\Phi_{2p}(p+1-j)-\partial_{s}\Phi_{2p}(p-j)\right)e^{\mathrm{i}j\theta}$
$\displaystyle=\sum_{j\in\mathbb{Z}}\left(\Phi_{2p-1}(p+1-j)-2\Phi_{2p-1}(p-j)+\Phi_{2p-1}(p-j-1)\right)e^{\mathrm{i}j\theta}$
$\displaystyle=-(2-2\cos\theta)\sum_{j\in\mathbb{Z}}\Phi_{2p-1}(p+1-j)e^{\mathrm{i}(j+1)\theta}$
$\displaystyle=-(2-2\cos\theta)\sum_{j\in\mathbb{Z}}\Phi_{2p-1}(p-j)e^{\mathrm{i}j\theta}=-(2-2\cos\theta)^{p+1}\sum_{j\in\mathbb{Z}}\frac{1}{(\theta+2j\pi)^{2p}},$

where in the last equality we used the same procedure as in the proof for
$M_{p}(\theta)$. Similarly, we can iteratively make the same steps as before
to show that

$\displaystyle\partial_{s}^{2p}\Phi_{2p+1}(p+1)+\sum_{j=1}^{p}\partial_{s}^{2p}\Phi_{2p+1}(p+1-j)(e^{-\mathrm{i}j\theta}+e^{\mathrm{i}j\theta})$
$\displaystyle=(-1)^{p}(2-2\cos\theta)^{p}\sum_{j\in\mathbb{Z}}\Phi_{1}(1-j)e^{\mathrm{i}j\theta}$
$\displaystyle=(-1)^{p}(2-2\cos\theta)^{p}.\qed$

######  Remark 5.2.

We observe that for $p\geq 1$, we have obtained
$B_{p}(\theta)=-(2-2\cos\theta)M_{p-1}(\theta)$, where we set
$M_{0}(\theta)\equiv 1$.

######  Remark 5.3.

The statement of Proposition 5.1 exhibits a close connection to [11, Lemma 6]
and [11, Lemma 7], where the authors investigate the symbols of matrices
associated with IgA discretizations. Specifically, the $p^{\operatorname{th}}$
symbol is $e_{p}(\theta)=f_{p}(\theta)/g_{p}(\theta)$, as detailed in, for
example, [6, Equation 15], where $f_{p}(\theta)=-B_{p}(\theta)$ and
$g_{p}(\theta)=M_{p}(\theta)$. For an extensive discussion on how symbols are
utilized to derive analytical spectral properties for IgA discretization of
boundary value problems, see [12] and the references provided therein.

Now, we determine $q^{\mathbf{K}^{p}(\rho,\delta)}(\pm 1)$ by computing
$M_{p}(0),B_{p}(0),D_{p}(0)$ and $M_{p}(\pi),B_{p}(\pi),D_{p}(\pi)$.

######  Corollary 5.4.

For $p\geq 1$, let $M_{p},B_{p},D_{p}$ be the functions defined in Proposition
5.1. It holds true that $M_{p}(0)=1$, and $B_{p}(0)=D_{p}(0)=0$. Furthermore,
we have

$\displaystyle M_{p}(\pi)=\frac{2}{\pi^{2(p+1)}}(2^{2(p+1)}-1)$
$\displaystyle\zeta(2(p+1)),\quad
B_{p}(\pi)=-2^{2}M_{p-1}(\pi)=-\frac{2^{3}}{\pi^{2p}}(2^{2p}-1)\zeta(2p)$
$\displaystyle\quad D_{p}(\pi)=(-1)^{p}2^{2p},$

where $\zeta$ is the Riemann zeta function. Therefore, we have
$q^{\mathbf{K}^{p}(\rho,\delta)}(1)=\rho$ and

$q^{\mathbf{K}^{p}(\rho,\delta)}(-1)=(-1)^{p}\left(-\frac{2^{3}}{\pi^{2p}}(2^{2p}-1)\zeta(2p)+\rho\left(\frac{2}{\pi^{2(p+1)}}(2^{2(p+1)}-1)\zeta\left(2(p+1)\right)+\delta
2^{2p}\right)\right).$

###### Proof.

The limits as $\theta\to 0$ and the values at $\theta=\pi$ are readily
obtained from the explicit expressions of $M_{p},B_{p},D_{p}$ in Proposition
5.1:

$\displaystyle\lim_{\theta\to 0}M_{p}(\theta)=1,\quad\lim_{\theta\to
0}B_{p}(\theta)=0,\quad\lim_{\theta\to 0}D_{p}(\theta)=0;$ (30) $\displaystyle
M_{p}(\pi)=2^{2(p+1)}\sum_{j\in\mathbb{Z}}$
$\displaystyle\frac{1}{(\pi+2j\pi)^{2(p+1)}},\quad
B_{p}(\pi)=-2^{2(p+1)}\sum_{j\in\mathbb{Z}}\frac{1}{(\pi+2j\pi)^{2p}},\quad
D_{p}(\pi)=(-1)^{p}2^{2p}.$ (31)

It remains to explicitly calculate the two infinite sums. We have, for $m\geq
1$,

$\displaystyle\sum_{j\in\mathbb{Z}}\frac{1}{(\pi+2j\pi)^{2m}}$
$\displaystyle=\frac{1}{\pi^{2m}}\sum_{j\in\mathbb{Z}}\frac{1}{(1+2j)^{2m}}=\frac{1}{\pi^{2m}}\sum_{j=0}^{\infty}\frac{1}{(1+2j)^{2m}}+\frac{1}{\pi^{2m}}\sum_{j=1}^{\infty}\frac{1}{(2j-1)^{2m}}$
$\displaystyle=\frac{2}{\pi^{2m}}\sum_{j=1}^{\infty}\frac{1}{(2j-1)^{2m}}=\frac{1}{\pi^{2m}}\frac{(2^{2m}-1)}{2^{2m-1}}\zeta(2m),$

where for the last equality, we used [13, Equation 0.233.5]. The expressions
for $q^{\mathbf{K}^{p}(\rho,\delta)}(\pm 1)$ follow from (29). ∎

######  Remark 5.5.

The cardinal splines satisfy the property of partition for unity (see [4,
Theorem 4.3]), meaning $\sum_{j\in\mathbb{Z}}\Phi_{p}(s-j)=1$, and
consequently $\sum_{j\in\mathbb{Z}}\partial_{s}\Phi_{p}(s-j)=0$. Therefore,
$q^{\mathbf{K}^{p}(\rho,\delta)}(1)$ can be computed directly without
Proposition 5.1 as

$q^{\mathbf{K}^{p}(\rho,\delta)}(1)=\sum_{j\in\mathbb{Z}}\partial_{s}^{2}\Phi_{2p+1}(p+1-j)+\rho\left(\sum_{j\in\mathbb{Z}}\Phi_{2p+1}(p+1-j)+\delta(-1)^{p}\sum_{j\in\mathbb{Z}}\partial_{s}^{2p}\Phi_{2p+1}(p+1-j)\right)=\rho.$

Recalling Proposition 4.7, we establish a condition on the parameters $\rho$
and $\delta$ that is necessary for the weak well-conditioning of the family of
matrices $\\{\mathbf{K}_{n}^{p}(\rho,\delta)\\}_{n}$, for each $p\geq 1$.

######  Proposition 5.6.

Let $p\geq 1$ and define the quantities

$\delta_{p}=-\frac{M_{p}(\pi)}{(-1)^{p}D_{p}(\pi)}=-\frac{(2^{2(p+1)}-1)}{2^{2p-1}\pi^{2(p+1)}}\zeta\left(2(p+1)\right),\quad\rho_{p}=-\frac{B_{p}(\pi)}{M_{p}(\pi)}=4\pi^{2}\frac{(2^{2p}-1)}{(2^{2(p+1)}-1)}\frac{\zeta(2p)}{\zeta(2(p+1))}$
(32)

with $M_{p},B_{p}$ and $D_{p}$ as in Proposition 5.1. Then the following
statements are valid.

  * •

If the family of matrices $\\{\mathbf{K}_{n}^{p}(\rho,\delta)\\}_{n}$ is
weakly well-conditioned for all $\rho>0$, then $\delta\leq\delta_{\rho}$.

  * •

If the family of matrices $\\{\mathbf{K}^{p}_{n}(\rho,0)\\}_{n}$ is weakly
well-conditioned, then $\rho\leq\rho_{p}$.

###### Proof.

The result directly follows by Proposition 4.7, and Corollaries 5.4 and 4.2.
Specifically, if the family of matrices
$\\{\mathbf{K}^{p}_{n}(\rho,\delta)\\}_{n}$ is weakly well-conditioned, then
we have
$(-1)^{p}q^{\mathbf{K}^{p}(\rho,\delta)}(-1)q^{\mathbf{K}^{p}(\rho,\delta)}(1)\leq
0$, which is equivalent to

$-B_{p}(\pi)+\rho(M_{p}(\pi)+\delta(-1)^{p}D_{p}(\pi))\leq 0$

or explicitly,

$-\frac{2^{3}}{\pi^{2p}}(2^{2p}-1)\zeta(2p)+\rho\left(\frac{2}{\pi^{2(p+1)}}(2^{2(p+1)}-1)\zeta\left(2(p+1)\right)+\delta
2^{2p}\right)\leq 0.\qed$

######  Remark 5.7.

Note that $\zeta(m)\to 1$ as $m\to+\infty$. This can observed from the
inequality

$1\leq\zeta(m)=\sum_{n=1}^{\infty}n^{-m}=1+\sum_{n=2}^{\infty}n^{-m}\leq
1+\int_{1}^{\infty}x^{-m}\differential x=1+\frac{1}{m-1}.$

Consequently, we deduce $\rho_{p}\to\pi^{2}$ and
$\delta_{p}\sim-8/\pi^{2}\,\pi^{-2p}\to 0$ as $p\to\infty$. These results are
consistent with those obtained in [9], where the trend
$\delta_{p}\approx-10^{-p}$ was observed. Furthermore, the constants
$\rho_{1}$ and $\delta_{1}$ are exactly the same as those obtained in [25].
Some explicit values of $\rho_{p}$ and $\delta_{p}$ are reported in Table 1.

$p$ | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
$\rho_{p}$ | 12 | 10 | $\nicefrac{{168}}{{17}}$ | $\nicefrac{{306}}{{31}}$ | $\nicefrac{{2349}}{{238}}$ | $\nicefrac{{7797}}{{790}}$
$\delta_{p}$ | $-\nicefrac{{1}}{{12}}$ | $-\nicefrac{{1}}{{120}}$ | $-\nicefrac{{17}}{{20160}}$ | $-\nicefrac{{5}}{{58529}}$ | $-\nicefrac{{2}}{{231067}}$ | $-\nicefrac{{1}}{{1140271}}$
Table 1: Values of $\rho_{p}$ and $\delta_{p}$ as defined in (32).

Finally, for all $p\geq 1$, we have $\rho_{p}>\pi^{2}$. Indeed, this is
equivalent to the inequality

$\frac{\zeta(2p+2)}{\zeta(2p)}<4\frac{2^{2p}-1}{2^{2(p+1)}-1},$ (33)

which is always true as shown in [23, Theorem 1.1] (see also [23, Equation
2.1]).

We finally establish that the conditions stated in Proposition 5.6 are also
sufficient for ensuring the weak well-conditioning of the family
$\\{\mathbf{K}^{p}_{n}(\rho,\delta)\\}_{n}$. This fact relies primarily on
technical estimates established in [5, 6], within the framework of symbols-
based spectral analysis for IgA discretizations.

In the following lemma, we investigate when $q^{\mathbf{K}^{p}(\rho,\delta)}$
has exactly two zeros on the unitary disk.

######  Lemma 5.8.

For $p\geq 1$, let $\delta_{p}$ and $\rho_{p}$ be defined as in (32). Define
the function

$F_{p}(\theta;\rho,\delta)=e^{-\mathrm{i}p\theta}q^{\mathbf{K}^{p}(\rho,\delta)}(e^{\mathrm{i}\theta})=B_{p}(\theta)+\rho(M_{p}(\theta)+\delta(-1)^{p}D_{p}(\theta))$

with $M_{p},B_{p}$ and $D_{p}$ as in Proposition 5.1. Then, the following
statements hold true.

  * •

The equation (in terms of $\theta$) $F_{p}(\theta;\rho,\delta)=0$ has exactly
two simple zeros in $[-\pi,\pi]$ for all $\rho>0$ if and only if
$\delta\leq\delta_{p}$;

  * •

The equation (in terms of $\theta$) $F_{p}(\theta;\rho,0)=0$ has exactly two
simple zeros in $[-\pi,\pi]$ if and only if $\rho\leq\rho_{p}$.

###### Proof.

Let $p\geq 1$. We study the equation only in the interval $[0,\pi]$ since
$F_{p}(\theta;\rho,\delta)$ is even in $\theta$. From (30) and the definition
in Proposition 5.1, we observe that $M_{p}(\theta)>0$ for all
$\theta\in[0,\pi]$. Hence, the ratio $F_{p}(\theta;\rho,\delta)/M_{p}(\theta)$
is well defined for $\theta\in[0,\pi]$. Our objective is to show that for all
$p\geq 1$, $\rho>0$ and $\delta\leq 0$, the function
$F_{p}(\theta;\rho,\delta)/M_{p}(\theta)$ is strictly decreasing in $\theta$
for $\theta\in(0,\pi)$. We begin by computing

$\frac{d}{d\theta}\frac{F_{p}(\theta;\rho,\delta)}{M_{p}(\theta)}=\frac{d}{d\theta}\frac{B_{p}(\theta)}{M_{p}(\theta)}+\rho\delta(-1)^{p}\left(\frac{D^{\prime}_{p}(\theta)M_{p}(\theta)-D_{p}(\theta)M_{p}^{\prime}(\theta)}{(M_{p}(\theta))^{2}}\right).$

It is shown in [6, Theorem 2] that
$\frac{d}{d\theta}\left[B_{p}(\theta)/M_{p}(\theta)\right]<0$ for all $p\geq
1$ and $\theta\in(0,\pi)$. Additionally, we calculate
$(-1)^{p}D_{p}^{\prime}(\theta)=2p\sin\theta(2-2\cos\theta)^{p-1}>0$ for all
$\theta\in(0,\pi)$. Therefore, combining [5, Lemma A.2.] that tell us that
$M_{p}^{\prime}(\theta)<0$ for all $p\geq 1$ and $\theta\in(0,\pi)$, the fact
that $M_{p}(\theta)>0$, and $(-1)^{p}D_{p}(\theta)>0$, we conclude that for
all $\rho>0$ and $\delta<0$,
$\frac{d}{d\theta}[F_{p}(\theta;\rho,\delta)/M_{p}(\theta)]<0$ for all
$\theta\in(0,\pi)$. This confirms that
$F_{p}(\theta;\rho,\delta)/M_{p}(\theta)$ is strictly decreasing in $\theta$
for $\rho>0$ and $\delta\leq 0$.

We can state that $F_{p}(\theta;\rho,\delta)$ has exactly one zero in
$[0,\pi]$ if and only if $F_{p}(0;\rho,\delta)\geq 0$ and
$F_{p}(\pi;\rho,\delta)\leq 0$. Both quantities have been computed in
Corollary 5.4, leading to the desired result. ∎

We conclude this section with the main result of the paper.

######  Theorem 5.9.

For $p\geq 1$, let $\delta_{p}$ and $\rho_{p}$ be the quantities defined in
(32). The following statements hold true.

  * •

The family of matrices $\\{\mathbf{K}_{n}^{p}(\rho,\delta)\\}_{n}$ is weakly
well-conditioned for all $\rho>0$ if and only $\delta\leq\delta_{p}$.

  * •

The family $\\{\mathbf{K}^{p}_{n}(\rho,0)\\}_{n}$ is weakly well-conditioned
if and only if $\rho\leq\rho_{p}$.

###### Proof.

The result follows from Lemmas 5.8 and 4.4. ∎

##  6 Alternative IgA stabilizations

In this section, we briefly analyze other possible stabilization terms, and
show that, in a certain sense, the one proposed in [9] turns out to be the
most effective.

Let us consider, for $k\geq 0$ and a suitable $\delta\in\mathbb{R}$, the
perturbed bilinear form $a_{\mu,h}^{p,k}:S_{0,\bullet,h}^{p}\times
S_{\bullet,0,h}^{p}\to\mathbb{R}$ defined as

$a^{p,k}_{\mu,h}(u_{h},v_{h})=-(\partial_{t}u_{h},\partial_{t}v_{h})_{L^{2}(0,T)}+\mu(u_{h},v_{h})_{L^{2}(0,T)}+\mu\delta
h^{2k}(\partial_{t}^{k}u_{h},\partial_{t}^{k}v_{h})_{L^{2}(0,T)}.$ (34)

Clearly, when $k=p$, this coincides with the stabilized bilinear form
$a_{\mu,h}^{p}$ defined in (5), and for $k\geq p+1$ the stabilizing term is
null, reducing to the unstabilized bilinear form as in (4). Moreover, since we
are ultimately employing FEM, the consistency term should be at least of order
$p+1$ to ensure optimal convergence properties in both the energy and
$L^{2}$-norms, that is $2k\geq p+1$.

######  Remark 6.1.

Another possibility, mantaining the symmetric structure of the numerical
scheme, could have been to include the more general stabilizing term:

$\mu\delta
h^{k_{1}+k_{2}}\left((\partial_{t}^{k_{1}}u_{h},\partial_{t}^{k_{2}}v_{h})_{L^{2}(0,T)}+(\partial_{t}^{k_{2}}u_{h},\partial_{t}^{k_{1}}v_{h})_{L^{2}(0,T)}\right),$

for some $k_{1},k_{2}\geq 0$. However, note that if $k_{1}+k_{2}\equiv 0\
(\bmod\,2)$, then after integration by parts, this term coincides with the one
added in (34) with $k=(k_{1}+k_{2})/2$. Conversely, if $k_{1}+k_{2}\equiv 1\
(\bmod\,2)$, recalling (12), this term vanishes when computed on central
splines.

Again for $\rho=\mu h^{2}$ and $n=N+p-1$, the associated linear system reads:

$\mathbf{K}_{\mu,h}^{p,k}(\delta)=-\mathbf{B}_{h}^{p}+\mu\mathbf{M}_{h}^{p}+\mu\delta
h^{2k}\mathbf{D}_{h}^{p,k},\quad\mathbf{K}_{n}^{p,k}(\rho,\delta)=h\mathbf{K}_{\mu,h}^{p,k}(\delta),$

with $\mathbf{M}_{h}^{p}$ and $\mathbf{B}_{h}^{p}$ as in (8), and the new
matrices

$\mathbf{D}^{p,k}_{h}[\ell,j]=\int_{0}^{T}\partial_{t}^{k}\varphi_{j}^{p}(t)\partial_{t}^{k}\varphi^{p}_{\ell-1}(t)\differential
t$

for $j,\ell=1,\ldots N+p-1$. Similarly as in Proposition 3.2, this new matrix
has the same structure as $\mathbf{D}_{h}^{p}$, with central non-vanishing
elements expressed as

$h^{2k-1}\mathbf{D}^{p,k}_{h}[\ell,\ell-1\pm
j]=(-1)^{k}\partial_{s}^{2k}\Phi_{2p+1}(p+1-j)$

for $j=0,\ldots,p$ and $\ell=2p+1,\ldots,N-p$. Note also that

$\mathbf{D}_{h}^{p,0}=\mathbf{M}^{p}_{h},\quad\mathbf{D}_{h}^{p,1}=\mathbf{B}^{p}_{h}\quad\text{and}\quad\mathbf{D}_{h}^{p,p}=\mathbf{D}_{h}^{p}.$

Proposition 5.1 and Corollary 5.4 can be readily generalized for $0\leq k\leq
p$.

######  Proposition 6.2.

For $0\leq k\leq p$ and $\theta\in[-\pi,\pi]$, the following relation is valid

$\partial_{s}^{2k}\Phi_{2p+1}(p+1)+\sum_{j=1}^{p}\partial_{s}^{2k}\Phi_{2p+1}(p+1-j)(e^{-\mathrm{i}j\theta}+e^{\mathrm{i}j\theta})=D^{k}_{p}(\theta),$

with $D_{p}^{k}:[-\pi,\pi]\to\mathbb{R}$ defined as
$D_{p}^{k}(\theta)=(-1)^{k}(2-2\cos\theta)^{k}M_{p-k}(\theta)$ with $M_{p}$ as
in Proposition 5.1. More explicitly,

$\displaystyle
D_{p}^{k}(\theta)=(-1)^{k}(2-2\cos\theta)^{p+1}\sum_{j\in\mathbb{Z}}\frac{1}{(\theta+2j\pi)^{2(p+1-k)}}.$

In particular, we can evaluate $D_{p}^{0}(0)=1$ and $D_{p}^{k}(0)=0$ for
$k\geq 1$, and

$D_{p}^{k}(\pi)=(-1)^{k}2^{2(k+1)}M_{p-k}(\pi)=(-1)^{k}2^{2k+1}\frac{2^{2(p+1-k)}-1}{\pi^{2(p+1-k)}}\zeta\left(2(p+1-k)\right),$
(35)

for $0\leq k\leq p$.

######  Remark 6.3.

Given that $\zeta(2)=\pi^{2}/6$, the value $D_{p}^{p}(\pi)$ is coherent with
the previously obtained $D_{p}(\pi)=(-1)^{p}2^{2p}$.

The polynomial associated with the family of matrices
$\\{\mathbf{K}_{n}^{p,k}(\rho,\delta)\\}_{n}$ is expressed for
$\theta\in[-\pi,\pi]$ as

$q^{\mathbf{K}^{p,k}(\rho,\delta)}(e^{\mathrm{i}\theta})=e^{\mathrm{i}p\theta}\left(B_{p}(\theta)+\rho\left(M_{p}(\theta)+\delta(-1)^{k}D_{p}^{k}(\theta)\right)\right).$

Generalizing Theorem 5.9, we obtain the main result of this section.

######  Theorem 6.4.

For $p\geq 1$ and $1\leq k\leq p$, let us define

$\delta_{p}^{k}=-\frac{M_{p}(\pi)}{(-1)^{k}D_{p}^{k}(\pi)}=-\frac{1}{2^{2k}\pi^{2k}}\frac{(2^{2(p+1)}-1)}{(2^{2(p+1-k)}-1)}\frac{\zeta\left(2(p+1)\right)}{\zeta\left(2(p+1-k)\right)}.$
(36)

Then, the family of matrices $\\{\mathbf{K}_{n}^{p,k}(\rho,\delta)\\}_{n}$ is
weakly well-conditioned for all $\rho>0$ if and only
$\delta\leq\delta_{p}^{k}$.

###### Proof.

Let us consider the function
$F_{p}^{k}(\theta;\rho,\delta)=e^{-\mathrm{i}p\theta}q^{\mathbf{K}^{p,k}(\rho,\delta)}(e^{\mathrm{i}\theta})$.
We aim to show that the equation (in terms of $\theta$)
$F^{k}_{p}(\theta;\rho,\delta)=0$ has exactly one simple zero in $[0,\pi]$ for
all $\rho>0$ if and only if $\delta\leq\delta_{p}^{k}$. Establishing this fact
leads to the conclusion by Lemma 4.4.

As in the proof of Lemma 5.8, for all $p\geq 1$, $\rho>0$, and $\delta\leq 0$,
we show that the function $F_{p}^{k}(\theta;\rho,\delta)/M_{p}(\theta)$ is
strictly decreasing in $\theta$. Explicitly, we can express it as

$\frac{F_{p}^{k}(\theta;\rho,\delta)}{M_{p}(\theta)}=\frac{B_{p}(\theta)}{M_{p}(\theta)}+\rho+\rho\delta(-1)^{k}\frac{D_{p}^{k}(\theta)}{M_{p}(\theta)}.$

By virtue of [6, Theorem 2], showing that the last term is a decreasing
functions suffices. From Proposition 6.2,

$(-1)^{k}\frac{D^{k}_{p}(\theta)}{M_{p}(\theta)}=(2-2\cos\theta)^{k}\frac{M_{p-k}(\theta)}{M_{p}(\theta)}=\prod_{j=0}^{k-1}(2-2\cos\theta)\frac{M_{p-k+j}(\theta)}{M_{p-k+j+1}(\theta)}.$

Each factor here is a decreasing function (as established in [5, Lemma A.2.])
and positive, hence the product is decreasing as well. Consequently,
$F_{p}^{k}(\theta;\rho,\delta)$ has only one simple zero in $[0,\pi]$ if and
only if $F_{p}^{k}(0;\rho,\delta)\geq 0$ and $F_{p}^{k}(\pi;\rho,\delta)\leq
0$. This holds true for all $\rho>0$ if and only if
$\delta\leq\delta_{p}^{k}$, given that $F_{p}^{k}(0;\rho,\delta)=\rho$. The
explicit expression of $\delta_{p}^{k}$ is a consequence of (35). ∎

######  Remark 6.5.

We note two distinct trends of $\delta_{p}^{k}$ as $p\to\infty$ depending on
the nature of $k$:

  * •

for $k\equiv k(p)$ such that $p+1-k(p)\to\infty$ as $p\to\infty$, then
$\delta_{p}^{k(p)}\to-\pi^{-2k(p)}$;

  * •

for $k\equiv k(p)$ such that $p+1-k(p)\to M\in\mathbb{N}$ as $p\to\infty$,
then $\delta_{p}^{k(p)}\sim-C_{M}\pi^{-2p}\to 0$ , as $p\to\infty$ with

$C_{M}=\frac{1}{\pi^{(2-2M)}(1-2^{-2M})\zeta(2M)}.$

If $k(p)=p$, as in the stabilization of [9], i.e., $M=1$, we obtain the
already mentioned constant $C_{1}=8/\pi^{2}$. Moreover, recalling (33) we
conclude that $\\{C_{M}\\}_{M}$ is an increasing sequence, with
$\nicefrac{{C_{M+1}}}{{C_{M}}}>\pi^{2}$, implying that $C_{M}>8\pi^{2M-4}$.

We conclude by highlighting the features of the stabilized bilinear form
$a_{h,\mu}^{p,p}=a_{h,\mu}^{p}$ (see (5)):

  * •

due to the presence of the parameter $\mu$ only as a monomial at the first
power, this bilinear form can be easily adapted as a space–time bilinear form
$a_{h}^{p}$ (see (6);

  * •

among the possible stabilized bilinear forms $a_{h,\mu}^{p,k}$ (see (34)) with
$0\leq k\leq p$, $a_{h,\mu}^{p}$ guarantees the optimal consistency error in
terms of the power of $h$;

  * •

by considering $a_{h,\mu}^{p,k}$ with $0\leq k\leq p$, and $k(p)$ a function
of $p$ such that $k(p)\to p+1-M$ as $p\to\infty$ for some $M\in\mathbb{N}$,
then the stabilized bilinear form $a_{h,\mu}^{p}$ guarantees the minimum
magnitude for the optimal stabilizing parameter $\delta_{p}^{k}\sim-
C_{M}\pi^{-2p}$ in terms of the constant $C_{M}$.

##  7 Numerical results

In this section, we provide two numerical tests to validate the sharpness of
the quantities $\rho_{p}$ and $\delta_{p}$ defined in (32).

In Figures 4-4 we show the condition numbers of the matrices
$\mathbf{K}_{n}^{p}(\rho,0)$ defined in (14) with $n=1000$ and
$p\in\\{1,2,3,4,5\\}$ by varying the parameter $\rho$. With the vertical lines
we also show the expected values of $\rho_{p}$ as detailed in Table 1.

Figure 3: Spectral condition numbers $\kappa_{2}(\mathbf{K}^{p}_{n}(\rho,0))$
in semi-logarithmic scale, with $n=1000$ by varying $\rho\in[8,13]$, with
$p\in\\{1,2,3\\}$.

Figure 4: Spectral condition numbers $\kappa_{2}(\mathbf{K}^{p}_{n}(\rho,0))$
in semi-logarithmic scale, with $n=2000$ by varying $\rho\in[9.84,9.92]$, with
$p\in\\{4,5\\}$.

In Figure 5, we show the condition numbers of the matrices
$\mathbf{K}_{n}^{p}(\rho,-|\delta|)$ again with $n=1000$,
$p\in\\{1,2,3,4,5,6\\}$, and $\rho=20000$ by varying the absolute value
$|\delta|$. We also show with vertical lines the expected values
$|\delta_{p}|$ as detailed in Table 1.

Figure 5: Spectral condition numbers
$\kappa_{2}(\mathbf{K}^{p}_{n}(\rho,-|\delta|))$ in logarithmic scale, with
$n=1000$ and $\rho=20000$, by varying $|\delta|\in[10^{-8},0.3]$ and
$p\in\\{1,2,3,4,5,6\\}$.

In all cases, Theorem 5.6 results to be sharp.

##  8 Conclusion

In this paper, we have examined the conditioning behaviour of a family of
matrices related to conforming space–time discretizations of the wave
equation, which uses splines of maximum regularity in time. We have provided a
theoretical explanation of the results presented in [9], and obtained precise
estimates on two crucial parameters: the CFL condition for the unstabilized
space–time bilinear form, and the optimal penalty parameter for the stabilized
method. Our analysis is based on results from numerical linear algebra and
properties of symbols associated with spline discretizations. Finally, we have
shown how the stabilization term proposed in [9] is optimal from several
points of view, and that the presented results are sharp.

##  9 Acknowledegments

This research was supported by the Austrian Science Fund (FWF) project
10.55776/F65 and and the Austrian Science Fund (FWF) project 10.55776/P33477.

## References

  * [1] P. Amodio and L. Brugnano.  The conditioning of Toeplitz band matrices.  Math. Comput. Modelling, 23(10):29–42, 1996.
  * [2] P. Bignardi and A. Moiola.  A space-time continuous and coercive formulation for the wave equation.  arXiv 2312.07268, 2023.
  * [3] A. Bressan and E. Sande.  Approximation in FEM, DG and IGA: a theoretical comparison.  Numer. Math., 143(4):923–942, 2019.
  * [4] C. K. Chui.  An introduction to wavelets, volume 1 of Wavelet Analysis and its Applications.  Academic Press, Inc., Boston, MA, 1992.
  * [5] M. Donatelli, C. Garoni, C. Manni, S. Serra-Capizzano, and H. Speleers.  Symbol-based multigrid methods for Galerkin B-spline isogeometric analysis.  SIAM J. Numer. Anal., 55(1):31–62, 2017.
  * [6] S.-E. Ekström, I. Furci, C. Garoni, C. Manni, S. Serra-Capizzano, and H. Speleers.  Are the eigenvalues of the B-spline isogeometric analysis approximation of $-\Delta u=\lambda u$ known in almost closed form?  Numer. Linear Algebra Appl., 25(5):e2198, 34, 2018.
  * [7] J. A. Evans, Y. Bazilevs, I. Babuška, and T. J. R. Hughes.  $n$-widths, sup-infs, and optimality ratios for the $k$-version of the isogeometric finite element method.  Comput. Methods Appl. Mech. Engrg., 198(21–26):1726–1741, 2009\.
  * [8] S. Fraschini.  Stability of space–time isogeometric methods for wave propagation problems.  arXiv 2030.15460, 2023.
  * [9] S. Fraschini, G. Loli, A. Moiola, and G. Sangalli.  An unconditionally stable space–time isogeometric method for the acoustic wave equation.  arXiv 2303.07268, 2023.
  * [10] S. Fraschini, A. Moiola, and G. Sangalli.  Stability of space–time isogeometric methods for wave propagation problems.  In WAVES 2022 - The 15th International Conference on Mathematical and Numerical Aspects of Wave Propagation, pages 286–287, 2022\.
  * [11] C. Garoni, C. Manni, F. Pelosi, S. Serra-Capizzano, and H. Speleers.  On the spectrum of stiffness matrices arising from isogeometric analysis.  Numer. Math., 127(4):751–799, 2014.
  * [12] C. Garoni, H. Speleers, S.-E. Ekström, A. Reali, S. Serra-Capizzano, and T. J. R. Hughes.  Symbol-based analysis of finite element and isogeometric B-spline discretizations of eigenvalue problems: exposition and review.  Arch. Comput. Methods Eng., 26(5):1639–1690, 2019.
  * [13] I. S. Gradshteyn and I. M. Ryzhik.  Table of integrals, series, and products.  Elsevier/Academic Press, Amsterdam, eighth edition, 2015.  Translated from the Russian, Translation edited and with a preface by Daniel Zwillinger and Victor Moll.
  * [14] T. J. R. Hughes, J. A. Cottrell, and Y. Bazilevs.  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.  Comput. Methods Appl. Mech. Engrg., 194(39-41):4135–4195, 2005\.
  * [15] T. J. R. Hughes and G. M. Hulbert.  Space–time finite element methods for elastodynamics: formulations and error estimates.  Comput. Methods Appl. Mech. Engrg., 66(3):339–363, 1988.
  * [16] T. J. R. Hughes, A. Reali, and G. Sangalli.  Duality and unified analysis of discrete approximations in structural dynamics and wave propagation: comparison of $p$-method finite elements with $k$-method NURBS.  Comput. Methods Appl. Mech. Engrg., 197(49–50):4104–4124, 2008\.
  * [17] L.-M. Imbert-Gérard, A. Moiola, and P. Stocker.  A space-time quasi-Trefftz DG method for the wave equation with piecewise-smooth coefficients.  Math. Comp., 92(341):1211–1249, 2023.
  * [18] P. Lakatos.  On zeros of reciprocal polynomials.  Publicationes mathematicae, 61:645–661, 11 2002.
  * [19] R. Löscher, O. Steinbach, and M. Zank.  Numerical results for an unconditionally stable space-time finite element method for the wave equation.  In Domain Decomposition Methods in Science and Engineering XXVI, pages 625–632. Springer, 2023.
  * [20] J. C. Mason and D. C. Handscomb.  Chebyshev polynomials.  Chapman & Hall/CRC, Boca Raton, FL, 2003.
  * [21] A. Moiola and I. Perugia.  A space-time Trefftz discontinuous Galerkin method for the acoustic wave equation in first-order formulation.  Numer. Math., 138(2):389–435, 2018.
  * [22] P. Monk and G. R. Richter.  A discontinuous Galerkin method for linear symmetric hyperbolic systems in inhomogeneous media.  J. Sci. Comput., 22/23:443–477, 2005.
  * [23] F. Qi.  A double inequality for the ratio of two non-zero neighbouring Bernoulli numbers.  J. Comput. Appl. Math., 351:1–5, 2019.
  * [24] T. J. Rivlin.  Chebyshev polynomials.  Pure and Applied Mathematics (New York). John Wiley & Sons, Inc., New York, second edition, 1990.  From approximation theory to algebra and number theory.
  * [25] O. Steinbach and M. Zank.  A stabilized space–time finite element method for the wave equation.  In: Advanced Finite Element Methods with Applications. Selected papers from the 30th Chemnitz FEM Symposium 2017, (T. Apel, U. Langer, A. Meyer, O. Steinbach eds., Lecture Notes in Computational Science and Engineering, pages 315–342, 2019.
  * [26] O. Steinbach and M. Zank.  Coercive space–time finite element methods for initial boundary value problems.  Electron. Trans. Numer. Anal., 52:154–194, 2020.
  * [27] R. Vieira.  Polynomials with symmetric zeros.  In Cheon Seoung Ryoo, editor, Polynomials, chapter 4. IntechOpen, Rijeka, 2019.
  * [28] M. Zank.  Inf-sup stable space–time methods for time-dependent partial differential equations.  Verlag d. Technischen Universität Graz, 2020.
  * [29] M. Zank.  Higher-order space–time continuous galerkin methods for the wave equation.  In 14th WCCM-ECCOMAS Congress 2020, volume 700, 2021.