悬臂梁的博客

但行好事 莫问前程

0%

特征问题的解法——向量迭代法(1)

发表于 更新于 分类于 ，

前言

在工程中，我们经常需要对模型进行模态分析，以确定模型的模态频率和振型。

模态分析本质上是求解特征值问题 \[ {\mathbf K \phi} = \lambda {\mathbf M \phi} \tag{1} \] 特别是求解最小特征值 \(\lambda_1, \dots, \lambda_p\) 及其特征向量 \(\phi_1, \dots, \phi_p\) 。

常用的特征值求解算法可分为四类，对应于求解算法利用的基本性质。

第一类是向量迭代法(Vector Iteration Method)，又称为幂法(Power Iteration Method)，利用的基本性质是 \[ {\mathbf K \phi_i} = \lambda_i{\mathbf M \phi_i} \tag{2} \]

第二类为变换法，利用的基本性质是 \[ {\mathbf{\Phi^T K \Phi}} = {\mathbf{\Lambda}} \tag{3} \] \[ {\mathbf{\Phi^T M \Phi}} = {\mathbf{I}} \tag{4} \] 其中， \({\rm \Phi} = \left[ {\rm \phi_i, \dots, \phi_n} \right]\) 和 \({\rm \Lambda} = {\rm diag(\lambda_i)}, i=1,\dots,n\) 。

第三类为多项式迭代法，利用的基本性质是 \[ p(\lambda_i) = 0 \tag{5} \] 其中， \[ p(\lambda) = {\rm{det}(\mathbf{K} - \lambda \mathbf{M})} \tag{6} \]

第四类求解算法利用特征多项式的 Sturm 序列性质 \[ p(\lambda) = {\rm{det}(\mathbf{K} - \lambda \mathbf{M})} \tag{7} \] \[ p^{(r)}(\lambda^{(r)}) = {\rm{det}(\mathbf{K}^{(r)} - \lambda^{(r)}\mathbf{M}^{(r)})}; \quad r=1,\dots,n-1 \tag{8} \] 其中， \(p^{(r)}(\lambda^{(r)})\) 是对应于\({\rm K \phi = \lambda M \phi}\) 的第 \(r\) 个相伴约束问题的特征多项式。

在四类求解方法中，每一类都提出了许多算法。

在学习具体算法前，需要认识到，求解算法应具有迭代性质， 因为求解特征问题 \({\rm K \phi = \lambda M \phi}\) 就等价于计算多项式 \(p(\lambda)\) 的根，它的阶等于 \({\rm K}\) 和 \({\rm M}\) 的阶。

接下来开始介绍第一类方法：向量迭代法。而在介绍向量迭代法之前，先尝试理解其基本思想。

向量迭代法

基本思想

假设向量 \({\rm A}\) 有 \(n\) 个线性无关的特征向量 \(\mathit{x_1, \dots, x_n}\) ，且相应的特征值满足 \[ |\lambda_1| > |\lambda_2| \geq \dots \geq |\lambda_n| \]

给定一个任意向量 \(\mathit{v}_0\) ，假设 \[ \mathit{v}_0 = \alpha_1 \mathit{x}_1 + \dots + \alpha_n \mathit{x}_n \]

将 \({\rm A}\) 作用于该向量，有 \[ \mathbf{A} \mathit{v}_0 = \alpha_1 \lambda_1 \mathit{x}_1 + \alpha_2 \lambda_2 \mathit{x}_2 + \dots + \alpha_n \lambda_n \mathit{x}_n \\ \mathbf{A}^2 \mathit{v}_0 = \alpha_1 \lambda_1^2 \mathit{x}_1 + \alpha_2 \lambda_2^2 \mathit{x}_2 + \dots + \alpha_n \lambda_n^2 \mathit{x}_n \]

且一般地， \[ \mathbf{A}^k \mathit{v}_0 = \alpha_1 \lambda_1^k \mathit{x}_1 + \alpha_2 \lambda_2^k \mathit{x}_2 + \dots + \alpha_n \lambda_n^k \mathit{x}_n \tag{9} \]

若定义 \[ \mathit{v}_k = \mathbf{A}^k \mathit{v}_0, \quad k=1, 2, \dots \]

则 \[ \frac{1}{\lambda_1^k} \mathit{v}_k = \alpha_1 \mathit{x}_1 + \alpha_2 (\frac{\lambda_2}{\lambda_1})^k \mathit{x}_2 + \dots + \alpha_2 (\frac{\lambda_n}{\lambda_1})^k \mathit{x}_n \tag{10} \]

由于 \[ \left| \frac{\lambda_i}{\lambda_1} \right| < 1, \quad, i=2,3,\dots,n \]

由此得到： \[ \frac{1}{\lambda_1^k} \mathit{v}_k \rightarrow \alpha_1 \mathit{x}_1, \quad k \rightarrow \infty \]

因此，若 \(\alpha_1 \neq 0\) ，则序列 \(\{ (1/\lambda_1^k) \mathit{v}_k \}\) 收敛到 \(\mathbf{A}\) 的特征向量 \(\alpha_1 \mathit{x}_1\)。

当然，由于 \(\lambda_1\) 是未知的，所以无法计算 \(\{ (1/\lambda_1^k) \mathit{v}_k \}\) 。 但好在不需要将序列 \(\{ \mathit{v}_k \}\) 用 \(1/\lambda_1^k\) 进行缩放。

这就是正迭代法的思想，能计算出最大的特征值。

而模态分析需要计算的是最小特征值，需要使用的是逆迭代法，即需要将 \(\mathbf{A}^{-1}\) 作用在向量上。

逆迭代和正迭代

向量迭代法所考虑的基本关系式是 \[ {\mathbf K \phi} = \lambda {\mathbf M \phi} \tag{1} \]

选择 \(\phi\) 的一个向量 \(\mathbf{x}_1\) ，对 \(\lambda\) 设定一个值，令 \(\lambda = 1\) 。于是，可以计算公式 (1) 的右手边，即可以计算 \[ \mathbf{R}_1 = (1)\mathbf{M x_1} \tag{11} \] 由于 \(\mathbf{x}_1\) 是任意向量，一般不满足 \(\mathbf{K x_1 = R_1}\) 。考虑静平衡方程 $$

; $$ 其中， \(\mathbf{x_2}\) 是位移解，对应于作用力 \(\mathbf{R}\) 。 由于我们是使用迭代法求解特征向量，所以我们可以直观地认为，\(\mathbf{x_2}\) 是比 \(\mathbf{x_1}\) 更好的近似特征向量。 通过反复进行此迭代，可以得到越来越好的近似特征向量。

上述过程即为逆迭代的基础。

而在正迭代中，迭代过程是反向的，即在第一步中，计算 \(\mathbf{R_1 = K x_1}\) ， 然后通过求解 \(\mathbf{M x_2 = R_1}\) ，得到改进的近似特征向量 \(\mathbf{x_2}\) 。

基本算法描述

首先考虑逆迭代法中所用的基本方程。

在求解中，假设初始迭代向量为 \(\mathbf{x}_1\)，在每个迭代步骤 \(k = 1, 2,\dots\) 时，计算 \[ \begin{align} \mathbf{K \bar{x}}_{k+1} &= \mathbf{M x}_k \tag{13} \\ \mathbf{x}_{k+1} &= \frac{\mathbf{\bar{x}}_{k+1}}{(\mathbf{\bar{x}}_{k+1}^T \mathbf{M} \mathbf{\bar{x}}_{k+1})^{1/2}} \tag{14} \end{align} \] 其中，只要 \(\mathbf{x}_1\) 和 \(\phi_i\) 不与 \(\mathbf{M}\) 正交，即 \(\mathbf{x_1^T M \phi_1} \neq 0\) ，就有 \[ 当 k \rightarrow \infty, \quad \mathbf{x}_{k+1} \rightarrow \phi_1 \]

迭代法的基本步骤是求解方程 (13)，每求解一次，我们就能够得到比前一次迭代向量 \(\mathbf{x}_{k}\) 更接近特征向量的新向量 \(\mathbf{x}_{k+1}\) 。 式 (14) 的计算对新的迭代向量进行缩放，使其与 \(\mathbf{M}\) 的加权长度为 \(1\) ，即令 \(\mathbf{x}_{k+1}\) 满足质量正交条件 \[ \mathbf{x_{k+1}^{T} M x_{k+1}} = 1 \tag{15} \] 当然，也可以将向量 \(\mathbf{x}_{k+1}\) 缩放到其他长度。

更有效的算法

式 (13)、(14) 是基本的逆迭代算法，但在实际实现中，以下方法更加有效。

假设 \(\mathbf{y}_1 = \mathbf{M x_1}\) ，计算 \(k = 1, 2, \dots\) ， \[ \begin{align} \mathbf{K \bar{x}}_{k+1} &= \mathbf{y}_k \tag{16} \\ \mathbf{\bar{y}}_{k+1} &= \mathbf{M \bar{x}}_{k+1} \tag{17} \\ \rho(\mathbf{\bar{x}}_{k+1}) &= \frac{\mathbf{\bar{x}}_{k+1}^T \mathbf{y}_k}{\mathbf{\bar{x}}_{k+1}^T \mathbf{\bar{y}}_{k+1}} \tag{18} \\ \mathbf{y}_{k+1} &= \frac{\mathbf{\bar{y}}_{k+1}}{(\mathbf{\bar{x}_{k+1}^T \mathbf{\bar{y}}_{k+1}})^{1/2}} \tag{19} \end{align} \] 其中，只要 \(\mathbf{y}_1^T \phi_1 \neq 0\) ，有 \[ 当 k \rightarrow \infty 时，\mathbf{y}_{k+1} \rightarrow \mathbf{M \phi_1} 和 \phi(\mathbf{\bar{x}}_{k+1}) \rightarrow \lambda_1 \] 本算法对 \(\mathbf{y}_k\) 进行迭代而非对 \(\mathbf{x}_k\) 进行迭代。 式 (18) 中，我们得到了特征值 \(\lambda_1\) 的近似值，其以瑞利商 \(\rho(\mathbf{\bar{x}}_{k+1})\) 的形式给出。 \(\lambda_1\) 的近似值可以用于评估迭代的收敛情况。 若当前的 \(\lambda_1\) 的近似值用 \(\lambda_1^{(k+1)}\) 表示，即 \(\lambda_1^{(k+1)} = \rho(\mathbf{\bar{x}}_{k+1})\)， 那么我们通过式 (20) 判断迭代是否收敛： \[ \frac{\left| \lambda_1^{(k+1)} - \lambda_1^{(k)} \right|}{\lambda_1^{(k+1)}} \leq tol \tag{20} \] 当特征值 \(\lambda_1\) 需要 \(2s\) 的精度时，\(tol\) 应该是 \(10^{-2s}\) 或更小。

代码实现和测试

实现

使用 Python 语言对基本算法和更高效的算法进行实现，代码如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
import numpy as np


def inverse_iteration_method1(k: np.ndarray, m: np.ndarray):
    n = k.shape[0]
    x = np.ones((n, 1))
    x_pre = x
    i_iter = 0
    while True:
        i_iter += 1
        x = np.linalg.solve(k, m @ x)
        norm = x.transpose() @ m @ x
        x /= np.sqrt(norm)
        if np.all(abs(x - x_pre) < 1e-6):
            break
        x_pre = x
    print(f'方法1：迭代{i_iter}次')
    print(f'特征值：\n{1/np.sqrt(norm)}')
    print(f'特征向量：\n{x}')


def inverse_iteration_method2(k: np.ndarray, m: np.ndarray):
    n = k.shape[0]
    x = np.ones((n, 1))
    y = m @ x
    rho_pre = 0.0
    n_iter = 0
    while True:
        n_iter += 1
        x = np.linalg.solve(k, y)
        y_bar = m @ x
        rho = x.transpose() @ y / (x.transpose() @ y_bar)
        if abs(rho - rho_pre) / rho < 1e-6:
            break
        y = y_bar / np.sqrt(x.transpose() @ y_bar)
        rho_pre = rho
    print(f'方法2：迭代{n_iter}次')
    print(f'特征值：\n{rho}')
    print(f'特征向量：\n{x / np.sqrt(x.transpose() @ y_bar)}')

测试

使用文献 [1] 中的例子进行测试。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
def test1():
    k = np.array([
        [2, -1, 0, 0],
        [-1, 2, -1, 0],
        [0, -1, 2, -1],
        [0, 0, -1, 1]
    ])
    m = np.array([
        [0, 0, 0, 0],
        [0, 2, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 1]
    ])
    inverse_iteration_method1(k, m)


def test2():
    k = np.array([
        [2, -1, 0, 0],
        [-1, 2, -1, 0],
        [0, -1, 2, -1],
        [0, 0, -1, 1]
    ])
    m = np.array([
        [0, 0, 0, 0],
        [0, 2, 0, 0],
        [0, 0, 0, 0],
        [0, 0, 0, 1]
    ])
    inverse_iteration_method2(k, m)

程序输出如下。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
方法1：迭代8次
特征值：
[[0.14644661]]
特征向量：
[[0.25000003]
 [0.50000006]
 [0.60355338]
 [0.70710669]]
方法2：迭代5次
特征值：
[[0.14644661]]
特征向量：
[[0.25000638]
 [0.50001275]
 [0.60355075]
 [0.70708874]]

最小特征值和特征向量的精确解为 \(\lambda_1 = \frac{1}{2} - \frac{\sqrt{2}}{4}\) 和 \(\phi_1 = [\frac{1}{4} \quad \frac{1}{2} \quad \frac{1+\sqrt{2}}{4} \quad \frac{\sqrt{2}}{2}]^T\)，可以看出，两种逆迭代算法正确算出了最小特征值和特征向量，而后一种算法所需的迭代次数更少。

计算出最小特征值后，剩余特征值的计算方法将在下一篇文章中介绍。

参考文献

[1]. Bathe K J. Finite element procedures[M]. Klaus-Jurgen Bathe, 2006.

[2]. 史蒂文・J.利昂著STEVENJ.LEON.线性代数[M].机械工业出版社,2015.