悬臂梁的博客

但行好事 莫问前程

0%

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

发表于 更新于 分类于 ，

前言

上一篇文章中介绍了计算特征值的逆迭代算法，我们知道基本逆迭代算法收敛于 \(\lambda_1\) 和 \(\Phi_1\)，基本正迭代算法收敛于 \(\lambda_n\) 和 \(\Phi_n\)。现在假定我们已经计算出了其中一个特征值对，如 \((\lambda_k, \Phi_k)\)，该怎么求解出其它的特征值对呢？

为了确保迭代过程不再收敛于 \((\lambda_k, \Phi_k)\)，需要收缩 (Deflation)矩阵或迭代向量。

收缩

矩阵收缩

矩阵收缩的基本思想是，求一个正交矩阵 \(\mathbf{H}\)，使得 \(\mathbf{HKH^{-1}}\) 为一个形如 \[ \left[ \begin{array}{c|cc} \lambda_1 & 0 & \cdots & 0 \\ \hline 0 & \\ \vdots & & \mathbf{K_1}\\ 0 \end{array} \tag{1} \right] \] 的矩阵。由于 \(\mathbf{K}\) 和 \(\mathbf{HKH^{-1}}\) 是相似的，故它们有相同的特征多项式。因此，若 \(\mathbf{HKH^{-1}}\) 形如 (1)，则 \[ \rm{det} (\mathbf{K - \lambda I}) = \rm{det} (\mathbf{HKH^{-1} -\lambda I}) = \mathbf{(\lambda_1 - \lambda)} \rm{det} (\mathbf{K_1 - \lambda I}) \] 由此可得，\(\mathbf{K}\) 的其余 \(n-1\) 个特征值是 \(\mathbf{K_1}\) 的特征值。

而一旦使用 \(\mathbf{K_1}\) 算出了第二个所求特征对后，就可以对 \(\mathbf{K_1}\) 重复该收缩过程，直到算出所有需要求解的特征值和特征向量。

应指出，矩阵 \(\mathbf{H}\) 不是唯一的，因此可以使用各种方法构建一个适当的变换矩阵。 另外，由于刚度矩阵 \(\mathbf{K}\) 是带状的，所以，该变换不应破坏其带状结构。

向量收缩与Gram-Schmidt方法

为了得到其它特征对，我们也可以选择不收缩矩阵，而是收缩迭代向量。

向量收缩的基本思想是，逆迭代或正迭代过程中，为了使迭代向量收敛于所求的特征向量，该迭代向量应不与特征向量正交。 反之，如果迭代向量正交于已计算出的特征向量，那么我们就消除了迭代向量收敛于这些特征向量的可能性， 即，迭代向量将收敛于其它特征向量。

Gram-Schmidt方法是一个被广泛应用的向量正交化方法。该方法可用于求解广义特征值问题 \(\mathbf{K \Phi = \lambda M \Phi}\)。

假设我们已经利用逆迭代算法求解出特征向量 \(\Phi_1,\Phi_2,\cdots,\Phi_m\)，我们需要 \(\mathbf{x_1}\) 与这些特征向量 \(\mathbf{M}\) 正交化。

在Gram-Schmidt正交化中，特征向量 \(\Phi_1,\Phi_2,\cdots,\Phi_m\) 的 \(\mathbf{M}\) 正交向量 \(\mathbf{\widetilde{x}_1}\) 的计算方法为： \[ \begin{align} \mathbf{\widetilde{x}_1} &= \mathbf{x}_1 - \sum_{i-1}^{m}{\alpha_i \Phi_i} \tag{2} \\ \alpha_i &= \mathbf{\Phi_i^T M x_1}; \quad i=1,\cdots,m \tag{3} \\ \end{align} \]

在逆迭代算法中，我们现在就可以使用 \(\mathbf{\widetilde{x}_1}\) 作为初始迭代向量。然后，只要 \(\mathbf{x_1^T M \Phi_{m+1}} \neq 0\)，迭代向量就将收敛于 \(\Phi_{m+1}\) 和 \({\lambda_{m+1}}\)。

代码实现和测试

在上一篇文章中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
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
import scipy as sp
import numpy as np


def inverse_iteration_method1(k: np.matrix, m: np.matrix):
    n = k.shape[0]
    eig_val = np.zeros(n)
    eig_vec = np.matrix(np.zeros((n, n)))

    eps = 1e-6
    for i in range(n):
        x = np.matrix(np.ones(n).reshape(n, 1))
        x_pre = x
        while True:
            for j in range(i):
                alpha = eig_vec[:, j].transpose() * m * x
                x -= alpha[0, 0] * eig_vec[:, j]
            x = np.linalg.solve(k, m * x)
            norm = x.transpose() * m * x
            x /= np.sqrt(norm)
            if np.all(abs(x - x_pre) < eps):
                eig_val[i] = 1/np.sqrt(norm[0, 0])
                eig_vec[:, i] = x
                break
            x_pre = x
    return eig_val, eig_vec


def inverse_iteration_method2(k: np.matrix, m: np.matrix):
    n = k.shape[0]
    eig_val = np.zeros(n)
    eig_vec = np.matrix(np.zeros((n, n)))

    eps = 1e-6
    rho_pre = 0.0
    for i in range(n):
        x = np.matrix(np.ones(n).reshape(n, 1))
        y = m * x
        while True:
            x = np.linalg.solve(k, y)
            for j in range(i):
                alpha = eig_vec[:, j].transpose() * m * x
                x -= alpha[0, 0] * eig_vec[:, j]
            y_bar = m * x
            rho = x.transpose() * y / (x.transpose() * y_bar)
            if abs(rho - rho_pre) / rho < eps:
                eig_val[i] = rho[0, 0]
                eig_vec[:, i] = x / np.sqrt(x.transpose() * y_bar)
                break
            y = y_bar / np.sqrt(x.transpose() * y_bar)
            rho_pre = rho
    return eig_val, eig_vec


def test():
    k = np.matrix(np.array([
        [5, -4, 1, 0],
        [-4, 6, -4, 1],
        [1, -4, 6, -4],
        [0, 1, -4, 5]
    ]))
    m = np.matrix(np.array([
        [2, 0, 0, 0],
        [0, 2, 0, 0],
        [0, 0, 1, 0],
        [0, 0, 0, 1]
    ]))
    n = k.shape[0]
    eig_vals, eig_vecs = sp.linalg.eig(k, m)[:2]
    eig_vals = eig_vals.real
    argsort = np.argsort(eig_vals)
    eig_vals = [eig_vals[argsort[i]] for i in range(n)]
    eig_vecs_copy = eig_vecs.copy()
    for i in range(n):
        eig_vecs[:, i] = eig_vecs_copy[:, argsort[i]]
    print('SciPy库：')
    print(eig_vals)
    print(eig_vecs)

    eig_vals, eig_vecs = inverse_iteration_method1(k, m)
    print('方法1：')
    print(f'特征值：\n{eig_vals}')
    print(f'特征向量：\n{eig_vecs}')

    eig_vals, eig_vecs = inverse_iteration_method2(k, m)
    print('方法2：')
    print(f'特征值：\n{eig_vals}')
    print(f'特征向量：\n{eig_vecs}')

程序输出如下： 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
SciPy库：
[0.09653732854936428, 1.3914654511583402, 4.373549554582955, 10.638447665709336]
[[-0.38576651  0.50215561 -0.55097093  0.11195372]
 [-0.61138815  0.14033458  0.52342732 -0.26606756]
 [-0.59120182 -0.55197537  0.02916661  0.75798869]
 [-0.35758794 -0.65074507 -0.64931054 -0.58491672]]
方法1：
特征值：
[ 0.09653733  1.39146545  4.37354955 10.63844769]
特征向量：
[[ 0.31262955  0.44526604  0.43866944 -0.10756245]
 [ 0.49547586  0.12443555 -0.41674032  0.25563076]
 [ 0.4791166  -0.48944215 -0.02322167 -0.72825456]
 [ 0.28979328 -0.5770219   0.51696616  0.56197132]]
方法2：
特征值：
[ 0.09653733  1.39146547  4.37354974 10.63844707]
特征向量：
[[ 0.31263514  0.44529746  0.43861285 -0.10764853]
 [ 0.49547613  0.12439816 -0.41670087  0.25571251]
 [ 0.47911154 -0.48944852 -0.023328   -0.72825001]
 [ 0.28978869 -0.57698412  0.51712097  0.56186986]]

可以看出，程序正确计算出了矩阵的全部特征值和特征向量。（特征向量因归一化方式不同，存在数字上的差异。若以相同的规则缩放，结果将相同）

参考文献

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

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

[3]. 朱伯芳.有限单元法原理与应用-第3版[M].中国水利水电出版社,2009.