Spectral Distribution and Numerical Methods for Rational Eigenvalue Problems

Abstract

Rational eigenvalue problems (REPs) have important applications in engineering applications and have attracted more and more attention in recent years. Based on the theory of low-rank modification, we discuss the spectral properties and distribution of the symmetric rational eigenvalue problems, and present two numerical iteration methods for the above REPs. Numerical experiments demonstrate the effectiveness of our proposed methods.

1. Introduction

We consider the following rational eigenvalue problem (REP)

$$R(\lambda )x=0,$$

(1)

where $R(\lambda )\in {\mathcal{R}}^{n\times n}$ is a matrix rational function with respect to the scalar parameter $\lambda \in \mathcal{R}$ and the details for the degree of $R(\lambda )$ can be seen in [1]. As we know, $\lambda$ is an eigenvalue of the problem (1) if and only if it satisfies the characteristic equation $\det (R(\lambda ))=0$ where $\det (·)$ denotes the determinant of its following matrix. The nonzero vector $x\in {\mathcal{R}}^n$ and the two-tuple $(\lambda ,x)$ are called as the corresponding eigenvector of $\lambda$ and an eigenpair of the REP (1), respectively. The REP arises in a wide variety of applications including vibration of fluid–solid structures [2], optimization of acoustic emissions of high-speed trains [3], free vibration of plates with elastically attached masses [4], free vibrations of a structure with a viscoelastic constitutive relation describing the behavior of a material [5,6], electronic structure calculations of quantum dots [7,8], and so on.

More precisely, in this paper we consider that $R(\lambda )$ is shown as follows:

$$R\left(\lambda \right)=P\left(\lambda \right)+\sum _{i=1}^N\frac{s_i\left(\lambda \right)}{q_i\left(\lambda \right)}{B}_i,$$

(2)

where $P\left(\lambda \right)={\lambda }^d{A}_d+\cdots +\lambda A_1+A_0$ is a matrix polynomial, $s_i(\lambda )$ and $q_i(\lambda )$ are coprime scalar polynomials of degrees $m_i$ and $n_i$ with $m_i<n_i$, respectively, and $A_i,B_i\in {\mathcal{R}}^{n\times n}$ are all constant matrices.

At present, there are mainly three types of numerical methods to compute the eigenvalues of the REP (1). The first type of numerical method is to solve the REP via a brute-force approach. That is to say, multiply the both sides of (2) by all scalar polynomials $q_i(\lambda )$, which results in a $d+{\sum }_{i=1}^N{n}_i$ order polynomial eigenvalue problem (PEP). Nevertheless, these methods are not as efficient as required for the large-scale problems, especially when the term ${\sum }_{i=1}^N{n}_i$ is not small enough. Moreover, the corresponding PEP would have more extra eigenvalues than the original REP (1). The second type of numerical method is to linearize the REP into a PEP with some specific tricks. For example, Su and Bai [9] presented a linearization-based method by converting the REP into a well-studied PEP and preserved the structures and properties of the original REP. Dopico and González-Pizarro [10] proposed a compact rational Krylov method for the large-scale REP. The third type of numerical method is to treat the REP as the general nonlinear eigenvalue problems, and solve them via a nonlinear eigensolver, see, e.g., [11,12,13,14,15,16,17,18,19,20,21,22,23,24]. Although the abovementioned methods can solve the REP well, they ignore some structures and properties of the original rational eigenvalue problems. To overcome this disadvantage, it is necessary to study spectral properties and distribution of the REP first, and then try to put forward some effective numerical methods according to these properties. As far as we know, there are engineering applications that require the computation of only some of the eigenvalues lying within an interval [5]. Therefore, in this paper we focus on some numerical methods to compute eigenpairs of the REP in an interval.

The rest of the paper is organized as follows. Section 2 briefly introduces some preliminary results. Section 3 discusses the spectral properties and the distribution of the rational eigenvalues. In Section 4, we develop two numerical methods for solving the symmetric REP based on the spectral properties. Some numerical examples are given to show the effectiveness of the proposed methods in Section 5. Finally, we give some concluding remarks in Section 6.

For convenience, we use the following notations: I denotes the identity matrix of suitable size. $e_j$ denotes the jth column of the identity matrix I. The superscript T denotes the transpose of a vector or a matrix, respectively. $∥·∥$ denotes the Euclidean norm of a vector or a matrix. ${∥·∥}_F$ denotes the Frobenius norm of a matrix. $\langle x,y\rangle$ denotes the inner product of vector x and vector y.

2. Preliminaries

In this section, to facilitate the theoretical analysis and further obtain the main results for rational eigenvalue problems (1) and (2), the following useful assumptions are addressed.

Hypothesis 1 (H1).

Coefficient matrices $A_i$ are the symmetric positive definite matrices with $i=1,\dots ,d$ and $A_0$ is symmetric.

Hypothesis 2 (H2).

Matrices $B_i$ are the low-rank symmetric positive semidefinite matrices with $i=1,2,\dots ,N$.

Hypothesis 3 (H3).

Rational functions $r_i(\lambda )={\displaystyle \frac{s_i\left(\lambda \right)}{q_i\left(\lambda \right)}}$ are monotonically increasing functions with respect to parameter λ on the intervals separated by the zeros of polynomials $q_i\left(\lambda \right)$ where $i=1,2,\dots ,N$.

Remark 1.

The following example provides an intuitive illustration for assumptionH3.

Example 1.

Assume that $r_i(\lambda )={\displaystyle \frac{a_1}{{\lambda }_1-\lambda }}+{\displaystyle \frac{a_2}{{\lambda }_2-\lambda }}+\cdots +{\displaystyle \frac{a_k}{{\lambda }_k-\lambda }}$, where ${\lambda }_1<{\lambda }_2<\cdots <{\lambda }_k$ and $a_1,a_2,\dots ,a_k>0$. We can easily get that $r_i^{\prime }(\lambda )={\displaystyle \frac{a_1}{{({\lambda }_1-\lambda )}^2}}+{\displaystyle \frac{a_2}{{({\lambda }_2-\lambda )}^2}}+\cdots +{\displaystyle \frac{a_k}{{({\lambda }_k-\lambda )}^2}}>0$ where $i=1,2,\dots ,k$. Therefore, we have that $r_i(\lambda )$ are monotonically increasing functions with respect to parameter λ on the intervals $(-\infty ,{\lambda }_1),({\lambda }_i,{\lambda }_{i+1}),({\lambda }_k,+\infty )$ where $i=1,2,\dots ,k-1$.

For the REP (1) and (2), in this paper we only discuss the eigenvalue distribution on the positive semi-real axis, namely $J\subseteq {\mathcal{R}}^{+}$. Assume that all zeros of $q_i(x)$ with $i=1,2,\dots ,N$ on the real positive semiaxis are arranged in the following order:

$$infJ\triangleq {\sigma }_0<{\sigma }_1\le \dots \le {\sigma }_L<{\sigma }_{L+1}<\stackrel{\Delta }{=}supJ.$$

Set $J_l=\left({\sigma }_{l-1},{\sigma }_l\right)$ where $l=1,\dots ,L+1$. Then we have $J={\displaystyle \bigcup _{l=1}^{L+1}}{J}_l={\mathcal{R}}^{+}\backslash {\left\{{\sigma }_i\right\}}_{i=1}^L$ and $J_l\cap J_k=\varphi$ if $l\ne k$.

As long as the matrices $A_i$ and $B_j$ are symmetrical for all $i=0,1,\dots ,d$ and $j=1,2,\dots ,N$, we can define the Rayleigh functional $p(x)$ for the REP. That is to say, if $p(x)$ satisfies the following equation

$$\langle R\left(p\left(x\right)\right)x,x\rangle =0,$$

(3)

$p(x)$ is the Rayleigh functional of $R(\lambda )$. Notice that in the linear case $R(\lambda )=\lambda I-A$, it is exactly the Rayleigh quotient. Let $f\left(\lambda ,x\right)=\langle R\left(\lambda \right)x,x\rangle$, then $p\left(x\right)$ is a root of $f\left(\lambda ,x\right)$ = 0. Because $A_i$ is positive definite with $i=1,\dots ,d$, $B_i$ is positive semidefinite and $r_i\left(\lambda \right)$ is monotonically increasing function, it is easy to verify that $f\left(\lambda ,x\right)$ is a monotone increasing function on the intervals separated by the zeros of polynomials $q_i\left(\lambda \right)$. Therefore,

$$\left(\lambda -p\left(x\right)\right)f\left(\lambda ,x\right)>0,$$

(4)

where $\lambda \ne p\left(x\right)$.

For each fixed $\lambda \in J$, we consider the following standard eigenvalue problem (SEP):

$$R\left(\lambda \right)x=\mu x.$$

(5)

We can easily see that if $\lambda$ is an eigenvalue of the REP (1) and (2), $\mu =0$ is an eigenvalue of the above SEP (5). Conversely, it is also true. Therefore, the eigenvalue of the REP (1) and (2) can be characterized by the zero eigenvalue of the SEP (5). For the standard eigenvalue problem, we have the minmax principle

$${\mu }_j=\underset{V\in H_j}{sup}\underset{x\in V\backslash \left\{0\right\}}{min}{\displaystyle \frac{\langle R\left(\lambda \right)x,x\rangle }{\langle x,x\rangle }},$$

where $H_j$ represents the set of Hilbert subspaces with dimension j of ${\mathcal{R}}^n$. Similarly, we have the minmax principle of the REP

$${\lambda }_j=\underset{V\in H_j,V\cap \mathcal{D}\ne \varphi }{min}\underset{x\in V\cap \mathcal{D}}{sup}p\left(x\right),$$

(6)

where $\mathcal{D}$ denotes the domain of the Rayleigh functional $p\left(x\right)$ which satisfies (3). For a more detailed discussion, see, e.g., [17,18].

3. Spectral Properties and Distribution of the REPs

According to the assumption (H2), we know that the matrices $B_i$ are low-rank. Hence, the REP (2) can be regarded as a low-rank perturbation of the following PEP [25]

$$\tilde{P}\left(\lambda \right)x=0,$$

(7)

where

$$\tilde{P}\left(\lambda \right)=P\left(\lambda \right)+\sum _{i=1}^N{\displaystyle \frac{s_i\left(\kappa \right)}{q_i\left(\kappa \right)}}{B}_i,$$

with $\kappa \in J_l$ and $l=1,2,\dots ,L+1$.

For any $l\in \left\{1,2,\dots ,L+1\right\}$ and any value $\kappa$ in the interval $J_l$, the REP (1) and (2) has the corresponding PEP (7). Assume that the PEP (7) has an eigenvalue $\mu$ in the interval $J_l$. Then we will prove that the eigenvalue $\lambda$ of the REP (1) and (2) is between $\kappa$ and $\mu$. Before giving this theorem, we first show some lemmas.

Lemma 1.

Assume that $s\left(\lambda \right)$ and $q\left(\lambda \right)$ represent any pair of scalar polynomials $s_i\left(\lambda \right)$ and $q_i\left(\lambda \right)$ with $i=1,2,\dots ,N$, respectively. For any ${\lambda }_1,{\lambda }_2\in J_l$ with $l=1,2,\dots ,L+1$, we have

$$s\left({\lambda }_1\right)q\left({\lambda }_2\right)-q\left({\lambda }_1\right)s\left({\lambda }_2\right)=({\lambda }_1-{\lambda }_2)g({\lambda }_1,{\lambda }_2),$$

where $g\left({\lambda }_1,{\lambda }_2\right)>0$ with $a_i$ and $b_j$ standing for the coefficients of $s(\lambda )$ and $q(\lambda )$, respectively.

Proof. 

Without losing generality, we let $s\left(\lambda \right)={\displaystyle \sum _{i=0}^m}{a}_i{\lambda }^i$ and $q\left(\lambda \right)={\displaystyle \sum _{j=0}^n}{b}_j{\lambda }^j$. Hence,

$$s\left({\lambda }_1\right)q\left({\lambda }_2\right)-q\left({\lambda }_1\right)s\left({\lambda }_2\right)=\left(\sum _{i=0}^m{a}_i{\lambda }_1^i\right)\left(\sum _{j=0}^n{b}_j{\lambda }_2^j\right)-\left(\sum _{i=0}^m{a}_i{\lambda }_2^i\right)\left(\sum _{j=0}^n{b}_j{\lambda }_1^j\right)=\sum _{i=0}^m\sum _{j=0}^n{a}_i{b}_j\left({\lambda }_1^i{\lambda }_2^j-{\lambda }_1^j{\lambda }_2^i\right).$$

Let

$$g({\lambda }_1,{\lambda }_2)=\sum _{i\ne j}sgn(i-j)a_i{b}_j{\left({\lambda }_1{\lambda }_2\right)}^{min(i,j)}\sum _{k=0}^{\left|i-j\right|-1}{\lambda }_1^k{\lambda }_2^{\left|i-j\right|-k-1},$$

and we have $s\left({\lambda }_1\right)q\left({\lambda }_2\right)-q\left({\lambda }_1\right)s\left({\lambda }_2\right)=\left({\lambda }_1-{\lambda }_2\right)g\left({\lambda }_1,{\lambda }_2\right)$.

Let $r\left(\lambda \right)={\displaystyle \frac{s\left(\lambda \right)}{q\left(\lambda \right)}}$. Then

$$r\left({\lambda }_1\right)-r\left({\lambda }_2\right)=\frac{s\left({\lambda }_1\right)}{q\left({\lambda }_1\right)}-\frac{s\left({\lambda }_2\right)}{q\left({\lambda }_2\right)}=\frac{s\left({\lambda }_1\right)q\left({\lambda }_2\right)-q\left({\lambda }_1\right)s\left({\lambda }_2\right)}{q\left({\lambda }_1\right)q\left({\lambda }_2\right)}=\frac{{\lambda }_1-{\lambda }_2}{q\left({\lambda }_1\right)q\left({\lambda }_2\right)}g\left({\lambda }_1,{\lambda }_2\right).$$

Because ${\lambda }_1,{\lambda }_2\in J_l$, then $q\left({\lambda }_1\right)q\left({\lambda }_2\right)>0$. Based on the assumption $\mathbf{H}\mathbf{3}$ we have $g\left({\lambda }_1,{\lambda }_2\right)>0$. Thus, the conclusion holds true.    □

For the PEP (7), the Rayleigh functional $R_{\kappa }\left(x\right)$ should satisfy $\langle \tilde{P}\left(R_{\kappa }\left(x\right)\right)x,x\rangle =0$, namely,

$$\langle P\left(R_{\kappa }\left(x\right)\right)x,x\rangle =-\sum _{i=1}^N{\displaystyle \frac{s_i\left(\kappa \right)}{q_i\left(\kappa \right)}}\langle B_{i}x,x\rangle .$$

(8)

Lemma 2.

Assume that $\kappa \in J_l$ and there exists a vector $x\in {\mathcal{H}}^1$, such that $R_{\kappa }\left(x\right)\in J_l$ with $l=1,2,\dots ,L+1$ where ${\mathcal{H}}^1=\left\{\left.x\in \mathcal{H}\right|\langle x,x\rangle =1\right\}$. Then $x\in {\mathcal{D}}_l$ and

$$min\left(\kappa ,R_{\kappa }\left(x\right)\right)\le p_l\left(x\right)\le max\left(\kappa ,R_{\kappa }\left(x\right)\right),$$

where ${\mathcal{D}}_l$ is the domain of the Rayleigh functional $p_l\left(x\right)$ which satisfies (3) in $J_l$.

Proof. 

With the fact that

$$f\left(R_{\kappa }\left(x\right),x\right)=\langle P\left(R_{\kappa }\left(x\right)\right)x,x\rangle +\sum _{i=1}^N{\displaystyle \frac{s_i\left(R_{\kappa }\left(x\right)\right)}{q_i\left(R_{\kappa }\left(x\right)\right)}}\langle B_{i}x,x\rangle$$

and the relation (8), we have

$$\begin{array}{cc}\hfill f\left(R_{\kappa }\left(x\right),x\right)& =-{\displaystyle \sum _{i=1}^N}\frac{s_i\left(\kappa \right)}{q_i\left(\kappa \right)}\langle B_{i}x,x\rangle +{\displaystyle \sum _{i=1}^N}{\displaystyle \frac{s_i\left(R_{\kappa }\left(x\right)\right)}{q_i\left(R_{\kappa }\left(x\right)\right)}}\langle B_{i}x,x\rangle \hfill \\ \hfill & ={\displaystyle \sum _{i=1}^N}\frac{\left(R_{\kappa }\left(x\right)-\kappa \right)g_i\left(\kappa ,R_{\kappa }\left(x\right)\right)}{q_i\left(\kappa \right)q_i\left(R_{\kappa }\left(x\right)\right)}\langle B_{i}x,x\rangle .\hfill \end{array}$$

It follows from Lemma 1 and the assumption (H2) that $f\left(R_{\kappa }\left(x\right),x\right)\le 0$ if $\kappa \ge R_{\kappa }\left(x\right)$ and $f\left(R_{\kappa }\left(x\right),x\right)\ge 0$ if $\kappa \le R_{\kappa }\left(x\right)$. Similarly, it can be proved that $f\left(\kappa ,x\right)\ge 0$ if $\kappa \ge R_{\kappa }\left(x\right)$ and $f\left(\kappa ,x\right)\le 0$ if $\kappa \le R_{\kappa }\left(x\right)$. Finally, because $f\left(\lambda ,x\right)$ is continuous, we have $x\in {\mathcal{D}}_l$ and $min\left(\kappa ,R_{\kappa }\left(x\right)\right)\le p_l\left(x\right)\le max\left(\kappa ,R_{\kappa }\left(x\right)\right)$, which completes the proof.    □

Theorem 1.

Assume that (H1)–(H3) hold, and let $\kappa \in J_l$ with $l=1,2,\dots ,L+1$. Assume that ${\mu }_j\in J_l$ and $J_l$ contains the j-th eigenvalue of the PEP (7). Then the REP (1) and (2) has a corresponding eigenvalue ${\lambda }_j\in J_l$ and $min(\kappa ,{\mu }_j)\le {\lambda }_j\le max(\kappa ,{\mu }_j)$.

Proof. 

We first show that there exists a subspace $V\in {\mathcal{H}}_j$, such that

$$V\cap {\mathcal{D}}_l\ne ⌀\mathrm{and}\underset{x\in V\cap {\mathcal{D}}_l}{sup}{p}_l\left(x\right)\le max(\kappa ,{\mu }_j).$$

In fact, suppose that there exist $W\in {\mathcal{H}}_j$ and $w\in W\backslash \left\{0\right\}$, such that

$${\mu }_j=\underset{V\in H_j,V\cap {\mathcal{D}}_l\ne \varphi }{min}\underset{x\in V\cap {\mathcal{D}}_l}{max}{R}_{\kappa }\left(x\right)=\underset{x\in W\backslash \left\{0\right\}}{max}{R}_{\kappa }\left(x\right)=R_{\kappa }\left(w\right).$$

Because $\kappa \in J_l$ and ${\mu }_j\in J_l$, we have $w\in {\mathcal{D}}_l$ from Lemma 2. Then $W\cap {\mathcal{D}}_l\ne ⌀$. For any $x\in W\backslash \left\{0\right\}$, we have $R_{\kappa }\left(x\right)\le {\mu }_j$. Therefore, $\langle P\left(R_{\kappa }\left(x\right)\right)x,x\rangle \le \langle P({\mu }_j)x,x\rangle$. That is,

$$\langle P({\mu }_j)x,x\rangle +\sum _{i=1}^N{\displaystyle \frac{s_i\left(\kappa \right)}{q_i\left(\kappa \right)}}\langle B_{i}x,x\rangle \ge 0.$$

Set $\delta =max(\kappa ,{\mu }_j)$. Then it is easy to obtain that $\langle P\left(\delta \right)x,x\rangle +{\displaystyle \sum _{i=1}^N}{\displaystyle \frac{s_i\left(\delta \right)}{q_i\left(\delta \right)}}\langle B_{i}x,x\rangle \ge 0$. That is, $f\left(\delta ,x\right)\ge 0$. It follows from (4) that $p_l\left(x\right)\le \delta$ for any $x\in {\mathcal{D}}_l\cap W$. Hence, $\underset{x\in W\cap D_l}{sup}{p}_l\left(x\right)\le \delta =max(\kappa ,{\mu }_j)$.

In the following, we show that for any $V\in {\mathcal{H}}_j$, if $V\cap \mathcal{D}\ne ⌀$, we have

$$\underset{x\in V\cap {\mathcal{D}}_l}{sup}{p}_l\left(x\right)\ge min(\kappa ,{\mu }_j).$$

We prove this result by reduction to absurdity. Suppose that there exists $V\in {\mathcal{H}}_j$ such that $V\cap \mathcal{D}\ne ⌀$, but $\underset{x\in V\cap {\mathcal{D}}_l}{sup}{p}_l(x)<min(\kappa ,{\mu }_j)$. Let $x_V\in V$ such that $R_{\kappa }\left(x_V\right)=\underset{x\in V}{max}{R}_{\kappa }\left(x\right)$. Then we have $R_{\kappa }\left(x_V\right)\notin J_l$. Actually, if $R_{\kappa }\left(x_V\right)\in J_l$, it is easy to get $x_V\in {\mathcal{D}}_l$ and $p_l\left(x_V\right)\ge min\left(\kappa ,R_{\kappa }\left(x_V\right)\right)$.

That is,

$$\underset{x\in V\cap {\mathcal{D}}_l}{sup}{p}_l\left(x\right)\ge p_l\left(x_V\right)\ge min\left(\kappa ,R_{\kappa }\left(x_V\right)\right)\ge min(\kappa ,{\mu }_j)$$

contradicting the fact that $\underset{x\in V\cap {\mathcal{D}}_l}{sup}{p}_l\left(x\right)<min(\kappa ,{\mu }_j)$. Hence, $R_{\kappa }\left(x_V\right)\notin J_l$.

Set $\delta =min(\kappa ,{\mu }_j)\le min\left(\kappa ,R_{\kappa }\left(x_V\right)\right)$. Then we have $\delta \le {\mu }_j\le R_{\kappa }\left(x_V\right)$. Moreover, because $f\left(\delta ,x_V\right)=\langle P\left(\delta \right)x_V,x_V\rangle +{\displaystyle \sum _{i=1}^N}{\displaystyle \frac{s_i\left(\delta \right)}{q_i\left(\delta \right)}}\langle B_i{x}_V,x_V\rangle$, it is easy to obtain

$$f\left(\delta ,x_V\right)\le \langle P\left(R_{\kappa }\left(x_V\right)\right)x_V,x_V\rangle -\sum _{i=1}^K\frac{s_i\left(\kappa \right)}{q_i\left(\kappa \right)}\langle B_i{x}_V,x_V\rangle =0.$$

For any $x\in V\cap {\mathcal{D}}_l$, we let $w\left(t\right)=tx+\left(1-t\right)x_V$ and $\phi \left(t\right)=f\left(\delta ,w\left(t\right)\right)$. Then $\phi \left(0\right)=f\left(\delta ,x_V\right)\le 0$. Because $p_l\left(x\right)<\delta$, it follows from (4) that $\phi \left(1\right)=f\left(\delta ,x\right)>f\left(p_l\left(x\right),x\right)=0$. There exists $\tilde{t}\in [0,1)$ such that $f\left(\delta ,w\left(\tilde{t}\right)\right)=0$. Thus, $w\left(\tilde{t}\right)\in V\cap {\mathcal{D}}_l$ and $p_l\left(w\left(\tilde{t}\right)\right)=\delta =min(\kappa ,{\mu }_j)$ which conflicts with the assumption.

To summarise, we have

$${\lambda }_j=\underset{V\in {\mathcal{H}}_j,V\cap {\mathcal{D}}_l\ne \varphi }{inf}\underset{x\in V\cap {\mathcal{D}}_l}{sup}{p}_l\left(v\right)\in J_l,\mathrm{and}min\left(\kappa ,{\mu }_j\right)\le {\lambda }_j\le max\left(\kappa ,{\mu }_j\right),$$

which completes the proof.    □

Actually, eigenvalue ${\mu }_j$ of the REP (1) and (2) is the function with respect to $\kappa$, which implies that ${\mu }_j={\mu }_j\left(\kappa \right)$. The following theorem will elaborate the continuity of the function ${\mu }_j\left(\kappa \right)$.

Theorem 2.

Assume that H1–H3 hold, then ${\mu }_j\left(\kappa \right)$ is a continuous decreasing function with respect to κ.

Proof. 

Let $\tilde{P}\left(\lambda ,\kappa \right)\stackrel{\Delta }{=}\tilde{P}\left(\lambda \right)$. For any $\epsilon >0$, there exists $\delta >0$ such that

$$∥\tilde{P}\left(\lambda ,{\kappa }_1\right)-\tilde{P}\left(\lambda ,{\kappa }_2\right)∥\le \sum _{i=1}^N\left|\frac{s_i\left({\kappa }_1\right)}{q_i\left({\kappa }_1\right)}-\frac{s_i\left({\kappa }_2\right)}{q_i\left({\kappa }_2\right)}\right|∥B_i∥=\left|{\kappa }_1-{\kappa }_2\right|\sum _{i=1}^N\left|\frac{g\left({\kappa }_1,{\kappa }_2\right)}{q_i\left({\kappa }_1\right)q_i\left({\kappa }_2\right)}\right|∥B_i∥<\epsilon$$

when $\left|{\kappa }_1-{\kappa }_2\right|<\delta$ with $\delta =\delta \left(\epsilon \right)>0$ small enough. Because the eigenvalue is a continuous function with respect to the elements of its matrix [26], we have that ${\mu }_j\left(\kappa \right)$ is a continuous function with respect to $\kappa$.

It follows from (8) that $R_{\kappa }\left(x\right)$ is the root of the following polynomial equation

$${\lambda }^d+a_{d-1}{\lambda }^{d-1}+\cdots +a_1\lambda +a_0=0,$$

where $a_i=\langle A_{i}x,x\rangle >0$ with $i=1,2,\dots ,d-1$ and $a_0=\langle A_{0}x,x\rangle +{\displaystyle \sum _{i=1}^N}{\displaystyle \frac{s_i\left(\kappa \right)}{q_i\left(\kappa \right)}}\langle B_{i}x,x\rangle .$ With fixed x, $a_i$ remains unchanged where $i=1,2,\dots ,d-1$. Moreover, from the assumption (H3), we know that $a_0$ is an increasing function of $\kappa$. Then we can easily prove that $R_{\kappa }\left(x\right)$ is a decreasing function of $\kappa$. Finally, through the minmax principle (6), we can conclude that ${\mu }_j\left(\kappa \right)$ is a continuous decreasing function of $\kappa$, which completes the proof.    □

Theorems 1 and 2 show that if the REP (1) and (2) have an eigenvalue ${\lambda }_j\in J_l$, there must be such a value $\kappa$ and one eigenvalue ${\mu }_j$ of the PEP (7) in the interval $J_l$. Conversely, if the REP (1) and (2) have no eigenvalues in the interval $J_l$, the PEP (7) will not have any eigenvalues in this interval even if the values $\kappa$ in the interval $J_l$ are taken all over. Therefore, there exists a one-to-many relationship between ${\lambda }_j$ and ${\mu }_j$.

On the other hand, suppose that there exists ${\kappa }_0\in J_l$, then the PEP (7) has two unequal eigenvalues ${\mu }_j\left({\kappa }_0\right)\ne {\tilde{\mu }}_j\left({\kappa }_0\right)\in J_l$. If ${\lambda }_j$ and ${\tilde{\lambda }}_j$ are fixed points of ${\mu }_j\left({\kappa }_0\right)$ and ${\tilde{\mu }}_j\left({\kappa }_0\right)$, respectively, we have ${\lambda }_j\ne {\tilde{\lambda }}_j$. In fact, because the eigenvalue is a continuous function with respect to the elements of its matrix, the multiplicity of the original eigenvalue will not change with the change of $\kappa$. That is, ${\mu }_j\left(\kappa \right)\ne {\tilde{\mu }}_j\left(\kappa \right)$ holds. Note that if $\kappa ={\lambda }_j$, we have ${\lambda }_j\ne {\tilde{\lambda }}_j$. Therefore, there is one-to-one correlation between ${\lambda }_j$ and ${\mu }_j$.

To summarise, we can obtain the following theorems.

Theorem 3.

There is a one-to-one relationship between the eigenvalue ${\lambda }_j$ of the REP (1) and (2) and the eigenvalue ${\mu }_j$ of the PEP (7) in $J_l$ with $l=1,2,\dots ,L+1$.

Theorem 4.

Assume that ${\sigma }_{l-1}<{\kappa }_1<{\kappa }_2<{\sigma }_l$. Let $N\left(\kappa \right)=max\left\{n:{\mu }_n\left(\kappa \right)\le \kappa \right\}$, then there are exactly $N\left({\kappa }_2\right)-N\left({\kappa }_1\right)$ eigenvalues in the semi-interval $({\kappa }_1,{\kappa }_2]$.

4. Numerical Methods for Solving the REPs

In this section, based on the above spectral distribution of the REP we discuss the numerical methods for solving the REP (1) and (2). Given a $\kappa \in J_l$, we can find an eigenvalue $\mu$ of the PEP (7). Here, how to select the next new value $\kappa$ is the key to propose the novel numerical algorithms. Because $\mu \left(\kappa \right)$ is a continuous decreasing function of $\kappa$ and $\lambda$ is a fixed point of $\mu \left(\kappa \right)$, we can choose the newest $\kappa$ by a certain fixed-point algorithm. For simplicity, we first consider to choose $\kappa$ via dichotomy as follows:

$$\kappa :=\frac{\kappa +\mu }{2}.$$

(9)

Therefore, we derive the following numerical method (Algorithm 1) for solving the REP (1) and (2).   

Algorithm 1: Dichotomy iteration method for the REP (1) and (2)

Input: rational matrix function $R\left(\lambda \right)$, the target point $\sigma \in J_l$ and the tolerance $\tau$. Output: the approximate eigenvalue $\lambda \in J_l$ closest to $\sigma$. Set $\kappa =\sigma \in J_l$; Construct the polynomial eigenvalue problem (7); Solve the PEP (7) to get eigenpair $\left(\mu ,x\right)$ where $\mu$ is closest to $\sigma$. If no such $\mu$ is present, then stop; otherwise; Compute $rn=∥R\left(\mu \right)x∥/\phantom{rn=∥R\left(\mu \right)x∥\left({\displaystyle \sum _{i=0}^d}{\left|\mu \right|}^i{∥A_i∥}_F+{\displaystyle \sum _{i=1}^N}\left|\frac{s_i\left(\mu \right)}{q_i\left(\mu \right)}\right|{∥E_i∥}_F\right)}\left({\displaystyle \sum _{i=0}^d}{\left|\mu \right|}^i{∥A_i∥}_F+{\displaystyle \sum _{i=1}^N}\left|{\displaystyle \frac{s_i\left(\mu \right)}{q_i\left(\mu \right)}}\right|{∥B_i∥}_F\right)$. If $rn<\tau$, then stop, and output $\lambda =\mu$; otherwise; Update $\kappa$ by using (9) and go to step 2.

Remark 2.

In actual computation, for the small-scale PEP (7) the classical approach is to turn it into a generalized eigenvalue problem (GEP) via linearization, or solve it directly by the in-built function of Matlab. For the large-scale ones, we can adopt the partially orthogonal projection method [27] to solve it.

Remark 3.

Assume that x is an eigenvector of μ for the PEP (7). If $\lambda =\mu$, x is also the eigenvector corresponding to λ of the REP (1) and (2). Therefore, we can use the Rayleigh functional to accelerate κ as follows:

$$\kappa :=p_l\left(x\right).$$

(10)

Thus another numerical method (Algorithm 2) for solving the REP (1) and (2) can be summarized as follows.

Algorithm 2: Rayleigh functional iteration method for the REP (1) and (2)

Input: rational matrix function $R\left(\lambda \right)$, the target point $\sigma \in J_l$ and the tolerance $\tau$. Output: the approximate eigenvalue $\lambda \in J_l$ closest to $\sigma$. Set $\kappa =\sigma \in J_l$; Construct the polynomial eigenvalue problem (7); Solve the PEP (7) to get eigenpair $\left(\mu ,x\right)$ where $\mu$ is closest to $\sigma$. If no such $\mu$ is present, then stop; otherwise; Compute $rn=∥R\left(\mu \right)x∥/\phantom{rn=∥R\left(\mu \right)x∥\left({\displaystyle \sum _{i=0}^d}{\left|\mu \right|}^i{∥A_i∥}_F+{\displaystyle \sum _{i=1}^N}\left|\frac{s_i\left(\mu \right)}{q_i\left(\mu \right)}\right|{∥E_i∥}_F\right)}\left({\displaystyle \sum _{i=0}^d}{\left|\mu \right|}^i{∥A_i∥}_F+{\displaystyle \sum _{i=1}^N}\left|{\displaystyle \frac{s_i\left(\mu \right)}{q_i\left(\mu \right)}}\right|{∥B_i∥}_F\right)$. If $rn<tol$, then stop, and output $\lambda =\mu$; otherwise; Update $\kappa$ by using (10) and go to step 2.

5. Numerical Results

In this section, we report some numerical examples to show the effectiveness of the proposed Algorithms 1 and 2. All computations are performed under Matlab (version R2019a). In our examples, ${\lambda }_{\ast }$ is an exact eigenvalue of the REP (1) and (2), and $\lambda$ is an approximate eigenvalue computed by Algorithm 1 or Algorithm 2. CPU denotes the CPU time (in seconds) for computing an approximate solution, and Iter denotes the number of iteration steps. The stopping tolerance for the residual norm is chosen to be $\tau ={10}^{-12}$.

Example 2

([9]). We consider the following REP:

$$R\left(\lambda \right)x=\left(\lambda B-A-E-\frac{\sigma }{\lambda -\sigma }E\right)x=0,$$

(11)

where A and B are the positive definite tridiagonal matrices defined as

$$A=\frac{1}{h}\left(\begin{array}{cccc}2& -1& & \\ -1& \ddots & \ddots & \\ & \ddots & 2& -1\\ & & -1& 1\end{array}\right),B=\frac{h}{6}\left(\begin{array}{cccc}4& 1& & \\ 1& \ddots & \ddots & \\ & \ddots & 4& 1\\ & & 1& 2\end{array}\right)$$

with $h={\displaystyle \frac{1}{n}}$ and $E=e_n{e}_n^T$.

We can easily check that the above eigenvalue problem meets the assumptions H1–H3 in Section 2. Let $n=1000$ and $\sigma =1$. Here we are interested in computing all eigenvalues of the REP (11) in the interval $J=\left(0,100\right)$, and we can divide it into two small intervals such as $J_1=\left(0,1\right)$ and $J_2=\left(1,100\right)$.

It is easy to verify that ${\lambda }_{\ast }=0.45731832$, $4.48202582$, $24.21875011$ and $63.69036457$ are the eigenvalues of the REP (11) in the interval $\left(0,100\right)$ because the corresponding smallest singular values of $R({\lambda }_{\ast })$ are less than ${10}^{-9}$.

Through Theorem 4, we have that the numbers of eigenvalues of the REP (11) in $J_1$ and $J_2$ are 1, and 3, respectively. The above result is completely consistent with the actual distribution of eigenvalues for the REP (11).

In the following, we choose different $\kappa$ in $J_1$ and $J_2$ such as $\kappa =0.5,13,38,75$ and apply Algorithms 1 and 2 to compute all eigenvalues of the (11) in $J_1$ and $J_2$.

The numerical results for Algorithms 1 and 2 are reported in

Table 1. Numerical results of Example 2.

$\mathit{\kappa }=0.5$	Iter	λ	rn	CPU
Algorithm 1	7	0.457318355879899	$6.60\times {10}^{-13}$	0.57
Algorithm 2	3	0.457318325580995	$1.10\times {10}^{-15}$	0.29
$\kappa =13$
Algorithm 1	21	4.482026373417336	$8.36\times {10}^{-13}$	1.69
Algorithm 2	3	4.482025808828371	$1.19\times {10}^{-13}$	0.29
$\kappa =38$
Algorithm 1	23	24.218751627286714	$6.27\times {10}^{-13}$	1.89
Algorithm 2	2	24.218750015378188	$8.46\times {10}^{-13}$	0.18
$\kappa =75$
Algorithm 1	23	63.690365902995510	$5.45\times {10}^{-13}$	1.85
Algorithm 2	2	63.690364569627210	$1.34\times {10}^{-15}$	0.18

, which shows that the proposed methods are very useful and efficient to solve rational eigenvalue problems in one interval. Moreover, Algorithm 2 requires less CPU and iteration steps than Algorithm 1. Moreover, the numerical results remain the same when n and $\sigma$ take the other different values.

6. Conclusions

In this paper, the spectral distribution of one class of rational eigenvalue problems has been studied in detail, and two simple iterative methods for solving this kind of rational eigenvalue problems have been proposed based on the spectral distribution. Numerical examples show the efficiency of the new approaches. The spectral properties and distribution of the general rational eigenvalue problems are remaining for our future work.