固有频率和屈曲分析，Kx=λM特征值和特征向量求解(python,数值积分)

第二十八篇 广义特征值问题

通常在工程实践中，在特征值方程的右边会有一个额外的矩阵，导致会编程这种形式

比如，在固有频率问题和屈曲分析中，[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 np
import B
n=4;tol=1.0e-5;limit=100
lower=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