# A note on the spectral analysis of matrix-sequences via GLT momentary
symbols: from all-at-once solution of parabolic problems to distributed
fractional order matrices

Matthias Bolten  Bergische Universität Wuppertal
Department of Mathematics and Computer Science
Gaußstraße 20, 42119 Wuppertal,
Germany<EMAIL_ADDRESS>Sven-Erik Ekström  Uppsala University
Department of Information Technology
Lägerhyddsv. 2, hus 2, SE-751 05, Uppsala
Sweden<EMAIL_ADDRESS>Isabella Furci  Bergische Universität
Wuppertal
Department of Mathematics and Computer Science
Gaußstraße 20, 42119, Wuppertal
Germany<EMAIL_ADDRESS>Stefano Serra-Capizzano  Insubria University
Department of Science and High Technology
via Valleggio 11, 22100, Como
Italy<EMAIL_ADDRESS>Uppsala University
Department of Information Technology
Lägerhyddsv. 2, hus 2, SE-751 05, Uppsala
Sweden<EMAIL_ADDRESS>
###### Abstract.

The first focus of this paper is the characterization of the spectrum and the
singular values of the coefficient matrix stemming from the discretization
with space-time grid for a parabolic diffusion problem and from the
approximation of distributed order fractional equations. For this purpose we
will use the classical GLT theory and the new concept of GLT momentary
symbols. The first permits to describe the singular value or eigenvalue
asymptotic distribution of the sequence of the coefficient matrices, the
latter permits to derive a function, which describes the singular value or
eigenvalue distribution of the matrix of the sequence, even for small matrix-
sizes but under given assumptions. The note is concluded with a list of open
problems, including the use of our machinery in the study of iteration
matrices, especially those concerning multigrid-type techniques.

###### keywords:

Toeplitz matrices; Asymptotic distribution of eigenvalues; Numerical solution
of discretized equations for boundary value problems involving PDEs;
Fractional partial differential equations

###### 1991 Mathematics Subject Classification:

15B05; 34L20; 65N22; 35R11

##  1\. Introduction and notation

As well known, many practical applications require to solve numerically linear
systems of Toeplitz kind and of large dimensions. As a consequence a number of
iterative techniques, such as preconditioned Krylov methods, multigrid
procedures, and sophisticated combination of them have been designed (see [9,
24] and the references therein). Linear systems with Toeplitz coefficient
matrices of large dimension arise when dealing with the numerical solution of
(integro-)differential equations and of problems with Markov chains. More
recently, new examples of real world problems have emerged.  The first focus
of this paper is the characterization of the spectrum and the singular values
of the coefficient matrix stemming from the discretization with a space-time
grid for a parabolic diffusion problem. More specifically, we consider the
diffusion equation in one space dimension,

$u_{t}=u_{xx},\quad x\in(a,b),\ t\in[0,T],$

and we approximate our parabolic model problem on a rectangular space-time
grid consisting of $N_{t}$ time intervals and $N_{x}$ space intervals.

The second focus concerns the matrix-sequences involved with the
discretization of distributed order fractional differential equations (FDEs)
which have gained a lot of attention. Owing to the nonlocal nature of
fractional operators, independently of the locality of the approximation
methods, the matrix structures are dense and under assumptions of uniform
step-sizing and of constant coefficients in the involved operators, the
matrices are again of Toeplitz type (unilevel, or multilevel according to the
dimensionality of the considered domains).

When the fractional order is fixed, the spectral analysis of such matrices
(conditioning, extremal eigenvalues etc) can be performed, by exploiting the
well-established analysis of the spectral features of Toeplitz matrix-
sequences generated by Lebesgue integrable functions and the more recent
Generalized Locally Toeplitz (GLT) theory [20]; see for instance [14, 15].
However in the case of the numerical approximation of distributed-order
fractional operators, also the spectral analysis of the resulting matrices is
more involved. We recall that distributed-order FDEs can be interpreted as a
parallel distribution of derivatives of fractional orders, whose most
immediate application consists in the physical modeling of systems
characterized by a superposition of different processes operating in parallel.
As an example, we mention the application of fractional distributed-order
operators as a tool for accounting memory effects in composite materials [11]
or multi-scale effects [10]. For a detailed review on the topic we refer the
reader to [13].

In order to study the involved structured linear systems of both integral and
differential equations, we will use the classical theory of GLT matrix-
sequences [20, 21] and the new concept of GLT momentary symbols. The first
permits to describe the singular value or eigenvalue asymptotic distribution
of the sequence of the coefficient matrices, the latter permits to derive a
function, which describes the singular value or eigenvalue distribution of a
fixed matrix of the sequence, even for small matrix-sizes, but under given
assumptions.

This paper is organized as follows. The remaining part of this section is
devoted to definitions, notation, and to the necessary background for our
analysis: in particular we provide a formal definition of GLT momentary
symbols. Section 2 is devoted to setting up the problem and to derive the
relevant matrix structures. The distributional analysis both for the
eigenvalues and singular values is the main focus of Subsection 2.3, while
Section 3 contains similar results for specific matrix-structures with
generating function depending on the matrix-size and which arise in the
context of fractional differential equations with distributed orders. Section
4 contains conclusions and a list of open problems, including the use of our
machinery in the study of iteration matrices, especially those concerning
multigrid-type techniques.

###  1.1. Background and definitions

Throughout this paper, we will use the following notations. Let
${f}:G\to\mathbb{C}$ be a function belonging to $L^{1}(G)$, with
$G\subseteq\mathbb{R}^{\ell}$, $\ell\geq 1$, a measurable set. We denote by
$\\{A_{n}\\}_{n}$ the matrix-sequence, whose elements are given by the
matrices $A_{n}$ of dimension $n\times n$. Let $s,d\in\mathbb{N}$. Let
$\mathbf{n}=(n_{1},n_{2},\dots,n_{d})$ be a multi-index, we indicate by
$\\{A_{\mathbf{n}}\\}_{\mathbf{n}}$, the $d$-level $s\times s$ block matrix-
sequence, whose elements are the matrices $A_{\mathbf{n}}$ of size
$d=d(\mathbf{n},s)=sn_{1}n_{2}\cdots n_{d}$.

###  1.2. Toeplitz and circulant matrix-sequences

In the following we report the main background concerning the concepts of
Toeplitz and circulant matrices, for simplicity, in the scalar unilevel
setting. We only provide the generalization in the block multilevel case of
the results that will be exploited for the purpose of the paper.

######  Definition 1.1.

An $n\times n$ Toeplitz matrix $A_{n}$ is a matrix that has equal entries
along each diagonal, and can be written as

$A_{n}=\left[a_{i-j}\right]_{i,j=1}^{n}=\left[\begin{smallmatrix}a_{0}&a_{-1}&a_{-2}&\cdots&\cdots&a_{1-n}\vphantom{\ddots}\\\
a_{1}&\ddots&\ddots&\ddots&&\vdots\\\
a_{2}&\ddots&\ddots&\ddots&\ddots&\vdots\\\
\vdots&\ddots&\ddots&\ddots&\ddots&a_{-2}\\\
\vdots&&\ddots&\ddots&\ddots&a_{-1}\\\
a_{n-1}&\cdots&\cdots&a_{2}&a_{1}&a_{0}\vphantom{\ddots}\\\
&\end{smallmatrix}\right],\ \ \ a_{j}\in\mathbb{C},\ j=1-n,\ldots,n-1.$

In the following we focus on the two important sub-classes given by the
Toeplitz matrices $T_{n}(f)\in\mathbb{C}^{n\times n}$ and the circulant
matrices $C_{n}(f)\in\mathbb{C}^{n\times n}$, associated with a function $f$,
called the generating function.

######  Definition 1.2.

Given $f$ belonging to $L^{1}([-\pi,\pi])$ and periodically extended to the
whole real line, the matrix $T_{n}(f)$ is defined as

$T_{n}(f)=\left[\hat{f}_{i-j}\right]_{i,j=1}^{n},$

where

$\hat{f}_{k}\coloneqq\frac{1}{2\pi}\int_{-\pi}^{\pi}\\!\\!f(\theta)\,\mathrm{e}^{-k\mathbf{i}\theta}\\!\mathrm{d}\theta,\qquad
k\in\mathbb{Z},\qquad\mathbf{i}^{2}=-1,$ (1.1)

are the Fourier coefficients of $f$, and

$f(\theta)=\\!\\!\sum_{k=-\infty}^{\infty}\\!\\!\hat{f}_{k}\mathrm{e}^{k\mathbf{i}\theta},$

is the Fourier series of $f$.

######  Definition 1.3.

Let the Fourier coefficients of a given function ${f}\in L^{1}([-\pi,\pi])$ be
defined as in formula (1.1). Then, we can define the $n\times n$ circulant
matrix $C_{n}(f)$ associated with $f$, as

$C_{n}(f)=\\!\\!\\!\\!\\!\\!\sum_{j=-(n-1)}^{n-1}\\!\\!\\!\\!\\!\\!\hat{a}_{j}Z_{n}^{j}=\mathbb{F}_{n}D_{n}(f)\mathbb{F}_{n}^{*},$
(1.2)

where ∗ denotes the transpose conjugate, $Z_{n}$ is the $n\times n$ matrix
defined by

$\left(Z_{n}\right)_{ij}=\begin{cases}1,&\text{if }\mathrm{mod}(i-j,n)=1,\\\
0,&\text{otherwise}.\end{cases}$

Moreover,

$D_{n}(f)=\operatorname*{diag}\left(s_{n}(f(\theta_{j,n}^{c}))\right),\quad
j=1,\ldots,n,$

where

$\theta_{j,n}^{c}=\frac{(j-1)2\pi}{n},\quad j=1,\ldots,n,$ (1.3)

and $s_{n}(f(\theta))$ is the $n$th Fourier sum of $f$ given by

$s_{n}(f({\theta}))=\sum_{k=1-n}^{n-1}\hat{f}_{k}\mathrm{e}^{k\mathbf{i}\theta}.$

The matrix $\mathbb{F}_{n}$ is the so called Fourier matrix of order $n$,
given by

$\displaystyle(\mathbb{F}_{n})_{i,j}=\frac{1}{\sqrt{n}}\mathrm{e}^{\mathbf{i}(i-1)\theta_{j,n}^{c}},\quad
i,j=1,\ldots,n.$

In the case of the Fourier matrix, we have
$\mathbb{F}_{n}\mathbb{F}_{n}^{*}=\mathbb{I}_{n}$, that is $\mathbb{F}_{n}$ is
complex-symmetric and unitary, with $\mathbb{I}_{n}$ being the identity of
size $n$.

The proof of the second equality in (1.2), which implies that the columns of
the Fourier matrix $\mathbb{F}_{n}$ are the eigenvectors of $C_{n}(f)$, can be
found in [20, Theorem 6.4].

Note that, from the definition follows that if ${f}$ is a trigonometric
polynomial of fixed degree less than $n$, the entries of $D_{n}(f)$ are the
eigenvalues of $C_{n}(f)$, explicitly given by sampling the generating
function $f$ using the grid $\theta_{j,n}^{c}$.

$\begin{split}\lambda_{j}(C_{n}(f))&=f\left(\theta_{j,n}^{c}\right),\quad
j=1,\ldots,n,\\\
D_{n}(f)&=\operatorname*{diag}\left(f\left(\theta_{j,n}^{c}\right)\right),\quad
j=1,\ldots,n.\end{split}$

The type of domain (either one-dimensional $[-\pi,\pi]$ or d-dimensional
$[-\pi,\pi]^{d}$) and codomain (either the complex field or the space of
$s\times s$ complex matrices) of $f$ gives rise to different kinds of Toeplitz
matrices, see Table 1 for a complete overview.

Table 1. Different types of generating function and the associated Toeplitz matrix. Type of generating function | Associated Toeplitz matrix
---|---
univariate scalar | $f(\theta):[-\pi,\pi]\to\mathbb{C}$ | unilevel scalar | $T_{n}(f)\in\mathbb{C}^{n\times n}$
$d$-variate scalar | $f(\boldsymbol{\theta}):[-\pi,\pi]^{d}\to\mathbb{C}$ |  $d$-level scalar | $T_{\mathbf{n}}(f)\in\mathbb{C}^{d(\mathbf{n},1)\times d(\mathbf{n},1)}$
univariate matrix-valued | $\mathbf{f}(\theta):[-\pi,\pi]\to\mathbb{C}^{s\times s}$ | unilevel block | $T_{n}(\mathbf{f})\in\mathbb{C}^{d({n},s)\times d({n},s)}$
$d$-variate matrix-valued | $\mathbf{f}(\boldsymbol{\theta}):[-\pi,\pi]^{d}\to\mathbb{C}^{s\times s}$ |  $d$-level block | $T_{\mathbf{n}}(\mathbf{f})\in\mathbb{C}^{d(\mathbf{n},s)\times d(\mathbf{n},s)}$

In particular, we provide the definition of a $d$-level $s\times s$ block
Toeplitz matrices $T_{\mathbf{n}}(\mathbf{f})$ starting from $d$-variate
matrix-valued function
$\mathbf{f}:[-\pi,\pi]^{d}\rightarrow\mathbb{C}^{s\times s}$, with with
$\mathbf{f}\in L^{1}([-\pi,\pi]^{d})$.

######  Definition 1.4.

Given a function $\mathbf{f}:[-\pi,\pi]^{d}\rightarrow\mathbb{C}^{s\times s}$
its Fourier coefficients are given by

$\hat{\mathbf{f}}_{\mathbf{k}}\coloneqq\frac{1}{(2\pi)^{d}}\int_{[-\pi,\pi]^{d}}\mathbf{f}(\boldsymbol{\theta})\mathrm{e}^{-\mathbf{i}\left\langle{\mathbf{k}},\boldsymbol{\theta}\right\rangle}\mathrm{d}\boldsymbol{\theta}\in\mathbb{C}^{s\times
s},\qquad\mathbf{k}=(k_{1},\ldots,k_{d})\in\mathbb{Z}^{d},$

where $\boldsymbol{\theta}=(\theta_{1},\ldots,\theta_{d})$,
$\left\langle\mathbf{k},\boldsymbol{\theta}\right\rangle=\sum_{i=1}^{d}k_{i}\theta_{i}$,
and the integrals of matrices are computed elementwise. The associated
generating function an be defined via its Fourier series as

$\mathbf{f}(\boldsymbol{\theta})=\sum_{\mathbf{k}\in\mathbb{Z}^{d}}\hat{\mathbf{f}}_{\mathbf{k}}\mathrm{e}^{\mathbf{i}\left\langle{\mathbf{k}},\boldsymbol{\theta}\right\rangle}.$

The $d$-level $s\times s$ block Toeplitz matrix associated with $\mathbf{f}$
is the matrix of dimension $d({\mathbf{n},s})$, where
$\mathbf{n}=(n_{1},\ldots,n_{d})$, given by

$T_{\mathbf{n}}(\mathbf{f})=\sum_{\mathbf{e}-\mathbf{n}\leq\mathbf{k}\leq\mathbf{n}-\mathbf{e}}T_{n_{1}}(\mathrm{e}^{\mathbf{i}k_{1}\theta_{1}})\otimes\cdots\otimes
T_{n_{d}}(\mathrm{e}^{\mathbf{i}k_{d}\theta_{1}})\otimes\hat{\mathbf{f}}_{\mathbf{k}},$

where $\mathbf{e}$ is the vector of all ones and where
$\mathbf{s}\leq\mathbf{t}$ means that $s_{j}\leq t_{j}$ for any
$j=1,\ldots,d$.

######  Definition 1.5.

If $\textbf{n}\in\mathbb{N}^{d}$ and
$\textbf{a}:[0,1]^{d}\to\mathbb{C}^{s\times s}$, we define the n-th $d$-level
and $s\times s$ block diagonal sampling matrix as the following multilevel
block diagonal matrix of dimension $d(\textbf{n},s)$:

$D_{\textbf{n}}(\textbf{a})=\operatorname*{diag}_{\textbf{e}\leq\textbf{j}\leq\textbf{n}}\textbf{a}\left(\frac{\textbf{j}}{\textbf{n}}\right),$

where we recall that $\textbf{e}\leq\textbf{j}\leq\textbf{n}$ means that j
varies from e to n following the lexicographic ordering.

The following result provides an important relation between tensor products
and multilevel Toeplitz matrices.

######  Lemma 1.6.

[21] Let $f_{1},\dots,f_{d}\in L^{1}([-\pi,\pi])$,
$\mathbf{n}=(n_{1},n_{2},\dots,n_{d})\in\mathbb{N}^{d}$. Then,

$T_{n_{1}}(f_{1})\otimes\dots\otimes
T_{n_{d}}(f_{d})=T_{\mathbf{n}}(f_{1}\otimes\dots\otimes f_{d}),$

where the Fourier coefficients of $f_{1}\otimes\dots\otimes f_{d}$ are given
by

$(f_{1}\otimes\dots\otimes
f_{d})_{\mathbf{k}}=(f_{1})_{k_{1}}\dots(f_{d})_{k_{d}},\quad\mathbf{k}\in\mathbb{Z}^{d}.$

###  1.3. Asymptotic distributions

In this subsection we introduce the definition of asymptotic distribution in
the sense of the eigenvalues and of the singular values, first for a generic
matrix-sequence $\\{A_{n}\\}_{n}$, and then we report specific results
concerning the distributions of Toeplitz and circulant matrix-sequences.
Finally, we recall the notion of GLT algebra and we introduce a general notion
of GLT momentary symbols. We remind that in a more specific and limited
setting, the notion of momentary symbols is given in [6]: here we generalize
the definition in [6].

######  Definition 1.7.

[20, 21, 22, 30] Let $f,{\mathfrak{f}}:G\to\mathbb{C}$ be measurable
functions, defined on a measurable set $G\subset\mathbb{R}^{\ell}$ with
$\ell\geq 1$, $0<\mu_{\ell}(G)<\infty$. Let $\mathcal{C}_{0}(\mathbb{K})$ be
the set of continuous functions with compact support over
$\mathbb{K}\in\\{\mathbb{C},\mathbb{R}_{0}^{+}\\}$ and let
$\\{A_{\mathbf{n}}\\}_{\mathbf{n}}$, be a sequence of matrices with
eigenvalues $\lambda_{j}(A_{\mathbf{n}})$, $j=1,\ldots,{d_{\mathbf{n}}}$ and
singular values $\sigma_{j}(A_{\mathbf{n}})$, $j=1,\ldots,{d_{\mathbf{n}}}$.
Then,

  * •

The matrix-sequence $\\{A_{\mathbf{n}}\\}_{\mathbf{n}}$ is distributed as $f$
in the sense of the singular values, and we write

$\displaystyle\\{A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\sigma}f,$

if the following limit relation holds for all
$F\in\mathcal{C}_{0}(\mathbb{R}_{0}^{+})$:

$\lim_{\mathbf{n}\to\infty}\frac{1}{{d_{\mathbf{n}}}}\sum_{j=1}^{{d_{\mathbf{n}}}}F(\sigma_{j}(A_{\mathbf{n}}))=\frac{1}{\mu_{\ell}(G)}\int_{G}F({|f(\boldsymbol{\theta})|})\,\\!\mathrm{d}{\boldsymbol{\theta}}.$
(1.4)

The function $f$ is called the singular value symbol which describes
asymptotically the singular value distribution of the matrix-sequence
$\\{A_{\mathbf{n}}\\}_{\mathbf{n}}$.

  * •

The matrix-sequence $\\{A_{\mathbf{n}}\\}_{\mathbf{n}}$ is distributed as
$\mathfrak{f}$ in the sense of the eigenvalues, and we write

$\\{A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\lambda}\mathfrak{f},$

if the following limit relation holds for all
$F\in\mathcal{C}_{0}(\mathbb{C})$:

$\lim_{\mathbf{n}\to\infty}\frac{1}{{d_{\mathbf{n}}}}\sum_{j=1}^{{d_{\mathbf{n}}}}F(\lambda_{j}(A_{\mathbf{n}}))=\frac{1}{\mu_{\ell}(G)}\int_{G}\displaystyle
F({\mathfrak{f}(\boldsymbol{\theta})})\,\\!\mathrm{d}{\boldsymbol{\theta}}.$
(1.5)

The function $\mathfrak{f}$ is called the eigenvalue symbol which describes
asymptotically the eigenvalue distribution of the matrix-sequence
$\\{A_{\mathbf{n}}\\}_{\mathbf{n}}$.

######  Remark 1.8.

Note that, if $A_{\mathbf{n}}$ is normal for any $\mathbf{n}$ or at least
definitely, then $\\{A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\sigma}f$ and
$\\{A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\lambda}\mathfrak{f}$ imply that
$f=\mathfrak{f}$. Of course this is true for Hermitian Toeplitz matrix-
sequences as emphasized in Theorem 1.9 and Theorem 1.10.

Moreover, considering the case $d=1$, if $f$ (or $\mathfrak{f}$) is smooth
enough, then the informal interpretation of the limit relation (1.4) (or
(1.5)) is that, for $n$ sufficiently large, the $n$ singular values (or
eigenvalues) of $A_{n}$ can be approximated by a sampling of $|f(\theta)|$ (or
$\mathfrak{f}(\theta)$) on an equispaced grid of the interval $G$, up to the
presence of possibly $o(n)$ outliers. It is worthy to notice that in most of
the Toeplitz and PDE/FDE applications the number of actual outliers is often
limited to $O(1)$ with often a small number of outliers (see [20, 21, 3, 2]
and references therein).

The generalization of Definition 1.7 and Remark 1.8 to the block setting and
multilevel block setting can be found in [3, 2] and in the references therein.

In the case where the matrix-sequence is a Toeplitz matrix-sequence generated
by a function, the singular value distribution and the spectral distribution
have been well studied in the past few decades. At the beginning Szegő in [22]
showed that the eigenvalues of the Toeplitz matrix $T_{n}(f)$ generated by
real-valued $f\in L^{\infty}([-\pi,\pi])$ are asymptotically distributed as
$f$. Moreover, under the same assumption on $f$, Avram and Parter [1, 25]
proved that the singular values of $T_{n}(f)$ are distributed as $|f|$. This
result has been undergone many generalizations and extensions among the years
(see [20, 21, 3, 2] and the references therein).

The generalized Szegő theorem that describes the singular value and spectral
distribution of Toeplitz sequences generated by a scalar $f\in
L^{1}([-\pi,\pi])$ is given as follows [31].

######  Theorem 1.9.

Suppose $f\in L^{1}([-\pi,\pi])$. Let $T_{n}(f)$ be the Toeplitz matrix
generated by $f$. We have

$\\{T_{n}(f)\\}_{n}\sim_{\sigma}f.$

Moreover, if $f$ is real-valued almost everywhere (a.e.), then

$\\{T_{n}(f)\\}_{n}\sim_{\lambda}f.$

Tilli [29] generalized the proof to the block-Toeplitz setting and we report
the extension of the eigenvalue result to the case of multivariate Hermitian
matrix-valued generating functions.

######  Theorem 1.10.

Suppose $\mathbf{f}\in L^{1}([-\pi,\pi]^{d},s)$ with positive integers $d,s$.
Let $T_{\bf n}(\mathbf{f})$ be the Toeplitz matrix generated by $\mathbf{f}$.
We have

$\\{T_{\bf n}(\mathbf{f})\\}_{{\bf n}}\sim_{\sigma}~{}\mathbf{f}.$

Moreover, if $\mathbf{f}$ is a Hermitian matrix-valued function a.e., then,

$\\{T_{\bf n}(\mathbf{f})\\}_{{\bf n}}\sim_{\lambda}~{}\mathbf{f}.$

Concerning the circulant matrix-sequences, though the eigenvalues of a $C_{\bf
n}(\textbf{f})$ are explicitly known, a result like Theorem 1.9 and Theorem
1.10 does not hold for sequences $\left\\{C_{\bf n}(\textbf{f})\right\\}_{{\bf
n}}$ in general. Indeed, the Fourier sum of f converges to f under quite
restrictive assumptions (see [32]). In particular, if f belongs to the Dini-
Lipschitz class, then $\left\\{C_{\bf n}(\textbf{f})\right\\}_{{\bf
n}}\sim_{\lambda}\textbf{f}$, (see [19] for more relationships between
circulant sequences and spectral distribution results).

###  1.4. Matrix algebras

A part from the circulant algebra, introduced in Section 1.2, we recall that
other particular matrix algebras have interesting properties and can be
exploited for our purpose. In particular, we mention the well-known
$\tau$-algebras, see [7] and references therein. Here, we restrict the
analysis to the case of the matrix algebras ${\tau_{\varepsilon,\varphi}}$,
introduced in [7], where an element of the algebra is a matrix

$\displaystyle
T_{n,\epsilon,\varphi}(g)=\left[\begin{array}[]{ccccc}a+\varepsilon b&b\\\
b&a&b\\\ &\ddots&\ddots&\ddots\\\ &&b&a&b\\\ &&&b&a+\varphi
b\end{array}\right],$ (1.11)

We can associate to this matrix a function $g$ of the form
$g(\theta)=a+2b\cos\theta$. For some values of $\varepsilon$ and $\varphi$ the
exact eigenvalues of $T_{n,\epsilon,\varphi}(g)$ are given by sampling with
specific grids; for detailed examples see [6] and [18] for asymptotic results.

In Table 2 we provide the proper grids $\theta^{(\varepsilon,\varphi)}_{j,n}$
and $\Theta_{i,j,n}^{(\varepsilon,\varphi)}$ to give the exact eigenvalues and
eigenvectors respectively for $\varepsilon,\varphi\in\\{-1,0,1\\}$.

Table 2. Grids for $\tau_{\varepsilon,\varphi}$-algebras, $\varepsilon,\varphi\in\\{-1,0,1\\}$; $\theta_{j,n}^{(\varepsilon,\varphi)}$ and $\Theta_{j,n}^{(\varepsilon,\varphi)}$ are the grids used to compute the eigenvalues and eigenvectors, respectively. The standard naming convention (dst-* and dct-*) in parenthesis; see, e.g., [12, Appendix 1]. |  | $\theta_{j,n}^{(\varepsilon,\varphi)}$ |  | $\Theta_{i,j,n}^{(\varepsilon,\varphi)}$
---|---|---|---|---
$\varepsilon$ $\varphi$ |  -1 |  0 |  1 |  -1, 0, 1
-1 |  (dst-2) $\frac{j\pi}{n}$ |  (dst-6) $\frac{j\pi}{n+1/2}$ |  (dst-4) $\frac{(j-1/2)\pi}{n}$ | $(i-1/2)\theta_{j,n}^{(\varepsilon,\varphi)}$
0 |  (dst-5) $\frac{j\pi}{n+1/2}$ |  (dst-1) $\frac{j\pi}{n+1}$ |  (dst-7) $\frac{(j-1/2)\pi}{n+1/2}$ | $i\theta_{j,n}^{(\varepsilon,\varphi)}$
1 |  (dct-4) $\frac{(j-1/2)\pi}{n}$ |  (dct-8) $\frac{(j-1/2)\pi}{n+1/2}$ |  (dct-2) $\frac{(j-1)\pi}{n}$ | $(i-1/2)\theta_{j,n}^{(\varepsilon,\varphi)}+\frac{\pi}{2}$

Since all grids $\theta_{j,n}^{(\varepsilon,\varphi)}$ associated with
$\tau_{\varepsilon,\varphi}$-algebras where
$\varepsilon,\varphi\in\\{-1,0,1\\}$ are uniformly spaced grids, we know that

$\displaystyle\theta_{j,n}^{(1,1)}$ $\displaystyle<\theta_{j,n}^{(0,1)}$
$\displaystyle=\theta_{j,n}^{(1,0)}$ $\displaystyle<\theta_{j,n}^{(-1,1)}$
$\displaystyle=\theta_{j,n}^{(1,-1)}$ $\displaystyle<\theta_{j,n}^{(0,0)}$
$\displaystyle<\theta_{j,n}^{(-1,0)}$ $\displaystyle=\theta_{j,n}^{(0,-1)}$
$\displaystyle<\theta_{j,n}^{(-1,-1)},\qquad\qquad\forall j=1,\ldots,n.$

###  1.5. Theory of Generalized Locally Toeplitz (GLT) sequences

In this subsection we will introduce the main properties from the theory of
Generalized Locally Toeplitz (GLT) sequences and the practical features, which
are sufficient for our purposes, see [20, 21, 3, 2].

In particular, we consider the multilevel and block setting with $d$ being the
number of levels.

GLT1:

Each GLT sequence has a singular value symbol
$\mathbf{f}(\boldsymbol{\theta,x})$ which is measurable according to the
Lebesgue measure and according to the second item in Definition 1.7 with
$\ell=2d$. In addition, if the sequence is Hermitian, then the distribution
also holds in the eigenvalue sense.

We specify that a GLT sequence $\\{A_{\mathbf{n}}\\}_{\mathbf{n}}$ has GLT
symbol $\mathbf{f}(\boldsymbol{\theta,x})$ writing
$\\{A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\textsc{glt}}\mathbf{f}(\boldsymbol{\theta,x})$,
$(\boldsymbol{\theta,x})\in[-\pi,\pi]^{d}\times[0,1]^{d}$.

GLT2:

The set of GLT sequences form a $*$-algebra, i.e., it is closed under linear
combinations, products, inversion (whenever the symbol is singular, at most,
in a set of zero Lebesgue measure), and conjugation. Hence, we obtain the GLT
symbol of algebraic operations of a finite set of GLT sequences by performing
the same algebraic manipulations of the symbols of the considered GLT
sequences.

GLT3:

Every Toeplitz sequence $\\{T_{\mathbf{n}}(\mathbf{f})\\}_{n}$ generated by a
function $\mathbf{f}(\boldsymbol{\theta})$ belonging to
$L^{1}([-\pi,\pi]^{d})$ is a GLT sequence and with GLT symbol given by
$\mathbf{f}$. Every diagonal sampling sequence
$\\{D_{\mathbf{n}}(\mathbf{a})\\}_{n}$ generated by a Riemann integrable
function $\mathbf{a}(\boldsymbol{x})$, $\boldsymbol{x}\in[0,1]^{d}$ is a GLT
sequence and with GLT symbol given by $\mathbf{a}$.

GLT4:

Every sequence which is distributed as the constant zero in the singular value
sense is a GLT sequence with symbol $0$. In particular:

  * •:

every sequence in which the rank divided by the size tends to zero, as the
matrix-size tends to infinity;

  * •:

every sequence in which the trace-norm (i.e., sum of the singular values)
divided by the size tends to zero, as the matrix-size tends to infinity.

From a practical view-point, on the one hand, one of the main advantages for a
sequence of belonging to the GLT class is that, under certain hypotheses,
crucial spectral and singular value information can be derived using the
concept of GLT symbol. On the other hand, the above properties imply the
following important features of the GLT symbol. Given a sequence
$\\{A_{\mathbf{n}}\\}_{\mathbf{n}}$ obtained by algebraic operations of a
finite set of GLT sequences, the small-norm and low-rank terms which composes
the sequence should be neglected in the computation of the GLT symbol.
Consequently, it happens that for small matrix-sizes $n$, the approximations
may not be as accurate as it is desirable.

For this reason in [6] it has been introduced and exploited the concept of a
(singular value and spectral) “momentary symbols”, starting from a special
case of Toeplitz structures. Here we generalize the notion to that of “GLT
momentary symbols”: the construction stems from that of the symbol in the GLT
sense, but in practice the information of the small norm contributions is kept
in the symbol and this may lead to higher accuracy, at least in some
emblematic cases, when approximating the singular values and eigenvalues of
Toeplitz-like matrices, even for small dimensions.

###  1.6. The GLT momentary symbol sequence

For clarity in this subsection we consider the matrix-sequences in detail only
in the unilevel and scalar setting. We want to avoid a cumbersome notation,
but the ideas are extensible in a plain manner to the case where the involved
GLT symbols are also matrix-valued and multivariate, as briefly sketched.

As an example, we take the following second-order differential equation with
Dirichlet boundary conditions

$\begin{cases}-(a(x)u^{\prime}(x))^{\prime}+b(x)u^{\prime}(x)+c(x)u(x)=f(x),&x\in(0,1),\\\
u(0)=\alpha,\qquad u(1)=\beta.&\end{cases}$

The well-posedness of the previous diffusion-convection-advection problem
holds in the case where $a(x)\in C^{1}(0,1)$. Furthermore, the uniqueness and
existence of the solution are guaranteed in the case where $a(x)>0$, $c(x)\geq
0$ and with continuous functions $b(x),c(x)$ on $[0,1]$, with $f(x)\in
L^{2}([0,1])$ (see [8]). For a more exhaustive discussion regarding the
conditions of existence and uniqueness, even in the multidimensional case, we
refer to [27, 28] and references therein.

From a GLT viewpoint, we only require the following much weaker assumptions
that are

  * •

$a(x),c(x)$ are real-valued functions, continuous almost everywhere, defined
in $[0,1]$,

  * •

$b(x)$ is a real-valued function on $[0,1]$, such that $|b(x)x^{\alpha}|$ is
bounded for some $\alpha<3/2$,

while $f(x)$ is a general function.

We employ central second-order finite differences for approximating the given
equation. We define the stepsize $h=\frac{1}{n+1}$ and the points $x_{k}=kh$
for $k$ belonging to the interval $[0,n+1]$. Let $a_{k}:=a(x_{\frac{k}{2}})$
for any $k\in[0,2n+2]$ and set $b_{j}:=b(x_{j})$, $c_{j}:=c(x_{j})$,
$f_{j}:=f(x_{j})$ for every $j=0,\ldots,n+1$. We compute approximations
$u_{j}$ of the values $u(x_{j})$ for $j=1,\ldots,n$ by solving the following
linear system

$A_{n}\begin{pmatrix}u_{1}\\\ u_{2}\\\ \vdots\\\ u_{n-1}\\\
u_{n}\end{pmatrix}+B_{n}\begin{pmatrix}u_{1}\\\ u_{2}\\\ \vdots\\\ u_{n-1}\\\
u_{n}\end{pmatrix}+C_{n}\begin{pmatrix}u_{1}\\\ u_{2}\\\ \vdots\\\ u_{n-1}\\\
u_{n}\end{pmatrix}=h^{2}\begin{pmatrix}f_{1}+\frac{1}{h^{2}}a_{1}\alpha+\frac{1}{2h}b_{1}\alpha\\\
f_{2}\\\ \vdots\\\ f_{n-1}\\\
f_{n}+\frac{1}{h^{2}}a_{2n+1}\beta-\frac{1}{2h}b_{n}\beta\end{pmatrix},$
(1.12)

where

$A_{n}=\begin{pmatrix}a_{1}+a_{3}&-a_{3}&&&\\\ -a_{3}&a_{3}+a_{5}&-a_{5}&&\\\
&\ddots&\ddots&\ddots&\\\ &&-a_{2n-3}&a_{2n-3}+a_{2n-1}&-a_{2n-1}\\\
&&&-a_{2n-1}&a_{2n-1}+a_{2n+1}\end{pmatrix},$
$B_{n}=\frac{h}{2}\begin{pmatrix}0&b_{1}&&&\\\ -b_{2}&0&b_{2}&&\\\
&\ddots&\ddots&\ddots&\\\ &&-b_{n-1}&0&b_{n-1}\\\
&&&-b_{n}&0\end{pmatrix},\quad
C_{n}=h^{2}\operatorname*{diag}(c_{1},\ldots,c_{n}).$

In the case where $a(x)\equiv 1$ and $b(x)\equiv 1$, we find the basic
Toeplitz structures

$K_{n}=T_{n}(2-2\cos\theta)=\begin{pmatrix}2&-1&&&\\\ -1&2&-1&&\\\
&\ddots&\ddots&\ddots&\\\ &&-1&2&-1\\\ &&&-1&2\end{pmatrix},$
$H_{n}=T_{n}(\mathbf{i}\sin\theta)=\frac{1}{2}\begin{pmatrix}0&1&&&\\\
-1&0&1&&\\\ &\ddots&\ddots&\ddots&\\\ &&-1&0&1\\\ &&&-1&0\end{pmatrix},$

which are of importance since $A_{n}=D_{n}(a)K_{n}+E_{n}$ and
$B_{n}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h}D_{n}(b)H_{n}$
with $\\{E_{n}\\}_{n}\sim_{\sigma}0$, as an immediate check in [4] can show.
Therefore, by using the GLT axioms (as done in detail in [4]) we obtain

$\\{A_{n}\\}_{n}\sim_{\textsc{glt}}a(x)(2-2\cos\theta),\quad\left\\{\frac{1}{h}B_{n}\right\\}_{n}\sim_{\textsc{glt}}\mathbf{i}b(x)\sin\theta,\quad\left\\{\frac{1}{h^{2}}C_{n}\right\\}_{n}\sim_{\textsc{glt}}c(x).$

As a conclusion $\left\\{B_{n}\right\\}_{n}\sim_{\textsc{glt}}\text{0}$, $\
\left\\{C_{n}\right\\}_{n}\sim_{\textsc{glt}}0$ and hence, setting
$X_{n}=A_{n}+B_{n}+C_{n}$ the actual coefficient matrix of the linear system
in (1.12), again by the $*$-algebra structure of the GLT matrix-sequences, we
deduce

$\\{X_{n}\\}_{n}\sim_{\textsc{glt}}a(x)(2-2\cos\theta).$

Now, following [6], the idea is to consider not only the asymptotic setting,
but also the case of moderate sizes. As a consequence, for increasing the
precision of the evaluation of eigenvalues and singular values, we can
associate to

$X_{n}=A_{n}+B_{n}+C_{n}$

the specific symbol
$f_{n}(x,\theta)=a(x)(2-2\cos\theta)+h\mathbf{i}b(x)\sin\theta+h^{2}c(x)$.

We are now in position to give a formal definition of GLT momentary symbols.

######  Definition 1.11 (GLT momentary symbols).

Let $\\{X_{n}\\}_{n}$ be a matrix-sequence and assume that there exist matrix-
sequences $\\{A_{n}^{(j)}\\}_{n}$, scalar sequences $c_{n}^{(j)}$,
$j=0,\ldots,t$, and measurable functions $f_{j}$ defined over
$[-\pi,\pi]\times[0,1]$, $t$ nonnegative integer independent of $n$, such that

$\displaystyle\left\\{\frac{A_{n}^{(j)}}{c_{n}^{(j)}}\right\\}_{n}$
$\displaystyle\sim_{\textsc{glt}}$ $\displaystyle f_{j},$ $\displaystyle
c_{n}^{(0)}=1,$ $\displaystyle c_{n}^{(s)}=o(c_{n}^{(r)}),\ \ t\geq s>r,$
$\displaystyle\\{X_{n}\\}_{n}$ $\displaystyle=$
$\displaystyle\\{A_{n}^{(0)}\\}_{n}+\sum_{j=1}^{t}\\{A_{n}^{(j)}\\}_{n}.$
(1.13)

Then, by a slight abuse of notation,

$f_{n}=f_{0}+\sum_{j=1}^{t}c_{n}^{(j)}f_{j}$ (1.14)

is defined as the GLT momentary symbol for $X_{n}$ and $\\{f_{n}\\}$ is the
sequence of GLT momentary symbols for the matrix-sequence $\\{X_{n}\\}_{n}$.

Of course, in line with Section 1.5, the momentary symbol could be matrix-
valued with a number of variables equal to $2d$ and domain
$[-\pi,\pi]^{d}\times[0,1]^{d}$ if the basic matrix-sequences appearing in
Definition 1.11 are, up to proper scaling, matrix-valued and multilevel GLT
matrix-sequences. For example in the scalar $d$-variate setting relation
(1.14) takes the form

$f_{\textbf{n}}=\sum_{\textbf{j}=\textbf{0}}^{\textbf{t}}c_{\textbf{n}}^{(\textbf{j})}f_{\textbf{j}},$

which is a plain multivariate (possibly block) version of (1.14).

Clearly there is a link with the GLT theory stated in the next result.

######  Theorem 1.12.

Assume that the matrix-sequence $\\{X_{n}\\}_{n}$ satisfies the requirements
in Definition 1.11. Then $\\{X_{n}\\}_{n}$ is a GLT matrix sequence and the
GLT symbol $f_{0}$ of the main term $A_{n}^{(0)}$ is the GLT symbol of
$\\{X_{n}\\}_{n}$, that is, $\\{X_{n}\\}_{n}\sim_{\textsc{glt}}f_{0}$ and
$\lim_{n\to\infty}f_{n}=f_{0}$ uniformly on the definition domain.

The given definition of momentary symbols is inspired, as it is clear from the
initial example of diffusion-convection-advection equation, by the example of
approximated differential equations, where the presence of differential
operators of different orders induces, after a possible proper scaling, a
structure like that reported in (1.13).

The idea is that the momentary symbol can be used for giving a more precise
evaluation either of the spectrum or of the eigenvalues for moderate sizes of
the matrices and not only asymptotically. However, we should be aware that,
intrinsically, there is no general recipe especially for the eigenvalues. In
fact, as already proven in [26], a rank one perturbation of infinitesimal
spectral norm actually can change the spectra of matrix-sequences, sharing the
same GLT symbol and even sharing the same sequence of momentary symbols.

Example 1:

Take the matrices $T_{n}(e^{\mathbf{i}\theta})$ and
$X_{n}=T_{n}(e^{\mathbf{i}\theta})+e_{1}e_{n}^{T}c_{n}^{(1)}$ with
$c_{n}^{(1)}=n^{-\alpha}$, $\alpha>0$ any positive number independent of the
matrix-size $n$. By direct inspection
$\\{e_{1}e_{n}^{T}c_{n}^{(1)}\\}_{n}\sim_{\sigma}0$ and hence it is a GLT
matrix-sequence with zero symbol, independently of the parameter $\alpha$. If
we look at the GLT momentary symbols then they coincide with the GLT symbol
for both $\\{T_{n}(e^{\mathbf{i}\theta})\\}_{n}$ and $\\{X_{n}\\}_{n}$:
however while in the first case, the eigenvalues are all equal to zero, in the
second case they distribute asymptotically as the GLT symbol
$e^{\mathbf{i}\theta}$ (which is also the GLT momentary symbol for any $n$).

Example 2:

Take a positive function $a$ defined on $[0,1]$ and the matrices
$D_{n}(a)T_{n}(e^{\mathbf{i}\theta})$ and
$X_{n}=D_{n}(a)T_{n}(e^{\mathbf{i}\theta})+e_{1}e_{n}^{T}c_{n}^{(1)}$ with
$c_{n}^{(1)}=n^{-\alpha}$, $\alpha>0$ any positive number independent of the
matrix-size $n$. Since $\\{e_{1}e_{n}^{T}c_{n}^{(1)}\\}_{n}$ is a GLT matrix-
sequence with zero symbol, independently of the parameter $\alpha$, we deduce
that both $\\{D_{n}(a)T_{n}(e^{\mathbf{i}\theta})\\}_{n}$ and
$\\{X_{n}\\}_{n}$ share the same GLT symbol $a(x)e^{\mathbf{i}\theta}$ (which
is also the momentary symbol for any $n$). Again there is dramatic change:
while in the first case, the eigenvalues are all equal to zero, in the second
case they distribute asymptotically as the function
$\hat{a}e^{\mathbf{i}\theta}$, where $\hat{a}$ is the limit (if it exists) of
the geometric mean of sampling values present in $D_{n}(a)$, as $n$ tends to
infinity: since $n^{-\alpha/n}$ converges to $1$ independently of the
parameter $\alpha$ as $n$ tends to infinity, $\hat{a}$ will depend only on the
diagonal values of $D_{n}(a)$. As a conclusion the eigenvalue distributions do
not coincide with the GLT momentary symbols and this is a message that the
present tool could be not effective and even misleading, when very non-normal
matrices are considered.

In this setting it must be emphasized that the asymptotic eigenvalue
distribution is discontinuous with respect to the standard norms or metrics
widely considered in the context of matrix-sequences.

##  2\. All-at-once solution of parabolic problems

The aim of this section is that of describing as accurate as possible the
spectra and singular values of the structured linear system sequence stemming
by the space-time discretization for a parabolic diffusion problem. Then, we
consider the diffusion equation in one space dimension,

$u_{t}=u_{xx},\quad x\in(a,b),\ t\in[0,T],$

where we are prescribing $u$ at $t=0$ and imposing the periodicity condition
$u(x\pm(b-a),t)=u(x,t)$.

We approximate our parabolic model problem on a rectangular space-time grid
consisting of $N_{t}$ time intervals and $N_{x}$ space intervals. We obtain a
sequence of linear systems, in which the each component is of the form

$A_{\mathbf{n}}x=b,\quad A_{\mathbf{n}}=J_{N_{t}}\oplus
Q_{N_{x}}=J_{N_{t}}\otimes\mathbb{I}_{N_{x}}+\mathbb{I}_{N_{t}}\otimes
Q_{N_{x}}\in\mathbb{R}^{N\times N},\quad x,b\in\mathbb{R}^{N},$ (2.1)

where $N=N_{t}N_{x}$, $\mathbf{n}=(N_{t},N_{x}),$ $\mathbb{I}_{m}$ is the
identity matrix of size $m$, and the matrices $J_{N_{t}}$ and $Q_{N_{x}}$ come
from the discretization in time and space, respectively. In the following, we
describe the time and space discretization and, in particular, how this leads
to structured components of the matrix $A_{\mathbf{n}}$.

###  2.1. Time discretization

The principal ingredients of the time discretization are:

  * •

Choosing $N_{t}$ equispaced points in $[0,T]$ with stepsize $h_{t}=T/N_{t}$,
that is, $t_{j}=jh_{t}$, for $j=1,\ldots,N_{t}$.

  * •

Discretizing in time by standard Euler backwards.

Regarding notations, for the sake of simplicity, since we are considering a 2D
problem, the symbols will have as Fourier variable $(\theta,\xi)$ instead of
the standard choice $(\theta_{1},\theta_{2})$ indicated in the notations of
Section 1.2 (see Definition 1.4).

The resulting matrix is $J_{N_{t}}$, which has the following unilevel scalar
Toeplitz structure:

$J_{N_{t}}=\frac{1}{h_{t}}\begin{bmatrix}1\\\ -1&1\\\ &\ddots&\ddots\\\
&&-1&1\end{bmatrix}=\frac{1}{h_{t}}T_{N_{t}}(f_{J}),$ (2.2)

where $f_{J}$ is the generating function of the matrix-sequence
$\\{T_{n}(f_{J})\\}_{n}$ with

$f_{J}(\theta)=1-\mathrm{e}^{\mathbf{i}\theta}.$

###  2.2. Space discretization

The principal elements of the time discretization are:

  * •

Choosing $N_{x}$ equispaced points in $[a,b]$. Since we are considering
periodic boundary conditions, we have step size $h_{x}=(a-b)/N_{x}$ and
$x_{j}=h_{x}(j-1)$, for $j=1,\ldots,N_{x}$.

  * •

Discretizing in space using second order finite differences.

Consequently, the space discretization matrix will be the circulant matrix
$Q_{N_{x}}$ of the form:

$Q_{N_{x}}=\frac{1}{h_{x}^{2}}\begin{bmatrix}2&-1&&&-1\\\ -1&2&-1\\\
&\ddots&\ddots&\ddots\\\ &&-1&2&-1\\\
-1&&&-1&2\end{bmatrix}=\frac{1}{h_{x}^{2}}C_{N_{x}}(f_{Q}),$

where

$f_{Q}(\xi)=2-2\cos\xi$

is the generating function of the matrix.

Of course a different choice of Dirichlet boundary conditions would lead to
the standard discrete Laplacian $T_{N_{x}}(f_{Q})$: the analysis is equivalent
since also this matrix admits a well known diagonalization matrix, that is the
sine transform matrix of type I, which is real, orthogonal and symmetric.

###  2.3. Analysis of the coefficient matrix $A_{\mathbf{n}}$

We have seen that discretizing of the problem of interest for a sequence of
discretization parameters $h_{x}$ and $h_{t}$ leads to a sequence of linear
systems, whose approximation error tends to zero as the coefficient matrix-
size grows to infinity. The $\mathbf{n}$th coefficient matrix component is of
the form

$A_{\mathbf{n}}=\frac{1}{h_{t}}T_{N_{t}}(f_{J})\otimes\mathbb{I}_{N_{x}}+\mathbb{I}_{N_{t}}\otimes\frac{1}{h_{x}^{2}}C_{N_{x}}(f_{Q}).$
(2.3)

In order to design efficient solvers for the considered linear systems, it is
of crucial importance to know the spectral propriety of the matrix-sequence
$\\{A_{\mathbf{n}}\\}_{\mathbf{n}}$. Hence, this section is devoted to the
analysis of the structure of the matrix-sequence
$\\{{A}_{\mathbf{n}}\\}_{\mathbf{n}}$ in (2.1). In particular, we provide the
singular values and spectral analysis using algebraic tricks, the GLT theory,
and the concept of GLT momentary symbols.

###  2.4. GLT analysis of the coefficient sequence
$\\{{A}_{\mathbf{n}}\\}_{\mathbf{n}}$

The asymptotic spectral and singular value distribution, for the matrix-size
$d(\mathbf{n})$ sufficiently large, of the matrix-sequence
$\\{{A}_{\mathbf{n}}\\}_{\mathbf{n}}$ depend on how $h_{x}$ and $h_{t}$
approaches zero. Let $c_{h}\coloneqq h_{x}^{2}/h_{t}$, we have three different
cases to consider.

  * Case 1.

$\left[c_{h}\to\infty\right]:$ If $h_{t}\to 0$ faster than $C_{1}h_{x}^{2}$,
where $C_{1}$ is a constant, then we can consider the matrix

$h_{t}{A}_{\mathbf{n}}=T_{N_{t}}(f_{J})\otimes\mathbb{I}_{N_{x}}+\mathbb{I}_{N_{t}}\otimes\underbrace{\frac{h_{t}}{h_{x}^{2}}}_{c_{h}^{-1}\to
0}C_{N_{x}}(f_{Q}).$

Then, the sequence
$\\{h_{t}A_{\mathbf{n}}\\}_{\mathbf{n}}=\\{T_{N_{t}}(f_{J})\otimes\mathbb{I}_{N_{x}}+\mathbb{N}_{\mathbf{n}}\\}_{\mathbf{n}}$,
where $\mathbb{N}_{\mathbf{n}}$ is a small-norm matrix in the sense of the
item 2 of property GLT4, with $\|\mathbb{N}_{\mathbf{n}}\|<C_{2}$, $C_{2}$
constant. Consequently, from GLT4,
$\\{\mathbb{N}_{\mathbf{n}}\\}_{\mathbf{n}}$ is a matrix-sequence distributed
in the singular value sense as $0$, which implies that
$\\{\mathbb{N}_{\mathbf{n}}\\}_{\mathbf{n}}$ is zero-distribued in GLT sense
as described in GLT4. Moreover $f_{J}$ is a trigonometric polynomial, then
Theorem 1.10, properties GLT1-GLT4 and Lemma 1.6 imply that

$\\{h_{t}A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\textsc{glt}}f_{J}(\theta)\otimes
1+1\otimes 0=f_{A}^{(1)}(\theta,\xi).$

The $1$ present in $f_{J}(\theta)\otimes 1$ should be interpreted as
$1\mathrm{e}^{0\mathbf{i}\xi}$, and $1\otimes 0$ should be interpreted as
$1\mathrm{e}^{0\mathbf{i}\theta}\otimes 0\mathrm{e}^{0\mathbf{i}\xi}$. Hence,
the GLT symbol of the sequence $\\{h_{t}A_{\mathbf{n}}\\}_{\mathbf{n}}$ is the
bivariate function

$f_{A}^{(1)}(\theta,\xi)=f_{J}(\theta)=1-\mathrm{e}^{\mathbf{i}\theta},$

and it should be interpreted as the function $f_{J}(\theta)\otimes 1$, with is
constant in the second component.

From the property GLT1, the function $f_{A}^{(1)}(\theta,\xi)$ describes the
singular value distribution in the sense of relation (1.4). More in detail

$\\{h_{t}A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\textsc{glt},\sigma}1-\mathrm{e}^{\mathbf{i}\theta}.$

However, the matrix-sequence $\\{h_{t}A_{\mathbf{n}}\\}_{\mathbf{n}}$ is not
symmetric, hence the distribution does not hold in the eigenvalue sense (see
also Example 1 and Example 2 at the end of Section 1.6). Because of the
structure of $J_{N_{t}}=\frac{1}{h_{t}}T_{N_{t}}(f_{J})$ in equation (2.2) it
is straightforward to see that the asymptotic spectral distribution is given
by ${\mathfrak{f}({\theta},\xi)}=1,$ accordingly to relation (1.5), that is

$\\{h_{t}A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\lambda}1.$

  * Case 2.

$\left[c_{h}\to 0\right]:$ If $h_{x}^{2}\to 0$ faster than $C_{1}h_{t}$, where
$C_{1}$ is a constant, then we have

$h_{x}^{2}A_{\mathbf{n}}=\underbrace{\frac{h_{x}^{2}}{h_{t}}}_{c_{h}\to
0}T_{N_{t}}(f_{J})\otimes\mathbb{I}_{N_{x}}+\mathbb{I}_{N_{t}}\otimes
C_{N_{x}}(f_{Q}).$

Then, the sequence
$\\{h_{x}^{2}A_{\mathbf{n}}\\}_{\mathbf{n}}=\\{\mathbb{N}_{\mathbf{n}}+\mathbb{I}_{N_{t}}\otimes
C_{N_{x}}(f_{Q})\\}_{\mathbf{n}}$, where $\mathbb{N}_{\mathbf{n}}$ is a small-
norm matrix in the sense of the item 2 of property GLT4, with
$\|\mathbb{N}_{\mathbf{n}}\|<C_{2}$, $C_{2}$ constant. Then,
$\\{\mathbb{N}_{\mathbf{n}}\\}_{\mathbf{n}}$ is a matrix-sequence distributed
in the singular value sense, and consequently in the GLT sense, as $0$.
Moreover, ${f_{Q}}$ belongs to the Dini-Lipschitz class, consequently,
properties GLT2-GLT4, and Lemma 1.6 imply that

$\\{h_{x}^{2}A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\textsc{glt}}0\otimes
1+1\otimes f_{Q}(\xi)=1\otimes f_{Q}(\xi)=f_{A}^{(2)}(\theta,\xi),$

where the GLT symbol is given by

$f_{A}^{(2)}(\theta,\xi)=f_{Q}(\xi)=2-2\cos\xi.$

In this case the function $f_{A}^{(2)}(\theta,\xi)$ is a singular value symbol
for the sequence $\\{h_{x}^{2}A_{\mathbf{n}}\\}_{\mathbf{n}}$, and also an
eigenvalue symbol, since the matrices $C_{N_{x}}(f_{Q})$ are Hermitian for
each $N_{x}$. Hence we have

$\\{h_{x}^{2}A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\textsc{glt},\sigma,\lambda}2-2\cos\xi.$

  * Case 3.

$\left[c_{h}=c=\text{constant}\right]:$ The last case is when $h_{x}^{2}$ and
$h_{t}$ are proportional and related by the constant
$c_{h}=c=\frac{h_{x}^{2}}{h_{t}}$, independent of the various step-sizes. In
this setting we have

$h_{x}^{2}A_{\mathbf{n}}=\underbrace{\frac{h_{x}^{2}}{h_{t}}}_{c_{h}}T_{N_{t}}(f_{J})\otimes\mathbb{I}_{N_{x}}+\mathbb{I}_{N_{t}}\otimes
C_{N_{x}}(f_{Q}).$

Consequently, from GLT2, GLT3 and Lemma 1.6, the following relationship holds
when $c_{h}$ is a constant,

$\\{h_{x}^{2}A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\textsc{glt}}cf_{J}(\theta)\otimes
1+1\otimes f_{Q}(\xi)=f_{A}^{(3)}(\theta,\xi).$

From considerations analogous to the case 1 and 2 we have

$\\{h_{x}^{2}A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\textsc{glt},\sigma}c(1-\mathrm{e}^{\mathbf{i}\theta})+(2-2\cos\xi).$

Since the matrix $h_{x}^{2}A_{\mathbf{n}}$ is not Hermitian, the eigenvalue
symbol $\mathfrak{f}({\theta},\xi)$ cannot be directly derived by
$f_{A}^{(3)}(\theta,\xi)$ (see again the discussion in the examples after
Definition 1.11).

In this setting the situation is simple because the involved twolevel
structure can be simply block-diagonalized, while the use of the GLT momentary
symbol becomes useful in approximation the singular values of the sequence
$\\{h_{x}^{2}A_{\mathbf{n}}\\}_{\mathbf{n}}$.

###  2.5. Analysis of the coefficient matrix-sequence
$\\{{A}_{\mathbf{n}}\\}_{\mathbf{n}}$ by algebraic manipulations and GLT
momentary symbols

The first observation is that the matrix in (2.3) admits a perfect
decomposition which shows in evidence a lower triangular matrix, which is
similar to the original one and hence all the eigenvalues are known exactly.
In fact, by looking carefully at (2.3), we obtain that

$\frac{1}{h_{t}}T_{N_{t}}(f_{J})\otimes\mathbb{I}_{N_{x}}=\mathbb{I}_{N_{t}}\left[\frac{1}{h_{t}}T_{N_{t}}(f_{J})\right]\mathbb{I}_{N_{t}}\otimes\mathbb{F}_{N_{x}}\mathbb{I}_{N_{x}}\mathbb{F^{*}}_{N_{x}}$

and

$\mathbb{I}_{N_{t}}\otimes\frac{1}{h_{x}^{2}}C_{N_{x}}(f_{Q})=\mathbb{I}_{N_{t}}\mathbb{I}_{N_{t}}\mathbb{I}_{N_{t}}\otimes\mathbb{F}_{N_{x}}\frac{1}{h_{x}^{2}}D_{N_{x}}\mathbb{F^{*}}_{N_{x}},$

where $\mathbb{F}_{N_{x}}$ is the unitary Fourier matrix of size $N_{x}$,
$\mathbb{F^{*}}_{N_{x}}$ is its transpose conjugate and hence its inverse, and
$D_{N_{x}}$ is the diagonal matrix containing the eigenvalues of
$C_{N_{x}}(f_{Q})$ that is $f_{Q}(2\pi j/N_{x})=2-2\cos(2\pi j/N_{x})$,
$j=0,1,\ldots,N_{x}-1$.

Since $T_{N_{t}}(f_{J})$ is lower bidiagonal matrix with $1$ on the main
diagonal, it can be easily seen that the eigenvalues of $A_{\mathbf{n}}$ in
(2.3) are exactly

$\frac{1}{h_{t}}+\frac{1}{h_{x}^{2}}(2-2\cos(2\pi j/N_{x})),\ \ \
j=0,1,\ldots,N_{x}-1,$

each of them with multiplicity $N_{t}$. As a consequence, by taking a proper
normalization, the spectral radius $\rho(h_{x}^{2}A_{\mathbf{n}})$ will
coincide simply with $4+c_{h}$.

It is clear that, in this context, due to the high non-normality of the term
$T_{N_{t}}(f_{J})$, after proper scalings depending on $h_{t}$ and $h_{x}$,
the eigenvalues are a uniform sampling of a function which is not the GLT
symbol and is not the associated GLT momentary symbol. This is not surprising
given the discussion regarding the asymptotical behaviour of the matrix-
sequences reported in Example 1 and in Example 2, when discussing the
potential and the limitations of the notion of GLT momentary symbols.

Also in this setting, by imposing (quite artificial) periodic boundary
conditions in time, the term $T_{N_{t}}(f_{J})$ will change into
$C_{N_{t}}(f_{J})$ and magically a one-rank correction repeated $N_{t}$ times
to the matrix $A_{\mathbf{n}}$ will produce a new matrix with the same GLT and
momentary symbols as before: however in this case the eigenvalues will be
exactly the sampling of such functions. This is a further confirmation of the
delicacies of the eigenvalues that can have dramatic changes due to minimal
corrections, when we are in a context of higlhy non-normal matrices.

####  2.5.1. Singular values of $h_{x}^{2}A_{\mathbf{n}}$ (exact)

The singular values $\sigma_{1}(h_{x}^{2}A_{\mathbf{n}}),\dots,$
$\sigma_{d(\mathbf{n})}(h_{x}^{2}A_{\mathbf{n}})$ of the matrix
$h_{x}^{2}A_{\mathbf{n}}$ are given the square root of the eigenvalues of the
Hermitian matrix $h_{x}^{4}A_{\mathbf{n}}A_{\mathbf{n}}^{\textsc{t}}$. Hence,
in order to provide exactly $\sigma_{i}(h_{x}^{2}A_{\mathbf{n}})$,
$i=1\dots,d(\mathbf{n})$, we are interested at the spectrum of the matrix

$\begin{split}&h_{x}^{4}A_{\mathbf{n}}A_{\mathbf{n}}^{\textsc{t}}=\\\
&\left[\begin{array}[]{cccccccccc}\tilde{Q}_{N_{x}}^{2}&-c_{h}\tilde{Q}_{N_{x}}\\\
-c_{h}\tilde{Q}_{N_{x}}&\tilde{Q}_{N_{x}}^{2}+c_{h}^{2}\mathbb{I}_{N_{x}}&-c_{h}\tilde{Q}_{N_{x}}\\\
&-c_{h}\tilde{Q}_{N_{x}}&\tilde{Q}_{N_{x}}^{2}+c_{h}^{2}\mathbb{I}_{N_{x}}&\ddots&\\\
&&&&\\\ &&\ddots&\ddots&-c_{h}\tilde{Q}_{N_{x}}\\\
&&&-c_{h}\tilde{Q}_{N_{x}}&\tilde{Q}_{N_{x}}^{2}+c_{h}^{2}\mathbb{I}_{N_{x}}\end{array}\right],\end{split}$

where $\tilde{Q}_{N_{x}}=C_{N_{x}}+c_{h}\mathbb{I}_{N_{x}}$. Note that
$h_{x}^{4}A_{\mathbf{n}}A_{\mathbf{n}}^{\textsc{t}}$ is not a pure block-
tridiagonal Toeplitz because of the missing constant $c_{h}^{2}$ in the block
in the top left corner. However, for each fixed $N_{t}$ and $N_{x}$, the
matrix $\tilde{Q}_{N_{x}}$ is a circulant matrix with generating function
$f_{\tilde{Q}_{N_{x}}}(\xi)=2-2\cos\xi+c_{h}$, which is also its GLT momentary
symbol. Thus we infer that
$h_{x}^{4}A_{\mathbf{n}}A_{\mathbf{n}}^{\textsc{t}}$ is similar to a matrix
$X_{d(\mathbf{n})}$, whose explicit expression is reported below

$\begin{split}&h_{x}^{4}A_{\mathbf{n}}A_{\mathbf{n}}^{\textsc{t}}\sim
X_{d(\mathbf{n})}=\\\
&\left[\begin{array}[]{cccccccccc}D_{\tilde{Q}}^{2}&-c_{h}D_{\tilde{Q}}\\\
-c_{h}D_{\tilde{Q}}&D_{\tilde{Q}}^{2}+c_{h}^{2}\mathbb{I}_{N_{x}}&-c_{h}D_{\tilde{Q}}\\\
&-c_{h}D_{\tilde{Q}}&D_{\tilde{Q}}^{2}+c_{h}^{2}\mathbb{I}_{N_{x}}&\ddots\\\
&&\ddots&\ddots&-c_{h}D_{\tilde{Q}}\\\
&&&-c_{h}D_{\tilde{Q}}&D_{\tilde{Q}}^{2}+c_{h}^{2}\mathbb{I}_{N_{x}}\end{array}\right],\end{split}$

with
$D_{\tilde{Q}}=\operatorname*{diag}_{\ell=1,\dots,N_{x}}\left(f_{\tilde{Q}_{N_{x}}}(\xi_{\ell,N_{x}})\right)$.
Consequently we study the spectrum of $X_{d(\mathbf{n})}$ to attain formulas
for the exact singular values of $h_{x}^{2}A_{\mathbf{n}}$. Let us consider a
permutation matrix $P$ such transforms $X_{d(\mathbf{n})}$ into an
$N_{t}\times N_{t}$ block diagonal matrix $PX_{d(\mathbf{n})}P^{\textsc{t}}$,
which has on the main diagonal, for $k=1,\ldots,N_{x}$, blocks of the form

$\begin{split}&\left(PX_{d(\mathbf{n})}P^{\textsc{t}}\right)_{i,j=(k-1)N_{t}+1}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}kN_{t}}}=\\\
&\left[\begin{array}[]{cccccccccc}C_{k}^{2}&-c_{h}C_{k}\\\
-c_{h}C_{k}&C_{k}^{2}+c_{h}^{2}&-c_{h}C_{k}\\\
&-c_{h}C_{k}&C_{k}^{2}+c_{h}^{2}&\ddots\\\ &&\ddots&\ddots&-c_{h}C_{k}\\\
&&&-c_{h}C_{k}&C_{k}^{2}+c_{h}^{2}\end{array}\right],\end{split}$ (2.4)

where $C_{k}=D_{\tilde{Q}}(k,k)=f_{\tilde{Q}_{N_{x}}}(\xi_{k,N_{x}})$.

Hence, the union of the eigenvalues of all blocks of
$PX_{d(\mathbf{n})}P^{\textsc{t}}$ is equivalent to the full spectrum of
$h_{x}^{4}A_{\mathbf{n}}A_{\mathbf{n}}^{\textsc{t}}$. These local eigenvalue
problems can be solved analytically (or numerically) independently from each
other. For example for $N_{t}=2$ we have for every $k=1,\ldots,N_{x}$ the
characteristic equation

$\left|\begin{array}[]{cccccccccc}\left(f_{\tilde{Q}_{N_{x}}}(\xi_{k,N_{x}})\right)^{2}-\lambda&-c_{h}(f_{\tilde{Q}_{N_{x}}}(\xi_{k,N_{x}}))\\\
-c_{h}(f_{\tilde{Q}_{N_{x}}}(\xi_{k,N_{x}}))&\left(f_{\tilde{Q}_{N_{x}}}(\xi_{k,N_{x}})\right)^{2}+c_{h}^{2}-\lambda\end{array}\right|=0.$

Thus, we have as singular value the union for $k=1,\dots,N_{x}$ of the
quantities

$\begin{split}\sigma^{(1)}(k,c_{h})&=\sqrt{\frac{2(f_{\tilde{Q}_{N_{x}}}(\xi_{k,N_{x}}))^{2}+c_{h}^{2}}{2}-\frac{c_{h}}{2}\sqrt{4(f_{\tilde{Q}_{N_{x}}}(\xi_{k,N_{x}}))^{2}+c_{h}^{2}}},\\\
\sigma^{(1)}(k,c_{h})&=\sqrt{\frac{2(f_{\tilde{Q}_{N_{x}}}(\xi_{k,N_{x}}))^{2}+c_{h}^{2}}{2}+\frac{c_{h}}{2}\sqrt{4(f_{\tilde{Q}_{N_{x}}}(\xi_{k,N_{x}}))^{2}+c_{h}^{2}}}.\end{split}$

Clearly, solving the characteristic equation for $k=1,\ldots,N_{x}$ becomes
more and more complex as $N_{t}$ grows. Hence, in the next section provide two
possible approximations given by the GLT theory and by the GLT momentary
formulations.

####  2.5.2. Singular values of $h_{x}^{2}A_{\mathbf{n}}$ (approximation) via
GLT momentary symbols

For case (2), in Section 2.4 and in Section 2.5, we have already shown that

$\\{h_{x}^{2}A_{\mathbf{n}}\\}_{\mathbf{n}}\sim_{\sigma}f_{A}^{(2)}(\theta,\xi)=2-2\cos\xi.$

On the other hand, the subsequent sequence
$\\{f_{\mathbf{n}}^{(2)}\\}_{\textbf{n}}$ with

$f_{\mathbf{n}}^{(2)}(\theta,\xi)=c_{h}(1-\mathrm{e}^{\mathbf{i}\theta})+(2-2\cos\xi),$

is the sequence of GLT momentary functions.

Remark 1.8 suggests to exploit these relations in order to obtain a better
approximation of the singular value of $h_{x}^{2}A_{\mathbf{n}}$, with respect
to the information obtained by the pure GLT symbol. In the following, we
compute the quantities $|f_{A}^{(2)}(\theta,\xi)|$ and
$|f_{\mathbf{n}}^{(2)}(\theta,\xi)|$ using the specific grid described below

$\theta_{j,N_{t}}=\frac{j\pi}{N_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}+1}},\,j=1,\dots,N_{t}\qquad\xi_{\ell,N_{x}}=\frac{2\pi(\ell-1)}{N_{x}}\,\ell=1,\dots,N_{x}.$
(2.5)

We can observe in Figure 1 that the singular value of
$h_{x}^{2}A_{\mathbf{n}}$ (blue circles) are well approximated by the
samplings of $|f_{\mathbf{n}}^{(2)}(\theta,\xi)|$ on the grid (2.5) (red
stars). The approximation by using $|f_{A}^{(2)}(\theta,\xi)|$, instead, is
good when $c_{h}$ is small, see the top panel of Figure 1 for $N_{t}=2$ and
$N_{x}=10$, but it tends to become a substantially less accurate approximation
otherwise, see the bottom panel of Figure 1 and Figure 2 where $N_{t}=10$ and
$N_{x}=10$.

Figure 1. Singular values of $h_{x}^{2}A_{\mathbf{n}}$ and samplings of
$|f_{A}^{(2)}(\theta,\xi)|$ and $|f_{\mathbf{n}}^{(2)}(\theta,\xi)|$ on the
grid (2.5), for $N_{t}=2$ and $N_{x}=10$ (top) and $N_{t}=10$ and $N_{x}=10$
(bottom).  Figure 2. Singular values and samplings of
$|f_{A}^{(2)}(\theta,\xi)|$ and $|f_{\mathbf{n}}^{(2)}(\theta,\xi)|$ for
$N_{t}=10$ and $N_{x}=10$ on the grid in (2.5).

####  2.5.3. 2-norm of $h_{x}^{2}A_{\mathbf{n}}$ (approximation)

In the following we are interested in providing a bound for the $2$-norm of
the matrix $h_{x}^{2}A_{\mathbf{n}}$. By definition it is given by
$\|h_{x}^{2}A_{\mathbf{n}}\|_{2}=\max_{j=1\dots,d(\mathbf{n})}|\sigma_{j}(h_{x}^{2}A_{\mathbf{n}})|$.
From the previous section we know that it can be computed by making the square
root of the maximum eigenvalue of the block in (2.4), corresponding to
$\xi_{N_{x}/2+1,N_{x}}=\pi$. Since the
$\max_{k}\left(f_{\tilde{Q}_{N_{x}}}(\xi_{k,N_{x}})\right)=4+c_{h}$, we are
interested in estimate the maximum eigenvalue of

$\left[\begin{array}[]{cccccccccc}(4+c_{h})^{2}&-c_{h}(4+c_{h})\\\
-c_{h}(4+c_{h})&(4+c_{h})^{2}+c_{h}^{2}&-c_{h}(4+c_{h})\\\
&-c_{h}(4+c_{h})&(4+c_{h})^{2}+c_{h}^{2}&\ddots\\\
&&\ddots&\ddots&-c_{h}(4+c_{h})\\\
&&&-c_{h}(4+c_{h})&(4+c_{h})^{2}+c_{h}^{2}\end{array}\right].$ (2.6)

For this purpose we exploit the concept of
$\tau_{\varepsilon,\varphi}$-algebras of Subsection 1.4.

In our case

$a=(4+c_{h})^{2}+c_{h}^{2},\quad b=-c_{h}(4+c_{h}),$

and the matrix belongs to the $\tau_{\frac{c_{h}}{4+c_{h}},0}$-algebra, since
the element with indices $i,j=1$ is $a+(c_{h}/(4+c_{h}))b$.

Hence, we have $g(\theta)=(4+c_{h})^{2}+c_{h}^{2}-2c_{h}(4+c_{h})\cos\theta$,
which coincides with the eigenvalue symbol $\mathfrak{g}_{\mathbf{n}}$ of the
matrix (2.6). Due to the Interlacing theorem, see [5] and the specific
relation between algebras [6], the following relationships can be derived

$\displaystyle\underbrace{\mathfrak{g}_{\mathbf{n}}\left(\frac{\pi(N_{t}-1/2)}{N_{t}+1/2}\right)}_{\mathrm{max}(\lambda_{j}(T_{N_{t},1,0}(g)))}<\underbrace{\|h_{x}^{2}A_{\mathbf{n}}\|_{2}^{2}}_{\mathrm{max}(\lambda_{j}(T_{N_{t},c_{h}/(4+c_{h}),0}(g)))}<$
(2.7) $\displaystyle\underbrace{\mathfrak{g}_{\mathbf{n}}\left(\frac{\pi
N_{t}}{N_{t}+1}\right)}_{\mathrm{max}(\lambda_{j}(T_{N_{t},0,0}(g)))}<\underbrace{\mathfrak{g}_{\mathbf{n}}(\pi)}_{\mathrm{max}(\lambda_{j}(T_{N_{t},-1,-1}(f)))=\mathrm{max}(\mathfrak{g}_{\mathbf{n}})}.$

As a consequence, good upper and lower bounds for the 2-norm of
$h_{x}^{2}A_{\mathbf{n}}$ are reported in the following set of inequalities

$\sqrt{\mathfrak{g}_{\mathbf{n}}\left(\frac{\pi(N_{t}-1/2)}{N_{t}+1/2}\right)}<\|h_{x}^{2}A_{\mathbf{n}}\|_{2}<\sqrt{\mathfrak{g}_{\mathbf{n}}\left(\frac{\pi
N_{t}}{N_{t}+1}\right)}.$

In Table 3 we present approximations of the 2-norm of
$h_{x}^{2}A_{\mathbf{n}}$ using the grid sampling from $\tau_{1,0}$ (lower
bound), $\tau_{0,0}$ (upper bound), and $\tau_{-1,-1}$. Note that the sampling
on the latter grid is equivalent to do the sampling of the singular value
momentary symbols $f_{\mathbf{n}}^{(2)}(\theta,\xi)$ at their maximum point.
The two-norm $\|h_{x}^{2}A_{\mathbf{n}}\|_{2}$ is computed numerically. We see
that the 2-norm is well described by the two bounds given above, as $N_{t}$
increases. Hence, for this type of examples, the GLT momentary symbols
provide, at least for moderate sizes, a more precise alternative to the pure
GLT symbol.

Table 3. Approximations of $2$-norm for different $N_{t}$ and $c_{h}$. The maximum bound is $\sqrt{\max\mathfrak{g}_{\mathbf{n}}}=4+2c_{h}$. $N_{t}$ | $c_{h}$ | $\sqrt{\mathfrak{g}_{\mathbf{n}}\left(\frac{\pi(N_{t}-1/2)}{N_{t}+1/2}\right)}$ | $\|h_{x}^{2}A_{\mathbf{n}}\|_{2}$ | $\sqrt{\mathfrak{g}_{\mathbf{n}}\left(\frac{\pi N_{t}}{N_{t}+1}\right)}$ | $4+2c_{h}$
---|---|---|---|---|---
1 | 1/8 | 4.06394205 | 4.12500000 | 4.12689350 | 4.25
10 | 1/8 | 4.24460651 | 4.24505679 | 4.24508270 | 4.25
100 | 1/8 | 4.24994073 | 4.24994128 | 4.24994131 | 4.25
1000 | 1/8 | 4.24999940 | 4.24999940 | 4.24999940 | 4.25
1 | 1 | 4.58257569 | 5.00000000 | 5.09901951 | 6.00
10 | 1 | 5.96286240 | 5.96511172 | 5.96614865 | 6.00
100 | 1 | 5.99959287 | 5.99959555 | 5.99959689 | 6.00
1000 | 1 | 5.99999589 | 5.99999589 | 5.99999590 | 6.00
1 | 8 | 10.58300524 | 12.00000000 | 14.42220510 | 20.00
10 | 8 | 19.78560029 | 19.78964627 | 19.80461186 | 20.00
100 | 8 | 19.99765486 | 19.99765952 | 19.99767802 | 20.00
1000 | 8 | 19.99997634 | 19.99997634 | 19.99997636 | 20.00

##  3\. The case of approximations of distributed order differential operators
via asymptotic expansion and GLT momentary symbols

In this last section we focus on the matrix-sequences arising from the
numerical approximation of distributed-order operators. In detail, such a
procedure consists of two steps:

  1. (1)

Employ a quadrature formula to discretize the distributed-order operator into
a multi-term constant-order fractional derivative;

  2. (2)

discretize each constant-order fractional derivative.

In particular, we focus on the case where the matrices under consideration
take the form

$\frac{h^{\alpha_{\ell}}}{\Delta{\alpha}}\mathcal{T}_{n}=c_{\ell}T_{n}(g_{\alpha_{\ell}})+c_{\ell-1}h^{\Delta\alpha}T_{n}(g_{\alpha_{\ell-1}})+\dots+c_{1}h^{\Delta\alpha({\ell-1})}T_{n}(g_{\alpha_{1}}),$
(3.1)

where $\ell$ is a positive integer, all the coefficients $c_{j}$ are positive,
independent of $n$, and contained in a specific positive range
$[c_{*},c^{*}]$. Moreover, $\Delta\alpha=\frac{1}{\ell}$ and all the terms
$1<\alpha_{1}<\alpha_{2}\cdots<\alpha_{\ell}<2$ are positive and defined by
$a_{k}=1+\left(k-\frac{1}{2}\right)\Delta\alpha$, $k=1,\dots,\ell$. More
importantly all the functions $g_{\alpha_{\ell}}$ are globally continuous,
monotonically increasing in the interval $[0,\pi]$ and even in the whole
definition domain $[-\pi,\pi]$.

The goal is to exploit the notion of GLT momentary symbols and use it in
combination with the asymptotic expansions derived in a quite new research
line (see [16] and references there reported), in order to have a very precise
description of the spectrum of such matrices.

Indeed, under specific hypotheses on the generating function $f$, and fixing
an integer $\nu\geq 0$, it is possible to give an accurate description of the
eigenvalues of $T_{n}(f)$ via the following asymptotic expansion:

$\displaystyle\lambda_{j}(T_{n}(f))=w_{0}(\theta_{j,n})+hw_{1}(\theta_{j,n})+h^{2}w_{2}(\theta_{j,n})+\ldots+h^{\nu}w_{\nu}(\theta_{j,n})+E_{j,n,\nu},$

where the eigenvalues of $T_{n}(f)$ are arranged in ascending order,
$h=\frac{1}{n}$, and $\theta_{j,n}=\frac{j\pi}{n}=j\pi h$ for $j=1,\ldots,n$,
$E_{j,n,\nu}=O(h^{\nu+1})$ is the error. Moreover,
$\\{w_{k}\\}_{k=1,2,\ldots}$ is a sequence of functions from $[0,\pi]$ to
$\mathbb{R}$. The idea of such procedure is that a numerical approximation of
the value $w_{k}(\theta_{j,n})$ can be obtained by fast interpolation-
extrapolation algorithms (see [17] and references therein). In particular,
choosing $\nu$ proper grids $\theta_{j,n_{1}}$, $\theta_{j,n_{2}}$, $\dots$
$\theta_{j,n_{\nu}}$ with $n>>n_{\nu}>\dots>n_{1}$ an approximation of the
quantities $\tilde{w_{k}}({\theta}_{j,n})\approx w_{k}({\theta}_{j,n})$ can be
obtained. In the Hermitian case, we find that $\tilde{w}_{0}$ coincides with
the generating function.

Concerning the example in (3.1), the idea is to link the functions
$\tilde{w}^{c_{i}g_{\alpha_{i}}}_{k}$, $k=1,\dots,\nu$, associated with each
$c_{i}g_{\alpha_{i}}$ with the GLT momentary symbols of
$\frac{h^{\alpha_{\ell}}}{\Delta{\alpha}}\mathcal{T}_{n}$. Precisely, for
$j=1,\dots,n$, for a fixed $\nu$, we approximate the eigenvalues of
$\frac{h^{\alpha_{\ell}}}{\Delta{\alpha}}\mathcal{T}_{n}$ by

$\begin{split}&\lambda_{j}\left(\frac{h^{\alpha_{\ell}}}{\Delta{\alpha}}\mathcal{T}_{n}\right)\approx\\\
&c_{\ell}g_{\alpha_{\ell}}(\theta_{j,n})+\sum_{t=1}^{\nu}h^{t}\left(\tilde{w}_{t}^{{\alpha_{\ell}}}(\theta_{j,n})+\sum_{i=\ell-1}^{1}\tilde{w}_{t-1}^{{\alpha_{i}}}(\theta_{j,n})h^{-(1-\Delta\alpha(\ell-i))}\right),\end{split}$
(3.2)

where, for sake of notation, we denoted by $\tilde{w}_{t}^{\alpha_{i}}$ the
approximation of the $t$-th asymptotic expansion coefficient associated with
$c_{i}g_{\alpha_{i}}$ and the term $\tilde{w}_{0}^{\alpha_{i}}$ coincides with
the evaluations of $c_{i}g_{\alpha_{i}}$.

We highlight that the terms in brackets on the right-hand side of the equality
act as possible asymptotic expansion coefficients associated with the GLT
momentary symbols $g_{n}$ of
$\frac{h^{\alpha_{\ell}}}{\Delta{\alpha}}\mathcal{T}_{n}$. Note that formula
(3.2) can be rewritten in compact form as

$\begin{split}\lambda_{j}\left(\frac{h^{\alpha_{\ell}}}{\Delta{\alpha}}\mathcal{T}_{n}\right)\approx\sum_{t=1}^{\nu}h^{t}\left(\tilde{w}_{t}^{{\alpha_{\ell}}}(\theta_{j,n})+\sum_{i=\ell}^{1}\tilde{w}_{t-1}^{{\alpha_{i}}}(\theta_{j,n})h^{-(1-\Delta\alpha(\ell-i))}\right).\end{split}$

Hence, it is easy to see that the GLT momentary symbols correspond to take
$\nu$ of the asymptotic expansion equal to 1.

In the following we consider the cases where $\ell=2$ and $\ell=5$ and
$\ell=n$ as in Section 4 of [23] and confirming at least numerically the
conjecture in (3.2) for a fixed $\nu=4$.

###  3.1. Examples:

For $\ell=2$, $\Delta\alpha$ is $\frac{1}{2}$ and the matrix in (3.1) becomes

$2h^{\alpha_{2}}\mathcal{T}_{n}=c_{2}T_{n}(g_{\alpha_{2}})+c_{1}h^{\frac{1}{2}}T_{n}(g_{\alpha_{1}}),$

where $\alpha_{1}=\frac{5}{4}$ and $\alpha_{2}=\frac{7}{4}$.  Exploiting the
procedure based on formula (3.2) with $\nu=4$, we compute an approximation of
the eigenvalues of
$c_{2}T_{n}(g_{\alpha_{2}})+c_{1}h^{\frac{1}{2}}T_{n}(g_{\alpha_{1}})$ by

$\begin{split}c_{2}g_{\alpha_{2}}({{\theta}_{j,n}})+&h\left[\tilde{w}^{{\alpha_{2}}}_{1}({\tilde{\theta}_{j,n}})+h^{-\frac{1}{2}}c_{1}g_{\alpha_{1}}({{\theta}_{j,n}})\right]+\\\
&\sum_{t=2}^{\nu}h^{t}\left[\tilde{w}^{{\alpha_{2}}}_{t}({{\theta}_{j,n}})+h^{-\frac{1}{2}}\tilde{w}^{{\alpha_{1}}}_{t-1}({{\theta}_{j,n}})\right],\end{split}$

for $j=1,\dots,n$. We consider the cases where $n=100,500,1000$ using an
initial grid with $n_{1}=10$ points and we compare the aforementioned
approximations with those obtained by the evaluations of GLT and GLT momentary
symbols associated with the sequence described by the matrices in (3.1). In
Figure 3 we can observe that the approximation of the spectrum obtained
computing the evaluations of the GLT momentary symbols is better with respect
to that provided by the evaluations $c_{\ell}g_{\alpha_{\ell}}(\theta_{j,n})$.
Moreover, the error of the approximation significantly reduces for almost all
the eigenvalues combining the notions of GLT momentary symbols with the
asymptotic expansion described before, see Figure 4. Note that the particular
shape of the asymptotic expansion error depends on the fact that in
correspondence with the grid points $\theta_{j,n_{t}}$, $t=1,\dots,\nu$ the
quantities $\tilde{w}_{t}^{\alpha_{i}}$ are calculated exactly by the
extrapolation-interpolation procedure. Moreover, note that the accuracy of the
approximation via the combination of GLT momentary symbols and spectral
asymptotic expansion seems to decrease corresponding to the maximum
eigenvalue. Actually, this behavior is expected by the theory of the
asymptotic expansion. Indeed, it is a consequence of the fact that the
involved symbols are not trigonometric polynomials and in particular they
become non-smooth when periodically extended on the real line.

Figure 3. Approximation of the eigenvalues of
$\frac{h^{\alpha_{\ell}}}{\Delta\alpha}\mathcal{T}_{n}$, $\ell=2$ by the
samplings of the GLT and GLT momentary symbols, and making use of the
momentary asymptotic expansion (MAE) with $\nu=4$ for $n=100,500,1000$ with an
initial grid of $n_{1}=10$ points.

Figure 4. Absolute errors of the approximation of the eigenvalues of
$\frac{h^{\alpha_{\ell}}}{\Delta\alpha}\mathcal{T}_{n}$, $\ell=2$ by the
samplings of the GLT and GLT momentary symbols, and making use of the
momentary asymptotic expansion (MAE) with $\nu=4$ for $n=100,500,1000$ with an
initial grid of $n_{1}=10$ points.

Following the analogous procedure, we consider the case where $\ell=5$ and
$\ell=n$ which are associated with $\Delta\alpha=\frac{1}{5}$ and
$\Delta\alpha=\frac{1}{n}$, respectively. In Figures 5 and 7 we plot the
approximations of the eigenvalues given by the three presented strategies for
$\ell=5$ and $\ell=n$.  Again, we obtain numerical confirmation that the
combination of the notions of GLT momentary symbols and asymptotic expansion
provides accurate results even for moderate sizes, as confirmed by the error
plots in Figures 6 and 8, for $n=100,500,1000$.  The good outcome of the
presented numerical tests gives ground for a finer analysis of the spectral
features of the matrices considered in the case left open in [23], which
arises when the integral partition width is asymptotic to the adopted
discretization step. That is, when in formula (3.1) we take $\alpha_{j}=jh$,
$j=1,\ldots,\ell$, $\ell=n$.

Moreover, efficient and fast algorithms which exploit the concept of momentary
symbols can be studied for computing the singular values and eigenvalues of
$T_{n}(f)$ with its possible block, and variable coefficients generalizations
and this will be investigated in the future.

Figure 5. Approximation of the eigenvalues of
$\frac{h^{\alpha_{\ell}}}{\Delta\alpha}\mathcal{T}_{n}$, $\ell=5$ by the
samplings of the GLT and GLT momentary symbols, and making use of the
momentary asymptotic expansion (MAE) with $\nu=4$ for $n=100,500,1000$ with an
initial grid of $n_{1}=10$ points.

Figure 6. Absolute errors of the approximation of the eigenvalues of
$\frac{h^{\alpha_{\ell}}}{\Delta\alpha}\mathcal{T}_{n}$, $\ell=5$ by the
samplings of the GLT and GLT momentary symbols, and making use of the
momentary asymptotic expansion (MAE) with $\nu=4$ for $n=100,500,1000$ with an
initial grid of $n_{1}=10$ points.

Figure 7. Approximation of the eigenvalues of
$\frac{h^{\alpha_{\ell}}}{\Delta\alpha}\mathcal{T}_{n}$, $\ell=n$ by the
samplings of the GLT and GLT momentary symbols, and making use of the
momentary asymptotic expansion (MAE) with $\nu=4$ for $n=100,500,1000$ with an
initial grid of $n_{1}=10$ points.

Figure 8. Absolute errors of the approximation of the eigenvalues of
$\frac{h^{\alpha_{\ell}}}{\Delta\alpha}\mathcal{T}_{n}$, $\ell=n$ by the
samplings of the GLT and GLT momentary symbols, and making use of the
momentary asymptotic expansion (MAE) with $\nu=4$ for $n=100,500,1000$ with an
initial grid of $n_{1}=10$ points.

##  4\. Concluding remarks

The main focus of this paper has been the characterization of the spectrum and
the singular values of the coefficient matrix stemming from the approximation
with space-time grid for a parabolic diffusion problem and from the
approximation of distributed order fractional equations. For this purpose we
employed the classical GLT theory and the new concept of GLT momentary
symbols. The first has permitted to describe the singular value or eigenvalue
asymptotic distribution of the sequence of the coefficient matrices. The
latter has permitted to derive a function, able to describe the singular value
or eigenvalue distribution of the matrix of the sequence, even for small
matrix-sizes, but under given assumptions. In particular, we exploited the
notion of GLT momentary symbols and we used it in combination with the
interpolation-extrapolation algorithms based on the spectral asymptotic
expansion of the involved matrices.

Many questions remain and below we list open problems to be considered in
future researches.

  * •

More examples of the use of GLT momentary symbols in a non-Toeplitz setting;

  * •

The application of GLT momentary symbol in a pure Toeplitz setting, but of
very involved nature, like that expressed in relation (3.1). The use of GLT
momentary symbol for the analysis of efficient iterative solvers, also of
multigrid type, of linear systems as those appearing in (2.1), also with the
inclusion of variable coefficients.

### Acknowledgment

This work was partially supported by INdAM-GNCS. Moreover, the work of
Isabella Furci was also supported by the Young Investigator Training Program
2020 (YITP 2019) promoted by ACRI.

## References

  * [1] F. Avram.  On bilinear forms in Gaussian random variables and Toeplitz matrices.  Probab. Theory Related Fields, 79(1):37–45, 1988.
  * [2] G. Barbarino, C. Garoni, and S. Serra-Capizzano.  Block generalized locally Toeplitz sequences: theory and applications in the multidimensional case.  Electron. Trans. Numer. Anal., 53:113–216, 2020.
  * [3] G. Barbarino, C. Garoni, and S. Serra-Capizzano.  Block generalized locally Toeplitz sequences: theory and applications in the unidimensional case.  Electron. Trans. Numer. Anal., 53:28–1112, 2020.
  * [4] G. Barbarino and S. Serra-Capizzano.  Non-Hermitian perturbations of Hermitian matrix-sequences and applications to the spectral analysis of the numerical approximation of partial differential equations.  Numer. Linear Algebra Appl., 27(3):e2286, 31, 2020.
  * [5] R. Bhatia.  Matrix Analysis.  Springer-Verlag, New York, 1997.
  * [6] M. Bolten, S-E. Ekström, and I. Furci.  Momentary Symbol: Spectral Analysis of Structured Matrices.  https://arxiv.org/abs/2010.06199.
  * [7] E. Bozzo and C. Di Fiore.  On the use of certain matrix algebras associated with discrete trigonometric transforms in matrix displacement decomposition.  SIAM. J. Matrix Anal. Appl., 16:312–326, 1995.
  * [8] H. Brezis.  Functional Analysis, Sobolev Spaces and Partial Differential Equations.  Springer, New York, 2011.
  * [9] R. Chan and X. Jin.  An introduction to iterative Toeplitz solvers.  Society for Industrial and Applied Mathematics (SIAM), 2007.
  * [10] G. Calcagni.  Towards multifractional calculus.  Front. Phys., 6:58, 2018.
  * [11] M. Caputo and M. Fabrizio.  The kernel of the distributed order fractional derivatives with an application to complex materials.  Fractal Fract., 1(1):13, 2017.
  * [12] T. Ceccherini-Silberstein, F. Scarabotti, and F. Tolli.  Harmonic Analysis on Finite Groups.  Cambridge University Press, 2008.
  * [13] W. Ding, S. Patnaik, S. Sidhardh, and F. Semperlotti.  Applications of distributed-order fractional operators: A review.  Entropy, 23(1):110, 2021.
  * [14] M. Donatelli, M. Mazza, and S. Serra-Capizzano.  Spectral analysis and structure preserving preconditioners for fractional diffusion equations.  J. Comput. Phys., 307:262–279, 2016.
  * [15] M. Donatelli, M. Mazza, and S. Serra-Capizzano.  Spectral analysis and multigrid methods for finite volume approximations of space-fractional diffusion equations.  SIAM J. Sci. Comput., 40(6):A4007–A4039, 2018.
  * [16] S.-E. Ekström, I. Furci, and S. Serra-Capizzano.  Exact formulae and matrix-less eigensolvers for block banded symmetric Toeplitz matrices.  BIT, 58:937–968, 2018.
  * [17] S.-E. Ekström, and C. Garoni,  A matrix-less and parallel interpolation–extrapolation algorithm for computing the eigenvalues of preconditioned banded symmetric Toeplitz matrices.  Numer. Algor., 80, 819–848, 2019.
  * [18] S.-E. Ekström, C. Garoni, A. Jozefiak, and J. Perla.  Eigenvalues and Eigenvectors of Tau Matrices with Applications to Markov Processes and Economics.  Linear Algebra Appl., 627, 41–71, 2021.
  * [19] C. Estatico and S. Serra-Capizzano.  Superoptimal approximation for unbounded symbols.  Linear Algebra Appl., 428(2-3):564–585, 2008.
  * [20] C. Garoni and S. Serra-Capizzano.  The Theory of Generalized Locally Toeplitz Sequences: Theory and Applications, Vol. I. Springer Monographs in Mathematics, Berlin, 2017.
  * [21] C. Garoni and S. Serra-Capizzano.  The Theory of Generalized Locally Toeplitz Sequences: Theory and Applications, Vol. II. Springer Monographs in Mathematics, Berlin, 2018.
  * [22] U. Grenander and G. Szegő.  Toeplitz forms and their applications.  Chelsea Publishing Co., New York, second edition, 1984.
  * [23] M. Mazza, S. Serra-Capizzano, and M. Usman.  Symbol-based preconditioning for Riesz distributed-order space-fractional diffusion equations.  Electr. Trans. Num. Anal., 54:499–513, 2021.
  * [24] M. Ng.  Iterative methods for Toeplitz systems.  Oxford University Press, 2004.
  * [25] S. V. Parter.  On the distribution of the singular values of Toeplitz matrices.  Linear Algebra Appl., 80:115–130, 1986.
  * [26] S. Serra-Capizzano, D. Bertaccini, and G. H. Golub.  How to deduce a proper eigenvalue cluster from a proper singular value cluster in the nonnormal case.  SIAM J. Matrix Anal. Appl., 27(1):82–86, 2005.
  * [27] A. Tesei and M. A. Pozio.  On the uniqueness of bounded solutions to singular parabolic problems.  Discrete Contin. Dyn. Syst., 13:117–137–, 2005.
  * [28] A. Tesei and F. Punzo.  Uniqueness of solutions to degenerate elliptic problems with unbounded coefficients.  Ann. I. H. Poincaré, 26:2001–2024, 2009.
  * [29] P. Tilli.  A note on the spectral distribution of Toeplitz matrices.  Linear and Multilinear Algebra, 45(2-3):147–159, 1998.
  * [30] E. Tyrtyshnikov and N. Zamarashkin.  Spectra of multilevel Toeplitz matrices: advanced theory via simple matrix relationships.  Linear Algebra Appl., 270:15–27, 1998.
  * [31] N. Zamarashkin and E. Tyrtyshnikov.  Distribution of the eigenvalues and singular numbers of Toeplitz matrices under weakened requirements on the generating function.  Mat. Sb., 188(3):83–92, 1997.
  * [32] A. Zygmund.  Trigonometric Series.  Cambridge University Press, Cambridge, 1959.