1 Introduction↩︎

Isogeometric Galerkin methods based on spline spaces with maximal smoothness on uniform grids provide an excellent approximation of almost the entire spectrum of the Laplace operator, in contrast to what is observed when \(C^0\) finite element methods are used [1]. These spline discretizations, however, suffer from a small number of spurious eigenvalues, called outliers. The number of outliers increases with the polynomial degree but, in the univariate case, is independent of the number of degrees of freedom.

Although limited in number, outliers degrade the overall quality of the discretization, affect the choice of time steps in explicit dynamics, and jeopardize the accurate approximation of high frequencies. It is well known that the presence of outliers is related to the treatment of boundary conditions, which induces low-rank perturbations in the otherwise clean algebraic structure of the matrices involved in the discretization [2], [3].

Recently, spline subspaces have been developed that are theoretically proven to be outlier-free, while still preserving the good approximation accuracy of the entire spectrum. They are obtained by imposition of special boundary conditions and are optimal in the sense of Kolmogorov \(n\)-widths [3].

The efficient use of such outlier-free spaces relies on the ability to represent their elements in terms of a basis with good computational properties. Since the spaces of interest are subspaces of classical spline spaces of maximal smoothness whose elements satisfy additional homogeneous boundary conditions, it is natural to seek bases that coincide with B-splines, except for a few basis elements affected by the new boundary conditions. For any degree \(p\) and for the most common boundary conditions (Dirichlet, Neumann, and mixed), a basis of the corresponding optimal space with this structure has been proposed in [3], [4], which can be easily constructed in terms of cardinal B-splines; see also [5]–[7] for some related work.

It turns out that the bases considered in [3], [4] additionally possess remarkable spectral properties: the eigenvectors of the corresponding stiffness matrix all coincide with the eigenvectors of the corresponding mass matrix, i.e., these matrices are simultaneously diagonalizable. Moreover, the corresponding mass and stiffness matrices admit explicitly known closed-form expressions for their eigenvalues. Therefore, such bases are called outlier-free. The above properties follow from the fact that the considered matrices exhibit a Toeplitz-minus-Hankel or Toeplitz-plus-Hankel structure [6]; see also [8]. This spectral knowledge is relevant not only from a theoretical point of view but also because it enables the construction of efficient solvers for the associated linear systems.

In this paper, we characterize all bases of optimal spline spaces – for any degree \(p\) and any type of boundary condition (Dirichlet, Neumann, and mixed) – that satisfy the following properties (see Definition 1):

• the basis elements have the same support structure as those of the bases in [3], [4];

• the corresponding mass and stiffness matrices are simultaneously diagonalizable.

We prove that, for each optimal space, any basis enjoying these properties can be obtained from the corresponding basis in [3], [4] by an orthogonal matrix; see Theorem 5. This implies that the corresponding mass and stiffness matrices share the same eigenvalues; see Corollary 2. In other words, all the selected bases are outlier-free.

The remainder of the paper is divided into five sections. Section 2 collects notation and preliminary material on splines and graphs, while Section 3 briefly summarizes the state of the art on outlier-free optimal spline spaces and their bases. Section 4 contains the main result of the paper (Theorem 5) and its proof. Section 5 presents a numerical optimization procedure that can be used to construct outlier-free bases with no pre-knowledge about any of them. We end in Section 6 with some final remarks.

2 Preliminaries and Notation↩︎

In this section, we introduce some notation and basic results related to splines and graphs, of interest in the rest of the paper.

2.1 Splines of Maximal Smoothness↩︎

We start by summarizing some results about spline spaces and their bases. For more details, we refer to [9] and references therein.

Let \(1\leq p \in \mathbb{N}\) and let \(\boldsymbol{\tau} \mathrel{\vcenter{:}}= \{0 = \tau_1 < \cdots < \tau_{m+1} = 1\}\) be a partition of the interval \([0,1]\) in \(m\) elements. The spline space of degree \(p\) and maximal smoothness on the partition \(\boldsymbol{\tau}\) is defined by \[\mathbb{S}_{p,\boldsymbol{\tau}} \mathrel{\vcenter{:}}= \left\{ s \in C^{p-1}([0,1]) : s_{|[\tau_i,\tau_{i+1}]} \in \mathbb{P}_p,\;i=1,\ldots,m \right\},\] where \(\mathbb{P}_p\) is the space of polynomials of degree at most \(p\). This space has dimension \(m+p\).

Let \(N \in \mathbb{N}\) such that \(N > p \geq 1\), and let \(\boldsymbol{\xi}\) be a non-decreasing sequence of real values, called knots, \[\boldsymbol{\xi} \mathrel{\vcenter{:}}= \{\xi_1 \le \xi_2 \le \cdots \le \xi_{N+p+1}\}.\] For \(j=1,\ldots,N\), the \(j\)-th B-spline of degree \(p\) associated with the knots \(\boldsymbol{\xi}\) is defined recursively by \[\mathcal{B}_{j,p,\boldsymbol{\xi}}(x) \mathrel{\vcenter{:}}= \frac{x - \xi_j}{\xi_{j+p} - \xi_j} \, \mathcal{B}_{j,p-1,\boldsymbol{\xi}}(x) \;+\; \frac{\xi_{j+p+1} - x}{\xi_{j+p+1} - \xi_{j+1}} \, \mathcal{B}_{j+1,p-1,\boldsymbol{\xi}}(x),\] with \[\mathcal{B}_{j,0,\boldsymbol{\xi}} \mathrel{\vcenter{:}}= \begin{cases} 1, & \text{if } x \in [\xi_j, \xi_{j+1}), \\ 0, & \text{otherwise}. \end{cases}\] Here fractions with zero denominator are set equal to zero and we adopt the convention \[\mathcal{B}_{j,p,\boldsymbol{\xi}}(\xi_{N+p+1})\mathrel{\vcenter{:}}=\lim_{x \to \xi_{N+p+1}^{-}}\mathcal{B}_{j,p,\boldsymbol{\xi}}(x).\] The \(N=m+p\) B-splines of degree \(p\) associated with the knots \[\label{osaitu} \{\underbrace{\tau_1 = \cdots = \tau_1}_{p+1} < \tau_2 < \cdots <\tau_m< \underbrace{\tau_{m+1} = \cdots = \tau_{m+1}}_{p+1}\}\tag{1}\] form a basis of \(\mathbb{S}_{p,\boldsymbol{\tau}}\), which is denoted by \(\mathscr{B}_{p,\boldsymbol{\tau}}\) and referred to as the open-knot B-spline basis.

Finally, the cardinal B-spline of degree \(p\) is defined recursively by \[\mathcal{N}_p(x) \mathrel{\vcenter{:}}= \frac{x}{p}\,\mathcal{N}_{p-1}(x) + \frac{p+1-x}{p}\,\mathcal{N}_{p-1}(x-1),\] with \[\mathcal{N}_0(x) \mathrel{\vcenter{:}}= \begin{cases} 1, & \text{if } x \in [0,1),\\ 0, & \text{otherwise}. \end{cases}\] For \(t \in \mathbb{Q}\), we set \(\mathcal{N}_{t,p} \mathrel{\vcenter{:}}= \mathcal{N}_p(x - t)\), i.e., a rational translate of the cardinal B-spline of degree \(p\).

2.2 Graph Theory↩︎

We now collect some basic facts on graphs, to be used later in the paper. For more details, we refer to [10] and references therein.

An undirected graph is a pair \(G \mathrel{\vcenter{:}}= (V,E)\), where \(V \mathrel{\vcenter{:}}= \{1,\ldots,n\}\) is the set of nodes (or vertices) and \(E \subseteq \big\{\{u,v\}\subseteq V: u\neq v\big\}\) is the set of edges (unordered pairs of nodes). We assume graphs to be simple (no loops and no multiple edges).

With a given graph, one can associate a matrix encoding the node–edge relations. Let \(G \mathrel{\vcenter{:}}= (V,E)\) be an undirected graph with \(|V|=n\) and \(|E|=m\). After fixing an arbitrary orientation for each edge \(e\in E\) (which is only used to define consistent signs), the incidence matrix is the matrix \(B\in\mathbb{R}^{n\times m}\), whose columns are indexed by edges and rows by nodes, defined by \[B_{v,e} \mathrel{\vcenter{:}}= \begin{cases} -1, & \text{if v is the \emph{head} of e},\\ +1, & \text{if v is the \emph{tail} of e},\\ 0, & \text{otherwise}. \end{cases}\]

An undirected graph \(G\) is said to be connected if for every pair of nodes \(u,v\in V\) there exists a path joining them. In general, \(G\) can be decomposed into \(c(G)\) maximal connected components, which are pairwise disjoint on the node set.

The following classical theorem characterizes the rank of the incidence matrix in terms of the connected components of the graph.

3 Laplace Eigenvalue Problem: Isogeometric Discretizations↩︎

We consider the eigenvalue problem associated with the 1D Laplace operator, \[\label{Laplace} -u''= \omega^2 u, \quad \text{in } (0,1),\tag{2}\] and the following standard boundary conditions:

• Dirichlet boundary conditions (also referred to as fixed or type 0 boundary conditions), \[\label{Dirichlet} u(0)=u(1)=0;\tag{3}\]

• Neumann boundary conditions (also referred to as free or natural or type 1 boundary conditions), \[\label{Neumann} u'(0)=u'(1)=0;\tag{4}\]

• a combination of the previous ones (also referred to as mixed or type 2 boundary conditions), \[\label{Miste} u(0)=u'(1)=0.\tag{5}\]

The non-trivial exact solutions of 2 subject to one of the boundary conditions 3 –5 form a numerable set of trigonometric functions, respectively, \[\begin{align} {3} u_l(x) &\mathrel{\vcenter{:}}= \sin(\omega_l x),\quad &\omega_l &\mathrel{\vcenter{:}}= l \pi, \quad &l &= 1,2,\ldots \\ u_l(x) &\mathrel{\vcenter{:}}= \cos(\omega_l x),\quad &\omega_l &\mathrel{\vcenter{:}}= l \pi, \quad &l &= 0,1,2,\ldots \\ u_l(x) &\mathrel{\vcenter{:}}= \sin(\omega_l x),\quad &\omega_l &\mathrel{\vcenter{:}}= (l-1/2)\pi, \quad &l &=1,2,\ldots \end{align}\]

The weak form of problem 2 reads as follows. Find non-trivial \(u \in \mathbb{V}\) and \(\omega^2 \in \mathbb{R}\) such that \[\int_0^1 u'(x)\, v'(x)\, \textrm{d}x = \omega^2 \int_0^1 u(x)\, v(x)\, \textrm{d}x \quad \forall v \in \mathbb{V},\] where \(\mathbb{V}\) is selected in accordance with the chosen boundary conditions:

• Dirichlet boundary conditions: \(\mathbb{V}=\{v \in \mathbb{H}^1(0,1) : v(0) = v(1) = 0 \}\);

• Neumann boundary conditions: \(\mathbb{V}=\mathbb{H}^1(0,1)\);

• mixed boundary conditions: \(\mathbb{V}=\{v \in \mathbb{H}^1(0,1) : v(0) = 0 \}\).

Following the Galerkin approach, we consider a finite-dimensional subspace \(\mathbb{V}_h\) of \(\mathbb{V}\) spanned by a basis \(\{\varphi_1, \ldots, \varphi_{n_h} \}\) and we look for approximate values \(\omega_h\) to \(\omega\) by solving \[K_h \mathbf{u}_h = (\omega_h)^2M_h\mathbf{u}_h,\] where the stiffness matrix \(K_h\) and the mass matrix \(M_h\) consist of the elements \[\label{stiff-mass} K_{h,i,j} \mathrel{\vcenter{:}}= \int_0^1 \varphi_j'(x)\varphi_i'(x)\, \textrm{d}x, \quad M_{h,i,j} \mathrel{\vcenter{:}}= \int_0^1\varphi_j(x)\varphi_i(x)\, \textrm{d}x,\tag{6}\] with \(i,j = 1,\ldots,n_h\). For \(l = 1, \ldots,n_h\), an approximation of the frequency \(\omega_l\) is given by the square root of the \(l\)-th eigenvalue of \(M_h^{-1}K_h\), denoted by \(\omega_{h,l}\). Here we assume that those eigenvalues are given in ascending order. Similarly, an approximation of the eigenfunction \(u_l\) is obtained by considering \[u_{h,l}(x) \mathrel{\vcenter{:}}= \sum_{i=1}^{n_h} u_{h,l,i}\varphi_i(x),\] where \(\mathbf{u}_{h,l} \mathrel{\vcenter{:}}= (u_{h,l,1}, \ldots, u_{h,l,n_h})\) is the \(l\)-th eigenvector of \(M_h^{-1}K_h\), properly normalized. More information on this eigenvalue problem can be found in [11].

In the classical isogeometric approach, the finite-dimensional subspace \(\mathbb{V}_h\) is selected as a proper subspace of the spline space \(\mathbb{S}_{p,\boldsymbol{\tau}}\) in accordance with the chosen boundary conditions:

• Dirichlet boundary conditions: \(\mathbb{V}_h=\{s \in \mathbb{S}_{p,\boldsymbol{\tau}} : s(0) = s(1) = 0 \}\);

• Neumann boundary conditions: \(\mathbb{V}_h= \mathbb{S}_{p,\boldsymbol{\tau}}\);

• mixed boundary conditions: \(\mathbb{V}_h=\{s \in \mathbb{S}_{p,\boldsymbol{\tau}} : s(0) = 0 \}\).

The above choices, with a uniform partition \(\boldsymbol{\tau}\), produce a very good approximation of the continuous spectrum, with increasing accuracy as the degree \(p\) increases [1]. However, the scheme still suffers from some numerical artifacts: a small portion of frequencies is poorly approximated. These spurious approximations are called outliers; see [3] and references therein.

3.1 Outlier-Free Spline Spaces↩︎

The outliers in the isogeometric eigenvalue approximation of the Laplace operator can be removed by taking special spline subspaces as the discretization space \(\mathbb{V}_h\). These subspaces turn out to be optimal in the sense of Kolmogorov \(n\)-widths for appropriate function classes identified by problem 2 with boundary conditions 3 –5 ; see [3], [5]. More precisely, let us consider the following \(n\)-dimensional spline spaces:

• Dirichlet boundary conditions: \[\label{S95p950} \mathbb{S}_{p,0} \mathrel{\vcenter{:}}= \bigl\{ s \in \mathbb{S}_{p,\boldsymbol{\tau}_{p,0}} : \partial^\alpha s(0) = \partial^\alpha s(1) = 0, \; 0 \le \alpha \le p, \;\alpha \;\text{even} \bigr\},\tag{7}\] where \[\boldsymbol{\tau}_{p,0} \mathrel{\vcenter{:}}= \begin{cases} \,\bigl\{0, \frac{1}{n+1}, \frac{2}{n+1}, \ldots, \frac{n}{n+1}, 1\bigr\}, & p \text{ odd},\\[3pt] \,\bigl\{0, \frac{1/2}{n+1}, \frac{3/2}{n+1}, \ldots, \frac{n+1/2}{n+1}, 1\bigr\}, & p \text{ even}; \end{cases}\]

• Neumann boundary conditions: \[\label{S95p951} \mathbb{S}_{p,1} \mathrel{\vcenter{:}}= \bigl\{ s \in \mathbb{S}_{p,\boldsymbol{\tau}_{p,1}} : \partial^\alpha s(0) = \partial^\alpha s(1) = 0, \; 0 \le \alpha \le p, \;\alpha \;\text{odd} \bigr\},\tag{8}\] where \[\boldsymbol{\tau}_{p,1} \mathrel{\vcenter{:}}= \begin{cases} \,\bigl\{0, \frac{1/2}{n}, \frac{3/2}{n}, \ldots, \frac{n-1/2}{n}, 1\bigr\}, & p \text{ odd},\\[3pt] \,\bigl\{0, \frac{1}{n}, \frac{2}{n}, \ldots, \frac{n-1}{n}, 1\bigr\}, & p \text{ even}; \end{cases}\]

• mixed boundary conditions: \[\label{S95p952} \begin{align} \mathbb{S}_{p,2} \mathrel{\vcenter{:}}= \bigl\{ s \in \mathbb{S}_{p,\boldsymbol{\tau}_{p,2}} : \partial^{\alpha_0} s(0) = \partial^{\alpha_1} s(1) = 0,\;& \;0 \le \alpha_0,\alpha_1 \le p,\;\\ & \;\alpha_0 \;\text{even}, \alpha_1 \;\text{odd} \bigr\}, \end{align}\tag{9}\] where \[\boldsymbol{\tau}_{p,2} \mathrel{\vcenter{:}}= \begin{cases} \,\bigl\{0, \frac{2}{2n+1}, \frac{4}{2n+1}, \ldots, \frac{2n}{2n+1}, 1\bigr\}, & p \text{ odd},\\[3pt] \,\bigl\{0, \frac{1}{2n+1}, \frac{3}{2n+1}, \ldots, \frac{2n-1}{2n+1}, 1\bigr\}, & p \text{ even}. \end{cases}\]

When taking \(\mathbb{V}_h=\mathbb{S}_{p,i}\), \(i=0,1,2\) in the Galerkin approach, accurate approximations are obtained for all the first \(n\) eigenvalues and eigenfunctions of the Laplace operator, with boundary conditions of type \(0,1,2\), respectively [3]. For this reason, such spaces are called outlier-free.

It is worth mentioning that the eigenvalues of \(M^{-1}_h K_h\) do not depend on the choice of the basis in the considered discretization space. On the contrary, those of the mass matrix \(M_h\) and of the stiffness matrix \(K_h\) clearly depend on the selected basis. In the next section, we present special bases for the spaces \(\mathbb{S}_{p,i}\), \(i=0,1,2\), such that the corresponding mass and stiffness matrices enjoy particularly nice spectral properties.

3.2 Outlier-Free Spline Bases↩︎

The effective use of the spaces \(\mathbb{S}_{p,i}\) is clearly tied to the possibility of representing their elements in terms of basis functions with pleasing properties, such as local support and non-negativity. In this regard, the following bases were proposed in [3], [4] (see also [5], [6]):

• for \(n \geq p + 1\), the basis \(\mathscr{E}_{p,0}\) for the space \(\mathbb{S}_{p,0}\) is given by the \(n\) functions \[\label{Base95Eliseo} \mathcal{E}_{j,p,0}(x) \mathrel{\vcenter{:}}= \begin{cases} -\mathcal{N}_{q-j,p}(y) + \mathcal{N}_{q+j,p}(y), & j \in \left\{ 1, \ldots, \left\lfloor \frac{p}{2} \right\rfloor \right\}, \\[3pt] \mathcal{N}_{q+j,p}(y), & j \in \left\{ \left\lfloor \frac{p}{2} \right\rfloor + 1, \ldots, n - \left\lfloor \frac{p}{2} \right\rfloor \right\}, \\[3pt] \mathcal{N}_{q+j,p}(y) - \mathcal{N}_{q + 2(n+1) - j,p}(y), & j \in \left\{ n + 1 - \left\lfloor \frac{p}{2} \right\rfloor, \ldots, n \right\}, \end{cases}\tag{10}\] where \(q \mathrel{\vcenter{:}}= \left\lfloor \frac{p}{2} \right\rfloor - p\), \(y(x) \mathrel{\vcenter{:}}= (n+1)x\) if \(p\) odd, \(y(x) \mathrel{\vcenter{:}}= (n+1)x + \frac{1}{2}\) if \(p\) even;

• for \(n \geq 2p - 2\lfloor \frac{p}{2} \rfloor + 1\), the basis \(\mathscr{E}_{p,1}\) for the space \(\mathbb{S}_{p,1}\) is given by the \(n\) functions \[\label{Base95Eliseo952} \mathcal{E}_{j,p,1}(x) \mathrel{\vcenter{:}}= \begin{cases} \mathcal{N}_{-\left\lfloor \frac{p}{2} \right\rfloor - j,p}(y) + \mathcal{N}_{-\left\lfloor \frac{p}{2} \right\rfloor + j - 1,p}(y), & j \in \left\{ 1, \ldots, -q \right\}, \\[3pt] \mathcal{N}_{-\left\lfloor \frac{p}{2} \right\rfloor + j - 1,p}(y), & j \in \left\{ -q + 1, \ldots, n + q \right\}, \\[3pt] \mathcal{N}_{-\left\lfloor \frac{p}{2} \right\rfloor + j - 1,p}(y) + \mathcal{N}_{-\left\lfloor \frac{p}{2} \right\rfloor + 2n - j,p}(y), & j \in \left\{ n + q + 1, \ldots, n \right\}, \end{cases}\tag{11}\] where \(q \mathrel{\vcenter{:}}= \left\lfloor \frac{p}{2} \right\rfloor - p\), \(y(x) \mathrel{\vcenter{:}}= nx + \frac{1}{2}\) if \(p\) odd, \(y(x) \mathrel{\vcenter{:}}= nx\) if \(p\) even;

• for \(n \geq p + 1\), the basis \(\mathscr{E}_{p,2}\) for the space \(\mathbb{S}_{p,2}\) is given by the \(n\) functions \[\label{Base95Eliseo953} \mathcal{E}_{j,p,2}(x) \mathrel{\vcenter{:}}= \mathcal{N}_{q+j,p}(y) + \gamma_j \mathcal{N}_{q-k_j,p}(y), \qquad j = 1, \ldots, n,\tag{12}\] where \[(\gamma_j, k_j) \mathrel{\vcenter{:}}= \begin{cases} (-1, j), & j \in \left\{ 1, \ldots, \left\lfloor \frac{p}{2} \right\rfloor \right\}, \\[3pt] (0, \cdot), & j \in \left\{ \left\lfloor \frac{p}{2} \right\rfloor + 1, \ldots, n + q \right\}, \\[3pt] (1, j - (2n + 1)), & j \in \left\{ n + q + 1, \ldots, n \right\}, \end{cases}\] and \(q \mathrel{\vcenter{:}}= \left\lfloor \frac{p}{2} \right\rfloor - p\), \(y(x) \mathrel{\vcenter{:}}= \frac{2n + 1}{2}\, x\) if \(p\) odd, \(y(x) \mathrel{\vcenter{:}}= \frac{2n + 1}{2}\, x + \frac{1}{2}\) if \(p\) even.

By construction, the basis functions in 10 –12 inherit the good computational properties of cardinal B-splines (and so of B-splines). Furthermore, they possess some remarkable spectral properties as detailed in the following.

Given a basis \(\mathscr{B}\) of one of the discretization spaces 7 –9 , we denote by \(M_{\mathscr{B}}\) and \(K_{\mathscr{B}}\) the corresponding mass and stiffness matrices. Recall that \(1\leq p \in \mathbb{N}\). For \(r=0,1\) and \(\theta \in [0,\pi]\), let \[\label{simboli} g^r_p(\theta) \mathrel{\vcenter{:}}= (-1)^r \mathcal{N}^{(2r)}_{2p+1}(p+1) + 2(-1)^r \sum_{k=1}^{p} \mathcal{N}^{(2r)}_{2p+1}(p+1-k)\cos(k\theta).\tag{13}\] It is known that \(g^0_p(\theta)>0\) for \(\theta \in[0,\pi]\), and \(g^1_p(\theta)>0\) for \(\theta \in(0,\pi]\); see [2]. Moreover, their ratio \[\label{monotonia} e_p(\theta) \mathrel{\vcenter{:}}= \frac{g^1_p(\theta)}{g^0_p(\theta)}\tag{14}\] is monotone increasing on \([0,\pi]\); see [15]. The following theorems were shown in [6] for the bases \(\mathscr{E}_{p,i}\), \(i = 0,1,2\), introduced above.

Theorems 2 to 4 and the monotonicity of the function 14 directly lead to the following corollary.

The relevance of the above mentioned spectral properties of the bases in 10 –12 is twofold, involving two different meanings of the notion outlier-free.

1. For \(i=0,1,2\), the matrices \(M_{\mathscr{E}_{p,i}}\) and \(K_{\mathscr{E}_{p,i}}\) admit a closed-form description of their spectra, i.e., all their eigenvalues are known samples, up to a scaling, of the functions \(g_p^0\) and \(g_p^1\) in 13 , which are the symbols of \(M_{\mathscr{E}_{p,i}}\) and \(K_{\mathscr{E}_{p,i}}\), respectively [2]. For this reason, the basis \(\mathscr{E}_{p,i}\) is called outlier-free. Note that the functions \(g_p^0\) and \(g_p^1\) are independent of the index \(i\), i.e., of the type of boundary conditions.

2. For \(i=0,1,2\), the matrices \(M_{\mathscr{E}_{p,i}}\) and \(K_{\mathscr{E}_{p,i}}\) are simultaneously diagonalizable, i.e., the eigenvectors of the stiffness matrix all coincide with the eigenvectors of the mass matrix – and so with those of \(M^{-1}_{\mathscr{E}_{p,i}} K_{\mathscr{E}_{p,i}}\). This property ensures that the generalized eigenvalues of \(K_{\mathscr{E}_{p,i}}\mathbf{u} = \lambda M_{\mathscr{E}_{p,i}}\mathbf{u}\) agree with the ratios \({\lambda_j \bigl( K_{\mathscr{E}_{p,i}} \bigr)}/ {\lambda_j \bigl( M_{\mathscr{E}_{p,i}} \bigr)}\). By taking into account the known explicit expressions of \({\lambda_j \bigl( M_{\mathscr{E}_{p,i}} \bigr)}\) and \({\lambda_j \bigl( K_{\mathscr{E}_{p,i}} \bigr)}\), it implies that the corresponding space \(\mathbb{S}_{p,i}\) is outlier-free as discussed in [6].

This raises the question whether there are other bases of the space \(\mathbb{S}_{p,i}\), \(i=0,1,2\), enjoying the same (or very similar) spectral properties and still maintaining a B-spline-like structure (with properties such as local support and non-negativity). This will be investigated in the next section.

4 Outlier-Free and B-Spline-Like Bases for Outlier-Free Spaces↩︎

This section contains the main result of the paper, which is stated in Theorem 5. We show that any basis \(\mathscr{B}_{p,i}\) of the space \(\mathbb{S}_{p,i}\), \(i=0,1,2\), such that

1. the elements of \(\mathscr{B}_{p,i}\) have the same support structure as those of \(\mathscr{E}_{p,i}\),

2. the corresponding mass and stiffness matrices are simultaneously diagonalizable,

can be obtained from the basis \(\mathscr{E}_{p,i}\) by a change of basis identified by an orthogonal matrix \(R_{\mathscr{B}_{p,i}}\). This implies (see Corollary 2) that the mass and stiffness matrices with respect to the basis \(\mathscr{B}_{p,i}\) share the same spectrum as the corresponding matrices with respect to the basis \(\mathscr{E}_{p,i}\), and thus they are outlier-free as well.

4.1 Main Result↩︎

Before we are able to formulate our main result, we need a precise definition of the family of bases of \(\mathbb{S}_{p,i}\) we are interested in. To this end, we start by recalling a well-known linear algebra result; see, e.g., [19].

We observe that for a basis \(\mathscr{B}\), the corresponding mass matrix \(M_{\mathscr{B}}\) and stiffness matrix \(K_{\mathscr{B}}\) are symmetric and thus diagonalizable by orthogonal transformations. This means that, according to Lemma 1, they are simultaneously diagonalizable if and only if their commutator vanishes, namely \[[M_{\mathscr{B}}, K_{\mathscr{B}}] \mathrel{\vcenter{:}}= M_{\mathscr{B}} K_{\mathscr{B}} - K_{\mathscr{B}} M_{\mathscr{B}} = 0.\] In this perspective, we consider bases \(\mathscr{B}_{p,i}\) of \(\mathbb{S}_{p,i}\) satisfying \([M_{\mathscr{B}_{p,i}}, K_{\mathscr{B}_{p,i}}] = 0\).

Furthermore, we want that the basis elements have the same support structure as those of \(\mathscr{E}_{p,i}\). More precisely, from 10 –12 , we deduce the following technical lemma.

The open-knot B-spline basis \(\mathscr{B}_{p,\boldsymbol{\tau}_{p,i}}\) is the basis with the most compact support of its elements for the space \(\mathbb{S}_{p,\boldsymbol{\tau}_{p,i}}\), which induces a pronounced sparsity structure in the mass and stiffness matrices. Due to the block structure of the matrix \(A_{\mathscr{E}_{p,i}}\), the basis \(\mathscr{E}_{p,i}\) mimics these features and modifies only those elements of the open-knot B-spline basis that are affected by the additional boundary conditions characterizing the space \(\mathbb{S}_{p,i}\); see 7 –9 . It is therefore natural to consider other bases of the space \(\mathbb{S}_{p,i}\) that preserve the same sparsity and support structure, i.e., bases that can be expressed in terms of the open-knot B-spline basis through a matrix having the same block structure as the matrix \(A_{\mathscr{E}_{p,i}}\) in ?? .

We are now ready to define the family of bases of \(\mathbb{S}_{p,i}\) of interest.

Our main result relates every element of this family to the basis \(\mathscr{E}_{p,i}\).

The proof of the theorem is provided in Section 4.2. We immediately deduce the following corollary.

Note that, vice versa, whenever the basis \(\mathscr{B}_{p,i}\) is obtained from \(\mathscr{E}_{p,i}\) by an orthogonal matrix, then 15 implies that \(M_{\mathscr{B}_{p,i}} K_{\mathscr{B}_{p,i}} - K_{\mathscr{B}_{p,i}} M_{\mathscr{B}_{p,i}} = 0\).

4.2 Proof of Theorem 5↩︎

This section is devoted to the proof of Theorem 5, which makes use of Chebyshev polynomials, graph theory, and few other preliminary lemmas. For the sake of brevity, we will only report the proof for the case \(\mathbb{S}_{p,0}\). The proof strategy for the cases \(\mathbb{S}_{p,1}\) and \(\mathbb{S}_{p,2}\) is similar, with minor technical changes. For further details, we refer to [20].

We start with a series of lemmas. In the first lemma, we show a trigonometric identity.

Let \(\operatorname{im}(A)\) and \(\ker(A)\) denote the image and kernel of a matrix \(A\). We now look at a special matrix \(C\), of interest later in the proof of Theorem 5, and state some of its properties.

For a given basis \(\mathscr{B}_{p,0} \in \mathfrak{F}_{p,0}\), the next lemma relates the structure of the matrix \(A_{\mathscr{B}_{p,0}}\) and the structure of the change-of-basis matrix \(R_{\mathscr{B}_{p,0}}\).

We now build a matrix \(H\) linked to the structure of the matrix \(R_{\mathscr{B}_{p,0}}\) in ?? and determine its rank by exploiting graph theory. To this end, we define a set of index pairs, which identifies some zero entries of the first row and of the \(\left\lfloor \frac{p}{2}\right\rfloor\)-th column of \(R_{\mathscr{B}_{p,0}}\).

Finally, our last lemma determines the rank of the matrix \(HC\), where \(H\) and \(C\) have been described above.

We have now at our disposal all the ingredients for the proof of Theorem 5.

5 Numerical Procedure↩︎

In this section, we briefly present a numerical procedure for finding B-spline-like and outlier-free bases for the optimal spaces \(\mathbb{S}_{p,i}\), without assuming any knowledge about the bases \(\mathscr{E}_{p,i}\), \(i = 0,1,2\). This is motivated by the prospect of constructing similar bases – if any – for optimal spaces for the isogeometric discretization of biharmonic and polyharmonic eigenvalue problems [21], where no direct explicit constructions analogous to the bases \(\mathscr{E}_{p,i}\) are available.

We focus on the case \(\mathbb{S}_{p,0}\) with \(p\) odd. The other cases can be addressed in a similar manner.

1. Let \(A_{\mathscr{B}_{p,0}}\) be a candidate representation matrix for the new basis \(\mathscr{B}_{p,0}\) in terms of the open-knot B-spline basis \(\mathscr{B}_{p,\boldsymbol{\tau}_{p,0}}\), with the following block structure: \[\label{A95B} A_{\mathscr{B}_{p,0}} = \begin{bmatrix} A_{11} & 0 & 0 \\ 0 & I_{r\times r} & 0 \\ 0 & 0 & A_{33} \end{bmatrix}, \quad A_{11}, A_{33} \in \mathbb{R}^{\left\lfloor \frac{p}{2} \right\rfloor \times p}.\tag{30}\]

2. To ensure that \(\mathscr{B}_{p,0}\) is a basis of the space \(\mathbb{S}_{p,0}\), the basis elements must satisfy the homogeneous boundary conditions specified in 7 . Compute the even derivatives of the open-knot B-splines at the end points and set up the kernel of these conditions. Then, based on this kernel, determine the matrices \(A_{11}\) and \(A_{33}\) in terms of a set of parameters \(X\), whose number only depends on the degree \(p\).

3. Construct the mass and stiffness matrices \(M_{\mathscr{B}_{p,0}}\) and \(K_{\mathscr{B}_{p,0}}\) in terms of the mass and stiffness matrices of the open-knot B-spline basis \(\mathscr{B}_{p,\boldsymbol{\tau}_{p,0}}\), using the representation matrix \(A_{\mathscr{B}_{p,0}}\). Then, minimize (numerically) the Frobenius norm of their commutator \(\|M_{\mathscr{B}_{p,0}} K_{\mathscr{B}_{p,0}} - K_{\mathscr{B}_{p,0}} M_{\mathscr{B}_{p,0}}\|_{F}\) with respect to the above parameters \(X\), imposing (linear) constraints on them that ensure the entries of \(A_{11}\) and \(A_{33}\) to be non-negative.

4. Build \(A_{\mathscr{B}_{p,0}^*}\) using the computed minimized parameter values \(X^*\) and obtain the new basis \(\mathscr{B}_{p,0}^* = A_{\mathscr{B}_{p,0}^*}\,\mathscr{B}_{p,\boldsymbol{\tau}_{p,0}}\).

The solution space of the minimization problem in step 3 is non-empty and the optimal solution is non-unique. The obtained solution will depend on the chosen optimization method and can be possibly steered through initialization feed or additional constraints. This numerical procedure finds justification in the previously presented theoretical result (Section 4.1). Indeed, from steps 1 to 3, we deduce that the basis \(\mathscr{B}_{p,0}^*\) belongs to the family \(\mathfrak{F}_{p,0}\) specified in Definition 1. From Corollary 2, we conclude that the obtained basis is outlier-free. Furthermore, the non-negativity constraints in step 3 ensure that the resulting basis functions are non-negative, being non-negative combinations of (open-knot) B-splines.

We implemented the above numerical procedure and tested it out for the case \(p=5\) and \(n=7\). We obtained a basis \(\mathscr{B}_{p,0}^*\) identified by the representation matrix \(A_{\mathscr{B}_{p,0}^*}\) as in 30 , where \[A_{11}^* = \begin{bmatrix} \;0 & \;0.1675 & \;0.5024 & \;0.8142 & \;0.0859\;\\ \;0 & \;0.0023 & \;0.0069 & \;0.1305 & \;0.9963\;\\ \end{bmatrix},\] and \(A_{33}^*\) has the same columns but in reversed order. We remark that the matrix \(A_{11}^0\) in ?? , corresponding to the basis \(\mathscr{E}_{p,0}\), is given by \[A_{11}^0 = \begin{bmatrix} \;0 & \;1/6 & \;1/2 & \;4/5 & \;0\;\\ \;0 & \;1/60 & \;1/20 & \;1/5 & \;1\;\\ \end{bmatrix},\] and it holds \[A_{11}^* = R_{11}\, A_{11}^0, \quad R_{11} = \begin{bmatrix} \; 0.9963 & \;0.0859\;\\ \;-0.0859 & \;0.9963\;\\ \end{bmatrix}.\] The matrix \(R_{11}\) is clearly a rotation matrix (with angle \(-0.086\)). This is in agreement with Theorem 5, taking into account the matrix structure ?? . The corresponding bases \(\mathscr{E}_{p,0}\) and \(\mathscr{B}_{p,0}^*\) are depicted in Figure 5. Since \(r=3\) in 30 , the three central functions in both bases are standard (cardinal) B-splines. We close this numerical example by observing that any rotation matrix with an angle \(\phi\) such that \(0 \leq -\phi \leq \arctan(1/10)\) would result in a matrix \(A_{11}^*\) with non-negative entries.

6 Conclusion↩︎

The use of optimal spline spaces in the isogeometric Galerkin method for the eigenvalue problem associated with the univariate Laplace operator – subject to any classical boundary conditions: Dirichlet, Neumann, and mixed – is a powerful approximation tool, which is theoretically proven to be outlier-free [3], [6]. While the absence of outliers is an intrinsic property of the considered discretization spaces \(\mathbb{S}_{p,i}\), \(i=0,1,2\), stemming from their optimality, the choice of the basis is of crucial importance from a practical point of view.

We have characterized the bases of these optimal spaces that enjoy a B-spline-like support structure and whose mass and stiffness matrices are simultaneously diagonalizable. It turns out that such bases are orthogonally equivalent to the outlier-free bases \(\mathscr{E}_{p,i}\) proposed in [3], [4], and therefore share their important spectral properties derived in [6]. As a consequence, all the selected bases are also outlier-free.

The obtained results can be reasonably extended to the spaces proposed in [14] (see also [7], [12]), which are closely related to the optimal spline spaces considered here and still provide outlier-free discretizations for the eigenvalue problem associated with the univariate Laplace operator [6].

Beyond their relevance for the Laplace operator, the obtained results and the related numerical procedure suggest a possible extension to the construction of similar bases for optimal spaces for the isogeometric discretization of biharmonic and polyharmonic eigenvalue problems [21]. Indeed, while a family of optimal spline spaces has been identified for the approximation of the eigenvalues of a polyharmonic operator of any order [21], the construction of suitable bases for such spaces remains an open question. Even in the biharmonic case, a direct attempt to construct a basis for the corresponding optimal spaces that reproduces the elegant expression of the basis \(\mathscr{E}_{p,i}\) in terms of cardinal B-splines fails. From this perspective, the results of this paper suggest to seek such a basis numerically by imposing a B-spline-like structure and minimizing a suitable norm of the commutator of the corresponding mass and stiffness matrices, in order to identify bases for which these matrices are simultaneously diagonalizable. This is a promising direction for future research.

References↩︎