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

对于标准特征值方程

由特征值问题编程基础可知，对于任何非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+1big=0for 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=jif abs(big)<1.0e-20:breakden=a[nr-1,nr-1]-a[nc-1,nc-1]if abs(den)<1.0e-20:alpha=pi/4.0if hold<0:alpha=-alphaelse:alpha=math.atan(2.0*hold/den)/2.0ct=math.cos(alpha);st=math.sin(alpha);p[:]=0for i in range(1,n+1):p[i-1,i-1]=1.0p[nr-1,nr-1]=ct;p[nc-1,nc-1]=ct;p[nr-1,nc-1]=-st;p[nc-1,nr-1]=sta=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:breakeold[:,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()

checkit

def checkit(loads,oldlds,tol):
#检查多个未知数的收敛neq=loads.shape[0]big=0.0converged=Truefor 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=Falsecheckit=convergedreturn  checkit

eliminate

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

终端输出结果如下：

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

雅可比旋转（Jacobi法）求对称矩阵的特征值和特征向量
雅可比方法该方法是求解对称矩阵全部特征值和特征向量的一种方法,它基于以下结论: ①任何实对称矩阵A可以通过正交相似变换成对角型,即存在正交矩阵Q,使得 QTAQ=diag(λ1,λ2,-,λn)Q^ ...
【Python】Numpy--np.linalg.eig()求对称矩阵的特征值和特征向量
[Python]Numpy–np.linalg.eig()求对称矩阵的特征值和特征向量文章目录 [Python]Numpy--np.linalg.eig()求对称矩阵的特征值和特征向量 1. 介绍 ...
c语言求矩阵特征值的程序,如何用C语言编写求对称矩阵的特征值和特征向量的程序编写对称矩阵的特征值和特征向量,其中矩阵用二维数组保存.特征向量要求有大到小放到数组里....
优质解答 //数值计算程序-特征值和特征向量 // //约化对称矩阵为三对角对称矩阵 //利用Householder变换将n阶实对称矩阵约化为对称三对角矩阵 //a-长度为n*n的数组,存放n阶实对称 ...
雅可比旋转求解对称二维矩阵的特征值和特征向量
问题描述: 给定一个矩阵,如下: A=[a11a21a12a22] A=\begin{bmatrix} a_{11}&a_{12}\\ a_{21}& a_{22} \end{bmat ...
对称矩阵到三对角矩阵的Lanczos推导(python,数值积分)
第三十二篇 Lanczos转化到三对角形式在之前的篇章里,有许多求解线性方程的迭代方法,如最陡下降法,可以通过向量乘法和各种简单的向量运算,简化为一个单个矩阵的循环.将矩阵化为三对角形式的Lancz ...
雅可比（Jacobi）迭代法解线性方程组的Matlab实现
雅可比(Jacobi)迭代法解线性方程组的Matlab实现代码运行代码迭代法解线性方程组的基本思想是构造一串收敛到解的序列,即建立一种从已有近似解计算新的近似解的规则,有不同的计算规则得到不同 ...
python最小化打开exe_如何用python使GoAgent窗口打开后自动最小化以及关闭之前的py.exe窗口...
python:3.4 goagent:3.1.22-33 写了一个想在ipv4/6之间切换的小脚本 path=r'D:\Documents\Downloads\Downloads\goagent-go ...
QR分解求矩阵全部特征值
QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式将A=A1化成相似的上三角阵,从而求出矩阵A的全部特征值. QR方法的计算步骤如下: 下面就依次进行介绍. 一. 将一般矩阵化为上H ...
numpy求矩阵的特征值与特征向量(np.linalg.eig函数详解)
numpy求矩阵的特征值与特征向量(np.linalg.eig) 语法 np.linalg.eig(a) 功能 Compute the eigenvalues and right eigenvecto ...

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

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

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

最新文章

热门文章

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

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

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

最新文章

热门文章

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