雅可比主对角线(Jacobi diagonalization)化求对称矩阵的特征值(python,数值积分)

第三十篇 雅可比主对角线化求对称矩阵的特征值

对于标准特征值方程

由特征值问题编程基础可知，对于任何非0解矩阵[P]，标准方程可以转化为具有相同特征值的方程

其中

这种转换技术的关键核心在于[A *]的特征值比原始[A]的特征值更容易找到。
然而，如果[A]是对称的，则变换后的矩阵不太可能保持对称。很容易证明下面变换

将保持[A∗]的对称性。为了使上面方程中给出的特征值[A∗]与[A]的特征值相同，必须把这两个性质综合起来

这种类型的矩阵称为“正交矩阵”，具有这种性质的矩阵称为“旋转矩阵”。

将此变换应用于下面的矩阵，我们有


其中很明显[A∗]对于任何α值都是对称的。这种情况，明显可以选择一个α值使得[A∗]成为一个对角矩阵，因为如果这样，对角线就是特征值。下面的情况，非对角线项将被消除

得出，tan α = 1和α = π/4，给出sin α = cos α = 1/√2。
得到的变换矩阵是

即[A]的特征值分别为3和1。
对于大于2 × 2的矩阵[A]，变换矩阵[P]必须通过在其他主对角线上放1，在所有非对角线上放0来“填充”，被消去的行和列选择上面的矩阵。例如，如果[A]是4 × 4，则变换矩阵可以选择6种形式中的一种，这取决于要消去初始矩阵中的哪些非对角项，例如:

上面第一个矩阵经过[P]T [A][P]变换后将原矩阵[A]中的a12和a21项消去，而第二个矩阵将消去a24和a42项。1和0的作用是让[A]的其他行和列保持不变。这意味着在一次转换中变为零的非对角线项在随后的转换中会恢复为非零值(尽管通常是很“小”的值)，因此正如期望的那样，该方法是迭代的。
这种类型的迭代的最早形式称为“雅可比对角化”，它通过消除每一次迭代剩余的“最大的”非对角项连续进行迭代。
对于任何对称矩阵[A]，得到广义方程为

得到[A∗]形式的非对角线项为

求α使这一项等于零

因此

因此，为了建立一个雅可比对角化的简单程序，必须在[a]中搜索“最大的”非对角线项，并找到它所在的行和列。“旋转矩阵”α按照之前的方法构建。矩阵[P]可以使用一个numpy库的transpose转化，然后矩阵乘积形成方程的[A *]。重复这个过程，直到[A∗]的主对角线在可接受的公差内收敛到[A]的特征值为止。
计算实例：
使用雅可比对角化去估算下面对称矩阵的特征值

下面的结果保留到小数点后四位，但实际计算的精确度更高。
第一次迭代
最大的非主对角项为a23 = a32 = −9.0,因此根据之前方程

第一次转换矩阵将包含下面的项

因此

通过转化矩阵得到

详细的数值为

最后

第二次迭代
最大的非主对角项为a12 = a21 = −7.7782，因此

第二次转化矩阵将包括下面的项

因此，

同上面一样，矩阵乘积将等于

可以看到，虽然位置(2,3)和(3,2)不再为零，但与初始矩阵中的值相比足够“小”。随着迭代的进行，旋转角度αk→0，变换矩阵[Pk]→[I]和变换矩阵[Ak]趋向于一个对角矩阵，特征值在对角线上。
对于这个例子，经过六次迭代，容差为1.0e-5
得到

因此[A]的特征值λ = 0.4659, 20.9681，−0.9340。特征向量将通过将每个特征值代入求线性方程的解。
程序如下
分为一个主程序和两个子程序，分别为判断收敛的子程序checkit，和高斯消元求特征向量的子程序eliminate。详情可参照之前文章的线性方程求解部分
主程序：

#雅可比主对角线化求对称矩阵的特征值
import numpy as np
import B
import math
n=3;tol=1.0e-5;limit=100
enew=np.zeros((n,1))
eold=np.zeros((n,1))
p=np.zeros((n,n))
a1=np.zeros((n,n))
a=np.array([[10,5,6],[5,20,4],[6,4,30]],dtype=np.float)
a2=a
pi=math.acos(-1)
x=np.zeros((n,1))
x=np.ones((3,1),dtype=np.float)
print('雅可比主对角线化求对称矩阵的特征值')
print('矩阵A')
print(a[:])
print('前几次迭代值')
iters=0;eold[:]=0
while(True):
    iters=iters+1
    big=0
    for i in range(1,n+1):
        for j in range(i+1,n+1):
            if abs(a[i-1,j-1]>big):
                big=abs(a[i-1,j-1]);hold=a[i-1,j-1];nr=i;nc=j
    if abs(big)<1.0e-20:
        break
    den=a[nr-1,nr-1]-a[nc-1,nc-1]
    if abs(den)<1.0e-20:
        alpha=pi/4.0
        if hold<0:
            alpha=-alpha
    else:
        alpha=math.atan(2.0*hold/den)/2.0
    ct=math.cos(alpha);st=math.sin(alpha);p[:]=0
    for i in range(1,n+1):
        p[i-1,i-1]=1.0
    p[nr-1,nr-1]=ct;p[nc-1,nc-1]=ct;p[nr-1,nc-1]=-st;p[nc-1,nr-1]=st
    a=np.dot(np.dot(np.transpose(p),a),p)
    if iters<5:
        for i in range(1,n+1):
            for j in range(1,n+1):
                print('{:13.4e}'.format(a[i-1,j-1]),end='')
            print(end='\n')
        print(end='\n')
    for i in range(1,n+1):
        enew[i-1,0]=a[i-1,i-1]
    if B.checkit(enew,eold,tol) or iters==limit:
        break
    eold[:,0]=enew[:,0]
print('迭代到收敛次数',iters)
print('最后的转化矩阵A')
for i in range(1,n+1):
    for j in range(1,n+1):
        print('{:13.4e}'.format(a[i-1,j-1]),end='')
    print(end='\n')
for i in range(1,n+1):
    a1[:]=a2[:]
    for j in range(1,n+1):
        a1[j-1,j-1]=a1[j-1,j-1]-a[i-1,i-1]
    x[:]=0;a1[i-1,i-1]=1.0e20;x[i-1]=1.0e20;x[:]=B.eliminate(a1,x)
    l2=np.linalg.norm(x)
    print('特征值','{:13.4e}'.format(a[i-1,i-1]))
    print('特征向量')
    for i in range(1,n+1):
        print('{:13.4e}'.format(x[i-1,0]/l2),end=' ')
    print()
    
    

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

checkit
1

def checkit(loads,oldlds,tol):
#检查多个未知数的收敛
  neq=loads.shape[0]
  big=0.0
  converged=True
  for i in range(1,neq+1):
    if abs(loads[i-1,0])>big:
      big=abs(loads[i-1,0])
  for i in range(1,neq+1):
    if abs(loads[i-1,0]-oldlds[i-1,0])/big>tol:
      converged=False
  checkit=converged
  return  checkit
1
2
3
4
5
6
7
8
9
10
11
12
13

eliminate
1

def eliminate(a,b):
  n=a.shape[0]
##确定主对角线最大值
  for i in range(1,n):
    big=abs(a[i-1,i-1]);ihold=i
    for j in range(i+1,n+1):
      if abs(a[j-1,i-1])>big:
        big=abs(a[j-1,i-1]); ihold=j
    if ihold!=i:
      for j in range(i,n+1):
        hold=a[i-1,j-1]; a[i-1,j-1]=a[ihold-1,j-1]; a[ihold-1,j-1]=hold
      hold=b[i-1,0]; b[i-1,0]=b[ihold-1,0]; b[ihold-1,0]=hold
##消元阶段
    for j in range(i+1,n+1):
      fac=a[j-1,i-1]/a[i-1,i-1]
      for l in range(i,n+1):
        a[j-1,l-1]=a[j-1,l-1]-a[i-1,l-1]*fac
      b[j-1]=b[j-1]-b[i-1]*fac
##从后迭代
  for i in range(n,0,-1):
    hold=0.0
    for l in range(i+1,n+1):
      hold=hold+a[i-1,l-1]*b[l-1]
    b[i-1]=(b[i-1]-hold)/a[i-1,i-1]
  return b
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

终端输出结果如下：