参考链接 https://wenku.baidu.com/view/ee7ecbeca98271fe910ef9fc.html?from=search
幂迭代算法：
逆幂迭代算法：
实际在使用中A可以先进行LU分解
无论是幂迭代，还是逆幂迭代，只能求出最大和最小特征值与对应的特征向量。
具体问题
python实现
1 import numpy as np
2 import matplotlib.pyplot as plt
3
4 eigenval1_ =
[]
5 eigenval2_ =
[]
6 eigenvec1_ =
[]
7 eigenvec2_ =
[]
8 cnt_ =
[]
9 x0 = np.array([1 / np.sqrt(2), 1 / np.sqrt(2
)])
10 x0.resize((2, 1
))
11 dis = 10 ** -12
12
13
14 def normalized_power_iteration(A):
15 v0 =
x0
16 distance = 1
17 times = 0.0
18
19 while distance >
dis:
20 v1 =
np.dot(A, v0)
21 v1 = v1 /
np.linalg.norm(v1)
22 distance = np.linalg.norm(v1 -
v0)
23 v0 =
v1
24 times += 1.0
25 eig = np.linalg.norm(np.dot(A, v1)) /
np.linalg.norm(v1)
26 return (v1, eig, times)
27
28
29 def normalized_inverse_power_iteration(A):
30 # 2x2矩阵 进行LU分解，在这里进行手动的LU分解求逆，因为题目找那个说了是一个2x2的矩阵，实际上当矩阵的维数比较高的时候，应该用LU分解然后解线性方程组来求解。
31 L = np.array([1, 0, A[1, 0] / A[0, 0], 1
])
32 L.resize((2, 2
))
33 LI = np.array([1, 0, -L[1, 0], 1
])
34 LI.resize((2, 2
))
35
36 U = np.array([A[0, 0], A[0, 1], 0, A[1, 1] - A[0, 1] * A[1, 0] /
A[0, 0]])
37 U.resize((2, 2
))
38 UI = np.array([1 / U[0, 0], -U[0, 1] / (U[0, 0] * U[1, 1]), 0, 1 / U[1, 1
]])
39 UI.resize((2, 2
))
40 lu =
np.dot(UI, LI)
41
42 v0 =
x0
43 distance = 1
44 times =
0
45
46 while distance >
dis:
47 u1 =
np.dot(lu, v0)
48 m1 =
np.linalg.norm(u1)
49 v1 = u1 /
m1
50 distance = np.linalg.norm(v1 -
v0)
51 v0 =
v1
52 times += 1
53 eig = np.linalg.norm(np.dot(lu, v1)) /
np.linalg.norm(v1)
54 return (v1, 1 /
eig, times)
55
56
57 for i
in np.arange(0, n):
58 (vec, val, round1) =
normalized_power_iteration(As[i, :, :])
59 eigenval1_.append(val)
60 eigenvec1_.append(vec[0, 0])
61 eigenvec1_.append(vec[1
, 0])
62 cnt_.append(round1)
63
64 (vec, val, round2) =
normalized_inverse_power_iteration(As[i, :, :])
65 eigenvec2_.append(vec[0, 0])
66 eigenvec2_.append(vec[1
, 0])
67 eigenval2_.append(val)
68
69 eigenval1 =
np.array(eigenval1_)
70 eigenvec1 =
np.array(eigenvec1_)
71 eigenvec1.resize((n, 2
))
72
73 cnt =
np.array(cnt_)
74
75 eigenval2 =
np.array(eigenval2_)
76 eigenvec2 =
np.array(eigenvec2_)
77 eigenvec2.resize((n, 2
))
78
79 plt.xlabel(
"bilv")
80 plt.ylabel(
"round times")
81 plt.title(
"result")
82 plt.show()
转载于:https://www.cnblogs.com/oldBook/p/9927168.html
