可对角化矩阵特征值分解扰动问题的快速求解方法

A Fast Computation Method for Eigenvalue Decomposition Perturbation Problems of Diagonalizable Matrix

摘要

针对特征值扰动计算的传统方法收敛速度慢的问题，提出了一种求解特征值扰动问题的快速迭代算法.首先，通过矩阵变换将初始矩阵的特征值扰动问题转化为对角矩阵的特征值扰动问题.然后，提出了一种快速迭代算法求解扰动参数，同时对算法的收敛性进行分析，并将其与基于摄动级数展开法导出的方法进行对比. 再次，采用逐一求解特征值并进行矩阵降阶的策略，有效降低运算量.最后，通过2个算例分别展示算法的计算过程及其在结构模态参数追踪方面的应用效果.

Abstract

A fast iterative algorithm for solving the eigenvalue perturbation problem is proposed in this paper for solving the problem of slow convergence of traditional methods for eigenvalue perturbation calculation. Firstly， the eigenvalue perturbation problem of the initial matrix is transformed into the eigenvalue perturbation problem of the diagonal matrix by matrix transformation. Then， a fast iterative algorithm is proposed to solve the perturbation parameter. The convergence of the algorithm is analyzed and compared with the method derived based on the perturbation series expansion method， and the strategy of solving the eigenvalues one by one and reducing the order of the matrix is adopted to effectively reduce the computation cost. Finally， two examples are used to show the calculation process of the algorithm and its application in the tracking of modal parameters of vibration structures.

Graphical abstract

关键词

特征值分解 / 特征值扰动 / 摄动级数展开法 / 可对角化矩阵 / 收敛性分析

Key words

eigenvalue decomposition / eigenvalue perturbation / perturbation series expansion method / diagonalizable matrix / convergence analysis

特征值分解不仅是矩阵理论的重要部分，也在工程计算中起着重要作用^［1-2］. 特征值和特征向量具有明确的物理意义，例如，在结构振动分析领域，特征值对应着结构的频率和阻尼比，而特征向量对应着结构的振型^［3-6］. 现代计算机技术的发展使得求解高阶矩阵的特征值问题成为可能，自20世纪50年代以来，科学家们提出了大量的数值方法计算矩阵的特征值分解^［1］，其中包括幂法、QR分解法、Jacobi方法等基础方法，以及基于Krylov子空间算法的Arnoldi算法和Lanczos算法^［1-2］. 当前实用的特征值求解器都是基于上述方法开发的^［7-8］.

谢和虎^［9-10］将大型矩阵的特征值问题转化为多重网格有限元空间上的边值问题，其核心在于通过降维减少运算量. 张宁等^［8］提出了一种广义共轭梯度算法用于求解特征值问题，并开发了可以并行计算的软件包. 尽管这些方法能有效地进行特征值分解计算，但在面对超大数值模型时，单次计算成本仍然很高. 结构优化设计和动态模型修正需要连续观察特征值变化，若模型每次修改都重新计算特征值分解，则运算量将无法承受^［11］. 另外，在结构实际运营期间，其物理参数（如刚度、质量等）有可能不断变化，要及时计算和追踪特征值和特征向量的变化^［12］. 然而，随着分析模型愈发庞大，矩阵的阶数显著增加，特征值分解计算的效率问题日渐突出.

当矩阵元素发生微小变化时，其特征值分解相应产生扰动. 此时，初始矩阵的已知特征值分解结果可作为特征值扰动计算的先验信息.基于矩阵摄动理论的方法是特征值扰动计算的重要方法，通过将特征值和特征向量进行幂级数展开以逼近精确结果^{［5， 13-18］}. 特征值分解的一阶扰动为线性项，在灵敏度分析中有重要作用，使用较为广泛^［19-20］. 然而，忽略高阶级数会造成累计误差，而采用高阶级数又会降低计算效率. 若摄动展开级数发散，则会导致计算错误. 为开发稳健及高效的特征值扰动计算方法，有必要优化算法收敛性能.

本文旨在提出一种求解可对角化矩阵受扰动后的特征值分解问题的高效算法. 首先将初始矩阵的特征值扰动问题转化为对角矩阵的特征值扰动问题. 引入一个可逆矩阵，使受扰动的矩阵相似于一个分块上三角矩阵，并用迭代算法求解该可逆矩阵的未知参数. 随后，逐个计算扰动特征值并不断降低待处理矩阵的阶数. 通过计算各特征值对应的收敛指标，采用优先处理孤立特征值并降阶的策略来加快计算.

1 特征值扰动问题

设初始矩阵\mathit{A}\in {\mathbf{ℂ}}^{n\times n}，存在n个孤立的特征值，且可作如下特征值分解：

\mathit{A}=\mathit{P}\mathit{\Lambda }{\mathit{P}}^{-\mathrm{1}}（1）

式中：\mathit{\Lambda }为对角矩阵，其对角元素为特征值；\mathit{P}矩阵的各列为特征向量. 在\mathit{\Lambda }和\mathit{P}皆可逆时，称初始矩阵\mathit{A}可对角化于\mathit{\Lambda }. 当初始矩阵受加法扰动{\mathit{A}}_{\mathrm{p}}时，相应扰动后特征值分解为：

\mathit{A}+{\mathit{A}}_{\mathrm{p}}=\tilde{\mathit{P}}\tilde{\mathit{\Lambda }}{\tilde{\mathit{P}}}^{-\mathrm{1}}（2）

此即关于矩阵\mathit{A}的特征值扰动问题，目标是在已知\mathit{\Lambda }、\mathit{P}和{\mathit{P}}^{-\mathrm{1}}的条件下求解\tilde{\mathit{\Lambda }}、\tilde{\mathit{P}}和{\tilde{\mathit{P}}}^{-\mathrm{1}}. 令：

\mathit{R}={\mathit{P}}^{-\mathrm{1}}{\mathit{A}}_{\mathrm{p}}\mathit{P}（3）

即可将一般矩阵的特征值扰动问题转化为对角矩阵的特征值扰动问题：

\mathit{\Lambda }+\mathit{R}=\mathit{Q}\tilde{\mathit{\Lambda }}{\mathit{Q}}^{-\mathrm{1}}（4）

式中：\mathit{Q}={\mathit{P}}^{-\mathrm{1}}\tilde{\mathit{P}}. 求出对角矩阵的特征值分解扰动结果，就可得扰动后特征值\tilde{\mathit{\Lambda }}，并可计算出初始矩阵\mathit{A}的扰动后特征向量：

\tilde{\mathit{P}}=\mathit{P}\mathit{Q}（5）

特征值扰动问题求解难度受\mathit{\Lambda }、\mathit{P}和{\mathit{A}}_{\mathrm{p}}的共同影响. 根据特征值分解灵敏度分析，当\mathit{\Lambda }中存在密集分布的特征值，其特征值分解对扰动非常敏感，可能导致摄动展开级数发散^［21］. 因此有必要分析\mathit{\Lambda }中特征值（在复平面上）的分布，预先挑选出易于计算的孤立特征值和不易计算的密集特征值. 对于扰动计算无法快速收敛的特征值，可采用分步计算方法^{［4， 22］}. 分步计算的思路是把扰动项\mathit{R}分为M步，取步长为\mu =\mathrm{1}/M，每步在当前矩阵基础上添加扰动项\mu \mathit{R}并进行扰动计算.

2 特征值扰动计算方法

在特征值扰动问题中，设扰动量为一阶小量，扰动后的特征值总是接近原矩阵的特征值.根据这一特点，本节给出单个特征值扰动的计算方法.将 式（4）中的\mathit{\Lambda }和\mathit{R}写作分块矩阵形式：

\mathit{\Lambda }=\left[\begin{array}{cc}\lambda & \mathrm{0}\\ \mathrm{0}& {\mathit{\Lambda }}_{\mathrm{d}}\end{array}\right],\mathit{R}=\left[\begin{array}{cc}r& {\mathit{r}}_{\mathrm{h}}\\ {\mathit{r}}_{\mathrm{c}}& {\mathit{R}}_{\mathrm{d}}\end{array}\right]（6）

式中：\lambda为\mathit{\Lambda }的第一个对角元素，{\mathit{\Lambda }}_{\mathrm{d}}为n-\mathrm{1}阶对角矩阵；r为\mathit{R}的第一个对角元素；{\mathit{r}}_{\mathrm{c}}为n-\mathrm{1}维列向量，{\mathit{r}}_{\mathrm{h}}为n-\mathrm{1}维行向量，{\mathit{R}}_{\mathrm{d}}为n-\mathrm{1}阶方阵. 引入可逆矩阵：

{\mathit{Q}}_{\mathrm{1}}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathit{x}& {\mathit{I}}_{n-\mathrm{1}}\end{array}\right],{\mathit{Q}}_{\mathrm{1}}^{-\mathrm{1}}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ -\mathit{x}& {\mathit{I}}_{n-\mathrm{1}}\end{array}\right]（7）

若\mathit{x}满足方程：

{\mathit{r}}_{\mathrm{c}}-\lambda \mathit{x}-r\mathit{x}+{\mathit{\Lambda }}_{\mathrm{d}}\mathit{x}+{\mathit{R}}_{\mathrm{d}}\mathit{x}-\mathit{x}{\mathit{r}}_{\mathrm{h}}\mathit{x}=\mathrm{0}（8）

结合式（6）和式（7）可得：

{\mathit{Q}}_{\mathrm{1}}^{-\mathrm{1}}\left(\mathit{\Lambda }+\mathit{R}\right){\mathit{Q}}_{\mathrm{1}}=\left[\begin{array}{cc}\lambda +r+{\mathit{r}}_{\mathrm{h}}\mathit{x}& {\mathit{r}}_{\mathrm{h}}\\ \mathrm{0}& {\mathit{\Lambda }}_{\mathrm{d}}+{\mathit{R}}_{\mathrm{d}}-\mathit{x}{\mathit{r}}_{\mathrm{h}}\end{array}\right]（9）

进而可得：

\mathit{\Lambda }+\mathit{R}={\mathit{Q}}_{\mathrm{1}}\left[\begin{array}{cc}\lambda +r+{\mathit{r}}_{\mathrm{h}}\mathit{x}& {\mathit{r}}_{\mathrm{h}}\\ \mathrm{0}& {\mathit{\Lambda }}_{\mathrm{d}}+{\mathit{R}}_{\mathrm{d}}-\mathit{x}{\mathit{r}}_{\mathrm{h}}\end{array}\right]{\mathit{Q}}_{\mathrm{1}}^{-\mathrm{1}}（10）

即\mathit{\Lambda }+\mathit{R}相似于一个分块上三角矩阵. 只要求解出\mathit{x}，就得到了一个扰动特征值：

\tilde{\lambda }=\lambda +r+{\mathit{r}}_{\mathrm{h}}\mathit{x}（11）

因此，可逆矩阵{\mathit{Q}}_{\mathrm{1}}的引入使特征值扰动问题转化为了关于\mathit{x}的非线性方程求解问题. 为此，本文提出一种新的求解该方程的迭代算法，并将其与经典的摄动级数展开法导出的结果进行对比.

2.1 迭代算法

在小扰动情况下，式（8）中关于\mathit{x}的非线性项\mathit{x}{\mathit{r}}_{\mathrm{h}}\mathit{x}为三阶小量，因此\left[\left(\lambda +r\right)\mathit{I}-{\mathit{\Lambda }}_{\mathrm{d}}-\mathit{R}\right]\mathit{x}\approx {\mathit{r}}_{\mathrm{c}}. 令

\mathit{G}={\left[\left(\lambda +r\right)\mathit{I}-{\mathit{\Lambda }}_{\mathrm{d}}-{\mathit{R}}_{\mathrm{d}}\right]}^{-\mathrm{1}}（12）

构造如下迭代算法：

\left\{\begin{array}{l}{\mathit{x}}_k=\mathit{G}{\mathit{r}}_{\mathrm{c}},k=\mathrm{1}\\ {\mathit{x}}_k={\mathit{x}}_{\mathrm{1}}-\mathit{G}\left({\mathit{r}}_{\mathrm{h}}{\mathit{x}}_{k-\mathrm{1}}\right){\mathit{x}}_{k-\mathrm{1}},k>\mathrm{1}\end{array}\right.（13）

式中：k=\mathrm{1},\mathrm{2},\cdots ,k_{\mathrm{m}\mathrm{a}\mathrm{x}},\mathrm{其中}{k}_{\mathrm{m}\mathrm{a}\mathrm{x}}为迭代次数上限.\mathit{G}矩阵在迭代过程中保持不变，仅需计算一次.根据 式（11），第k步得到的扰动特征值为：

\tilde{\lambda }(k)=\lambda +r+{\mathit{r}}_{\mathrm{h}}{\mathit{x}}_k（14）

为分析该迭代算法的收敛性，设每步迭代计算 \mathit{x}的改变量为：

\left\{\begin{array}{l}{\mathit{y}}_k={\mathit{x}}_k,k=\mathrm{1}\\ {\mathit{y}}_k={\mathit{x}}_k-{\mathit{x}}_{k-\mathrm{1}},k>\mathrm{1}\end{array}\right.（15）

可得

{\mathit{y}}_{k+\mathrm{1}}=\mathit{G}{\mathit{\Gamma }}_k{\mathit{y}}_k,k>\mathrm{1}（16）

式中：{\mathit{\Gamma }}_k=-{\mathit{r}}_{\mathrm{h}}^{\mathrm{T}}{\mathit{x}}_k\mathit{I}-{\mathit{x}}_{k-\mathrm{1}}{\mathit{r}}_{\mathrm{h}}^{\mathrm{T}}. 理论上，只要向量{\mathit{y}}_k的范数不断减小，即可保证算法收敛. 根据诱导范数性质，有

‖{\mathit{y}}_{k+\mathrm{1}}‖\le {‖\mathit{G}{\mathit{\Gamma }}_k‖}_{\mathrm{2}}\bullet ‖{\mathit{y}}_k‖\le {‖\mathit{G}‖}_{\mathrm{2}}\bullet {‖{\mathit{\Gamma }}_k‖}_{\mathrm{2}}\bullet ‖{\mathit{y}}_k‖（17）

式中：‖\bullet ‖表示欧氏范数（即向量长度），{‖\bullet ‖}_{\mathrm{2}}表示矩阵2范数. 只要控制\mathit{G}{\mathit{\Gamma }}_k矩阵2范数小于1，即可保证算法收敛. \mathit{G}的2范数受到特征值分布的影响，在小扰动情况下：

{‖\mathit{G}‖}_{\mathrm{2}}\approx \underset{i\ne \mathrm{1}}{\mathrm{m}\mathrm{a}\mathrm{x}}\frac{\mathrm{1}}{\left|\mathit{\Lambda }(i,i)-\lambda \right|}（18）

即{‖\mathit{G}‖}_{\mathrm{2}}取决于\lambda与其他特征值的最小间距，这表明较为孤立的特征值的扰动计算更容易收敛. 由于在算法执行时并无必要计算每次迭代对应的范数，仅需考虑{\mathit{\Gamma }}_{\mathrm{2}}对应的收敛指标，该指标仅与\mathit{G}、{\mathit{r}}_{\mathrm{c}}和{\mathit{r}}_{\mathrm{h}}有关. 因此，结合式（18），定义以下便于计算的收敛指标：

\rho =\mathrm{2}‖{\mathit{r}}_{\mathrm{c}}‖\bullet ‖{\mathit{r}}_{\mathrm{h}}‖\bullet \underset{i\ne \mathrm{1}}{\mathrm{m}\mathrm{a}\mathrm{x}}\frac{\mathrm{1}}{{\left|\mathit{\Lambda }(i,i)-\lambda \right|}^{\mathrm{2}}}（19）

该指标可在调用迭代程序之前进行计算，用于预测算法的收敛性. \rho值越小，则迭代算法越容易收敛，因此可根据\rho值将待处理的特征值进行排序，优先处理易收敛的特征值. 按照指标\rho由小到大排列，确定特征值处理顺序向量\mathit{p}，并将调整了索引顺序的矩阵\mathit{\Lambda }(\mathit{p},\mathit{p})和\mathit{R}\left(\mathit{p},\mathit{p}\right)输入迭代算法. 而在迭代子程序内部通过监测特征值的变化以判断收敛性，一旦算法达到预定精度或迭代次数达到上限，即停止迭代.

根据式（14），可知各次迭代对应的特征值改变量为：

{\Delta }_{\mathrm{1}}=r+{\mathit{r}}_{\mathrm{h}}{\mathit{x}}_{\mathrm{1}}

{\Delta }_k={\mathit{r}}_{\mathrm{h}}{\mathit{y}}_k,k>\mathrm{1}（20）

本文将{\mathrm{\Delta }}_k作为精度判断指标，当时停止迭代. 算例中取ϵ=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}，迭代步数上限为20. 该算法的流程图如图1所示，当算法收敛时输出满足精度要求的{\tilde{\lambda }}_k和{\mathit{x}}_k;算法不收敛时则输出\tilde{\lambda }=\lambda和\mathit{x}=\mathrm{0}，此时{\mathit{Q}}_{\mathrm{1}}为单位矩阵. 在该迭代算法不收敛时，可采用分步法进行扰动计算\mathit{x}.

2.2 摄动级数展开法

根据摄动理论^［11］，将扰动矩阵\mathit{\Lambda }+\epsilon \mathit{R}对应的\mathit{x}展开为\epsilon的幂级数：

\mathit{x}=\epsilon {\mathit{y}}_{\mathrm{1}}+{\epsilon }^{\mathrm{2}}{\mathit{y}}_{\mathrm{2}}+{\epsilon }^{\mathrm{3}}{\mathit{y}}_{\mathrm{3}}+\cdots（21）

式中\epsilon是一个小参数. 将该摄动解代入式（9）可得：

\left(\lambda \mathit{I}-{\mathit{\Lambda }}_{\mathrm{d}}\right)\mathit{x}=\epsilon {\mathit{r}}_{\mathrm{c}}-\epsilon \left({\mathit{R}}_{\mathrm{d}}-r\mathit{I}\right)\mathit{x}-\epsilon \mathit{x}{\mathit{r}}_{\mathrm{h}}\mathit{x}（22）

通过比较\epsilon的同次幂系数，可归纳出{\mathit{y}}_k的递推公式：

{\mathit{y}}_{\mathrm{1}}={\left(\lambda \mathit{I}-{\mathit{\Lambda }}_{\mathrm{d}}\right)}^{-\mathrm{1}}{\mathit{r}}_{\mathrm{c}}

{\mathit{y}}_k={\left(\lambda \mathit{I}-{\mathit{\Lambda }}_{\mathrm{d}}\right)}^{-\mathrm{1}}\left[{\mathit{R}}_{\mathrm{d}}{\mathit{y}}_{k-\mathrm{1}}-{\sum }_{l=\mathrm{1}}^{k-\mathrm{1}}{∆}_{k-l}{\mathit{y}}_l\right],k>\mathrm{1}（23）

式中：{∆}_k为特征值改变量，按下式计算：

{∆}_{\mathrm{1}}=r

{∆}_k={\mathit{r}}_{\mathrm{h}}{\mathit{y}}_{k-\mathrm{1}},k>\mathrm{1}（24）

根据式（11）和式（21），可得到第k步对应的特征值结果：

\tilde{\lambda }\left(k\right)=\lambda +{∆}_{\mathrm{1}}+\cdots +{∆}_k（25）

式（21）中取\epsilon =\mathrm{1}，将各步迭代得到的{\mathit{y}}_k叠加，可得到\mathit{x}的渐近解：

\begin{array}{l}{\mathit{x}}_k={\mathit{x}}_{\mathrm{1}}+{\left(\lambda \mathit{I}-\mathit{B}\right)}^{-\mathrm{1}}\left[{\mathit{R}}_{\mathrm{d}}{\mathit{x}}_{k-\mathrm{1}}-{\sum }_{l=\mathrm{1}}^{k-\mathrm{1}}{∆}_{k-l}{\mathit{x}}_l\right],\\ k>\mathrm{2}\end{array}（26）

该公式计算结果与经典摄动法求解特征值扰动的结果一致^［11］. 高阶摄动解可以通过递归的方式计算，在计算第k阶的结果时，需要采用前k-\mathrm{1}阶的摄动解. 当所需递归步数增加，数据存储量亦将增多. 由式（23）可见，与本文提出的迭代算法类似，摄动级数展开法的收敛性也受到初始矩阵特征值分布间距的影响，距其他特征值较远的特征值对应的扰动计算更容易收敛. 在迭代算法和摄动级数展开法皆收敛时，两者都能得出准确的结果. 迭代算法每步的修正量为三阶小量，而摄动级数展开法每步修正量为二阶小量，因而本文提出的迭代算法收敛速度更快.

3 矩阵降阶方法

采用迭代算法可求解单一特征值的扰动问题，当矩阵阶数n较大时，需进行n次循环求解全部扰动后的特征值. 为提高算法效率，本节提出一种矩阵降阶的方法以提高特征值扰动的计算速度. 根据式（8）和式（10），右侧矩阵中的{\mathit{\Lambda }}_{\mathrm{d}}可作为新的对角矩阵，{\mathit{R}}_{\mathrm{d}}-\mathit{x}{\mathit{r}}_{\mathrm{h}}^{}作为新的扰动项. 因此，可采用第2节的迭代算法求解{\mathit{\Lambda }}_{\mathrm{d}}(\mathrm{1,1})对应的特征值扰动结果，并得到新的分块对角矩阵. 每次调用迭代算法后待求解矩阵降低1阶，第i次调用迭代算法求解出的{\mathit{x}}^{(i)}\in {\mathbf{ℂ}}^{n-i}. 重复该过程，可得到\mathit{\Lambda }+\mathit{R}的相似上三角矩阵：

{\mathit{T}}_{\mathrm{r}}={\mathit{V}}^{-\mathrm{1}}\left(\mathit{\Lambda }+\mathit{R}\right)\mathit{V}（27）

其中：

\mathit{V}={\prod }_{i=\mathrm{1}}^{n-\mathrm{1}}{\mathit{W}}_i,{\mathit{V}}^{-\mathrm{1}}={\prod }_{i=\mathrm{1}}^{n-\mathrm{1}}{\mathit{W}}_{n-i}^{-\mathrm{1}}（28）

{\mathit{W}}_i=\left[\begin{array}{cc}{\mathit{I}}_{i-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& {\mathit{Q}}_i\end{array}\right]， {\mathit{W}}_i^{-\mathrm{1}}=\left[\begin{array}{cc}{\mathit{I}}_{i-\mathrm{1}}& \mathrm{0}\\ \mathrm{0}& {\mathit{Q}}_i^{-\mathrm{1}}\end{array}\right]（29）

{\mathit{Q}}_i=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ {\mathit{x}}^{(i)}& {\mathit{I}}_{n-i}\end{array}\right]， {\mathit{Q}}_i^{-\mathrm{1}}=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ -{\mathit{x}}^{(i)}& {\mathit{I}}_{n-i}\end{array}\right]（30）

式中：{\mathit{I}}_{i-\mathrm{1}}和{\mathit{I}}_{n-i}分别为i-\mathrm{1}阶和n-i阶单位矩阵，当i=\mathrm{1}时，{\mathit{W}}_{\mathrm{1}}={\mathit{Q}}_{\mathrm{1}}；矩阵{\mathit{Q}}_i由第i次调用迭代算法得到. 此处列出式（28）~式（30）是便于更好地理解矩阵的结构，实际计算中无须存储{\mathit{W}}_i和{\mathit{Q}}_i矩阵. 事实上，\mathit{V}是对角线元素皆为1的下三角矩阵，且其对角线以下的元素为：

\mathit{V}\left(i+\mathrm{1}:n,i\right)={\mathit{x}}^{(i)}（31）

图2给出了循环计算多个特征值扰动问题的流程图. 计算时需调用2.1节介绍的迭代算法，每循环一次就可以更新\mathit{V}和{\mathit{V}}^{-\mathrm{1}}的部分元素，并使特征值扰动问题下降1阶. 当完成n-\mathrm{1}次循环后，可用 式（27）计算{\mathit{T}}_{\mathrm{r}}，并输出结果.

上三角矩阵{\mathit{T}}_{\mathrm{r}}的特征值分解容易计算：

{\mathit{T}}_{\mathrm{r}}=\mathit{\Psi }\tilde{\mathit{\Lambda }}{\mathit{\Psi }}^{-\mathrm{1}}（32）

事实上，矩阵{\mathit{T}}_{\mathrm{r}}的对角线元素即为其特征值，也是矩阵扰动后的特征值. 结合式（4）、式（26）及 式（31），可得

\tilde{\mathit{P}}=\mathit{P}\mathit{V}\mathit{\Psi }（33）

即可得到初始矩阵扰动后的特征值分解：

\mathit{A}+{\mathit{A}}_{\mathrm{p}}=\left(\mathit{P}\mathit{V}\mathit{\Psi }\right)\tilde{\mathit{\Lambda }}\left({\mathit{\Psi }}^{-\mathrm{1}}{\mathit{V}}^{-\mathrm{1}}{\mathit{P}}^{-\mathrm{1}}\right)（34）

对于n阶矩阵的扰动问题，若按上述方法调用n-1次迭代算法皆收敛，可精确地求解出所有特征值扰动结果. 若成功调用迭代算法i次，则解决了i个特征值的扰动计算，原矩阵的扰动问题降阶为n-i阶矩阵的扰动问题. 结合降阶方法的迭代算法在求解孤立特征值扰动问题时易收敛，同时降低问题的阶数，能有效减少计算量.

4 算例

4.1 含密集特征值矩阵的扰动分析

为验证本文所提算法，以一个包含密集特征值的4阶矩阵为例进行扰动计算. 矩阵\mathit{\Lambda }和\mathit{R}取值如下：

\mathit{\Lambda }=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[\mathrm{0.1}+\mathrm{1.0}i\mathrm{1.1}+\mathrm{0.1}i\mathrm{1.0}+\mathrm{1.0}i\mathrm{1.05}+\mathrm{1.05}i\right]

\mathit{R}=\left[\begin{array}{l}-\mathrm{0.09}+\mathrm{0.20}i\mathrm{0.12}-\mathrm{0.07}i\mathrm{0.02}+\mathrm{0.16}i\mathrm{0.11}+\mathrm{0.07}i\\ -\mathrm{0.07}-\mathrm{0.70}i\mathrm{0.41}+\mathrm{0.28}i\mathrm{0.32}-\mathrm{0.05}i\mathrm{0.13}-\mathrm{0.47}i\\ \mathrm{0.19}+\mathrm{0.02}i-\mathrm{0.31}-\mathrm{0.10}i\mathrm{0.10}-\mathrm{0.18}i\mathrm{0.04}+\mathrm{0.12}i\\ \mathrm{0.29}-\mathrm{0.31}i\mathrm{0.01}+\mathrm{0.38}i-\mathrm{0.06}-\mathrm{0.10}i-\mathrm{0.41}+\mathrm{0.06}i\end{array}\right]

扰动矩阵\mathit{R}中元素的实部和虚部数值都是从均值为0、方差为0.2的正态分布中抽样得到的. 矩阵\mathit{\Lambda }的第1、2个特征值与其他特征值在复平面上相距（相对）较远；第3、4个特征值相互邻近. 扰动前后的特征值分布如图3所示，图中Re（\lambda）和Im（\lambda）分别为复数的实部和虚部.表1给出了各特征值对应的收敛指标\rho ,{\lambda }_{\mathrm{1}}和{\lambda }_{\mathrm{2}}对应的收敛指标较小，而由于{\lambda }_{\mathrm{3}}和{\lambda }_{\mathrm{4}}距离较近，故其收敛指标较大. 本例中\mathit{p}=[\mathrm{1,2},\mathrm{3,4}]，矩阵行、列顺序不需调整.

首先用摄动级数展开法和本文提出的迭代算法分别计算特征值{\lambda }_{\mathrm{1}}的扰动结果.表2中列出了两种算法前4步迭代计算得到的\tilde{\lambda }(k)和\left|{\mathrm{\Delta }}_k\right|的结果.由 图4可知，两种算法对应的\left|{\mathrm{\Delta }}_k\right|随迭代步数增加而减小，两种算法得到的特征值变化量{\mathrm{\Delta }}_k随迭代步数增加而变小，因此两种算法都收敛.本文提出的迭代算法收敛速度更快，第4步迭代就可得到相对误差为0.01%的特征值结果；而摄动法第4步迭代得到的特征值结果还包含0.23%的相对误差.特征值相对误差计算公式为：

{\delta }_{\mathrm{E}}=\frac{\left|\tilde{\lambda }\left(k\right)-{\tilde{\lambda }}_{\mathrm{E}}\right|}{\left|{\tilde{\lambda }}_{\mathrm{E}}\right|}\times \mathrm{100}\mathrm{\%}（35）

式中：{\tilde{\lambda }}_{\mathrm{E}}为扰动后的精确特征值. 以作为收敛判断条件，本文迭代算法仅需10步即可达到精度要求，而摄动级数展开法需43步才能达到这一精度要求.

在求解出{\mathit{x}}^{(\mathrm{1})}之后，即可采用第4节中的算法继续计算第2个特征值扰动问题对应的{\mathit{x}}^{(\mathrm{2})}. 由于第1、2个特征值与其他特征值相距较远，迭代算法收敛较快. 然而在计算时{\mathit{x}}^{(\mathrm{3})}时，由于第3、第4个特征值相邻，迭代算法收敛困难. 此时受扰矩阵已降为2阶，受扰矩阵和扰动矩阵分别为：

{\mathit{\Lambda }}^{(\mathrm{3})}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left[\begin{array}{cc}\mathrm{1.0}+\mathrm{1.0}i& \mathrm{1.05}+\mathrm{1.05}i\end{array}\right]

{\mathit{R}}^{\left(\mathrm{3}\right)}=\left[\begin{array}{cc}\mathrm{0.132}\mathrm{9}+\mathrm{0.016}\mathrm{7}\mathrm{i}& \mathrm{0.297}\mathrm{3}+\mathrm{0.188}\mathrm{3}i\\ \mathrm{0.079}\mathrm{2}-\mathrm{0.122}\mathrm{0}i& -\mathrm{0.393}\mathrm{6}-\mathrm{0.089}\mathrm{2}i\end{array}\right]

利用分步方法计算该2阶矩阵扰动问题，取分步数M=\mathrm{3}.三个分步调用迭代算法得到的\left|{\mathrm{\Delta }}_k\right|变化如图5所示，结果表明特征值改变量都随迭代次数下降，各分步计算都收敛. 不过由于特征值初始间距较小，第1分步收敛速度相对较慢，需要19次迭代才能达到精度要求；后两个分步收敛速度快，分别经过4次和3次迭代即可达到精度要求.

在{\mathit{x}}^{(\mathrm{1})}、{\mathit{x}}^{(\mathrm{2})}和{\mathit{x}}^{(\mathrm{3})}都计算出后，根据式（26）~式（30）即可得到矩阵\mathit{V}、{\mathit{V}}^{-\mathrm{1}}和{\mathit{T}}_{\mathrm{r}}.进而对上三角矩阵进行特征值分解得到\mathit{\Psi }，即可实现\mathit{\Lambda }+\mathit{R}的特征值分解.限于篇幅，本文仅列出\mathit{V}的结果（因为\mathit{V}的对角线以下元素包含了{\mathit{x}}^{(\mathrm{1})}、{\mathit{x}}^{(\mathrm{2})}和{\mathit{x}}^{(\mathrm{3})}的值）：

\mathit{V}=\left[\begin{array}{cccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -\mathrm{0.147}\mathrm{7}+\mathrm{0.223}\mathrm{2}i& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ -\mathrm{0.156}\mathrm{4}-\mathrm{0.005}\mathrm{1}i& -\mathrm{0.043}\mathrm{9}-\mathrm{0.482}\mathrm{0}i& \mathrm{1}& \mathrm{0}\\ -\mathrm{0.325}\mathrm{7}+\mathrm{0.386}\mathrm{7}i& -\mathrm{0.211}\mathrm{7}+\mathrm{0.166}\mathrm{8}i& \mathrm{0.135}\mathrm{6}-\mathrm{0.222}\mathrm{1}i& \mathrm{1}\end{array}\right]

上述含密集特征值矩阵扰动问题的计算展示了算法的计算过程和算法的有效性，同时表明了本文迭代算法的收敛速度快于传统的摄动级数展开法.

4.2 具有时变参数的结构模态分析

如图6所示的平面桁架结构模型，包含20个节点及46根杆件. 采用集中质量模型，各节点集中质量为m=\mathrm{10}\mathrm{k}\mathrm{g}，不考虑杆件重量.各杆件截面积皆为A=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}{\mathrm{m}}^{\mathrm{2}}，材料弹性模量为E=\mathrm{206}\mathrm{G}\mathrm{P}\mathrm{a}. 竖向和横向杆件长1.0 m，斜向撑杆长\sqrt[]{\mathrm{2}}m，桁架总跨度为9 m.桁架节点1受水平和竖向约束，节点19受竖向约束.

利用有限元法可获得该桁架结构的质量矩阵\mathit{M}和整体刚度矩阵\mathit{K}，为便于分析，采用瑞利阻尼假定\mathit{C}=\alpha \mathit{M}+\beta \mathit{K}，取\alpha =\mathrm{0.1},\beta =\mathrm{0.000}\mathrm{1}. 利用结构参数可得到结构的状态空间模型系统矩阵：

{\mathit{A}}_c=\left[\begin{array}{cc}\mathrm{0}& \mathit{I}\\ -{\mathit{M}}^{-\mathrm{1}}\mathit{K}& -{\mathit{M}}^{-\mathrm{1}}\mathit{C}\end{array}\right]（36）

对于欠阻尼系统，系统矩阵的特征值为复数，包含了结构模态频率和阻尼信息. 对{\mathit{A}}_c进行特征值分解得到的特征值为：

{\lambda }_{\mathrm{c}i},{\lambda }_{\mathrm{c}i}^{*}=-{\xi }_i{\omega }_i\pm i{\omega }_i\sqrt[]{\mathrm{1}-{\xi }_i^{\mathrm{2}}}（37）

式中：{\lambda }_{\mathrm{c}i}^{*}为{\lambda }_{\mathrm{c}i}的共轭复数，这说明系统矩阵的特征值以共轭形式成对出现；{\omega }_i和{\xi }_i分别表示第i阶模态振动的固有频率和阻尼比.表3中列出了该结构的前5阶模态频率、阻尼比以及系统矩阵的特征值.

实际工程经常出现结构的模态参数随时间发生变化的现象，为追踪特征值和特征向量的变化，若采用对不同时间点对应的系统矩阵进行特征值分解的方法，一方面要进行多次特征值求解，另一方面无法保证不同时间点系统矩阵特征值的顺序一一对应. 采用本文算法进行结构模态参数追踪，只需要在初始时刻计算一次系统矩阵的特征值分解，再依次对不同时间的系统矩阵求解特征值扰动问题.对于图6中的桁架结构，假设连接节点1和3的杆件刚度随时间发生扰动，其值为

∆E\left(\mathit{t}\right)=E\bullet \mathrm{0.2}\mathrm{s}\mathrm{i}\mathrm{n}\left(\mathrm{2}\mathrm{\pi }\mathit{t}\right)（38）

式中：\mathit{t}=∆t\bullet [\mathrm{1,2},\cdots ,\mathrm{50}]，∆t=\mathrm{0.02}，总共分析50个时间步.此时，可用本文算法追踪振动系统特征值的扰动.图7给出了本文算法得出的不同时间桁架结构前5阶模态频率受时变参数影响的结果.可见，杆件的刚度变化会导致多阶模态频率发生变化，且总体上各阶频率随杆件刚度增加而增加.同时，对比各阶频率变化幅度，可发现各阶频率对特定杆件刚度变化的敏感性不同，这与结构动力参数敏感性分析的结论相符. 上述结果表明，采用本文算法可高效地跟踪计算时变结构的模态参数.

5 结 论

本文针对摄动级数展开法求解特征值扰动问题收敛较慢的问题，提出了一种可对角化矩阵特征值分解扰动问题的快速求解算法. 首先，将一般矩阵的扰动问题转化为对角矩阵的扰动问题，并利用一对可逆矩阵计算扰动后矩阵的相似上三角矩阵. 其次，给出了迭代算法的递推公式及详细的流程图，并与摄动级数展开法进行了比较. 再次，提出了可快速计算的收敛性评价指标，以便挑选并优先处理孤立特征值. 本文算法的另一特点在于结合了矩阵降阶策略，使需要采用分步计算求解的特征值的扰动问题的计算量大幅降低. 算例结果展示了所提算法收敛速度快于传统摄动级数展开法；并利用所提算法计算了一榀具有时变刚度构件的桁架结构的时变模态频率，验证了算法在追踪时变结构模态参数方面的应用潜力.

参考文献

[1] GOLUB G H， VAN DER VORST H A ．Eigenvalue computation in the 20th century［J］．Journal of Computational and Applied Mathematics，2000，123（1/2）：35-65．

[2] 朱晓临 ．数值分析［M］．2版．合肥：中国科学技术大学出版社，2014．

[3] ZHU X L ．Numerical analysis［M］．2nd ed．Hefei：University of Science and Technology of China Press， 2014．（in Chinese）

[4] 王树青，田晓洁 ．结构振动测试与模态识别［M］．青岛：中国海洋大学出版社，2021．

[5] WANG S Q， TIAN X J ．Structural vibration testing and modal identification［M］．Qingdao：China Ocean University Press，2021．（in Chinese）

[6] LÜ Z H， SHAO C， FENG Z D ．High accuracy step-by-step perturbation method of the generalized eigenvalue problem and its application to the parametric history analysis of structural vibration characteristics［J］．Journal of Sound and Vibration，1993，164（3）：459-469．

[7] YU L， CHENG L， YAM L H，et al ．Application of eigenvalue perturbation theory for detecting small structural damage using dynamic responses［J］．Composite Structures，2007，78（3）：402-409．

[8] CHEN S H， MA L， MENG G W，et al ．An efficient method for evaluating the natural frequencies of structures with uncertain-but-bounded parameters［J］．Computers & Structures，2009， 87（9/10）：582-590．

[9] SAAD Y ．Numerical methods for large eigenvalue problems［M］．Rev．ed．Philadelphia，PA：Society for Industrial and Applied Mathematics，2011．

[10] 张宁，李瑜，谢和虎，等 ．一种求解特征值问题的广义共轭梯度算法［J］．中国科学（数学），2021，51（8）：1297-1320．

[11] ZHANG N， LI Y， XIE H H，et al ．A generalized conjugate gradient method for eigenvalue problems［J］．Scientia Sinica （Mathematica），2021，51（8）：1297-1320．（in Chinese）

[12] 谢和虎 ．非线性特征值问题的多重网格算法［J］．中国科学（数学），2015，45（8）：1193-1204．

[13] XIE H H ．A multigrid method for nonlinear eigenvalue problems［J］．Scientia Sinica （Mathematica），2015，45（8）：1193-1204．（in Chinese）

[14] 谢和虎 ．子空间扩展算法及其应用［J］．数值计算与计算机应用，2020，41（3）：169-191．

[15] XIE H H ．An augmented subspace method and its applications［J］．Journal on Numerical Methods and Computer Applications，2020，41（3）：169-191．（in Chinese）

[16] 邱志平，姜南 ．标准特征值问题扰动分析的精确方法［J］．科学通报，2017，62（）：3400-3408．

[17] QIU Z P， JIANG N ．An exact solution method for perturbation analysis of standard eigenvalue problems［J］．Chinese Science Bulletin，2017，62（Sup.2）：3400-3408．（in Chinese）

[18] XIA Y Y， FRISWELL M ．Efficient solution of the fuzzy eigenvalue problem in structural dynamics［J］．Engineering Computations，2014，31（5）：864-878．

[19] 吴志刚，高强 ．Hamilton系统特征值问题的摄动方法及其应用［J］．振动工程学报，2004，17（1）：7-10．

[20] WU Z G， GAO Q ．Perturbation method for eigenvalue problems of Hamiltonian systems with applications［J］．Journal of Vibration Engineering，2004，17（1）：7-10．（in Chinese）

[21] CHEN S H， WU X M， YANG Z J ．Eigensolution reanalysis of modified structures using epsilon-algorithm［J］．International Journal for Numerical Methods in Engineering，2006，66（13）：2115-2130．

[22] YANG X W， CHEN S H， WU B S ．Eigenvalue reanalysis of structures using perturbations and padé approximation［J］．Mechanical Systems and Signal Processing，2001，15（2）：257-263．

[23] LIU J K， XU W H， CAI C W ．A universal matrix perturbation technique for complex modes［J］．Applied Mathematics and Mechanics，2001，22（3）：361-367．

[24] 刘济科，徐伟华，蔡承武 ．一种通用的复模态矩阵摄动法［J］．应用数学和力学，2001，22（3）：314-320．

[25] LIU J K， XU W H， CAI C W ．A universal matrix perturbation technique for complex modes［J］．Applied Mathematics and Mechanics，2001，22（3）：314-320．（in Chinese）

[26] 刘济科，孙慧 ．广义特征值摄动问题的一种改进的通用方法［J］．暨南大学学报（自然科学与医学版），2005，26（1）：39-42．

[27] LIU J K， SUN H ．An improved method for general perturbation eigenvalue problems［J］．Journal of Jinan University （Natural Science & Medicine Edition），2005，26（1）：39-42．（in Chinese）

[28] GREENBAUM A， LI R C， OVERTON M L ．First-order perturbation theory for eigenvalues and eigenvectors［J］．SIAM Review，2020，62（2）：463-482．

[29] PANDA S， TRIPURA T， HAZRA B ．First-order error-adapted eigen perturbation for real-time modal identification of vibrating structures［J］．Journal of Vibration and Acoustics，2021，143（5）：051001．

[30] 孙继广 ．矩阵扰动分析［M］．2版．北京：科学出版社，2001．

[31] SUN J G ．Matrix perturbation analysis［M］．2nd ed．Beijing：Science Press，2001．（in Chinese）

[32] GÜNEL S， ZORAL E Y ．Parametric history analysis of resonance problems via step-by-step eigenvalue perturbation technique［J］．IET Microwaves，Antennas & Propagation，2010，4（4）：466．

基金资助

国家自然科学基金资助项目(52178283)

国家自然科学基金资助项目(52378298)

National Natural ScienceFoundation of China(52178283)

National Natural ScienceFoundation of China(52378298)

安徽省杰出青年基金项目(2208085J20)

Natural Science Fund for Distinguished Young Scholars of Anhui Province(2208085J20)