第二十八篇 广义特征值问题
通常在工程实践中，在特征值方程的右边会有一个额外的矩阵，导致会编程这种形式
比如，在固有频率问题和屈曲分析中，[K]是系统的“刚度矩阵”，[M]是系统的“质量”或“几何”矩阵。
通过重新排列上面的方程，可以写出任意一个等价的特征值方程，
本程序对应最下面方程的最大特征值1/λ，其倒数为前一个方程的最小特征值λ。
对最开始方程进行向量迭代，让λ = 1，并猜测右边的{x}0。矩阵与向量的乘积会得到
通过求解线性方程组得到{x}∗1的新估计
当新的{x}∗1被计算出来时，它可以通过除以“最大”分量来达到正交化，从而得到{x}1，并从带回开始方程，重复这个过程直至收敛。由于在整个迭代过程中[K]矩阵不变，通过求得[K]的[L][U]因子，可以更加有效地进行迭代过程。应用之后，就是在每次迭代中从前和从后替换计算{x}∗i，详情可以参看之前的三篇文章，移位取逆迭代，移位向量迭代，向量迭代。
程序如下：
其中有一个主程序，四个子程序，分别为检查收敛的子程序checkit，因式分解的子程序lufac，从前迭代的子程序subfor，从后迭代的子程序subbac。详情可以参看LDLT分解高斯消元
主程序：
#Kx=λMx的向量迭代 import numpy as npimport Bn=4;tol=1.0e-5;limit=100lower=np.zeros((n,n))upper=np.zeros((n,n))k=np.array([[8,4,-24,0],[4,16,0,4],[-24,0,192,24],[0,4,24,8]],dtype=np.float)m=np.array([[0.06667,-0.01667,-0.1,0],[-0.01667,0.1333,0,-0.01667],[-0.1,0,4.8,0.1],[0,-0.01667,0.1,0.06667]],dtype=np.float)x1=np.zeros((n,1))x=np.ones((4,1),dtype=np.float)print('矩阵K')print(k[:])print('矩阵M')print(m[:])print('初始猜测值',x[:,0])B.lufac(k,lower,upper)print('前几次迭代值')iters=0while(True):iters=iters+1x1=np.dot(m,x)B.subfor(lower,x1)B.subbac(upper,x1)big=0.0for i in range(1,n+1):if abs(x1[i-1,0])>abs(big):big=x1[i-1,0]x1[:,0]=x1[:,0]/bigif B.checkit(x1,x,tol)==True or iters==limit:breakx[:,0]=x1[:,0]if iters<5:for i in range(1,n+1):print("{:13.4e}".format(x[i-1,0]),end=" ")print(end="\n")l2=np.linalg.norm(x1)x1[:,0]=x1[:,0]/l2print('迭代到收敛的次数',iters)print('M转秩*K的“最小”特征值',"{:13.4e}".format(1.0/big))print('对应的特征向量')for i in range(1,n+1):print("{:13.4e}".format(x1[i-1,0]),end=" ")
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
lufac
def lufac(a,lower,upper):n=a.shape[0]upper[0,:]=a[0,:]for i in range(1,n+1):lower[i-1,i-1]=1.0for k in range(1,n):if abs(upper[k-1,k-1])>1.0e-10:for i in range(k+1,n+1):#下三角分解for j in range(1,i):total=0for l in range(1,j):total=total-lower[i-1,l-1]*upper[l-1,j-1]lower[i-1,j-1]=(a[i-1,j-1]+total)/upper[j-1,j-1]#上三角分解for j in range(1,n+1):total=0for l in range(1,i):total=total-lower[i-1,l-1]*upper[l-1,j-1]upper[i-1,j-1]=a[i-1,j-1]+totalelse:print('有0向量在第',k,'行')break
subfor
def subfor(a,b):#一个下三角的从前迭代法n=a.shape[0]for i in range(1,n+1):total=b[i-1]if i>1:for j in range(1,i):total=total-a[i-1,j-1]*b[j-1]b[i-1]=total/a[i-1,i-1]
subbac
def subbac(a,b):#一个下三角的从后迭代法n=a.shape[0]for i in range(n,0,-1):total=b[i-1]if i<n:for j in range(i+1,n+1):total=total-a[i-1,j-1]*b[j-1]b[i-1]=total/a[i-1,i-1]
终端输出结果如下
