Python 振动分析 迭代法计算高阶特征值及特征向量

参考书 :

<<振动分析>> 张准 汪凤泉 编著 东南大学出版社 ISBN 7-80123-583-4

参考章节 : 4.6.2 和 4.6.3

<<数值分析>> 崔瑞彩 谢伟松 天津大学出版社 ISBN 7-5618-1366-X

参考章节 : 3.1

参考资料:

<<交替使用幂法和降阶法求解矩阵全部特征值>> 下载地址:https://pan.baidu.com/s/1fmNMnS8zyaMv4B_6jd7rnQ

2018-03-11 笔记

1 求解器代码

import numpy as np

class EigenValueModule(object) :
    K=np.mat([[1]])     #K矩阵
    M=np.mat([[1]])     #M矩阵
    s=1                 #前s阶
    Lambdas=[]
    Xs=[]
    def GetMax(self,Mat):
        position=np.argmax(Mat)
        row,column=Mat.shape
        row,column=divmod(position, column)
        Result=Mat[row,column]
        return Result
    def IsValid(self,LambdaLast,LambdaNext):
        #用于判断迭代出的特征值是否满足要求
        #要求 Abs(LambdaLast-LambdaNext)/LambdaLast < Epsilon
        myResult=False
        Epsilon=1e-3                    #0.001
        Ratio=abs(LambdaLast-LambdaNext)/LambdaLast
        if  Ratio<Epsilon :
            myResult=True
        return myResult
    def ReduceOrder(self,A,X):
        a=A[0]
        X=X/X[0,0]
        aX=X*a
        myResult=A-aX
        myResult=np.delete(myResult,0,0)
        myResult=np.delete(myResult,0,1)
        return myResult
    def solve(self):
        self.Lambdas=[]
        self.Xs=[]
        #导入外部变量
        K=self.K
        M=self.M
        s=self.s
        #初始化
        R=K.I
        A=R*M
        #进行计算
        for i in range(s):
            X=self.ChildSolve(A)
            A=self.ReduceOrder(A,X)
    def solveA(self,A):
        self.Lambdas=[]
        self.Xs=[]
        s=self.s
        #进行计算
        for i in range(s):
            X=self.ChildSolve(A)
            print(A)
            A=self.ReduceOrder(A,X)

    def ChildSolve(self,A):                                     #单个计算过程
        n=len(A)
        X=np.mat(np.ones((1,n))).T                              #创建试算向量
        XMax=self.GetMax(X)
        LambdaLast=XMax                                         #第一个Lambda
        Y=X/XMax
        X=A*Y
        XMax=self.GetMax(X)             
        LambdaNext=XMax                                         #下一个Lambda
        while self.IsValid(LambdaLast,LambdaNext)!=True :       #Lambda不符合要求
            LambdaLast=LambdaNext                               #更新Lambda
            Y=X/XMax
            X=A*Y
            XMax=self.GetMax(X)
            LambdaNext=XMax                                     #更新Lambda
        self.Lambdas.append(LambdaNext)
        self.Xs.append(X)
        return X

 

2 调试代码

　　调试代码1

 

import numpy as np
from SolveEigenValue import EigenValueModule
Figure=EigenValueModule()
Figure.M=np.mat([[1,0,0],[0,1,0],[0,0,1]])
Figure.K=np.mat([[2,-1,0],[-1,2,-1],[0,-1,1]])
Figure.s=3
Figure.solve()
print(Figure.Lambdas)

 

　　 调试代码2

 

import numpy as np
from SolveEigenValue import EigenValueModule
Figure=EigenValueModule()
Figure.s=3
A=np.mat([[1,7,8,9,3],[0,2,5,8,6],[0,2,4,2,2],[0,0,1,6,5],[0,0,0,0,9]])
Figure.solveA(A)
print(Figure.Lambdas)