Regularized Normalization Methods for Solving Linear and Nonlinear Eigenvalue Problems

Abstract

To solve linear and nonlinear eigenvalue problems, we develop a simple method by directly solving a nonhomogeneous system obtained by supplementing a normalization condition on the eigen-equation for the uniqueness of the eigenvector. The novelty of the present paper is that we transform the original homogeneous eigen-equation to a nonhomogeneous eigen-equation by a normalization technique and the introduction of a simple merit function, the minimum of which leads to a precise eigenvalue. For complex eigenvalue problems, two normalization equations are derived utilizing two different normalization conditions. The golden section search algorithms are employed to minimize the merit functions to locate real and complex eigenvalues, and simultaneously, we can obtain precise eigenvectors to satisfy the eigen-equation. Two regularized normalization methods can accelerate the convergence speed for two extensions of the simple method, and a derivative-free fixed-point Newton iterative scheme is developed to compute real eigenvalues, the convergence speed of which is ten times faster than the golden section search algorithm. Newton methods are developed for solving two systems of nonlinear regularized equations, and the efficiency and accuracy are significantly improved. Over ten examples demonstrate the high performance of the proposed methods. Among them, the two regularization methods are better than the simple method.

1. Introduction

It is well known that the Rayleigh quotient [1,2]:

$$R\left(\mathbf{x}\right):=\frac{{\mathbf{x}}^{\mathrm{T}}\mathbf{A}\mathbf{x}}{{\mathbf{x}}^{\mathrm{T}}\mathbf{x}}$$

(1)

can be used to determine the real eigenvalue of a symmetric matrix $\mathbf{A}$.

In this paper, we derive a simple normalized condition solver to obtain the eigenvalues of

$$\mathbf{A}\mathbf{x}=\lambda \mathbf{x},$$

(2)

where $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ is a given matrix, $\mathbf{x}\in {\mathbb{R}}^n$ is an unknown vector, and $\lambda$ is an unknown eigenvalue in the standard linear eigen-equation. When $\mathbf{A}$ is not a symmetric matrix, the Rayleigh quotient (1) cannot be used to determine the eigenvalues. Liu et al. [3] developed a new quotient to determine the eigenvalues of Equation (2).

As noticed by Liu et al. [4], it is hard to directly determine the eigenvalue and eigenvector from Equation (2) by a numerical method. In fact, from $\mathbf{A}\mathbf{x}-\lambda \mathbf{x}=\mathbf{0}$, we always obtain $\mathbf{x}=\mathbf{0}$ by a numerical method since the right-hand side is a zero vector. In [4], a new strategy to overcome this difficulty uses the variable transformation to a new nonhomogeneous linear system. It possesses a nonzero external excitation term on the right-hand side, such that one can obtain a nonzero eigenvector when the eigen-parameter $\lambda$ is an eigenvalue. We are going to propose a simple method to nonhomogenize the eigen-equation to obtain a nonhomogeneous linear system, and then it is easy to find the eigenvalue and eigenvector by using the minimization technique.

The standard free vibration model of elastic structural elements is

$$\mathbf{M}\ddot{\mathbf{q}}\left(t\right)+\mathbf{C}\dot{\mathbf{q}}\left(t\right)+\mathbf{K}\mathbf{q}\left(t\right)=\mathbf{0},$$

(3)

which by $\mathbf{q}\left(t\right)=e^{\lambda t}\mathbf{x}$ renders a quadratically nonlinear eigenvalue problem [5]:

$$({\lambda }^2\mathbf{M}+\lambda \mathbf{C}+\mathbf{K})\mathbf{x}=\mathbf{0}.$$

(4)

A lot of applications and solvers of quadratic eigenvalue problems have been proposed, e.g., stability analysis of time-delay systems [6], free vibrations of fluid–solids structures [7], a modified second-order Arnoldi method [8], the inexact residual iteration method [9], the homotopy perturbation technique [10], electromagnetic wave propagation and analysis of an acoustic fluid contained in a cavity with absorbing walls [11], and a friction-induced vibration problem under variability [12]. In addition, several applications and solvers of generalized eigenvalue problems have been addressed, e.g., the block Arnoldi-type contour integral spectral projection method [13], small-sample statistical condition estimation [14], matrix perturbation methods [15], the overlapping finite element method [16], the complex HZ method [17], the context of sensor selection [18], and a generalized Arnoldi method [19].

As done in [4], we can take

$$\mathbf{y}=\lambda \mathbf{x},$$

(5)

and combine Equations (5) and (4) to obtain

$$\left[\begin{array}{cc}{\mathbf{0}}_n& {\mathbf{I}}_n\\ -\mathbf{K}& -\mathbf{C}\end{array}\right]\left[\begin{array}{c}\mathbf{x}\\ \mathbf{y}\end{array}\right]=\lambda \left[\begin{array}{cc}{\mathbf{I}}_n& {\mathbf{0}}_n\\ {\mathbf{0}}_n& \mathbf{M}\end{array}\right]\left[\begin{array}{c}\mathbf{x}\\ \mathbf{y}\end{array}\right].$$

(6)

Upon defining

$$\mathbf{X}:=\left[\begin{array}{c}\mathbf{x}\\ \mathbf{y}\end{array}\right],\mathbf{A}:=\left[\begin{array}{cc}{\mathbf{0}}_n& {\mathbf{I}}_n\\ -\mathbf{K}& -\mathbf{C}\end{array}\right],\mathbf{B}:=\left[\begin{array}{cc}{\mathbf{I}}_n& {\mathbf{0}}_n\\ {\mathbf{0}}_n& \mathbf{M}\end{array}\right],$$

(7)

Equation (6) becomes a generalized eigenvalue problem for the $2n$-vector $\mathbf{X}$:

$$\mathbf{A}\mathbf{X}=\lambda \mathbf{B}\mathbf{X},$$

(8)

where $\mathbf{A},\mathbf{B}\in {\mathbb{R}}^{2n\times 2n}$. Equation (8) is used to determine the eigen-pair $(\lambda ,\mathbf{X})$, which is a linear eigen-equation associated with the pencil $\mathbf{A}-\lambda \mathbf{B}$, where $\lambda$ is an eigen-parameter. A main drawback of this argumentation is that the dimension is raised doubly from n to $2n$.

Equation (8) can be written as

$$(\mathbf{A}-\lambda \mathbf{B})\mathbf{X}=\mathbf{0}.$$

(9)

Because the right-hand side is a zero vector, solving it by the numerical method we can only obtain the trivial solution $\mathbf{X}=\mathbf{0}$. To avoid this situation, Liu et al. [4] introduced an external excitation method by letting $\mathbf{X}=\mathbf{Y}+{\mathbf{X}}_0$, such that

$$(\mathbf{A}-\lambda \mathbf{B})\mathbf{Y}=-(\mathbf{A}-\lambda \mathbf{B}){\mathbf{X}}_0.$$

(10)

Solving this equation for $\mathbf{Y}$, then the eigenvector $\mathbf{X}=\mathbf{Y}+{\mathbf{X}}_0$ is obtained. However, it is a problem how to select the proper exciting vector ${\mathbf{X}}_0$. The basic idea is the transformation from a homogeneous Equation (9) to a nonhomogeneous Equation (10). We need to determine if we can develop a simpler method to realize such a type transformation but without introducing an extra exciting vector ${\mathbf{X}}_0$, which is an interesting problem. The present paper attempts to make this type transformation very easy, which is the main motivation and the major novelty: to realize this transformation by a simple normalization technique. The present idea is simpler than that in [4], so the new technique is called a simple method and is introduced in Section 2.

Nonlinear eigenvalue problems are important and find a lot of real applications in engineering and applied fields. Betcke et al. [20] collected 52 nonlinear eigenvalue problems in the form of a MATLAB toolbox, which contains problems from models of real-life applications as well as ones constructed specifically to have particular properties. Recently, El-Guide et al. [21] presented two approximation methods for computing eigenfrequencies and eigenmodes of large-scale nonlinear eigenvalue problems resulting from boundary element method solutions of some types of acoustic eigenvalue problems in 3D space. We extend Equation (4) to a general nonlinear eigenvalue problem [20]:

$$\mathbf{N}\left(\lambda \right)\mathbf{x}=\mathbf{0},$$

(11)

which is a nonlinear eigen-equation of $\lambda$ used to solve the eigen-pair $(\lambda ,\mathbf{x})$, where $\mathbf{N}\left(\lambda \right)\in {\mathbb{R}}^{n\times n}$. Equation (11) is a nonlinear eigenvalue problem because $\mathbf{N}\left(\lambda \right)$ is a nonlinear matrix function of the eigen-parameter $\lambda$. In Equation (9), $\mathbf{N}\left(\lambda \right)=\mathbf{A}-\lambda \mathbf{B}$ is a linear matrix function of $\lambda$, so that it is a linear eigenvalue problem.

Most numerical methods that deal with the nonlinear eigenvalue problems are Newton-type methods [22,23,24,25]. In [26], some available solution techniques for nonlinear eigenvalue problems using the Jacobi–Davidson, Arnoldi and the rational Krylov methods were presented. Zhou [27] used the Leray–Schauder fixed-point theorem to acquire the existence of positive solutions of a nonlinear eigenvalue problem. El-Ajou [28] demonstrated the general exact and numerical solutions of four significant matrix fractional differential equations, and a new computational skill was applied for obtaining the general solutions of the nonlinear issue in the Caputo sense. Jadamba et al. [29] addressed the nonlinear inverse problem of estimating the stochasticity of a random parameter in stochastic partial differential equations by using a regularized projected stochastic gradient scheme. Later, Harcha et al. [30] tackled the nonlinear eigenvalue problem with the p-Laplacian fractional involving singular weights and obtained the nonexistence of solutions by utilizing a typical version of Picone’s identity.

The nonlinear eigenvalue problem is a great challenge for developing efficient and accurate methods [31]. Even for polynomial nonlinear eigenvalue problems, the linearizations to the linear eigenvalue problems in a larger space are quite complicated and are in general not unique. The present paper intends to overcome these challenges, wherein we will directly solve the nonlinear eigenvalue problem in its nonhomogeneous form by incorporating a normalization condition in the original space.

In this paper, we will encounter a problem to solve a nonlinear equation $f\left(x\right)=0$, but the explicit function of $f\left(x\right)$ is not available. The Newton method for iteratively solving $f\left(x\right)=0$ is given by

$$x_{n+1}=x_n-\frac{f\left(x_n\right)}{f^{\prime }\left(x_n\right)},n=0,1,\dots ,$$

(12)

which needs to carry out a point-wise derivative $f^{\prime }\left(x_n\right)$ in the iteration. For some problems, the explicit function $f\left(x\right)$ might not be available, and this induces great difficulty in using the Newton method to solve the nonlinear equation. To overcome this inefficiency, Liu [32] derived a derivative-free iterative scheme based on a new splitting technique:

$$x_{n+1}=x_n-\frac{f\left(x_n\right)}{a+bf\left(x_n\right)},n=0,1,\dots ,$$

(13)

where a and b are constants. In Section 6.3, we will develop a derivative-free fixed-point Newton method to determine a and b. With regard to the derivative-free fixed-point Newton methods, one can refer to [33] and references therein.

In addition to the derivative-free fixed-point Newton method and the minimization techniques, we will also develop the Newton method for the nonlinear equations system by incorporating the normalization condition into the eigen-equation. Arnoldi [34] proposed that the method of minimized iterations was recommended as a rapid means for determining a small number of the larger eigenvalues and modal columns of a large matrix. After that, many iterative methods were surveyed in [35] at length. Argyros et al. [36] addressed a semilocal analysis of the Newton–Kurchatov method for solving nonlinear equations involving the splitting of an operator. They also acquired weaker sufficient semilocal convergence criteria and tighter error estimates than in earlier works. Argyros and Shakhno [37] employed local convergence of the combined Newton–Kurchatov method for solving Banach-space-valued equations. Further, they also mentioned that these modifications of earlier conditions resulted in tighter convergence analysis and more precise information on the location of the solution.

This paper develops several simple approaches, including two regularization methods, to solve nonlinear eigenvalue problems. The contributions and innovation points of this paper are given as follows:

• When solving nonlinear eigenvalue problems, they can be transformed into minimization problems regardless of real and complex eigenvalues.
• For solving linear and nonlinear eigenvalue problems, this paper presents normalization techniques to create new nonhomogeneous systems and merit functions.
• Two simple regularization methods are combined with the Newton iteration method, which results in very fast convergence to solve nonlinear eigenvalue problems.
• Using the derivative-free fixed-point Newton method to directly solve the regularized scalar equation for nonlinear eigenvalue problems, we can can quickly obtain high-precision eigenvalues.

The remainder of the paper is arranged as follows. In Section 2, we consider a normalization condition for the uniqueness of the eigenvector and derive a simple nonhomogeneous linear system to minimize the residual of the eigen-equation by using the 1D golden section search algorithm (1D GSSA) to determine the real eigenvalue, which results in a simple method (SM). Some examples of linear eigenvalue problems in Section 3 exhibit the advantages of the present methodology of the SM to find the approximate solution of Equation (2). A simple method (SM) of the nonlinear eigenvalue problem (11) is presented in Section 4, which is combined with the golden section search algorithm to be a stable solver of eigenvalues and eigenvectors. For complex eigenvalue problems, we propose two normalization equations with nonhomogeneous terms that appear on the right-hand side. Section 5 displays some examples of nonlinear eigenvalue problems solved by the SM and GSSA. In Section 6, we discuss two simple regularization methods and provide a derivative-free fixed-point Newton method for quickly finding the real eigenvalues. The combination of Newton’s method and regularized equations is carried out in this section. Finally, the conclusions are drawn in Section 7.

2. A Simple Method for Standard Linear Eigenvalue Problems

We can observe that $\mathbf{x}\ne \mathbf{0}$ in Equation (2) is not unique because $\beta \mathbf{x},\beta \ne 0$ is also a solution if $\mathbf{x}$ is a solution of Equation (2). Therefore, for the uniqueness of the eigenvector of Equation (2), an extra normalization condition

$${\mathbf{b}}^{\mathrm{T}}\mathbf{x}=1$$

(14)

can be imposed, where $\mathbf{b}$ is a given nonzero vector. For example, if we take $\mathbf{b}={(1,0,\dots ,0)}^{\mathrm{T}}$, then the first component of $\mathbf{x}$ is normalized to $x_1=1$. If $b_J=1$ for a given J and $b_k=0,k\ne J$, we normalize the J-th component of $\mathbf{x}$ to be $x_J=1$.

Theorem 1.

If $\mathbf{x}\in {\mathbb{R}}^n$ in Equation (2) is constrained by a normalization condition (14) for the uniqueness of $\mathbf{x}$, we can derive a nonhomogeneous system to determine $\mathbf{x}$:

$$({\mathbf{A}}^T\mathbf{A}-\lambda \mathbf{A}-\lambda {\mathbf{A}}^T+{\lambda }^2{\mathbf{I}}_n+\mathbf{b}{\mathbf{b}}^T)\mathbf{x}=\mathbf{b}.$$

(15)

Proof. 

Equations (2) and (14) yield an over-determined system:

$$\left[\begin{array}{c}\mathbf{A}-\lambda {\mathbf{I}}_n\\ {\mathbf{b}}^{\mathrm{T}}\end{array}\right]\mathbf{x}=\left[\begin{array}{c}{\mathbf{0}}_n\\ 1\end{array}\right].$$

(16)

Multiplying both sides by

$${\left[\begin{array}{c}\mathbf{A}-\lambda {\mathbf{I}}_n\\ {\mathbf{b}}^{\mathrm{T}}\end{array}\right]}^{\mathrm{T}}=\left[\begin{array}{cc}{\mathbf{A}}^{\mathrm{T}}-\lambda {\mathbf{I}}_n& \mathbf{b}\end{array}\right],$$

we have

$$\left[\begin{array}{cc}{\mathbf{A}}^{\mathrm{T}}-\lambda {\mathbf{I}}_n& \mathbf{b}\end{array}\right]\left[\begin{array}{c}\mathbf{A}-\lambda {\mathbf{I}}_n\\ {\mathbf{b}}^{\mathrm{T}}\end{array}\right]\mathbf{x}=\left[\begin{array}{cc}{\mathbf{A}}^{\mathrm{T}}-\lambda {\mathbf{I}}_n& \mathbf{b}\end{array}\right]\left[\begin{array}{c}{\mathbf{0}}_n\\ 1\end{array}\right].$$

(17)

Expanding this, we can prove Equation (15). □

Remark 1.

If the coefficient matrix $\mathbf{A}$ is highly ill-conditioned, it is better to replace Equation (15) with

$$(\mathbf{A}-\lambda {\mathbf{I}}_n+\mathbf{b}{\mathbf{b}}^T)\mathbf{x}=\mathbf{b},$$

(18)

the condition number of which is reduced to one-half. Equation (18) is easily derived by adding ${\mathbf{b}}^T\mathbf{x}\mathbf{b}$ on both sides of Equation (2) and taking into account Equation (14) with ${\mathbf{b}}^T\mathbf{x}=1$.

If $\lambda$ is an eigenvalue and $\mathbf{x}\ne \mathbf{0}$ is the corresponding eigenvector, $\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel =0$ by Equation (2). In other cases, $\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel >0\forall \mathbf{x}\ne \mathbf{0}$. As a consequence, $\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel \ge 0\forall \mathbf{x}\in {\mathbb{R}}^n$. For a given $\lambda \in \mathbb{R}$, if $\mathbf{x}$ is solved from Equation (18), then we can determine the correct value of $\lambda$ by minimizing the following merit function:

$$\underset{\lambda \in [a,b]}{min}f\left(\lambda \right):=\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel \ge 0,$$

(19)

where $\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel$ denotes the Euclidean norm of $\mathbf{A}\mathbf{x}-\lambda \mathbf{x}$, and $\lambda \in [a,b]$ is a real eigenvalue to be sought.

The numerical procedures for solving Equation (2) are given by (i) selecting $[a,b]$ and $\mathbf{b}$, (ii) solving Equation (18) for each required ${\lambda }_i\in [a,b]$, and (iii) applying the one-dimensional golden section search algorithm (1D GSSA) to Equation (19) to pick up the eigenvalue.

Remark 2.

This method, involving merely two Equations: (18) and (19), is the simplest method to find the real eigenvalues of Equation (2) and is thus labeled as a simple method (SM). When a SM is used to solve the nonhomogeneous linear system in Equation (18) for each eigen-parameter λ, we can easily compute the real eigenvalue λ and then the corresponding eigenvector $\mathbf{x}$ with the aid of Equation (19).

3. Examples of Linear Eigenvalue Problems

In order to assess the performance of the newly developed SM, we test some linear eigenvalue problems.

Example 1.

We first demonstrate the case

$$\mathbf{A}=\left[\begin{array}{ccc}6& 2& -2\\ 2& 5& 0\\ -2& 0& 7\end{array}\right],$$

(20)

and the exact eigenvalues are $\lambda \left(\mathbf{A}\right)=\{3,6,9\}$.

Although this example is very simple, we adopt it to test the accuracy and efficiency of the proposed SM since the exact eigenvalues are known. By applying the SM, we take $J=1$ and $\epsilon ={10}^{-15}$ as used in the 1D GSSA. We plot $f\left(\lambda \right)$ with respect to the eigen-parameter over an interval, as shown in Figure 1, for which three minimal points are the corresponding eigenvalues $\lambda =3,6,9$, the corresponding eigenvectors of which are given as follows:

$${\mathbf{x}}^{\left(1\right)}=\left[\begin{array}{c}1\\ -1\\ 0.5\end{array}\right],{\mathbf{x}}^{\left(2\right)}=\left[\begin{array}{c}1\\ 2\\ 2\end{array}\right],{\mathbf{x}}^{\left(3\right)}=\left[\begin{array}{c}1\\ 0.5\\ -1\end{array}\right].$$

(21)

Table 1. Example 1 solved by SM, listing EE, the error $\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel$, and NI.

Exact $\mathit{\lambda }$	3	6	9
$[a,b]$	[2.5,4]	[5,7]	[8,11]
EE	0	0	0
$\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel$	$9.62\times {10}^{-17}$	$2.48\times {10}^{-16}$	0
NI	73	74	90

lists some results obtained by SM, where EE means the error of the eigenvalue and NI denotes the number of iterations.

Example 2.

We consider

$$\mathbf{A}=\left[\begin{array}{cccc}-1.6407& 1.0814& 1.2014& 1.1539\\ 1.0814& 4.1573& 7.4035& -1.0463\\ 1.2014& 7.4035& 2.7890& -1.5737\\ 1.1539& -1.0463& -1.5737& 8.6944\end{array}\right],$$

(22)

where $\lambda \left(\mathbf{A}\right)=\{-4,-2,8,12\}$.

In the SM, we take $\mathbf{b}={(1,1,0,0)}^T$ and $\epsilon ={10}^{-15}$ for the 1D GSSA. We plot $f\left(\lambda \right)$ in Figure 2, the four minimal points of which are the corresponding eigenvalues $\lambda =-4,-2,8,12$.

Table 2. Example 2 solved by SM, listing EE, the error $\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel$, and NI.

Exact $\mathit{\lambda }$	−4	−2	8	12
$[a,b]$	[−5,−3]	[−3,−1]	[7,9]	[11,13]
Numerical $\lambda$	−4.0000734	−1.9999357	7.99995888	12.0000501
$\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel$	$6.78\times {10}^{-16}$	$4.55\times {10}^{-16}$	$4.70\times {10}^{-15}$	$2.28\times {10}^{-15}$
NI	74	75	73	72

lists some results obtained by SM. The eigenvectors corresponding to $\lambda =\{-4,-2,8,12\}$ are given as follows:

$$\begin{array}{c}{\mathbf{x}}^{\left(1\right)}=\left[\begin{array}{c}0.134043\\ 0.865957\\ -0.982563\\ -0.062616\end{array}\right],{\mathbf{x}}^{\left(2\right)}=\left[\begin{array}{c}1.219790\\ -0.219790\\ -0.017379\\ -0.155674\end{array}\right],\hfill \\ {\mathbf{x}}^{\left(3\right)}=\left[\begin{array}{c}0.335959\\ 0.664041\\ 0.526805\\ 1.636074\end{array}\right],{\mathbf{x}}^{\left(4\right)}=\left[\begin{array}{c}0.091303\\ 0.908697\\ 0.855567\\ -0.663053\end{array}\right].\hfill \end{array}$$

(23)

Example 3.

A nonsymmetric matrix is given by

$$\mathbf{A}=\left[\begin{array}{cccc}3.5488& 15.593& 8.5775& -4.0123\\ 2.3595& 24.526& 14.596& -5.8157\\ 0.089953& 27.599& 21.483& -5.8415\\ 1.9227& 55.667& 39.717& -10.558\end{array}\right],$$

(24)

where $\lambda \left(\mathbf{A}\right)=\{1,2,6,30\}$.

With $\mathbf{b}={(1,1,1,1)}^{\mathrm{T}}$, we plot $f\left(\lambda \right)$ in Figure 3, the four minimal points of which are the corresponding eigenvalues $\lambda =1,2,6,30$.

Table 3. Example 3 solved by SM, listing EE, the error $\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel$, and NI.

Exact $\mathit{\lambda }$	1	2	6	30
$[a,b]$	[0,2]	[1.5,2.5]	[5.5,7.5]	[29,31]
Numerical $\lambda$	1.0002159	1.9997765	6.0002064	29.999601
$\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel$	$8.01\times {10}^{-15}$	$1.19\times {10}^{-15}$	$1.28\times {10}^{-15}$	$2.55\times {10}^{-15}$
NI	75	73	74	71

lists some results obtained by SM. The eigenvectors corresponding to $\lambda =\{1,2,6,30\}$ are given as follows:

$$\begin{array}{c}{\mathbf{x}}^{\left(1\right)}=\left[\begin{array}{c}2.320260\\ -2.054856\\ 2.309442\\ -1.574846\end{array}\right],{\mathbf{x}}^{\left(2\right)}=\left[\begin{array}{c}0.271484\\ 0.133857\\ -0.009676\\ 0.604335\end{array}\right],\hfill \\ {\mathbf{x}}^{\left(3\right)}=\left[\begin{array}{c}-0.281715\\ -0.120998\\ 0.542044\\ 0.860669\end{array}\right],{\mathbf{x}}^{\left(4\right)}=\left[\begin{array}{c}0.107796\\ 0.174668\\ 0.238818\\ 0.478718\end{array}\right].\hfill \end{array}$$

(25)

Example 4.

Consider the Frank nonsymmetric matrix [38,39,40]:

$$\mathbf{A}=\left[\begin{array}{ccccc}1& 1& 1& \dots & 1\\ 1& 2& 2& \dots & 2\\ 0& 2& 3& \dots & 3\\ ⋮& 0& \ddots & \dots & ⋮\\ 0& 0& \dots & n-1& n\end{array}\right].$$

(26)

For $n=30$, the largest eigenvalue is given by

$$\lambda =96.200622293285.$$

(27)

In the SM, we take $J=10$ and $\epsilon ={10}^{-15}$ for the 1D GSSA. The function $f\left(\lambda \right)$ is plotted in Figure 4, for which the minimal points are the last nine eigenvalues. Through 103 iterations of the 1D GSSA, $\lambda =96.20062229328501$ and $\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel =3.22\times {10}^{-16}$ are obtained.

Example 5.

Let $\mathbf{A}=\left[a_{ij}\right],i,j=1,\dots ,n$. In Equation (2), we take the Hilbert matrix:

$$a_{ij}=\frac{1}{i+j-1},i,j=1,\dots ,n.$$

(28)

Since the Hilbert matrix is highly ill-conditioned, we take Equation (18) instead of Equation (15) to compute the eigenvalue. In Figure 5, with $\mathbf{b}={\mathbf{1}}_n,n=7$, we plot $f\left(\lambda \right)$ with respect to the eigen-parameter over an interval, for which the seven minimal points are the corresponding eigenvalues. There exists the largest eigenvalue between 1.5 and 2.

In

Table 4. For Example 5, the largest eigenvalues of the Hilbert matrices with n = 7, 8, 9, 10.

n	7	8	9	10
Present	1.66088533892693	1.69593899692195	1.72588266090185	1.75191967026518
[41]	1.66088533892693	1.69593899692195	1.72588266090185	1.75191967026518
[3]	1.66088533892658	1.69593899692256	1.72588266090195	1.75191967026518

, we list the largest eigenvalues for $n=7,\dots ,10$ obtained in [41] using the cyclic Jacobi method [42] and in [3] using the external excitation method.

Due to the highly ill-conditioned nature of the Hilbert matrix with $n=100$, it is a quite difficult linear eigenvalue problem. For this problem, we take $J=1$ to compute the largest eigenvalue, which is given by $\lambda =2.182696097757424$. The SM converges very fast with 69 iterations under $\epsilon ={10}^{-15}$, and the error of the eigen-equation is $\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel =2.4\times {10}^{-15}$. Notice that the smallest eigenvalue of the Hilbert matrix with a large n is very difficult to compute since it is very close to zero. However, for $n=100$ and $b_{49}=b_{50}=1$, we can obtain the smallest eigenvalue $6.18\times {10}^{-30}$, whose error is $\parallel \mathbf{A}\mathbf{x}-\lambda \mathbf{x}\parallel =1.41\times {10}^{-16}$.

4. A Simple Method for Nonlinear Eigenvalue Problems

If $\mathbf{x}$ is an eigenvector of Equation (11), then with $\alpha \ne 0$, $\alpha \mathbf{x}$ is also an eigenvector, which means that the eigenvector of Equation (11) is not unique. Therefore, we can impose on Equation (11) an extra normalization condition (14).

Theorem 2.

If $\mathbf{x}\in {\mathbb{R}}^n$ in Equation (11) is imposed by a normalization condition (14) for the uniqueness of $\mathbf{x}\in {\mathbb{R}}^n$, we can derive a nonhomogeneous equation system to determine $\mathbf{x}$:

$$[\mathbf{N}\left(\lambda \right)+\mathbf{b}{\mathbf{b}}^T]\mathbf{x}=\mathbf{b}.$$

(29)

Proof. 

Equation (29) is easily derived by adding ${\mathbf{b}}^{\mathrm{T}}\mathbf{x}\mathbf{b}$ on both sides of Equation (11):

$$\mathbf{N}\left(\lambda \right)\mathbf{x}+\mathbf{b}{\mathbf{b}}^{\mathrm{T}}\mathbf{x}={\mathbf{b}}^{\mathrm{T}}\mathbf{x}\mathbf{b},$$

(30)

which, with ${\mathbf{b}}^{\mathrm{T}}\mathbf{x}=1$ by Equation (14) being taken on the right-hand side, yields

$$\mathbf{N}\left(\lambda \right)\mathbf{x}+\mathbf{b}{\mathbf{b}}^{\mathrm{T}}\mathbf{x}=\mathbf{b}.$$

(31)

Thus, we prove Equation (29). □

The numerical procedures for determining the real eigenvalues of Equation (11) are summarized as follows: (i) Select $[a,b]$ and $\mathbf{b}$. (ii) Solve Equation (29) for each required ${\lambda }_i\in [a,b]$. (iii) Apply the 1D GSSA to Equation (19) with $f\left(\lambda \right)=\parallel \mathbf{N}\left(\lambda \right)\mathbf{x}\parallel$ for picking up the eigenvalue. This method is the simplest method to find the real eigenvalues of Equation (11), which is labeled as a simple method (SM).

Theorem 3.

If $\mathbf{x}\in {\mathbb{C}}^n$ in Equation (11) is imposed by a normalization condition (14) for the uniqueness of $\mathbf{x}$, and λ is a complex eigenvalue, we can derive a nonhomogeneous equations system to determine $\mathbf{x}$:

$$\left[\begin{array}{cc}{\mathbf{N}}_1+\mathbf{b}{\mathbf{b}}^T& -{\mathbf{N}}_2\\ {\mathbf{N}}_2& {\mathbf{N}}_1+\mathbf{b}{\mathbf{b}}^T\end{array}\right]\left[\begin{array}{c}\mathbf{u}\\ \mathbf{v}\end{array}\right]=\left[\begin{array}{c}\mathbf{b}\\ {\mathbf{0}}_n\end{array}\right],$$

(32)

where

$$\lambda ={\lambda }_R+i{\lambda }_I,\mathbf{N}={\mathbf{N}}_1+i{\mathbf{N}}_2,\mathbf{x}=\mathbf{u}+i\mathbf{v}.$$

(33)

Proof. 

Inserting Equation (33) into Equation (29) yields

$$({\mathbf{N}}_1+\mathbf{b}{\mathbf{b}}^T+i{\mathbf{N}}_2)(\mathbf{u}+i\mathbf{v})=\mathbf{b}.$$

(34)

Equating the real and imaginary parts of Equation (34), we have

$$\begin{array}{ccc}\hfill & & ({\mathbf{N}}_1+\mathbf{b}{\mathbf{b}}^T)\mathbf{u}-{\mathbf{N}}_2\mathbf{v}=\mathbf{b},\hfill \\ & & {\mathbf{N}}_2\mathbf{u}+({\mathbf{N}}_1+\mathbf{b}{\mathbf{b}}^T)\mathbf{v}={\mathbf{0}}_n.\hfill \end{array}$$

(35)

which can be recast to that in Equation (32). □

When $\mathbf{x}$ is solved from Equation (32), we can employ the following minimization:

$$\underset{({\lambda }_R,{\lambda }_I)\in [a,b]\times [c,d]}{min}f({\lambda }_R,{\lambda }_I):=\parallel \mathbf{N}({\lambda }_R,{\lambda }_I)\mathbf{x}\parallel \ge 0$$

(36)

to determine the complex eigenvalue.

For the complex eigenvalue problem, we can also derive another normalization equation.

Theorem 4.

For $\mathbf{x}=\mathbf{u}+i\mathbf{v}\in {\mathbb{C}}^n$ in Equation (11), it is imposed by a normalization condition:

$${\mathbf{c}}^T\mathbf{y}=1,\mathbf{y}:=\left[\begin{array}{c}\mathbf{u}\\ \mathbf{v}\end{array}\right]$$

(37)

for the uniqueness of $\mathbf{x}$, and λ is a complex eigenvalue. We can derive a nonhomogeneous equation system to determine $\mathbf{x}=\mathbf{u}+i\mathbf{v}$:

$$(\mathbf{D}+\mathbf{c}{\mathbf{c}}^T)\mathbf{y}=\mathbf{c},$$

(38)

where $\mathbf{c}$ is a $2n$-dimensional constant vector and

$$\mathbf{D}:=\left[\begin{array}{cc}{\mathbf{N}}_1& -{\mathbf{N}}_2\\ {\mathbf{N}}_2& {\mathbf{N}}_1\end{array}\right].$$

(39)

Proof. 

Inserting Equation (33) into Equation (11) and equating the real and imaginary parts, we can derive

$$\mathbf{D}\mathbf{y}=\mathbf{0}.$$

(40)

Adding $\mathbf{c}{\mathbf{c}}^{\mathrm{T}}\mathbf{y}$ on both sides yields

$$(\mathbf{D}+\mathbf{c}{\mathbf{c}}^{\mathrm{T}})\mathbf{y}=\mathbf{c}{\mathbf{c}}^{\mathrm{T}}\mathbf{y},$$

(41)

which by using Equation (37) on the right-hand side generates Equation (38). □

Therefore, the numerical procedures for determining the complex eigenvalue of a nonlinear eigenvalue problem (11) are depicted as follows: (i) Select $[a,b]\times [c,d]$ and $\mathbf{b}$. (ii) For each $({\lambda }_R,{\lambda }_I)\in [a,b]\times [c,d]$ required, we solve Equation (32). (iii) Apply the two-dimensional GSSA (2D GSSA) to Equation (36). With regard to the two-dimensional golden section search algorithm, one may refer to [43].

Remark 3.

Even the proofs of Theorems 2–4 are simple and straightforward; they are crucial for the development of the proposed numerical methods for effectively and accurately solving nonlinear eigenvalue problems.

5. Examples of Nonlinear Eigenvalue Problems

Example 6.

To demonstrate the new idea in Equation (29), we consider a generalized eigenvalue problem $\mathbf{A}\mathbf{x}=\lambda \mathbf{B}\mathbf{x}$ endowed by [44]:

$$\mathbf{A}=\left[\begin{array}{ccccc}2& 3& 4& 5& 6\\ 4& 4& 5& 6& 7\\ 0& 3& 6& 7& 8\\ 0& 0& 2& 8& 9\\ 0& 0& 0& 1& 10\end{array}\right],\mathbf{B}=\left[\begin{array}{ccccc}1& -1& -1& -1& -1\\ 0& 1& -1& -1& -1\\ 0& 0& 1& -1& -1\\ 0& 0& 0& 1& -1\\ 0& 0& 0& 0& 1\end{array}\right].$$

(42)

By using the SM with $\mathbf{b}={(1,0,0,0,0)}^{\mathrm{T}}$ and $\epsilon ={10}^{-15}$ for the 1D GSSA, $f\left(\lambda \right)$ is plotted in Figure 6, for which five eigenvalues appear as minimums. With $[a,b]=[-1,0]$, we obtain $\lambda =-0.1873528931969768$, NI = 73 and $\parallel \mathbf{N}\mathbf{x}\parallel =1.07\times {10}^{-15}$. Those eigenvalues are listed as follows:

$$\begin{array}{ccc}\hfill & & {\lambda }_1=-0.1873528931969768,{\lambda }_2=1.313278952662423,\hfill \\ & & {\lambda }_3=5.537956370847891,{\lambda }_4=12.0896928530668,{\lambda }_5=21.24642471661987.\hfill \end{array}$$

(43)

Example 7.

To display the advantage of Equation (32), we consider a standard eigenvalue problem with

$$\mathbf{A}=\left[\begin{array}{cc}1& -1\\ 1& 2\end{array}\right],\mathbf{B}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],$$

(44)

which possesses the complex eigenvalues:

$$\lambda =\frac{3}{2}\pm \frac{i\sqrt{3}}{2}.$$

(45)

By using the SM with $\mathbf{b}={(1,1)}^{\mathrm{T}}$ and $\epsilon ={10}^{-15}$ for the 2D GSSA and with $[a,b]\times [c,d]=[1.2,1.9]\times [0.5,0.9]$, we can obtain NI = 73 and $\parallel \mathbf{N}\mathbf{x}\parallel =7.22\times {10}^{-16}$, and the error of the eigenvalue is $4.44\times {10}^{-16}$. The SM is very accurate in finding the complex eigenvalue.

On the other hand, using the SM in Theorem 4 with $\mathbf{c}={(0,1,0,0)}^T$ and $\epsilon ={10}^{-15}$ for the 2D GSSA and with $[a,b]\times [c,d]=[1.2,1.9]\times [0.5,0.9]$, we can obtain NI = 73 and $\parallel \mathbf{N}\mathbf{x}\parallel =1.67\times {10}^{-16}$, and the error of the obtained eigenvalue is zero.

Example 8.

We consider

$$\mathbf{N}\left(\lambda \right)\mathbf{x}=({\lambda }^2\mathbf{M}+\sqrt{\lambda }\mathbf{C}+\mathbf{K})\mathbf{x}=\mathbf{0},$$

(46)

where

$$\mathbf{M}=\left[\begin{array}{ccc}0& 6& 0\\ 0& 6& 0\\ 0& 0& 1\end{array}\right],\mathbf{C}=\left[\begin{array}{ccc}1& -6& 0\\ 2& -7& 0\\ 0& 0& 0\end{array}\right],\mathbf{K}={\mathbf{I}}_3.$$

(47)

By using the SM with $\mathbf{b}={(1,1,1)}^{\mathrm{T}}$ and $\epsilon ={10}^{-15}$, we plot $f\left(\lambda \right)$ in Figure 7a with respect to the eigen-parameter over an interval, for which two real eigenvalues are $\lambda =0.04080314176866112$ and $\lambda =0.7425972620277184$, where $\parallel \mathbf{N}\mathbf{x}\parallel =1.22\times {10}^{-15}$ and $\parallel \mathbf{N}\mathbf{x}\parallel =2.53\times {10}^{-15}$ are obtained, respectively.

Let $\sqrt{\lambda }=\mu$, and we can derive

$${\mathbf{N}}_1=[{({\mu }_R^2-{\mu }_I^2)}^2-4{\mu }_R^2{\mu }_I^2]\mathbf{M}+{\mu }_R\mathbf{C}+\mathbf{K},{\mathbf{N}}_2=4{\mu }_R{\mu }_I({\mu }_R^2-{\mu }_I^2)\mathbf{M}+{\mu }_I\mathbf{C}.$$

(48)

There are a total of 24 eigenvalues, as shown in Figure 7b, with respect to ${\lambda }_R={\mu }_R^2-{\mu }_I^2$ and ${\lambda }_I=2{\mu }_R{\mu }_I$.

Example 9.

As an application, we consider a time-delay linear system of first-order ordinary differential equations:

$$\dot{\mathbf{q}}\left(t\right)=\mathbf{A}\mathbf{q}\left(t\right)+\mathbf{B}\mathbf{q}(t-1),$$

(49)

where $\mathbf{B}\mathbf{q}(t-1)$ is a time-delay external force. Inserting $\mathbf{q}\left(t\right)=e^{\lambda t}\mathbf{x}$ into Equation (49) renders

$$\lambda e^{\lambda t}\mathbf{x}=\mathbf{A}{e}^{\lambda t}\mathbf{x}+e^{-\lambda }{e}^{\lambda t}\mathbf{B}\mathbf{x}.$$

(50)

By canceling $e^{\lambda t}$ on both sides, we obtain a time-delay nonlinear eigenvalue problem:

$$\mathbf{A}\mathbf{x}+e^{-\lambda }\mathbf{B}\mathbf{x}=\lambda \mathbf{x},$$

(51)

where $\mathbf{A},\mathbf{B}\in {\mathbb{R}}^{n\times n}$. The eigenvalues of this system are very important as they reflect the stability of the time-delay system.

We take $n=3$ and consider [25]

$$\mathbf{N}\left(\lambda \right)\mathbf{x}=[\mathbf{A}+exp(-\lambda )\mathbf{B}-\lambda {\mathbf{I}}_3]\mathbf{x}=\mathbf{0},$$

(52)

where

$$\mathbf{A}=\left[\begin{array}{ccc}0& 1& 0\\ 0& 0& 1\\ -a_3& -a_2& -a_1\end{array}\right],\mathbf{B}=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ -b_3& -b_2& -b_1\end{array}\right].$$

(53)

This describes a time-delay system.

We take $a_1=1.5,a_2=1,a_3=0.5$ and $b_1=0.3,b_2=0.2,b_3=0.1$, and there exist four complex eigenvalues. Through some manipulations, we find that

$${\mathbf{N}}_1({\lambda }_R,{\lambda }_I)=\mathbf{A}+exp(-{\lambda }_R)cos{\lambda }_I\mathbf{B}-{\lambda }_R{\mathbf{I}}_3,{\mathbf{N}}_2({\lambda }_R,{\lambda }_I)=-exp(-{\lambda }_R)sin{\lambda }_I\mathbf{B}-{\lambda }_I{\mathbf{I}}_3.$$

(54)

We apply the SM to solve this problem with $\mathbf{b}={(1,1,1)}^{\mathrm{T}}$ and $\epsilon ={10}^{-15}$. In $[a,b]\times [c,d]=[-1.5,1.5]\times [0,2]$, we obtain NI = 76 and $\parallel \mathbf{N}\mathbf{x}\parallel =9.76\times {10}^{-16}$. The complex eigenvalue is $\lambda =-3.208498325636312\pm i0.6608850649667832$. In $[a,b]\times [c,d]=[-1.5,-1.4]\times [0.9,1.1]$, we obtain NI = 75 and $\parallel \mathbf{N}\mathbf{x}\parallel =1.05\times {10}^{-15}$. The complex eigenvalue is $\lambda =-1.422926059141501\pm i1.035178169828797$.

By taking the parameters of $a_i,b_i,i=1,2,3$ as those listed in [25,45], there exists a double non-semisimple eigenvalue $3\pi i$. With $[a,b]\times [c,d]=[-{10}^{-15},{10}^{-15}]\times [9.42,9.45]$, we obtain NI = 62 and $\parallel \mathbf{N}\mathbf{x}\parallel =2\times {10}^{-15}$. The eigenvalue obtained is very close to $3\pi i$, with an error of $4.72\times {10}^{-8}$.

Example 10.

Consider

$$\mathbf{N}\left(\lambda \right)={\lambda }^2{\mathbf{A}}_2+\lambda {\mathbf{A}}_1+{\mathbf{A}}_0,$$

(55)

where

$${\mathbf{A}}_2=\left[\begin{array}{ccc}-4& 3& 12\\ -17& -11& 0\\ 1& -1& 3\end{array}\right],{\mathbf{A}}_1=\left[\begin{array}{ccc}2& -6& 1\\ -2& 22& 11\\ 7& -1& 1\end{array}\right],{\mathbf{A}}_0=\left[\begin{array}{ccc}-16& -4& 7\\ -14& 7& 13\\ 6& 8& 7\end{array}\right].$$

(56)

We can derive

$${\mathbf{N}}_1({\lambda }_R,{\lambda }_I)=({\lambda }_R^2-{\lambda }_I^2){\mathbf{A}}_2+{\lambda }_R{\mathbf{A}}_1+{\mathbf{A}}_0,{\mathbf{N}}_2({\lambda }_R,{\lambda }_I)=2{\lambda }_R{\lambda }_I{\mathbf{A}}_2+{\lambda }_I{\mathbf{A}}_1.$$

(57)

As shown in Figure 8, there exist two minima of $f\left(\lambda \right)$ for this quadratic eigenvalue problem having two real eigenvalues. We apply the SM to solve this problem with $\mathbf{b}={(0,1,0)}^T$ and $\epsilon ={10}^{-15}$. The first real eigenvalue is $\lambda =-0.2328574586400297$, and we obtain NI = 71 and $\parallel \mathbf{N}\mathbf{x}\parallel =3.97\times {10}^{-15}$. The second real eigenvalue is $\lambda =2.355885632295363$, and we obtain NI = 72 and $\parallel \mathbf{N}\mathbf{x}\parallel =8.33\times {10}^{-15}$. Then, by applying the 2D GSSA to solve this problem with $[a,b]\times [c,d]=[-1,1]\times [0.1,2]$, we obtain NI = 75 and $\parallel \mathbf{N}\mathbf{x}\parallel =2.95\times {10}^{-14}$, and the complex eigenvalues are $\lambda =-0.188835959350602\pm i1.06014959301131$.

Similarly, we adopt the second normalization equation in Theorem 4 with $\mathbf{c}={(1,0,0,0,0,0)}^T$, and by applying the 2D GSSA to solve this problem with $[a,b]\times [c,d]=[-1,1]\times [0.1,2]$, we obtain NI = 75 and $\parallel \mathbf{N}\mathbf{x}\parallel =5.73\times {10}^{-14}$, and the complex eigenvalues are $\lambda =-0.188835959350602\pm i1.06014959301131$, which are equal to those obtained above with the SM derived from Theorem 3.

Example 11

(From [25]).

$$\mathbf{N}\left(\lambda \right)={\lambda }^3{\mathbf{A}}_3+{\lambda }^2{\mathbf{A}}_2+{\mathbf{A}}_0,$$

(58)

where

$${\mathbf{A}}_3=\left[\begin{array}{ccc}-4& 3& 12\\ -17& -11& 0\\ 1& -1& 3\end{array}\right],{\mathbf{A}}_2=\left[\begin{array}{ccc}2& -6& 1\\ -2& 22& 11\\ 7& -1& 1\end{array}\right],{\mathbf{A}}_0=\left[\begin{array}{ccc}-16& -4& 7\\ -14& 7& 13\\ 6& 8& 7\end{array}\right],$$

(59)

we can derive

$${\mathbf{N}}_1({\lambda }_R,{\lambda }_I)=({\lambda }_R^3-3{\lambda }_R{\lambda }_I^2){\mathbf{A}}_3+({\lambda }_R^2-{\lambda }_I^2){\mathbf{A}}_2+{\mathbf{A}}_0,{\mathbf{N}}_2({\lambda }_R,{\lambda }_I)=(3{\lambda }_R^2{\lambda }_I-{\lambda }_I^3){\mathbf{A}}_3+2{\lambda }_R{\lambda }_I{\mathbf{A}}_2.$$

(60)

We adopt the second normalization Equation (38) in Theorem 4 to solve this nonlinear eigenvalue problem. By applying the 2D GSSA to solve this problem with $[a,b]\times [c,d]=[0.01,0.03]\times [0.4,0.5]$ and with $\mathbf{c}={(0,0,0,0,0,1)}^T$, NI = 68, $\parallel \mathbf{N}\mathbf{x}\parallel =3.73\times {10}^{-15}$, and $\lambda =0.02570242595103074+i0.4701394321627313$ can be obtained.

6. A Derivative-Free Newton Method and Regularizations for the Simple Method

6.1. Two Regularization Methods

In the first regularization method (FRM), we take

$$\mathbf{b}=\alpha \mathbf{d},$$

(61)

where $\alpha \ne 0$ is a regularization parameter and $\mathbf{d}$ is a constant vector. Inserting Equation (61) into Equations (29) and (14) generates the first regularization equation:

$$[\mathbf{N}\left(\lambda \right)+{\alpha }^2\mathbf{d}{\mathbf{d}}^{\mathrm{T}}]\mathbf{x}=\alpha \mathbf{d},\alpha {\mathbf{d}}^{\mathrm{T}}\mathbf{x}=1.$$

(62)

In the second regularization method (SRM), we consider another normalization condition instead of Equation (14):

$${\mathbf{d}}^{\mathrm{T}}\mathbf{x}=\alpha .$$

(63)

If $\alpha =1$, we recover to Equation (14). Then, as done in the proof of Theorem 2, we can derive the second regularization equation:

$$[\mathbf{N}\left(\lambda \right)+\mathbf{d}{\mathbf{d}}^{\mathrm{T}}]\mathbf{x}=\alpha \mathbf{d}.$$

(64)

Equations (62) and (64) are different as the ${\alpha }^2\mathbf{d}{\mathbf{d}}^{\mathrm{T}}$ in the first one becomes $\mathbf{d}{\mathbf{d}}^{\mathrm{T}}$ in the second one. The second regularization method is simpler than the first regularization method.

6.2. Newton Iterative Methods

The regularized Equation (62) consists of a system of nonlinear equations for $\mathbf{x}$ and $x_m:=\lambda ,m=n+1$, the Jacobian matrix of which at the k-th step is

$${\mathbf{J}}_k=\left[\begin{array}{cc}\mathbf{N}\left(x_m^k\right)+{\alpha }^2\mathbf{d}{\mathbf{d}}^T& {\mathbf{N}}^{\prime }\left(x_m^k\right){\mathbf{x}}^k\\ \alpha {\mathbf{d}}^T& 0\end{array}\right].$$

(65)

Thus, the Newton iterative method together with the FRM is given by

$$\left[\begin{array}{c}{\mathbf{x}}^{k+1}\\ x_m^{k+1}\end{array}\right]=\left[\begin{array}{c}{\mathbf{x}}^k\\ x_m^k\end{array}\right]-{\left[\begin{array}{cc}\mathbf{N}\left(x_m^k\right)+{\alpha }^2\mathbf{d}{\mathbf{d}}^T& {\mathbf{N}}^{\prime }\left(x_m^k\right){\mathbf{x}}^k\\ \alpha {\mathbf{d}}^T& 0\end{array}\right]}^{-1}\left[\begin{array}{c}\mathbf{N}\left(x_m^k\right){\mathbf{x}}^k+{\alpha }^2\mathbf{d}·{\mathbf{x}}^k\mathbf{d}-\alpha \mathbf{d}\\ \alpha \mathbf{d}·{\mathbf{x}}^k-1\end{array}\right].$$

(66)

To construct the Newton method, letting $x_{n+1}=\lambda$ and using Equations (64) and (63) yields the following nonlinear regularized equations with dimension $m=n+1$:

$$\begin{array}{ccc}\hfill & & [\mathbf{N}\left(x_m\right)+\mathbf{d}{\mathbf{d}}^T]\mathbf{x}=\alpha \mathbf{d},\hfill \end{array}$$

(67)

$$\begin{array}{ccc}\hfill & & d_1{x}_1+\dots +d_n{x}_n-\alpha =0.\hfill \end{array}$$

(68)

At the k-th step, the Jacobian matrix reads as

$${\mathbf{J}}_k=\left[\begin{array}{cc}\mathbf{N}\left(x_m^k\right)+\mathbf{d}{\mathbf{d}}^T& {\mathbf{N}}^{\prime }\left(x_m^k\right){\mathbf{x}}^k\\ {\mathbf{d}}^T& 0\end{array}\right].$$

(69)

Then, the Newton iterative method together with the SRM is given by

$$\left[\begin{array}{c}{\mathbf{x}}^{k+1}\\ x_m^{k+1}\end{array}\right]=\left[\begin{array}{c}{\mathbf{x}}^k\\ x_m^k\end{array}\right]-{\left[\begin{array}{cc}\mathbf{N}\left(x_m^k\right)+\mathbf{d}{\mathbf{d}}^T& {\mathbf{N}}^{\prime }\left(x_m^k\right){\mathbf{x}}^k\\ {\mathbf{d}}^T& 0\end{array}\right]}^{-1}\left[\begin{array}{c}\mathbf{N}\left(x_m^k\right){\mathbf{x}}^k+\mathbf{d}·{\mathbf{x}}^k\mathbf{d}-\alpha \mathbf{d}\\ \mathbf{d}·{\mathbf{x}}^k-\alpha \end{array}\right].$$

(70)

The iteration is terminated if a given convergence criterion $\epsilon ={10}^{-15}$ is fulfilled.

As usual, the Newton iterative method needs to invert the Jacobian matrix at each iteration step, which may result in much computational time for a large-scale eigenvalue problem. In order to save the computational time, we return to solve the scalar equation $\parallel \mathbf{N}\left(\lambda \right)\mathbf{x}\parallel =0$ given below.

6.3. A Derivative-Free Fixed-Point Newton Method

The NI of 1D GSSA is usually over 70, as shown by Examples 6–10 in Section 5. To reduce the computational burden, a derivative-free fixed-point Newton method (DFFPNM) for solving a scalar equation $F\left(\lambda \right)=0$ can be derived as follows [46,47]. Mathematically speaking, solving Equation (11) is equivalent to solving

$$F\left(\lambda \right):=\parallel \mathbf{N}\left(\lambda \right)\mathbf{x}\parallel =0,$$

(71)

which, however, after inserting the solution $\mathbf{x}$ obtained from Equation (62) or Equation (64), is a highly nonlinear and implicit function of $\lambda$.

Suppose that ${\lambda }^{*}$ is root with $F\left({\lambda }^{*}\right)=0$. In order to get rid of the derivative term in the Newton method, we consider

$$F^{\prime }\left({\lambda }^n\right)=F^{\prime }\left({\lambda }^{*}\right)+F^{″}\left({\lambda }^{*}\right)({\lambda }^n-{\lambda }^{*})+\frac{1}{2}{F}^{‴}\left({\lambda }^{*}\right){({\lambda }^n-{\lambda }^{*})}^2+\cdots .$$

(72)

Neglecting the higher-order terms and inserting it into the Newton iterative scheme, we have

$${\lambda }^{n+1}={\lambda }^n-\frac{F\left({\lambda }^n\right)}{F^{\prime }\left({\lambda }^{*}\right)+F^{″}\left({\lambda }^{*}\right)({\lambda }^n-{\lambda }^{*})},$$

(73)

which, combining with

$$F\left({\lambda }^n\right)=F\left({\lambda }^{*}\right)+F^{\prime }\left({\lambda }^{*}\right)({\lambda }^n-{\lambda }^{*})=F^{\prime }\left({\lambda }^{*}\right)({\lambda }^n-{\lambda }^{*}),$$

(74)

yields

$${\lambda }^{n+1}={\lambda }^n-\frac{F\left({\lambda }^n\right)}{a+bF\left({\lambda }^n\right)},$$

(75)

where

$$a=F^{\prime }\left({\lambda }^{*}\right),b=\frac{F^{″}\left({\lambda }^{*}\right)}{F^{\prime }\left({\lambda }^{*}\right)}.$$

(76)

Determining a and b by a fixed-point estimation, the first step is to choose two initial guesses ${\lambda }^0$ and ${\lambda }^2$. Then, we take ${\lambda }^1=({\lambda }^0+{\lambda }^2)/2$. As the approximations of a and b in Equation (76), we can evaluate them by using the technique of finite difference:

$$a=\frac{F\left({\lambda }^2\right)-F\left({\lambda }^0\right)}{{\lambda }^2-{\lambda }^0},b=\frac{1}{a}\frac{F\left({\lambda }^2\right)-2F\left({\lambda }^1\right)+F\left({\lambda }^0\right)}{{({\lambda }^1-{\lambda }^0)}^2}.$$

(77)

The resulting iterative scheme (75) together with a and b above is a derivative-free fixed-point Newton method (DFFPNM).

6.4. Numerical Tests

We employ Example 6 in Section 5 to demonstrate the effectiveness of these two regularization methods. In the first regularization method (FRM), we take $\alpha =4$ and $\mathbf{d}={(1,0,0,0,0)}^{\mathrm{T}}$, and in the second regularization method (SRM), we take $\alpha =0.1$ and $\mathbf{d}={(1,0,0,0,0)}^{\mathrm{T}}$.

Comparing the two curves in Figure 9 to Figure 6, we see it is easier to find the zero points of $F\left(\lambda \right)=0$ by the Newton-like method.

Now, we solve Example 6 in Section 5 by the first regularization method (FRM).

Table 5. Results of Example 6 solved by FRM and DFFPNM.

$\mathit{\lambda }$	$\mathit{\alpha }$	$\parallel \mathbf{N}\mathbf{x}\parallel$	NI	$[{\mathit{\lambda }}^0, {\mathit{\lambda }}^2]$
−0.1873528931969755	4	$5.93\times {10}^{-16}$	9	[−0.18,−0.17]
1.313278952662422	4	$5.57\times {10}^{-16}$	12	[1.2,1.3]
5.537956370847892	4	$5\times {10}^{-16}$	7	[5.51,5.52]
12.0896928530668	4	$9.81\times {10}^{-16}$	8	[11.8,11.9]
21.24642471661986	5	$1.11\times {10}^{-16}$	10	[21.1,21.2]

lists the eigenvalue, $\alpha$, error $\parallel \mathbf{N}\mathbf{x}\parallel$, number of iterations (NI), and $[{\lambda }^0,{\lambda }^2]$ results used in the DFFPNM.

Next we solve Example 6 by the second regularization method (SRM).

Table 6. Results of Example 6 solved by SRM and DFFPNM.

$\mathit{\lambda }$	$\mathit{\alpha }$	$\parallel \mathbf{N}\mathbf{x}\parallel$	NI	$[{\mathit{\lambda }}^0, {\mathit{\lambda }}^2]$
−0.1873528931969756	0.1	$2.77\times {10}^{-16}$	14	[−0.18,−0.17]
1.313278952662424	0.1	$6.07\times {10}^{-16}$	13	[1.32,1.33]
5.537956370847892	0.1	$1.85\times {10}^{-16}$	12	[5.54,5.55]
12.0896928530668	0.1	$4.04\times {10}^{-16}$	11	[12.09,12.1]
21.24642471661986	0.1	$8.89\times {10}^{-17}$	7	[21.247,21.248]

lists the eigenvalue, $\alpha$, $\parallel \mathbf{N}\mathbf{x}\parallel$, NI, and $[{\lambda }^0,{\lambda }^2]$ results used in the DFFPNM.

In

Table 7. Results of Example 6 solved by the Newton method and FRM.

$\mathit{\lambda }$	${\mathit{c}}_0$	$\parallel \mathbf{N}\mathbf{x}\parallel$	NI
−0.1873528931969773	−0.2	$1.0\times {10}^{-15}$	5
1.313278952662422	1.5	$8.05\times {10}^{-16}$	6
5.537956370847892	5	$4.34\times {10}^{-16}$	7
12.0896928530668	12	$4.58\times {10}^{-16}$	6
21.24642471661986	22	$1.43\times {10}^{-17}$	7

, we list the results computed from the Newton method together with the FRM, for which the initial guess is ${\mathbf{x}}^0=\mathbf{1}$, ${\lambda }^0=c_0$, and $\alpha =3$.

In

Table 8. Results of Example 6 solved by the Newton method and SRM.

$\mathit{\lambda }$	${\mathit{c}}_0$	$\parallel \mathbf{N}\mathbf{x}\parallel$	NI
−0.1873528931969765	−0.2	$9.6\times {10}^{-18}$	5
1.313278952662422	1.5	$2.1\times {10}^{-17}$	6
5.537956370847892	5	$7.81\times {10}^{-18}$	7
12.0896928530668	12	$1.59\times {10}^{-17}$	5
21.24642471661986	22	$2.78\times {10}^{-17}$	8

, we list the results computed from the Newton method together with the SRM, for which the initial guess is ${\mathbf{x}}^0=\mathbf{1}$, ${\lambda }^0=c_0$, and $\alpha =0.01$.

We solve Example 8 by the FRM with a fixed $\mathbf{d}={(1,1,1)}^{\mathrm{T}}$ and two values of $\alpha =1$ and $\alpha =2$. The two curves are compared in Figure 10: it can be seen that the curve with $\alpha =2$ is better than that with $\alpha =1$.

Table 9. Results of Example 8 solved by FRM and DFFPNM.

$\mathit{\lambda }$	$\mathit{\alpha }$	$\parallel \mathbf{N}\mathbf{x}\parallel$	NI	$[{\mathit{\lambda }}^0,{\mathit{\lambda }}^2]$
0.04080314176866105	1	$4.02\times {10}^{-16}$	7	[0.041,0.042]
0.7425972620277158	1	$8.56\times {10}^{-16}$	19	[0.72,0.721]
0.04080314176866114	2	$4.42\times {10}^{-16}$	5	[0.041,0.042]
0.7425972620277167	2	$7.89\times {10}^{-16}$	13	[0.72,0.721]

reveals that $\alpha =2$ converges faster than $\alpha =1$.

Here we apply the FRM to solve Example 10 with $\mathbf{d}={(0,1,0)}^{\mathrm{T}}$ and $\epsilon ={10}^{-15}$. For $\alpha =5$, the first real eigenvalue is found to be $\lambda =-0.2328574586400297$, and NI = 7 and $\parallel \mathbf{N}\mathbf{x}\parallel =3.14\times {10}^{-16}$ are obtained. Then, with $\mathbf{d}={(1,1,1)}^{\mathrm{T}}$ and $\alpha =10$, the second real eigenvalue is given by $\lambda =2.355885632295364$, and we obtain NI = 6 and $\parallel \mathbf{N}\mathbf{x}\parallel =9.89\times {10}^{-16}$. The two curves of $F\left(\lambda \right)$ are compared in Figure 11. By comparing the SM together with the 1D GSSA, NI can save ten times.

When we apply the SRM to solve Example 10 with $\alpha ={10}^{-10}$, $\mathbf{d}={(0,1,0)}^T$, and $\epsilon ={10}^{-15}$, we obtain $\lambda =-0.2328576719480232$, NI = 10, and $\parallel \mathbf{N}\mathbf{x}\parallel =4.65\times {10}^{-16}$. For the second real eigenvalue, we obtain $\lambda =2.355885667978917$, NI = 7, and $\parallel \mathbf{N}\mathbf{x}\parallel =1.95\times {10}^{-16}$.

Example 12.

This example is Example 6.2 of [48], and we have a quadratic eigenvalue problem (4) with

$$\begin{array}{ccc}\hfill & & \mathbf{M}=c_{11}{\mathbf{I}}_m\otimes \tilde{\mathbf{M}}+c_{12}\tilde{\mathbf{M}}\otimes {\mathbf{I}}_m,\hfill \\ & & \mathbf{C}=c_{21}{\mathbf{I}}_m\otimes \tilde{\mathbf{C}}+c_{22}\tilde{\mathbf{C}}\otimes {\mathbf{I}}_m,\hfill \\ & & \mathbf{K}=c_{31}{\mathbf{I}}_m\otimes \tilde{\mathbf{K}}+c_{32}\tilde{\mathbf{K}}\otimes {\mathbf{I}}_m,\hfill \end{array}$$

(78)

where

$$\begin{array}{ccc}\hfill & & \tilde{\mathbf{M}}=\frac{1}{6}(4{\mathbf{I}}_m+\mathbf{B}+{\mathbf{B}}^T),\hfill \\ & & \tilde{\mathbf{C}}=\mathbf{B}-{\mathbf{B}}^T,\hfill \\ & & \tilde{\mathbf{K}}=\mathbf{B}+{\mathbf{B}}^T-2{\mathbf{I}}_m,\hfill \\ & & \mathbf{B}=\left[\begin{array}{cccc}0& 0& \dots & 0\\ 1& 0& \dots & 0\\ & \ddots & \ddots & \\ \dots & 0& 1& 0\end{array}\right].\hfill \end{array}$$

(79)

By taking $m=5$, we have $n=25$, and we take $c_{11}=1$, $c_{12}=1.3$, $c_{21}=0.1$, $c_{22}=1.1$, $c_{31}=1$, and $c_{32}=1.2$. In

Table 10. Results of Example 12 solved by the combination of FRM and the Newton method.

$\mathit{\lambda }$	${\mathit{c}}_0$	$\parallel \mathbf{N}\mathbf{x}\parallel$	NI
0.6726432606416950	0.6	$8.45\times {10}^{-16}$	7
0.9866442954682593	1	$6.63\times {10}^{-16}$	7
1.068910201481607	1.06	$7.15\times {10}^{-16}$	7

, we list the results computed from the Newton method and FRM, where $\mathbf{d}={(1,0,\dots ,0)}^T$ and the initial guess is ${\mathbf{x}}^0=\mathbf{1}$, ${\lambda }^0=c_0$, and $\alpha =2.5$.

In

Table 11. Results of Example 12 solved by the combination of SRM and the Newton method.

$\mathit{\lambda }$	${\mathit{c}}_0$	$\parallel \mathbf{N}\mathbf{x}\parallel$	NI
0.672643260641695	0.6	$2.12\times {10}^{-17}$	7
0.9866442954682593	1	$2.2\times {10}^{-17}$	6
1.068910201481607	1.06	$1.02\times {10}^{-17}$	7

, we list the results computed from the Newton method and SRM, where $\mathbf{d}={(1,0,\dots ,0)}^T$ and the initial guess is ${\mathbf{x}}^0=\mathbf{1}$, ${\lambda }^0=c_0$, and $\alpha =0.01$.

Example 13.

As a practical application, we consider a five-story shear building with [49]

$$\begin{array}{ccc}\hfill & & \mathbf{M}=\left[\begin{array}{ccccc}140& 0& 0& 0& 0\\ 0& 120& 0& 0& 0\\ 0& 0& 120& 0& 0\\ 0& 0& 0& 120& 0\\ 0& 0& 0& 0& 100\end{array}\right]kip/g,\hfill \\ & & \mathbf{K}=\left[\begin{array}{ccccc}800& -400& 0& 0& 0\\ -400& 600& -200& 0& 0\\ 0& -200& 400& -200& 0\\ 0& 0& -200& 300& -100\\ 0& 0& 0& -100& 100\end{array}\right]kip/in,\hfill \end{array}$$

(80)

and $\mathbf{C}=\mathbf{0}$. By inserting $\mathbf{q}=e^{i\omega t}\mathbf{x}$ into Equation (3), we can obtain a nonlinear eigenvalue problem:

$$({\omega }^2\mathbf{M}-\mathbf{K})\mathbf{x}=\mathbf{0}.$$

(81)

For the design of engineering structures, knowing the frequencies ω of the free vibration modes $\mathbf{x}$ is of utmost importance.

In

Table 12. Results of Example 13 solved by the combination of FRM and the Newton method.

$\mathit{\omega }$	${\mathit{c}}_0$	$\parallel \mathbf{N}\mathbf{x}\parallel$	NI
0.451662663134403	0.4	$6.32\times {10}^{-15}$	6
1.093583306689083	1.1	$3.63\times {10}^{-14}$	6
1.597324416019366	1.7	$8.64\times {10}^{-14}$	8
2.206998531216506	2.3	$4.26\times {10}^{-13}$	8
2.953880097579612	3	$6.44\times {10}^{-13}$	11

, we list the results computed from the Newton method and FRM, where $\mathbf{d}=\mathbf{1}$ and the initial guess is ${\mathbf{x}}^0=\mathbf{1}$, ${\omega }^0=c_0$, and $\alpha =2$. It can be seen that all $\parallel \mathbf{N}\mathbf{x}\parallel$, where $\mathbf{N}={\omega }^2\mathbf{M}-\mathbf{K}$ are very small.

In

Table 13. Results of Example 13 solved by the combination of SRM and the Newton method.

$\mathit{\omega }$	${\mathit{c}}_0$	$\parallel \mathbf{N}\mathbf{x}\parallel$	NI
0.4516626631344028	0.4	$5.58\times {10}^{-15}$	6
1.093583306689083	1	$1.96\times {10}^{-15}$	7
1.597324416019366	1.5	$6.28\times {10}^{-16}$	7
2.206998531216506	2.1	$4.36\times {10}^{-15}$	9
2.953880097579611	2.9	$1.83\times {10}^{-15}$	7

, we list the results computed from the Newton method and SRM, where $d_1=1,d_j=0,j=2,\dots ,n$ and the initial guess is ${\mathbf{x}}^0=\mathbf{1}$, ${\omega }^0=c_0$, and $\alpha =0.01$. The corresponding five modes of free vibration are plotted in Figure 12, wherein all the first components are normalized to one. It can be seen that all $\parallel \mathbf{N}\mathbf{x}\parallel$ are very small, which indicates the high accuracy of the proposed Newton method based on the SRM, which is slightly more accurate than the FRM in Table 12.

Indeed, the regularization parameter $\alpha$ controls the convergence speed and accuracy. In

Table 14. Example 13 solved by the combination of SRM and the Newton method, listing NI and $\parallel \mathbf{N}\mathbf{x}\parallel$ for different values of $\alpha$.

α	0.01	0.5	1	4	4.8	6
NI	7	7	8	10	12	$>1000$
$\parallel \mathbf{N}\mathbf{x}\parallel$	$1.96\times {10}^{-15}$	$7.11\times {10}^{-14}$	$1.17\times {10}^{-13}$	$2.05\times {10}^{-13}$	$2.54\times {10}^{-13}$	$5.21\times {10}^{-13}$

, we list NI and $\parallel \mathbf{N}\mathbf{x}\parallel$ for different values of $\alpha$. When $\alpha \ge 6$, it does not converge within 1000 steps, and the accuracy is reduced to $5.21\times {10}^{-13}$. The best value is $\alpha =0.01$.

When $\alpha =1$, the normalization condition (63) in the SRM recovers to the normalization condition (14) in the SM. However, $\alpha =1$ is not the best one, as shown in Table 14. When the proper values of $\alpha$ are taken in the FRM and SRM, they are better than the SM.

7. Conclusions

Fast and accurate iterative solution methods of linear eigenvalue problems and nonlinear eigenvalue problems were studied in the paper. We transformed the original homogeneous eigen-equation to a nonhomogeneous linear system by imposing an extra normalization condition for the uniqueness of the eigenvector. Over the range given, the curve of the merit function is constructed, for which the real eigenvalues are local minimums. In the merit function, the vector variable is solved from the newly derived nonhomogeneous linear system. Real eigenvalues by using the 1D golden section search algorithm and complex eigenvalues by using the 2D golden section search algorithm can be obtained quite fast. Very accurate eigenvalues and eigenvectors—as reflected by the very small errors to satisfy the eigen-equation with orders of ${10}^{-15}$ and ${10}^{-16}$—were available merely through a few iterations, and the computations of merit functions are quite saving. Complex eigenvalue problems are more difficult than real eigenvalue problems. In Theorems 3 and 4, we explored two normalization equations for complex eigenvalue problems. For real eigenvalue problems, two regularization methods were constructed, which upon combination with the derivative-free fixed-point Newton method can find the real eigenvalues about ten times more quickly than using the 1D GSSA in the simple methods. Moreover, the combination of Newton’s iterative technique with two regularized normalization methods confirmed that using suitable values for the regularization parameters can enhance the convergence speed and also the accuracy of the solutions. Compared to the two Newton methods, which need to invert the Jacobian matrices at each iteration step, the derivative-free fixed-point Newton method only solves a scalar equation; hence, it can save much computational time and have the same good convergence rate as the two Newton methods, for which the number of iterations for the tested examples are few.