Block Generalized Locally Toeplitz Sequences: From the Theory to the Applications

Abstract

The theory of generalized locally Toeplitz (GLT) sequences is a powerful apparatus for computing the asymptotic spectral distribution of matrices $A_n$ arising from virtually any kind of numerical discretization of differential equations (DEs). Indeed, when the mesh fineness parameter n tends to infinity, these matrices $A_n$ give rise to a sequence ${\left\{A_n\right\}}_n$, which often turns out to be a GLT sequence or one of its “relatives”, i.e., a block GLT sequence or a reduced GLT sequence. In particular, block GLT sequences are encountered in the discretization of systems of DEs as well as in the higher-order finite element or discontinuous Galerkin approximation of scalar DEs. Despite the applicative interest, a solid theory of block GLT sequences has been developed only recently, in 2018. The purpose of the present paper is to illustrate the potential of this theory by presenting a few noteworthy examples of applications in the context of DE discretizations.

1. Introduction

The theory of generalized locally Toeplitz (GLT) sequences stems from Tilli’s work on locally Toeplitz (LT) sequences [1] and from the spectral theory of Toeplitz matrices [2,3,4,5,6,7,8,9,10,11,12]. It was then carried forward in [13,14,15,16], and was recently extended by Barbarino [17]. This theory is a powerful apparatus for computing the asymptotic spectral distribution of matrices arising from the numerical discretization of continuous problems, such as integral equations (IEs) and, especially, differential equations (DEs). The experience reveals that virtually any kind of numerical methods for the discretization of DEs gives rise to structured matrices $A_n$ whose asymptotic spectral distribution, as the mesh fineness parameter n tends to infinity, can be computed through the theory of GLT sequences. We refer the reader to ([13] Section 10.5), ([14] Section 7.3), and [15,16,18] for applications of the theory of GLT sequences in the context of finite difference (FD) discretizations of DEs; to ([13] Section 10.6), ([14] Section 7.4), and [16,18,19] for the finite element (FE) case; to [20] for the finite volume (FV) case; to ([13] Section 10.7), ([14] Sections 7.5–7.7), and [21,22,23,24,25,26] for the case of isogeometric analysis (IgA) discretizations, both in the collocation and Galerkin frameworks; and to [27] for a further recent application to fractional DEs. We also refer the reader to ([13] Section 10.4) and [28,29] for a look at the GLT approach for sequences of matrices arising from IE discretizations.

It is worth emphasizing that the asymptotic spectral distribution of DE discretization matrices, whose computation is the main objective of the theory of GLT sequences, is not only interesting from a theoretical viewpoint, but can also be used for practical purposes. For example, it is known that the convergence properties of mainstream iterative solvers, such as multigrid and preconditioned Krylov methods, strongly depend on the spectral features of the matrices to which they are applied. The spectral distribution can then be exploited to design efficient solvers of this kind and to analyze/predict their performance. In this regard, we recall that noteworthy estimates on the superlinear convergence of the conjugate gradient method obtained by Beckermann and Kuijlaars in [30] are closely related to the asymptotic spectral distribution of the considered matrices. Furthermore, in the context of Galerkin and collocation IgA discretizations of elliptic DEs, the spectral distribution computed through the theory of GLT sequences in a series of recent papers [21,22,23,24,25] was exploited in [31,32,33] to devise and analyze optimal and robust multigrid solvers for IgA linear systems.

In the very recent work [34], starting from the original intuition by the third author ([16] Section 3.3), the theory of block GLT sequences has been developed in a systematic way as an extension of the theory of GLT sequences. Such an extension is of the utmost importance in practical applications. In particular, it provides the necessary tools for computing the spectral distribution of block structured matrices arising from the discretization of systems of DEs ([16] Section 3.3) and from the higher-order finite element or discontinuous Galerkin approximation of scalar/vectorial DEs [35,36,37]. The purpose of this paper is to illustrate the potential of the theory of block GLT sequences [34] and of its multivariate version—which combines the results of [34] with the “multivariate technicalities” from [14]—by presenting a few noteworthy examples of applications. Actually, the present paper can be seen as a necessary completion of the purely theoretical work [34].

The paper is organized as follows. In Section 2, we report a summary of the theory of block GLT sequences. In Section 3, we focus on the FD discretization of a model system of univariate DEs; through the theory of block GLT sequences, we compute the spectral distribution of the related discretization matrices. In Section 4, we focus on the higher-order FE approximation of the univariate diffusion equation; again, we compute the spectral distribution of the associated discretization matrices through the theory of block GLT sequences. In Section 5, we summarize the multivariate version of the theory of block GLT sequences, also known as the theory of multilevel block GLT sequences. In Section 6, we describe the general GLT approach for computing the spectral distribution of matrices arising from the discretization of systems of partial differential equations (PDEs). In Section 7, we focus on the B-spline IgA approximation of a bivariate variational problem for the curl–curl operator, which is of interest in magnetostatics; through the theory of multilevel block GLT sequences, we compute the spectral distribution of the related discretization matrices. Final considerations are collected in Section 8.

2. The Theory of Block GLT Sequences

In this section, we summarize the theory of block GLT sequences, which was originally introduced in ([16] Section 3.3) and has been recently revised and systematically developed in [34].

Sequences of Matrices and Block Matrix-Sequences. Throughout this paper, a sequence of matrices is any sequence of the form ${\left\{A_n\right\}}_n$, where $A_n$ is a square matrix of size $d_n$ and $d_n\to \infty$ as $n\to \infty$. Let $s\ge 1$ be a fixed positive integer independent of n; an s-block matrix-sequence (or simply a matrix-sequence if s can be inferred from the context or we do not need/want to specify it) is a special sequence of matrices ${\left\{A_n\right\}}_n$ in which the size of $A_n$ is $d_n=sn$.

Singular Value and Eigenvalue Distribution of a Sequence of Matrices. Let ${\mu }_k$ be the Lebesgue measure in ${\mathbb{R}}^k$. Throughout this paper, all the terminology from measure theory (such as “measurable set”, “measurable function”, “a.e.”, etc.) is referred to the Lebesgue measure. A matrix-valued function $f:D\subseteq {\mathbb{R}}^k\to {\mathbb{C}}^{r\times r}$ is said to be measurable (resp., continuous, Riemann-integrable, in $L^p\left(D\right)$, etc.) if its components $f_{\alpha \beta }:D\to \mathbb{C},\alpha ,\beta =1,\dots ,r$, are measurable (resp., continuous, Riemann-integrable, in $L^p\left(D\right)$, etc.). We denote by $C_c\left(\mathbb{R}\right)$ (resp., $C_c\left(\mathbb{C}\right)$) the space of continuous complex-valued functions with bounded support defined on $\mathbb{R}$ (resp., $\mathbb{C}$). If $A\in {\mathbb{C}}^{m\times m}$, the singular values and the eigenvalues of A are denoted by ${\sigma }_1\left(A\right),\dots ,{\sigma }_m\left(A\right)$ and ${\lambda }_1\left(A\right),\dots ,{\lambda }_m\left(A\right)$, respectively.

Definition 1.

Let ${\left\{A_n\right\}}_n$ be a sequence of matrices, with $A_n$ of size $d_n$, and let $f:D\subset {\mathbb{R}}^k\to {\mathbb{C}}^{r\times r}$ be a measurable function defined on a set D with $0<{\mu }_k\left(D\right)<\infty$.

• We say that ${\left\{A_n\right\}}_n$ has a (asymptotic) singular value distribution described by f, and we write ${\left\{A_n\right\}}_n{\sim }_{\sigma }f$, if

  $$\underset{n\to \infty }{\lim }\frac{1}{d_n}\sum _{i=1}^{d_n}F\left({\sigma }_i\left(A_n\right)\right)=\frac{1}{{\mu }_k\left(D\right)}{\int }_D\frac{{\sum }_{i=1}^{r}F\left({\sigma }_i\left(f\left(\mathbf{x}\right)\right)\right)}{r}\mathrm{d}\mathbf{x},\forall F\in C_c\left(\mathbb{R}\right).$$

  (1)
  In this case, f is referred to as a singular value symbol of ${\left\{A_n\right\}}_n$.
• We say that ${\left\{A_n\right\}}_n$ has a (asymptotic) spectral (or eigenvalue) distribution described by f, and we write ${\left\{A_n\right\}}_n{\sim }_{\lambda }f$, if

  $$\underset{n\to \infty }{\lim }\frac{1}{d_n}\sum _{i=1}^{d_n}F\left({\lambda }_i\left(A_n\right)\right)=\frac{1}{{\mu }_k\left(D\right)}{\int }_D\frac{{\sum }_{i=1}^{r}F\left({\lambda }_i\left(f\left(\mathbf{x}\right)\right)\right)}{r}\mathrm{d}\mathbf{x},\forall F\in C_c\left(\mathbb{C}\right).$$

  (2)
  In this case, f is referred to as a spectral (or eigenvalue) symbol of ${\left\{A_n\right\}}_n$.

If ${\left\{A_n\right\}}_n$ has both a singular value and an eigenvalue distribution described by f, we write ${\left\{A_n\right\}}_n{\sim }_{\sigma ,\lambda }f$.

We note that Definition 1 is well-posed because the functions $\mathbf{x}\mapsto {\sum }_{i=1}^{r}F\left({\sigma }_i\left(f\left(\mathbf{x}\right)\right)\right)$ and $\mathbf{x}\mapsto {\sum }_{i=1}^{r}F\left({\lambda }_i\left(f\left(\mathbf{x}\right)\right)\right)$ are measurable ([34] Lemma 2.1). Whenever we write a relation such as ${\left\{A_n\right\}}_n{\sim }_{\sigma }f$ or ${\left\{A_n\right\}}_n{\sim }_{\lambda }f$, it is understood that f is as in Definition 1; that is, f is a measurable function defined on a subset D of some ${\mathbb{R}}^k$ with $0<{\mu }_k\left(D\right)<\infty$, and f takes values in ${\mathbb{C}}^{r\times r}$ for some $r\ge 1$.

Remark 1.

The informal meaning behind the spectral distribution (2) is the following: assuming that f possesses r Riemann-integrable eigenvalue functions ${\lambda }_i\left(f\left(\mathbf{x}\right)\right)$, $i=1,\dots ,r$, the eigenvalues of $A_n$, except possibly for $o\left(d_n\right)$ outliers, can be subdivided into r different subsets of approximately the same cardinality; and, for n large enough, the eigenvalues belonging to the ith subset are approximately equal to the samples of the ith eigenvalue function ${\lambda }_i\left(f\left(\mathbf{x}\right)\right)$ over a uniform grid in the domain D. For instance, if $k=1$, $d_n=nr$, and $D=[a,b]$, then, assuming we have no outliers, the eigenvalues of $A_n$ are approximately equal to

$${\lambda }_i\left(f\left(a+j\frac{b-a}{n}\right)\right),j=1,\dots ,n,i=1,\dots ,r,$$

for n large enough; similarly, if $k=2$, $d_n=n^{2}r$, and $D=[a_1,b_1]\times [a_2,b_2]$, then, assuming we have no outliers, the eigenvalues of $A_n$ are approximately equal to

$${\lambda }_i\left(f\left(a_1+j_1\frac{b_1-a_1}{n},a_2+j_2\frac{b_2-a_2}{n}\right)\right),j_1,j_2=1,\dots ,n,i=1,\dots ,r,$$

for n large enough; and so on for $k\ge 3$. A completely analogous meaning can also be given for the singular value distribution (1).

Remark 2.

Let $D=[a_1,b_1]\times \cdots \times [a_k,b_k]\subset {\mathbb{R}}^k$ and let $f:D\to {\mathbb{C}}^{r\times r}$ be a measurable function possessing r real-valued Riemann-integrable eigenvalue functions ${\lambda }_i\left(f\left(\mathbf{x}\right)\right)$, $i=1,\dots ,r$. Compute for each $\rho \in \mathbb{N}$ the uniform samples

$${\lambda }_i\left(f\left(a_1+j_1\frac{b_1-a_1}{\rho },\dots ,a_k+j_k\frac{b_k-a_k}{\rho }\right)\right),j_1,\dots ,j_k=1,\dots ,\rho ,i=1,\dots ,r,$$

sort them in non-decreasing order and put them in a vector $({\varsigma }_1,{\varsigma }_2,\dots ,{\varsigma }_{r{\rho }^k})$. Let ${\varphi }_{\rho }:[0,1]\to \mathbb{R}$ be the piecewise linear non-decreasing function that interpolates the samples $({\varsigma }_0={\varsigma }_1,{\varsigma }_1,{\varsigma }_2,\dots ,{\varsigma }_{r{\rho }^k})$ over the nodes $(0,\frac{1}{r{\rho }^k},\frac{2}{r{\rho }^k},\dots ,1)$, i.e.,

$$\left\{\begin{array}{c}{\varphi }_{\rho }\left(\frac{i}{r{\rho }^k}\right)={\varsigma }_i,i=0,\dots ,r{\rho }^k,\hfill \\ {\varphi }_{\rho } linear on \left[\frac{i}{r{\rho }^k},\frac{i+1}{r{\rho }^k}\right] for i=0,\dots ,r{\rho }^k-1.\hfill \end{array}\right.$$

Suppose ${\varphi }_{\rho }$ converges in measure over $[0,1]$ to some function ϕ as $\rho \to \infty$ (this is always the case in real-world applications). Then,

$${\int }_0^{1}F\left(\varphi \left(t\right)\right)\mathrm{d}t=\frac{1}{{\mu }_k\left(D\right)}{\int }_D\frac{{\sum }_{i=1}^{r}F\left({\lambda }_i\left(f\left(\mathbf{x}\right)\right)\right)}{r}\mathrm{d}\mathbf{x},\forall F\in C_c\left(\mathbb{C}\right).$$

(3)

This result can be proved by adapting the argument used in ([13] solution of Exercise 3.1). The function ϕ is referred to as the canonical rearranged version of f. What is interesting about ϕ is that, by (3), if ${\left\{A_n\right\}}_n{\sim }_{\lambda }f$ then ${\left\{A_n\right\}}_n{\sim }_{\lambda }\varphi$, i.e., if f is a spectral symbol of ${\left\{A_n\right\}}_n$ then the same is true for ϕ. Moreover, ϕ is a univariate scalar function and hence it is much easier to handle than f. According to Remark 1, assuming that ϕ is Riemann-integrable, if we have ${\left\{A_n\right\}}_n{\sim }_{\lambda }f$ (and hence also ${\left\{A_n\right\}}_n{\sim }_{\lambda }\varphi$), then, for n large enough, the eigenvalues of $A_n$, with the possible exception of $o\left(d_n\right)$ outliers, are approximately equal to the samples of ϕ over a uniform grid in $[0,1]$.

The next two theorems are useful tools for computing the spectral distribution of sequences formed by Hermitian or perturbed Hermitian matrices. For the related proofs, we refer the reader to ([38] Theorem 4.3) and ([39] Theorem 1.1). In the following, the conjugate transpose of the matrix A is denoted by $A^{*}$. If $A\in {\mathbb{C}}^{m\times m}$ and $1\le p\le \infty$, we denote by ${\parallel A\parallel }_p$ the Schatten p-norm of A, i.e., the p-norm of the vector $({\sigma }_1\left(A\right),\dots ,{\sigma }_m\left(A\right))$. The Schatten ∞-norm ${\parallel A\parallel }_{\infty }$ is the largest singular value of A and coincides with the spectral norm $\parallel A\parallel$. The Schatten 1-norm ${\parallel A\parallel }_1$ is the sum of the singular values of A and is often referred to as the trace-norm of A. The Schatten 2-norm ${\parallel A\parallel }_2$ coincides with the Frobenius norm of A. For more on Schatten p-norms, see [40].

Theorem 1.

Let ${\left\{X_n\right\}}_n$ be a sequence of matrices, with $X_n$ Hermitian of size $d_n$, and let ${\left\{P_n\right\}}_n$ be a sequence such that $P_n\in {\mathbb{C}}^{d_n\times {\delta }_n}$, $P_n^{*}{P}_n=I_{{\delta }_n}$, ${\delta }_n\le d_n$ and ${\delta }_n/d_n\to 1$ as $n\to \infty$. Then, ${\left\{X_n\right\}}_n{\sim }_{\sigma ,\lambda }\kappa$ if and only if ${\left\{P_n^{*}{X}_n{P}_n\right\}}_n{\sim }_{\sigma ,\lambda }\kappa$.

Theorem 2.

Let ${\left\{X_n\right\}}_n$ and ${\left\{Y_n\right\}}_n$ be sequences of matrices, with $X_n$ and $Y_n$ of size $d_n$. Assume that:

• the matrices $X_n$ are Hermitian and ${\left\{X_n\right\}}_n{\sim }_{\lambda }\kappa$;
• $\parallel Y_n{\parallel }_2=o\left(\sqrt{d_n}\right)$;

then ${\{X_n+Y_n\}}_n{\sim }_{\lambda }\kappa$.

Block Toeplitz Matrices. Given a function $f:[-\pi ,\pi ]\to {\mathbb{C}}^{s\times s}$ in $L^1\left([-\pi ,\pi ]\right)$, its Fourier coefficients are denoted by

$$f_k=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }f\left(\theta \right){\mathrm{e}}^{-\mathrm{i}k\theta }\mathrm{d}\theta \in {\mathbb{C}}^{s\times s},k\in \mathbb{Z},$$

where the integrals are computed componentwise. The nth block Toeplitz matrix generated by f is defined as

$$T_n\left(f\right)={\left[f_{i-j}\right]}_{i,j=1}^n\in {\mathbb{C}}^{sn\times sn}.$$

It is not difficult to see that all the matrices $T_n\left(f\right)$ are Hermitian when f is Hermitian a.e.

Block Diagonal Sampling Matrices. For $n\in \mathbb{N}$ and $a:[0,1]\to {\mathbb{C}}^{s\times s}$, we define the block diagonal sampling matrix $D_n\left(a\right)$ as the diagonal matrix

$$D_n\left(a\right)=\underset{i=1,\dots ,n}{\mathrm{diag}}a\left(\frac{i}{n}\right)=\left[\begin{array}{cccc}a\left(\frac{1}{n}\right)& & & \\ & a\left(\frac{2}{n}\right)& & \\ & & \ddots & \\ & & & a\left(1\right)\end{array}\right]\in {\mathbb{C}}^{sn\times sn}.$$

Zero-Distributed Sequences. A sequence of matrices ${\left\{Z_n\right\}}_n$ such that ${\left\{Z_n\right\}}_n{\sim }_{\sigma }0$ is referred to as a zero-distributed sequence. Note that, for any $r\ge 1$, ${\left\{Z_n\right\}}_n{\sim }_{\sigma }0$ is equivalent to ${\left\{Z_n\right\}}_n{\sim }_{\sigma }{O}_r$ (throughout this paper, $O_m$ and $I_m$ denote the $m\times m$ zero matrix and the $m\times m$ identity matrix, respectively). Proposition 1 provides an important characterization of zero-distributed sequences together with a useful sufficient condition for detecting such sequences. Throughout this paper, we use the natural convention $1/\infty =0$.

Proposition 1.

Let ${\left\{Z_n\right\}}_n$ be a sequence of matrices, with $Z_n$ of size $d_n$.

• ${\left\{Z_n\right\}}_n$ is zero-distributed if and only if $Z_n=R_n+N_n$ with $\mathrm{rank}\left(R_n\right)/d_n\to 0$ and $\parallel N_n\parallel \to 0$.
• ${\left\{Z_n\right\}}_n$ is zero-distributed if there exists a $p\in [1,\infty ]$ such that $\parallel Z_n{\parallel }_p/{\left(d_n\right)}^{1/p}\to 0$.

Approximating Classes of Sequences. The notion of approximating classes of sequences (a.c.s.) is the fundamental concept on which the theory of block GLT sequences is based.

Definition 2.

Let ${\left\{A_n\right\}}_n$ be a sequence of matrices, with $A_n$ of size $d_n$, and let ${\left\{{\left\{B_{n,m}\right\}}_n\right\}}_m$ be a sequence of sequences of matrices, with $B_{n,m}$ of size $d_n$. We say that ${\left\{{\left\{B_{n,m}\right\}}_n\right\}}_m$ is an approximating class of sequences (a.c.s.) for ${\left\{A_n\right\}}_n$ if the following condition is met: for every m there exists $n_m$ such that, for $n\ge n_m$,

$$A_n=B_{n,m}+R_{n,m}+N_{n,m},\mathrm{rank}\left(R_{n,m}\right)\le c\left(m\right)d_n,\parallel N_{n,m}\parallel \le \omega \left(m\right),$$

where $n_m,c\left(m\right),\omega \left(m\right)$ depend only on m and ${\displaystyle \underset{m\to \infty }{\lim }c\left(m\right)=\underset{m\to \infty }{\lim }\omega \left(m\right)=0}$.

Roughly speaking, ${\left\{{\left\{B_{n,m}\right\}}_n\right\}}_m$ is an a.c.s. for ${\left\{A_n\right\}}_n$ if, for large m, the sequence ${\left\{B_{n,m}\right\}}_n$ approximates ${\left\{A_n\right\}}_n$ in the sense that $A_n$ is eventually equal to $B_{n,m}$ plus a small-rank matrix (with respect to the matrix size $d_n$) plus a small-norm matrix. It turns out that, for each fixed sequence of positive integers $d_n$ such that $d_n\to \infty$, the notion of a.c.s. is a notion of convergence in the space

$$\mathcal{E}=\{{\left\{A_n\right\}}_n:A_n\in {\mathbb{C}}^{d_n\times d_n}\}.$$

More precisely, there exists a pseudometric $d_{\mathrm{a}.\mathrm{c}.\mathrm{s}.}$ in $\mathcal{E}$ such that ${\left\{{\left\{B_{n,m}\right\}}_n\right\}}_m$ is an a.c.s. for ${\left\{A_n\right\}}_n$ if and only if $d_{\mathrm{a}.\mathrm{c}.\mathrm{s}.}({\left\{B_{n,m}\right\}}_n,{\left\{A_n\right\}}_n)\to 0$ as $m\to \infty$. We therefore use the convergence notation ${\left\{B_{n,m}\right\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\left\{A_n\right\}}_n$ to indicate that ${\left\{{\left\{B_{n,m}\right\}}_n\right\}}_m$ is an a.c.s. for ${\left\{A_n\right\}}_n$. A useful criterion to identify an a.c.s. is provided in the next proposition ([13] Corollary 5.3).

Proposition 2.

Let ${\left\{A_n\right\}}_n$ be a sequence of matrices, with $A_n$ of size $d_n$, let ${\left\{{\left\{B_{n,m}\right\}}_n\right\}}_m$ be a sequence of sequences of matrices, with $B_{n,m}$ of size $d_n$, and let $p\in [1,\infty ]$. Suppose that for every m there exists $n_m$ such that, for $n\ge n_m$,

$$\parallel A_n-B_{n,m}{\parallel }_p\le \epsilon (m,n){\left(d_n\right)}^{1/p},$$

where ${\displaystyle \underset{m\to \infty }{\lim }\underset{n\to \infty }{\lim \sup } \epsilon (m,n)=0}$. Then, ${\left\{B_{n,m}\right\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\left\{A_n\right\}}_n$.

If $X\in {\mathbb{C}}^{m_1\times m_2}$ and $Y\in {\mathbb{C}}^{{\ell }_1\times {\ell }_2}$ are any two matrices, the tensor (Kronecker) product of X and Y is the $m_1{\ell }_1\times m_2{\ell }_2$ matrix defined as follows:

$$X\otimes Y={\left[x_{ij}Y\right]}_{\begin{array}{c}i=1,\dots ,m_1\\ j=1,\dots ,m_2\end{array}}=\left[\begin{array}{ccc}{x}_{11}Y& \cdots & x_{1m_2}Y\\ ⋮& & ⋮\\ x_{m_{1}1}Y& \cdots & x_{m_1{m}_2}Y\end{array}\right].$$

We recall that the tensor product operation ⊗ is associative and bilinear. Moreover,

$$\begin{array}{cc}\hfill \parallel X\otimes Y\parallel & =\parallel X\parallel \parallel Y\parallel ,\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill \mathrm{rank}(X\otimes Y)& =\mathrm{rank}\left(X\right)\mathrm{rank}\left(Y\right),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill {(X\otimes Y)}^T& =X^T\otimes Y^T.\hfill \end{array}$$

(6)

Finally, if $X_1,X_2$ can be multiplied and $Y_1,Y_2$ can be multiplied, then

$$(X_1\otimes Y_1)(X_2\otimes Y_2)=\left(X_1{X}_2\right)\otimes \left(Y_1{Y}_2\right).$$

(7)

Lemma 1.

For $i,j=1,\dots ,s$, let ${\left\{A_{n,ij}\right\}}_n$ be a sequence of matrices and suppose that ${\left\{B_{n,ij}^{\left(m\right)}\right\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\left\{A_{n,ij}\right\}}_n$. Then,

$${\left[B_{n,ij}^{\left(m\right)}\right]}_{i,j=1}^s\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\left[A_{n,ij}\right]}_{i,j=1}^s.$$

Proof. 

Let $E_{ij}$ be the $s\times s$ matrix having 1 in position $(i,j)$ and 0 elsewhere. Note that

$${\left[A_{n,ij}\right]}_{i,j=1}^s=\sum _{i,j=1}^s{E}_{ij}\otimes A_{n,ij},{\left[B_{n,ij}^{\left(m\right)}\right]}_{i,j=1}^s=\sum _{i,j=1}^s{E}_{ij}\otimes B_{n,ij}^{\left(m\right)}.$$

(8)

Since ${\{B_{n,ij}^{\left(m\right)}\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\{A_{n,ij}\}}_n$, it is clear from (4), (5) and the definition of a.c.s. that

$$E_{ij}\otimes B_{n,ij}^{\left(m\right)}\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{E}_{ij}\otimes A_{n,ij},i,j=1,\dots ,s.$$

(9)

Now, if ${\{B_{n,m}^{\left[k\right]}\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\{A_n^{\left[k\right]}\}}_n$ for $k=1,\dots ,K$ then ${\{{\sum }_{k=1}^K{B}_{n,m}^{\left[k\right]}\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\{{\sum }_{k=1}^K{A}_n^{\left[k\right]}\}}_n$ (this is an obvious consequence of the definition of a.c.s.). Thus, the thesis follows from (8) and (9). ☐

Block GLT Sequences. Let $s\ge 1$ be a fixed positive integer. An s-block GLT sequence (or simply a GLT sequence if s can be inferred from the context or we do not need/want to specify it) is a special s-block matrix-sequence ${\left\{A_n\right\}}_n$ equipped with a measurable function $\kappa :[0,1]\times [-\pi ,\pi ]\to {\mathbb{C}}^{s\times s}$, the so-called symbol. We use the notation ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ to indicate that ${\left\{A_n\right\}}_n$ is a GLT sequence with symbol $\kappa$. The symbol of a GLT sequence is unique in the sense that if ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ and ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\varsigma$ then $\kappa =\varsigma$ a.e. in $[0,1]\times [-\pi ,\pi ]$. The main properties of s-block GLT sequences proved in [34] are listed below. If A is a matrix, we denote by $A^{†}$ the Moore–Penrose pseudoinverse of A (recall that $A^{†}=A^{-1}$ whenever A is invertible). If $f_m,f:D\subseteq {\mathbb{R}}^k\to {\mathbb{C}}^{r\times r}$ are measurable matrix-valued functions, we say that $f_m$ converges to f in measure (resp., a.e., in $L^p\left(D\right)$, etc.) if ${\left(f_m\right)}_{\alpha \beta }$ converges to $f_{\alpha \beta }$ in measure (resp., a.e., in $L^p\left(D\right)$, etc.) for all $\alpha ,\beta =1,\dots ,r$.

GLT 1.

If ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ then ${\left\{A_n\right\}}_n{\sim }_{\sigma }\kappa$. If, moreover, each $A_n$ is Hermitian, then ${\left\{A_n\right\}}_n{\sim }_{\lambda }\kappa$.

GLT 2.

We have:

• ${\left\{T_n\left(f\right)\right\}}_n{\sim }_{\mathrm{GLT}}\kappa (x,\theta )=f\left(\theta \right)$ if $f:[-\pi ,\pi ]\to {\mathbb{C}}^{s\times s}$ is in $L^1\left([-\pi ,\pi ]\right)$;
• ${\left\{D_n\left(a\right)\right\}}_n{\sim }_{\mathrm{GLT}}\kappa (x,\theta )=a\left(x\right)$ if $a:[0,1]\to {\mathbb{C}}^{s\times s}$ is Riemann-integrable;
• ${\left\{Z_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa (x,\theta )=O_s$ if and only if ${\left\{Z_n\right\}}_n{\sim }_{\sigma }0$.

GLT 3.

If ${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ and ${\left\{B_n\right\}}_n{\sim }_{\mathrm{GLT}}\varsigma$, then:

• ${\left\{A_n^{*}\right\}}_n{\sim }_{\mathrm{GLT}}{\kappa }^{*}$;
• ${\{\alpha A_n+\beta B_n\}}_n{\sim }_{\mathrm{GLT}}\alpha \kappa +\beta \varsigma$ for all $\alpha ,\beta \in \mathbb{C}$;
• ${\left\{A_n{B}_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa \varsigma$;
• ${\{A_n^{†}\}}_n{\sim }_{\mathrm{GLT}}{\kappa }^{-1}$ provided that $\kappa$ is invertible a.e.

GLT 4.

${\left\{A_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ if and only if there exist s-block GLT sequences ${\left\{B_{n,m}\right\}}_n{\sim }_{\mathrm{GLT}}{\kappa }_m$ such that ${\left\{B_{n,m}\right\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\left\{A_n\right\}}_n$ and ${\kappa }_m\to \kappa$ in measure.

Remark 3.

The reader might be astonished by the fact that we have talked so far about block GLT sequences without defining them. Actually, we intentionally avoided giving a definition for two reasons. First, the definition is rather cumbersome as it requires introducing other related (and complicated) concepts such as “block LT operators” and “block LT sequences”. Second, from a practical viewpoint, the definition is completely useless because everything that can be derived from it can also be derived fromGLT 1–GLT 4(and in a much easier way). The reader who is interested in the formal definition of block GLT sequences can find it in ([34] Section 5) along with the proofs of properties GLT 1–GLT 4.

3. FD Discretization of a System of DEs

Consider the following system of DEs:

$$\left\{\begin{array}{ccc}\hfill -a_{11}\left(x\right)u_1^{″}\left(x\right)+a_{12}\left(x\right)u_2^{\prime }\left(x\right)& =f_1\left(x\right),\hfill & \hfill x\in (0,1),\\ \hfill a_{21}\left(x\right)u_1^{\prime }\left(x\right)+a_{22}\left(x\right)u_2\left(x\right)& =f_2\left(x\right),\hfill & \hfill x\in (0,1),\\ \hfill u_1\left(0\right)=0,u_1\left(1\right)=0,\\ \hfill u_2\left(0\right)=0,u_2\left(1\right)=0.\end{array}\right.$$

(10)

In this section, we consider the classical central FD discretization of (10). Through the theory of block GLT sequences, we show that the corresponding sequence of (normalized) FD discretization matrices enjoys a spectral distribution described by a $2\times 2$ matrix-valued function. We remark that the number 2, which identifies the matrix space ${\mathbb{C}}^{2\times 2}$ where the spectral symbol takes values, coincides with the number of equations that compose the system (10). In what follows, we use the following notation:

$$\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccc}{\beta }_j& {\alpha }_j& {\gamma }_j\end{array}\right]=\left[\begin{array}{ccccc}{\alpha }_1& {\gamma }_1& & & \\ {\beta }_2& {\alpha }_2& {\gamma }_2& & \\ & \ddots & \ddots & \ddots & \\ & & {\beta }_{n-1}& {\alpha }_{n-1}& {\gamma }_{n-1}\\ & & & {\beta }_n& {\alpha }_n\end{array}\right].$$

The parameters ${\alpha }_j,{\beta }_j,{\gamma }_j$ may be either scalars or $s\times s$ blocks for some $s>1$, in which case the previous matrix is a block tridiagonal matrix.

3.1. FD Discretization

Let $n\ge 1$, and set $h=\frac{1}{n+1^{\phantom{0}}}$ and $x_j=jh$ for $j=0,\dots ,n+1$. Using the classical central FD schemes $(-1,2,-1)$ and $\frac{1}{2^{\phantom{0}}}(-1,0,1)$ for the discretization of, respectively, the (negative) second derivative and the first derivative, for each $j=1,\dots ,n$ we obtain the following approximations:

$$\begin{array}{cc}\hfill {\left.[-a_{11}\left(x\right)u_1^{″}\left(x\right)+a_{12}\left(x\right)u_2^{\prime }\left(x\right)]\right|}_{x=x_j}& \approx a_{11}\left(x_j\right)\frac{-u_1\left(x_{j+1}\right)+2u_1\left(x_j\right)-u_1\left(x_{j-1}\right)}{h^2}\hfill \\ & +a_{12}\left(x_j\right)\frac{u_2\left(x_{j+1}\right)-u_2\left(x_{j-1}\right)}{2h},\hfill \\ \hfill {\left.[a_{21}\left(x\right)u_1^{\prime }\left(x\right)+a_{22}\left(x\right)u_2\left(x\right)]\right|}_{x=x_j}& \approx a_{21}\left(x_j\right)\frac{u_1\left(x_{j+1}\right)-u_1\left(x_{j-1}\right)}{2h}+a_{22}\left(x_j\right)u_2\left(x_j\right).\hfill \end{array}$$

This means that the nodal values of the solutions $u_1,u_2$ of (10) satisfy approximately the equations

$$\begin{array}{cc}\hfill a_{11}\left(x_j\right)\left[-u_1\left(x_{j+1}\right)+2u_1\left(x_j\right)-u_1\left(x_{j-1}\right)\right]+\frac{h}{2}{a}_{12}\left(x_j\right)\left[u_2\left(x_{j+1}\right)-u_2\left(x_{j-1}\right)\right]& =h^2{f}_1\left(x_j\right),\hfill \\ \hfill \frac{1}{2}{a}_{21}\left(x_j\right)\left[u_1\left(x_{j+1}\right)-u_1\left(x_{j-1}\right)\right]+h{a}_{22}\left(x_j\right)u_2\left(x_j\right)& =h{f}_2\left(x_j\right),\hfill \end{array}$$

for $j=1,\dots ,n$. We then approximate the solution $u_1$ (resp., $u_2$) by the piecewise linear function that takes the value $u_{1,j}$ (resp., $u_{2,j}$) at $x_j$ for all $j=0,\dots ,n+1$, where $u_{1,0}=u_{1,n+1}=u_{2,0}=u_{2,n+1}=0$ and the vectors ${\mathbf{u}}_1={(u_{1,1},\dots ,u_{1,n})}^T$ and ${\mathbf{u}}_2={(u_{2,1},\dots ,u_{2,n})}^T$ solve the linear system

$$\begin{array}{cc}\hfill a_{11}\left(x_j\right)\left[-u_{1,j+1}+2u_{1,j}-u_{1,j-1}\right]+\frac{h}{2}{a}_{12}\left(x_j\right)\left[u_{2,j+1}-u_{2,j-1}\right]& =h^2{f}_1\left(x_j\right),j=1,\dots ,n,\hfill \\ \hfill \frac{1}{2}{a}_{21}\left(x_j\right)\left[u_{1,j+1}-u_{1,j-1}\right]+h{a}_{22}\left(x_j\right)u_{2,j}& =h{f}_2\left(x_j\right),j=1,\dots ,n.\hfill \end{array}$$

(11)

This linear system can be rewritten in matrix form as follows:

$$A_n\left[\begin{array}{c}{\mathbf{u}}_1\\ {\mathbf{u}}_2\end{array}\right]=\left[\begin{array}{c}{h}^2{\mathbf{f}}_1\\ h{\mathbf{f}}_2\end{array}\right],$$

(12)

where ${\mathbf{f}}_1={\left[f_1\left(x_j\right)\right]}_{j=1}^n,{\mathbf{f}}_2={\left[f_2\left(x_j\right)\right]}_{j=1}^n$,

$$\begin{array}{cc}\hfill A_n& =\left[\begin{array}{cc}{K}_n\left(a_{11}\right)& h{H}_n\left(a_{12}\right)\\ H_n\left(a_{21}\right)& h{M}_n\left(a_{22}\right)\end{array}\right]=\left[\begin{array}{cc}{K}_n\left(a_{11}\right)& H_n\left(a_{12}\right)\\ H_n\left(a_{21}\right)& M_n\left(a_{22}\right)\end{array}\right]\left[\begin{array}{cc}{I}_n& O_n\\ O_n& h{I}_n\end{array}\right],\hfill \end{array}$$

(13)

and

$$\begin{array}{cc}\hfill K_n\left(a_{11}\right)& =\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}-a_{11}\left(x_j\right)& |& 2a_{11}\left(x_j\right)& |& -a_{11}\left(x_j\right)\end{array}\right]=\left(\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{11}\left(x_j\right)\right)T_n(2-2\cos \theta ),\hfill \\ \hfill H_n\left(a_{12}\right)& =\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}-\frac{1}{2^{\phantom{0}}}{a}_{12}\left(x_j\right)& |& 0& |& \frac{1}{2^{\phantom{0}}}{a}_{12}\left(x_j\right)\end{array}\right]=\left(\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{12}\left(x_j\right)\right)T_n(-\mathrm{i}\sin \theta ),\hfill \\ \hfill H_n\left(a_{21}\right)& =\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}-\frac{1}{2^{\phantom{0}}}{a}_{21}\left(x_j\right)& |& 0& |& \frac{1}{2^{\phantom{0}}}{a}_{21}\left(x_j\right)\end{array}\right]=\left(\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{21}\left(x_j\right)\right)T_n(-\mathrm{i}\sin \theta ),\hfill \\ \hfill M_n\left(a_{22}\right)& =\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{22}\left(x_j\right).\hfill \end{array}$$

In view of (13), the linear system (12) is equivalent to

$$B_n\left[\begin{array}{c}{\mathbf{v}}_1\\ {\mathbf{v}}_2\end{array}\right]=\left[\begin{array}{c}{h}^2{\mathbf{f}}_1\\ h{\mathbf{f}}_2\end{array}\right],$$

(14)

where ${\mathbf{v}}_1={\mathbf{u}}_1,{\mathbf{v}}_2=h{\mathbf{u}}_2$, and

$$B_n=\left[\begin{array}{cc}{K}_n\left(a_{11}\right)& H_n\left(a_{12}\right)\\ H_n\left(a_{21}\right)& M_n\left(a_{22}\right)\end{array}\right].$$

(15)

Let $v_{1,1},\dots ,v_{1,n}$ and $v_{2,1},\dots ,v_{2,n}$ be the components of ${\mathbf{v}}_1$ and ${\mathbf{v}}_2$, respectively. When writing the linear system (11) in the form (14), we are implicitly assuming the following.

• The unknowns are sorted as follows:

  $$\left[\begin{array}{c}{\left[v_{1,j}\right]}_{j=1,\dots ,n}\\ {\left[v_{2,j}\right]}_{j=1,\dots ,n}\end{array}\right]=\left[\begin{array}{c}{v}_{1,1}\\ v_{1,2}\\ ⋮\\ v_{1,n}\\ v_{2,1}\\ v_{2,2}\\ ⋮\\ v_{2,n}\end{array}\right].$$

  (16)
• The equations are sorted as follows, in accordance with the ordering (16) for the unknowns:

  $$\left[\frac{{\left[a_{11}\left(x_j\right)\left[-v_{1,j+1}+2v_{1,j}-v_{1,j-1}\right]+\frac{1}{2^{\phantom{0}}}{a}_{12}\left(x_j\right)\left[v_{2,j+1}-v_{2,j-1}\right]=h^2{f}_1\left(x_j\right)\right]}_{j=1,\dots ,n}}{{\left[\frac{1}{2^{\phantom{0}}}{a}_{21}\left(x_j\right)\left[u_{1,j+1}-u_{1,j-1}\right]+a_{22}\left(x_j\right)v_{2,j}=h{f}_2\left(x_j\right)\right]}_{j=1,\dots ,n}}\right].$$

  (17)

Suppose we decide to change the ordering for both the unknowns and the equations. More precisely, suppose we opt for the following orderings.

• The unknowns are sorted as follows:

  $${\left[\begin{array}{c}{v}_{1,j}\\ v_{2,j}\end{array}\right]}_{j=1,\dots ,n}=\left[\begin{array}{c}{v}_{1,1}\\ v_{2,1}\\ v_{1,2}\\ v_{2,2}\\ ⋮\\ v_{1,n}\\ v_{2,n}\end{array}\right].$$

  (18)
• The equations are sorted as follows, in accordance with the ordering (18) for the unknowns:

  $${\left[\begin{array}{c}{a}_{11}\left(x_j\right)\left[-v_{1,j+1}+2v_{1,j}-v_{1,j-1}\right]+\frac{1}{2^{\phantom{0}}}{a}_{12}\left(x_j\right)\left[v_{2,j+1}-v_{2,j-1}\right]=h^2{f}_1\left(x_j\right)\\ \frac{1}{2^{\phantom{0}}}{a}_{21}\left(x_j\right)\left[v_{1,j+1}-v_{1,j-1}\right]+a_{22}\left(x_j\right)v_{2,j}=h{f}_2\left(x_j\right)\end{array}\right]}_{j=1,\dots ,n}.$$

  (19)

The matrix $C_n$ associated with the linear system (11) assuming the new orderings (18) and (19) is the $2\times 2$ block tridiagonal matrix given by

$$C_n=\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}-a_{11}\left(x_j\right)& -\frac{1}{2^{\phantom{0}}}{a}_{12}\left(x_j\right)\\ -\frac{1}{2^{\phantom{0}}}{a}_{21}\left(x_j\right)& 0\end{array}& |& \begin{array}{cc}2a_{11}\left(x_j\right)& 0\\ 0& a_{22}\left(x_j\right)\end{array}& |& \begin{array}{cc}-a_{11}\left(x_j\right)& \frac{1}{2^{\phantom{0}}}{a}_{12}\left(x_j\right)\\ \frac{1}{2^{\phantom{0}}}{a}_{21}\left(x_j\right)& 0\end{array}\end{array}\right].$$

(20)

The matrix $C_n$ is similar to $B_n$. Indeed, by permuting both rows and columns of $B_n$ according to the permutation $1,n+1,2,n+2,\dots ,n,2n$ we obtain $C_n$. More precisely, let ${\mathbf{e}}_1,\dots ,{\mathbf{e}}_n$ and ${\tilde{\mathbf{e}}}_1,\dots ,{\tilde{\mathbf{e}}}_{2n}$ be the vectors of the canonical basis of ${\mathbb{R}}^n$ and ${\mathbb{R}}^{2n}$, respectively, and let ${\mathsf{\Pi }}_n$ be the permutation matrix associated with the permutation $1,n+1,2,n+2,\dots ,n,2n$, that is,

$${\mathsf{\Pi }}_n=\left[\begin{array}{c}{\tilde{\mathbf{e}}}_1^T\\ {\tilde{\mathbf{e}}}_{n+1}^T\\ {\tilde{\mathbf{e}}}_2^T\\ {\tilde{\mathbf{e}}}_{n+2}^T\\ ⋮\\ {\tilde{\mathbf{e}}}_n^T\\ {\tilde{\mathbf{e}}}_{2n}^T\end{array}\right]=\left[\begin{array}{c}{I}_2\otimes {\mathbf{e}}_1^T\\ I_2\otimes {\mathbf{e}}_2^T\\ ⋮\\ I_2\otimes {\mathbf{e}}_n^T\end{array}\right].$$

(21)

Then, $C_n={\mathsf{\Pi }}_n{B}_n{\mathsf{\Pi }}_n^T$.

3.2. GLT Analysis of the FD Discretization Matrices

The main result of this section (Theorem 3) shows that ${\left\{C_n\right\}}_n$ is a block GLT sequence whose spectral distribution is described by a $2\times 2$ matrix-valued symbol, which is obtained by replacing the matrix-sequences ${\left\{K_n\left(a_{11}\right)\right\}}_n$, ${\left\{H_n\left(a_{12}\right)\right\}}_n$, ${\left\{H_n\left(a_{21}\right)\right\}}_n$, ${\left\{M_n\left(a_{22}\right)\right\}}_n$ appearing in the expression (15) of $B_n$ with the corresponding symbols $a_{11}\left(x\right)(2-2\cos \theta )$, $-\mathrm{i}{a}_{12}\left(x\right)\sin \theta$, $-\mathrm{i}{a}_{21}\left(x\right)\sin \theta$, $a_{22}\left(x\right)$. In this regard, we note that, assuming for instance $a_{11},a_{12},a_{21},a_{22}\in C\left([0,1]\right)$, we have

$$\begin{array}{cc}\hfill {\left\{K_n\left(a_{11}\right)\right\}}_n& {\sim }_{\mathrm{GLT}}{a}_{11}\left(x\right)(2-2\cos \theta ),\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill {\left\{H_n\left(a_{12}\right)\right\}}_n& {\sim }_{\mathrm{GLT}}-\mathrm{i}{a}_{12}\left(x\right)\sin \theta ,\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill {\left\{H_n\left(a_{21}\right)\right\}}_n& {\sim }_{\mathrm{GLT}}-\mathrm{i}{a}_{21}\left(x\right)\sin \theta ,\hfill \end{array}$$

(24)

$$\begin{array}{cc}\hfill {\left\{M_n\left(a_{22}\right)\right\}}_n& {\sim }_{\mathrm{GLT}}{a}_{22}\left(x\right).\hfill \end{array}$$

(25)

To prove (22), it suffices to observe that

$$\begin{array}{cc}\hfill \parallel K_n\left(a_{11}\right)-D_n\left(a_{11}\right)T_n(2-2\cos \theta )\parallel & \le ∥\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{11}\left(x_j\right)-D_n\left(a_{11}\right)∥\parallel T_n(2-2\cos \theta )\parallel \hfill \\ & =\underset{j=1,\dots ,n}{\max }\left|a_{11}\left(x_j\right)-a_{11}\left(\frac{j}{n}\right)\right|\parallel T_n(2-2\cos \theta )\parallel \le 4{\omega }_{a_{11}}\left(h\right),\hfill \end{array}$$

where ${\omega }_{a_{11}}(·)$ is the modulus of continuity of $a_{11}$. Since ${\omega }_{a_{11}}\left(h\right)\to 0$ as $n\to \infty$, it follows from Proposition 1 that ${\{K_n\left(a_{11}\right)-D_n\left(a_{11}\right)T_n(2-2\cos \theta )\}}_n{\sim }_{\sigma }0$, and so GLT 2 and GLT 3 immediately yield (22). The relations (23)–(25) are proved in the same way.

Theorem 3.

Suppose that $a_{11},a_{12},a_{21},a_{22}\in C\left([0,1]\right)$. Then,

$$\begin{array}{cc}\hfill {\left\{C_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa (x,\theta )& =\left[\begin{array}{cc}{a}_{11}\left(x\right)(2-2\cos \theta )& -\mathrm{i}{a}_{12}\left(x\right)\sin \theta \\ -\mathrm{i}{a}_{21}\left(x\right)\sin \theta & a_{22}\left(x\right)\end{array}\right]\hfill \end{array}$$

(26)

and

$${\left\{C_n\right\}}_n{\sim }_{\sigma }\kappa (x,\theta ).$$

(27)

If, moreover, $a_{21}=-a_{12}$, then we also have

$${\left\{C_n\right\}}_n{\sim }_{\lambda }\kappa (x,\theta ).$$

(28)

Proof. 

From (20), we have

$$\begin{array}{cc}\hfill C_n& =\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}-a_{11}\left(x_j\right)& -\frac{1}{2^{\phantom{0}}}{a}_{12}\left(x_j\right)\\ -\frac{1}{2^{\phantom{0}}}{a}_{21}\left(x_j\right)& 0\end{array}& |& \begin{array}{cc}2a_{11}\left(x_j\right)& 0\\ 0& a_{22}\left(x_j\right)\end{array}& |& \begin{array}{cc}-a_{11}\left(x_j\right)& \frac{1}{2^{\phantom{0}}}{a}_{12}\left(x_j\right)\\ \frac{1}{2^{\phantom{0}}}{a}_{21}\left(x_j\right)& 0\end{array}\end{array}\right]\hfill \\ & =\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}-a_{11}\left(x_j\right)& 0\\ 0& 0\end{array}& |& \begin{array}{cc}2a_{11}\left(x_j\right)& 0\\ 0& 0\end{array}& |& \begin{array}{cc}-a_{11}\left(x_j\right)& 0\\ 0& 0\end{array}\end{array}\right]\hfill \\ & +\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}0& -\frac{1}{2^{\phantom{0}}}{a}_{12}\left(x_j\right)\\ 0& 0\end{array}& |& \begin{array}{cc}0& 0\\ 0& 0\end{array}& |& \begin{array}{cc}0& \frac{1}{2^{\phantom{0}}}{a}_{12}\left(x_j\right)\\ 0& 0\end{array}\end{array}\right]\hfill \\ & +\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}0& 0\\ -\frac{1}{2^{\phantom{0}}}{a}_{21}\left(x_j\right)& 0\end{array}& |& \begin{array}{cc}0& 0\\ 0& 0\end{array}& |& \begin{array}{cc}0& 0\\ \frac{1}{2^{\phantom{0}}}{a}_{21}\left(x_j\right)& 0\end{array}\end{array}\right]\hfill \\ & +\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}0& 0\\ 0& 0\end{array}& |& \begin{array}{cc}0& 0\\ 0& a_{22}\left(x_j\right)\end{array}& |& \begin{array}{cc}0& 0\\ 0& 0\end{array}\end{array}\right]\hfill \\ & =\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{11}\left(x_j\right)I_2·\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}-1& 0\\ 0& 0\end{array}& |& \begin{array}{cc}2& 0\\ 0& 0\end{array}& |& \begin{array}{cc}-1& 0\\ 0& 0\end{array}\end{array}\right]\hfill \\ & +\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{12}\left(x_j\right)I_2·\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}0& -\frac{1}{2^{\phantom{0}}}\\ 0& 0\end{array}& |& \begin{array}{cc}0& 0\\ 0& 0\end{array}& |& \begin{array}{cc}0& \frac{1}{2^{\phantom{0}}}\\ 0& 0\end{array}\end{array}\right]\hfill \\ & +\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{21}\left(x_j\right)I_2·\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}0& 0\\ -\frac{1}{2^{\phantom{0}}}& 0\end{array}& |& \begin{array}{cc}0& 0\\ 0& 0\end{array}& |& \begin{array}{cc}0& 0\\ \frac{1}{2^{\phantom{0}}}& 0\end{array}\end{array}\right]\hfill \\ & +\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{22}\left(x_j\right)I_2·\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}0& 0\\ 0& 0\end{array}& |& \begin{array}{cc}0& 0\\ 0& 1\end{array}& |& \begin{array}{cc}0& 0\\ 0& 0\end{array}\end{array}\right]\hfill \\ & =\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{11}\left(x_j\right)I_2·T_n\left((2-2\cos \theta )E_{11}\right)\hfill \\ & +\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{12}\left(x_j\right)I_2·T_n\left((-\mathrm{i}\sin \theta )E_{12}\right)\hfill \\ & +\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{21}\left(x_j\right)I_2·T_n\left((-\mathrm{i}\sin \theta )E_{21}\right)\hfill \\ & +\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{22}\left(x_j\right)I_2·T_n\left(E_{22}\right),\hfill \end{array}$$

(29)

where $E_{pq}$ is the $2\times 2$ matrix having 1 in position $(p,q)$ and 0 elsewhere. It is clear that, for every $p,q=1,2$,

$$∥\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{pq}\left(x_j\right)I_2-D_n\left(a_{pq}{I}_2\right)∥\le {\omega }_{a_{pq}}\left(h\right)\to 0$$

as $n\to \infty$; hence, by Proposition 1, GLT 2 and GLT 3,

$${\left\{\underset{j=1,\dots ,n}{\mathrm{diag}}{a}_{pq}\left(x_j\right)I_2\right\}}_n{\sim }_{\mathrm{GLT}}a\left(x\right)I_2,p,q=1,2.$$

Consequently, the decomposition (29), GLT 2 and GLT 3 imply (26), which in turn implies (27) by GLT 1. It only remains to prove (28) in the case where $a_{21}=-a_{12}$. In this case, we have

$$C_n=\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}-a_{11}\left(x_j\right)& -\frac{1}{2^{\phantom{1}}}{a}_{12}\left(x_j\right)\\ \frac{1}{2^{\phantom{1}}}{a}_{12}\left(x_j\right)& 0\end{array}& |& \begin{array}{cc}2a_{11}\left(x_j\right)& 0\\ 0& a_{22}\left(x_j\right)\end{array}& |& \begin{array}{cc}-a_{11}\left(x_j\right)& \frac{1}{2^{\phantom{1}}}{a}_{12}\left(x_j\right)\\ -\frac{1}{2^{\phantom{1}}}{a}_{12}\left(x_j\right)& 0\end{array}\end{array}\right].$$

Consider the symmetric approximation of $C_n$ given by

$${\tilde{C}}_n=\underset{j=1,\dots ,n}{\mathrm{tridiag}}\left[\begin{array}{ccccc}\begin{array}{cc}-a_{11}\left(x_{j-1}\right)& -\frac{1}{2^{\phantom{1}}}{a}_{12}\left(x_{j-1}\right)\\ \frac{1}{2^{\phantom{1}}}{a}_{12}\left(x_{j-1}\right)& 0\end{array}& |& \begin{array}{cc}2a_{11}\left(x_j\right)& 0\\ 0& a_{22}\left(x_j\right)\end{array}& |& \begin{array}{cc}-a_{11}\left(x_j\right)& \frac{1}{2^{\phantom{1}}}{a}_{12}\left(x_j\right)\\ -\frac{1}{2^{\phantom{1}}}{a}_{12}\left(x_j\right)& 0\end{array}\end{array}\right].$$

It is not difficult to see that $\parallel C_n-{\tilde{C}}_n\parallel \to 0$ as $n\to \infty$ by invoking the inequality

$$\parallel X\parallel \le \sqrt{\left({\max }_{i=1,\dots ,n}\sum _{j=1}^n\left|x_{ij}\right|\right)\left({\max }_{j=1,\dots ,n}\sum _{i=1}^n\left|x_{ij}\right|\right)},X\in {\mathbb{C}}^{n\times n};$$

(30)

see, e.g., ([13] Section 2.4.1). Therefore:

• in view of the decomposition ${\tilde{C}}_n=C_n+({\tilde{C}}_n-C_n)$, we have ${\left\{{\tilde{C}}_n\right\}}_n{\sim }_{\mathrm{GLT}}\kappa (x,\theta )$ by (26), Proposition 1, GLT 2 and GLT 3, so in particular ${\left\{{\tilde{C}}_n\right\}}_n{\sim }_{\lambda }\kappa (x,\theta )$ by GLT 1 as ${\tilde{C}}_n$ is symmetric;
• $\parallel C_n-{\tilde{C}}_n{\parallel }_2\le \sqrt{n}\parallel C_n-{\tilde{C}}_n\parallel =o\left(\sqrt{2n}\right)$ as $n\to \infty$.

Thus, (28) follows from Theorem 2. ☐

Example 1.

Suppose that $a_{11},a_{12},a_{21},a_{22}\in C\left([0,1]\right)$ and $a_{21}=-a_{12}$, so that ${\left\{C_n\right\}}_n{\sim }_{\lambda }\kappa (x,\theta )$ by Theorem 3. The eigenvalue functions of $\kappa (x,\theta )$ are given by

$${\lambda }_{1,2}\left(\kappa (x,\theta )\right)=\frac{a_{11}\left(x\right)(2-2\cos \theta )+a_{22}\left(x\right)\pm \sqrt{{(a_{11}\left(x\right)(2-2\cos \theta )-a_{22}\left(x\right))}^2+4{\left(a_{12}\left(x\right)\sin \theta \right)}^2}}{2}$$

and are continuous on $[0,1]\times [-\pi ,\pi ]$. Let ϕ be the canonical rearranged version of $\kappa (x,\theta )$ obtained as the limit of the piecewise linear functions ${\varphi }_{\rho }$, according to the construction in Remark 2. Figure 1 shows the graph of ϕ and the eigenvalues ${\lambda }_1,\dots ,{\lambda }_{2n}$ of $C_n$ for $a_{11}\left(x\right)=2+\cos \left(\pi x\right)$, $a_{12}\left(x\right)=-a_{21}\left(x\right)={\mathrm{e}}^{-x}\sin \left(\pi x\right)$, $a_{22}\left(x\right)=2x+\sin \left(\pi x\right)$ and $n=40$. The graph of ϕ has been obtained by plotting the graph of ${\varphi }_{\rho }$ corresponding to a large value of ρ. The eigenvalues of $C_n$, which turn out to be real, although $C_n$ is not symmetric, have been sorted in non-decreasing order and placed at the points $(t_q,{\lambda }_q)$ with $t_q=\frac{q}{2n}$, $q=1,\dots ,2n$. We clearly see from the figure an excellent agreement between ϕ and the eigenvalues of $C_n$, as predicted by Remark 2. In particular, we observe no outliers in this case.

4. Higher-Order FE Discretization of the Diffusion Equation

Consider the diffusion equation

$$\left\{\begin{array}{cc}-{\left(a\left(x\right)u^{\prime }\left(x\right)\right)}^{\prime }=f\left(x\right),\hfill & x\in (0,1),\hfill \\ u\left(0\right)=u\left(1\right)=0.\hfill \end{array}\right.$$

(31)

In this section, we consider the higher-order FE discretization of (31). Through the theory of block GLT sequences, we show that the corresponding sequence of (normalized) FE discretization matrices enjoys a spectral distribution described by a $(p-k)\times (p-k)$ matrix-valued function, where p and k represent, respectively, the degree and the smoothness of the piecewise polynomial functions involved in the FE approximation. Note that this result represents a remarkable argument in support of ([35] Conjecture 2).

4.1. FE Discretization

The weak form of (31) reads as follows: find $u\in H_0^1\left([0,1]\right)$ such that

$${\int }_0^{1}a\left(x\right)u^{\prime }\left(x\right)w^{\prime }\left(x\right)\mathrm{d}x={\int }_0^{1}f\left(x\right)w\left(x\right)\mathrm{d}x,\forall w\in H_0^1\left([0,1]\right).$$

In the FE method, we fix a set of basis functions $\{{\phi }_1,\dots ,{\phi }_N\}\subset H_0^1\left([0,1]\right)$ and we look for an approximation of the exact solution in the space $\mathcal{W}=\mathrm{span}({\phi }_1,\dots ,{\phi }_N)$ by solving the following discrete problem: find $u_{\mathcal{W}}\in \mathcal{W}$ such that

$${\int }_0^{1}a\left(x\right)u_{\mathcal{W}}^{\prime }\left(x\right)w^{\prime }\left(x\right)\mathrm{d}x={\int }_0^{1}f\left(x\right)w\left(x\right)\mathrm{d}x,\forall w\in \mathcal{W}.$$

Since $\{{\phi }_1,\dots ,{\phi }_N\}$ is a basis of $\mathcal{W}$, we can write $u_{\mathcal{W}}={\sum }_{j=1}^N{u}_j{\phi }_j$ for a unique vector $\mathbf{u}={(u_1,\dots ,u_N)}^T$. By linearity, the computation of $u_{\mathcal{W}}$ (i.e., of $\mathbf{u}$) reduces to solving the linear system

$$A\mathbf{u}=\mathbf{f},$$

where $\mathbf{f}={\left({\int }_0^{1}f\left(x\right){\phi }_1\left(x\right)\mathrm{d}x,\dots ,{\int }_0^{1}f\left(x\right){\phi }_N\left(x\right)\mathrm{d}x\right)}^T$ and A is the stiffness matrix,

$$\begin{array}{c}\hfill A={\left[{\int }_0^{1}a\left(x\right){\phi }_j^{\prime }\left(x\right){\phi }_i^{\prime }\left(x\right)\mathrm{d}x\right]}_{i,j=1}^N.\end{array}$$

(32)

4.2. p-Degree $C^k$ B-spline Basis Functions

Following the higher-order FE approach, the basis functions ${\phi }_1,\dots ,{\phi }_N$ will be chosen as piecewise polynomials of degree $p\ge 1$. More precisely, for $p,n\ge 1$ and $0\le k\le p-1$, let $B_{1,[p,k]},\dots ,B_{n(p-k)+k+1,[p,k]}:\mathbb{R}\to \mathbb{R}$ be the B-splines of degree p and smoothness $C^k$ defined on the knot sequence

$$\begin{array}{c}\hfill \{{\tau }_1,\dots ,{\tau }_{n(p-k)+p+k+2}\}=\left\{\underset{p+1}{\underbrace{0,\dots ,0}},\underset{p-k}{\underbrace{\frac{1}{n},\dots ,\frac{1}{n}}},\underset{p-k}{\underbrace{\frac{2}{n},\dots ,\frac{2}{n}}},\dots ,\underset{p-k}{\underbrace{\frac{n-1}{n},\dots ,\frac{n-1}{n}}},\underset{p+1}{\underbrace{1,\dots ,1}}\right\}.\end{array}$$

(33)

We collect here a few properties of $B_{1,[p,k]},\dots ,B_{n(p-k)+k+1,[p,k]}$ that we shall use in this paper. For the formal definition of B-splines, as well as for the proof of the properties listed below, see [41,42].

• The support of the ith B-spline is given by

  $$\mathrm{supp}\left(B_{i,[p,k]}\right)=[{\tau }_i,{\tau }_{i+p+1}],i=1,\dots ,n(p-k)+k+1.$$

  (34)
• Except for the first and the last one, all the other B-splines vanish on the boundary of $[0,1]$, i.e.,

  $$B_{i,[p,k]}\left(0\right)=B_{i,[p,k]}\left(1\right)=0,i=2,\dots ,n(p-k)+k.$$

  (35)
• $\{B_{1,[p,k]},\dots ,B_{n(p-k)+k+1,[p,k]}\}$ is a basis for the space of piecewise polynomial functions of degree p and smoothness $C^k$, that is,

  $${\mathcal{V}}_{n,[p,k]}=\left\{v\in C^k{\left([0,1]\right):v|}_{\left[\frac{i}{n},\frac{i+1}{n}\right]}\in {\mathbb{P}}_p\mathrm{for} \mathrm{all}i=0,\dots ,n-1\right\},$$

  where ${\mathbb{P}}_p$ is the space of polynomials of degree $\le p$. Moreover, $\{B_{2,[p,k]},\dots ,B_{n(p-k)+k,[p,k]}\}$ is a basis for the space

  $${\mathcal{W}}_{n,[p,k]}=\{w\in {\mathcal{V}}_{n,[p,k]}:w\left(0\right)=w\left(1\right)=0\}.$$
• The B-splines form a non-negative partition of unity over $[0,1]$:

  $$\begin{array}{cc}\hfill B_{i,[p,k]}& \ge 0\mathrm{over}\mathbb{R},i=1,\dots ,n(p-k)+k+1,\hfill \end{array}$$

  (36)

  $$\begin{array}{cc}\hfill \sum _{i=1}^{n(p-k)+k+1}{B}_{i,[p,k]}& =1\mathrm{over}[0,1].\hfill \end{array}$$

  (37)
• The derivatives of the B-splines satisfy

  $$\sum _{i=1}^{n(p-k)+k+1}\left|B_{i,[p,k]}^{\prime }\right|\le c_{p}n\mathrm{over}[0,1],$$

  (38)

  where $c_p$ is a constant depending only on p. Note that the derivatives $B_{i,[p,k]}^{\prime }$ may not be defined at some of the grid points $0,\frac{1}{n},\frac{2}{n},\dots ,\frac{n-1}{n},1$ in the case of $C^0$ smoothness ($k=0$). In (38), it is assumed that the undefined values are excluded from the summation.
• All the B-splines, except for the first $k+1$ and the last $k+1$, are uniformly shifted-scaled versions of $p-k$ fixed reference functions ${\beta }_{1,[p,k]},\dots ,{\beta }_{p-k,[p,k]}$, namely the first $p-k$ B-splines defined on the reference knot sequence

  $$\underset{p-k}{\underbrace{0,\dots ,0}},\underset{p-k}{\underbrace{1,\dots ,1}},\dots ,\underset{p-k}{\underbrace{\eta ,\dots ,\eta }},\eta =⌈\frac{p+1}{p-k}⌉.$$
  In formulas, setting

  $$\nu =⌈\frac{k+1}{p-k}⌉,$$

  (39)

  for the B-splines $B_{k+2,[p,k]},\dots ,B_{k+1+(n-\nu )(p-k),[p,k]}$, we have

  $$\begin{array}{c}\hfill B_{k+1+(p-k)(r-1)+q,[p,k]}\left(x\right)={\beta }_{q,[p,k]}(nx-r+1),r=1,\dots ,n-\nu ,q=1,\dots ,p-k.\end{array}$$

  (40)

  We point out that the supports of the reference B-splines ${\beta }_{q,[p,k]}$ satisfy

  $$\mathrm{supp}\left({\beta }_{1,[p,k]}\right)\subseteq \mathrm{supp}\left({\beta }_{2,[p,k]}\right)\subseteq \dots \subseteq \mathrm{supp}\left({\beta }_{p-k,[p,k]}\right)=[0,\eta ].$$

  Figure 2 and Figure 3 show the graphs of the B-splines $B_{1,[p,k]},\dots ,B_{n(p-k)+k+1,[p,k]}$ for the degree $p=3$ and the smoothness $k=1$, and the graphs of the associated reference B-splines ${\beta }_{1,[p,k]},{\beta }_{2,[p,k]}$.

The basis functions ${\phi }_1,\dots ,{\phi }_N$ are defined as follows:

$$\begin{array}{cc}\hfill {\phi }_i& =B_{i+1,[p,k]},i=1,\dots ,n(p-k)+k-1.\hfill \end{array}$$

(41)

In particular, with the notations of Section 4.1, we have $N=n(p-k)+k-1$ and $\mathcal{W}={\mathcal{W}}_{n,[p,k]}$.

4.3. GLT Analysis of the Higher-Order FE Discretization Matrices

The stiffness matrix (32) resulting from the choice of the basis functions as in (41) will be denoted by $A_{n,[p,k]}\left(a\right)$,

$$A_{n,[p,k]}\left(a\right)={\left[{\int }_0^{1}a\left(x\right)B_{j+1,[p,k]}^{\prime }\left(x\right)B_{i+1,[p,k]}^{\prime }\left(x\right)\mathrm{d}x\right]}_{i,j=1}^{n(p-k)+k-1}.$$

(42)

The main result of this section (Theorem 4) gives the spectral distribution of the normalized sequence ${\left\{n^{-1}{A}_{n,[p,k]}\left(a\right)\right\}}_n$. The proof of Theorem 4 is entirely based on the theory of block GLT sequences and it is therefore referred to as “GLT analysis”. It also requires the following lemma, which provides an approximate construction of the matrix $A_{n,[p,k]}\left(1\right)$ corresponding to the constant-coefficient case where $a\left(x\right)=1$ identically. In view of what follows, define the $(p-k)\times (p-k)$ blocks

$$K_{[p,k]}^{[\ell ]}={\left[{\int }_{\mathbb{R}}{\beta }_{j,[p,k]}^{\prime }\left(t\right){\beta }_{i,[p,k]}^{\prime }(t-\ell )\mathrm{d}t\right]}_{i,j=1}^{p-k},\ell \in \mathbb{Z},$$

(43)

and the $(p-k)\times (p-k)$ matrix-valued function ${\kappa }_{[p,k]}:[-\pi ,\pi ]\to {\mathbb{C}}^{(p-k)\times (p-k)}$,

$${\kappa }_{[p,k]}\left(\theta \right)=\sum _{\ell \in \mathbb{Z}}{K}_{[p,k]}^{[\ell ]}{\mathrm{e}}^{\mathrm{i}\ell \theta }=K_{[p,k]}^{\left[0\right]}+\sum _{\ell >0}\left(K_{[p,k]}^{[\ell ]}{\mathrm{e}}^{\mathrm{i}\ell \theta }+{\left(K_{[p,k]}^{[\ell ]}\right)}^T{\mathrm{e}}^{-\mathrm{i}\ell \theta }\right).$$

(44)

Due to the compact support of the reference functions ${\beta }_{1,[p,k]},\dots ,{\beta }_{p-k,[p,k]}$, there is only a finite number of nonzero blocks $K_{[p,k]}^{[\ell ]}$ and, consequently, the series in (44) is actually a finite sum.

Lemma 2.

Let $p,n\ge 1$ and $0\le k\le p-1$. Define ${\tilde{A}}_{n,[p,k]}\left(1\right)$ as the principal submatrix of $A_{n,[p,k]}\left(1\right)$ of size $(n-\nu )(p-k)$ corresponding to the indices $k+1,\dots ,k+(n-\nu )(p-k)$, where $\nu =\lceil (k+1)/(p-k)\rceil$ as in (39). Then, ${\tilde{A}}_{n,[p,k]}\left(1\right)=n{T}_{n-\nu }\left({\kappa }_{[p,k]}\right)$.

Proof. 

By (34) and (40), for all $r,R=1,\dots ,n-\nu$ and $q,Q=1,\dots ,p-k$ we have

$$\begin{array}{cc}\hfill {\left({\tilde{A}}_{n,[p,k]}\left(1\right)\right)}_{(p-k)(r-1)+q,(p-k)(R-1)+Q}& ={\int }_0^1{B}_{k+1+(p-k)(R-1)+Q,[p,k]}^{\prime }\left(x\right)B_{k+1+(p-k)(r-1)+q,[p,k]}^{\prime }\left(x\right)\mathrm{d}x\hfill \\ & ={\int }_{\mathbb{R}}{B}_{k+1+(p-k)(R-1)+Q,[p,k]}^{\prime }\left(x\right)B_{k+1+(p-k)(r-1)+q,[p,k]}^{\prime }\left(x\right)\mathrm{d}x\hfill \\ & =n^2{\int }_{\mathbb{R}}{\beta }_{Q,[p,k]}^{\prime }(nx-R+1){\beta }_{q,[p,k]}^{\prime }(nx-r+1)\mathrm{d}x\hfill \\ & =n{\int }_{\mathbb{R}}{\beta }_{Q,[p,k]}^{\prime }\left(y\right){\beta }_{q,[p,k]}^{\prime }(y-r+R)\mathrm{d}y\hfill \end{array}$$

and

$${\left(T_{n-\nu }\left({\kappa }_{[p,k]}\right)\right)}_{(p-k)(r-1)+q,(p-k)(R-1)+Q}={\left(K_{[p,k]}^{[r-R]}\right)}_{q,Q}={\int }_{\mathbb{R}}{\beta }_{Q,[p,k]}^{\prime }\left(y\right){\beta }_{q,[p,k]}^{\prime }(y-r+R)\mathrm{d}y,$$

which completes the proof. ☐

Theorem 4.

Let $a\in L^1\left([0,1]\right)$, $p\ge 1$ and $0\le k\le p-1$. Then, ${\left\{n^{-1}{A}_{n,[p,k]}\left(a\right)\right\}}_n{\sim }_{\sigma ,\lambda }a\left(x\right){\kappa }_{[p,k]}\left(\theta \right)$.

Proof. 

The proof consists of four steps. Throughout this proof, we use the following notation.

• $\nu =\lceil (k+1)/(p-k)\rceil$ as in (39).
• For every square matrix A of size $n(p-k)+k-1$, we denote by $\tilde{A}$ the principal submatrix of A corresponding to the row and column indices $i,j=k+1,\dots ,k+(n-\nu )(p-k)$.
• $P_{n,[p,k]}$ is the $\left(n\right(p-k)+k-1)\times (n-\nu )(p-k)$ matrix having $I_{(n-\nu )(p-k)}$ as the principal submatrix corresponding to the row and column indices $i,j=k+1,\dots ,k+(n-\nu )(p-k)$ and zeros elsewhere. Note that $P_{n,[p,k]}^T{P}_{n,[p,k]}=I_{(n-\nu )(p-k)}$ and $P_{n,[p,k]}^{T}A{P}_{n,[p,k]}=\tilde{A}$ for every square matrix A of size $n(p-k)+k-1$.

Step 1. Consider the linear operator $A_{n,[p,k]}(·):L^1\left([0,1]\right)\to {\mathbb{R}}^{\left(n\right(p-k)+k-1)\times \left(n\right(p-k)+k-1)}$,

$$A_{n,[p,k]}\left(g\right)={\left[{\int }_0^{1}g\left(x\right)B_{j+1,[p,k]}^{\prime }\left(x\right)B_{i+1,[p,k]}^{\prime }\left(x\right)\mathrm{d}x\right]}_{i,j=1}^{n(p-k)+k-1}.$$

The next three steps are devoted to show that

$${\left\{P_{n,[p,k]}^T\left(n^{-1}{A}_{n,[p,k]}\left(g\right)\right)P_{n,[p,k]}\right\}}_n={\left\{n^{-1}{\tilde{A}}_{n,[p,k]}\left(g\right)\right\}}_n{\sim }_{\mathrm{GLT}}g\left(x\right){\kappa }_{[p,k]}\left(\theta \right),\forall g\in L^1\left([0,1]\right).$$

(45)

Once this is done, the theorem is proven. Indeed, from (45), we immediately obtain the relation ${\left\{P_{n,[p,k]}^T\left(n^{-1}{A}_{n,[p,k]}\left(a\right)\right)P_{n,[p,k]}\right\}}_n{\sim }_{\mathrm{GLT}}a\left(x\right){\kappa }_{[p,k]}\left(\theta \right)$. We infer that ${\left\{P_{n,[p,k]}^T\left(n^{-1}{A}_{n,[p,k]}\left(a\right)\right)P_{n,[p,k]}\right\}}_n{\sim }_{\sigma ,\lambda }a\left(x\right){\kappa }_{[p,k]}\left(\theta \right)$ by GLT 1 and ${\left\{n^{-1}{A}_{n,[p,k]}\left(a\right)\right\}}_n{\sim }_{\sigma ,\lambda }a\left(x\right){\kappa }_{[p,k]}\left(\theta \right)$ by Theorem 1.

Step 2. We first prove (45) in the constant-coefficient case where $g\left(x\right)=1$ identically. In this case, by Lemma 2, we have $n^{-1}{\tilde{A}}_{n,[p,k]}\left(1\right)=T_{n-\nu }\left({\kappa }_{[p,k]}\right)$. Hence, the desired relation ${\{n^{-1}{\tilde{A}}_{n,[p,k]}\left(1\right)\}}_n{\sim }_{\mathrm{GLT}}{\kappa }_{[p,k]}\left(\theta \right)$ follows from GLT 2.

Step 3. Now we prove (45) in the case where $g\in C\left(\right[0,1\left]\right)$. Let

$$Z_{n,[p,k]}\left(g\right)=n^{-1}{\tilde{A}}_{n,[p,k]}\left(g\right)-n^{-1}{D}_{n-\nu }\left(g{I}_{p-k}\right){\tilde{A}}_{n,[p,k]}\left(1\right).$$

By (33), (34) and (38), for all $r,R=1,\dots ,n-\nu$ and $q,Q=1,\dots ,p-k$, we have

$$\begin{array}{c}|{\left(n{Z}_{n,[p,k]}\left(g\right)\right)}_{(p-k)(r-1)+q,(p-k)(R-1)+Q}|\hfill \\ =\left|{\left({\tilde{A}}_{n,[p,k]}\left(g\right)\right)}_{(p-k)(r-1)+q,(p-k)(R-1)+Q}-{\left(D_{n-\nu }\left(g{I}_{p-k}\right){\tilde{A}}_{n,[p,k]}\left(1\right)\right)}_{(p-k)(r-1)+q,(p-k)(R-1)+Q}\right|\hfill \\ =\left|{\int }_0^1\left[g\left(x\right)-g\left(\frac{r}{n-\nu }\right)\right]B_{k+1+(p-k)(R-1)+Q,[p,k]}^{\prime }\left(x\right)B_{k+1+(p-k)(r-1)+q,[p,k]}^{\prime }\left(x\right)\mathrm{d}x\right|\hfill \\ =\left|{\int }_{{\tau }_{k+1+(p-k)(r-1)+q}}^{{\tau }_{k+1+(p-k)(r-1)+q+p+1}}\left[g\left(x\right)-g\left(\frac{r}{n-\nu }\right)\right]B_{k+1+(p-k)(R-1)+Q,[p,k]}^{\prime }\left(x\right)B_{k+1+(p-k)(r-1)+q,[p,k]}^{\prime }\left(x\right)\mathrm{d}x\right|\hfill \\ \le c_p^2{n}^2{\int }_{(r-1)/n}^{(r+p)/n}\left|g\left(x\right)-g\left(\frac{r}{n-\nu }\right)\right|\mathrm{d}x\le c_p^2(p+1)n{\omega }_g\left(\frac{\nu +p}{n}\right),\hfill \end{array}$$

where ${\omega }_g(·)$ is the modulus of continuity of g and the last inequality is justified by the fact that the distance of the point $r/(n-\nu )$ from the interval $\left[\right(r-1)/n,(r+p)/n]$ is not larger than $(\nu +p)/n$. It follows that each entry of $Z_{n,[p,k]}\left(g\right)$ is bounded in modulus by $C_p{\omega }_g(1/n)$, where $C_p$ is a constant depending only on p. Moreover, by (34), the matrix $Z_{n,[p,k]}\left(g\right)$ is banded with bandwidth bounded by a constant $w_p$ depending only on p. Thus, by (30), $\parallel Z_{n,[p,k]}\left(g\right)\parallel \le w_p{C}_p{\omega }_g(1/n)\to 0$ as $n\to \infty$, and so ${\left\{Z_{n,[p,k]}\left(g\right)\right\}}_n$ is zero-distributed by Proposition 1. Since

$$n^{-1}{\tilde{A}}_{n,[p,k]}\left(g\right)=n^{-1}{D}_{n-\nu }\left(g{I}_{p-k}\right){\tilde{A}}_{n,[p,k]}\left(1\right)+Z_{n,[p,k]}\left(g\right),$$

we conclude that ${\{n^{-1}{\tilde{A}}_{n,[p,k]}\left(g\right)\}}_n{\sim }_{\mathrm{GLT}}g\left(x\right){\kappa }_{[p,k]}\left(\theta \right)$ by GLT 2, GLT 3 and Step 2.

Step 4. Finally, we prove (45) in the general case where $g\in L^1\left([0,1]\right)$. By the density of $C\left(\right[0,1\left]\right)$ in $L^1\left([0,1]\right)$, there exist functions $g_m\in C\left([0,1]\right)$ such that $g_m\to g$ in $L^1\left([0,1]\right)$. By Step 3,

$${\left\{n^{-1}{\tilde{A}}_{n,[p,k]}\left(g_m\right)\right\}}_n{\sim }_{\mathrm{GLT}}{g}_m\left(x\right){\kappa }_{[p,k]}\left(\theta \right).$$

(46)

Moreover,

$$g_m\left(x\right){\kappa }_{[p,k]}\left(\theta \right)\to g\left(x\right){\kappa }_{[p,k]}\left(\theta \right)\mathrm{in} \mathrm{measure}.$$

(47)

We show that

$${\left\{n^{-1}{\tilde{A}}_{n,[p,k]}\left(g_m\right)\right\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\left\{n^{-1}{\tilde{A}}_{n,[p,k]}\left(g\right)\right\}}_n.$$

(48)

Once this is done, the thesis (45) follows immediately from GLT 4. To prove (48), we recall that

$${\parallel X\parallel }_1\le \sum _{i,j=1}^N\left|x_{ij}\right|,X\in {\mathbb{C}}^{N\times N};$$

(49)

see, e.g., ([13] Section 2.4.3). By (38), we obtain

$$\begin{array}{cc}\hfill \parallel {\tilde{A}}_{n,[p,k]}\left(g\right)-{\tilde{A}}_{n,[p,k]}\left(g_m\right){\parallel }_1& \le \sum _{i,j=1}^{n(p-k)+k-1}\left|{\int }_0^1\left[g\left(x\right)-g_m\left(x\right)\right]B_{j+1,[p,k]}^{\prime }\left(x\right)B_{i+1,[p,k]}^{\prime }\left(x\right)\mathrm{d}x\right|\hfill \\ & \le {\int }_0^1\left|g\left(x\right)-g_m\left(x\right)\right|\sum _{i,j=1}^{n(p-k)+k-1}|B_{j+1,[p,k]}^{\prime }\left(x\right)\left|\right|B_{i+1,[p,k]}^{\prime }\left(x\right)|\mathrm{d}x\hfill \\ & \le c_p^2{n}^2{\parallel g-g_m\parallel }_{L^1}.\hfill \end{array}$$

Thus, the a.c.s. convergence (48) follows from Proposition 2. ☐

Remark 4.

By following step by step the proof of Theorem 4, we can give an alternative (much simpler) proof of ([36] Theorem A.6) based on the theory of block GLT sequences.

5. The Theory of Multilevel Block GLT Sequences

As illustrated in Section 3 and Section 4, the theory of block GLT sequences allows the computation of the singular value and eigenvalue distribution of block structured matrices arising from the discretization of univariate DEs. In order to cope with multivariate DEs, i.e., PDEs, we need the multivariate version of the theory of block GLT sequences, also known as the theory of multilevel block GLT sequences. The present section is devoted to a careful presentation of this theory, which is obtained by combining the results of [34] with the necessary technicalities for tackling multidimensional problems [14].

Multi-Index Notation. The multi-index notation is an essential tool for dealing with sequences of matrices arising from the discretization of PDEs. A multi-index $\mathit{i}\in {\mathbb{Z}}^d$, also called a d-index, is simply a (row) vector in ${\mathbb{Z}}^d$; its components are denoted by $i_1,\dots ,i_d$.

• $\mathbf{0},\mathbf{1},\mathbf{2},\dots$ are the vectors of all zeros, all ones, all twos, etc. (their size will be clear from the context).
• For any d-index $\mathit{m}$, we set $N\left(\mathit{m}\right)={\prod }_{j=1}^d{m}_j$ and we write $\mathit{m}\to \infty$ to indicate that $\min \left(\mathit{m}\right)\to \infty$.
• If $\mathit{h},\mathit{k}$ are d-indices, $\mathit{h}\le \mathit{k}$ means that $h_r\le k_r$ for all $r=1,\dots ,d$.
• If $\mathit{h},\mathit{k}$ are d-indices such that $\mathit{h}\le \mathit{k}$, the multi-index range $\mathit{h},\dots ,\mathit{k}$ is the set $\{\mathit{j}\in {\mathbb{Z}}^d:\mathit{h}\le \mathit{j}\le \mathit{k}\}$. We assume for this set the standard lexicographic ordering:

  $${\left[\dots {\left[{\left[(j_1,\dots ,j_d)\right]}_{j_d=h_d,\dots ,k_d}\right]}_{j_{d-1}=h_{d-1},\dots ,k_{d-1}}\dots \right]}_{j_1=h_1,\dots ,k_1}.$$

  (50)
  For instance, in the case $d=2$, the ordering is the following: $(h_1,h_2)$, $(h_1,h_2+1)$, $\dots ,$$(h_1,k_2)$, $(h_1+1,h_2)$, $(h_1+1,h_2+1)$, $\dots ,$$(h_1+1,k_2)$, $\dots \dots \dots ,$$(k_1,h_2)$, $(k_1,h_2+1)$, $\dots ,$$(k_1,k_2)$.
• When a d-index $\mathit{j}$ varies over a multi-index range $\mathit{h},\dots ,\mathit{k}$ (this is sometimes written as $\mathit{j}=\mathit{h},\dots ,\mathit{k}$), it is understood that $\mathit{j}$ varies from $\mathit{h}$ to $\mathit{k}$ following the specific ordering (50). For instance, if $\mathit{m}\in {\mathbb{N}}^d$ and if we write $\mathbf{x}={\left[x_{\mathit{i}}\right]}_{\mathit{i}=\mathbf{1}}^{\mathit{m}}$, then $\mathbf{x}$ is a vector of size $N\left(\mathit{m}\right)$ whose components $x_{\mathit{i}},\mathit{i}=\mathbf{1},\dots ,\mathit{m}$, are ordered in accordance with (50): the first component is $x_{\mathbf{1}}=x_{(1,\dots ,1,1)}$, the second component is $x_{(1,\dots ,1,2)}$, and so on until the last component, which is $x_{\mathit{m}}=x_{(m_1,\dots ,m_d)}$. Similarly, if $X={\left[x_{\mathit{i}\mathit{j}}\right]}_{\mathit{i},\mathit{j}=\mathbf{1}}^{\mathit{m}}$, then X is a $N\left(\mathit{m}\right)\times N\left(\mathit{m}\right)$ matrix whose components are indexed by two d-indices $\mathit{i},\mathit{j}$, both varying from $\mathbf{1}$ to $\mathit{m}$ according to the lexicographic ordering (50).
• Given $\mathit{h},\mathit{k}\in {\mathbb{Z}}^d$ with $\mathit{h}\le \mathit{k}$, the notation ${\sum }_{\mathit{j}=\mathit{h}}^{\mathit{k}}$ indicates the summation over all $\mathit{j}$ in $\mathit{h},\dots ,\mathit{k}$.
• Operations involving d-indices that have no meaning in the vector space ${\mathbb{Z}}^d$ must be interpreted in the componentwise sense. For instance, $\mathit{i}\mathit{j}=(i_1{j}_1,\dots ,i_d{j}_d),\mathit{i}/\mathit{j}=(i_1/j_1,\dots ,i_d/j_d)$, etc.

Multilevel Block Matrix-Sequences. Given $d,s\ge 1$, a d-level s-block matrix-sequence (or simply a matrix-sequence if d and s can be inferred from the context or we do not need/want to specify them) is a sequence of matrices of the form $\left\{A_{\mathit{n}}\right\}{}_n$, where:

• n varies in some infinite subset of $\mathbb{N}$;
• $\mathit{n}=\mathit{n}\left(n\right)$ is a d-index in ${\mathbb{N}}^d$ which depends on n and satisfies $\mathit{n}\to \infty$ as $n\to \infty$;
• A_n is a square matrix of size $N\left(\mathit{n}\right)s$.

Multilevel Block Toeplitz Matrices. Given a function $f:{[-\pi ,\pi ]}^d\to {\mathbb{C}}^{s\times s}$ in $L^1\left({[-\pi ,\pi ]}^d\right)$, its Fourier coefficients are denoted by

$$f_{\mathit{k}}=\frac{1}{{\left(2\pi \right)}^d}{\int }_{{[-\pi ,\pi ]}^d}f\left(\mathit{\theta }\right){\mathrm{e}}^{-\mathrm{i}\mathit{k}·\mathit{\theta }}\mathrm{d}\mathit{\theta }\in {\mathbb{C}}^{s\times s},\mathit{k}\in {\mathbb{Z}}^d,$$

where $\mathit{k}·\mathit{\theta }=k_1{\theta }_1+\dots +k_d{\theta }_d$ and the integrals are computed componentwise. For $\mathit{n}\in {\mathbb{N}}^d$, the nth multilevel block Toeplitz matrix generated by f is defined as

$$T_{\mathit{n}}\left(f\right)={\left[f_{\mathit{i}-\mathit{j}}\right]}_{\mathit{i},\mathit{j}=\mathbf{1}}^{\mathit{n}}\in {\mathbb{C}}^{N\left(\mathit{n}\right)s\times N\left(\mathit{n}\right)s}.$$

It is not difficult to see that the map $f\mapsto T_{\mathit{n}}\left(f\right)$ is linear. Moreover, it can be shown that

$$T_{\mathit{n}}(f^{*}{)=(T_{\mathit{n}}\left(f\right))}^{*},$$

(51)

where the transpose conjugate function $f^{*}$ is defined by $f^{*}\left(\mathit{\theta }\right)={\left(f\left(\mathit{\theta }\right)\right)}^{*}$; in particular, all the matrices $T_{\mathit{n}}\left(f\right)$ are Hermitian whenever f is Hermitian a.e. We also recall that, if $\mathit{n}\in {\mathbb{N}}^d$ and $f_1,f_2,\dots ,f_d:[-\pi ,\pi ]\to \mathbb{C}$ belong to $L^1\left([-\pi ,\pi ]\right)$, then

$$T_{n_1}\left(f_1\right)\otimes T_{n_2}\left(f_2\right)\otimes \cdots \otimes T_{n_d}\left(f_d\right)=T_{\mathit{n}}\left(f\right),$$

(52)

where $f:{[-\pi ,\pi ]}^d\to \mathbb{C}$ is defined by $f\left(\mathit{\theta }\right)=f\left({\theta }_1\right)f\left({\theta }_2\right)\cdots f\left({\theta }_d\right)$; see, e.g., ([14] Lemma 3.3).

Multilevel Block Diagonal Sampling Matrices. For $\mathit{n}\in {\mathbb{N}}^d$ and $a:{[0,1]}^d\to {\mathbb{C}}^{s\times s}$, we define the multilevel block diagonal sampling matrix $D_{\mathit{n}}\left(a\right)$ as the block diagonal matrix

$$D_{\mathit{n}}\left(a\right)=\underset{\mathit{i}=\mathbf{1},\dots ,\mathit{n}}{\mathrm{diag}}a\left(\frac{\mathit{i}}{\mathit{n}}\right)\in {\mathbb{C}}^{N\left(\mathit{n}\right)s\times N\left(\mathit{n}\right)s}.$$

Multilevel Block GLT Sequences. Let $d,s\ge 1$ be fixed positive integers. A d-level s-block GLT sequence (or simply a GLT sequence if d and s can be inferred from the context or we do not need/want to specify them) is a special d-level s-block matrix-sequence $\left\{A_{\mathit{n}}\right\}{}_n$ equipped with a measurable function $\kappa :{[0,1]}^d\times {[-\pi ,\pi ]}^d\to {\mathbb{C}}^{s\times s}$, the so-called symbol. We use the notation $\left\{A_{\mathit{n}}\right\}{}_n{\sim }_{\mathrm{GLT}}\kappa$ to indicate that $\left\{A_{\mathit{n}}\right\}{}_n$ is a GLT sequence with symbol $\kappa$. The symbol of a GLT sequence is unique in the sense that if $\left\{A_{\mathit{n}}\right\}{}_n{\sim }_{\mathrm{GLT}}\kappa$ and $\left\{A_{\mathit{n}}\right\}{}_n{\sim }_{\mathrm{GLT}}\varsigma$ then $\kappa =\varsigma$ a.e. in ${[0,1]}^d\times {[-\pi ,\pi ]}^d$. The main properties of d-level s-block GLT sequences are listed below.

GLT 1.

If $\left\{A_{\mathit{n}}\right\}{}_n{\sim }_{\mathrm{GLT}}\kappa$ then $\left\{A_{\mathit{n}}\right\}{}_n{\sim }_{\sigma }\kappa$. If, moreover, each A_n is Hermitian then $\left\{A_{\mathit{n}}\right\}{}_n{\sim }_{\lambda }\kappa$.

GLT 2.

We have:

• $\left\{T_{\mathit{n}}\right(f){\}}_n{\sim }_{\mathrm{GLT}}\kappa (\mathbf{x},\mathit{\theta })=f\left(\mathit{\theta }\right)$ if $f:{[-\pi ,\pi ]}^d\to {\mathbb{C}}^{s\times s}$ is in $L^1\left({[-\pi ,\pi ]}^d\right)$;
• $\left\{D_{\mathit{n}}\right(a){\}}_n{\sim }_{\mathrm{GLT}}\kappa (\mathbf{x},\mathit{\theta })=a\left(\mathbf{x}\right)$ if $a:{[0,1]}^d\to {\mathbb{C}}^{s\times s}$ is Riemann-integrable;
• $\left\{Z_{\mathit{n}}\right\}n{\sim }_{\mathrm{GLT}}\kappa (\mathbf{x},\mathit{\theta })=O_s$ if and only if $\left\{Z_{\mathit{n}}\right\}{}_n{\sim }_{\sigma }0$.

GLT 3.

If $\left\{A_{\mathit{n}}\right\}{}_n{\sim }_{\mathrm{GLT}}\kappa$ and $\left\{B_{\mathit{n}}\right\}{}_n{\sim }_{\mathrm{GLT}}\varsigma$ then:

• ${\left\{A_{\mathit{n}}^{*}\right\}}_n{\sim }_{\mathrm{GLT}}{\kappa }^{*}$;
• $\{\alpha A_{\mathit{n}}+\beta B_{\mathit{n}}\}{}_n{\sim }_{\mathrm{GLT}}\alpha \kappa +\beta \varsigma$ for all $\alpha ,\beta \in \mathbb{C}$;
• ${\left\{A_{\mathit{n}}{B}_{\mathit{n}}\right\}}_n{\sim }_{\mathrm{GLT}}\kappa \varsigma$;
• ${\{A_{\mathit{n}}^{†}\}}_n{\sim }_{\mathrm{GLT}}{\kappa }^{-1}$ provided that $\kappa$ is invertible a.e.

GLT 4.

$\left\{A_{\mathit{n}}\right\}{}_n{\sim }_{\mathrm{GLT}}\kappa$ if and only if there exist GLT sequences ${\left\{B_{\mathit{n},m}\right\}}_n{\sim }_{\mathrm{GLT}}{\kappa }_m$ such that ${\left\{B_{\mathit{n},m}\right\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}\{A_{\mathit{n}}\}{}_n$ and ${\kappa }_m\to \kappa$ in measure.

6. Discretizations of Systems of PDEs: The General GLT Approach

In this section, we outline the main ideas of a multidimensional block GLT analysis for general discretizations of PDE systems. What we are going to present here is then a generalization of what is shown in Section 3. We begin by proving a series of auxiliary results. In the following, given $\mathit{n}\in {\mathbb{N}}^d$ and $s\ge 1$, we denote by ${\mathsf{\Pi }}_{\mathit{n},s}$ the permutation matrix given by

$${\mathsf{\Pi }}_{\mathit{n},s}=\left[\begin{array}{c}{I}_s\otimes {\mathbf{e}}_{\mathbf{1}}^T\\ I_s\otimes {\mathbf{e}}_{\mathbf{2}}^T\\ ⋮\\ I_s\otimes {\mathbf{e}}_{\mathit{n}}^T\end{array}\right]=\sum _{\mathbf{k}=\mathbf{1}}^{\mathit{n}}{\mathbf{e}}_k\otimes I_s\otimes {\mathbf{e}}_{\mathit{k}}^T,$$

(53)

where ${\mathbf{e}}_{\mathit{i}},\mathit{i}=\mathbf{1},\dots ,\mathit{n},$ are the vectors of the canonical basis of ${\mathbb{R}}^{N\left(\mathit{n}\right)}$, which, for convenience, are indexed by a d-index $\mathit{i}=\mathbf{1},\dots ,\mathit{n}$ instead of a linear index $i=1,\dots ,N\left(\mathit{n}\right)$. Note that ${\mathsf{\Pi }}_{n,2}$ coincides with the matrix ${\mathsf{\Pi }}_n$ in (21).

Lemma 3.

Let $\mathit{n}\in {\mathbb{N}}^d$, let $f_{ij}:{[-\pi ,\pi ]}^d\to \mathbb{C}$ be in $L^1\left({[-\pi ,\pi ]}^d\right)$ for $i,j=1,\dots ,s$, and set $f={\left[f_{ij}\right]}_{i,j=1}^s$. The block matrix $T_{\mathit{n}}=\left[T_{\mathit{n}}\right(f_{ij}){]}_{i,j=1}^s$ is similar via the permutation (53) to the multilevel block Toeplitz matrix $T_{\mathit{n}}(f)$, that is, ${\mathsf{\Pi }}_{\mathit{n},s}{T}_{\mathit{n}}{\mathsf{\Pi }}_{\mathit{n},s}^T=T_{\mathit{n}}(f)$.

Proof. 

Let $E_{ij}$ be the $s\times s$ matrix having 1 in position $(i,j)$ and 0 elsewhere. Since $T_{\mathit{n}}={\sum }_{i,j=1}^s{E}_{ij}\otimes T_{\mathit{n}}(f_{ij})$ and $T_{\mathit{n}}(f)={\sum }_{i,j=1}^s{T}_{\mathit{n}}(f_{ij}{E}_{ij})$ by the linearity of the map $T_{\mathit{n}}(·)$, it is enough to show that

$${\mathsf{\Pi }}_{\mathit{n},s}(E\otimes T_{\mathit{n}}(g)){\mathsf{\Pi }}_{\mathit{n},s}^T=T_{\mathit{n}}(gE),\forall g\in L^1\left({[-\pi ,\pi ]}^d\right),\forall E\in {\mathbb{C}}^{s\times s}.$$

By (6) and (7),

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}_{\mathit{n},s}(E\otimes T_{\mathit{n}}(g)){\mathsf{\Pi }}_{\mathit{n},s}^T& =\left[\sum _{\mathit{k}=\mathbf{1}}^{\mathit{n}}{\mathbf{e}}_{\mathit{k}}\otimes I_s\otimes {\mathbf{e}}_{\mathit{k}}^T\right](E\otimes T_{\mathit{n}}(g))\left[\sum _{\ell =\mathbf{1}}^{\mathit{n}}{\mathbf{e}}_{\ell }^T\otimes I_s\otimes {\mathbf{e}}_{\ell }\right]\hfill \\ & =\sum _{\mathit{k},\ell =1}^{\mathit{n}}({\mathbf{e}}_{\mathit{k}}\otimes I_s\otimes {\mathbf{e}}_{\mathit{k}}^T)(E\otimes T_{\mathit{n}}(g))({\mathbf{e}}_{\ell }^T\otimes I_s\otimes {\mathbf{e}}_{\ell })\hfill \\ & =\sum _{\mathit{k},\ell =\mathbf{1}}^{\mathit{n}}{\mathbf{e}}_{\mathit{k}}{\mathbf{e}}_{\ell }^T\otimes E\otimes {\mathbf{e}}_{\mathit{k}}^T{T}_{\mathit{n}}(g){\mathbf{e}}_{\ell }=\sum _{\mathit{k},\ell =\mathbf{1}}^{\mathit{n}}{\mathbf{e}}_{\mathit{k}}{\mathbf{e}}_{\ell }^T\otimes \left(T_{\mathit{n}}\right(g){)}_{\mathit{k}\ell }E=T_{\mathit{n}}\left(gE\right),\hfill \end{array}$$

as required. ☐

Lemma 4.

Let $\mathit{n}\in {\mathbb{N}}^d$, let $a_{ij}:{[0,1]}^d\to \mathbb{C}$ for $i,j=1,\dots ,s$, and set $a={\left[a_{ij}\right]}_{i,j=1}^s$. The block matrix $D_{\mathit{n}}=\left[D_{\mathit{n}}\right(a_{ij}){]}_{i,j=1}^s$ is similar via the permutation (53) to the multilevel block diagonal sampling matrix $D_{\mathit{n}}(a)$, that is, ${\mathsf{\Pi }}_{\mathit{n},s}{D}_{\mathit{n}}{\mathsf{\Pi }}_{\mathit{n},s}^T=D_{\mathit{n}}(a)$.

Proof. 

With obvious adaptations, it is the same as the proof of Lemma 3. ☐

We recall that a d-variate trigonometric polynomial is a finite linear combination of the d-variate Fourier frequencies ${\mathrm{e}}^{\mathrm{i}\mathit{k}·\theta },\mathit{k}\in {\mathbb{Z}}^d$.

Theorem 5.

For $i,j=1,\dots ,s$, let ${\left\{A_{\mathit{n},ij}\right\}}_n$ be a d-level 1-block GLT sequence with symbol ${\kappa }_{ij}:{[0,1]}^d\times {[-\pi ,\pi ]}^d\to \mathbb{C}$. Set $A_{\mathit{n}}={\left[A_{\mathit{n},ij}\right]}_{i,j=1}^s$ and $\kappa ={\left[{\kappa }_{ij}\right]}_{i,j=1}^s$. Then, the matrix-sequence ${\left\{{\mathsf{\Pi }}_{\mathit{n},s}{A}_{\mathit{n}}{\mathsf{\Pi }}_{\mathit{n},s}^T\right\}}_n$ is a d-level s-block GLT sequence with symbol κ.

Proof. 

The proof consists of the following two steps.

Step 1. We first prove the theorem under the additional assumption that $A_{\mathit{n},ij}$ is of the form

$$A_{\mathit{n},ij}=\sum _{\ell =1}^{L_{ij}}{D}_{\mathit{n}}(a_{\ell ,ij})T_{\mathit{n}}(f_{\ell ,ij}),$$

(54)

where $L_{ij}\in \mathbb{N},a_{\ell ,ij}:{[0,1]}^d\to \mathbb{C}$ is Riemann-integrable, and $f_{\ell ,ij}:{[-\pi ,\pi ]}^d\to \mathbb{C}$ belongs to $L^1\left({[-\pi ,\pi ]}^d\right)$. Note that the symbol of ${\left\{A_{\mathit{n},ij}\right\}}_n$ is

$${\kappa }_{ij}(\mathbf{x},\mathit{\theta })=\sum _{\ell =1}^{L_{ij}}{a}_{\ell ,ij}\left(\mathbf{x}\right)f_{\ell ,ij}\left(\mathit{\theta }\right).$$

By setting $L={\max }_{i,j=1,\dots ,s}{L}_{ij}$ and by adding zero matrices of the form $D_{\mathit{n}}(0\left)T_{\mathit{n}}\right(0)$ in the summation (54) whenever $L_{ij}<L$, we can assume, without loss of generality, that

$$\begin{array}{cc}\hfill A_{\mathit{n},ij}& =\sum _{\ell =1}^L{D}_{\mathit{n}}(a_{\ell ,ij})T_{\mathit{n}}(f_{\ell ,ij}),\hfill \\ \hfill {\kappa }_{ij}(\mathbf{x},\mathit{\theta })& =\sum _{\ell =1}^L{a}_{\ell ,ij}\left(\mathbf{x}\right)f_{\ell ,ij}\left(\mathit{\theta }\right),\hfill \end{array}$$

with L independent of $i,j$. Let $E_{ij}$ be the $s\times s$ matrix having 1 in position $(i,j)$ and 0 elsewhere. Then,

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}_{\mathit{n},s}{A}_{\mathit{n}}{\mathsf{\Pi }}_{\mathit{n},s}^T& =\sum _{\ell =1}^L{\mathsf{\Pi }}_{\mathit{n},s}{\left[D_{\mathit{n}}(a_{\ell ,ij}\left)T_{\mathit{n}}\right(f_{\ell ,ij})\right]}_{i,j=1}^s{\mathsf{\Pi }}_{\mathit{n},s}^T\hfill \\ & =\sum _{\ell =1}^L{\mathsf{\Pi }}_{\mathit{n},s}\left[\sum _{i,j=1}^s(E_{ij}\otimes D_{\mathit{n}}(a_{\ell ,ij})\left)(E_{ij}\otimes T_{\mathit{n}}(f_{\ell ,ij})\right)\right]{\mathsf{\Pi }}_{\mathit{n},s}^T\hfill \\ & =\sum _{\ell =1}^L\sum _{i,j=1}^s{\mathsf{\Pi }}_{\mathit{n},s}(E_{ij}\otimes D_{\mathit{n}}(a_{\ell ,ij})\left){\mathsf{\Pi }}_{\mathit{n},s}^T{\mathsf{\Pi }}_{\mathit{n},s}(E_{ij}\otimes T_{\mathit{n}}(f_{\ell ,ij})\right){\mathsf{\Pi }}_{\mathit{n},s}^T.\hfill \end{array}$$

By Lemmas 3 and 4,

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}_{\mathit{n},s}(E_{ij}\otimes D_{\mathit{n}}(a_{\ell ,ij})){\mathsf{\Pi }}_{\mathit{n},s}^T& =D_{\mathit{n}}(a_{\ell ,ij}{E}_{ij}),\hfill \\ \hfill {\mathsf{\Pi }}_{\mathit{n},s}(E_{ij}\otimes T_{\mathit{n}}(f_{\ell ,ij})){\mathsf{\Pi }}_{\mathit{n},s}^T& =T_{\mathit{n}}(f_{\ell ,ij}{E}_{ij}).\hfill \end{array}$$

It follows that

$${\mathsf{\Pi }}_{\mathit{n},s}{A}_{\mathit{n}}{\mathsf{\Pi }}_{\mathit{n},s}^T=\sum _{\ell =1}^L\sum _{i,j=1}^s{D}_{\mathit{n}}(a_{\ell ,ij}{E}_{ij}\left)T_{\mathit{n}}\right(f_{\ell ,ij}{E}_{ij}).$$

Thus, by GLT 2 and GLT 3, ${\left\{{\mathsf{\Pi }}_{\mathit{n},s}{A}_{\mathit{n}}{\mathsf{\Pi }}_{\mathit{n},s}^T\right\}}_n$ is a d-level s-block GLT sequence with symbol

$$\kappa (\mathbf{x},\theta )=\sum _{\ell =1}^L\sum _{i,j=1}^s{a}_{\ell ,ij}\left(\mathbf{x}\right)f_{\ell ,ij}\left(\mathit{\theta }\right)E_{ij}={\left[{\kappa }_{ij}(\mathbf{x},\mathit{\theta })\right]}_{i,j=1}^s.$$

Step 2. We now prove the theorem in its full generality. Since ${\left\{A_{\mathit{n},ij}\right\}}_n{\sim }_{\mathrm{GLT}}{\kappa }_{ij}$, by ([14] Theorem 5.6) there exist functions $a_{\ell ,ij}^{\left(m\right)},f_{\ell ,ij}^{\left(m\right)},\ell =1,\dots ,L_{ij}^{\left(m\right)}$, such that

• $a_{\ell ,ij}^{\left(m\right)}\in C^{\infty }\left({[0,1]}^d\right)$ and $f_{\ell ,ij}^{\left(m\right)}$ is a d-variate trigonometric polynomial,
• ${\kappa }_{ij}^{\left(m\right)}(\mathbf{x},\mathit{\theta })={\sum }_{\ell =1}^{L_{ij}^{\left(m\right)}}{a}_{\ell ,ij}^{\left(m\right)}\left(\mathbf{x}\right)f_{\ell ,ij}^{\left(m\right)}\left(\mathit{\theta }\right)\to {\kappa }_{ij}(\mathbf{x},\mathit{\theta })$ a.e.;
• ${\left\{A_{\mathit{n},ij}^{\left(m\right)}={\sum }_{\ell =1}^{L_{ij}^{\left(m\right)}}{D}_{\mathit{n}}(a_{\ell ,ij}^{\left(m\right)})T_{\mathit{n}}(f_{\ell ,ij}^{\left(m\right)})\right\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\left\{A_{\mathit{n},ij}\right\}}_n$.

Set $A_{\mathit{n}}^{\left(m\right)}={\left[A_{\mathit{n},ij}^{\left(m\right)}\right]}_{i,j=1}^s$ and ${\kappa }^{\left(m\right)}(\mathbf{x},\mathit{\theta })={\left[{\kappa }_{ij}^{\left(m\right)}(\mathbf{x},\mathit{\theta })\right]}_{i,j=1}^s$. We have:

• ${\left\{{\mathsf{\Pi }}_{\mathit{n},s}{A}_{\mathit{n}}^{\left(m\right)}{\mathsf{\Pi }}_{\mathit{n},s}^T\right\}}_n{\sim }_{\mathrm{GLT}}{\kappa }^{\left(m\right)}$ by Step 1;
• ${\kappa }^{\left(m\right)}\to \kappa$ a.e. (and hence also in measure);
• ${\left\{{\mathsf{\Pi }}_{\mathit{n},s}{A}_{\mathit{n}}^{\left(m\right)}{\mathsf{\Pi }}_{\mathit{n},s}^T\right\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}{\left\{{\mathsf{\Pi }}_{\mathit{n},s}{A}_{\mathit{n}}{\mathsf{\Pi }}_{\mathit{n},s}^T\right\}}_n$ because ${\left\{A_{\mathit{n}}^{\left(m\right)}\right\}}_n\stackrel{\mathrm{a}.\mathrm{c}.\mathrm{s}.}{⟶}\{A_{\mathit{n}}\}{}_n$ by Lemma 1.

We conclude that ${\left\{{\mathsf{\Pi }}_{\mathit{n},s}{A}_{\mathit{n}}{\mathsf{\Pi }}_{\mathit{n},s}^T\right\}}_n{\sim }_{\mathrm{GLT}}\kappa$ by GLT 4. ☐

Now, suppose we have a system of linear PDEs of the form

$$\left\{\begin{array}{c}\sum _{j=1}^s{\mathcal{L}}_{1j}{u}_j\left(\mathbf{x}\right)=f_1\left(\mathbf{x}\right),\\ \sum _{j=1}^s{\mathcal{L}}_{2j}{u}_j\left(\mathbf{x}\right)=f_2\left(\mathbf{x}\right),\\ ⋮\\ \sum _{j=1}^s{\mathcal{L}}_{sj}{u}_j\left(\mathbf{x}\right)=f_s\left(\mathbf{x}\right),\end{array}\right.$$

(55)

where $\mathbf{x}\in {(0,1)}^d$. The matrices A_n resulting from any standard discretization of (55) are parameterized by a d-index $\mathit{n}=(n_1,\dots ,n_d)$, where $n_i$ is related to the discretization step $h_i$ in the ith direction, and $n_i\to \infty$ if and only if $h_i\to 0$ (usually, $h_i\approx 1/n_i$). By choosing each $n_i$ as a function of a unique discretization parameter $n\in \mathbb{N}$, as it normally happens in practice where the most natural choice is $n_i=n$ for all $i=1,\dots ,d$, we see that $\mathit{n}=\mathit{n}\left(n\right)$ and, consequently, $\left\{A_{\mathit{n}}\right\}{}_n$ is a (d-level) matrix-sequence. Moreover, it turns out that, after a suitable normalization that we ignore in this discussion—the normalization we are talking about is the analog of the normalization that we have seen in Section 3, which allowed us to pass from the matrix $A_n$ in (13) to the matrix $B_n$ in (15)—, A_n has the following block structure:

$$A_{\mathit{n}}={\left[A_{\mathit{n},ij}\right]}_{i,j=1}^s,$$

where $A_{\mathit{n},ij}$ is the (normalized) matrix arising from the discretization of the differential operator ${\mathcal{L}}_{ij}$. In the simplest case where the operators ${\mathcal{L}}_{ij}$ have constant coefficients and we use equispaced grids in each direction, the matrix $A_{\mathit{n},ij}$ takes the form

$$A_{\mathit{n},ij}=T_{\mathit{n}}(f_{ij})+Z_{\mathit{n},ij},$$

where $f_{ij}$ is a d-variate trigonometric polynomial, while the perturbation $Z_{\mathit{n},ij}$ is usually a low-rank correction due to boundary conditions and, in any case, we have ${\left\{Z_{\mathit{n},ij}\right\}}_n{\sim }_{\sigma }0$. Hence,

$${\left\{A_{\mathit{n},ij}\right\}}_n{\sim }_{\mathrm{GLT}}{f}_{ij}$$

by GLT 2 and GLT 3, and it follows from Theorem 5 that

$${\left\{{\mathsf{\Pi }}_{\mathit{n},s}{A}_{\mathit{n}}{\mathsf{\Pi }}_{\mathit{n},s}^T\right\}}_n{\sim }_{\mathrm{GLT}}{\left[f_{ij}\right]}_{i,j=1}^s.$$

In the case where the operators ${\mathcal{L}}_{ij}$ have variable coefficients, the matrix $A_{\mathit{n},ij}$ usually takes the form

$$A_{\mathit{n},ij}=\sum _{\ell =1}^{L_{ij}}{D}_{\mathit{n}}(a_{\ell ,ij})T_{\mathit{n}}(f_{\ell ,ij})+Z_{\mathit{n},ij},$$

where $L_{ij}\in \mathbb{N},f_{\ell ,ij}$ is a d-variate trigonometric polynomial, ${\left\{Z_{\mathit{n},ij}\right\}}_n{\sim }_{\sigma }0$, and the functions $a_{\ell ,ij}:{[0,1]}^d\to \mathbb{R},\ell =1,\dots ,L_{ij},$ are related to the coefficients of ${\mathcal{L}}_{ij}$ (for example, in Section 3, while proving (22), we have seen that $K_n\left(a_{11}\right)$, which plays there the same role as the matrix $A_{\mathit{n},11}$ here, is equal to $D_n\left(a_{11}\right)T_n(2-2\cos \theta )+Z_n$ for some zero-distributed sequence ${\left\{Z_n\right\}}_n$). Hence,

$${\left\{A_{\mathit{n},ij}\right\}}_n{\sim }_{\mathrm{GLT}}{\kappa }_{ij}(\mathbf{x},\mathit{\theta })=\sum _{\ell =1}^{L_{ij}}{a}_{\ell ,ij}\left(\mathbf{x}\right)f_{\ell ,ij}\left(\mathit{\theta }\right)$$

by GLT 2 and GLT 3, and it follows from Theorem 5 that

$${\left\{{\mathsf{\Pi }}_{\mathit{n},s}{A}_{\mathit{n}}{\mathsf{\Pi }}_{\mathit{n},s}^T\right\}}_n{\sim }_{\mathrm{GLT}}{\left[{\kappa }_{ij}\right]}_{i,j=1}^s.$$

7. B-Spline IgA Discretization of a Variational Problem for the Curl–Curl Operator

For any function $\mathbf{u}(x_1,x_2)={[u_1(x_1,x_2),u_2(x_1,x_2)]}^T$, defined over some open set $\mathsf{\Omega }\subseteq {\mathbb{R}}^2$ and taking values in ${\mathbb{R}}^2$, the curl operator is formally defined as follows:

$$(\nabla \times \mathbf{u})(x_1,x_2)=\frac{\partial u_2}{\partial x_1}(x_1,x_2)-\frac{\partial u_1}{\partial x_2}(x_1,x_2),(x_1,x_2)\in \mathsf{\Omega }.$$

Clearly, this definition has meaning when the components $u_1,u_2$ belong to $H^1\left(\mathsf{\Omega }\right)$, so that their partial derivatives exist in the Sobolev sense. Now, let $\mathsf{\Omega }={(0,1)}^2$, set

$$\begin{array}{cc}\hfill {\left(L^2\left(\mathsf{\Omega }\right)\right)}^2& =\{\mathbf{u}:\mathsf{\Omega }\to {\mathbb{R}}^2:u_1,u_2\in L^2\left(\mathsf{\Omega }\right)\},\hfill \\ \hfill H(\mathrm{curl},\mathsf{\Omega })& =\{\mathbf{u}\in {\left(L^2\left(\mathsf{\Omega }\right)\right)}^2:\nabla \times \mathbf{u}\mathrm{exists} \mathrm{in} \mathrm{the} \mathrm{Sobolev} \mathrm{sense},\nabla \times \mathbf{u}\in L^2\left(\mathsf{\Omega }\right)\},\hfill \end{array}$$

and consider the following variational problem: find $\mathbf{u}\in H(\mathrm{curl},\mathsf{\Omega })$ such that

$$(\nabla \times \mathbf{u},\nabla \times \mathbf{v})=(\mathbf{f},\mathbf{v}),\forall \mathbf{v}\in H(\mathrm{curl},\mathsf{\Omega }),$$

(56)

where $\mathbf{f}(x_1,x_2)={[f_1(x_1,x_2),f_2(x_1,x_2)]}^T$ is a vector field in ${\left(L^2\left(\mathsf{\Omega }\right)\right)}^2$ and

$$\begin{array}{cc}\hfill (\nabla \times \mathbf{u},\nabla \times \mathbf{v})& ={\int }_{\mathsf{\Omega }}(\nabla \times \mathbf{u})(x_1,x_2)(\nabla \times \mathbf{v})(x_1,x_2)\mathrm{d}{x}_1\mathrm{d}{x}_2,\hfill \\ \hfill (\mathbf{f},\mathbf{v})& ={\int }_{\mathsf{\Omega }}\left[f_1(x_1,x_2)v_1(x_1,x_2)+f_2(x_1,x_2)v_2(x_1,x_2)\right]\mathrm{d}{x}_1\mathrm{d}{x}_2.\hfill \end{array}$$

Variational problems of the form of (56) arise in important applications, such as time harmonic Maxwell’s equations and magnetostatic problems. In this section, we consider a so-called compatible B-spline IgA discretization of (56); see [43] for details. We show that the corresponding sequence of discretization matrices enjoys a spectral distribution described by a $2\times 2$ matrix-valued function whose determinant is zero everywhere. The results of this section have already been obtained in [38], but the derivation presented here is entirely based on the theory of multilevel block GLT sequences and turns out to be simpler and more lucid than that in [38]. For simplicity, throughout this section, the B-splines $B_{i,[p,p-1]}$, $i=1,\dots ,n+p$, and the associated reference B-spline ${\beta }_{1,[p,p-1]}$, are denoted by $B_{i,\left[p\right]}$, $i=1,\dots ,n+p$, and ${\beta }_{\left[p\right]}$, respectively. The function ${\beta }_{\left[p\right]}$ is the so-called cardinal B-spline of degree p over the knot sequence $\{0,1,\dots ,p+1\}$. In view of the following, we recall from [42] and ([23] Lemma 4) that the cardinal B-spline ${\beta }_{\left[q\right]}$ is defined for all degrees $q\ge 0$, belongs to $C^{q-1}\left(\mathbb{R}\right)$, and satisfies the following properties:

$$\mathrm{supp}\left({\beta }_{\left[q\right]}\right)=[0,q+1]$$

(57)

for $q\ge 1$,

$${\beta }_{\left[q\right]}^{\prime }\left(t\right)={\beta }_{[q-1]}\left(t\right)-{\beta }_{[q-1]}(t-1),$$

(58)

for $t\in \mathbb{R}$ and $q\ge 1$, and

$${\int }_{\mathbb{R}}{\beta }_{\left[q_1\right]}^{\left(r_1\right)}\left(\tau \right){\beta }_{\left[q_2\right]}^{\left(r_2\right)}(\tau +t)\mathrm{d}\tau ={(-1)}^{r_1}{\beta }_{[q_1+q_2+1]}^{(r_1+r_2)}(q_1+1+t)={(-1)}^{r_2}{\beta }_{[q_1+q_2+1]}^{(r_1+r_2)}(q_2+1-t)$$

(59)

for $t\in \mathbb{R}$ and $q_1,q_2,r_1,r_2\ge 0$. Moreover, property (40) in the case $k=p-1$ simplifies to

$$B_{i,\left[p\right]}\left(x\right)={\beta }_{\left[p\right]}(nx-i+p+1),i=p+1,\dots ,n.$$

(60)

7.1. Compatible B-Spline IgA Discretization

Let $\mathit{n}=(n_1,n_2)\in {\mathbb{N}}^2$, let $p\ge 2$, and define the space

$$\begin{array}{cc}\hfill {\mathcal{V}}_{\mathit{n},\left[p\right]}(\mathrm{curl},\mathsf{\Omega })=\mathrm{span}\left\{\left[\begin{array}{c}{B}_{i_1,[p-1]}\left(x_1\right)B_{i_2,\left[p\right]}\left(x_2\right)\\ B_{j_1,\left[p\right]}\left(x_1\right)B_{j_2,[p-1]}\left(x_2\right)\end{array}\right]:\right.i_1& =1,\dots ,n_1+p-1,i_2=1,\dots ,n_2+p,\hfill \\ \hfill j_1& =1,\dots ,n_1+p,j_2=1,\dots ,n_2+p-1\left.\phantom{\left[\begin{array}{c}{B}_{i_1,[p-1]}\left(x_1\right)B_{i_2,\left[p\right]}\left(x_2\right)\\ B_{j_1,\left[p\right]}\left(x_1\right)B_{j_2,[p-1]}\left(x_2\right)\end{array}\right]}\right\}.\hfill \end{array}$$

(61)

Following a compatible B-spline approach [43], we look for an approximation of the solution in the space ${\mathcal{V}}_{\mathit{n},\left[p\right]}(\mathrm{curl},\mathsf{\Omega })$ by solving the following discrete problem: find ${\mathbf{u}}_{\mathcal{V}}\in {\mathcal{V}}_{\mathit{n},\left[p\right]}(\mathrm{curl},\mathsf{\Omega })$ such that

$$(\nabla \times {\mathbf{u}}_{\mathcal{V}},\nabla \times \mathbf{v})=(\mathbf{f},\mathbf{v}),\forall \mathbf{v}\in {\mathcal{V}}_{\mathit{n},\left[p\right]}(\mathrm{curl},\mathsf{\Omega }).$$

After choosing a suitable ordering on the basis functions of ${\mathcal{V}}_{\mathit{n},\left[p\right]}(\mathrm{curl},\mathsf{\Omega })$ displayed in (61), by linearity the computation of ${\mathbf{u}}_{\mathcal{V}}$ reduces to solving a linear system whose coefficient matrix is given by

$$A_{\mathit{n},\left[p\right]}=\left[\begin{array}{cc}{A}_{\mathit{n},\left[p\right],11}& A_{\mathit{n},\left[p\right],12}\\ A_{\mathit{n},\left[p\right],21}& A_{\mathit{n},\left[p\right],22}\end{array}\right]=\left[\begin{array}{cc}{M}_{n_1,[p-1]}\otimes K_{n_2,\left[p\right]}& -H_{n_1,\left[p\right]}\otimes {\left(H_{n_2,\left[p\right]}\right)}^T\\ -{\left(H_{n_1,\left[p\right]}\right)}^T\otimes H_{n_2,\left[p\right]}& K_{n_1,\left[p\right]}\otimes M_{n_2,[p-1]}\end{array}\right],$$

where

$$\begin{array}{cc}\hfill {\left(M_{n,[p-1]}\right)}_{ij}& ={\int }_0^1{B}_{j,[p-1]}\left(x\right)B_{i,[p-1]}\left(x\right)\mathrm{d}x,i,j=1,\dots ,n+p-1,\hfill \\ \hfill {\left(K_{n,\left[p\right]}\right)}_{ij}& ={\int }_0^1{B}_{j,\left[p\right]}^{\prime }\left(x\right)B_{i,\left[p\right]}^{\prime }\left(x\right)\mathrm{d}x,i,j=1,\dots ,n+p,\hfill \\ \hfill {\left(H_{n,\left[p\right]}\right)}_{ij}& ={\int }_0^1{B}_{j,\left[p\right]}^{\prime }\left(x\right)B_{i,[p-1]}\left(x\right)\mathrm{d}x,i=1,\dots ,n+p-1,j=1,\dots ,n+p.\hfill \end{array}$$

Note that $M_{n,[p-1]}$ is a square matrix of size $n+p-1,K_{n,\left[p\right]}$ is a square matrix of size $n+p$, while $H_{n,\left[p\right]}$ is a rectangular matrix of size $(n+p-1)\times (n+p)$.

7.2. GLT Analysis of the B-Spline IgA Discretization Matrices

In the main result of this section (Theorem 6), assuming that $\mathit{n}=n\nu$ for a fixed vector $\nu$, we show that the spectral distribution of the sequence ${\left\{A_{\mathit{n},\left[p\right]}\right\}}_n$ is described by a $2\times 2$ matrix-valued function whose determinant is zero everywhere (Remark 5). To prove Theorem 6, some preliminary work is necessary. We first note that, in view of the application of Theorem 5, the matrix $A_{\mathit{n},\left[p\right]}$ has an unpleasant feature: the anti-diagonal blocks $A_{\mathit{n},\left[p\right],12}$ and $A_{\mathit{n},\left[p\right],21}$ are not square and the square diagonal blocks $A_{\mathit{n},\left[p\right],11}$ and $A_{\mathit{n},\left[p\right],22}$ do not have the same size whenever $n_1\ne n_2$. Let us then introduce the nicer matrix

$${\tilde{A}}_{\mathit{n},\left[p\right]}=\left[\begin{array}{cc}{\tilde{A}}_{\mathit{n},\left[p\right],11}& {\tilde{A}}_{\mathit{n},\left[p\right],12}\\ {\tilde{A}}_{\mathit{n},\left[p\right],21}& {\tilde{A}}_{\mathit{n},\left[p\right],22}\end{array}\right]=\left[\begin{array}{cc}{\tilde{M}}_{n_1,[p-1]}\otimes K_{n_2,\left[p\right]}& -{\tilde{H}}_{n_1,\left[p\right]}\otimes {({\tilde{H}}_{n_2,\left[p\right]})}^T\\ -{({\tilde{H}}_{n_1,\left[p\right]})}^T\otimes {\tilde{H}}_{n_2,\left[p\right]}& K_{n_1,\left[p\right]}\otimes {\tilde{M}}_{n_2,[p-1]}\end{array}\right],$$

where ${\tilde{M}}_{n,[p-1]}$ and ${\tilde{H}}_{n,\left[p\right]}$ are square matrices of size $n+p$ given by

$${\tilde{M}}_{n,[p-1]}=\left[\begin{array}{cccc}& & & 0\\ & M_{n,[p-1]}& & ⋮\\ & & & 0\\ 0& \cdots & 0& 0\end{array}\right],{\tilde{H}}_{n,\left[p\right]}=\left[\begin{array}{}\\ & H_{n,\left[p\right]}& \\ \\ 0& \cdots & 0\end{array}\right].$$

Each block ${\tilde{A}}_{\mathit{n},\left[p\right],ij}$ of the matrix ${\tilde{A}}_{\mathit{n},\left[p\right]}$ is now a square block of size $(n_1+p)(n_2+p)$. Moreover,

$$M_{n,[p-1]}=P_{n,\left[p\right]}{\tilde{M}}_{n,[p-1]}{\left(P_{n,\left[p\right]}\right)}^T,H_{n,\left[p\right]}=P_{n,\left[p\right]}{\tilde{H}}_{n,\left[p\right]},$$

where the matrix

$$P_{n,\left[p\right]}=\left[\begin{array}{cccc}& & & 0\\ & I_{n+p-1}& & ⋮\\ & & & 0\end{array}\right]$$

satisfies $P_{n,\left[p\right]}{\left(P_{n,\left[p\right]}\right)}^T=I_{n+p-1}$. By (6) and (7),

$$\begin{array}{cc}\hfill A_{\mathit{n},\left[p\right],11}& =(P_{n_1,\left[p\right]}\otimes I_{n_2+p}){\tilde{A}}_{\mathit{n},\left[p\right],11}{(P_{n_1,\left[p\right]}\otimes I_{n_2+p})}^T,\hfill \\ \hfill A_{\mathit{n},\left[p\right],12}& =(P_{n_1,\left[p\right]}\otimes I_{n_2+p}){\tilde{A}}_{\mathit{n},\left[p\right],12}{(I_{n_1+p}\otimes P_{n_2,\left[p\right]})}^T,\hfill \\ \hfill A_{\mathit{n},\left[p\right],21}& =(I_{n_1+p}\otimes P_{n_2,\left[p\right]}){\tilde{A}}_{\mathit{n},\left[p\right],21}{(P_{n_1,\left[p\right]}\otimes I_{n_2+p})}^T,\hfill \\ \hfill A_{\mathit{n},\left[p\right],22}& =(I_{n_1+p}\otimes P_{n_2,\left[p\right]}){\tilde{A}}_{\mathit{n},\left[p\right],22}{(I_{n_1+p}\otimes P_{n_2,\left[p\right]})}^T,\hfill \end{array}$$

and so

$$A_{\mathit{n},\left[p\right]}=P_{\mathit{n},\left[p\right]}{\tilde{A}}_{\mathit{n},\left[p\right]}{\left(P_{\mathit{n},\left[p\right]}\right)}^T,P_{\mathit{n},\left[p\right]}=\left[\begin{array}{cc}{P}_{n_1,\left[p\right]}\otimes I_{n_2+p}& \\ & I_{n_1+p}\otimes P_{n_2,\left[p\right]}\end{array}\right].$$

(62)

In view of the application of Theorem 1, we note that

$$\begin{array}{c}{P}_{\mathit{n},\left[p\right]}\in {\mathbb{R}}^{[(n_1+p-1)(n_2+p)+(n_1+p)(n_2+p-1)]\times 2(n_1+p)(n_2+p)},\hfill \end{array}$$

(63)

$$\begin{array}{c}{P}_{\mathit{n},\left[p\right]}{\left(P_{\mathit{n},\left[p\right]}\right)}^T=I_{(n_1+p-1)(n_2+p)+(n_1+p)(n_2+p-1)},\hfill \end{array}$$

(64)

$$\begin{array}{c}\underset{\mathit{n}\to \infty }{\lim }\frac{(n_1+p-1)(n_2+p)+(n_1+p)(n_2+p-1)}{2(n_1+p)(n_2+p)}=\underset{\mathit{n}\to \infty }{\lim }\left[\frac{n_1+p-1}{2(n_1+p)}+\frac{n_2+p-1}{2(n_2+p)}\right]=1.\hfill \end{array}$$

(65)

Lemma 5.

Let $p\ge 2$ and $n\ge 1$. Then,

$$\begin{array}{cccc}\hfill n^{-1}{K}_{n,\left[p\right]}& =T_{n+p}\left(f_p\right)+Q_{n,\left[p\right]},\hfill & \hfill \mathrm{rank}\left(Q_{n,\left[p\right]}\right)& \le 4p,\hfill \end{array}$$

(66)

$$\begin{array}{cccc}\hfill {\tilde{H}}_{n,\left[p\right]}& =T_{n+p}\left(g_p\right)+R_{n,\left[p\right]},\hfill & \hfill \mathrm{rank}\left(R_{n,\left[p\right]}\right)& \le 4p,\hfill \end{array}$$

(67)

$$\begin{array}{cccc}\hfill n{\tilde{M}}_{n,[p-1]}& =T_{n+p}\left(h_p\right)+S_{n,\left[p\right]},\hfill & \hfill \mathrm{rank}\left(S_{n,\left[p\right]}\right)& \le 4p,\hfill \end{array}$$

(68)

where

$$\begin{array}{cc}\hfill f_p\left(\theta \right)& =\sum _{k\in \mathbb{Z}}-{\beta }_{[2p+1]}^{\prime \prime }(p+1-k){\mathrm{e}}^{\mathrm{i}k\theta },\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill g_p\left(\theta \right)& =\sum _{k\in \mathbb{Z}}-{\beta }_{\left[2p\right]}^{\prime }(p-k){\mathrm{e}}^{\mathrm{i}k\theta },\hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill h_p\left(\theta \right)& =\sum _{k\in \mathbb{Z}}{\beta }_{[2p-1]}(p-k){\mathrm{e}}^{\mathrm{i}k\theta },\hfill \end{array}$$

(71)

and we note that the three series are actually finite sums because of (57).

Proof. 

For every $i,j=p+1,\dots ,n$, since $[-i+p+1,n-i+p+1]\supseteq [0,p+1]=\mathrm{supp}\left({\beta }_{\left[p\right]}\right)$ and $[-i+p,n-i+p]\supseteq [0,p]=\mathrm{supp}\left({\beta }_{[p-1]}\right)$, by (59) and (60) we obtain

$$\begin{array}{cc}\hfill {(K_{n,\left[p\right]})}_{ij}& ={\int }_0^1{B}_{j,\left[p\right]}^{\prime }\left(x\right)B_{i,\left[p\right]}^{\prime }\left(x\right)\mathrm{d}x=n^2{\int }_0^1{\beta }_{\left[p\right]}^{\prime }(nx-j+p+1){\beta }_{\left[p\right]}^{\prime }(nx-i+p+1)\mathrm{d}x\hfill \\ & =n{\int }_{-i+p+1}^{n-i+p+1}{\beta }_{\left[p\right]}^{\prime }(\tau +i-j){\beta }_{\left[p\right]}^{\prime }\left(\tau \right)\mathrm{d}\tau =n{\int }_{\mathbb{R}}{\beta }_{\left[p\right]}^{\prime }\left(\tau \right){\beta }_{\left[p\right]}^{\prime }(\tau +i-j)\mathrm{d}\tau \hfill \\ & =-n{\beta }_{[2p+1]}^{″}(p+1+i-j)=-n{\beta }_{[2p+1]}^{″}(p+1-i+j),\hfill \\ \hfill {({\tilde{H}}_{n,\left[p\right]})}_{ij}& ={\int }_0^1{B}_{j,\left[p\right]}^{\prime }\left(x\right)B_{i,[p-1]}\left(x\right)\mathrm{d}x=n{\int }_0^1{\beta }_{\left[p\right]}^{\prime }(nx-j+p+1){\beta }_{[p-1]}(nx-i+p)\mathrm{d}x\hfill \\ & ={\int }_{-i+p}^{n-i+p}{\beta }_{\left[p\right]}^{\prime }(\tau +i-j+1){\beta }_{[p-1]}\left(\tau \right)\mathrm{d}\tau ={\int }_{\mathbb{R}}{\beta }_{[p-1]}\left(\tau \right){\beta }_{\left[p\right]}^{\prime }(\tau +i-j+1)\mathrm{d}\tau \hfill \\ & ={\beta }_{\left[2p\right]}^{\prime }(p+i-j+1)=-{\beta }_{\left[2p\right]}^{\prime }(p-i+j),\hfill \\ \hfill {({\tilde{M}}_{n,[p-1]})}_{ij}& ={\int }_0^1{B}_{j,[p-1]}\left(x\right)B_{i,[p-1]}\left(x\right)\mathrm{d}x={\int }_0^1{\beta }_{[p-1]}(nx-j+p){\beta }_{[p-1]}(nx-i+p)\mathrm{d}x\hfill \\ & =n^{-1}{\int }_{-i+p}^{n-i+p}{\beta }_{[p-1]}(\tau +i-j){\beta }_{[p-1]}\left(\tau \right)\mathrm{d}\tau =n^{-1}{\int }_{\mathbb{R}}{\beta }_{[p-1]}\left(\tau \right){\beta }_{[p-1]}(\tau +i-j)\mathrm{d}\tau \hfill \\ & =n^{-1}{\beta }_{[2p-1]}(p+i-j)=n^{-1}{\beta }_{[2p-1]}(p-i+j).\hfill \end{array}$$

Thus,

$$\begin{array}{cc}\hfill {\left[{\left(n^{-1}{K}_{n,\left[p\right]}\right)}_{ij}\right]}_{i,j=p+1}^n& ={[-{\beta }_{[2p+1]}^{\prime \prime }(p+1-i+j)]}_{i,j=p+1}^n=T_{n-p}\left(f_p\right),\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill [({\tilde{H}}_{n,\left[p\right]}{{)}_{ij}]}_{i,j=p+1}^n& ={[-{\beta }_{\left[2p\right]}^{\prime }(p-i+j)]}_{i,j=p+1}^n=T_{n-p}\left(g_p\right),\hfill \end{array}$$

(73)

$$\begin{array}{cc}\hfill [(n{\tilde{M}}_{n,[p-1]}{{)}_{ij}]}_{i,j=p+1}^n& ={\left[{\beta }_{[2p-1]}(p-i+j)\right]}_{i,j=p+1}^n=T_{n-p}\left(h_p\right).\hfill \end{array}$$

(74)

It follows from (72) that the principal submatrix of $n^{-1}{K}_n^{\left[p\right]}-T_{n+p}\left(f_p\right)$ corresponding to the row and column indices $i,j=p+1,\dots ,n$ is the zero matrix, which implies (66). Similarly, (73) and (74) imply (67) and (68), respectively. ☐

Theorem 6.

Let $p\ge 2$, let $\nu =({\nu }_1,{\nu }_2)\in {\mathbb{Q}}^2$ be a vector with positive components, and assume that $\mathit{n}=n\nu$ (it is understood that n varies in the infinite subset of $\mathbb{N}$ such that $\mathit{n}=n\mathit{\nu }\in {\mathbb{N}}^2$). Then,

$${\left\{A_{\mathit{n},\left[p\right]}\right\}}_n{\sim }_{\sigma ,\lambda }\kappa \left(\mathit{\theta }\right)=\left[\begin{array}{cc}{\displaystyle \frac{{\nu }_2}{{\nu }_1}}{h}_p\left({\theta }_1\right)f_p\left({\theta }_2\right)& -g_p\left({\theta }_1\right)\overline{g_p\left({\theta }_2\right)}\\ -\overline{g_p\left({\theta }_1\right)}{g}_p\left({\theta }_2\right)& {\displaystyle \frac{{\nu }_1}{{\nu }_2}}{f}_p\left({\theta }_1\right)h_p\left({\theta }_2\right)\end{array}\right].$$

Proof. 

The thesis follows immediately from Theorem 1 and (62)–(65) as soon as we have proven that

$${\left\{{\tilde{A}}_{\mathit{n},\left[p\right]}\right\}}_n{\sim }_{\sigma ,\lambda }\kappa \left(\mathit{\theta }\right).$$

(75)

We show that

$${\left\{{\tilde{A}}_{\mathit{n},\left[p\right],ij}\right\}}_n{\sim }_{\mathrm{GLT}}{\kappa }_{ij}\left(\mathit{\theta }\right),i,j=1,2.$$

(76)

Once this is done, the thesis (75) follows immediately from Theorem 5 and GLT 1 as the matrix ${\tilde{A}}_{\mathit{n},\left[p\right]}$ is symmetric. Actually, we only prove (76) for $(i,j)=(1,2)$ because the proof for the other pairs of indices $(i,j)$ is conceptually the same. Setting $\mathit{p}=(p,p)$ and keeping in mind the assumption $\mathit{n}=n\nu$, by Lemma 5 and Equations (5), (51) and (52), we have

$$\begin{array}{cc}\hfill {\tilde{A}}_{\mathit{n},\left[p\right],12}& =-{\tilde{H}}_{n_1,\left[p\right]}\otimes {\left({\tilde{H}}_{n_2,\left[p\right]}\right)}^T=-(T_{n_1+p}\left(g_p\right)+R_{n_1,\left[p\right]})\otimes {(T_{n_2+p}\left(g_p\right)+R_{n_2,\left[p\right]})}^T\hfill \\ & =-(T_{n_1+p}\left(g_p\right)+R_{n_1,\left[p\right]})\otimes (T_{n_2+p}\left(\overline{g_p}\right)+{\left(R_{n_2,\left[p\right]}\right)}^T)\hfill \\ & =-(T_{\mathit{n}+\mathit{p}}\left(g_p\left({\theta }_1\right)\overline{g_p\left({\theta }_2\right)}\right)+T_{n_1+p}\left(g_p\right)\otimes {\left(R_{n_2,\left[p\right]}\right)}^T+R_{n_1,\left[p\right]}\otimes {\left({\tilde{H}}_{n_2,\left[p\right]}\right)}^T)\hfill \\ & =T_{\mathit{n}+\mathit{p}}\left({\kappa }_{12}\right)+V_{\mathit{n},\left[p\right]},\hfill \end{array}$$

where $\mathrm{rank}\left(V_{\mathit{n},\left[p\right]}\right)\le 4p(n_1+p)+4p(n_2+p)$. Thus, ${\left\{V_{\mathit{n},\left[p\right]}\right\}}_n{\sim }_{\sigma }0$ by Proposition 1, and (76) (for $(i,j)=(1,2)$) follows from GLT 2 and GLT 3. ☐

Remark 5.

Using (58), it is not difficult to see that the functions $f_p\left(\theta \right)$ and $g_p\left(\theta \right)$ in (69) and (70) can be expressed in terms of $h_p\left(\theta \right)$ as follows:

$$f_p\left(\theta \right)=(2-2\cos \theta )h_p\left(\theta \right),g_p\left(\theta \right)=({\mathrm{e}}^{-\mathrm{i}\theta }-1)h_p\left(\theta \right).$$

Therefore, the $2\times 2$ matrix-valued function $\kappa \left(\mathit{\theta }\right)$ appearing in Theorem 6 can be simplified as follows:

$$\kappa \left(\mathit{\theta }\right)=\frac{1}{{\nu }_1{\nu }_2}{h}_p\left({\theta }_1\right)h_p\left({\theta }_2\right)\left[\begin{array}{c}{\nu }_2({\mathrm{e}}^{\mathrm{i}{\theta }_2}-1)\\ -{\nu }_1({\mathrm{e}}^{\mathrm{i}{\theta }_1}-1)\end{array}\right]\left[\begin{array}{cc}{\nu }_2({\mathrm{e}}^{-\mathrm{i}{\theta }_2}-1)& -{\nu }_1({\mathrm{e}}^{-\mathrm{i}{\theta }_1}-1)\end{array}\right].$$

In particular, $\det \left(\kappa \right(\mathit{\theta }\left)\right)=0$ for all θ. According to the informal meaning behind the spectral distribution ${\left\{A_{\mathit{n},\left[p\right]}\right\}}_n{\sim }_{\lambda }\kappa \left(\mathit{\theta }\right)$ reported in Remark 1, this means that, for large n, one half of the eigenvalues of $A_{\mathit{n},\left[p\right]}$ are approximately zero and one half is given by a uniform sampling over ${[-\pi ,\pi ]}^2$ of

$$\mathrm{trace}\left(\kappa \left(\mathit{\theta }\right)\right)=\frac{1}{{\nu }_1{\nu }_2}{h}_p\left({\theta }_1\right)h_p\left({\theta }_2\right)\left[{\nu }_1^2(2-2\cos {\theta }_1)+{\nu }_2^2(2-2\cos {\theta }_1)\right].$$

8. Conclusions

We have illustrated through specific examples the applicative interest of the theory of block GLT sequences and of its multivariate version, thus bringing to completion the purely theoretical work [34]. It should be said, however, that the theory of GLT sequences is still incomplete. In particular, besides filling in the details of the theory of multilevel block GLT sequences—the results of Section 5 have been obtained as a combination of the results in [14,34], but formal proofs of them are still missing and will be the subject of a future paper—, it will be necessary to develop the theory of the so-called reduced GLT sequences, as explained in ([13] Chapter 11).