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)