大型龍門銑床主軸滑枕有限元分析說(shuō)明書(shū),大型,龍門,銑床,主軸,有限元分析,說(shuō)明書(shū),仿單 基于逆冪法的解決本地修改系統(tǒng)本征解的一種數(shù)值方法 K. Krishnapillai, R.Jones 摘要：快速再分析修改系統(tǒng)本征解的問(wèn)題具有相當(dāng)大的實(shí)際意義，現(xiàn)在已經(jīng)開(kāi)發(fā)出了用于計(jì)算特征值修改系統(tǒng)的幾種方法。本文提出的正是這種基于逆功率的方法，僅修改部分使用自由度。這種做法不論規(guī)模的修改和使用的濃縮程度的修改部分，都能使精確解的特征值盡快找到。提出的方法的優(yōu)勢(shì)是比較研究解決方案的逆功率方法的若干數(shù)值例子。這一做法將有助于改進(jìn)系統(tǒng)的自由度大和修改系統(tǒng)小的問(wèn)題。 關(guān)鍵詞：局部修改系統(tǒng) 逆冪法 特征值 特征向量 再分析 1 引言 振動(dòng)工程包括從機(jī)械震動(dòng)到電氣振蕩和擺動(dòng)的大結(jié)構(gòu)的廣泛的各種領(lǐng)域。動(dòng)態(tài)分析的特征值和特征向量在自由振動(dòng)的這種系統(tǒng)（動(dòng)力學(xué)）是一項(xiàng)基本的分支工程，而且大量的數(shù)值計(jì)算方法已經(jīng)被提出了有效地執(zhí)行這種分析。最近取得的進(jìn)展有限元法（ FEM ）和大型計(jì)算機(jī)，使很多分析家治療更大程度的自由（自由度） ，而且越來(lái)越多的軟件正在開(kāi)發(fā)，以適應(yīng)這些因素。但是，如果任何參數(shù)的變化（例如，形狀，材料，初始條件或環(huán)境條件）在先前分析系統(tǒng)，整個(gè)系統(tǒng)都必須再分析。特征值分析是一個(gè)更艱巨的任務(wù)，比較相應(yīng)的靜態(tài)分析，后者的計(jì)算時(shí)間一般是前者的幾倍甚至幾十倍。如果一個(gè)系統(tǒng)進(jìn)行多次修改，因?yàn)榉治鲂抻喓蟮南到y(tǒng)一樣需要時(shí)間，勞動(dòng)力和金錢，就像分析原始系統(tǒng)一樣，這可以使這一進(jìn)程相當(dāng)費(fèi)時(shí)。 “再分析”是用來(lái)執(zhí)行部分分析（而不是全面分析），以便由于設(shè)計(jì)修改或其他一些變化使部分矩陣被更改時(shí)而獲得有效數(shù)據(jù)。如果可以找到采用先前計(jì)算特征值和特征向量表演再分析的方便方法，這將是非常有用的解決涉及部分改變矩陣的特征值問(wèn)題的。再分析方法在許多文件中得到高度重視，這些文件涉及變化的特征值在本地修改系統(tǒng)中的問(wèn)題[1-9]。許多這些文件都使用微擾理論[1，4，7]。Fox and Kapoor [1]和Rogers [7]應(yīng)用一階微分方程的特征值和特征向量的設(shè)計(jì)變量。這種擾動(dòng)的方法是使設(shè)計(jì)變量的初始值假定不變而有相對(duì)較小的調(diào)整時(shí)的是有效解決方案。但是，如果改變?cè)O(shè)計(jì)變量大甚至是增加，這些方法一般是提供低精度的擾動(dòng)。 除Hirai等人[2]證實(shí)了只使用程序的修改部分以獲取準(zhǔn)確的特征值和特征向量證明以外，還有Parazzola等人[ 6，10 ]采用此方法做為特征值和特征向量的阻尼系統(tǒng)的理論解決方案。這些方法在確切的解決辦法的基礎(chǔ)上使用特征值和特征向量的修改制度來(lái)確定一個(gè)基本公式具有相同程度的矩陣代表后修改的制度。他們不論規(guī)模的修改都是有效的。如果特征值可再分析問(wèn)題使用這種凝聚方程，這將意味著，矩陣可以簡(jiǎn)化程度較低，這是一個(gè)減少計(jì)算時(shí)間非常有效的方法。這基本方程是非線性的，它可以用來(lái)解決一個(gè)矩陣方程的組成是否合理的問(wèn)題。但解決這類問(wèn)題的試驗(yàn)和錯(cuò)誤是非常耗時(shí)的，而且通過(guò)這種過(guò)程無(wú)法顯示順序特征。此外，很少有前景的實(shí)現(xiàn)穩(wěn)定的計(jì)算過(guò)程中使用Newton—Raphson法或其他方法，其中的差別系數(shù)為當(dāng)?shù)亟鉀Q辦法和初步估計(jì)而設(shè)置為適當(dāng)?shù)闹?，其解決辦法必須合理地接近實(shí)際的解決辦法（如果不是合理地接近實(shí)際的解決辦法則其進(jìn)程不可預(yù)測(cè)）。即使大致范圍的解決方案已經(jīng)確定，很難找到穩(wěn)定的解決方案，利用逐次逼近的方法，如反向線性插值或多項(xiàng)式逼近。使用這些方法可能提供一些解決辦法，但眾所周知，數(shù)值計(jì)算可以提供解決方案的一個(gè)隨機(jī)秩序。當(dāng)一個(gè)解決方案已經(jīng)被發(fā)現(xiàn)使用，例如，Newton—Raphson法，通貨緊縮是建立和重復(fù)的程序，以確定下一個(gè)特征。因此，傳統(tǒng)的方法解決這些非線性方程組有一些根本性的弱點(diǎn)。 Kashiwagi等人對(duì)Hirai等人[11–16]出版的方程壓縮版本進(jìn)行了系統(tǒng)的研究[ 11-16日] 。他們已經(jīng)解決了在當(dāng)?shù)亟?jīng)修改的制度中找到所有特征值的這個(gè)問(wèn)題，提出了找到低度系統(tǒng)為特征的組合Durand–Kerner[ 12 ]法與Newton法[ 16 ]和合理的功能[ 13 ]的可靠的解決方案。他們還表明，Sturm序列對(duì)某些壓縮的強(qiáng)烈非線性方程組[ 14,15 ] 是有用的 ，討論了如何確定該地區(qū)的Sturm 序列的特征值是否存在問(wèn)題，并說(shuō)明了使用Sturm序列一分為二的解決方案。Sturm序列方法是計(jì)算特征值常用的傳統(tǒng)計(jì)算方法，一般，它最適合已轉(zhuǎn)化為三對(duì)角矩陣的矩陣 [ 9 ] 。轉(zhuǎn)換的過(guò)程中，需要進(jìn)行大量的計(jì)算，一般以三對(duì)角矩陣形式和代表半數(shù)以上的特征值來(lái)計(jì)算確定。由Kashiwagi等人制定的Sturm序列法的版本需要查明所有的特征值和特征向量的修改制度，但并不需要轉(zhuǎn)化成三對(duì)角矩陣的形式，所以它是這一領(lǐng)域獨(dú)特的貢獻(xiàn)。 許多特征解決問(wèn)題，尋求公正的幾個(gè)特征值及其相應(yīng)的來(lái)自低度系統(tǒng)的載體。逆功率職能是這種系統(tǒng)用迭代方法來(lái)獲得解決方案的一個(gè)例子 [ 17 ] 。逆冪法是一個(gè)獲得解決方案的特征值分析的基本方法;許多分析都是根據(jù)它得到的。例如，子空間迭代[ 18,19 ]是一個(gè)確定本征值非常常用的方法。因此，我們必須早日完成工具箱的方法來(lái)確定特征值的本地修改的制度，但目前還沒(méi)有一種已提議的基于逆功率的職能。 本文提出的正是這種僅修改部分使用自由度的基于逆功率的方法。這種做法不論使用的濃縮程度的修改部分的規(guī)模都能盡快找到精確解的特征值。當(dāng)程度的矩陣低，計(jì)算時(shí)間很短，特別是當(dāng)逆冪函數(shù)相結(jié)合轉(zhuǎn)變的起源，所有的特征值和特征向量的修改系統(tǒng)從最小的本征解開(kāi)始對(duì)前幾個(gè)本征解（ 10個(gè)或更少的實(shí)際系統(tǒng)）相對(duì)較小的系統(tǒng)是否有效必須被了解。以下各節(jié)中描述的理論和算法證明了數(shù)值求解典型特征值問(wèn)題這一有效的做法。 2 本地修改系統(tǒng)的理論與算法的逆功率方法 本節(jié)描述修改系統(tǒng)的逆功率方法的理論和什么是本地修改系統(tǒng)轉(zhuǎn)向逆功率方法的理論。這些理論使只是濃縮版修改部分的本征解完全確定。 2.1 本地修改系統(tǒng)的逆功率方法 一般特征值問(wèn)題如下，假設(shè)一個(gè)n×n實(shí)對(duì)稱矩陣A和一個(gè)積極的實(shí)對(duì)稱矩陣B組： （1） 從最小的數(shù)值和是相應(yīng)的特征值特征向量開(kāi)始，是隨著那里的特征值。A和B有以下幾種形式： （2） （3） 如果是模式矩陣的特征向量包含的問(wèn)題，則用等式（1）表示 （4） 然后用下列關(guān)系表示： （5） （6） 這里，Ⅰ是n × n的矩陣，由等式 （ 5 ）可得 （7） 通常由有限元分析的自由結(jié)構(gòu)振動(dòng)和屈曲有兩個(gè)特征值。在自由振動(dòng)問(wèn)題中，A是剛度矩陣K，B是質(zhì)量矩陣M。在屈曲中， A是剛度矩陣K，B是幾何剛度矩陣KG。 根據(jù)當(dāng)?shù)氐男薷模?A是由A+A替代，B是由B+B替代的。讓我們假設(shè)A和B如下： （8） （9） 如果零要素從A變到B ，我們會(huì)得到和： （10） （11） 這是化簡(jiǎn)后的部分（m×m矩陣，在等式（8）和（9）中假設(shè)m=2） 。 我們考慮一個(gè)m × n的布爾矩陣（使等式（12）中m = 2） ： （12） 和由下式給出 （13） （14） 一般等式（1）表示的特征值問(wèn)題可以用來(lái)表示上述化簡(jiǎn)的當(dāng)?shù)鼐仃嚕? （15） 等式（15）可重寫(xiě)布爾矩陣形式： （16） 讓我們重寫(xiě)上述矩陣，分別用P代表n × n矩陣，Q代表n ×m矩陣，R代表m× n矩陣，代表m×m對(duì)稱矩陣。然后，我們可以寫(xiě)成 （17） 的逆矩陣是 （18） 這里，是修改部分元素對(duì)應(yīng)的m×m矩陣 。 用逆功率的方法，如果后面的特征值試驗(yàn)值為然后等式（16）和（18）可得 （19） 我們又可以寫(xiě)成 （20） 是下面等式（25）中的；是前迭代中獲得的近似特征向量?？紤]到式中的 （等式（9）中）可得 （21） 這里，和是未修改部分的矢量 （22） （23） 然后，由下式給出 （24） 這里 （25） 因此，我們可以用數(shù)值本身的表示來(lái)尋找該特征向量。同樣在等式（24）中我們只需要只用簡(jiǎn)明修改部分的自由度的化簡(jiǎn)特征向量： （26） （27） 一旦等式（26）中的被找出和保存，則由等式（24）所得出的自由度就可以計(jì)算出近似的特征向量。 2.2 本地修改系統(tǒng)的逆功率轉(zhuǎn)換法 讓我們指定當(dāng)前特征值和下一個(gè)特征值為。冪函數(shù)的收斂速度由||/||給出。當(dāng)原特征值為時(shí)，逆冪函數(shù)的收斂速度隨原值的改變由|－|/|－|給出，其值是一個(gè)比未修改的逆冪函數(shù)較低的值。 因此，選擇適當(dāng)?shù)闹翟谠俜治龅倪^(guò)程中收斂可能會(huì)被大大提高。在本節(jié)中，我們得出一個(gè)本地修改系統(tǒng)的轉(zhuǎn)換逆功率方法。等式（28）的基本方程具有跟一個(gè)完整的系統(tǒng)逆功率轉(zhuǎn)換方法相同的格式，因此預(yù)計(jì)將對(duì)本地修改系統(tǒng)是有效地。 此處重述等式（1）代表本地修改系統(tǒng)逆功率轉(zhuǎn)換方法作為一般特征值問(wèn)題： （28） 等式（7）也可被重新寫(xiě)為： （29） 我們按照上一節(jié)中給出的相同的方法： (30) (31) 這確保我們說(shuō)明地點(diǎn)時(shí)尋找到近似的特征向量。而且，由等式（30）可得到濃縮修改部分基礎(chǔ)上的近似特征向量： (32) (33) 一旦用等式（32）算出并保存，近似特征向量可以用等式（30）給予的自由度計(jì)算得到。 2.3 初始近似特征向量和初始近似特征值 這種重新分析的方法需要以往修改系統(tǒng)的特征向量和特征值的知識(shí)。我們討論如何將這些數(shù)據(jù)可用于生成隨著i的初始近似特征值和特征向量，我們由下式開(kāi)始： (34) 考慮到正?；臈l件為： 我們發(fā)現(xiàn)： (35) 表明如果我們考慮到，可得： (36) 初始的近似特征向量只占簡(jiǎn)明修改部分： (37) 此外，考慮到初始近似特征值為： 并考慮本地化的和，我們發(fā)現(xiàn) (38) 2.4 特征向量正?；吞卣髦? 我們現(xiàn)在描述如何隨著i近似特征向量正常化以及如何計(jì)算近似特征值。我們指定 (39) 和條件 由等式（24）、（30）和的本地化，我們可以得出 （40） 然后 （41） 并且初始近似特征向量根據(jù)濃縮修改部分 （42） 我們指定初始近似特征值為 并且由等式（24）、（30）和和本地化，我們可以得到為： （43） 這操作要求用代替 2.5 提取已取得的特征值 當(dāng)發(fā)現(xiàn)高階特征值，有必要提取使用逆冪函數(shù)已經(jīng)發(fā)現(xiàn)的本征解。革蘭氏Schmidt正交是用來(lái)配合等式（6）、（24）、（30）和的本地化： （44） 并且： （45） 方程僅由簡(jiǎn)明修改部分得到。如果我們定義的組成部分為，那么我們可以重寫(xiě)如下： （46） 2.6 算法 本節(jié)說(shuō)明一個(gè)根據(jù)擬議法創(chuàng)建算法的例子： （1）K＝0（初始值）。和版本的特征值，特征向量和m×m矩陣（m：一定程度的修改部分）可以得出，并且mq（數(shù)目的本征解）初始近似特征向量和特征值可以使用等式（37）和（38）算出?？梢杂玫仁剑?5）算出并保存。得出的特征值按升序重新排列。 （2）K＝K＋1（計(jì)算K的特征值）。 1、當(dāng)上一步計(jì)算的近似特征值和本步驟中的近似特征值之差的絕對(duì)值除以本步驟近似特征值所除得的值大于一些指定的值（在實(shí)驗(yàn)中用計(jì)算）則用逆功率近似法；當(dāng)小于時(shí)則用轉(zhuǎn)向逆功率近似法。在目前的數(shù)值試驗(yàn)中，轉(zhuǎn)向的起源（是上一步驟中的近似特征值）。 2、簡(jiǎn)明特征向量和值可由等式（32）和（33）計(jì)算得出。 3、（K－1）組的特征向量可由等式（45）計(jì)算得出。 4、簡(jiǎn)明近似特征向量的規(guī)范化使用等式（42）和（43），并可以計(jì)算出近似值。 5、收斂測(cè)試；如果計(jì)算值被收斂，則用等式（30）所得的自由度計(jì)算特征值，否則，該程序返回到步驟1。 （3）如果沒(méi)有得到指定的一些特征值，則該程序返回到步驟（2）。 3 數(shù)值實(shí)驗(yàn) 圖1. 平面框架模型（ A型） 本實(shí)驗(yàn)使用的是如圖1顯示的利用有限元法求解包含關(guān)節(jié)強(qiáng)度的結(jié)構(gòu)框架的一般特征值問(wèn)題。當(dāng)上一步的近似特征值與目前步驟中的近似特征值所確定的差的絕對(duì)值除以目前步驟中的近似特征值小于時(shí)，每一個(gè)特征值則被決定。所有的數(shù)值計(jì)算都用雙精度。計(jì)算機(jī)使用的是一個(gè)富士通FMV（賽揚(yáng)處理器2.40GHz，RAM768MB，Windows XP系統(tǒng)，F(xiàn)ujitsu Fortran和C Academic Package Ver.3）。由鋼筋混凝土制成的柱子和梁都具有楊氏模量和密度。柱子和梁的橫截面分別是800mm×800mm和400mm×800mm（寬×高）。柱子的長(zhǎng)度都是3000mm（落地式地板高度），梁的跨度都是6000mm。結(jié)構(gòu)中的每個(gè)節(jié)點(diǎn)有三個(gè)自由度（X方向的位移u ， Y方向的位移V和傾斜角度）。剛度矩陣A和質(zhì)量系數(shù)矩陣B中的元素值與變量的角度傾斜是對(duì)應(yīng)的，并允許是不同的層面。矩陣B表示分布式質(zhì)量，而不是集中質(zhì)量。因?yàn)橥ǔＴ诠こ虇?wèn)題中的需要，前18個(gè)特征值的計(jì)算可以用兩個(gè)建議的方法和逆冪函數(shù)的方法。表1給出了數(shù)值實(shí)驗(yàn)中的參數(shù)。為計(jì)算五個(gè)級(jí)別的自由度n，如圖1所顯示的修改了地點(diǎn)的①和②；這些列被刪除。當(dāng)只有列①被移動(dòng)時(shí)標(biāo)明“m=3”，當(dāng)列①②都被移動(dòng)時(shí)則標(biāo)明“m=6”。 表1 例舉參數(shù) 類型 mr nr n A 8 8 216 B 12 12 468 C 16 16 816 D 20 20 1260 E 24 24 1800 F 28 28 2436 mr:存儲(chǔ)數(shù)目，nr:跨度，n:自由度 表2 所得特征值的舉例（類型 A,m=3） 序號(hào) 提出的方法 逆冪法 1 0.272888641053 0.272888641053 2 2.712802279429 2.712802279430 3 4.132181003310 4.132181003310 4 8.949844231526 8.949844231525 5 17.110318347381 17.110318347409 6 18.455201910472 18.455201911020 7 18.694790497310 18.694790496682 8 20.124657349311 20.124657349413 9 21.272968144586 21.272968144488 10 22.391939618166 22.391939618438 11 23.161264697183 23.161264696928 12 27.052631104578 27.052631104671 13 28.734377712798 28.734377712946 14 29.783422706474 29.783422706224 15 34.877946938428 34.877946941330 16 35.081364019585 35.081364016670 17 40.495367536585 40.495367536877 18 43.004770987907 43.004770987659 圖2.計(jì)算時(shí)間為最低特征值18（m=3） 圖3. 迭代次數(shù)（類型E，m= 3 ） 圖4. 計(jì)算時(shí)間為18最低特征值（m=6） 圖5.迭代次數(shù)（類型E，m=3） 表2列舉了一個(gè)所取得結(jié)果的例子。很顯然，本文提出的這種方法能正確計(jì)算出特征值，因?yàn)樗鼈兣c其他計(jì)算方法的預(yù)測(cè)是一致的。圖 2和圖4提供了這些結(jié)果和逆冪函數(shù)方法所預(yù)測(cè)的比較（帶系數(shù)矩陣）和CPU所需要的時(shí)間。在m=3時(shí)，用逆功率方法解決本地系統(tǒng)修改部分所需要的時(shí)間大約為該逆功率方法解決整個(gè)系統(tǒng)所需時(shí)間的1/25，而且轉(zhuǎn)向逆功率方法解決本地系統(tǒng)修改部分所需要的時(shí)間大約為逆功率方法解決整個(gè)系統(tǒng)所需時(shí)間的1/100。在m = 6時(shí) ，用逆功率方法解決本地系統(tǒng)修改部分所需要的時(shí)間大約為逆功率方法解決整個(gè)系統(tǒng)所需時(shí)間的一半，而且轉(zhuǎn)向逆功率方法解決本地系統(tǒng)修改部分所需要的時(shí)間大約為逆功率方法解決整個(gè)系統(tǒng)所需要時(shí)間的1/30。 圖3和圖5比較了三種方法所需要的迭代次數(shù)。在m=3時(shí)，用逆功率方法解決本地系統(tǒng)修改部分的收斂迭代次數(shù)遠(yuǎn)小于用逆功率方法解決整個(gè)系統(tǒng)的收斂迭代次數(shù)，而在m=6時(shí)，幾乎沒(méi)有區(qū)別這兩種方法。本地修改部分的逆功率轉(zhuǎn)向方法需要比其他兩種方法更少的迭代次數(shù)，保證使用這種方法時(shí)有較低的計(jì)算。 4 總結(jié) 本文提出的理論研究，其目的是利用以前獲得的數(shù)據(jù)用于未修改系統(tǒng)，使修改系統(tǒng)再分析獲得其新的特征值和特征向量，是一種穩(wěn)定且有效的分析方法。這種方法的主要特點(diǎn)是經(jīng)過(guò)修改可以得到一個(gè)壓縮的低階矩陣，而不論修改的規(guī)模，這矩陣為修改后的特征值和特征向量提供了一個(gè)準(zhǔn)確的解決方案。這個(gè)較低階矩陣為大大縮短計(jì)算時(shí)間提供了解決方案。 在一個(gè)部分修改系統(tǒng)特征值分析的基本方程基礎(chǔ)上產(chǎn)生了逆功率法，并且這種方法的有效性和穩(wěn)定性證明了數(shù)值分析。在這一分析，也有人認(rèn)為，修改部分的轉(zhuǎn)向逆功率分析特別有效。它比修改部分和整個(gè)系統(tǒng)的未轉(zhuǎn)向逆功率分析更急劇地減少迭代而收斂。然而，數(shù)量的計(jì)算主要是確定原系統(tǒng)的自由度，特征值和特征向量，而且代表一定修改程度的矩陣中的數(shù)據(jù)必須在采用這種方法之前就已獲得。預(yù)計(jì)該方法也可與斯特姆序列或二等分法結(jié)合進(jìn)一步提高穩(wěn)定性。未來(lái)的研究需要提供這種算法的進(jìn)一步驗(yàn)證和處理理論方面的裁斷。 參考文獻(xiàn) [1] R.L. Fox, M.P. Kapoor, Rates of eigenvalues and eigenvectors, AIAA J. 6 (1968) 2426–2429. [2] I. Hirai, T. Yoshimura, K. Takamura, On a direct eigensolution of locally modified structure, Int. J. Numer. Methods Eng. 6 (1973) 441–442. [3] S.A. Jasbir, Survey of structural reanalysis technique, in: ASCE ST4 Proceedings Paper, vol. 12056, 1976, pp. 783–802. [4] U. Kirsch, Approximate structural reanalysis for optimization along a line, Int. J. Numer. Methods Eng. 18 (1982) 635–651. [5] L. Kitis, W.D. Pilky, I. Hirai, Reanalysis, in: Finite Element Handbook, McGraw- Hill, New York, 1987. [6] A.B. Parazzola, Vibration of locally modified mechanical and structural systems, Doctoral Thesis, University of Virginia, 1981. [7] L.C. Rogers, Derivatives of eigenvalues and eigenvectors, AIAA J. 8 (1970) 943–945. [8] B.P. Wang, W.D. Pilkey, Efficient reanalysis of locally modified structures, in: Proceedings of the First Chautaqua on Finite Element Modeling, Schaeffer Analysis, 1980. [9] J.H. Wilkinson, The Algebraic Eigenvalue Problem, Oxford University Press, Oxford, 1965. [10] B.P. Wang, A.B. Parazzola, W.D. Pilkey, Reanalysis Modal Synthesis and Dynamic Design in State-of-the-Art Surveys of Finite Element Methods, ASME, 1983. [11] I. Hirai, M. Kashiwagi, Derivatives of eigenvectors of locally modified structure, Int. J. Numer. Methods Eng. 11 (1977) 1769–1773. [12] M. Kashiwagi, A numerical method for eigensolution of locally modified systems by Durand–Kerner like method, J. Struct. Eng. 417 (1990) 71–77. [13] M. Kashiwagi, I. Hirai, T. Katayama, A numerical method for eigensolution of locally modified systems by a rational function approximation, J. Struct. Eng. 37A (1991) 271–278. [14] M. Kashiwagi, I. Hirai, S. Ohwaki, W.D. Pilkey, Stable eigensolution of locally modified systems based on the Sturm sequence property, Finite Elem. Ana. Des. 9 (1991) 133–139. [15] M. Kashiwagi, I. Hirai, T. Katayama, W.D. Pilkey, Sturm sequence recurrence formula of locally modified systems, Finite Elem. Anal. Des. 12 (1993) 31–36. [16] T. Katayama, S. Miyamura, M. Kashiwagi, I. Hirai, Eigensolution of locally modified systems by Newton-like method, J. Struct. Eng. 38A (1992) 303–310. [17] M. Kashiwagi, An eigensolution method of large sparse symmetric matrices by the CG method, Trans. Jpn. Soc Ind. Appl. Math. 15 (2005) 29–43. [18] K.J. Bathe, Finite Element Procedures in Engineering Analysis, Prentice-Hall, Englewood Cliffs, NJ, 1982. [19] B. Nour-omid, B.N. Parlett, R.L. Taylor, Lanczos versus subspace iteration for solution of eigenvalue problems, Int. J. Numer. Method Eng. 19 (1983) 859–971.