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

胡志祥 1，杨其东 1，黄潇 2，贺文宇 1?; HU Zhixiang1，YANG Qidong1，HUANG Xiao2，HE Wenyu1?

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

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

- ORCID：
胡志祥 ¹
- ORCID：
杨其东 ¹
- ORCID：
黄潇 ²
- ORCID：
贺文宇 ¹
✉

1. 合肥工业大学土木与水利工程学院,安徽合肥 230009； 2. 安徽建筑大学土木工程学院,安徽合肥 230601

中图分类号： O324

最近更新：2024-07-29

DOI： 10.16339/j.cnki.hdxbzkb.2024073

摘要

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

关键词

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

特征值分解不仅是矩阵理论的重要部分，也在工程计算中起着重要作用^［

1-2］. 特征值和特征向量具有明确的物理意义，例如，在结构振动分析领域，特征值对应着结构的频率和阻尼比，而特征向量对应着结构的振型^{［参考文献 3-6}3-6］. 现代计算机技术的发展使得求解高阶矩阵的特征值问题成为可能，自20世纪50年代以来，科学家们提出了大量的数值方法计算矩阵的特征值分解^{［参考文献 1

百度学术}1］，其中包括幂法、QR分解法、Jacobi方法等基础方法，以及基于Krylov子空间算法的Arnoldi算法和Lanczos算法^{［参考文献 1-2}1-2］. 当前实用的特征值求解器都是基于上述方法开发的^{［参考文献 7-8}7-8］.

谢和虎^［

9-10］将大型矩阵的特征值问题转化为多重网格有限元空间上的边值问题，其核心在于通过降维减少运算量. 张宁等^{［参考文献 8

百度学术}8］提出了一种广义共轭梯度算法用于求解特征值问题，并开发了可以并行计算的软件包. 尽管这些方法能有效地进行特征值分解计算，但在面对超大数值模型时，单次计算成本仍然很高. 结构优化设计和动态模型修正需要连续观察特征值变化，若模型每次修改都重新计算特征值分解，则运算量将无法承受^{［参考文献 11

百度学术}11］. 另外，在结构实际运营期间，其物理参数（如刚度、质量等）有可能不断变化，要及时计算和追踪特征值和特征向量的变化^{［参考文献 12

百度学术}12］. 然而，随着分析模型愈发庞大，矩阵的阶数显著增加，特征值分解计算的效率问题日渐突出.

当矩阵元素发生微小变化时，其特征值分解相应产生扰动. 此时，初始矩阵的已知特征值分解结果可作为特征值扰动计算的先验信息.基于矩阵摄动理论的方法是特征值扰动计算的重要方法，通过将特征值和特征向量进行幂级数展开以逼近精确结果^［

5， 13-18］. 特征值分解的一阶扰动为线性项，在灵敏度分析中有重要作用，使用较为广泛^{［参考文献 19-20}19-20］. 然而，忽略高阶级数会造成累计误差，而采用高阶级数又会降低计算效率. 若摄动展开级数发散，则会导致计算错误. 为开发稳健及高效的特征值扰动计算方法，有必要优化算法收敛性能.

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

1 特征值扰动问题

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

A = P Λ P^{- 1}

（1）

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

A + A_{p} = \tilde{P} \tilde{Λ} {\tilde{P}}^{- 1}

（2）

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

R = P^{- 1} A_{p} P

（3）

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

Λ + R = Q \tilde{Λ} Q^{- 1}

（4）

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

\tilde{P} = P Q

（5）

特征值扰动问题求解难度受 $Λ$ 、 $P$ 和 $A_{p}$ 的共同影响. 根据特征值分解灵敏度分析，当 $Λ$ 中存在密集分布的特征值，其特征值分解对扰动非常敏感，可能导致摄动展开级数发散^［

21］. 因此有必要分析

Λ

中特征值（在复平面上）的分布，预先挑选出易于计算的孤立特征值和不易计算的密集特征值. 对于扰动计算无法快速收敛的特征值，可采用分步计算方法^{［参考文献 4

百度学术}4， 22］. 分步计算的思路是把扰动项

R

分为

M

步，取步长为

μ = 1 / M

，每步在当前矩阵基础上添加扰动项

μ R

并进行扰动计算.

2 特征值扰动计算方法

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

Λ = [\begin{matrix} λ & 0 \\ 0 & Λ_{d} \end{matrix}], R = [\begin{matrix} r & r_{h} \\ r_{c} & R_{d} \end{matrix}]

（6）

式中： $λ$ 为 $Λ$ 的第一个对角元素， $Λ_{d}$ 为 $n - 1$ 阶对角矩阵； $r$ 为 $R$ 的第一个对角元素； $r_{c}$ 为 $n - 1$ 维列向量， $r_{h}$ 为 $n - 1$ 维行向量， $R_{d}$ 为 $n - 1$ 阶方阵. 引入可逆矩阵：

Q_{1} = [\begin{matrix} 1 & 0 \\ x & I_{n - 1} \end{matrix}], Q_{1}^{- 1} = [\begin{matrix} 1 & 0 \\ - x & I_{n - 1} \end{matrix}]

（7）

若 $x$ 满足方程：

r_{c} - λ x - r x + Λ_{d} x + R_{d} x - x r_{h} x = 0

（8）

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

Q_{1}^{- 1} (Λ + R) Q_{1} = [\begin{matrix} λ + r + r_{h} x & r_{h} \\ 0 & Λ_{d} + R_{d} - x r_{h} \end{matrix}]

（9）

进而可得：

Λ + R = Q_{1} [\begin{matrix} λ + r + r_{h} x & r_{h} \\ 0 & Λ_{d} + R_{d} - x r_{h} \end{matrix}] Q_{1}^{- 1}

（10）

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

\tilde{λ} = λ + r + r_{h} x

（11）

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

2.1 迭代算法

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

G = {[(λ + r) I - Λ_{d} - R_{d}]}^{- 1}

（12）

构造如下迭代算法：

\{\begin{array}{l} x_{k} = G r_{c}, k = 1 \\ x_{k} = x_{1} - G (r_{h} x_{k - 1}) x_{k - 1}, k > 1 \end{array}

（13）

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

\tilde{λ} (k) = λ + r + r_{h} x_{k}

（14）

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

\{\begin{array}{l} y_{k} = x_{k}, k = 1 \\ y_{k} = x_{k} - x_{k - 1}, k > 1 \end{array}

（15）

可得

y_{k + 1} = G Γ_{k} y_{k}, k > 1

（16）

式中： $Γ_{k} = - r_{h}^{T} x_{k} I - x_{k - 1} r_{h}^{T}$ . 理论上，只要向量 $y_{k}$ 的范数不断减小，即可保证算法收敛. 根据诱导范数性质，有

‖y_{k + 1}‖ \leq {‖G Γ_{k}‖}_{2} \cdot ‖y_{k}‖ \leq {‖G‖}_{2} \cdot {‖Γ_{k}‖}_{2} \cdot ‖y_{k}‖

（17）

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

{‖G‖}_{2} \approx \underset{i \neq 1}{m a x} \frac{1}{|Λ (i, i) - λ|}

（18）

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

ρ = 2 ‖r_{c}‖ \cdot ‖r_{h}‖ \cdot \underset{i \neq 1}{m a x} \frac{1}{{|Λ (i, i) - λ|}^{2}}

（19）

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

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

Δ_{1} = r + r_{h} x_{1}

Δ_{k} = r_{h} y_{k}, k > 1

（20）

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

图1 迭代算法流程图

Fig.1 Flowchart of the iterative algorithm

2.2 摄动级数展开法

根据摄动理论^［

11］，将扰动矩阵

Λ + ε R

对应的

x

展开为

ε

的幂级数：

x = ε y_{1} + ε^{2} y_{2} + ε^{3} y_{3} + \dots

（21）

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

(λ I - Λ_{d}) x = ε r_{c} - ε (R_{d} - r I) x - ε x r_{h} x

（22）

通过比较 $ε$ 的同次幂系数，可归纳出 $y_{k}$ 的递推公式：

y_{1} = {(λ I - Λ_{d})}^{- 1} r_{c}

y_{k} = {(λ I - Λ_{d})}^{- 1} [R_{d} y_{k - 1} - \sum_{l = 1}^{k - 1} ∆_{k - l} y_{l}], k > 1

（23）

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

∆_{1} = r

∆_{k} = r_{h} y_{k - 1}, k > 1

（24）

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

\tilde{λ} (k) = λ + ∆_{1} + \dots + ∆_{k}

（25）

式（21）中取 $ε = 1$ ，将各步迭代得到的 $y_{k}$ 叠加，可得到 $x$ 的渐近解：

\begin{array}{l} x_{k} = x_{1} + {(λ I - B)}^{- 1} [R_{d} x_{k - 1} - \sum_{l = 1}^{k - 1} ∆_{k - l} x_{l}], \\ k > 2 \end{array}

（26）

该公式计算结果与经典摄动法求解特征值扰动的结果一致^［

11］. 高阶摄动解可以通过递归的方式计算，在计算第

k

阶的结果时，需要采用前

k - 1

阶的摄动解. 当所需递归步数增加，数据存储量亦将增多. 由式（23）可见，与本文提出的迭代算法类似，摄动级数展开法的收敛性也受到初始矩阵特征值分布间距的影响，距其他特征值较远的特征值对应的扰动计算更容易收敛. 在迭代算法和摄动级数展开法皆收敛时，两者都能得出准确的结果. 迭代算法每步的修正量为三阶小量，而摄动级数展开法每步修正量为二阶小量，因而本文提出的迭代算法收敛速度更快.

3 矩阵降阶方法

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

T_{r} = V^{- 1} (Λ + R) V

（27）

其中：

V = \prod_{i = 1}^{n - 1} W_{i}, V^{- 1} = \prod_{i = 1}^{n - 1} W_{n - i}^{- 1}

（28）

W_{i} = [\begin{matrix} I_{i - 1} & 0 \\ 0 & Q_{i} \end{matrix}]

，

W_{i}^{- 1} = [\begin{matrix} I_{i - 1} & 0 \\ 0 & Q_{i}^{- 1} \end{matrix}]

（29）

Q_{i} = [\begin{matrix} 1 & 0 \\ x^{(i)} & I_{n - i} \end{matrix}]

，

Q_{i}^{- 1} = [\begin{matrix} 1 & 0 \\ - x^{(i)} & I_{n - i} \end{matrix}]

（30）

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

V (i + 1 : n, i) = x^{(i)}

（31）

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

图2 计算 $V$ 、 $V^{- 1}$ 和 $T_{r}$ 的流程图

Fig.2 The flowchart for computing $V$ 、 $V^{- 1}$ and $T_{r}$

上三角矩阵 $T_{r}$ 的特征值分解容易计算：

T_{r} = Ψ \tilde{Λ} Ψ^{- 1}

（32）

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

\tilde{P} = P V Ψ

（33）

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

A + A_{p} = (P V Ψ) \tilde{Λ} (Ψ^{- 1} V^{- 1} P^{- 1})

（34）

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

4 算例

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

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

Λ =

d i a g [0.1 + 1.0 i 1.1 + 0.1 i 1.0 + 1.0 i 1.05 + 1.05 i]

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

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

图3 特征值扰动前后位置变化

Fig.3 Distribution of the eigenvalues before and after perturbation

表1 各特征值扰动计算迭代算法的收敛性指标

Tab.1 Convergence index of iterative algorithm for perturbation of eigenvalues

特征值	最小距离	$ρ$
$λ_{1}$	1	0.42
$λ_{2}$	1	0.95
$λ_{3}$	0.07	60.55
$λ_{4}$	0.07	121.07

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

δ_{E} = \frac{|\tilde{λ} (k) - {\tilde{λ}}_{E}|}{|{\tilde{λ}}_{E}|} \times 100 %

（35）

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

表2 摄动法与迭代算法计算

λ_{1}

扰动结果对比

Tab.2 Comparision of Perturbation series expansion method with Iterative algorithm for

λ_{1}

步数	摄动级数展开法			迭代算法
步数	$\tilde{λ} (k)$	δ_E/%	$\|Δ_{k}\|$	$\tilde{λ} (k)$	δ_E/%	$\|Δ_{k}\|$
1	-0.090 0+1.200 0i	5.99	0.668 6	-0.165 4+1.235 1i	0.71	0.642 8
2	-0.151 8+1.253 0i	1.77	0.303 2	-0.155 4+1.230 8i	0.17	0.073 4
3	-0.145 5+1.236 1i	1.02	0.200 1	-0.157 7+1.232 0i	0.04	0.018 0
4	-0.154 6+1.232 6i	0.23	0.084 1	-0.157 2+1.237i	0.01	0.004 1

图4 $|Δ_{k}|$ 随迭代步数增加的变化

Fig.4 Trend of $|Δ_{k}|$ as iterative step number increase

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

Λ^{(3)} = d i a g [\begin{matrix} 1.0 + 1.0 i & 1.05 + 1.05 i \end{matrix}]

R^{(3)} =

[\begin{matrix} 0.132 9 + 0.016 7 i & 0.297 3 + 0.188 3 i \\ 0.079 2 - 0.122 0 i & - 0.393 6 - 0.089 2 i \end{matrix}]

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

图5 分步计算 $λ_{3}$ 特征值扰动时调用迭代算法 $|Δ_{k}|$ 的变化

Fig.5 Trend of $|Δ_{k}|$ for computing $λ_{3}$ using iterative algorithm

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

V =

[\begin{matrix} 1 & 0 & 0 & 0 \\ - 0.147 7 + 0.223 2 i & 1 & 0 & 0 \\ - 0.156 4 - 0.005 1 i & - 0.043 9 - 0.482 0 i & 1 & 0 \\ - 0.325 7 + 0.386 7 i & - 0.211 7 + 0.166 8 i & 0.135 6 - 0.222 1 i & 1 \end{matrix}]

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

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

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

图6 平面桁架结构布置（单位：m）

Fig.6 Layout of the plane truss structure（unit：m）

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

A_{c} = [\begin{matrix} 0 & I \\ - M^{- 1} K & - M^{- 1} C \end{matrix}]

（36）

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

λ_{c i}, λ_{c i}^{*} = - ξ_{i} ω_{i} \pm i ω_{i} \sqrt[]{1 - ξ_{i}^{2}}

（37）

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

表3 平面桁架前5阶模态频率及阻尼比

Tab.3 The first 5 modal frequency and the eigenvalues of the plane truss

阶数	$ω_{i}$ /（rad·s^-1）	$ξ_{i}$	$λ_{c i}, λ_{c i}^{*}$
1	80.173 6	0.010 2	-0.821 4 $\pm$ 80.169 4i
2	217.471 4	0.013 2	-2.864 7 $\pm$ 217.452 5i
3	312.366 7	0.017 2	-5.378 6 $\pm$ 312.320 4i
4	549.511 2	0.028 4	-15.598 1 $\pm$ 549.289 8i
5	697.006 4	0.035 6	-24.791 0 $\pm$ 696.565 4i

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

∆ E (t) = E \cdot 0.2 s i n (2 π t)

（38）

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

（a）第1阶

（b）第2阶

（c）第3阶

（d）第4阶

（e）第5阶

图7 桁架结构前5阶模态频率变化曲线

Fig.7 Modal frequency variation curves of the first 5^th orders for the truss

5 结论

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

参考文献

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． [百度学术]

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

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

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

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

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． [百度学术]

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． [百度学术]

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． [百度学术]

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

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

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） [百度学术]

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

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

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

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） [百度学术]

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

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） [百度学术]

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

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

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） [百度学术]

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． [百度学术]

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． [百度学术]

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． [百度学术]

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

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） [百度学术]

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

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） [百度学术]

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

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． [百度学术]

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

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

GÜ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． [百度学术]

作者稿件一经被我刊录用，如无特别声明，即视作同意授予我刊论文整体的全部复制传播的权利，包括但不限于复制权、发行权、信息网络传播权、广播权、表演权、翻译权、汇编权、改编权等著作使用权转让给我刊，我刊有权根据工作需要，允许合作的数据库、新媒体平台及其他数字平台进行数字传播和国际传播等。特此声明。

关闭

首页

期刊简介

编委会

作者中心

下载中心

学术道德

常见问题

版权声明

联系我们

English

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

摘要

关键词

1 特征值扰动问题

2 特征值扰动计算方法

2.1 迭代算法

2.2 摄动级数展开法

3 矩阵降阶方法

4 算例

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

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

5 结论

参考文献

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

摘要

关键词

1 特征值扰动问题

2 特征值扰动计算方法

2.1 迭代算法

2.2 摄动级数展开法

3 矩阵降阶方法

4 算例

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

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

5 结 论

参考文献

5 结论