【前言】子空間迭代法可同時(shí)求解幾個(gè)極端特征值和相應(yīng)的特征向量，但它有收斂較慢，運(yùn)算量較大，積累誤差的缺點(diǎn)；隨后，人們對(duì)其作了進(jìn)一步的研究，出現(xiàn)了預(yù)處理子空間迭代法，這種方法的運(yùn)算量較之子空間迭代法有所減小，但還是比較大，另外，當(dāng)要求的特征值與不要求的特征值之間的間隔較大時(shí)，預(yù)處理子空間迭代法收斂也比較慢。因此需要有一種更高效的方法來求解，今天給大家?guī)淼氖翘m索斯法(Lanczos法)振型求解思路方法。
關(guān)于子空間迭代法可以觀看下列文章：
【JY】振型求解之子空間迭代
【JY|STR】求解器之三維結(jié)構(gòu)振型分析
目前科學(xué)計(jì)算、理論研究、科學(xué)實(shí)驗(yàn)并列為當(dāng)今世界科學(xué)活動(dòng)的三種主要方式。許多科學(xué)和工程領(lǐng)域如果沒有科學(xué)計(jì)算不可能有一流的研究成果。而在工程領(lǐng)域，矩陣計(jì)算具有舉足輕重的地位。
在土木工程上，矩陣的特征值具有廣泛的應(yīng)用
，如建筑結(jié)構(gòu)中的固有頻率分析問題以及鋼結(jié)構(gòu)的穩(wěn)定性，建筑結(jié)構(gòu)的振動(dòng)問題，大型橋梁的顫振問題等等，都與矩陣的特征值密切相關(guān)。一個(gè)復(fù)雜結(jié)構(gòu)的動(dòng)力分析和穩(wěn)定性分析可轉(zhuǎn)化為大型稀疏矩陣的特征值/特征向量問題。
求解結(jié)構(gòu)振型和頻率實(shí)質(zhì)上是對(duì)以下式子的特征值和特征周期進(jìn)行求解。但是，根據(jù)伽羅瓦理論，對(duì)于次數(shù)大于四的多項(xiàng)式求根不存在一種通用的方法。換句話說，對(duì)于大型矩陣并無法直接求解得到通法解析解，從而得到特征值和特征向量。于是，人們通過尋求其他的求解途徑，提出了許多很好的思想方法和具體算法，促進(jìn)了該學(xué)科的不斷發(fā)展。
本期介紹一種針對(duì)稀疏矩陣高效的特征值、特征向量計(jì)算方法——Lanczos迭代法。Lanczos方法存儲(chǔ)量個(gè)大，計(jì)算速度較快，而且適用于求解矩陣的極端特征值。由于在實(shí)際問題中，矩陣的階數(shù)雖然很大，但往往只是需要求其依模最大或最小的特征值，因此Lanczos方法得以廣泛的應(yīng)用，特別是適合求大型稀疏對(duì)稱矩陣的部分特征值。
對(duì)于Lanczos方法的本質(zhì)
：將實(shí)對(duì)稱矩陣轉(zhuǎn)換成為三對(duì)角矩陣（稀疏矩陣），從而便于計(jì)算機(jī)儲(chǔ)存和后續(xù)的計(jì)算。（這樣就把一個(gè)求實(shí)對(duì)稱矩陣的特征值問題轉(zhuǎn)化為求一個(gè)三對(duì)角矩陣的特征值問題。）
有興趣的小伙伴也可以了解下另外兩種方法將一個(gè)實(shí)對(duì)稱矩陣轉(zhuǎn)化為三對(duì)角矩陣的方法Givens法及Householdes。
假
定該結(jié)構(gòu)特征矩陣是大型
、稀疏的實(shí)對(duì)稱矩陣，求其幾個(gè)最大或最小的特征值的問題，可以用Lanczos算法解決，它的原理是先產(chǎn)生一個(gè)三對(duì)角陣T，于是問題轉(zhuǎn)化為求一個(gè)三對(duì)角矩陣的特征值，變得相對(duì)簡便。
因此，我們的目的是將矩陣A轉(zhuǎn)化為T，并且只要求得T的特征值和特征向量，則可以通過一定關(guān)系得到原結(jié)構(gòu)特征矩陣的特征值和特征向量，也即可以得到結(jié)構(gòu)振型和頻率。
首先對(duì)于一般結(jié)構(gòu)來說，均需要進(jìn)行初始向量的預(yù)設(shè)進(jìn)行迭代，但是大部分振型都難以提前預(yù)估，我們可以聽下威爾金森(Wilkinson)的建議預(yù)先不必猜測，而統(tǒng)一初始的向量為：
然后進(jìn)行預(yù)設(shè)需求的振型數(shù)量為 i≤n (結(jié)構(gòu)特征矩陣維度)進(jìn)行求解，以下流程圖對(duì)于自己編寫Lanczos方法中的T矩陣集成有較好的理解幫助。
通過上述的迭代方法，并將求得的α和β帶入矩陣T中，即可將結(jié)構(gòu)特征矩陣是大型、稀疏的實(shí)對(duì)稱矩陣化為
一個(gè)三對(duì)角陣T。
那么結(jié)構(gòu)的頻率和振型與該三對(duì)角矩陣T的關(guān)系是什么呢？
接下來討論下T矩陣和結(jié)構(gòu)的頻率和振型的關(guān)系，繼上面的公式推導(dǎo)得到三對(duì)角矩陣T，由Gram-Schmidt正交化條件得到：
將Gram-Schmidt正交化條件帶入上式中可得到：
考慮列向量間的正交性，并將得到的上述公式往會(huì)帶入可得到：
得到這個(gè)標(biāo)注了三星的重要級(jí)公式！其中標(biāo)記紅框部分是左乘部分。若將紅框這么一畫，就變成了：
再通過上述這樣一個(gè)化解，原結(jié)構(gòu)和變換后的T的頻率和振型的關(guān)系一目了然。若令原結(jié)構(gòu)的特征向量（振型）為：
也就是求解T矩陣的完全解為，則該完全解的特征向量與原結(jié)構(gòu)的特征向量(振型)的關(guān)系式之間相差一個(gè)系數(shù)矩陣X，而特征值(頻率)是原結(jié)構(gòu)的倒數(shù)。因此取Lanczos向量的個(gè)數(shù)i等于系統(tǒng)的階數(shù)n，則Lanczos變換法給出了原特征值問題的精確解。
一般情況下我們?nèi)?i ≤ n 。由于組裝結(jié)構(gòu)特征矩陣時(shí)候，剛度求逆變換了一次，相當(dāng)于進(jìn)行了一次逆迭代，因此Lanczos變換法給出了系統(tǒng)低階特征值的很好的近似解。值得注意的是，為了保證求解精度，在迭代過程中可以用Gram-Schmidt對(duì)迭代向量進(jìn)行重新正交化，并采用移軸法可提高其效率。
至此，該問題變?yōu)榍笠粋€(gè)i階三對(duì)角陣T的特征值和特征向量。
對(duì)于求該降階的特征值、特征向量，有許多方法如：比較高效的QR分解法等等，這里就不再多贅述，小伙伴們可以自行編程嘗試。
建立一個(gè)無樓板的簡易框架結(jié)構(gòu)，體驗(yàn)一把Abaqus的蘭索斯算法，并對(duì)比下自編的Lanczos。模型如下：
求解前4階模態(tài)
利用Abaqus 直接提取該結(jié)構(gòu)的質(zhì)量矩陣、剛度矩陣、阻尼矩陣(非必要)，一個(gè)簡單的結(jié)構(gòu)密密麻麻的數(shù)，感興趣的小伙伴可以試一試：
并將這三個(gè)矩陣放入自行編寫的動(dòng)力計(jì)算器（由于該動(dòng)力計(jì)算器尚未完善，所以暫不公開使用，小伙伴們按這個(gè)思路也可以自己編寫得到屬于自己的求解器。）
可以看得到結(jié)果非常接近，誤差來源可能來源于計(jì)算機(jī)計(jì)算的舍入誤差。而且計(jì)算速度非常高效，因此針對(duì)稀疏矩陣高效的特征值、特征向量的計(jì)算方法可以采用可存儲(chǔ)量個(gè)大，計(jì)算速度較快的Lanczos迭代法。
