# Outlier-free isogeometric discretizations for Laplace eigenvalue problems: closed-form eigenvalue and eigenvector expressions

Noureddine Lamsahel, Carla Manni, Ahmed Ratnani, Stefano Serra-Capizzano, Hendrik Speleers

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

## ► To cite this version:

Noureddine Lamsahel, Carla Manni, Ahmed Ratnani, Stefano Serra-Capizzano, Hendrik Speleers. Outlier-free isogeometric discretizations for Laplace eigenvalue problems: closed-form eigenvalue and eigenvector expressions. Numerische Mathematik, 2025, 10.1007/s00211-025-01470-x. http://dx.doi.org/10.1007/s00211-025-01470-x

HAL Id: hal-05127663

https://hal.science/hal-05127663v1

Submitted on 24 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Outlier-free isogeometric discretizations for Laplace eigenvalue problems: closed-form eigenvalue and eigenvector expressions

Noureddine Lamsahel¹,² · Carla Manni³ · Ahmed Ratnani¹ · Stefano Serra-Capizzano⁴,⁵ · Hendrik Speleers³

Received: 26 April 2024 / Revised: 2 March 2025 / Accepted: 20 April 2025
© The Author(s) 2025

## Abstract

We derive explicit closed-form expressions for the eigenvalues and eigenvectors of the matrices resulting from isogeometric Galerkin discretizations based on outlier-free spline subspaces for the Laplace operator, under different types of homogeneous boundary conditions on bounded intervals. For optimal spline subspaces and specific reduced spline spaces, represented in terms of B-spline-like bases, we show that the corresponding mass and stiffness matrices exhibit a Toeplitz-minus-Hankel or Toeplitz-plus-Hankel structure. Such matrix structure holds for any degree *p* and implies that the eigenvalues are an explicitly known sampling of the spectral symbol of the Toeplitz part. Moreover, by employing tensor-product arguments, we extend the closed-form property of the eigenvalues and eigenvectors to a *d*-dimensional box. As a side result, we have an algebraic confirmation that the considered optimal and reduced spline spaces are indeed outlier-free.

**Mathematics Subject Classification** 34L16 · 41A15 · 65D07 · 65F15 · 65L15 · 65L60

## 1 Introduction

Isogeometric analysis (IGA) is a well-established paradigm for the numerical treatment of differential problems, introduced in [29], with the aim of unifying computer aided design (CAD) and finite element analysis (FEA). The leading idea in IGA is to use the same representation tools for the design as well as for the analysis (in an isoparametric environment), providing a true design-through-analysis methodology [8, 29]. In its original form, IGA adopts the core functions of CAD systems, i.e., B-splines and their rational counterpart (NURBS), to approximate the solution of the considered differential problems. The isogeometric approach shows important advantages over classical $C^0$ FEA. In particular, the inherently high smoothness of B-splines and

<footnote>Extended author information available on the last page of the article</footnote>
NURBS allows for a higher accuracy per degree of freedom; see [6, 34, 35] and references therein.

Spectral analysis typically provides a global picture of the eigenvalues and eigen-
functions of a given (discretized) differential operator. The study of the spectral error is
of major interest in applications because, in practice, computations are not performed
in the limit of mesh refinement, and for a large class of boundary and initial-value
problems the total discretization error on a given mesh can be recovered from its
spectral error; see, e.g., [30]. Primarily second order eigenvalue problems have been
addressed and related IGA discretizations have been investigated by several authors;
see [9, 27, 32] and references therein.

The explicit error estimates given in [34, 35] for spline spaces of different smooth-
ness show that discretizations in spline spaces of maximal smoothness on uniform
grids are preferable. This theoretical finding totally agrees with the numerical obser-
vations: maximally smooth spline discretizations of degree p on uniform grids lead
to a prime discrete spectrum, which almost entirely converges to the continuous one
for increasing p and does not present the well-known artifacts of high-order C⁰ FEA,
the so-called optical branches. From a spectral and matrix theoretical point of view
[2, 26, 27], it is clearly understood that the number of optical branches amounts to
(p − k)d − 1, with k being the imposed smoothness and d the dimensionality of the
differential problem. The previous formula indicates that these spectral artifacts dis-
appear only in the case of maximal smoothness, i.e., k = p − 1, while in all other
cases their number is exponential in the dimensionality d.

Notwithstanding maximally smooth spline spaces on uniform grids are an excellent
choice for addressing eigenvalue problems, they are not completely satisfactory: a
very small portion of the eigenvalues are poorly approximated and the corresponding
computed values are much larger than the exact ones. These spurious values are usually
referred to as outliers [9]. The number of outliers increases with the degree *p* and, for
fixed *p*, it is independent of the degrees of freedom for univariate problems, while a
“thin layer” of outliers appears in the multivariate setting. Discretization techniques
that remove those outliers are often called outlier-free discretizations. Accurate outlier-
free discretizations are important both for their superior description of the spectrum of
the continuous operator and for their beneficial effects in various applicative contexts
such as (explicit) dynamics.

A solid analytical answer to the problem of outlier removal from isogeometric Galerkin discretizations of second order eigenvalue problems has been recently given in [32] by considering subspaces of maximally smooth splines of any degree that are optimal spaces in the sense of $L^2$ Kolmogorov $n$-widths for certain function classes related to the considered eigenvalue problems [19, 20]. These optimal spline subspaces are simply identified by additional homogeneous boundary conditions and by knots that are uniform (or shifted uniform) and are equipped with a B-spline-like basis. For any degree, they give rise to a complete removal of the outliers without loss of accuracy in the approximation of all eigenfunctions.

The fact that outliers are related to the treatment of boundary conditions has been
observed already for a long time. Interesting empirical procedures for outlier removal
based on this intuition have been presented in [12, 28], where the authors exploit the
imposition of suitable homogeneous additional boundary conditions in a similar (but
slightly different) manner to the above mentioned optimal spline subspaces. More precisely, in [12] a penalization of specific high-order derivatives near the boundary is proposed, while in [28] the same high-order derivatives at the boundary are strongly set equal to zero, by using suitable reduced spline spaces as trial spaces; see also [38, 39]. Although there is a clear numerical evidence that the strong imposition of the considered additional boundary conditions removes the outliers without affecting the accuracy of the approximation for the global spectrum, a proper theoretical foundation of the whole process is missing.

The knowledge of the spectrum of matrices arising from a given discretization technique for partial differential equations, not necessarily confined to eigenvalue differential problems, plays also an important role in designing efficient and robust solvers for the resulting linear systems. It is the foundation for the convergence analysis of (preconditioned) Krylov methods and for the theoretical analysis of effective preconditioners and fast multigrid/multi-iterative solvers.

The spectral behavior of matrices arising from isogeometric Galerkin methods based on standard B-splines and NURBS, with boundary conditions in strong or weak form, has been analyzed in a series of recent papers; see, e.g., [22, 23] and references therein. The analysis typically exploits the theory of generalized locally Toeplitz (GLT) sequences [24, 25] as well as properties of the bases used to represent the discretization spaces. More precisely, it has been shown that sequences of such matrices corresponding to different levels of mesh refinement admit an eigenvalue distribution, which can be compactly described by a function called (GLT) symbol. This symbol has a canonical structure incorporating the approximation technique, the geometry of the domain, and the coefficients of the principal terms of the differential problem. The eigenvalues of the considered matrices can be approximated by a uniform sampling of the symbol and, when the size of the matrices increases, such an approximation converges to the true spectrum, up to few outliers. As application, the spectral information obtained from the symbol has been utilized to design efficient preconditioners and multigrid solvers [17], but the presence of the outliers dirties and weakens the elegant symbol information.

We want to emphasize that none of the above results provides closed-form expressions of the discrete eigenvalues and the corresponding eigenvectors for the considered (full or reduced) spline discretization spaces, unless in few specific cases of low degree [18]. However, it has been observed in [11] that the discretization matrices for the standard full spline spaces are almost Toeplitz-minus-Hankel or Toeplitz-plus-Hankel matrices and that, if they would have exactly this structure, explicit closed-form expressions for the corresponding eigenpairs could be realized. Such structured matrices belong to certain $\tau$ matrix algebras, which have been widely studied in the literature; see, e.g., [3, 5, 14] and references therein.

In this paper, we continue the spectral study of isogeometric Galerkin discretizations in optimal spline subspaces for the Laplace operator, under different types of boundary conditions, and we assume that the considered subspaces are represented in terms of the B-spline-like basis from [32] (see also [15]). We derive explicit closed-form expressions for the eigenvalues and eigenvectors of the corresponding mass and stiffness matrices, for any degree. These expressions show that the eigenvalues are a simple sampling of the GLT symbol [21], with no outliers. Moreover, when fixing the
type of boundary conditions, both the mass and the stiffness matrix exhibit the same Toeplitz-minus-Hankel or Toeplitz-plus-Hankel structure, and thus they belong to the same $\tau$ matrix algebra [5]. This allows us to deduce also explicit closed-form expressions for the approximated spectrum of the Laplace operator, under the addressed boundary conditions. In this way, we get an algebraic confirmation that optimal spline subspaces lead to outlier-free discretizations.

Furthermore, we extend our analysis to isogeometric Galerkin discretizations based on the reduced spline spaces considered in [28] (see also [38, 39]) and represented in terms of the B-spline-like basis from [32]. It turns out that the corresponding stiffness and mass matrices exhibit again a Toeplitz-minus/plus-Hankel structure. This gives rise to explicit closed-form expressions for their eigenvalues and eigenvectors. As a consequence, we can conclude that such spaces are outlier-free as well.

## 1.1 Main results

The eigenvalues of a square matrix A of size n are denoted by $\lambda_j(A)$, $j = 1, \dots, n$, and the corresponding eigenvectors by $u_j(A)$ with components $u_{i,j}(A)$, $i, j = 1, \dots, n$. For $\theta \in [0, \pi]$, $r \in \{0, 1\}$, and $p \ge 1$, we define

$$g_p^r(\theta) = (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1) + 2(-1)^r \sum_{k=1}^{p} \mathcal{N}_{2p+1}^{(2r)}(p+1-k) \cos(k\theta),$$

where $\mathcal{N}_p$ is the cardinal B-spline of degree $p \in \mathbb{N}$, given by (16) and (17). Then, our main results can be summarized as follows.

1. We derive explicit closed-form expressions for the eigenvalues and eigenvectors of the mass matrix $M_n^{p,b}$ and the stiffness matrix $K_n^{p,b}$ when considering the B-spline-like basis of the univariate optimal spline subspaces from [32] of dimension $n$, for any degree $p \ge 1$ and three different types of boundary conditions, indicated with the index $b = 0, 1, 2$.

* Dirichlet boundary conditions (2), indicated with $b = 0$: The matrices $M_n^{p,0}$ and $K_n^{p,0}$ are symmetric and centrosymmetric. Under the assumption $n \ge \max\{p+1, p+\lfloor \frac{p}{2} \rfloor - 1\}$, Theorem 11 and Corollary 13 reveal that for $j = 1, \dots, n$,

$$\lambda_j(M_n^{p,0}) = (n+1)^{-1} g_p^0 \left( \frac{j\pi}{n+1} \right), \quad \lambda_j(K_n^{p,0}) = (n+1) g_p^1 \left( \frac{j\pi}{n+1} \right),$$

and for $i, j = 1, \dots, n$,

$$u_{i,j}(M_n^{p,0}) = u_{i,j}(K_n^{p,0}) = \sqrt{\frac{2}{n+1}} \sin \left( \frac{ij\pi}{n+1} \right).$$

* Neumann boundary conditions (3), indicated with $b = 1$: The matrices $M_n^{p,1}$ and $K_n^{p,1}$ are symmetric and centrosymmetric. Under the assumption $n \ge \max\{p+1, p+\lfloor \frac{p}{2} \rfloor - 1\}$, Theorem 11 and Corollary 13 reveal that for $j = 1, \dots, n$,
max {2p - ⌊p/2⌋, 2p - 2 ⌊p/2⌋ + 1}, Theorem 17 and Corollary 18 reveal that
for j = 1, ..., n,

$$
\lambda_j (M_n^{p,1}) = n^{-1} g_p^0 \left( \frac{(j-1)\pi}{n} \right), \quad \lambda_j (K_n^{p,1}) = n g_p^1 \left( \frac{(j-1)\pi}{n} \right),
$$

and for i, j = 1, ..., n,

$$
\begin{align*}
u_{i,j}(M_n^{p,1}) &= u_{i,j}(K_n^{p,1}) = \sqrt{\frac{2}{n}} c_j \cos\left(\frac{(j-1)\pi}{n}\left(i - \frac{1}{2}\right)\right), \\
c_j &= \begin{cases} \frac{1}{\sqrt{2}}, & j=1, \\ 1, & j \ge 2. \end{cases}
\end{align*}
$$

* Mixed boundary conditions (4), indicated with $b=2$: The matrices $M_n^{p,2}$ and $K_n^{p,2}$ are only symmetric. Under the assumption $n \ge \max\{p+1, p+\lfloor\frac{p}{2}\rfloor\}$, Theorem 23 and Corollary 24 reveal that for $j=1, \dots, n$,

$$
\lambda_j (M_n^{p,2}) = \left(\frac{2n+1}{2}\right)^{-1} g_p^0 \left(\frac{(2j-1)\pi}{2n+1}\right),
$$

$$
\lambda_j (K_n^{p,2}) = \left(\frac{2n+1}{2}\right) g_p^1 \left(\frac{(2j-1)\pi}{2n+1}\right),
$$

and for i, j = 1, ..., n,

$$
u_{i,j}(M_n^{p,2}) = u_{i,j}(K_n^{p,2}) = \sqrt{\frac{4}{2n+1}} \sin \left( \frac{i(2j-1)\pi}{2n+1} \right).
$$

We also identify a new Hankel matrix for which we derive the explicit eigen-
values and eigenvectors in Theorem 21.

2. The above results for optimal spline subspaces are expanded to other outlier-free spline spaces. For simplicity of presentation, we concentrate exclusively on Dirichlet boundary conditions but similar results hold for other boundary conditions as well.

* We derive explicit closed-form expressions for the eigenvalues and eigenvec-
tors of the mass matrix $\overline{M}_n^{p,0}$ and the stiffness matrix $\overline{K}_n^{p,0}$, when considering
the B-spline-like basis of the univariate reduced spline spaces from [32] of
dimension $n$, for any even degree $p \ge 1$ and Dirichlet boundary conditions.
The matrices $\overline{M}_n^{p,0}$ and $\overline{K}_n^{p,0}$ are symmetric and centrosymmetric. Under
the assumption $n \ge p + \frac{p}{2}$, Theorem 28 and Corollary 29 reveal that for
$j = 1, \dots, n$,

$$
\lambda_j (\overline{M}_n^{p,0}) = n^{-1} g_p^0 \left( \frac{j\pi}{n} \right), \quad \lambda_j (\overline{K}_n^{p,0}) = n g_p^1 \left( \frac{j\pi}{n} \right),
$$
and for $i, j = 1, \dots, n$,

$$ 
\begin{aligned}
u_{i,j}(\overline{M}_n^{p,0}) &= u_{i,j}(\overline{K}_n^{p,0}) = \sqrt{\frac{2}{n}} c_j \sin\left(\frac{j\pi}{n}\left(i - \frac{1}{2}\right)\right), \\
c_j &= \begin{cases} \frac{1}{\sqrt{2}}, & j = n, \\ 1, & j \le n - 1. \end{cases}
\end{aligned}
$$

For other degrees $p$ and other boundary conditions, the reduced spline spaces from [28] are either optimal spline subspaces, and hence covered by item 1, or the analysis is similar.

* We extend the closed-form property of the eigenvalues and eigenvectors to the multivariate tensor-product setting on a $d$-dimensional box; see Theorem 33.

In all cases, explicit closed-form expressions for the resulting approximations of the spectrum of the Laplace operator, under the addressed boundary conditions, can be immediately deduced; see Corollaries 14, 19, 25, and 30. Our main results confirm that the optimal spline subspaces from [32] and the reduced spline spaces from [28] are indeed outlier-free; see Remarks 16, 20, 26, and 31. Moreover, it follows that, for degree $p = 3$ (resp., $p = 2$) and under Dirichlet (resp., Neumann) boundary conditions, the penalization technique in [12], when the penalty parameter goes to infinity, exactly corresponds to the first (resp., second) optimal spline subspace discretization. Finally, it is also evident that the GLT symbols of the optimal spline subspaces and the above reduced spline spaces are identical and they are equal to the symbol of the related full spline space of the same degree.

## 1.2 Outline of the paper

In Sect. 2, we provide the necessary notations, definitions, and preliminary results crucial for our analysis. We start by introducing the eigenvalue problems of interest and by providing a brief overview of optimal spline subspaces designed for three classical types of boundary conditions (Dirichlet/Neumann/mixed) along with the construction of the corresponding B-spline-like bases. Afterwards, we collect from the literature some results about the closed form of the eigenvalues and eigenvectors for structured matrices of Toeplitz-minus-Hankel or Toeplitz-plus-Hankel form. In Sect. 3, we delve into the matrices resulting from isogeometric Galerkin discretizations based on optimal spline subspaces. This section is organized into three subsections, each dedicated to a specific type of boundary conditions. Adopting a coherent format across these subsections, we initiate the discussion by introducing a matrix notation designed to handle both mass and stiffness matrices simultaneously. Upon establishing the Toeplitz-minus-Hankel or Toeplitz-plus-Hankel structure of this matrix, we deduce a closed-form expression for its eigenvalues and eigenvectors. For mixed boundary conditions, a new class of Hankel matrices is identified and the expression for the associated eigenvalues and eigenvectors is derived. To enhance readability, the detailed proofs of the main results presented in this section are collected in Appendix A. In Sect. 4, our focus shifts to other outlier-free spline spaces. Firstly, we consider specific
reduced spline spaces, similar to the optimal spline subspaces. Following a parallel approach to that in Sect. 3, we derive closed-form expressions for the eigenvalues and eigenvectors of the corresponding mass and stiffness matrices, which allows us to conclude that these spaces are truly outlier-free. Then, we extend the closed-form expressions to spline discretizations defined on a bounded box of $\mathbb{R}^d$, $d \ge 1$. Finally, Sect. 5 summarizes the results, provides some further comments, and discusses some future work.

## 2 Model problems and preliminaries

In this section, we describe our eigenvalue problems of interest, considering three classical types of boundary conditions and for each of them we define a family of optimal spline subspaces. In [32], it was shown that these optimal spline subspaces lead to outlier-free discretizations and a B-spline-like basis was constructed, which we summarize here. As we will show later in the paper, the corresponding discretization matrices possess a Toeplitz-minus-Hankel or Toeplitz-plus-Hankel structure. Therefore, we present some known results about the closed form of the eigenvalues and eigenvectors of such matrices here as well.

### 2.1 Eigenvalue problems

We consider three one-dimensional eigenvalue problems, each corresponding to the Laplace operator, of the form

$$ \left\{ \begin{array}{ll} -u'' = \lambda u, & \text{in } (0, 1), \\ \text{BC}. & \end{array} \right. \qquad (1) $$

Here, BC serves as an abbreviation for boundary conditions and represents one of the following.

* Dirichlet boundary conditions:

$$ u(0) = u(1) = 0; \qquad (2) $$

* Neumann boundary conditions:

$$ u'(0) = u'(1) = 0; \qquad (3) $$

* Mixed boundary conditions:

$$ u(0) = u'(1) = 0. \qquad (4) $$
It is well established that the families of non-trivial exact solutions ($\lambda_k, u_k$) of (1) subject to the boundary conditions (2)–(4) are given by, respectively,

$$ \left\{ \left( (k\pi)^2, \sin(k\pi \cdot) \right), \quad k = 1, 2, \dots \right\}, \qquad (5) $$

$$ \left\{ \left( ((k-1)\pi)^2, \cos((k-1)\pi \cdot) \right), \quad k = 1, 2, \dots \right\}, \qquad (6) $$

$$ \left\{ \left( \left( \frac{(2k-1)\pi}{2} \right)^2, \sin \left( \frac{(2k-1)\pi}{2} \cdot \right) \right), \quad k = 1, 2, \dots \right\}. \qquad (7) $$

In [18, Section 3], it was illustrated how a closed-form expression for the approximate eigenvalues and eigenvectors can be achieved for Galerkin discretizations in standard full spline spaces related to the Dirichlet boundary conditions (2). However, the proposed approach is restricted to $p = 1, 2$. When $p$ is greater than or equal to 3, the mass and stiffness matrices are not in the same matrix algebra anymore, making the approach unfeasible. More recently, in [12, Section 5], it was analytically shown that a closed-form expression for the approximate eigenvalues and eigenvectors can be derived using another isogeometric discretization method, namely the penalization technique. Similar to the limitation in [18], the exact expression is confined to specific scenarios, namely $p = 3$ for Dirichlet boundary conditions (2) and $p = 2$ for Neumann boundary conditions (3).

One of the main objectives of this paper is to derive an explicit closed-form expression for the approximate eigenvalues and eigenvectors of the pairs in (5)–(7) when an optimal spline subspace (or other types of outlier-free spline spaces) is used as a Galerkin discretization space, for any degree $p \ge 1$.

**Remark 1** Other classical boundary conditions for the eigenvalue problem (1) are periodic boundary conditions. In this case, the discretization based on classical periodic B-splines on uniform grids gives rise to a circulant structure both for the stiffness and the mass matrix, and thus results in a full (explicit) knowledge of their spectrum; see, e.g., [39].

## 2.2 Optimal spline subspaces

We provide a concise overview of the optimal spline discretizations associated with (1) and presented in [32]. We start by defining the optimal spline subspaces, followed by the corresponding variational formulation of the discrete systems. Lastly, we describe the construction of a B-spline-like basis for each subspace in terms of cardinal B-splines. For more details on these spaces, we refer the reader to [32] and references therein.

Let $p, n_{el} \in \mathbb{N}^*$. Following [32] and using the same notation, we consider the vector of strictly increasing breakpoints

$$ \tau = (\tau_0 = 0, \tau_1, \dots, \tau_{n_{el}-1}, \tau_{n_{el}} = 1), \qquad (8) $$
where $n_{el}$ is the number of intervals of $[0, 1]$. We then introduce the standard spline space of maximal smoothness

$$ \mathbb{S}_{p, \tau} = \{ s \in \mathbb{C}^{p-1}([0, 1]) : s|_{[\tau_i, \tau_{i+1})} \in \mathbb{P}_p, i = 0, \dots, n_{el} - 1 \}, $$

where $\mathbb{P}_p$ is the space of polynomials of degree less than or equal to $p$. A comprehensive exposition of the space $\mathbb{S}_{p, \tau}$ and its B-spline basis can be found in [10, 31, 36].

Now, let us fix an integer number $\ell \in \{0, \dots, p\}$. In the following, we will focus on specific subspaces of $\mathbb{S}_{p, \tau}$ designed to selectively capture elements with certain vanishing derivatives at the boundary, precisely up to order $\ell$:

$$ \mathbb{S}_{p, \tau, 0}^{\ell} = \{ s \in \mathbb{S}_{p, \tau} : s^{(\beta)}(0) = s^{(\beta)}(1) = 0, 0 \le \beta \le \ell, \beta \text{ even} \}, \quad (9) $$

$$ \mathbb{S}_{p, \tau, 1}^{\ell} = \{ s \in \mathbb{S}_{p, \tau} : s^{(\beta)}(0) = s^{(\beta)}(1) = 0, 0 \le \beta \le \ell, \beta \text{ odd} \}, \quad (10) $$

$$ \mathbb{S}_{p, \tau, 2}^{\ell} = \{ s \in \mathbb{S}_{p, \tau} : s^{(\beta_0)}(0) = s^{(\beta_1)}(1) = 0, 0 \le \beta_0, \beta_1 \le \ell, \beta_0 \text{ even}, \beta_1 \text{ odd} \}. \quad (11) $$

It can be easily verified that the eigenfunctions in (5)–(7) of the problems (1) satisfy the boundary conditions in the spaces $\mathbb{S}_{p, \tau, b}^{\ell}$ for $b = 0, 1, 2$, respectively. The optimal spline subspaces of dimension $n$ ($n \in \mathbb{N}^*$) associated with the boundary conditions (2)–(4) are obtained by taking $\ell = p$ and specific breakpoints. More precisely,

$$ \mathbb{S}_{p,n,b}^{\text{opt}} = \mathbb{S}_{p, \tau_{p,n,b}^{\text{opt}}, b}^{p}, \quad b = 0, 1, 2, \quad (12) $$

where the breakpoints $\tau_{p,n,b}^{\text{opt}}, b = 0, 1, 2$ are specified as

$$ \tau_{p,n,0}^{\text{opt}} = \begin{cases} \left(0, \frac{1}{n+1}, \frac{2}{n+1}, \dots, \frac{n}{n+1}, 1\right), & p \text{ odd}, \\ \left(0, \frac{1}{n+1}, \frac{3}{n+1}, \dots, \frac{n+1}{n+1}, 1\right), & p \text{ even}, \end{cases} $$

$$ \tau_{p,n,1}^{\text{opt}} = \begin{cases} \left(0, \frac{1}{n}, \frac{3}{n}, \dots, \frac{n-1}{n}, 1\right), & p \text{ odd}, \\ \left(0, \frac{1}{n}, \frac{2}{n}, \dots, \frac{n-1}{n}, 1\right), & p \text{ even}, \end{cases} $$

$$ \tau_{p,n,2}^{\text{opt}} = \begin{cases} \left(0, \frac{2}{2n+1}, \frac{4}{2n+1}, \dots, \frac{2n}{2n+1}, 1\right), & p \text{ odd}, \\ \left(0, \frac{1}{2n+1}, \frac{3}{2n+1}, \dots, \frac{2n-1}{2n+1}, 1\right), & p \text{ even}. \end{cases} $$

In the Galerkin method, the approximate solutions of (1) with (2)–(4) are computed as follows:

$$ \begin{array}{l} \text{Find } (\lambda_{k,h}, u_{k,h}) \in \mathbb{R} \times \mathbb{S}_{p,n,b}^{\text{opt}} \text{ such that} \\ K(u_{k,h}, v) = \lambda_{k,h} M(u_{k,h}, v), \quad \forall v \in \mathbb{S}_{p,n,b}^{\text{opt}}, \end{array} \quad (13) $$
where

$$M(u_{k,h}, v) = \int_{0}^{1} u_{k,h}(x)v(x) \, dx,$$

$$K(u_{k,h}, v) = \int_{0}^{1} u'_{k,h}(x)v'(x) \, dx.$$

For well-posedness of the problem, the convergence of (13), and the outlier-free property of the spaces (12), we refer the reader to [32, Section 4].

Given a fixed $b \in \{0, 1, 2\}$, we have $\dim(\mathbb{S}^{\text{opt}}_{p,n,b}) = n$. Let $B_b^p = \{N_{1,b}^p, N_{2,b}^p, \dots, N_{n,b}^p\}$ denote a basis of the space $\mathbb{S}^{\text{opt}}_{p,n,b}$ (our basis of interest will be specified later in Sects. 2.2.1–2.2.3). The problem described by (13) can be reformulated as a discrete eigenvalue problem,

$$(M_n^{p,b})^{-1} K_n^{p,b} \mathbf{u}_{k,h} = \lambda_{k,h} \mathbf{u}_{k,h}, \quad (14)$$

where $\mathbf{u}_{k,h}$ represents the coefficient vector of $u_{k,h}$ in terms of the basis $B_b^p$. The matrices $M_n^{p,b}$ and $K_n^{p,b}$ are called the mass matrix and the stiffness matrix, respectively, and are defined by

$$M_n^{p,b} = \left[ \int_0^1 N_{i,b}^p(x) N_{j,b}^p(x) \, dx \right]_{i,j=1}^n, \\ K_n^{p,b} = \left[ \int_0^1 (N_{i,b}^p)'(x) (N_{j,b}^p)'(x) \, dx \right]_{i,j=1}^n. \quad (15)$$

We will conclude our summary of the optimal spline subspaces presented in [32] by introducing a B-spline-like basis for each space. To this end, we need to define the cardinal B-splines and their properties; see, e.g., [10, 31] for more details. Let $\mathcal{N}_p: \mathbb{R} \to \mathbb{R}$ be the cardinal B-spline of degree $p$ recursively defined as

$$\mathcal{N}_0(t) = \begin{cases} 1, & t \in [0, 1), \\ 0, & \text{otherwise,} \end{cases} \quad (16)$$

and

$$\mathcal{N}_p(t) = \frac{t}{p} \mathcal{N}_{p-1}(t) + \frac{p+1-t}{p} \mathcal{N}_{p-1}(t-1), \quad t \in \mathbb{R}, \ p \ge 1. \quad (17)$$

The function class $(\mathcal{N}_p)_{p \ge 0}$ possesses the following properties.

* Smoothness and local support:

$$\mathcal{N}_p \in C^{p-1}(\mathbb{R}), \quad \text{supp}(\mathcal{N}_p) = [0, p+1]. \quad (18)$$
• Symmetry:

$$ \mathcal{N}_p^{(r)} \left( \frac{p+1}{2} + t \right) = (-1)^r \mathcal{N}_p^{(r)} \left( \frac{p+1}{2} - t \right), \quad r \le p. \quad (19) $$

• The $L^2$ inner products:

$$ \begin{aligned} \int_{\mathbb{R}} \mathcal{N}_{p_1}^{(r_1)}(t) \mathcal{N}_{p_2}^{(r_2)}(t+\rho) \, dt &= (-1)^{r_1} \mathcal{N}_{p_1+p_2+1}^{(r_1+r_2)}(p_1+1+\rho) \\ &= (-1)^{r_2} \mathcal{N}_{p_1+p_2+1}^{(r_1+r_2)}(p_2+1-\rho), \end{aligned} \quad (20) $$

for every $\rho \in \mathbb{R}$ and every $r_1 \le p_1$ and $r_2 \le p_2$.

We are now ready to discuss the construction of the B-spline-like basis $B_b^p$ corresponding to $S_{p,n,b}^{\text{opt}}$ for $b = 0, 1, 2$ in the following subsections. For each basis we present its extraction formula from cardinal B-splines. For a more in-depth understanding, we refer the reader to [15, 32].

### 2.2.1 Dirichlet boundary conditions

Let us consider the first optimal spline subspace $S_{p,n,0}^{\text{opt}}$ corresponding to the boundary condition (2) of our eigenvalue problem (1). The B-spline-like basis of this space is denoted by $B_0^p = \{N_{i,0}^p, i = 1, \dots, n\}$ and its definition depends on the parity of $p$ as follows.

Let $m = 2n + 2$. In the case $p$ odd, the basis is computed as

$$ \begin{bmatrix} N_{1,0}^p \\ N_{2,0}^p \\ \vdots \\ N_{n,0}^p \end{bmatrix} = \left[ \underbrace{\cdots \left| L_{(n,m)}^0 \right| L_{(n,m)}^0}_{n+p+1} \mid \underbrace{\cdots \left| R_{(n,m)}^0 \right| R_{(n,m)}^0}_{n+p+1} \mid \cdots \right] \begin{bmatrix} N_{-p}^p \\ \vdots \\ N_0^p \\ N_1^p \\ \vdots \\ N_n^p \end{bmatrix}, \quad (21) $$

where the set $\{N_i^p, i = -p, \dots, n\}$ is derived from the cardinal B-spline $\mathcal{N}_p$ in (17) as

$$ N_i^p(x) = \mathcal{N}_p((n+1)x - i), \quad i = -p, \dots, n, \quad (22) $$

and the matrices $L_{(n,m)}^0$ and $R_{(n,m)}^0$ are composed of the identity matrix $I_n$ of size $n$, the exchange matrix $J_n$ of size $n$, and the zero column vector $O_{(n,1)}$ of length $n$ as

$$ L_{(n,m)}^0 = \left[ I_n \mid O_{(n,1)} \mid -J_n \mid O_{(n,1)} \right], \quad R_{(n,m)}^0 = \left[ O_{(n,1)} \mid -J_n \mid O_{(n,1)} \mid I_n \right]. $$
In the case $p$ even, the basis is computed as

$$ \begin{bmatrix} N_{1,0}^p \\ N_{2,0}^p \\ \vdots \\ N_{n,0}^p \end{bmatrix} = \left[ \underbrace{\cdots | L_{(n,m)}^0 | L_{(n,m)}^0}_{\frac{p}{2}+1} | I_n \underbrace{R_{(n,m)}^0 | R_{(n,m)}^0 | \cdots}_{\frac{p}{2}+1} \right] \begin{bmatrix} N_{-p}^p \\ \vdots \\ N_0^p \\ N_1^p \\ \vdots \\ N_{n+1}^p \end{bmatrix}, \quad (23) $$

where the set $\{N_i^p, i = -p, \dots, n+1\}$ is defined by

$$ N_i^p(x) = \mathcal{N}_p \left( (n+1)x - \left(i - \frac{1}{2}\right) \right), \quad i = -p, \dots, n+1. \quad (24) $$

### 2.2.2 Neumann boundary conditions

Now, we turn our attention to the second optimal spline subspace $S_{p,n,1}^{\text{opt}}$, which is associated with the boundary condition (3). The B-spline-like basis of this space is denoted by $B_1^p = \{N_{i,1}^p, i = 1, \dots, n\}$ and determined as follows, with a dependence on the parity of $p$.

In the case $p$ odd, the basis is computed as

$$ \begin{bmatrix} N_{1,1}^p \\ N_{2,1}^p \\ \vdots \\ N_{n,1}^p \end{bmatrix} = \left[ \underbrace{\cdots | J_n | I_n | J_n}_{\frac{p}{2}} | I_n \underbrace{J_n | I_n | J_n | \cdots}_{\frac{p}{2}+1} \right] \begin{bmatrix} N_{-p}^p \\ \vdots \\ N_0^p \\ N_1^p \\ \vdots \\ N_n^p \end{bmatrix}, \quad (25) $$

where the set $\{N_i^p, i = -p, \dots, n\}$ is given by

$$ N_i^p(x) = \mathcal{N}_p \left( nx - \left(i - \frac{1}{2}\right) \right), \quad i = -p, \dots, n. \quad (26) $$

In the case $p$ even, the basis is computed as

$$ \begin{bmatrix} N_{1,1}^p \\ N_{2,1}^p \\ \vdots \\ N_{n,1}^p \end{bmatrix} = \left[ \underbrace{\cdots | J_n | I_n | J_n}_{\frac{p}{2}} | I_n \underbrace{J_n | I_n | J_n | \cdots}_{\frac{p}{2}+1} \right] \begin{bmatrix} N_{-p}^p \\ \vdots \\ N_0^p \\ N_1^p \\ \vdots \\ N_{n-1}^p \end{bmatrix}, \quad (27) $$

<figure><img src="image_3.png" /></figure>
where the set {$N_i^p, i = -p, \dots, n-1$} is defined by

$$N_i^p(x) = \mathcal{N}_p(nx - i), \quad i = -p, \dots, n-1. \qquad (28)$$

### 2.2.3 Mixed boundary conditions

Finally, let us consider the last optimal spline subspace $S_{p,n,2}^{\text{opt}}$. The B-spline-like basis of this space is denoted by $B_2^p = \{N_{i,2}^p, i = 1, \dots, n\}$ and determined as follows; there is again a dependence on the parity of $p$.

Let $m = 2n + 1$. In the case $p$ odd, the basis is computed as

$$\begin{bmatrix} N_{1,2}^p \\ N_{2,2}^p \\ \vdots \\ N_{n,2}^p \end{bmatrix} = \left[ \underbrace{\left( \cdots \left| L_{(n,m)}^2 \right| - L_{(n,m)}^2 \right) \left| L_{(n,m)}^2 \right|}_{n+p+1} \quad \underbrace{\left( \cdots \left| R_{(n,m)}^2 \right| - R_{(n,m)}^2 \right) \left| R_{(n,m)}^2 \right|}_{n+p+1} \right] \begin{bmatrix} N_{-p}^p \\ \vdots \\ N_0^p \\ N_1^p \\ \vdots \\ N_n^p \end{bmatrix}, \qquad (29)$$

where the set {$N_i^p, i = -p, \dots, n$} is defined by

$$N_i^p(x) = \mathcal{N}_p \left( \frac{2n+1}{2} x - i \right), \quad i = -p, \dots, n, \qquad (30)$$

and the matrices $L_{(n,m)}^2$ and $R_{(n,m)}^2$ are defined by

$$L_{(n,m)}^2 = [-I_n \mid -J_n \mid O_{(n,1)}], \quad R_{(n,m)}^2 = [J_n \mid O_{(n,1)} \mid -I_n].$$

In the case $p$ even, the basis is computed as

$$\begin{bmatrix} N_{1,2}^p \\ N_{2,2}^p \\ \vdots \\ N_{n,2}^p \end{bmatrix} = \left[ \underbrace{\left( \cdots \left| L_{(n,m)}^2 \right| - L_{(n,m)}^2 \right) \left| L_{(n,m)}^2 \right|}_{n+p+1} \quad \underbrace{\left( \cdots \left| R_{(n,m)}^2 \right| - R_{(n,m)}^2 \right) \left| R_{(n,m)}^2 \right|}_{n+p+1} \right] \begin{bmatrix} N_{-p}^p \\ \vdots \\ N_0^p \\ N_1^p \\ \vdots \\ N_n^p \end{bmatrix}, \qquad (31)$$

where the set {$N_i^p, i = -p, \dots, n$} is given by

$$N_i^p(x) = \mathcal{N}_p \left( \frac{2n+1}{2} x - \left(i - \frac{1}{2}\right) \right), \quad i = -p, \dots, n. \qquad (32)$$
## 2.3 Structured matrices of Toeplitz-minus-Hankel or Toeplitz-plus-Hankel form

We describe three known results regarding closed-form expressions for eigenvalues and eigenvectors of structured matrices of Toeplitz-minus-Hankel or Toeplitz-plus-Hankel form. Further details and other types of matrices can be found in [11] (see also [4, Table 1.2]). All theorems presented here are taken from [11], adjusted to align more effectively with the notation of the GLT symbol presented in [24]. Additionally, various matrix algebras can be found in [3, 5, 14], illustrating the connections between these fields.

Suppose a given integer number $n \in \mathbb{N}^*$ and a given function $a$ belonging to the Lebesgue space $L^1([-\pi, \pi])$. The Toeplitz matrix of size $n$ generated by $a$ is defined by

$$ T_n(a) = [a_{i-j}]_{i,j=1}^n = \begin{bmatrix} a_0 & a_{-1} & a_{-2} & \cdots & a_{-(n-1)} \\ a_1 & \ddots & \ddots & \ddots & \vdots \\ a_2 & \ddots & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & \ddots & a_{-2} \\ \vdots & \ddots & \ddots & \ddots & a_{-1} \\ a_{n-1} & \cdots & \cdots & a_2 & a_1 & a_0 \end{bmatrix}, $$

where $a_k$ represents the $k$-th Fourier coefficient of $a$, given by

$$ a_k = \frac{1}{2\pi} \int_{-\pi}^{\pi} a(\theta) e^{-ik\theta} \, d\theta, \quad k \in \mathbb{Z}, $$

with $i^2 = -1$.

Let $p$ be a positive integer such that $n \ge p + 1$. Given the vector $\alpha = (\alpha_0, \alpha_1, \dots, \alpha_p) \in \mathbb{R}^{p+1}$, we define the symmetric Toeplitz matrix of size $n$ associated with $\alpha$ as follows

$$ T_n^\alpha = \begin{bmatrix} \alpha_0 & \alpha_1 & \cdots & \alpha_p & 0 \\ \alpha_1 & \alpha_0 & \alpha_1 & \cdots & \alpha_p \\ \vdots & \alpha_1 & \alpha_0 & \alpha_1 & \cdots & \alpha_p \\ \alpha_p & \cdots & \alpha_1 & \alpha_0 & \alpha_1 & \cdots & \alpha_p \\ 0 & \cdots & \alpha_p & \cdots & \alpha_1 & \alpha_0 & \alpha_1 \\ \alpha_p & \cdots & \alpha_1 & \alpha_0 & \alpha_1 & \cdots & \alpha_p \\ \vdots & \alpha_p & \cdots & \alpha_1 & \alpha_0 & \alpha_1 & \alpha_p \\ 0 & \cdots & \alpha_p & \cdots & \alpha_1 & \alpha_0 & \alpha_p \end{bmatrix}. \quad (33) $$
Through straightforward calculations, we can deduce the function that generates this matrix. More precisely, we have

$$T_n^\alpha = T_n(g_p^\alpha),$$

where

$$g_p^\alpha(\theta) = \alpha_0 + 2 \sum_{k=1}^{p} \alpha_k \cos(k\theta), \quad \forall \theta \in [0, \pi]. \qquad (34)$$

Before delving into the results, it is essential to introduce two distinct types of Hankel matrices of size n that are directly linked to the vector α:

$$H_n^{\alpha,1} = \begin{bmatrix} \alpha_1 & \alpha_2 & \alpha_3 & \cdots & \alpha_p & 0 \\ \alpha_2 & \alpha_3 & \alpha_4 & \cdots & \alpha_p & 0 \\ \alpha_3 & \alpha_4 & \alpha_5 & \cdots & \alpha_p & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \alpha_p & \alpha_1 & \alpha_2 & \cdots & \alpha_{p-1} & 0 \\ 0 & \alpha_p & \alpha_1 & \cdots & \alpha_{p-1} & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \alpha_p & \alpha_1 & \alpha_2 & \cdots & \alpha_{p-1} & 0 \\ \alpha_p & \alpha_1 & \alpha_2 & \cdots & \alpha_{p-1} & 0 \end{bmatrix}, \quad H_n^{\alpha,2} = \begin{bmatrix} \alpha_2 & \alpha_3 & \alpha_4 & \cdots & \alpha_p & 0 \\ \alpha_3 & \alpha_4 & \alpha_5 & \cdots & \alpha_p & 0 \\ \alpha_4 & \alpha_5 & \alpha_6 & \cdots & \alpha_p & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \alpha_p & \alpha_1 & \alpha_2 & \cdots & \alpha_{p-1} & 0 \\ 0 & \alpha_p & \alpha_1 & \cdots & \alpha_{p-1} & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\ \alpha_p & \alpha_1 & \alpha_2 & \cdots & \alpha_{p-1} & 0 \\ \alpha_p & \alpha_1 & \alpha_2 & \cdots & \alpha_{p-1} & 0 \end{bmatrix}. \qquad (35)$$

Note that the initial element of the Hankel matrix $H_n^{\alpha,1}$ corresponds to the second element of the Toeplitz matrix (33). Similarly, the Hankel matrix $H_n^{\alpha,2}$ begins with the third element of the same Toeplitz matrix.

Now we have the necessary background to introduce the following theorems; their proofs can be found in [11].

**Theorem 2** Let $n, p \in \mathbb{N}^*$ such that $n \ge p + 1$. Then, for $j = 1, \dots, n$, the j-th eigenpair of the matrix $T_n^\alpha - H_n^{\alpha,2}$ can be expressed as

$$\lambda_j (T_n^\alpha - H_n^{\alpha,2}) = g_p^\alpha \left( \frac{j\pi}{n+1} \right),$$

and

$$u_{i,j} (T_n^\alpha - H_n^{\alpha,2}) = \sqrt{\frac{2}{n+1}} \sin \left( \frac{ij\pi}{n+1} \right),$$

with $i = 1, \dots, n$. The function $g_p^\alpha$ is defined in (34).
**Theorem 3** Let $n, p \in \mathbb{N}^*$ such that $n \ge p + 1$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $T_n^\alpha + H_n^{\alpha, 1}$ can be expressed as

$$ \lambda_j (T_n^\alpha + H_n^{\alpha, 1}) = g_p^\alpha \left( \frac{(j-1)\pi}{n} \right), $$

and

$$ u_{i,j} (T_n^\alpha + H_n^{\alpha, 1}) = \sqrt{\frac{2}{n}} c_j \cos \left( \frac{(j-1)\pi}{n} \left( i - \frac{1}{2} \right) \right), \quad c_j = \begin{cases} \frac{1}{\sqrt{2}}, & j = 1, \\ 1, & j \ge 2, \end{cases} $$

with $i = 1, \dots, n$. The function $g_p^\alpha$ is defined in (34).

**Theorem 4** Let $n, p \in \mathbb{N}^*$ such that $n \ge p + 1$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $T_n^\alpha - H_n^{\alpha, 1}$ can be expressed as

$$ \lambda_j (T_n^\alpha - H_n^{\alpha, 1}) = g_p^\alpha \left( \frac{j\pi}{n} \right), $$

and

$$ u_{i,j} (T_n^\alpha - H_n^{\alpha, 1}) = \sqrt{\frac{2}{n}} c_j \sin \left( \frac{j\pi}{n} \left( i - \frac{1}{2} \right) \right), \quad c_j = \begin{cases} \frac{1}{\sqrt{2}}, & j = n, \\ 1, & j \le n-1, \end{cases} $$

with $i = 1, \dots, n$. The function $g_p^\alpha$ is defined in (34).

**Remark 5** Theorems 2, 3, and 4 precisely mirror those in [11, Theorems 2.1, 2.4, and 2.2], differing only in the normalization of the eigenvectors. Furthermore, an alternative version of Theorem 2 can be found in [3, Proposition 2.2], wherein the eigenvalues are provided implicitly.

**Remark 6** From Appendix B, it is clear that the matrices in Theorems 2, 3, and 4 belong to certain matrix algebras, namely

$$ T_n^\alpha - H_n^{\alpha, 2} \in \tau_n(0, 0), \quad T_n^\alpha + H_n^{\alpha, 1} \in \tau_n(1, 1), \quad T_n^\alpha - H_n^{\alpha, 1} \in \tau_n(-1, -1). $$

## 3 Eigenpairs of the optimal spline discretization matrices

Our objective is to establish closed-form expressions for the eigenvalues and eigenvectors of the matrices $M_n^{p,b}$ and $K_n^{p,b}$, $b = 0, 1, 2$, defined in (15) and arising from the discretization in the spaces $S_{p,n,b}^{\text{opt}}$ related to the boundary conditions (2)–(4). For the sake of unification, we introduce a general matrix $X^{p,r,b}$ as follows.

**Definition 7** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $r \le p$. Let $B_b^p = \{N_{i,b}^p, i = 1, \dots, n\}$ be the basis constructed in Sections 2.2.1–2.2.3 for $b = 0, 1, 2$,
respectively. We consider the matrix $X^{p,r,b} = [X_{i,j}^{p,r,b}]_{i,j=1}^n$, whose elements are given by

$$X_{i,j}^{p,r,b} = \int_0^1 (N_{i,b}^p)^{(r)}(x)(N_{j,b}^p)^{(r)}(x) \, dx, \quad i, j = 1, \dots, n,$$

for $b = 0, 1, 2$. In particular, in the cases $r = 0, 1$, we have $X^{p,0,b} = M_n^{p,b}$ and $X^{p,1,b} = K_n^{p,b}$.

In this section, by leveraging the properties (18)–(20) of cardinal B-splines, we show that the matrices in Definition 7 exhibit a Toeplitz-minus-Hankel or Toeplitz-plus-Hankel structure, for all the values of $p$. Then, the results from Section 2.3 lead to closed-form expressions for the eigenvalues and eigenvectors of our matrices of interest. The section is divided into three subsections, each corresponding to one of the boundary conditions (2)–(4). Since all three boundary cases can be treated using the same technique, we primarily detail the Dirichlet case. However, to enhance readability, we delegate the technical proofs to Appendix A. As a side result, for the mixed boundary conditions (4), we identify a new Hankel structure, not covered in [11].

## 3.1 The Dirichlet case

Let us consider the discretization of the eigenvalue problem (1) in the space $S_{p,n,0}^{\text{opt}}$ defined in (12) and related to the boundary conditions (2). We aim to show that the matrix $X^{p,r,0}$ in Definition 7 has a Toeplitz-minus-Hankel form. Before going to the main result, two essential lemmas are in order. The first lemma is a key ingredient in the proof and the second one contributes to its reduction.

**Lemma 8** Let $a, b, \rho \in \mathbb{R}$, $p \in \mathbb{N}^*$, and $r \in \mathbb{N}$ such that $r \le p$. Then,

$$\int_a^b \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+\rho) \, dt = \int_{p+1-b-\rho}^{p+1-a-\rho} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+\rho) \, dt.$$

**Proof** Fix $a, b, \rho \in \mathbb{R}$. By the symmetry property of cardinal B-splines (19) we have

$$\mathcal{N}_p^{(r)}(t) = (-1)^r \mathcal{N}_p^{(r)}(p+1-t),$$

and

$$\mathcal{N}_p^{(r)}(t+\rho) = (-1)^r \mathcal{N}_p^{(r)}(p+1-t-\rho).$$

Integrating the product of these functions and using a change of variables ($t \leftarrow p+1-t-\rho$) gives the stated result. □

**Lemma 9** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $r \le p$. Then, the matrix $X^{p,r,0}$ is both symmetric and centrosymmetric.
**Proof** From Definition 7 it is clear that $X^{p,r,0}$ is a symmetric matrix. Moreover, let

$$p_0 = p - 2 \lfloor \frac{p}{2} \rfloor,$$

and from the properties of cardinal B-splines (18)–(20) we deduce

$$N_i^p(x) = N_{n+1-p_0-p-i}^p(1-x), \quad i = -p, \dots, n+1-p_0, \quad x \in [0, 1],$$

where $N_i^p$ is defined in (22) and (24). Note that $p_0$ allows us to unify the cases of even and odd degree. Then, from the constructions (21) and (23) it follows

$$N_{i,0}^p(x) = N_{n+1-i,0}^p(1-x), \quad i = 1, \dots, n, \quad x \in [0, 1],$$

which implies that $X^{p,r,0}$ is a centrosymmetric matrix. □

**Remark 10** The symmetry and centrosymmetry of the matrix $X^{p,r,0}$ hold without any restriction on the size $n$. Therefore, the insights given in [7] become particularly valuable for studying the eigenvectors, eigenvalues, and the development of fast solvers.

In the next theorem, we establish that the matrix $X^{p,r,0}$ takes the form of $T_n^\alpha - H_n^{\alpha,2}$ for a specific vector $\alpha$, where $T_n^\alpha$ and $H_n^{\alpha,2}$ are defined in (33) and (35), respectively. From Theorem 2, it is evident that a closed-form expression for the eigenvalues and eigenvectors requires at least $n \ge p+1$. However, in our case, we need to assume $n \ge C(p) \ge p+1$ and we aim to determine the smallest possible value for $C(p)$ that gives us the structure $T_n^\alpha - H_n^{\alpha,2}$ without any perturbation.

**Theorem 11** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $n \ge \max\{p+1, p+\lfloor \frac{p}{2} \rfloor - 1\}$ and $r \le p$. Then,

$$X^{p,r,0} = (n+1)^{2r-1} (T_n^\alpha - H_n^{\alpha,2}), \quad (36)$$

where the elements of the vector $\alpha$ are given by

$$\alpha_k = (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1-k), \quad k = 0, \dots, p. \quad (37)$$

**Proof** For ease of reading, the detailed proof is moved to Section A.1 (see Appendix A) and we only sketch the idea here. Assuming $n \ge p+1$, we start by providing a simplified definition of the basis $B_0^p$, derived from (21) and (23). This allows us to obtain explicit expressions for the elements of the matrix $X^{p,r,0}$ by leveraging the properties of cardinal B-splines (18)–(20) and utilizing Lemma 8. We do this in two steps. In the first step, we show that the central part of the matrix corresponds to a Toeplitz form,

$$X_{i,j}^{p,r,0} = (n+1)^{2r-1} (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1-i-j), \\ i, j = \lfloor \frac{p}{2} \rfloor + 1, \dots, n - \lfloor \frac{p}{2} \rfloor.$$
In the second step, under the stated conditions on $n$ and $p$, we show that the first block of size $(\lfloor \frac{p}{2} \rfloor, n)$ corresponds to a Toeplitz-minus-Hankel form,

$$X_{i,j}^{p,r,0} = (n+1)^{2r-1} \left( (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i-j) - (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i+j) \right), \\ i = 1, \dots, \lfloor \frac{p}{2} \rfloor, \quad j = i, \dots, n.$$

Exploiting the symmetry and centrosymmetry of $X^{p,r,0}$ (see Lemma 9) completes the proof. □

**Remark 12** As a generalization of $X^{p,r,0}$, one can consider two indices $r_1, r_2 \in \mathbb{N}$ instead of just $r$ in Definition 7. Then, the same analysis can be followed to derive a more diverse range of Toeplitz-minus/plus-Hankel matrices that depend on the pair $(r_1, r_2)$. For brevity, we omit this analysis.

For $\theta \in [0, \pi]$ and $0 \le r \le p$, let us set

$$g_p^r(\theta) = (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1) + 2(-1)^r \sum_{k=1}^{p} \mathcal{N}_{2p+1}^{(2r)}(p+1-k) \cos(k\theta). \quad (38)$$

The next corollary follows directly from Theorems 2 and 11, taking into account that the function $g_p^r$ in (38) is equal to $g_p^\alpha$ in (34), with the elements of $\alpha$ specified in (37).

**Corollary 13** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $n \ge \max\{p+1, p+\lfloor \frac{p}{2} \rfloor - 1\}$ and $r \le p$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $X^{p,r,0}$ can be expressed as

$$\lambda_j(X^{p,r,0}) = (n+1)^{2r-1} g_p^r \left( \frac{j\pi}{n+1} \right),$$

and

$$u_{i,j}(X^{p,r,0}) = \sqrt{\frac{2}{n+1}} \sin \left( \frac{ij\pi}{n+1} \right),$$

with $i = 1, \dots, n$. The function $g_p^r$ is defined in (38).

We now state a corollary regarding the closed form of eigenvalues and eigenvectors of the discretized Laplace operator (14) with the Dirichlet boundary conditions (2). We recall that the matrices $M_n^{p,0}$ and $K_n^{p,0}$ defined in (15) represent a particular case of $X^{p,r,0}$ ($r = 0, 1$).

**Corollary 14** Let $p, n \in \mathbb{N}^*$ such that $n \ge \max\{p+1, p+\lfloor \frac{p}{2} \rfloor - 1\}$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $L_n^{p,0} = (M_n^{p,0})^{-1} K_n^{p,0}$ can be expressed
as

$$ \lambda_j (L_n^{p,0}) = (n+1)^2 \frac{g_p^1 \left( \frac{j\pi}{n+1} \right)}{g_p^0 \left( \frac{j\pi}{n+1} \right)}, $$

and

$$ u_{i,j} (L_n^{p,0}) = \sqrt{\frac{2}{n+1}} \sin \left( \frac{ij\pi}{n+1} \right), $$

with $i = 1, \dots, n$. The function $g_p^r$ with $r \in \{0, 1\}$ is defined in (38).

**Proof** The results follow directly from Corollary 13 by considering the cases $r = 0$ and $r = 1$. Additionally, it was shown in [21] that

$$ \left(\frac{4}{\pi^2}\right)^{p+1} \leq g_p^0(\theta) \leq 1, \quad \theta \in [0, \pi], $$

ensuring the validity of the division by $g_p^0 \left( \frac{j\pi}{n+1} \right)$, $j = 1, \dots, n$. $\square$

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

**Remark 15** The functions $g_p^1$ and $g_p^0$ are the so-called (GLT) symbols of the stiffness and mass matrices resulting from the Galerkin discretization in the full spline space $S_{p,\tau}$, where $\tau$ corresponds to a uniform grid. These symbols are independent of the boundary conditions. There is an extensive literature on their construction and on the analysis of their properties; see, e.g., [16, 18, 21] and references therein. The above results also show that $g_p^1$ and $g_p^0$ are the (GLT) symbols of the same matrices when the discretization is performed in optimal spline subspaces (or in reduced spline spaces; see Section 4.1) of the same degree.

**Remark 16** The explicit closed-form expressions, described in Corollary 14, for the approximated spectrum of the Laplace operator with Dirichlet boundary conditions give an algebraic confirmation that the optimal spline subspaces $S_{p,n,0}^{\text{opt}}$ lead to outlier-free discretizations. Indeed, it was shown in [18, Theorem 1] that a uniform sampling of the function

$$ e_p(\theta) = \frac{g_p^1(\theta)}{g_p^0(\theta)}, \quad \theta \in [0, \pi] $$

provides an approximation of the foreseen spectrum, with no outliers; see also [27]. In particular, by [18, Lemma 1] we have for $p \ge 1$ and $\theta \in [0, \pi]$,

$$ 0 \le \frac{e_p(\theta) - \theta^2}{\theta^2} \le \frac{4\pi(\pi - \theta)}{(2\pi - \theta)^2} \left(\frac{\theta}{2\pi - \theta}\right)^{2p} + 5 \left(\frac{\theta}{2\pi + \theta}\right)^{2p}. \quad (39) $$

<figure><img src="image_5.png" /></figure>
By Corollary 14 and setting $\theta_j = \frac{j\pi}{n+1}$, $j = 1, \dots, n$, we get

$$ \frac{\lambda_j (L_n^{p,0}) - (j\pi)^2}{(j\pi)^2} = \frac{(n+1)^2}{(j\pi)^2} (e_p(\theta_j) - (\theta_j)^2) = \frac{e_p(\theta_j) - (\theta_j)^2}{(\theta_j)^2}. $$

Thus, for all $j = 1, \dots, n$, the inequalities in (39) show that $\lambda_j (L_n^{p,0})$ is an approximation of the $j$-th eigenvalue of the Laplace operator with Dirichlet boundary conditions, see (5), whose relative error converges to zero as $p$ increases. Thus, the approximation is outlier-free.

## 3.2 The Neumann case

After addressing the Dirichlet boundary conditions (2) in Section 3.1, our attention shifts to the Neumann boundary conditions (3). Although the proof of the main result and the analysis in this section is similar to the Dirichlet case, we still include it with the necessary details in Appendix A for two reasons. Firstly, this serves as an illustration how to handle cases not explicitly covered in this work. Secondly, and more importantly, the analysis of the reduced spline space proposed in Section 4.1 can be viewed as a small variation of the Neumann scenario. Consequently, the results obtained here will straightforwardly imply the results of the reduced spline space.

Let us consider the discretization of the eigenvalue problem (1) in the space $S_{p,n,1}^{\text{opt}}$ defined in (12) and related to the boundary conditions (3). The matrix $X^{p,r,1}$ in Definition 7 exhibits the structure of $T_n^\alpha + H_n^{\alpha,1}$ for a specific vector $\alpha$, where $T_n^\alpha$ and $H_n^{\alpha,1}$ are given by (33) and (35), respectively. Following a proof akin to that of Lemma 9, it can be shown that $X^{p,r,1}$ is centrosymmetric. This property will be exploited in the proof of the next theorem, which shows the Toeplitz-plus-Hankel form of $X^{p,r,1}$.

**Theorem 17** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $n \ge \max\{2p - \lfloor \frac{p}{2} \rfloor, 2p - 2\lfloor \frac{p}{2} \rfloor + 1\}$ and $r \le p$. Then,

$$ X^{p,r,1} = n^{2r-1} (T_n^\alpha + H_n^{\alpha,1}), \quad (40) $$

where the elements of the vector $\alpha$ are given by (37).

**Proof** A similar proof strategy can be employed to the one for Theorem 11 in the Dirichlet case. The detailed steps can be found in Section A.2 (see Appendix A). □

By combining Theorems 3 and 17 we directly obtain a closed-form expression for the eigenvalues and eigenvectors of the matrix $X^{p,r,1}$.

**Corollary 18** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $n \ge \max\{2p - \lfloor \frac{p}{2} \rfloor, 2p - 2\lfloor \frac{p}{2} \rfloor + 1\}$ and $r \le p$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $X^{p,r,1}$ can be expressed as

$$ \lambda_j (X^{p,r,1}) = n^{2r-1} g_p^r \left( \frac{(j-1)\pi}{n} \right), $$
and

$$u_{i,j}(X^{p,r,1}) = \sqrt{\frac{2}{n}} c_j \cos\left(\frac{(j-1)\pi}{n}\left(i - \frac{1}{2}\right)\right), \quad c_j = \begin{cases} \frac{1}{\sqrt{2}}, & j = 1, \\ 1, & j \ge 2, \end{cases}$$

with $i = 1, \dots, n$. The function $g_p^r$ is defined in (38).

Thanks to Corollary 18 and the definition of $X^{p,r,1}$ ($r = 0, 1$), we conclude with the following result on the eigenpairs of the system (14) related to the basis $B_1^p$.

**Corollary 19** Let $p, n \in \mathbb{N}^*$ such that $n \ge \max\{2p - \lfloor \frac{p}{2} \rfloor, 2p - 2\lfloor \frac{p}{2} \rfloor + 1\}$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $L_n^{p,1} = (M_n^{p,1})^{-1} K_n^{p,1}$ can be expressed as

$$\lambda_j(L_n^{p,1}) = n^2 \frac{g_p^1\left(\frac{(j-1)\pi}{n}\right)}{g_p^0\left(\frac{(j-1)\pi}{n}\right)},$$

and

$$u_{i,j}(L_n^{p,1}) = \sqrt{\frac{2}{n}} c_j \cos\left(\frac{(j-1)\pi}{n}\left(i - \frac{1}{2}\right)\right), \quad c_j = \begin{cases} \frac{1}{\sqrt{2}}, & j = 1, \\ 1, & j \ge 2, \end{cases}$$

with $i = 1, \dots, n$. The function $g_p^r$ with $r \in \{0, 1\}$ is defined in (38).

**Proof** The argument of the proof is the same as that of Corollary 14. □

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

**Remark 20** Corollary 19 gives an algebraic confirmation that the optimal spline subspaces $S_{p,n,1}^{\text{opt}}$ lead to outlier-free discretizations for the spectrum of the Laplace operator with Neumann boundary conditions. This follows from (6) in the same way as described in Remark 16. In particular, by setting $\theta_j = \frac{(j-1)\pi}{n}$, $j = 2, \dots, n$, we get

$$\frac{\lambda_j(L_n^{p,1}) - ((j-1)\pi)^2}{((j-1)\pi)^2} = \frac{n^2}{((j-1)\pi)^2} (e_p(\theta_j) - (\theta_j)^2) = \frac{e_p(\theta_j) - (\theta_j)^2}{(\theta_j)^2}.$$

Moreover, $\lambda_1(L_n^{p,1})$ agrees with the exact value of the first eigenvalue in (6) because $g_p^1(0) = 0$.

### 3.3 The mixed case

The mixed boundary conditions (4) represent a combination of a Dirichlet boundary condition and a Neumann boundary condition. Since, by nature, these boundary conditions are not symmetric, the mass and stiffness matrices arising from this case

<figure><img src="image_7.png" /></figure>
are not centrosymmetric. However, these matrices can be described by combining the matrices of the previously studied boundary conditions and this gives rise to a new class of matrices. We also establish closed-form expressions for the eigenvalues and eigenvectors of this type of matrices. The same reasoning can be applied if we permute the mixed boundary conditions.

We define a new Hankel matrix, which combines the Hankel matrices of the Dirichlet case ($H_n^{\alpha,2}$) and the Neumann case ($H_n^{\alpha,1}$):

$$H_n^{\alpha,2,1} = \begin{bmatrix} -\alpha_2 & -\alpha_3 & -\alpha_4 & \cdots & -\alpha_p & 0 \\ -\alpha_3 & -\alpha_4 & \cdots & & -\alpha_p & \\ -\alpha_4 & \cdots & -\alpha_p & & & \\ \vdots & \ddots & & & & \vdots \\ -\alpha_p & & & & & 0 \\ 0 & & & & & \alpha_p \\ & & & & \alpha_p & \cdots & \alpha_3 \\ & & \alpha_p & \cdots & \alpha_3 & \alpha_2 & \alpha_1 \\ & & & \alpha_p & \cdots & \alpha_3 & \alpha_2 & \alpha_1 \end{bmatrix}$$

The next theorem is the key to the closed-form expression of the eigenpairs considered in this section. It can be proved using the same technique as done in [11, Theorem 2.1]; for brevity, we omit the details.

**Theorem 21** Let $n, p \in \mathbb{N}^*$ such that $n \ge p + 1$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $T_n^\alpha + H_n^{\alpha,2,1}$ can be expressed as

$$\lambda_j (T_n^\alpha + H_n^{\alpha,2,1}) = g_p^\alpha \left( \frac{(2j-1)\pi}{2n+1} \right),$$

and

$$u_{i,j} (T_n^\alpha + H_n^{\alpha,2,1}) = \sqrt{\frac{4}{2n+1}} \sin \left( \frac{i(2j-1)\pi}{2n+1} \right),$$

with $i = 1, \dots, n$. The function $g_p^\alpha$ is defined in (34).

**Remark 22** From Appendix B it is clear that the matrix $T_n^\alpha + H_n^{\alpha,2,1}$ falls within the matrix algebra $\tau_n(0, 1)$. This algebra is referred to as the algebra $\Gamma$ in [14] and initially introduced in [1].

Now, with all the necessary prerequisites, we present the main result on the Toeplitz-plus-Hankel form of the matrix $X^{p,r,2}$ in Definition 7.

**Theorem 23** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $n \ge \max\{p+1, p+\lfloor\frac{p}{2}\rfloor\}$ and $r \le p$. Then,

$$X^{p,r,2} = \left(\frac{2n+1}{2}\right)^{2r-1} (T_n^\alpha + H_n^{\alpha,2,1}),$$
where the elements of the vector $\alpha$ are given by (37).

**Proof** We can again employ a similar proof strategy as in the previous boundary cases (Theorems 11 and 17). However, due to the lack of centrosymmetry of the matrix $X^{p,r,2}$, we need to do some extra work in our proof. More precisely, we divide the second step into two parts: the first part is akin to the Dirichlet case, while the second part resembles the Neumann case. The detailed steps are provided in Section A.3 (see Appendix A). □

Since the structure considered in Theorem 23 is a specific case of Theorem 21, we directly arrive at the following result.

**Corollary 24** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $n \ge \max\{p+1, p+\lfloor\frac{p}{2}\rfloor\}$ and $r \le p$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $X^{p,r,2}$ can be expressed as

$$ \lambda_j (X^{p,r,2}) = \left(\frac{2n+1}{2}\right)^{2r-1} g_p^r \left(\frac{(2j-1)\pi}{2n+1}\right), $$

and

$$ u_{i,j} (X^{p,r,2}) = \sqrt{\frac{4}{2n+1}} \sin \left( \frac{i(2j-1)\pi}{2n+1} \right), $$

with $i = 1, \dots, n$. The function $g_p^r$ is defined in (38).

The closed form of the eigenpairs of the system (14) related to the basis $B_2^p$ follows from the above corollary, taking $r = 0, 1$.

**Corollary 25** Let $p, n \in \mathbb{N}^*$ such that $n \ge \max\{p+\lfloor\frac{p}{2}\rfloor, p+1\}$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $L_n^{p,2} = (M_n^{p,2})^{-1} K_n^{p,2}$ can be expressed as

$$ \lambda_j (L_n^{p,2}) = \left(\frac{2n+1}{2}\right)^2 \frac{g_p^1 \left(\frac{(2j-1)\pi}{2n+1}\right)}{g_p^0 \left(\frac{(2j-1)\pi}{2n+1}\right)}, $$

and

$$ u_{i,j} (L_n^{p,2}) = \sqrt{\frac{4}{2n+1}} \sin \left( \frac{i(2j-1)\pi}{2n+1} \right), $$

with $i = 1, \dots, n$. The function $g_p^r$ with $r \in \{0, 1\}$ is defined in (38).

**Remark 26** Corollary 25 gives an algebraic confirmation that the optimal spline subspaces $S_{p,n,2}^{\text{opt}}$ lead to outlier-free discretizations for the spectrum of the Laplace operator with mixed boundary conditions. This follows from (7) in the same way
as described in Remark 16. In particular, by setting $\theta_j = \frac{(2j-1)\pi}{2n+1}$, $j = 1, \dots, n$, we get

$$ \frac{4\lambda_j (L_n^{p,1}) - ((2j-1)\pi)^2}{((2j-1)\pi)^2} = \frac{(2n+1)^2}{((2j-1)\pi)^2} (e_p(\theta_j) - (\theta_j)^2) = \frac{e_p(\theta_j) - (\theta_j)^2}{(\theta_j)^2} $$

# 4 Other outlier-free spline spaces

Upon completing the derivation of closed-form expressions for the eigenpairs of the discretization matrices in optimal spline subspaces, we now aim to illustrate that a similar analysis can be seamlessly applied to other (univariate) reduced spline spaces, known in the outlier-free IGA literature. This reinforces the flexibility and applicability of our methodology for addressing various scenarios. Furthermore, we expand all our findings from the univariate to the multivariate tensor-product setting. For simplicity of presentation, in this section, we concentrate exclusively on Dirichlet boundary conditions but similar results hold for other boundary conditions as well.

## 4.1 Reduced spline spaces

Here, we address the reduced spline spaces considered in [28] (see also [32, 38, 39]). We start by defining these spaces in the Dirichlet case (2) and we recall the B-spline-like basis from [32]. Subsequently, we present the main result, which directly stems from the analysis provided for the optimal spline subspace addressed in Section 3.2.

Let us approximate the eigenvalue problem (1) under the Dirichlet boundary conditions (2) using the reduced spline space $\bar{S}_{p,n,0}$, given by

$$ \bar{S}_{p,n,0} = S_{p,\bar{\tau},0}^{p-1} $$

where $S_{p,\bar{\tau},0}^{p-1}$ is defined in (9), taking $\ell = p-1$, and $\bar{\tau}$ corresponds to the uniform grid

$$ \bar{\tau} = \left( 0, \frac{1}{n}, \frac{2}{n}, \dots, \frac{n-1}{n}, 1 \right). $$

It has been observed numerically in [28] that these spaces are outlier-free. Here we confirm these observations theoretically and actually provide closed-form expressions for the discrete eigenvalues and eigenvectors of the matrices arising from such discretizations.

We first remark that $S_{p,\bar{\tau},0}^{p} = S_{p,\bar{\tau},0}^{p-1}$ and $\bar{S}_{p,n+1,0} = S_{p,n,0}^{\text{opt}}$ for $p$ odd; see also [32, Remark 4.1]. Therefore, the results from Section 3.1 are valid here for $p$ odd and it suffices to restrict the analysis to the case of $p$ even. Hereafter, we assume that $p$ is even in this section, and we have $\dim(\bar{S}_{p,n,0}) = n$. From [32], the B-spline-like basis
$\overline{B}_0^p = \{ \overline{N}_{i,0}^p, i = 1, \dots, n \}$ is computed as

$$ \left[ \begin{array}{c} \overline{N}_{1,0}^p \\ \overline{N}_{2,0}^p \\ \vdots \\ \overline{N}_{n,0}^p \end{array} \right] = \left[ \left\{ \overline{L}_{(n,2n)} \left| \overline{L}_{(n,2n)} \right. \right\}_{n+p} \right] \left[ \begin{array}{c} N_{-p}^p \\ \vdots \\ N_0^p \\ N_1^p \\ \vdots \\ N_{n-1}^p \end{array} \right], \quad (41) $$

where the set $\{N_i^p, i = -p, \dots, n-1\}$ is derived from the cardinal B-spline $N_p$ in (17) as

$$ N_i^p(x) = N_p(nx - i), \quad i = -p, \dots, n-1, $$

and the matrices $\overline{L}_{(n,2n)}$ and $\overline{R}_{(n,2n)}$ are defined by

$$ \overline{L}_{(n,2n)} = [I_n \mid -J_n], \quad \overline{R}_{(n,2n)} = [-J_n \mid I_n]. $$

Let $\overline{M}_n^{p,0}$ and $\overline{K}_n^{p,0}$ denote the mass matrix and the stiffness matrix, respectively, related to the basis $\overline{B}_0^p$. For the sake of unification, similar to Definition 7, we introduce a general matrix $\overline{X}^{p,r,0}$.

**Definition 27** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $p$ even and $r \le p$. Let $\overline{B}_0^p = \{\overline{N}_{i,0}^p, i = 1, \dots, n\}$ be the basis constructed in (41). We consider the matrix $\overline{X}^{p,r,0} = [\overline{X}_{i,j}^{p,r,0}]_{i,j=1}^n$, whose elements are given by

$$ \overline{X}_{i,j}^{p,r,0} = \int_0^1 (\overline{N}_{i,0}^p)^{(r)}(x) (\overline{N}_{j,0}^p)^{(r)}(x) \, dx, \quad i, j = 1, \dots, n. $$

In particular, in the cases $r = 0, 1$, we have $\overline{X}^{p,0,0} = \overline{M}_n^{p,0}$ and $\overline{X}^{p,1,0} = \overline{K}_n^{p,0}$.

Drawing upon the proof of Lemma 9, one can check that the matrix $\overline{X}^{p,r,0}$ is centrosymmetric without any restriction on the size $n$. In the next theorem, we establish that $\overline{X}^{p,r,0}$ takes the form of $T_n^\alpha - H_n^{\alpha,1}$, where $T_n^\alpha$ and $H_n^{\alpha,1}$ are defined by (33) and (35), respectively. Then, by applying Theorem 4, we obtain a closed-form expression for the eigenvalues and eigenvectors.

**Theorem 28** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $p$ even, $n \ge p + \frac{p}{2}$, and $r \le p$. Then,

$$ \overline{X}^{p,r,0} = n^{2r-1} (T_n^\alpha - H_n^{\alpha,1}), \quad (42) $$
where the elements of the vector $\alpha$ are given by (37).

**Proof** Using the assumption $n \ge p + 1$, the expression (41) simplifies to

$$ \begin{bmatrix} \bar{N}_{1,0}^p \\ \bar{N}_{2,0}^p \\ \vdots \\ \bar{N}_{n,0}^p \end{bmatrix} = \left[ \underbrace{-J_n}_{n+p} \quad \underbrace{I_n}_{n+p} \quad \underbrace{-J_n}_{n+p} \right] \begin{bmatrix} N_{-p}^p \\ \vdots \\ N_0^p \\ N_1^p \\ \vdots \\ N_{n-1}^p \end{bmatrix}, $$

and we obtain

$$ \bar{N}_{i,0}^p = \begin{cases} -N_{q-i}^p + N_{q+i-1}^p, & i \in \{1, \dots, -q\}, \\ N_{q+i-1}^p, & i \in \{-q+1, \dots, n+q\}, \\ N_{q+i-1}^p - N_{q+2n-i}^p, & i \in \{n+q+1, \dots, n\}, \end{cases} \quad (43) $$

where

$$ q = \lfloor \frac{p}{2} \rfloor - p = -\frac{p}{2}. $$

Taking into account that $p$ is even and $\lfloor \frac{p}{2} \rfloor = \frac{p}{2} = -q$, we observe a strong similarity between (43) and (57) in the proof of Theorem 17 (see Section A.2). This similarity implies that we can follow exactly the same proof, with the only difference that for the reduced spline space we need to replace the general form of the basis functions in (58) with

$$ \bar{N}_{i,0}^p = N_{q+i-1}^p - N_{q-k_i}^p. $$

Therefore, the general form of the elements of $\bar{X}^{p,r,0}$ is obtained by replacing the sign ‘+’ in the two terms of (59) that have only $k_i$ or $k_j$ by ‘-’. Then, we can simply follow the same steps as in the proof of Theorem 17, under the assumption $n \ge p + \frac{p}{2}$, and we arrive at the identity (42). $\square$

By combining Theorems 4 and 28 we directly deduce a closed-form expression for the eigenvalues and eigenvectors of the matrix $\bar{X}^{p,r,0}$.

**Corollary 29** Let $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $p$ even, $n \ge p + \frac{p}{2}$, and $r \le p$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $\bar{X}^{p,r,0}$ can be expressed as

$$ \lambda_j (\bar{X}^{p,r,0}) = n^{2r-1} g_p^r \left( \frac{j\pi}{n} \right), $$
and

$$u_{i,j}(\bar{X}^{p,r,0}) = \sqrt{\frac{2}{n}} c_j \sin\left(\frac{j\pi}{n}\left(i - \frac{1}{2}\right)\right), \quad c_j = \begin{cases} \frac{1}{\sqrt{2}}, & j = n, \\ 1, & j \le n-1, \end{cases}$$

with $i = 1, \dots, n$. The function $g_p^r$ is defined in (38).

A simplified version of Corollary 29 (considering $r = 0, 1$) reveals the expression for eigenvalues and eigenvectors of the discretized Laplace eigenvalue problem (1) with Dirichlet boundary conditions (2).

**Corollary 30** Let $p, n \in \mathbb{N}^*$ such that $p$ even and $n \ge p + \frac{p}{2}$. Then, for $j = 1, \dots, n$, the $j$-th eigenpair of the matrix $\bar{L}_n^{p,0} = (\bar{M}_n^{p,0})^{-1} \bar{K}_n^{p,0}$ can be expressed as

$$\lambda_j(\bar{L}_n^{p,0}) = n^2 \frac{g_p^1\left(\frac{j\pi}{n}\right)}{g_p^0\left(\frac{j\pi}{n}\right)},$$

and

$$u_{i,j}(\bar{L}_n^{p,0}) = \sqrt{\frac{2}{n}} c_j \sin\left(\frac{j\pi}{n}\left(i - \frac{1}{2}\right)\right), \quad c_j = \begin{cases} \frac{1}{\sqrt{2}}, & j = n, \\ 1, & j \le n-1, \end{cases}$$

with $i = 1, \dots, n$. The function $g_p^r$ with $r \in \{0, 1\}$ is defined in (38).

**Remark 31** The explicit closed-form expressions, described in Corollary 30, for the approximated spectrum of the Laplace operator with Dirichlet boundary conditions show that also the reduced spline spaces $\bar{S}_{p,n,0}$ ($p$ even) lead to outlier-free discretizations. This follows from (5), using the same line of arguments as in Remark 16. Indeed, by setting $\theta_j = \frac{j\pi}{n}$, $j = 1, \dots, n$, we get

$$\frac{\lambda_j(\bar{L}_n^{p,0}) - (j\pi)^2}{(j\pi)^2} = \frac{n^2}{(j\pi)^2}(e_p(\theta_j) - (\theta_j)^2) = \frac{e_p(\theta_j) - (\theta_j)^2}{(\theta_j)^2}.$$

Thus, for all $j = 1, \dots, n$, from (39) we deduce that the relative error between $\lambda_j(\bar{L}_n^{p,0})$ and the $j$-th eigenvalue of the Laplace operator with Dirichlet boundary conditions converges to zero as $p$ increases, resulting in an outlier-free approximation.

**Remark 32** From Corollaries 30 and 14 we see that, when considering the spaces $\bar{S}_{p,n+1,0}$ and $\bar{S}_{p,n,0}^{\text{opt}}$ with $p$ even, they both yield $n$ identical eigenvalues for the discretized Laplace operator – note that the former space has an additional $(n+1)$-st eigenvalue. However, the corresponding eigenvectors differ consistently. For the

<figure><img src="image_8.png" /></figure>
reduced spline space, there is a shift of $-1/2$ compared to the optimal spline subspace. This shift of $-1/2$ is not introduced in the construction of the basis of the reduced spline space (41) but is present in the basis of the optimal spline subspace (23); see in particular (24).

## 4.2 Multidimensional extension

Now, we expand the closed-form expressions for eigenvalues and eigenvectors from the 1-dimensional to the $d$-dimensional context. We start by introducing the $d$-dimensional eigenvalue problem and outline the approximate variational formulation. We focus on the tensor product of optimal spline subspaces suited for the Dirichlet boundary conditions (2). Similar results can be obtained for the other previously mentioned boundary conditions and spaces, including the reduced spline spaces, and we refer the reader to the general setting presented in [18].

Let $d \in \mathbb{N}^*$, we consider the following $d$-dimensional Laplace eigenvalue problem:

$$ \begin{cases} -\Delta u = \lambda u, & \text{in } (0, 1)^d, \\ u = 0, & \text{on } \partial((0, 1)^d). \end{cases} \quad (44) $$

The non-trivial exact solutions $(\lambda_k, u_k)$ of (44) are given by

$$ \lambda_k = \sum_{s=1}^{d} (k_s \pi)^2, \quad u_k(\mathbf{x}) = \prod_{s=1}^{d} \sin(k_s \pi x_s), $$

for $\mathbf{k} = (k_1, \dots, k_d) \in (\mathbb{N}^*)^d$ and $\mathbf{x} = (x_1, \dots, x_d) \in [0, 1]^d$.

Given $\mathbf{n} = (n_1, \dots, n_d) \in (\mathbb{N}^*)^d$ and $\mathbf{p} = (p_1, \dots, p_d) \in (\mathbb{N}^*)^d$, let us consider the tensor-product space $\bigotimes_{s=1}^d S_{p_s, n_s, 0}^{\text{opt}}$ in the context of the Galerkin method to discretize (44). From Section 2.2.1, it is clear that

$$ N_{i,0}^{\mathbf{p}} = N_{i_1,0}^{p_1} \otimes \dots \otimes N_{i_d,0}^{p_d}, \quad \mathbf{i} = \mathbf{1}, \dots, \mathbf{n}, $$

forms a basis for this space, where $\mathbf{i} = (i_1, \dots, i_d) \in (\mathbb{N}^*)^d$ is a multi-index varying between $\mathbf{1} = (1, \dots, 1)$ and $\mathbf{n} = (n_1, \dots, n_d)$. The approximate solutions of (44) can be computed as follows:

$$ \begin{align} \text{Find } (\lambda_{\mathbf{k},\mathbf{h}}, u_{\mathbf{k},\mathbf{h}}) \in \mathbb{R} & \times \bigotimes_{s=1}^{d} S_{p_s, n_s, 0}^{\text{opt}} \quad \text{such that} \tag{45} \\ K(u_{\mathbf{k},\mathbf{h}}, v) = \lambda_{\mathbf{k},\mathbf{h}} M(u_{\mathbf{k},\mathbf{h}}, v), \quad \forall v & \in \bigotimes_{s=1}^{d} S_{p_s, n_s, 0}^{\text{opt}} \notag \end{align} $$
where

$$M(u_{k,h}, v) = \int_{(0,1)^d} u_{k,h}(\mathbf{x})v(\mathbf{x}) \, d\mathbf{x},$$

$$K(u_{k,h}, v) = \int_{(0,1)^d} \nabla u_{k,h}(\mathbf{x}) \cdot \nabla v(\mathbf{x}) \, d\mathbf{x}.$$

The problem described by (45) can be reformulated as a discrete eigenvalue problem,

$$ (M_n^{p,0})^{-1} K_n^{p,0} u_{k,h} = \lambda_{k,h} u_{k,h}. \quad (46) $$

The matrices $M_n^{p,0}$ and $K_n^{p,0}$ are called the mass matrix and the stiffness matrix, respectively, and are defined by

$$M_n^{p,0} = \left[ \int_{(0,1)^d} N_{i,0}^p(\mathbf{x}) N_{j,0}^p(\mathbf{x}) \, d\mathbf{x} \right]_{i,j=1}^n,$$

$$K_n^{p,0} = \left[ \int_{(0,1)^d} \nabla N_{i,0}^p(\mathbf{x}) \cdot \nabla N_{j,0}^p(\mathbf{x}) \, d\mathbf{x} \right]_{i,j=1}^n.$$

Note that

$$M_n^{p,0} = \bigotimes_{s=1}^{d} M_{n_s}^{p_s,0}, \quad K_n^{p,0} = \sum_{r=1}^{d} \left( \bigotimes_{s=1}^{r-1} M_{n_s}^{p_s,0} \right) \otimes K_{n_r}^{p_r,0} \otimes \left( \bigotimes_{s=r+1}^{d} M_{n_s}^{p_s,0} \right), \quad (47)$$

and

$$L_n^{p,0} = (M_n^{p,0})^{-1} K_n^{p,0} = \sum_{r=1}^{d} \left( \bigotimes_{s=1}^{r-1} I_{n_s} \right) \otimes \left( (M_{n_r}^{p_r,0})^{-1} K_{n_r}^{p_r,0} \right) \otimes \left( \bigotimes_{s=r+1}^{d} I_{n_s} \right). \quad (48)$$

The matrices $K_{n_r}^{p_r,0}$ and $M_{n_r}^{p_r,0}$ for a fixed $r \in \{1, \dots, d\}$ are precisely those defined in (15) and studied in detail in Section 3.1. For more details on the above constructions and arguments of this section, we refer the reader to [18, Section 6.1].

In the next theorem we provide a closed-form expression for the eigenvalues and eigenvectors of the system (46), which is a direct consequence of Corollary 13.

**Theorem 33** Let $d \in \mathbb{N}^*$ and $\mathbf{p}, \mathbf{n} \in (\mathbb{N}^*)^d$ such that $n_s \ge \max\{p_s + 1, p_s + \lfloor \frac{p_s}{2} \rfloor - 1\}$ for all $s \in \{1, \dots, d\}$. Then, for $\mathbf{j} = \mathbf{1}, \dots, \mathbf{n}$, the $\mathbf{j}$-th eigenpairs of the matrices $M_n^{p,0}$, $L_n^{p,0}$, and $K_n^{p,0}$ can be expressed as

$$\lambda_j (M_n^{p,0}) = \prod_{s=1}^{d} (n_s + 1)^{-1} \left[ g_{p_s}^0 \left( \frac{j_s \pi}{n_s + 1} \right) \right],$$

<figure><img src="image_9.png" /></figure>
$$ \lambda_j (L_n^{p,0}) = \sum_{r=1}^{d} (n_r + 1)^2 \left[ \frac{g_{p_r}^1}{g_{p_r}^0} \left( \frac{j_r \pi}{n_r + 1} \right) \right], $$

$$ \lambda_j (K_n^{p,0}) = \lambda_j (M_n^{p,0}) \lambda_j (L_n^{p,0}), $$

and

$$ \mathbf{u}_j (M_n^{p,0}) = \mathbf{u}_j (L_n^{p,0}) = \mathbf{u}_j (K_n^{p,0}) = \bigotimes_{s=1}^{d} \left[ \sqrt{\frac{2}{n_s + 1}} \sin \left( \frac{ij_s \pi}{n_s + 1} \right) \right]_{i=1}^{n_s}. $$

The function $g_p^r$ with $r \in \{0, 1\}$ is defined in (38).

**Proof** Based on Corollary 13, for any fixed $s \in \{1, \dots, d\}$, we have

$$ M_{n_s}^{p,s,0} = Q_{n_s} \left( (n_s + 1)^{-1} \operatorname{diag} \left[ g_{p_s}^0 \left( \frac{j\pi}{n_s + 1} \right) \right] \right) Q_{n_s}, $$

$$ K_{n_s}^{p,s,0} = Q_{n_s} \left( (n_s + 1) \operatorname{diag} \left[ g_{p_s}^1 \left( \frac{j\pi}{n_s + 1} \right) \right] \right) Q_{n_s}, $$

where $Q_{n_s}$ is a symmetric and orthogonal matrix containing the eigenvectors, given by

$$ Q_{n_s} = \sqrt{\frac{2}{n_s + 1}} \left[ \sin \left( \frac{ij_s \pi}{n_s + 1} \right) \right]_{i,j=1}^{n_s}. $$

Then, from (47) and (48), we deduce

$$ M_n^{p,0} = \bigotimes_{s=1}^{d} M_{n_s}^{p,s,0} \\ = \left( \bigotimes_{s=1}^{d} Q_{n_s} \right) \left( \bigotimes_{s=1}^{d} \left( (n_s + 1)^{-1} \operatorname{diag} \left[ g_{p_s}^0 \left( \frac{j\pi}{n_s + 1} \right) \right] \right) \right) \left( \bigotimes_{s=1}^{d} Q_{n_s} \right), $$

and

$$ L_n^{p,0} = \sum_{r=1}^{d} \left( \bigotimes_{s=1}^{r-1} I_{n_s} \right) \otimes \left( \left( M_{n_r}^{p_r,0} \right)^{-1} K_{n_r}^{p_r,0} \right) \otimes \left( \bigotimes_{s=r+1}^{d} I_{n_s} \right) \\ = \sum_{r=1}^{d} \left( \bigotimes_{s=1}^{r-1} I_{n_s} \right) \otimes \left\{ Q_{n_r} \left( (n_r + 1)^2 \operatorname{diag} \left[ \frac{g_{p_r}^1}{g_{p_r}^0} \left( \frac{j\pi}{n_r + 1} \right) \right] \right) Q_{n_r} \right\} \\ \otimes \left( \bigotimes_{s=r+1}^{d} I_{n_s} \right) $$
$$ = \left( \bigotimes_{s=1}^{d} Q_{n_s} \right) \left\{ \sum_{r=1}^{d} \left( \bigotimes_{s=1}^{r-1} I_{n_s} \right) \otimes \left( (n_r + 1)^2 \operatorname{diag} \left[ \frac{g_{p_r}^1}{g_{p_r}^0} \left( \frac{j\pi}{n_r + 1} \right) \right] \right) \right\} \left( \bigotimes_{s=r+1}^{d} I_{n_s} \right) \left\{ \bigotimes_{s=1}^{d} Q_{n_s} \right\}. $$

Finally, we have $K_n^{p,0} = M_n^{p,0} L_n^{p,0}$, which completes the proof. □

## 5 Conclusions, further comments, and future work

We have provided explicit closed-form expressions for the eigenvalues and eigenvectors of the mass and stiffness matrices resulting from isogeometric Galerkin discretizations in optimal spline subspaces and specific reduced spline spaces that are outlier-free for the Laplace operator, under different types of boundary conditions. Our methodology is fully algebraic: it is based on proving that the mass and stiffness matrices possess a Toeplitz-minus-Hankel or Toeplitz-plus-Hankel structure, for any degree $p$. This specific perturbation of Toeplitz matrices results in the elimination of outliers and allows us to deduce a closed form for eigenvalues and eigenvectors. In all the studied cases, the eigenvalues of the mass and stiffness matrices are an explicitly known (uniform) sampling of the spectral symbol of the associated Toeplitz matrices.

The results in Sect. 3 are not confined to the mass and stiffness matrices related to isogeometric Galerkin discretizations of zeroth and second order problems in the considered spline subspaces (the cases $r = 0, 1$). They also hold for higher values of $r$, i.e., for isogeometric Galerkin discretizations of polyharmonic problems of order $r$ (with $p \ge 2r - 1$) under the so-called Laplace boundary conditions; see [33, Section 2.1]. For these problems, the spaces in (12) turn out to be optimal as well; see [33, Section 3.1]. Such boundary conditions, however, are of not much practical interest; they result in solutions that coincide with those for $r = 1$. Outlier-free discretizations for polyharmonic problems with more interesting boundary conditions have been proposed in [33] and an algebraic investigation of the spectral properties of the resulting matrices is a captivating topic for future research.

A natural question is to what extent the proposed outlier-free discretizations can be fruitfully used for addressing general problems with non-homogeneous boundary behavior. It is clear that a plain discretization in outlier-free spline subspaces generally leads to a substantial loss of approximation power compared to the corresponding full spline space because of the additional homogeneous boundary conditions that characterize the subspaces of interest. However, for problems identified by sufficiently smooth data, a suitable data-correction process for the missing boundary derivatives has been proposed in [32, Section 5], which is analogous to the classical reduction from non-homogeneous to homogeneous Dirichlet boundary conditions. When coupled with this boundary data correction, the discretization in outlier-free spline subspaces achieves full approximation order both in the univariate and in the multivariate tensor-product case. It is noteworthy that the presented spectral results are highly valuable
for the development of powerful solvers in that context. For more details and insights, we direct readers to [18, Remark 1] and references therein.

Lastly, an important aspect is that the machinery used here for the outlier-free theory could extend to the case of coercive second order operators with variable coefficients. The key will be the nature of the linear positive operator (see [37], [24, Corollary 6.2, Section 6.3], [25, Corollary 3.2, Section 3.3] and references therein) of the underlying matrices and we are convinced that designing outlier-free approximating matrices in a more general setting is worth to be investigated in future research.

## A Proofs of main theorems

To simplify the notation, we define the parameters

$$q = \lfloor \frac{p}{2} \rfloor - p, \quad p_0 = p - 2 \lfloor \frac{p}{2} \rfloor,$$

which help in unifying the different basis constructions for even and odd degree $p$. Moreover, in the case $p = 1$, any matrix block with $\lfloor \frac{p}{2} \rfloor$ rows or columns is considered to be non-existent.

### A.1 Proof of Theorem 11

Fix $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $n \ge \max\{p+1, p+\lfloor \frac{p}{2} \rfloor - 1\}$ and $r \le p$. From the construction of the basis in (21) and (23), and using the assumption $n \ge p+1$, we deduce

$$ \begin{bmatrix} N_{1,0}^p \\ N_{2,0}^p \\ \vdots \\ N_{n,0}^p \end{bmatrix} = \left[ \underbrace{\left( \begin{array}{c} \lfloor \frac{p}{2} \rfloor + 1 \\ L_{(n, \lfloor \frac{p}{2} \rfloor + 1)}^0 \end{array} \right)}_{n+p+2-p_0} \quad \underbrace{\left( \begin{array}{c} \lfloor \frac{p}{2} \rfloor + 1 \\ R_{(n, \lfloor \frac{p}{2} \rfloor + 1)}^0 \end{array} \right)}_{n+p+2-p_0} \right] \begin{bmatrix} N_{-p}^p \\ \vdots \\ N_0^p \\ N_1^p \\ \vdots \\ N_{n+1-p_0}^p \end{bmatrix}, $$

where

$$ L_{(n, \lfloor \frac{p}{2} \rfloor + 1)}^0 = \begin{bmatrix} -J_{\lfloor \frac{p}{2} \rfloor} & O_{(\lfloor \frac{p}{2} \rfloor, 1)} \\ O_{(n-\lfloor \frac{p}{2} \rfloor, \lfloor \frac{p}{2} \rfloor)} & O_{(n-\lfloor \frac{p}{2} \rfloor, 1)} \end{bmatrix}, $$

$$ R_{(n, \lfloor \frac{p}{2} \rfloor + 1)}^0 = \begin{bmatrix} O_{(n-\lfloor \frac{p}{2} \rfloor, 1)} & O_{(n-\lfloor \frac{p}{2} \rfloor, \lfloor \frac{p}{2} \rfloor)} \\ O_{(\lfloor \frac{p}{2} \rfloor, 1)} & -J_{\lfloor \frac{p}{2} \rfloor} \end{bmatrix}. $$

We recall that $I_k$ is the identity matrix of size $k$, $J_k$ the exchange matrix of size $k$, and $O_{(k,l)}$ the zero matrix of the indicated size. Through a straightforward computation,
we arrive at

$$N_{i,0}^p = \begin{cases} -N_{q-i}^p + N_{q+i}^p, & i \in \{1, \dots, \lfloor \frac{p}{2} \rfloor\}, \\ N_{q+i}^p, & i \in \{\lfloor \frac{p}{2} \rfloor + 1, \dots, n - \lfloor \frac{p}{2} \rfloor\}, \\ N_{q+i}^p - N_{q+2(n+1)-i}^p, & i \in \{n+1 - \lfloor \frac{p}{2} \rfloor, \dots, n\}. \end{cases} \quad (49)$$

When we assume their restriction to the interval $[0, 1]$, the basis functions in (49) can be written in the following general form:

$$N_{i,0}^p = N_{q+i}^p - N_{q-k_i}^p, \quad i = 1, \dots, n, \quad (50)$$

where $k_i$ represents a certain index depending on the value of $i \in \{1, \dots, n\}$. Here we exploit the fact that $N_i^p$ has no support on $[0, 1]$ for very large positive or negative values of index $l$; see (18).

Now, we proceed to compute the elements of the matrix $X^{p,r,0}$; see Definition 7 with $b=0$. For fixed $i, j \in \{1, \dots, n\}$, we have

$$X_{i,j}^{p,r,0} = \int_0^1 (N_{i,0}^p)^{(r)}(x)(N_{j,0}^p)^{(r)}(x) \, dx \\ = \int_0^1 \left[ (N_{q+i}^p)^{(r)}(x) - (N_{q-k_i}^p)^{(r)}(x) \right] \left[ (N_{q+j}^p)^{(r)}(x) - (N_{q-k_j}^p)^{(r)}(x) \right] \, dx \\ = \int_0^1 (N_{q+i}^p)^{(r)}(x)(N_{q+j}^p)^{(r)}(x) \, dx + \int_0^1 (N_{q-k_i}^p)^{(r)}(x)(N_{q-k_j}^p)^{(r)}(x) \, dx \\ - \int_0^1 (N_{q+i}^p)^{(r)}(x)(N_{q-k_j}^p)^{(r)}(x) \, dx - \int_0^1 (N_{q-k_i}^p)^{(r)}(x)(N_{q+j}^p)^{(r)}(x) \, dx,$$

and using the definitions of $N_i^p$ in (22) and (24), we find

$$X_{i,j}^{p,r,0} = (n+1)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{n+1-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + (n+1)^{2r-1} \int_{k_j-q+\frac{1-p_0}{2}}^{n+1+k_j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+k_i-k_j) \, dt \\ - (n+1)^{2r-1} \int_{k_j-q+\frac{1-p_0}{2}}^{n+1+k_j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-k_j) \, dt \\ - (n+1)^{2r-1} \int_{k_i-q+\frac{1-p_0}{2}}^{n+1+k_i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-k_i-j) \, dt. \quad (51)$$

As stipulated earlier, the indices $k_i$ and $k_j$ depend on the value of $i$ and $j$, respectively. Thus, to compute $X^{p,r,0}$ it is imperative to address all the possible cases in (49). From Lemma 9 we know that $X^{p,r,0}$ is symmetric and centrosymmetric. This observation reduces the proof of (36) to two main steps.
**Step 1** We show that the block

$$ (X_{i,j}^{p,r,0})_{\lfloor \frac{p}{2} \rfloor + 1 \le i, j \le n - \lfloor \frac{p}{2} \rfloor} $$

is a Toeplitz matrix. Let us fix $i, j \in \{\lfloor \frac{p}{2} \rfloor + 1, \dots, n - \lfloor \frac{p}{2} \rfloor\}$. In this case, the expression of the matrix element in (51) simplifies to

$$ X_{i,j}^{p,r,0} = (n+1)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{n+1-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt. $$

By observing that

$$ n+1-i-q+\frac{1-p_0}{2} \ge p+1, \quad -i-q+\frac{1-p_0}{2} \le 0, $$

and by using (18), (20), we obtain

$$ \begin{aligned} X_{i,j}^{p,r,0} &= (n+1)^{2r-1} \int_{0}^{p+1} \mathcal{N}_{p}^{(r)}(t) \mathcal{N}_{p}^{(r)}(t+i-j) \, dt \\ &= (n+1)^{2r-1} (-1)^{r} \mathcal{N}_{2p+1}^{(2r)}(p+1+i-j). \end{aligned} $$

This concludes the proof of the first step.

**Step 2** We show that the block

$$ (X_{i,j}^{p,r,0})_{\substack{1 \le i \le \lfloor \frac{p}{2} \rfloor \\ i \le j \le n}} $$

has a Toeplitz-minus-Hankel form, where the Hankel part starts with the third element of the associated Toeplitz matrix. More precisely, for the corresponding indices $(i, j)$ we have

$$ X_{i,j}^{p,r,0} = (n+1)^{2r-1} \left( (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i-j) - (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i+j) \right). \quad (52) $$

Let us fix $i \in \{1, \dots, \lfloor \frac{p}{2} \rfloor\}$. In the following, we discuss the three different cases for $j$ arising from (49).

If $j \in \{i, \dots, \lfloor \frac{p}{2} \rfloor\}$, then we have $k_i = i$ and $k_j = j$ in (50). In this case, the expression of the matrix element in (51) reads as
$$X_{i,j}^{p,r,0} = (n+1)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{n+1-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + (n+1)^{2r-1} \int_{j-q+\frac{1-p_0}{2}}^{n+1+j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ - (n+1)^{2r-1} \int_{j-q+\frac{1-p_0}{2}}^{n+1+j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j) \, dt \\ - (n+1)^{2r-1} \int_{i-q+\frac{1-p_0}{2}}^{n+1+i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j) \, dt. \tag{53}$$

Since $n \ge p+1$, we have

$$n+1 \pm i-q+\frac{1-p_0}{2} \ge p+1, \quad n+1+j-q+\frac{1-p_0}{2} \ge p+1.$$

Together with (18), this implies that the equality in (53) is still valid if the upper bounds of the four integration intervals in (53) are replaced by $p+1$. By leveraging Lemma 8 and the identity

$$p+1+q-\frac{1-p_0}{2}=-q+\frac{1-p_0}{2},$$

we get

$$\int_{j-q+\frac{1-p_0}{2}}^{p+1} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt = \int_{-i+j}^{-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt, \tag{54}$$

and

$$\int_{i-q+\frac{1-p_0}{2}}^{p+1} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j) \, dt = \int_{i+j}^{j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j) \, dt. \tag{55}$$

By substituting (54)-(55) into (53) and by using (18), (20), we arrive at (52).

Now, let us consider $j \in \{\lfloor \frac{p}{2} \rfloor + 1, \dots, n - \lfloor \frac{p}{2} \rfloor\}$. In this case, the expression of the matrix element in (51) reads as

$$X_{i,j}^{p,r,0} = (n+1)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{n+1-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ - (n+1)^{2r-1} \int_{i-q+\frac{1-p_0}{2}}^{n+1+i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j) \, dt.$$

<figure><img src="image_10.png" /></figure>
By following the same reasoning as in the first case and by utilizing the implication

$$t \le \pm i - q + \frac{1 - p_0}{2} \implies t \ne i - j \le 0,$$

we obtain (52).

Lastly, we take $j$ in the range $\{n+1 - \lfloor \frac{p}{2} \rfloor, \dots, n\}$. In this case, the expression of the matrix element in (51) reads as

$$X_{i,j}^{p,r,0} = (n+1)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{n+1-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + (n+1)^{2r-1} \int_{-2(n+1)+j-q+\frac{1-p_0}{2}}^{-(n+1)+j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j+2(n+1)) \, dt \\ - (n+1)^{2r-1} \int_{-2(n+1)+j-q+\frac{1-p_0}{2}}^{-(n+1)+j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j+2(n+1)) \, dt \\ - (n+1)^{2r-1} \int_{i-q+\frac{1-p_0}{2}}^{n+1+i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j) \, dt. \quad (56)$$

Since $n \ge p+1$, we find

$$n + 1 \pm i - q + \frac{1 - p_0}{2} \ge p + 1, \quad -2(n + 1) + j - q + \frac{1 - p_0}{2} \le 0.$$

Together with (18), this implies that the equality in (56) is still valid if the upper bounds of the first and the last integration interval in (56) are replaced by $p+1$, and the lower bounds of the second and the third integration interval in (56) are replaced by zero. Furthermore, from $n \ge p + \lfloor \frac{p}{2} \rfloor - 1$ we deduce

$$t \ge 0 \implies t \pm i - j + 2(n + 1) \ge p + 1,$$

and thus the second and the third integral in (56) are both equal zero. Then, by following the same reasoning as in the previous two cases, we arrive at (52) again. This concludes the proof of the second step.

Finally, regrouping both steps and using the symmetry and centrosymmetry of $X^{p,r,0}$, we arrive at (36), and hence the proof of the theorem is complete.

## A.2 Proof of Theorem 17

The proof closely mirrors the one of Theorem 11 (see the previous section). Let us fix $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $r \le p$ and

$$n \ge \max \left\{ 2p - \lfloor \frac{p}{2} \rfloor, 2p - 2 \lfloor \frac{p}{2} \rfloor + 1 \right\} = \max \left\{ p + \lfloor \frac{p}{2} \rfloor + p_0, p + 1 + p_0 \right\}.$$
From the construction of the basis in (25) and (27), and using the assumption $n \ge p + 1 + p_0$, we have

$$ \begin{bmatrix} N_{1,1}^p \\ N_{2,1}^p \\ \vdots \\ N_{n,1}^p \end{bmatrix} = \left[ \begin{array}{c|c} \underbrace{\overbrace{J_n}^{\lfloor \frac{p}{2} \rfloor + p_0}}_{n+p+p_0} & \underbrace{\overbrace{I_n}^{\lfloor \frac{p}{2} \rfloor + p_0}}_{n+p+p_0} \\ \hline \end{array} \right] \begin{bmatrix} N_{-p}^p \\ \vdots \\ N_0^p \\ N_1^p \\ \vdots \\ N_{n+p_0-1}^p \end{bmatrix}. $$

Through a straightforward computation, we obtain

$$ N_{i,1}^p = \begin{cases} N_{-\lfloor \frac{p}{2} \rfloor - i}^p + N_{-\lfloor \frac{p}{2} \rfloor + i - 1}^p, & i \in \{1, \dots, -q\}, \\ N_{-\lfloor \frac{p}{2} \rfloor + i - 1}^p, & i \in \{-q + 1, \dots, n + q\}, \\ N_{-\lfloor \frac{p}{2} \rfloor + i - 1}^p + N_{-\lfloor \frac{p}{2} \rfloor + 2n - i}^p, & i \in \{n + q + 1, \dots, n\}. \end{cases} \quad (57) $$

When we assume their restriction to the interval $[0, 1]$, the basis functions in (57) can be written in the following general form:

$$ N_{i,1}^p = N_{-\lfloor \frac{p}{2} \rfloor + i - 1}^p + N_{-\lfloor \frac{p}{2} \rfloor - k_i}^p, \quad (58) $$

where $k_i$ represents again a certain index depending on the value of $i \in \{1, \dots, n\}$.

Now, we look into the elements of the matrix $X^{p,r,1}$; see Definition 7 with $b = 1$. Similar to the derivation of (51) and taking into account the definitions of $N_i^p$ in (26) and (28), we find

$$ \begin{aligned} X_{i,j}^{p,r,1} = {}& n^{2r-1} \int_{-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt \\ & + n^{2r-1} \int_{k_j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n+k_j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+k_i-k_j) \, dt \\ & + n^{2r-1} \int_{k_j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n+k_j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-(i-1+k_j)) \, dt \\ & + n^{2r-1} \int_{k_i+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n+k_i+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-(k_i+j-1)) \, dt. \end{aligned} \quad (59) $$

A key observation in the proof of Theorem 11 was the fact that the matrix $X^{p,r,0}$ is symmetric and centrosymmetric (Lemma 9). In the same way, one can show that the matrix $X^{p,r,1}$ is also symmetric and centrosymmetric. This reduces the proof of (40) to two main steps.
**Step 1** We show that the block

$$ (X_{i,j}^{p,r,1})_{-q+1 \le i, j \le n+q} $$

is a Toeplitz matrix. Let $i, j \in \{-q+1, \dots, n+q\}$. Then, we can eliminate all terms corresponding to $k_i$ and $k_j$ in (59), i.e.,

$$ X_{i,j}^{p,r,1} = n^{2r-1} \int_{-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt. $$

By observing that

$$ n - i + 1 + \lfloor \frac{p}{2} \rfloor + \frac{p_0}{2} \ge p + 1, \quad -i + 1 + \lfloor \frac{p}{2} \rfloor + \frac{p_0}{2} \le 0, $$

and by using (18), (20), we obtain

$$ X_{i,j}^{p,r,1} = n^{2r-1} \int_{0}^{p+1} \mathcal{N}_{p}^{(r)}(t) \mathcal{N}_{p}^{(r)}(t+i-j) \, dt = n^{2r-1} (-1)^{r} \mathcal{N}_{2p+1}^{(2r)}(p+1+i-j). $$

This concludes the proof of the first step.

**Step 2** We show that the block

$$ (X_{i,j}^{p,r,1})_{\substack{1 \le i \le -q \\ i \le j \le n}} $$

has a Toeplitz-plus-Hankel form, where the Hankel part starts with the second element of the associated Toeplitz matrix. More precisely, for the corresponding indices $(i, j)$ we have

$$ X_{i,j}^{p,r,1} = n^{2r-1} \left( (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i-j) + (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i+j-1) \right). \quad (60) $$

Let us fix $i \in \{1, \dots, -q\}$. In the following we discuss the three different cases for $j$ arising from (57).
In the first case, where we consider $j \in \{i, \dots, -q\}$, we have $k_i = i$ and $k_j = j$ in (58). Hence, the expression of the matrix element in (59) reads as

$$X_{i,j}^{p,r,1} = n^{2r-1} \int_{-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + n^{2r-1} \int_{j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n+j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + n^{2r-1} \int_{j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n+j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-i-j+1) \, dt \\ + n^{2r-1} \int_{i+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n+i+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-i-j+1) \, dt. \tag{61}$$

Since $n \ge p+1+p_0$, we have

$$n \pm i + \lfloor \frac{p}{2} \rfloor + \frac{p_0}{2} \ge p+1, \quad n+j + \lfloor \frac{p}{2} \rfloor + \frac{p_0}{2} \ge p+1.$$

Together with (18), this implies that the equality in (61) is still valid if the upper bounds of the four integration intervals in (61) are replaced by $p+1$. By exploiting Lemma 8 and the identity

$$p - \lfloor \frac{p}{2} \rfloor - \frac{p_0}{2} = \lfloor \frac{p}{2} \rfloor + \frac{p_0}{2},$$

we get

$$\int_{j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{p+1} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt \\ = \int_{-i+j}^{-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt, \tag{62}$$

and

$$\int_{i+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{p+1} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-i-j+1) \, dt \\ = \int_{i+j-1}^{j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-i-j+1) \, dt. \tag{63}$$

By substituting (62)–(63) into (61) and by using (18), (20), we infer (60).
Now, we fix $j \in \{-q+1, \dots, n+q\}$. In this case, the expression of the matrix element in (59) reads as

$$X_{i,j}^{p,r,1} = n^{2r-1} \int_{-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + n^{2r-1} \int_{i+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n+i+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-i-j+1) \, dt.$$

By following the same reasoning as in the first case and by utilizing the implication

$$t \le \pm i + \lfloor \frac{p}{2} \rfloor + \frac{p_0}{2} \implies t \ne i - j + 1 \le 0,$$

we obtain (60).

In the last case, where we consider $j \in \{n+q+1, \dots, n\}$, the expression of the matrix element in (59) reads as

$$X_{i,j}^{p,r,1} = n^{2r-1} \int_{-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n-i+1+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + n^{2r-1} \int_{j-2n+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{-n+j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j+2n) \, dt \\ + n^{2r-1} \int_{j-2n+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{-n+j+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-i-j+2n+1) \, dt \\ + n^{2r-1} \int_{i+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}}^{n+i+\lfloor \frac{p}{2} \rfloor + \frac{p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-i-j+1) \, dt. \tag{64}$$

Since $n \ge p+1+p_0$, we find

$$n \pm i + \lfloor \frac{p}{2} \rfloor + \frac{p_0}{2} \ge p+1, \quad j-2n + \lfloor \frac{p}{2} \rfloor + \frac{p_0}{2} \le 0.$$

Together with (18), this implies that the equality in (64) is still valid if the upper bounds of the first and the last integration interval in (64) are replaced by $p+1$, and the lower bounds of the second and the third integration interval in (64) are replaced by zero. Furthermore, from $n \ge p + \lfloor \frac{p}{2} \rfloor + p_0$ we deduce

$$t \ge 0 \implies \{t - i - j + 2n + 1, t + i - j + 2n\} \ge p + 1,$$

and thus the second and the third integral in (64) are both equal zero. Then, by following the same reasoning as in the previous two cases, we arrive at (60) again. The latter concludes the proof of the second step and consequently the proof of the theorem is complete.
### A.3 Proof of Theorem 23

Also here, the proof is similar to the one of Theorem 11 (see Section A.1). Let us fix $p, n \in \mathbb{N}^*$ and $r \in \mathbb{N}$ such that $n \ge \max\{p+1, p+\lfloor\frac{p}{2}\rfloor\}$ and $r \le p$. From the construction of the basis in (29) and (31), and using the assumption $n \ge p+1$, we have

$$ \begin{bmatrix} N_{1,2}^p \\ N_{2,2}^p \\ \vdots \\ N_{n,2}^p \end{bmatrix} = \left[ \overbrace{\left( \begin{bmatrix} -J_n & O_{(n,1)} & I_n \end{bmatrix} \middle| \begin{bmatrix} \lfloor \frac{p}{2} \rfloor + p_0 \\ \lfloor \frac{p}{2} \rfloor \\ \vdots \\ \lfloor \frac{p}{2} \rfloor \end{bmatrix} \right)}^{\lfloor \frac{p}{2} \rfloor} \right] \begin{bmatrix} N_{-p}^p \\ N_0^p \\ N_1^p \\ \vdots \\ N_n^p \end{bmatrix}, $$

and we can express the basis functions in the following general form:

$$ N_{i,2}^p = N_{q+i}^p + \gamma_i N_{q-k_i}^p, \quad i = 1, \dots, n, $$

where

$$ (\gamma_i, k_i) = \begin{cases} (-1, i), & i \in \{1, \dots, \lfloor \frac{p}{2} \rfloor\}, \\ (0, \cdot), & i \in \{\lfloor \frac{p}{2} \rfloor + 1, \dots, n+q\}, \\ (1, i - (2n+1)), & i \in \{n+q+1, \dots, n\}. \end{cases} \quad (65) $$

Now, we look into the elements of the matrix $X^{p,r,2}$; see Definition 7 with $b=2$. Similar to the derivation of (51) and taking into account the definitions of $N_i^p$ in (30) and (32), we find

$$ \begin{aligned} X_{i,j}^{p,r,2} = & \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt \\ & + \gamma_i \gamma_j \left(\frac{2n+1}{2}\right)^{2r-1} \int_{k_j-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}+k_j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+k_i-k_j) \, dt \\ & + \gamma_j \left(\frac{2n+1}{2}\right)^{2r-1} \int_{k_j-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}+k_j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-i-k_j) \, dt \\ & + \gamma_i \left(\frac{2n+1}{2}\right)^{2r-1} \int_{k_i-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}+k_i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-k_i-j) \, dt. \end{aligned} \quad (66) $$

In the preceding proofs, we relied on both symmetry and centrosymmetry of the matrices $X^{p,r,0}$ and $X^{p,r,1}$. However, the matrix $X^{p,r,2}$ lacks centrosymmetry and thus, we can only leverage the symmetry of the matrix to proceed. This introduces an additional step.
**Step 1** We show that the block

$$ (X_{i,j}^{p,r,2})_{\lfloor \frac{p}{2} \rfloor + 1 \le i, j \le n+q} $$

is a Toeplitz matrix. Let us fix $i, j \in \{\lfloor \frac{p}{2} \rfloor + 1, \dots, n+q\}$. In this case, we have $\gamma_i = \gamma_j = 0$ and the expression of the matrix element in (66) simplifies to

$$ X_{i,j}^{p,r,2} = \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{\frac{2n+1-i-q}{2}-\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt. $$

By observing that

$$ \frac{2n+1}{2} - i - q + \frac{1-p_0}{2} \ge p+1, \quad -i - q + \frac{1-p_0}{2} \le 0, $$

and by using (18), (20), we deduce

$$ \begin{aligned} X_{i,j}^{p,r,2} &= \left(\frac{2n+1}{2}\right)^{2r-1} \int_{0}^{p+1} \mathcal{N}_{p}^{(r)}(t) \mathcal{N}_{p}^{(r)}(t+i-j) \, dt \\ &= \left(\frac{2n+1}{2}\right)^{2r-1} (-1)^{r} \mathcal{N}_{2p+1}^{(2r)}(p+1+i-j). \end{aligned} $$

This concludes the proof of the first step.

**Step 2** We show that the block

$$ (X_{i,j}^{p,r,2})_{\substack{1 \le i \le \lfloor \frac{p}{2} \rfloor \\ i \le j \le n}} $$

has a Toeplitz-minus-Hankel form, where the Hankel part starts with the third element of the associated Toeplitz matrix. This form is the same as the Dirichlet case. More precisely, for the corresponding indices $(i, j)$ we have

$$ X_{i,j}^{p,r,2} = \left(\frac{2n+1}{2}\right)^{2r-1} \qquad (67) \\ \left((-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i-j) - (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i+j)\right). $$

Let us fix $i \in \{1, \dots, \lfloor \frac{p}{2} \rfloor\}$. In the following we discuss the three different cases for $j$ arising from (65).
Firstly, let us consider $j \in \{i, \dots, \lfloor \frac{p}{2} \rfloor\}$. In this case, we have $\gamma_i = \gamma_j = -1$, $k_i = i$, and $k_j = j$. Then, the expression of the matrix element in (66) reads as

$$X_{i,j}^{p,r,2} = \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + \left(\frac{2n+1}{2}\right)^{2r-1} \int_{j-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}+j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ - \left(\frac{2n+1}{2}\right)^{2r-1} \int_{j-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}+j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j) \, dt \\ - \left(\frac{2n+1}{2}\right)^{2r-1} \int_{i-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}+i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j) \, dt. \tag{68}$$

Since $n \ge p+1$, we have

$$\frac{2n+1}{2} \pm i - q + \frac{1-p_0}{2} \ge p+1, \quad \frac{2n+1}{2} + j - q + \frac{1-p_0}{2} \ge p+1.$$

Together with (18), this implies that the equality in (68) is still valid if the upper bounds of the four integration intervals in (68) are replaced by $p+1$. By using Lemma 8 and the identity

$$p+1+q-\frac{1-p_0}{2}=-q+\frac{1-p_0}{2}, \tag{69}$$

we get

$$\int_{j-q+\frac{1-p_0}{2}}^{p+1} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt = \int_{-i+j}^{-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt, \tag{70}$$

and

$$\int_{i-q+\frac{1-p_0}{2}}^{p+1} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-(i+j)) \, dt = \int_{i+j}^{j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-i-j) \, dt. \tag{71}$$

By substituting (70)–(71) into (68) and by using (18), (20), we arrive at (67).
Now, let us fix $j \in \{\lfloor \frac{p}{2} \rfloor + 1, \dots, n+q\}$ and so $\gamma_j = 0$. In this case, the expression of the matrix element in (66) reads as

$$X_{i,j}^{p,r,2} = \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ - \left(\frac{2n+1}{2}\right)^{2r-1} \int_{i-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}+i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j) \, dt.$$

By following the same reasoning as in the first case and by using the implication

$$t \le \pm i - q + \frac{1 - p_0}{2} \implies t \neq i - j \le 0,$$

we obtain (67).

The last case addresses $j \in \{n+q+1, \dots, n\}$, where we have $\gamma_j = 1$ and $k_j = j - (2n+1)$. Then, the expression of the matrix element in (66) reads as

$$X_{i,j}^{p,r,2} = \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ - \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-(2n+1)+j-q+\frac{1-p_0}{2}}^{-\frac{2n+1}{2}+j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j+2n+1) \, dt \\ + \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-(2n+1)+j-q+\frac{1-p_0}{2}}^{-\frac{2n+1}{2}+j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j+2n+1) \, dt \\ - \left(\frac{2n+1}{2}\right)^{2r-1} \int_{i-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}+i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j) \, dt. \quad (72)$$

Since $n \ge p+1$, we find

$$\frac{2n+1}{2} \pm i - q + \frac{1 - p_0}{2} \ge p+1, \quad -(2n+1) + j - q + \frac{1 - p_0}{2} \le 0.$$

Together with (18), this implies that the equality in (72) is still valid if the upper bounds of the first and the last integration interval in (72) are replaced by $p+1$, and the lower bounds of the second and the third integration interval in (72) are replaced by zero. Furthermore, from $n \ge p + \lfloor \frac{p}{2} \rfloor$ we deduce

$$t \ge 0 \implies t \pm i - j + 2n + 1 \ge p + 1,$$

and thus the second and the third integral in (72) are both equal zero. Then, by following the same reasoning as in the previous two cases, we infer (67). This concludes the proof of the second step.

Finally, we observe that the assumption $n \ge p + \lfloor \frac{p}{2} \rfloor$ employed at the end of the last case of the second step also implies that for $j \ge n+q+1$, the Hankel part is
zero, i.e.,

$$(-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i+j) = 0, \quad i=1, \dots, \lfloor \frac{p}{2} \rfloor, \quad j=n+q+1, \dots, n. \tag{73}$$

**Step 3** This additional step is necessary by the absence of the centrosymmetry property of $X_{i,j}^{p,r,2}$. We show that the block

$$(X_{i,j}^{p,r,2})_{\substack{n+1+q \le i \le n \\ 1 \le j \le n}}$$

has a Toeplitz-plus-Hankel form, where the Hankel part starts with the second element of the associated Toeplitz matrix. This form is the same as the Neumann case. More precisely, for the corresponding indices $(i, j)$ we have

$$X_{i,j}^{p,r,2} = \left(\frac{2n+1}{2}\right)^{2r-1} \left( (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i-j) + (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1-i-j+2n+1) \right). \tag{74}$$

As before, we fix $i \in \{n+q+1, \dots, n\}$ and we discuss the three different cases for $j$.

Firstly, we consider $j \in \{n+q+1, \dots, n\}$. In this case, we have $\gamma_i = \gamma_j = 1$, $k_i = i - (2n+1)$, and $k_j = j - (2n+1)$. Then, the expression of the matrix element in (66) reads as

$$X_{i,j}^{p,r,2} = \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-(2n+1)+j-q+\frac{1-p_0}{2}}^{-\frac{2n+1}{2}+j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-(2n+1)+j-q+\frac{1-p_0}{2}}^{-\frac{2n+1}{2}+j-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j+2n+1) \, dt \\ + \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-(2n+1)+i-q+\frac{1-p_0}{2}}^{-\frac{2n+1}{2}+i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t)\mathcal{N}_p^{(r)}(t-i-j+2n+1) \, dt. \tag{75}$$

Taking into account $n \ge p+1$, one can check that

$$\begin{aligned} & -i - q + \frac{1 - p_0}{2} \le 0, \quad -(2n+1) + j - q + \frac{1 - p_0}{2} \le 0, \\ & -(2n+1) + i - q + \frac{1 - p_0}{2} \le 0. \end{aligned}$$
Together with (18), this implies that the equality in (75) is still valid if the lower bounds of the four integration intervals in (75) are replaced by zero. By leveraging Lemma 8 and the identity (69), we get

$$ 
\begin{aligned}
& \int_{0}^{-\frac{2n+1}{2} + j - q + \frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t + i - j) \, dt \\
& = \int_{\frac{2n+1}{2} - i - q + \frac{1-p_0}{2}}^{p+1-i+j} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t + i - j) \, dt,
\end{aligned} 
\quad (76) $$

and

$$ 
\begin{aligned}
& \int_{0}^{-\frac{2n+1}{2} + i - q + \frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t - i - j + 2n + 1) \, dt \\
& = \int_{-\frac{2n+1}{2} + j - q + \frac{1-p_0}{2}}^{p+1-(2n+1)+i+j} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t - i - j + 2n + 1) \, dt.
\end{aligned} 
\quad (77) $$

By substituting equations (76)–(77) into (75) and by using (18), (20), we infer (74).

Now, let us consider the case $j \in \{\lfloor \frac{p}{2} \rfloor + 1, \dots, n+q\}$, where $\gamma_j = 0$. Then, the expression of the matrix element in (66) reads as

$$ X_{i,j}^{p,r,2} = \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-i-q+\frac{1-p_0}{2}}^{\frac{2n+1}{2}-i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t+i-j) \, dt \\ + \left(\frac{2n+1}{2}\right)^{2r-1} \int_{-(2n+1)+i-q+\frac{1-p_0}{2}}^{-\frac{2n+1}{2}+i-q+\frac{1-p_0}{2}} \mathcal{N}_p^{(r)}(t) \mathcal{N}_p^{(r)}(t-i-j+2n+1) \, dt. $$

By following the same reasoning as in the first case and by using the implication

$$ t \ge \frac{2n+1}{2} \pm i - q + \frac{1-p_0}{2} \implies t \ne i - j \ge p + 1, $$

we obtain (74).

For the last case, where $j \in \{1, \dots, \lfloor \frac{p}{2} \rfloor\}$, we exploit the symmetry of the matrix, the result of Step 2, and the observation in (73), to deduce that

$$ X_{i,j}^{p,r,2} = \left(\frac{2n+1}{2}\right)^{2r-1} (-1)^r \mathcal{N}_{2p+1}^{(2r)}(p+1+i-j). $$

Note that, since $n \ge p + \lfloor \frac{p}{2} \rfloor$, this expression is equal to (74) for $j \le \lfloor \frac{p}{2} \rfloor$. This concludes the proof of the third step and consequently the proof of the theorem is complete.
## B Matrix algebras

In this appendix, we give a compact overview of $\tau$ matrix algebras, focusing specifically on those encountered in our analysis. Additional algebras can be found in [5] and [14].

Following [5] (see also [18]), we consider the following class of real symmetric matrices

$$ T_n(\varepsilon, \varphi) = \begin{bmatrix} \varepsilon & 1 & 0 & \cdots & 0 \\ 1 & 0 & \ddots & \ddots & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ \vdots & \ddots & \ddots & 0 & 1 \\ 0 & \cdots & 0 & 1 & \varphi \end{bmatrix}, \quad \varepsilon, \varphi \in \{-1, 0, 1\}. $$

Then, the matrix algebras generated over $\mathbb{R}$ by $T_n(\varepsilon, \varphi)$ are defined as

$$ \tau_n(\varepsilon, \varphi) = \left\{ Q_n(\varepsilon, \varphi) \left( \operatorname{diag} [a_j] \right) Q_n^T(\varepsilon, \varphi) : (a_j)_{1 \le j \le n} \in \mathbb{R}^n \right\}, $$

where $Q_n(\varepsilon, \varphi)$ represents an orthogonal matrix determined through the diagonalization of $T_n(\varepsilon, \varphi)$:

$$ T_n(\varepsilon, \varphi) = Q_n(\varepsilon, \varphi) \Lambda(\varepsilon, \varphi) Q_n^T(\varepsilon, \varphi). $$

It can be observed from our closed-form expressions of eigenvalues and eigenvectors (see Sects. 3 and 4.1) that our mass and stiffness matrices fall within one of the following algebras.

* Consider $\varepsilon = \varphi = 0$. Then,

$$ Q_n(0, 0) = \sqrt{\frac{2}{n+1}} \left[ \sin \left( \frac{ij\pi}{n+1} \right) \right]_{i,j=1}^n, $$

$$ \Lambda(0, 0) = \operatorname{diag}_{j=1,...,n} \left[ 2 \cos \left( \frac{j\pi}{n+1} \right) \right]. $$

In this particular case, the matrix $Q_n(0, 0)$ is also symmetric.

* Consider $\varepsilon = \varphi = 1$. Then,

$$ Q_n(1, 1) = \sqrt{\frac{2}{n}} \left[ c_j \cos \left( \frac{(j-1)\pi}{n} \left(i - \frac{1}{2}\right) \right) \right]_{i,j=1}^n, \quad c_j = \begin{cases} \frac{1}{\sqrt{2}}, & j = 1, \\ 1, & j \ge 2, \end{cases} $$

$$ \Lambda(1, 1) = \operatorname{diag}_{j=1,...,n} \left[ 2 \cos \left( \frac{(j-1)\pi}{n} \right) \right]. $$
* Consider $\epsilon = 0$ and $\varphi = 1$. Then,

$$Q_n(0, 1) = \sqrt{\frac{4}{2n+1}} \left[ \sin \left( \frac{i(2j-1)\pi}{2n+1} \right) \right]_{i,j=1}^n, \\ \Lambda(0, 1) = \operatorname{diag}_{j=1,...,n} \left[ 2 \cos \left( \frac{(2j-1)\pi}{2n+1} \right) \right].$$

* Consider $\epsilon = \varphi = -1$. Then,

$$Q_n(-1, -1) = \sqrt{\frac{2}{n}} \left[ c_j \sin \left( \frac{j\pi}{n} \left(i - \frac{1}{2}\right) \right) \right]_{i,j=1}^n, \quad c_j = \begin{cases} \frac{1}{\sqrt{2}}, & j = n, \\ 1, & j \le n-1, \end{cases} \\ \Lambda(-1, -1) = \operatorname{diag}_{j=1,...,n} \left[ 2 \cos \left( \frac{j\pi}{n} \right) \right].$$

For additional properties of elements within the matrix algebras $\tau_n(\epsilon, \varphi)$ in terms of sine and cosine transforms, as well as the computational cost of the matrix-vector product $Q_n(\epsilon, \varphi)v$ for a given vector $v$ in a preconditioning setting, we refer the reader to [13] along with the references cited therein.

**Acknowledgements** C. Manni, S. Serra-Capizzano, and H. Speleers are members of the research group GNCS (Gruppo Nazionale per il Calcolo Scientifico) of INdAM (Istituto Nazionale di Alta Matematica). C. Manni and H. Speleers acknowledge the MUR Excellence Department Project MatMod@TOV (CUP E83C23000330006) awarded to the Department of Mathematics of the University of Rome Tor Vergata. They are also partially supported by a Project of Relevant National Interest (PRIN) under the National Recovery and Resilience Plan (PNRR) funded by the European Union – Next Generation EU (CUP E53D23017910001), by the Italian Research Center in High Performance Computing, Big Data and Quantum Computing (CUP E83C22003230001), and by a GNCS project (CUP E53C23001670001). The work of S. Serra-Capizzano is partially supported by a PRIN-PNRR project (CUP J53D23003780006) and by a GNCS project (CUP E53C23001670001). He is also funded via the European High-Performance Computing Joint Undertaking (JU) under grant agreement No 955701. The JU receives support from the European Union's Horizon 2020 research and innovation programme and Belgium, France, Germany, Switzerland. Furthermore, he is grateful for the support of the Laboratory of Theory, Economics and Systems - Department of Computer Science at Athens University of Economics and Business and to the "Como Lake Center for AstroPhysics" of the University of Insubria.

**Funding** Open access funding provided by Università degli Studi di Roma Tor Vergata within the CRUI-CARE Agreement.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
References

1. Bapat, R.B., Sunder, V.S.: On hypergroups of matrices. Linear and Multilinear Algebra **29**, 125–140 (1991)

2. Barbarino, G., Garoni, C., Serra-Capizzano, S.: Block generalized locally Toeplitz sequences: Theory and applications in the multidimensional case. Electronic Transactions on Numerical Analysis **53**, 113–216 (2020)

3. Bini, D., Capovani, M.: Spectral and computational properties of band symmetric Toeplitz matrices. Linear Algebra and its Applications **52**, 99–126 (1983)

4. Bolten, M., Ekström, S.-E., Furci, I., Serra-Capizzano, S.: A note on the spectral analysis of matrix sequences via GLT momentary symbols: From all-at-once solution of parabolic problems to distributed fractional order matrices. Electronic Transactions on Numerical Analysis **58**, 136–163 (2023)

5. Bozzo, E., Di Fiore, C.: On the use of certain matrix algebras associated with discrete trigonometric transforms in matrix displacement decomposition. SIAM Journal on Matrix Analysis and Applications **16**, 312–326 (1995)

6. Bressan, A., Sande, E.: Approximation in FEM, DG and IGA: A theoretical comparison. Numerische Mathematik **143**, 923–942 (2019)

7. Butler, P., Cantoni, A.: Eigenvalues and eigenvectors of symmetric centrosymmetric matrices. Linear Algebra and its Applications **13**, 275–288 (1976)

8. Cottrell, J.A., Hughes, T.J.R., Bazilevs, Y.: Isogeometric Analysis: Toward Integration of CAD and FEA. John Wiley & Sons, Chichester (2009)

9. Cottrell, J.A., Reali, A., Bazilevs, Y., Hughes, T.J.R.: Isogeometric analysis of structural vibrations. Computer Methods in Applied Mechanics and Engineering **195**, 5257–5296 (2006)

10. de Boor, C.: A Practical Guide to Splines, vol. 27. Springer-Verlag, New York (1978)

11. Deng, Q.: Analytical solutions to some generalized and polynomial eigenvalue problems. Special Matrices **9**, 240–256 (2021)

12. Deng, Q., Calo, V.M.: A boundary penalization technique to remove outliers from isogeometric analysis on tensor-product meshes. Computer Methods in Applied Mechanics and Engineering **383**, 113907 (2021)

13. Di Benedetto, F., Serra-Capizzano, S.: Optimal multilevel matrix algebra operators. Linear and Multilinear Algebra **48**, 35–66 (2000)

14. Di Fiore, C., Zellini, P.: Matrix decompositions using displacement rank and classes of commutative matrix algebras. Linear Algebra and its Applications **229**, 49–99 (1995)

15. Di Vona, E.: Kolmogorov $n$-width e spazi spline ottimi. Università degli Studi di Roma Tor Vergata, Tesi di Laura Magistrale (2019)

16. Donatelli, M., Garoni, C., Manni, C., Serra-Capizzano, S., Speleers, H.: Spectral analysis and spectral symbol of matrices in isogeometric collocation methods. Mathematics of Computation **85**, 1639–1680 (2016)

17. Donatelli, M., Garoni, C., Manni, C., Serra-Capizzano, S., Speleers, H.: Symbol-based multigrid methods for Galerkin B-spline isogeometric analysis. SIAM Journal on Numerical Analysis **55**, 31–62 (2017)

18. Ekström, S.-E., Furci, I., Garoni, C., Manni, C., Serra-Capizzano, S., Speleers, H.: Are the eigenvalues of the B-spline isogeometric analysis approximation of $-\Delta u = \lambda u$ known in almost closed form? Numerical Linear Algebra with Applications **25**, e2198 (2018)

19. Floater, M.S., Sande, E.: Optimal spline spaces of higher degree for $L^2$ $n$-widths. Journal of Approximation Theory **216**, 1–15 (2017)

20. Floater, M.S., Sande, E.: Optimal spline spaces for $L^2$ $n$-width problems with boundary conditions. Constructive Approximation **50**, 1–18 (2019)

21. Garoni, C., Manni, C., Pelosi, F., Serra-Capizzano, S., Speleers, H.: On the spectrum of stiffness matrices arising from isogeometric analysis. Numerische Mathematik **127**, 751–799 (2014)

22. Garoni, C., Manni, C., Pelosi, F., Speleers, H.: Spectral analysis of matrices resulting from isogeometric immersed methods and trimmed geometries. Computer Methods in Applied Mechanics and Engineering **400**, 115551 (2022)

23. Garoni, C., Manni, C., Serra-Capizzano, S., Speleers, H.: NURBS in isogeometric discretization methods: A spectral analysis. Numerical Linear Algebra with Applications **27**, e2318 (2020)

24. Garoni, C., Serra-Capizzano, S.: Generalized Locally Toeplitz Sequences: Theory and Applications, vol. 1. Springer, Cham (2017)
25. Garoni, C., Serra-Capizzano, S.: Generalized Locally Toeplitz Sequences: Theory and Applications, vol. 2. Springer, Cham (2018)

26. Garoni, C., Serra-Capizzano, S., Sesana, D.: Spectral analysis and spectral symbol of $d$-variate $Q_p$ Lagrangian FEM stiffness matrices. SIAM Journal on Matrix Analysis and Applications 36, 1100–1128 (2015)

27. Garoni, C., Speleers, H., Ekström, S.-E., Reali, A., Serra-Capizzano, S., Hughes, T.J.R.: Symbol-based analysis of finite element and isogeometric B-spline discretizations of eigenvalue problems: Exposition and review. Archives of Computational Methods in Engineering 26, 1639–1690 (2019)

28. Hiemstra, R.R., Hughes, T.J.R., Reali, A., Schillinger, D.: Removal of spurious outlier frequencies and modes from isogeometric discretizations of second- and fourth-order problems in one, two, and three dimensions. Computer Methods in Applied Mechanics and Engineering 387, 114115 (2021)

29. Hughes, T.J.R., Cottrell, J.A., Bazilevs, Y.: Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. Computer Methods in Applied Mechanics and Engineering 194, 4135–4195 (2005)

30. Hughes, T.J.R., Evans, J.A., Reali, A.: Finite element and NURBS approximations of eigenvalue, boundary-value, and initial-value problems. Computer Methods in Applied Mechanics and Engineering 272, 290–320 (2014)

31. Lyche, T., Manni, C., Speleers, H.: Foundations of spline theory: B-splines, spline approximation, and hierarchical refinement. In: Lyche, T., Manni, C., Speleers, H. (eds.) Splines and PDEs: From Approximation Theory to Numerical Linear Algebra. Lecture Notes in Mathematics, vol. 2219, pp. 1–76. Springer, Cham (2018)

32. Manni, C., Sande, E., Speleers, H.: Application of optimal spline subspaces for the removal of spurious outliers in isogeometric discretizations. Computer Methods in Applied Mechanics and Engineering 389, 114260 (2022)

33. Manni, C., Sande, E., Speleers, H.: Outlier-free spline spaces for isogeometric discretizations of biharmonic and polyharmonic eigenvalue problems. Computer Methods in Applied Mechanics and Engineering 417, 116314 (2023)

34. Sande, E., Manni, C., Speleers, H.: Sharp error estimates for spline approximation: Explicit constants, $n$-widths, and eigenfunction convergence. Mathematical Models and Methods in Applied Sciences 29, 1175–1205 (2019)

35. Sande, E., Manni, C., Speleers, H.: Explicit error estimates for spline approximation of arbitrary smoothness in isogeometric analysis. Numerische Mathematik 144, 889–929 (2020)

36. Schumaker, L.L.: Spline Functions: Basic Theory. Cambridge University Press, Cambridge (2007)

37. Serra-Capizzano, S.: Some theorems on linear positive operators and functionals and their applications. Computers & Mathematics with Applications 39, 139–167 (2000)

38. Sogn, J., Takacs, S.: Robust multigrid solvers for the biharmonic problem in isogeometric analysis. Computers & Mathematics with Applications 77, 105–124 (2019)

39. Takacs, S., Takacs, T.: Approximation error estimates and inverse inequalities for B-splines of maximum smoothness. Mathematical Models and Methods in Applied Sciences 26, 1411–1445 (2016)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Authors and Affiliations

Noureddine Lamsahel¹,² · Carla Manni³ · Ahmed Ratnani¹ ·
Stefano Serra-Capizzano⁴,⁵ · Hendrik Speleers³

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

Hendrik Speleers

speleers@mat.uniroma2.it

Noureddine Lamsahel

noureddine.lamsahel@um6p.ma, noureddine.lamsahel@etu.univ-littoral.fr
Carla Manni
manni@mat.uniroma2.it

Ahmed Ratnani
ahmed.ratnani@um6p.ma

Stefano Serra-Capizzano
s.serracapizzano@uninsubria.it, stefano.serra@it.uu.se

1. The UM6P Vanguard Center, Mohammed VI Polytechnic University, Hay Moulay Rachid, Lot 660, 43150 Benguerir, Morocco

2. Laboratory of Pure and Applied Mathematics, University of Littoral Côte d'Opale, Rue F. Buisson 50, 62228 Calais-Cedex, France

3. Department of Mathematics, University of Rome Tor Vergata, Via della Ricerca Scientifica 1, 00133 Rome, Italy

4. Department of Science and High Technology, University of Insubria, Via Valleggio 11, 22100 Como, Italy

5. Department of Information Technology, Uppsala University, Lägerhyddsv. 2, Hus 2, P.O. Box 337, SE-751 05 Uppsala, Sweden