Applications of mixed finite element method based on Bernstein polynomials in numerical solution of Stokes equations

1.   Introduction

The spectral properties of tridiagonal matrices is a well-studied topic for which a vast literature can be found (e.g. ^{[1,5,16,17,19,25,27,35]}, among others), and even formulae for the corresponding inverse of these matrices has also been discussed over the last decades of twentieth century (see ^[15] and references therein). Recently, taking advantage of basic properties of the Chebyshev polynomials, some authors have established localization theorems for the eigenvalues of real pentadiagonal and heptadiagonal symmetric Toeplitz matrices by expressing them as the zeros of explicit rational functions ^[12,32]. The eigenvalues of a special kind of heptadiagonal matrices were still derived in ^[26] by employing other methods, namely, determinant properties and recurrence relations.

In fact, the above-mentioned matrices are typical examples of a much more wider class called band matrices (see ^[30], page 13), and the idea of having explicit formulas to compute its eigenvalues, eigenvectors or establishing some other properties is both appealing and challenging by reason of their usefulness in many areas of science and engineering (see, for instance, ^{[4,10,11,14,20,24,33]}).

In order to give a contribution on this matter, we shall obtain the eigenvalues of the following $n \times n$ heptadiagonal matrix

as the zeros of explicit rational functions, also providing upper/lower bounds non-depending of any unknown parameter to each of them. Further, we shall compute eigenvectors for these sort of matrices at the expense of the prescribed eigenvalues. To accomplish these purposes, we will obtain an orthogonal block diagonalization for matrix (1.1) where each block is a sum of a diagonal matrix plus dyads, i.e.

where ${{\bf{u}}}_{j}, {{\bf{v}}}_{j}$, $j = 1, \ldots, m$ are $\kappa \times 1$ matrices, by exploiting the modification technique introduced by Fasino in ^[13] for matrices of the type (1.1). This key ingredient allows us to get formulas for the characteristic polynomial of ${{\bf{H}}}_{n}$ on one hand, and for the inverse of ${{\bf{H}}}_{n}$ on the other (assuming, of course, its nonsingularity). With the aim of getting expressions as explicit as possible, we will use not only results concerning the secular equation of diagonal matrices perturbed by the addition of rank-one matrices developed by Anderson in the nineties ^[2], but also a Miller's formula of the eighties for the inverse of the sum of matrices ^[29]. In section four of the paper, applications are given for the established results, showing its potential usage.

Since the class of matrices ${{\bf{H}}}_{n}$ includes the ones considered in ^[12] and ^[32], our statements will extend necessarily the results of these papers. Moreover, the current approach also points a way to achieve localization formulas for the eigenvalues of general symmetric quasi-Toeplitz matrices. In detail, the eigenvalues of any symmetric quasi-Toeplitz matrix enjoying a block diagonalization with diagonal elements of the form (1.2) are precisely the eigenvalues of each one of these diagonal blocks, which in turn can be located/computed by rational functions via Anderson's secular equation.

2.   Auxiliary tools

In this paper, $n$ is generally assumed to be an integer greater or equal to four and $a, b, c, d, \xi, \eta$ in (1.1) will be taken as real numbers; in fact, this last restriction can be discarded because the majority of forthcoming statements remain valid when $a, b, c, d, \xi, \eta$ are complex numbers. Moreover, ${{\bf{S}}}_{n}$ will be the $n \times n$ symmetric, involutory and orthogonal matrix defined by

Our first auxiliary result is an orthogonal diagonalization for the following $n \times n$ heptadiagonal symmetric matrix

Lemma 1. Let $a, b, c, d$ be real numbers and

If ${{\bf{\widehat{H}}}}_{n}$ is the $n \times n$ matrix (2.2), then

where ${\boldsymbol{\Lambda}}_{n} = {{\mathrm{diag}}} \left(\lambda_{1}, \ldots, \lambda_{n} \right)$, and ${{\bf{S}}}_{n}$ is the matrix (2.1).

Proof. Supposing the $n \times n$ matrix

it is a simple matter of routine to verify that

Using the spectral decomposition

where

(i.e. the $\ell$th column of ${{\bf{S}}}_{n}$), it follows

where ${\boldsymbol{\Lambda}}_{n} = {{\mathrm{diag}}} \left(\lambda_{1}, \ldots, \lambda_{n} \right)$, and ${{\bf{S}}}_{n}$ is the matrix (2.1). The proof is complete. □

The next statement is an orthogonal block diagonalization for matrices ${{\bf{H}}}_{n}$ of the form (1.1) and it extends Proposition 3.1 in ^[7], which is valid only for heptadiagonal symmetric Toeplitz matrices.

Lemma 2. Let $a, b, c, d, \xi, \eta$ be real numbers, $\lambda_{k}$, $k = 1, \ldots, n$ be given by (2.3) and ${{\bf{H}}}_{n}$ be the $n \times n$ matrix (1.1).

(a) If $n$ is even,

and

then

where ${{\bf{S}}}_{n}$ is the $n \times n$ matrix (2.1), ${{\bf{P}}}_{n}$ is the $n \times n$ permutation matrix defined by

and

(b) If $n$ is odd,

and

then

where ${{\bf{S}}}_{n}$ is the $n \times n$ matrix (2.1), ${{\bf{P}}}_{n}$ is the $n \times n$ permutation matrix defined by

and

Proof. Consider $a, b, c, d, \xi, \eta$ as real numbers, $\lambda_{k}$, $k = 1, \ldots, n$ given by (2.3) and ${{\bf{H}}}_{n}$ as the $n \times n$ matrix (1.1). Setting $\theta : = c + \xi - a$, $\vartheta : = d + \eta - b$,

and

we have from Lemma 1

where ${{\bf{S}}}_{n}$ is the $n \times n$ matrix (2.1), ${{\bf{\widehat{H}}}}_{n}$ is the $n \times n$ matrix (2.2),

and

Since $\left[{{\bf{G}}}_{n} \right]_{k, \ell} = \left[{{\bf{K}}}_{n} \right]_{k, \ell} = 0$ whenever $k + \ell$ is odd, we can permute rows and columns of ${\boldsymbol{\Lambda}}_{n} + {{\bf{G}}}_{n} + {{\bf{K}}}_{n}$ according to the permutation matrices (2.4c) and (2.5c), yielding: for $n$ even,

where ${{\bf{P}}}_{n}$ is the matrix (2.4c), ${\boldsymbol{\Upsilon}}_{\frac{n}{2}} = {{\mathrm{diag}}}(\lambda_{1}, \lambda_{3}, \ldots, \lambda_{n-1})$, ${\boldsymbol{\Delta}}_{\frac{n}{2}} = {{\mathrm{diag}}}(\lambda_{2}, \lambda_{4}, \ldots, \lambda_{n})$ and ${{\bf{x}}}$, ${{\bf{y}}}$ are given by (2.4a); for $n$ odd,

with ${{\bf{P}}}_{n}$ defined in (2.5c), ${\boldsymbol{\Upsilon}}_{\frac{n+1}{2}} = {{\mathrm{diag}}}(\lambda_{1}, \lambda_{3}, \ldots, \lambda_{n})$, ${\boldsymbol{\Delta}}_{\frac{n-1}{2}} = {{\mathrm{diag}}}(\lambda_{2}, \lambda_{4}, \ldots, \lambda_{n-1})$ and ${{\bf{v}}}$, ${{\bf{w}}}$ defined by (2.5a). The proof is complete. □

Remark 1. Let us point out that the decomposition for real heptadiagonal symmetric Toeplitz matrices established in Proposition 3.1 of ^[7] at the expense of the bordering technique is no more useful for matrices having the shape (1.1). As consequence, some results stated by these authors will be necessarily extended, particularly, the referred decomposition and a formula to compute the determinant of real heptadiagonal symmetric Toeplitz matrices (Corollary 3.1 of ^[7]).

3.   Main results

3.1.   Determinant of ${{\bf{H}}}_{n}$

The orthogonal block diagonalization established in Lemma 2 will lead us to an explicit formula for the determinant of the matrix ${{\bf{H}}}_{n}$.

Theorem 1. Let $a, b, c, d, \xi, \eta$ be real numbers, $\lambda_{k}$, $k = 1, \ldots, n$ be given by (2.3), $x_{k} = \sin \left(\frac{k \pi}{n+1} \right)$, $k = 1, \ldots, 2n$ and ${{\bf{H}}}_{n}$ the $n \times n$ matrix (1.1). If $\theta : = c + \xi - a$, $\vartheta : = d + \eta - b$ and

(a) $n$ is even, then

(b) $n$ is odd, then

Proof. Since both assertions can be proven in the same way, we only prove (a). Consider $a, b, c, d, \xi, \eta$ are real numbers, $x_{k} = \sin \left(\frac{k \pi}{n+1} \right)$, $k = 1, \ldots, 2n$, $\lambda_{k}$, $k = 1, \ldots, n$ as given by (2.3), $\theta : = c + \xi - a$, $\vartheta : = d + \eta - b$ and the notations used in Lemma 2. The determinant formula for block-triangular matrices (see ^[21], page 185) and Lemma 2 ensure $\det({{\bf{H}}}_{n}) = \det \left({\boldsymbol{\Phi}}_{\frac{n}{2}} \right) \det \left({\boldsymbol{\Psi}}_{\frac{n}{2}} \right)$. We shall first assume $\lambda_{k} \neq 0$ for all $k = 1, \ldots, n$,

and

Putting ${\boldsymbol{\Upsilon}}_{\frac{n}{2}} : = {{\mathrm{diag}}} \left(\lambda_{1}, \lambda_{3}, \ldots, \lambda_{n-1} \right)$ and ${\boldsymbol{\Delta}}_{\frac{n}{2}} : = {{\mathrm{diag}}} \left(\lambda_{2}, \lambda_{4}, \ldots, \lambda_{n} \right)$, we have

and

(see ^[29], pages 69 and 70), i.e.

Since both sides of (3.3) are polynomials in the variables $a, b, c, d, \xi, \eta$, conditions (3.1a)–(3.2c) as well as $\lambda_{k} \neq 0$ can be dropped, and (3.3) is valid more generally. □

Example 1. Suppose the following symmetric quasi-Toeplitz matrix

(when $\xi = a$, we have a pentadiagonal symmetric Toeplitz matrix). Let us notice that Theorem 3 of ^[12] cannot be employed to compute $\det({{\bf{T}}}_{n})$. However, according to our Theorem 1 we get (with $d = 0$, $\eta = b$ and $\vartheta = 0$)

where

and $x_{k} = \sin \left(\frac{k \pi}{n+1} \right)$, $k = 1, \ldots, 2n$. Moreover, if $\xi = a - c$ in ${{\bf{T}}}_{n}$, then $\det({{\bf{T}}}_{n})$ simply turns into $\lambda_{1} \lambda_{2} \ldots \lambda_{n}$ (let us stress that this includes the particular case $c = 0$, i.e. the determinant of tridiagonal symmetric Toeplitz matrices).

3.2.   Eigenvalue localization for ${{\bf{H}}}_{n}$

The following lemma will allows us to express the eigenvalues of key matrices in this paper as the zeros of explicit rational functions providing, additionally, explicit upper and lower bounds for each one. We will denote the Euclidean norm by $\left\lVert{\; \cdot \; } \right\rVert$.

Lemma 3. Let $a, b, c, d, \xi, \eta$ be real numbers and $\lambda_{k}$, $k = 1, \ldots, n$ be given by (2.3).

(a) If $n$ is even,

ⅰ. ${{\bf{x}}}, {{\bf{y}}}$ are given by (2.4a) and the eigenvalues of

are not of the form $a + 2b \cos \left[\frac{(2k-1)\pi}{n+1} \right] + 2c \cos \left[\frac{2(2k-1) \pi}{n+1} \right] + 2d \cos \left[\frac{3(2k-1) \pi}{n+1} \right]$, $k = 1, \ldots, \frac{n}{2}$, then the eigenvalues of (3.4a) are precisely the zeros of the rational function

Moreover, if $\mu_{1} \leqslant \mu_{2} \leqslant \ldots \leqslant \mu_{\frac{n}{2}}$ are the eigenvalues of (3.4a) and $\lambda_{\tau(1)} \leqslant \lambda_{\tau(3)} \leqslant \ldots \leqslant \lambda_{\tau(n-1)}$ are arranged in a nondecreasing order by some bijection $\tau$ defined in $\{1, 3, \ldots, n-1\}$, then

for each $k = 1, \ldots, \tfrac{n}{2}$.

ⅱ. ${{\bf{v}}}, {{\bf{w}}}$ are given by (2.4b) and the eigenvalues of

are not of the form $a + 2b \cos \left(\frac{2k\pi}{n+1} \right) + 2c \cos \left(\frac{4k \pi}{n+1} \right) + 2d \cos \left(\frac{6k \pi}{n+1} \right)$, $k = 1, \ldots, \frac{n}{2}$, then the eigenvalues of (3.5a) are precisely the zeros of the rational function

Furthermore, if $\nu_{1} \leqslant \nu_{2} \leqslant \ldots \leqslant \nu_{\frac{n}{2}}$ are the eigenvalues of (3.5a) and $\lambda_{\sigma(2)} \leqslant \lambda_{\sigma(4)} \leqslant \ldots \leqslant \lambda_{\sigma(n)}$ are arranged in a nondecreasing order by some bijection $\sigma$ defined in $\{2, 4, \ldots, n\}$, then

for every $k = 1, \ldots, \tfrac{n}{2}$.

(b) If $n$ is odd,

ⅰ. ${{\bf{x}}}, {{\bf{y}}}$ are given by (2.5a) and the eigenvalues of

are not of the form $a + 2b \cos \left[\frac{(2k-1)\pi}{n+1} \right] + 2c \cos \left[\frac{2(2k-1) \pi}{n+1} \right] + 2d \cos \left[\frac{3(2k-1) \pi}{n+1} \right]$, $k = 1, \ldots, \frac{n+1}{2}$, then the eigenvalues of (3.6a) are precisely the zeros of the rational function

Moreover, if $\mu_{1} \leqslant \mu_{2} \leqslant \ldots \leqslant \mu_{\frac{n+1}{2}}$ are the eigenvalues of (3.6a) and $\lambda_{\tau(1)} \leqslant \lambda_{\tau(3)} \leqslant \ldots \leqslant \lambda_{\tau(n)}$ are arranged in a nondecreasing order by some bijection $\tau$ defined in $\{1, 3, \ldots, n\}$, then

for any $k = 1, \ldots, \tfrac{n+1}{2}$.

ⅱ. ${{\bf{v}}}, {{\bf{w}}}$ are given by (2.5b) and the eigenvalues of

are not of the form $a + 2b \cos \left(\frac{2k\pi}{n+1} \right) + 2c \cos \left(\frac{4k \pi}{n+1} \right) + 2d \cos \left(\frac{6k \pi}{n+1} \right)$, $k = 1, \ldots, \frac{n-1}{2}$, then the eigenvalues of (3.7a) are precisely the zeros of the rational function

Furthermore, if $\nu_{1} \leqslant \nu_{2} \leqslant \ldots \leqslant \nu_{\frac{n-1}{2}}$ are the eigenvalues of (3.7a) and $\lambda_{\sigma(2)} \leqslant \lambda_{\sigma(4)} \leqslant \ldots \leqslant \lambda_{\sigma(n-1)}$ are arranged in a nondecreasing order by some bijection $\sigma$ defined in $\{2, 4, \ldots, n-1\}$, then

for all $k = 1, \ldots, \tfrac{n-1}{2}$.

Proof. Suppose real numbers $a, b, c, d, \xi, \eta$, $\lambda_{k}$, $k = 1, \ldots, n$ given by (2.3) and put $\theta : = c + \xi - a$, $\vartheta : = d + \eta - b$. We shall denote by $\mathcal{S}(k, m)$ the collection of all $k$-element subsets of $\{1, 2, \ldots, m \}$ written in increasing order; additionally, for any rectangular matrix ${{\bf{M}}}$, we shall indicate by $\det \left({{\bf{M}}}[I, J] \right)$ the minor determined by the subsets $I = \left\{i_{1} < i_{2} < \ldots < i_{k} \right\}$ and $J = \left\{j_{1} < j_{2} < \ldots < j_{k} \right\}$. Supposing $\theta \neq 0$,

and

Theorem 1 of ^[2] ensures that $\zeta$ is an eigenvalue of (3.4a) if, and only if,

provided that $\zeta$ is not an eigenvalue of ${{\mathrm{diag}}} \left(\lambda_{1}, \lambda_{3}, \ldots, \lambda_{n-1} \right)$. Since

we obtain (3.4b). Considering now $\theta = 0$ and setting

we still have that $\zeta$ is an eigenvalue of (3.4a) if, and only if,

assuming that $\zeta$ is not an eigenvalue of ${{\mathrm{diag}}} \left(\lambda_{1}, \lambda_{3}, \ldots, \lambda_{n-1} \right)$. Hence,

and (3.4b) is established. Let $\mu_{1} \leqslant \mu_{2} \leqslant \ldots \leqslant \mu_{\frac{n}{2}}$ be the eigenvalues of (3.4a) and $\lambda_{\tau(1)} \leqslant \lambda_{\tau(3)} \leqslant \ldots \leqslant \lambda_{\tau(n-1)}$ be arranged in a nondecreasing order by some bijection $\tau$ defined in $\{1, 3, \ldots, n-1\}$. Thus,

for each $k = 1, \ldots, \tfrac{n}{2}$ (see ^[23], page 242). Since the characteristic polynomial of $\theta {{\bf{x}}} {{\bf{x}}}^{\top} + \vartheta {{\bf{x}}} {{\bf{y}}}^{\top} + \vartheta {{\bf{y}}} {{\bf{x}}}^{\top}$ is

we have that its spectrum is

where $\alpha^{\pm} : = \frac{\theta \left\lVert{{{\bf{x}}}} \right\rVert^{2} + 2 \vartheta {{\bf{x}}}^{\top} {{\bf{y}}} \pm \sqrt{\left(\theta \left\lVert{{{\bf{x}}}} \right\rVert^{2} + 2 \vartheta {{\bf{x}}}^{\top} {{\bf{y}}}\right)^{2} - 4\vartheta^{2} \left[\left({{\bf{x}}}^{\top} {{\bf{y}}} \right)^{2} - \left\lVert{{{\bf{x}}}} \right\rVert^{2} \left\lVert{{{\bf{y}}}} \right\rVert^{2} \right]}}{2}$. From the identities

it follows $\lVert {{\bf{x}}} \rVert = \lVert {{\bf{y}}} \rVert = 1$ and ${{\bf{x}}}^{\top} {{\bf{y}}} = 0$. Hence, (3.8) and (3.9) yields (3.4c). The proofs of the remaining assertions are performed in the same way and so will be omitted. □

The next statement allows us to locate the eigenvalues of ${{\bf{H}}}_{n}$, providing also explicit bounds for each of them.

Theorem 2. Let $a, b, c, d, \xi, \eta$ be real numbers, $\lambda_{k}$, $k = 1, \ldots, n$ be given by (2.3) and ${{\bf{H}}}_{n}$ be the $n \times n$ matrix (1.1).

(a) If $n$ is even, the eigenvalues of ${\boldsymbol{\Phi}}_{\frac{n}{2}}$ in (2.4d) are not of the form $\lambda_{2k-1}$, $k = 1, \ldots, \frac{n}{2}$ and the eigenvalues of ${\boldsymbol{\Psi}}_{\frac{n}{2}}$ in (2.4e) are not of the form $\lambda_{2k}$, $k = 1, \ldots, \frac{n}{2}$, then the eigenvalues of ${{\bf{H}}}_{n}$ are precisely the zeros of the rational functions $f(t)$ and $g(t)$ given by (3.4b) and (3.5b), respectively. Moreover, if $\mu_{1} \leqslant \mu_{2} \leqslant \ldots \leqslant \mu_{\frac{n}{2}}$ are the zeros of $f(t)$ and $\nu_{1} \leqslant \nu_{2} \leqslant \ldots \leqslant \nu_{\frac{n}{2}}$ are the zeros of $g(t)$ (counting multiplicities in both cases), then $\mu_{k}$, $k = 1, \ldots, \tfrac{n}{2}$ and $\nu_{k}$, $k = 1, \ldots, \tfrac{n}{2}$ satisfy (3.4c) and (3.5c), respectively.

(b) If $n$ is odd, the eigenvalues of ${\boldsymbol{\Phi}}_{\frac{n+1}{2}}$ in (2.5d) are not of the form $\lambda_{2k-1}$, $k = 1, \ldots, \frac{n+1}{2}$ and the eigenvalues of ${\boldsymbol{\Psi}}_{\frac{n-1}{2}}$ in (2.5e) are not of the form $\lambda_{2k}$, $k = 1, \ldots, \frac{n-1}{2}$, then the eigenvalues of ${{\bf{H}}}_{n}$ are precisely the zeros of the rational functions $f(t)$ and $g(t)$ given by (3.6b) and (3.7b), respectively. Furthermore, if $\mu_{1} \leqslant \mu_{2} \leqslant \ldots \leqslant \mu_{\frac{n+1}{2}}$ are the zeros of $f(t)$ and $\nu_{1} \leqslant \nu_{2} \leqslant \ldots \leqslant \nu_{\frac{n-1}{2}}$ are the zeros of $g(t)$ (counting multiplicities in both cases), then $\mu_{k}$, $k = 1, \ldots, \tfrac{n+1}{2}$ and $\nu_{k}$, $k = 1, \ldots, \tfrac{n-1}{2}$ satisfy (3.6c) and (3.7c), respectively.

Proof. Suppose $a, b, c, d, \xi, \eta$ are real numbers and $\lambda_{k}$, $k = 1, \ldots, n$ as given by (2.3).

(a) According to Lemma 2 and the determinant formula for block-triangular matrices (see ^[21], page 185), the characteristic polynomial of ${{\bf{H}}}_{n}$ for $n$ even is

where ${\boldsymbol{\Phi}}_{\frac{n}{2}}$ and ${\boldsymbol{\Psi}}_{\frac{n}{2}}$ are given by (2.4d) and (2.4e), respectively, so that the thesis is a direct consequence of Lemma 2.

(b) For $n$ odd, we obtain

where ${\boldsymbol{\Phi}}_{\frac{n+1}{2}}$ and ${\boldsymbol{\Psi}}_{\frac{n-1}{2}}$ are given by (2.5d) and (2.5e), respectively. The conclusion follows from Lemma 2. □

From Geršgorin theorem (see ^[23], Theorem 6.1.1), it can also be stated that all eigenvalues of ${{\bf{H}}}_{n}$ $(n \geqslant 7)$ belong to $[h_{\min}, h_{\max}]$, where

and

Further, all eigenvalues of the $n \times n$ heptadiagonal symmetric Toeplitz matrix

are contained in the interval

where $\varphi(t) = a + 2b \cos(t) + 2c \cos(2t) + 2d \cos(3t)$, $-\pi \leqslant t \leqslant \pi$ (see ^[18], Theorem 6.1). As illustrated, eigenvalues of ${{\bf{H}}}_{n}$ and those of ${{\mathrm{hepta}}}_{n}(d, c, b, a, b, c, d)$ with $a = 0$, $b = -2$, $c = -1$, $d = 2$, $\xi = 9$, $\eta = -7$ are depicted in complex plane for increasing values of $n$.

A distinctive feature of the blue graphics is the existence of two outliers for ${{\bf{H}}}_{n}$, i.e. eigenvalues that do not belong to the interval $\left[-\frac{154}{27}, 7 \right]$, which seems to become just one as $n \rightarrow \infty$. This numerical experiment also reveals that as the matrix size increases, the spectrum of quasi-Toeplitz matrix ${{\bf{H}}}_{n}$ approaches the spectrum of Toeplitz matrix ${{\mathrm{hepta}}}_{n}(2, -1, -2, 1, -2, -1, 2)$ plus the outliers; this is the scenario that is consistent with the study presented in ^[6].

Remark 2. In ^[12] and ^[32], similar localization results were established for the eigenvalues of symmetric Toeplitz matrices (pentadiagonal and heptadiagonal, respectively). The referred papers make use of Chebyshev polynomials and their properties to earn rational functions with a more concise form. However, its statements do not cover the broader class of matrices (1.1).

3.3.   Eigenvectors of ${{\bf{H}}}_{n}$

The decomposition presented in Lemma 2 allows us also to compute eigenvectors for ${{\bf{H}}}_{n}$ in (1.1).

Theorem 3. Let $a, b, c, d, \xi, \eta$ be real numbers, $\lambda_{k}$, $k = 1, \ldots, n$ be given by (2.3) and ${{\bf{H}}}_{n}$ be the $n \times n$ matrix (1.1).

(a) If $n$ is even, ${{\bf{S}}}_{n}$ is the $n \times n$ matrix (2.1), ${{\bf{P}}}_{n}$ is the $n \times n$ permutation matrix (2.4c), the zeros $\mu_{1}, \ldots, \mu_{\frac{n}{2}}$ of (3.4b) are not of the form $\lambda_{2k-1}$, $k = 1, \ldots, \frac{n}{2}$, the zeros $\nu_{1}, \ldots, \nu_{\frac{n}{2}}$ of (3.5b) are not of the form $\lambda_{2k}$, $k = 1, \ldots, \frac{n}{2}$,

and $b \neq d + \eta$, then

is an eigenvector of ${{\bf{H}}}_{n}$ associated to $\mu_{k}$, $k = 1, \ldots, \frac{n}{2}$, and

is an eigenvector of ${{\bf{H}}}_{n}$ associated to $\nu_{k}$, $k = 1, \ldots, \frac{n}{2}$.

(b) If $n$ is odd, ${{\bf{S}}}_{n}$ is the $n \times n$ matrix (2.1), ${{\bf{P}}}_{n}$ is the $n \times n$ permutation matrix (2.4c), the zeros $\mu_{1}, \ldots, \mu_{\frac{n+1}{2}}$ of (3.6b) are not of the form $\lambda_{2k-1}$, $k = 1, \ldots, \frac{n+1}{2}$, the zeros $\nu_{1}, \ldots, \nu_{\frac{n-1}{2}}$ of (3.7b) are not of the form $\lambda_{2k}$, $k = 1, \ldots, \frac{n-1}{2}$,

and $b \neq d + \eta$, then

is an eigenvector of ${{\bf{H}}}_{n}$ associated to $\mu_{k}$, $k = 1, \ldots, \frac{n+1}{2}$, and

is an eigenvector of ${{\bf{H}}}_{n}$ associated to $\nu_{k}$, $k = 1, \ldots, \frac{n-1}{2}$.

Proof. Since both assertions can be proven in the same way, we only prove (a). Let $n$ be even. We can rewrite the matricial equation $(\mu_{k} {{\bf{I}}}_{n} - {{\bf{H}}}_{n}) {{\bf{q}}} = {{\bf{0}}}$ as

where ${{\bf{S}}}_{n}$ is the matrix (2.1), ${{\bf{P}}}_{n}$ is the permutation matrix (2.4c) and ${\boldsymbol{\Phi}}_{\frac{n}{2}}$ and ${\boldsymbol{\Psi}}_{\frac{n}{2}}$ are given by (2.4d) and (2.4e), respectively. Thus,

That is,

for $\alpha \neq 0$ (see ^[8], page 41), and

is a nontrivial solution of (3.12). Thus, choosing $\alpha = 1$, we conclude that (3.10a) is an eigenvector of ${{\bf{H}}}_{n}$ associated to the eigenvalue $\mu_{k}$. Similarly, from $(\nu_{k} {{\bf{I}}}_{n} - {{\bf{H}}}_{n}) {{\bf{q}}} = {{\bf{0}}}$, we have

and

for $\alpha \neq 0$, which is an eigenvector of ${{\bf{H}}}_{n}$ associated to the eigenvalue $\nu_{k}$. □

3.4.   Expression of ${{\bf{H}}}_{n}^{-1}$

The orthogonal block diagonalization presented in Lemma 2 and Miller's formula for the inverse of the sum of nonsingular matrices lead us to an explicit expression for the inverse of ${{\bf{H}}}_{n}$.

Theorem 4. Let $a, b, c, d, \xi, \eta$ be real numbers, $\lambda_{k}$, $k = 1, \ldots, n$ be given by (2.3) and ${{\bf{H}}}_{n}$ be the $n \times n$ matrix (1.1). If $\lambda_{k} \neq 0$ for every $k = 1, \ldots, n$, ${{\bf{H}}}_{n}$ is nonsingular and:

(a) $n$ is even, then

where ${{\bf{S}}}_{n}$ is the $n \times n$ matrix (2.1), ${{\bf{P}}}_{n}$ is the $n \times n$ permutation matrix (2.4c),

with ${\boldsymbol{\Upsilon}}_{\frac{n}{2}} : = {{\mathrm{diag}}} \left(\lambda_{1}, \lambda_{3}, \ldots, \lambda_{n-1} \right)$, ${{\bf{x}}}, {{\bf{y}}}$ given by (2.4a),

and

with ${\boldsymbol{\Delta}}_{\frac{n}{2}} : = {{\mathrm{diag}}} \left(\lambda_{2}, \lambda_{4}, \ldots, \lambda_{n} \right)$, ${{\bf{v}}}, {{\bf{w}}}$ given by (2.5a) and

(b) $n$ is odd, then

where ${{\bf{S}}}_{n}$ is the $n \times n$ matrix (2.1), ${{\bf{P}}}_{n}$ is the $n \times n$ permutation matrix (2.5c),

with ${\boldsymbol{\Upsilon}}_{\frac{n+1}{2}} : = {{\mathrm{diag}}} \left(\lambda_{1}, \lambda_{3}, \ldots, \lambda_{n} \right)$, ${{\bf{x}}}, {{\bf{y}}}$ given by (2.5a),

and

with ${\boldsymbol{\Delta}}_{\frac{n-1}{2}} : = {{\mathrm{diag}}} \left(\lambda_{2}, \lambda_{4}, \ldots, \lambda_{n-1} \right)$, ${{\bf{v}}}, {{\bf{w}}}$ in (2.5b),

Proof. Consider $a, b, c, d, \xi, \eta$ as real numbers, $\lambda_{k} \neq 0$, $k = 1, \ldots, n$ are given by (2.3) and ${{\bf{H}}}_{n}$ in (1.1) is nonsingular. Recall that if ${{\bf{H}}}_{n}$ is nonsingular, then $\rho$ and $\varrho$ in (3.13b) and (3.13d), respectively, are both nonzero. Setting $\theta : = c + \xi - a$, $\vartheta : = d + \eta - b$ and assuming that conditions (3.1a) and (3.1b) are satisfied (note that (3.1c) corresponds to $\rho \neq 0$), we have from the main result of ^[29] (see pages 69 and 70),

and

with ${\boldsymbol{\Upsilon}}_{\frac{n}{2}} : = {{\mathrm{diag}}} \left(\lambda_{1}, \lambda_{3}, \ldots, \lambda_{n-1} \right)$, ${{\bf{x}}}, {{\bf{y}}}$ given by (2.4a) and $\rho$ in (3.13b). In the same way, supposing (3.2a) and (3.2b) (observe that (3.2c) is $\varrho \neq 0$), we obtain

and

where ${\boldsymbol{\Delta}}_{\frac{n}{2}} : = {{\mathrm{diag}}} \left(\lambda_{2}, \lambda_{4}, \ldots, \lambda_{n} \right)$ and ${{\bf{v}}}, {{\bf{w}}}$ are given by (2.5a) and $\varrho$ in (3.13d). Since the nonsingularity of ${{\bf{H}}}_{n}$ and $\lambda_{k} \neq 0$, for all $k = 1, \ldots, n$ are sufficient for both sides of (3.15) and (3.16) to be well-defined, conditions (3.1a), (3.1b), (3.2a) and (3.2b) previously assumed can be dropped. Hence, the block diagonalization provided in (a) of Lemma 2 together with 8.5b of ^[21] (see page 88) establish the thesis in (a). The proof of (b) is analogous, so we will omit the details. □

4.   Applications

4.1.   Matrix derivative operator

It is well known that the fourth derivative can be computed through the following centered finite-formula

(see ^[9], page 556). Consider an interval $[a, b]$ $(a < b)$, a mesh of points $x_{k} = a + k h$, $k = 0, 1, \ldots, N$ where $h = (b - a)/N$ and a function $f\colon [a, b] \longrightarrow \mathbb{R}$, such that $f(a) = 0 = f(b)$. By setting

for some $\alpha\in \mathbb{R}$, the matrix operator corresponding to (4.1) for the fourth derivative is

A remarkable example that involves the fourth derivative is the ordinary differential equation that governs the deflection of a laterally loaded symmetrical beam of length $L$,

where ${{\mathrm{E}}}$ is the modulus of elasticity of the beam material, ${{\mathrm{I}}}(x)$ is the moment of inertia of the beam cross section and $q(x)$ is the distributed load. The ordinary differential equation (4.3) can be equipped with the boundary conditions $y(0) = 0 = y(L)$ (see, for instance, ^[22]).

The eigenvalues of derivative matrices are very useful. In fact, they can be compared with those of the exact (continuous) derivative operator to gauge the accuracy of the finite difference approximation. On the other hand, in the context of partial differential equations, the eigenvalues of the spatial operator is considered along with the stability diagram of the time-integration scheme to evaluate the stability of the numerical solution for the partial differential equation ^[3]. The statements of subsection 3.2 can be employed to locate (bound) the eigenvalues of (4.2).

Another example of a derivative matrix is

which appears in the discretization of the second-derivative operator via three-point centered finite-difference formula with Neumann boundary conditions $f'(x_{0}) = a$ and $f'(x_{N}) = b$ (see ^[3], pages 133 and 134). Our results can also be used to locate (bound) its eigenvalues by noticing that the eigenvalues of (4.4) and

are exactly the same.

4.2.   Beveridge-Nelson smoother in ARIMA(1, 1, 0) process

Consider $n$ pairs of observations $(x_{1}, y_{1}), (x_{2}, y_{2}), \ldots, (x_{n}, y_{n})$ such that

where $r$ is the regression function to be estimated. The estimator of $r(x)$ is usually denoted by $\widehat{r}(x)$ and called smoother. An estimator $\widehat{r}$ of $r$ is a linear smoother if, for each $x$, there exists a vector ${\boldsymbol{\varsigma}}(x) = (\varsigma_{1}(x), \ldots, \varsigma_{n}(x))^{\top}$ such that

Defining the vector of fitted values ${{\bf{\widehat{y}}}} = (\widehat{r}_{n}(x_{1}), \ldots, \widehat{r}_{n}(x_{n}))^{\top}$, it follows

where ${\boldsymbol{\Sigma}}$ is an $n \times n$ matrix whose $k^{{{\mathrm{th}}}}$ row is ${\boldsymbol{\varsigma}}(x_{k})^{\top}$, called the smoothing matrix and ${{\bf{y}}} = (y_{1}, \ldots, y_{n})^{\top}$ (see ^[34], page 66).

The eigendecomposition of the smoothing matrix ${\boldsymbol{\Sigma}}$ provides a useful characterization of the properties of a smoother. In fact, if ${\boldsymbol{\Sigma}} = \sum_{k = 1}^{n} \lambda_{k} {\boldsymbol{\sigma}}_{k} {\boldsymbol{\sigma}}_{k}^{\top}$ is the spectral decomposition of the smoothing matrix, where $\lambda_{k}$ are the ordered eigenvalues and ${\boldsymbol{\sigma}}_{k}$ the corresponding eigenvectors, we can meaningfully decompose the fit as ${{\bf{\widehat{y}}}} = \sum_{k = 1}^{n} \alpha_{k} \lambda_{k} {\boldsymbol{\sigma}}_{k}$, where the eigenvectors ${\boldsymbol{\sigma}}_{k}$ illustrate what sequences are preserved or compressed via a scalar multiplication and $\alpha_{k}$ are the specific coefficients of the projection of ${{\bf{y}}}$ onto the space spanned by the eigenvectors ${\boldsymbol{\sigma}}_{k}$, ${{\bf{y}}} = \sum_{k = 1}^{n} \alpha_{k} {\boldsymbol{\sigma}}_{k}$. Moreover, ${{\mathrm{tr}}}({\boldsymbol{\Sigma}}) = \sum_{k = 1}^{n} \lambda_{k}$ provides the number of degrees of freedom of a smoother, which is a measure of the equivalent number of parameters used to obtain the fit ${{\bf{\widehat{y}}}}$ that allows us to compare alternative filters according to their degree of smoothing (see ^[28] and the references therein).

The smoothing matrix associated to the Beveridge-Nelson smoother (see ^[31] for details) when the observed series is generated by an ${{\mathrm{ARIMA}}}$(1, 1, 0) model with $-1 < \phi < 0$ and (half) bandwidth filter $m = 1$ is the following tridiagonal matrix:

(see ^[28]). Since the matrices ${\boldsymbol{\Sigma}}$ and

share the same eigenvalues, we are able to locate (bound) the eigenvalues of ${\boldsymbol{\Sigma}}$ by using results of subsection 3.2. Moreover, at the expense of the prescribed eigenvalues, an eigendecomposition for ${\boldsymbol{\Sigma}}$ can also be obtained at the expense of statements in subsection 3.3.

5.   Conclusions

In this paper, a procedure to express the eigenvalues and associated eigenvectors of a symmetric heptadiagonal quasi-Toeplitz matrix was presented, as well as an explicit formula for its inverse. The proposed method allowed us to get rational functions to locate the eigenvalues and closed-form formulas to the corresponding eigenvectors for the class of matrices under analysis, which cannot be considered in recent works on this subject, but most of all leave an open door for additional statements on symmetric quasi-Toeplitz matrices in general. The numerical example provided to highlight the differences between the quasi-Toeplitz and Toeplitz cases also raised some open questions. Indeed, despite Geršgorin theorem leading us to an interval containing all eigenvalues of generic quasi-Toeplitz matrices, it would be interesting to have a more precise tool, as in the "pure" Toeplitz case. A method that could predict the number of outliers and its asymptotic behavior when $n$ tends to infinity would be also very welcome. Of course, another open problem closely related to the content of this paper would be the obtention of a block diagonalization for nonsymmetric quasi-Toeplitz matrices in the same spirit of Lemma 2.

Use of AI tools declaration

The author declares he has not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The author would like to thank Professor Yongjian Hu for the invitation to submit the manuscript, and also to anonymous referees for the careful reading of it as well as their very constructive comments, which greatly improved the final version of the paper.

This work is funded by national funds through the FCT - Fundação para a Ciência e a Tecnologia, I.P., under the scope of project UIDB/04035/2020 (GeoBioTec).

Conflict of interest

The author declares there is no conflict of interest.