有限元法基礎-7線性代數(shù)方程組的解法

1、第七章 線性代數(shù)方程組的解法 7.1 高斯消去法及其變化7.2 帶狀系數(shù)矩陣的直接法7.3 利用外存的直接法7.4 迭代解法7. 線性代數(shù)方程組的解法 本章要點本章要點lGaussGauss消去法和三角分解法的原理和算法步驟消去法和三角分解法的原理和算法步驟l二維等帶寬存儲和一維變帶寬存儲的特點二維等帶寬存儲和一維變帶寬存儲的特點l分塊解法的原理和實施方案分塊解法的原理和實施方案l幾種迭代解法幾種迭代解法 有限元法基礎7. 線性代數(shù)方程組的解法 關(guān)鍵概念關(guān)鍵概念高斯循環(huán)消去法高斯循環(huán)消去法 三角分解法三角分解法二維等帶寬存儲二維等帶寬存儲 一維變帶寬存儲一維變帶寬存儲分塊解法分塊解法 迭代解法

2、迭代解法 超松弛迭代法超松弛迭代法梯度法梯度法 共軛梯度法共軛梯度法 預條件共軛梯度法預條件共軛梯度法 有限元法基礎7. 線性代數(shù)方程組的解法 彈性力學的有限元方程為彈性力學的有限元方程為對于彈性（本構(gòu)關(guān)系線性）小變形（幾何方程線性）問對于彈性（本構(gòu)關(guān)系線性）小變形（幾何方程線性）問題題K與與q無關(guān)，無關(guān)，為常數(shù)矩陣，方程組為線性代數(shù)方程組。為常數(shù)矩陣，方程組為線性代數(shù)方程組。求解是有限元方程分析中費時最多的步驟。求解是有限元方程分析中費時最多的步驟。 有限元法基礎KqQ7. 線性代數(shù)方程組的解法 線性代數(shù)方程組的解法分為兩大類，即線性代數(shù)方程組的解法分為兩大類，即直接解法直接解法和和迭迭代法

3、代法。 直接法的特點是，事先可按規(guī)定的算法步驟計算出它直接法的特點是，事先可按規(guī)定的算法步驟計算出它所需要的算術(shù)運算操作數(shù)，直接給出最后的結(jié)果。所需要的算術(shù)運算操作數(shù)，直接給出最后的結(jié)果。 迭代法的特點是，首先假定初始解，然后按一定的算迭代法的特點是，首先假定初始解，然后按一定的算法進行迭代，在每次的迭代過程檢查解的誤差，通過多法進行迭代，在每次的迭代過程檢查解的誤差，通過多次迭代直至滿足解的精度要求。次迭代直至滿足解的精度要求。 有限元法基礎7. 線性代數(shù)方程組的解法 有限元法基礎直接解法直接解法 以高斯消去法為基礎，以高斯消去法為基礎， 求解效率高，適用于小于求解效率高，適用于小于 一定階

4、數(shù)的方程組，根據(jù)計算機和軟件的不同有所一定階數(shù)的方程組，根據(jù)計算機和軟件的不同有所不同，比如不同，比如1 1萬萬1010萬階方程組。萬階方程組。迭代解法迭代解法當方程組階數(shù)過高當方程組階數(shù)過高 時，由于計算機有效位數(shù)的限時，由于計算機有效位數(shù)的限制，直接解法中舍入制，直接解法中舍入 誤差的積累影響精度，誤差的積累影響精度， 采用采用迭代迭代 解法。解法。 7.1 高斯消去法及其變化形式 有限元法基礎一一. .高斯循序消去法高斯循序消去法對于對于n n階線性方程組階線性方程組1.1.消元消元 1112111212222212nnnnnnnnkkkaPkkkaPkkkaP *( 1)1111211

5、*( 1)22222*( 1)000nnnnnnnnnakkkPakkPakP 7.1 高斯消去法及其變化形式 有限元法基礎對于對于n階線性方程組，共需進行階線性方程組，共需進行n-1次消元：次消元：l第第m次消元：次消元：以第以第m-1次消元結(jié)果為基礎次消元結(jié)果為基礎第第m行元素為消元行，行元素為消元行， 為主元為主元僅對僅對m+1 n行元素進行元素進 行元，并將行元，并將m列元素中列元素中 m+1 n消為消為0 0 1mmmK7.1 高斯消去法及其變化形式 有限元法基礎l對對i行行m列（列（im）消元）消元, ,將將m列從列從m+1列的元素列的元素消為消為0 0 稱為高斯消去因子稱為高斯消

6、去因子7.1 高斯消去法及其變化形式 有限元法基礎l因此消元過程可以寫為因此消元過程可以寫為 最終的最終的 為上三角陣。其中為上三角陣。其中(0)(1)(0)(1),nnKLKPLP(1)nK2131321211210010,1nnnlllll LLL LL(1)(1)(1)(1)(1)(1),1,2,/,1,2,1,nnniiiinnniiijiidKinKKKinji in為對角陣單位上三角陣KDKD7.1 高斯消去法及其變化形式 有限元法基礎l因此因此 因為因為K K(0)(0)為對稱矩陣，所以為對稱矩陣，所以 (0)(1)nKLDK(1)nTKL(0)TKLDL三角分解法的基礎三角分解

7、法的基礎7.1 高斯消去法及其變化形式 有限元法基礎l特點特點原系數(shù)矩陣是對稱的，則每次消元后矩陣依然是對稱原系數(shù)矩陣是對稱的，則每次消元后矩陣依然是對稱的，只需存儲一半的矩陣的，只需存儲一半的矩陣消元結(jié)果中，消元結(jié)果中， 和和 中的第中的第i i行就是（行就是（i-1)i-1)次次消元的結(jié)果消元的結(jié)果 (1)nP(1)nK(1)(1)(1)(1),1,2,1,niniijijiiKKPPinji in7.1 高斯消去法及其變化形式 有限元法基礎載荷列陣載荷列陣 消元用到的元素都是矩陣消元用到的元素都是矩陣 中的元素，中的元素，因此，因此， 的消元過程隨時可進行，對于多載荷工況，可的消元過程隨

8、時可進行，對于多載荷工況，可以利用消元后的以利用消元后的 矩陣進行消元和回代求解。這樣可矩陣進行消元和回代求解。這樣可大量節(jié)省求解所需的計算時間。大量節(jié)省求解所需的計算時間。 這是直接法相對迭代法的一個優(yōu)點。這是直接法相對迭代法的一個優(yōu)點。 P(1)nKPK7.1 高斯消去法及其變化形式 有限元法基礎2.2.回代求解回代求解回代公式回代公式 (1)(1)(1)(1)(1)11,2,2,1()/nnnnnnnnnniiijjiij iPaKinnaPKaK 7.1 高斯消去法及其變化形式 有限元法基礎l例例: :用高斯消元法求下列矩陣的解用高斯消元法求下列矩陣的解7.1 高斯消去法及其變化形式

9、有限元法基礎44320/7( 20/7)5/32,45/615/7aaa 34232112/5( 16/5)2( 4)3,214/55aaaaaa 回代求解得：回代求解得：7.1 高斯消去法及其變化形式 有限元法基礎二二. . 三角分解法三角分解法由高斯消去法能得到對由高斯消去法能得到對 的三角分解的三角分解設設KTKLDL下三角陣下三角陣對角陣對角陣TDLS上三角陣上三角陣7.1 高斯消去法及其變化形式 有限元法基礎由代數(shù)方程由代數(shù)方程可分解為可分解為高斯消元法相當于高斯消元法相當于令令KaPK = LS單位下三角陣單位下三角陣上三角陣上三角陣LSaP11L LSaL P1L PVP P在消

10、元后的結(jié)果在消元后的結(jié)果7.1 高斯消去法及其變化形式 有限元法基礎l三角分解后的代數(shù)方程求解步驟三角分解后的代數(shù)方程求解步驟KaPTLDL a = P1V = L P1TL a = D V1Ta = L D V7.1 高斯消去法及其變化形式 有限元法基礎l三角分解的遞推公式三角分解的遞推公式 K K 中任意元素中任意元素 TKLDL7.1 高斯消去法及其變化形式 有限元法基礎按行分解按行分解 i i = =1 1 i i = =2 2 7.1 高斯消去法及其變化形式 有限元法基礎 i i = =3,4,n3,4,n 7.1 高斯消去法及其變化形式 有限元法基礎 按行分解存儲情況按行分解存儲情

11、況7.1 高斯消去法及其變化形式 有限元法基礎按列分解按列分解 j=1j=1j j =2,=2,3,n3,n 7.1 高斯消去法及其變化形式 有限元法基礎 按列分解存儲情況按列分解存儲情況7.1 高斯消去法及其變化形式 有限元法基礎 l關(guān)于三角分解法關(guān)于三角分解法稱為改進稱為改進CholeskiCholeski法，經(jīng)典方法法，經(jīng)典方法比高斯消去法效率更高比高斯消去法效率更高只是改變了高斯消去法的循環(huán)循序和存儲只是改變了高斯消去法的循環(huán)循序和存儲 TKUU 按行三角分解按行三角分解 Do 15 i = 1, n Do 15 j = 1, n Do 15 m = 1, i-1 K(i,j )= K

12、(i,j) K(m,i) * K(m,j) /K(m,m)15 continue高斯循環(huán)消去法高斯循環(huán)消去法 Do 15 m = 1, n-1 Do 15 i = m+1, n Do 15 j = i, n K(i,j )= K(i,j) K(m,i) * K(m,j) /K(m,m)15 continue7.2 帶狀系數(shù)矩陣的直接法 有限元法基礎 系數(shù)矩陣在計算機中的存儲方法系數(shù)矩陣在計算機中的存儲方法l 等帶寬存儲等帶寬存儲K的特點：的特點：對稱、帶狀、稀疏對稱、帶狀、稀疏7.2 帶狀系數(shù)矩陣的直接法 有限元法基礎 二維等帶寬存儲二維等帶寬存儲(nND)7.2 帶狀系數(shù)矩陣的直接法 有限元

13、法基礎 l相關(guān)節(jié)點：相關(guān)節(jié)點：所有與節(jié)點所有與節(jié)點i i共單元的節(jié)點共單元的節(jié)點稱為節(jié)點稱為節(jié)點i i的相關(guān)節(jié)點的相關(guān)節(jié)點如果節(jié)點如果節(jié)點j j是節(jié)點是節(jié)點i i的相關(guān)節(jié)點的相關(guān)節(jié)點則則如果不是相關(guān)節(jié)點如果不是相關(guān)節(jié)點則則0ijK0ijK 7.2 帶狀系數(shù)矩陣的直接法 有限元法基礎 l一維變列高存儲主對角線位置主對角線位置MM：1 1，2 2，4 4，6 6，1010，1212，1616，1818，2222j j列上最上面非零元素行號列上最上面非零元素行號在一維存儲中得位置在一維存儲中得位置( )M jji7.2 帶狀系數(shù)矩陣的直接法 有限元法基礎 l兩種存儲方式比較二維等帶寬存儲二維等帶寬存

14、儲一維變帶寬存儲一維變帶寬存儲占內(nèi)存較多占內(nèi)存較多乘除法計算量相對較多乘除法計算量相對較多編程簡單編程簡單尋址時間較少尋址時間較少占內(nèi)存較少占內(nèi)存較少乘除法計算量相對較少乘除法計算量相對較少程序編制復雜程序編制復雜尋址時間較多尋址時間較多變列高變列高找元素找元素7.2 帶狀系數(shù)矩陣的直接法 有限元法基礎 二.二維等帶寬的高斯消去法l工作三角形工作三角形 由于系數(shù)矩陣呈帶狀由于系數(shù)矩陣呈帶狀 每次消元只涉及包括每次消元只涉及包括 主元在內(nèi)的一個三角主元在內(nèi)的一個三角 形內(nèi)的元素，稱為工形內(nèi)的元素，稱為工作三角形作三角形 。 7.2 帶狀系數(shù)矩陣的直接法 有限元法基礎 l二維等帶寬高斯消去法公式7

15、.2 帶狀系數(shù)矩陣的直接法 二維等帶寬存儲二維等帶寬存儲(n(nND)ND)采用按行分解采用按行分解I=i , J=j-i+1I=i , J=j-i+1新的循環(huán)界：新的循環(huán)界：r=max(j-ND+1, i-1)r=max(j-ND+1, i-1)有限元法基礎 l二維等帶寬三角分解7.2 帶狀系數(shù)矩陣的直接法 有限元法基礎 三.一維變列高存儲高斯消去法l采用按列分解7.3 利用計算機外存的直接法 有限元法基礎 主要解決計算機內(nèi)存不足的問題，充分利用外存保存分解后的系數(shù)矩陣與未分解的系數(shù)矩陣，以達到小內(nèi)存算大問題的目的。7.3 利用計算機外存的直接法 有限元法基礎 一.高斯消去法的特點1）第m次

16、消元過程中，所涉及的元素僅在三角形次消元過程中，所涉及的元素僅在三角形工作區(qū)內(nèi)工作區(qū)內(nèi)消元行元素消元行元素7.3 利用計算機外存的直接法 有限元法基礎 2）在）在m次消元過程中，前面的元素不再參加消元，次消元過程中，前面的元素不再參加消元，后面的元素尚未參加消元。后面的元素尚未參加消元。3）在整個消元過程中，工作區(qū)自上向下運動。）在整個消元過程中，工作區(qū)自上向下運動。 為分塊解法奠定基礎為分塊解法奠定基礎7.3 利用計算機外存的直接法 有限元法基礎 二二.分塊解法分塊解法 設允許使用內(nèi)存為設允許使用內(nèi)存為NANA=NQ ND行數(shù)NAND要求 在每一分塊，在每一分塊，NQNQNDND行集成完畢，

17、可進行消元修行集成完畢，可進行消元修正，最后的正，最后的NDND行進入到下一塊系數(shù)矩陣一起集成，行進入到下一塊系數(shù)矩陣一起集成，消元修正。消元修正。7.3 利用計算機外存的直接法 有限元法基礎 l 分塊解法的特點分塊解法的特點在每一分塊中，系數(shù)矩陣的元素是先集成后消在每一分塊中，系數(shù)矩陣的元素是先集成后消元修正元修正從求解的全過程看，系數(shù)矩陣的集成和消元修從求解的全過程看，系數(shù)矩陣的集成和消元修正是交替進行正是交替進行7.3 利用計算機外存的直接法 有限元法基礎 分塊解法簡單框圖分塊解法簡單框圖7.3 利用計算機外存的直接法 有限元法基礎 三三. .波前法（波前法（Front MethodFr

18、ont Method） 高斯循序消去法和三角分解法的求解規(guī)模與高斯循序消去法和三角分解法的求解規(guī)模與帶寬帶寬NDND有關(guān)。有些情況下帶寬會很大，占用內(nèi)存有關(guān)。有些情況下帶寬會很大，占用內(nèi)存很大，限制了計算機的求解能力。很大，限制了計算機的求解能力。 波前法和分塊解法基本思想都是基于對高斯波前法和分塊解法基本思想都是基于對高斯消去法的再分析，由先集成后消元修正，發(fā)展到消去法的再分析，由先集成后消元修正，發(fā)展到集成和消元修正交替進行。集成和消元修正交替進行。 7.3 利用計算機外存的直接法 有限元法基礎 三三. .波前法（波前法（Front MethodFront Method） 高斯循序消去法和

19、三角分解法的求解規(guī)模與高斯循序消去法和三角分解法的求解規(guī)模與帶寬帶寬NDND有關(guān)。有些情況下帶寬會很大，占用內(nèi)存有關(guān)。有些情況下帶寬會很大，占用內(nèi)存很大，限制了計算機的求解能力。很大，限制了計算機的求解能力。 波前法和分塊解法基本思想都是基于對高斯波前法和分塊解法基本思想都是基于對高斯消去法的再分析，由先集成后消元修正，發(fā)展到消去法的再分析，由先集成后消元修正，發(fā)展到集成和消元修正交替進行。集成和消元修正交替進行。 7.3 利用計算機外存的直接法 有限元法基礎 l波前法的特點波前法的特點1 1）剛度矩陣）剛度矩陣K K和載荷矩陣和載荷矩陣P P不按自然編號進入內(nèi)存，不按自然編號進入內(nèi)存，而是按

20、計算時參加運算的順序排列而是按計算時參加運算的順序排列2 2）在內(nèi)存中保留盡可能少的一部分）在內(nèi)存中保留盡可能少的一部分K K和和P P中的元素中的元素3 3）完成消元修正的行保存在外存）完成消元修正的行保存在外存 求解的自然編號是節(jié)點順序求解的自然編號是節(jié)點順序 參加運算的順序是單元順序參加運算的順序是單元順序7.3 利用計算機外存的直接法 有限元法基礎 內(nèi)存內(nèi)存最大工作三角塊（波前區(qū)）最大工作三角塊（波前區(qū)）三角形直角邊為波前數(shù)三角形直角邊為波前數(shù)W最大波前數(shù)最大波前數(shù)W ND需要存儲大量信息，以用于恢復完成集成，消元需要存儲大量信息，以用于恢復完成集成，消元的自由度號的自由度號 7.3

21、利用計算機外存的直接法 有限元法基礎 l波前法的計算過程波前法的計算過程1 1）按單元順序掃描計算單元剛度矩陣，并送入內(nèi)）按單元順序掃描計算單元剛度矩陣，并送入內(nèi)存進行集成存進行集成2 2）檢查那些）檢查那些DOFDOF已完成集成，將集成完畢的已完成集成，將集成完畢的DOFDOF作為主元，對其他行、列進行消元修正作為主元，對其他行、列進行消元修正3 3）完成消元修正后，將主）完成消元修正后，將主DOFDOF行有關(guān)的行有關(guān)的K K和和P P中的中的元素移到外存元素移到外存4 4）重復）重復1 13 3步，將全部單元掃描完畢步，將全部單元掃描完畢5 5）按消元順序，由后向前依次回代求解）按消元順序

22、，由后向前依次回代求解 7.3 利用計算機外存的直接法 有限元法基礎 波前法一度是有限元研究者廣泛采用的方法波前法一度是有限元研究者廣泛采用的方法與分塊解法相比，波前法利用內(nèi)存更少與分塊解法相比，波前法利用內(nèi)存更少由于頻繁使用內(nèi)外存交換求解效率低由于頻繁使用內(nèi)外存交換求解效率低編程復雜，以效率換取求解規(guī)模編程復雜，以效率換取求解規(guī)模隨著計算機硬件的發(fā)展，目前已較少應用隨著計算機硬件的發(fā)展，目前已較少應用 7.4 迭代法 有限元法基礎 一一. .雅克比迭代法雅克比迭代法方程組為方程組為方程組非奇異，且方程組非奇異，且 0(1,2, )iiain7.4 迭代法 有限元法基礎 方程組可改寫為方程組可

23、改寫為雅克比迭代法雅克比迭代法 設初始解設初始解 迭代方程迭代方程 0(1,2, )ixin7.4 迭代法 有限元法基礎 為了便于編程，方程組可改寫為為了便于編程，方程組可改寫為精度檢查準則精度檢查準則 為允許誤差為允許誤差當系數(shù)矩陣為嚴格對角優(yōu)勢矩陣時，方法收斂當系數(shù)矩陣為嚴格對角優(yōu)勢矩陣時，方法收斂 1(0,1,2,)kkkxxxk1(1,2, )niiijjj iaain7.4 迭代法 有限元法基礎 例例1 1：用雅克比迭代法求解方程組：用雅克比迭代法求解方程組 Ax=b,Ax=b,其中其中相對誤差控制為相對誤差控制為初始解取為初始解取為精確解為精確解為經(jīng)過經(jīng)過3030次迭代可達到精度要

24、求。次迭代可達到精度要求。 61000,0,0,0Tx 1,2,3,4Tx 7.4 迭代法 有限元法基礎 例例2 2：用雅克比迭代法求解方程組：用雅克比迭代法求解方程組 Ax=b,Ax=b,其中其中相對誤差控制為相對誤差控制為初始解取為初始解取為精確解為精確解為雅克比迭代法不收斂。原因：不滿足嚴格對角優(yōu)勢。雅克比迭代法不收斂。原因：不滿足嚴格對角優(yōu)勢。 61000,0,0,0Tx 1,2,3,4Tx 7.4 迭代法 有限元法基礎 二二. .高斯賽德爾高斯賽德爾(Gauss(Gauss-Seidel)迭代法迭代法雅克比迭代式雅克比迭代式在計算在計算 時，時， 為已知量，將迭為已知量，將迭代式修改

25、為代式修改為 1kix1(1,2,1)kjxji7.4 迭代法 有限元法基礎 三三. .超松弛超松弛迭代法迭代法 逐次超松弛迭代是逐次超松弛迭代是G-C迭代的一種加速收斂算法，迭代的一種加速收斂算法，引入松弛因子引入松弛因子 111G-C迭代法迭代法低松弛迭代法低松弛迭代法超松弛迭代法超松弛迭代法7.4 迭代法 有限元法基礎 三三. .超松弛超松弛迭代法迭代法 逐次超松弛迭代是逐次超松弛迭代是G-C迭代的一種加速收斂算法，迭代的一種加速收斂算法，引入松弛因子引入松弛因子 111G-C迭代法迭代法低松弛迭代法低松弛迭代法超松弛迭代法超松弛迭代法7.4 迭代法 有限元法基礎 例例3 3：松弛因子：

26、松弛因子迭代法求解例迭代法求解例1 其余條件與例其余條件與例1 相同，考察松弛因子對收斂相同，考察松弛因子對收斂性的影響。性的影響。 7.4 迭代法 有限元法基礎 例例4 4：松弛因子：松弛因子迭代法求解例迭代法求解例2 其余條件與例其余條件與例2 相同，考察松弛因子對收斂相同，考察松弛因子對收斂性的影響。性的影響。 可見，當系數(shù)矩陣不是嚴格的對角優(yōu)勢，算法也可見，當系數(shù)矩陣不是嚴格的對角優(yōu)勢，算法也收斂。通常取收斂。通常取 。 1.27.4 迭代法 有限元法基礎 四四. .共軛梯度法共軛梯度法 共軛梯度法由共軛梯度法由Hestenes和和Stiefel提出，可以提出，可以看作是稱為最速下降法

27、的梯度法發(fā)展而來，看作是稱為最速下降法的梯度法發(fā)展而來，簡稱簡稱CG法（法（Conjugate Gradient Method）。）。 7.4 迭代法 有限元法基礎 1.1.梯度法梯度法 對與很多數(shù)學物理問題，如果方程是自伴隨的，對與很多數(shù)學物理問題，如果方程是自伴隨的，可等效為求解對應的二次泛函的極值問題?？傻刃榍蠼鈱亩畏汉臉O值問題。線性方程組線性方程組 對應的二次泛函對應的二次泛函在在xk處的梯度為處的梯度為記負梯度方向記負梯度方向 1( )2TTfxx Ax -b xAx = b( )kkfx xxAx -bxkkrb- Ax7.4 迭代法 有限元法基礎 設設 xk 為為 的一

28、個近似解，利用最速下降法的一個近似解，利用最速下降法改善近似值：改善近似值：在在 的方向移動的方向移動xk ，即，即再選擇再選擇 使使 取得極小值，即令取得極小值，即令 Ax = b-grad f()x1kkk=-grad f()xxxk1()kfx()/0kk kkdf+dxrTkkkTkkr rr Ar7.4 迭代法 有限元法基礎 l迭代公式迭代公式初始解初始解 x0 0 可任意選取?？扇我膺x取。實際計算中發(fā)現(xiàn)收斂速度不高。實際計算中發(fā)現(xiàn)收斂速度不高。10,1,2,kkk kkkTkkkTkk=+kxxrrb- Axr rr Ar7.4 迭代法 有限元法基礎 2.2.共軛梯度法共軛梯度法 定義：向量定義：向量 pi i 和和 pj j , ,若若稱為關(guān)于稱為關(guān)于A正交或共軛。正交或共軛。 使用共軛梯度方向搜索近似解使用共軛梯度方向搜索近似解xk1 ，可加快收，可加快收斂速度。斂速度。 0()Tijijp Ap近似解為近似解為共軛梯度方向共軛梯度方向選擇選擇 求二次泛函極值求二次泛函極值7.4 迭代法 有限元法基礎 1kkkk=+xxp110kkkkTkkprpp Apk7.4 迭代法 有限元法基礎 l共軛梯度法迭代公式共軛梯度法迭代公式111111kkkkkkkkTkkkTk