1、第七章 线性代数方程组的解法 7.1 高斯消去法及其变化7.2 带状系数矩阵的直接法7.3 利用外存的直接法7.4 迭代解法7.线性代数方程组的解法 本章要点本章要点lGaussGauss消去法和三角分解法的原理和算法步骤消去法和三角分解法的原理和算法步骤l二维等带宽存储和一维变带宽存储的特点二维等带宽存储和一维变带宽存储的特点l分块解法的原理和实施方案分块解法的原理和实施方案l几种迭代解法几种迭代解法 有限元法基础7.线性代数方程组的解法 关键概念关键概念高斯循环消去法高斯循环消去法 三角分解法三角分解法二维等带宽存储二维等带宽存储 一维变带宽存储一维变带宽存储分块解法分块解法 迭代解法迭代
2、解法 超松弛迭代法超松弛迭代法梯度法梯度法 共轭梯度法共轭梯度法 预条件共轭梯度法预条件共轭梯度法 有限元法基础7.线性代数方程组的解法 弹性力学的有限元方程为弹性力学的有限元方程为对于弹性（本构关系线性）小变形（几何方程线性）问对于弹性（本构关系线性）小变形（几何方程线性）问题题K与与q无关，无关，为常数矩阵，方程组为线性代数方程组。为常数矩阵，方程组为线性代数方程组。求解是有限元方程分析中费时最多的步骤。求解是有限元方程分析中费时最多的步骤。有限元法基础KqQ7.线性代数方程组的解法 线性代数方程组的解法分为两大类，即线性代数方程组的解法分为两大类，即直接解法直接解法和和迭迭代法代法。直接
3、法的特点是，事先可按规定的算法步骤计算出它直接法的特点是，事先可按规定的算法步骤计算出它所需要的算术运算操作数，直接给出最后的结果。所需要的算术运算操作数，直接给出最后的结果。迭代法的特点是，首先假定初始解，然后按一定的算迭代法的特点是，首先假定初始解，然后按一定的算法进行迭代，在每次的迭代过程检查解的误差，通过多法进行迭代，在每次的迭代过程检查解的误差，通过多次迭代直至满足解的精度要求。次迭代直至满足解的精度要求。有限元法基础7.线性代数方程组的解法 有限元法基础直接解法直接解法 以高斯消去法为基础，以高斯消去法为基础，求解效率高，适用于小于求解效率高，适用于小于 一定阶数的方程组，根据计算
4、机和软件的不同有所一定阶数的方程组，根据计算机和软件的不同有所不同，比如不同，比如1 1万万1010万阶方程组。万阶方程组。迭代解法迭代解法当方程组阶数过高当方程组阶数过高 时，由于计算机有效位数的限时，由于计算机有效位数的限制，直接解法中舍入制，直接解法中舍入 误差的积累影响精度，误差的积累影响精度，采用采用迭代迭代 解法。解法。7.1 高斯消去法及其变化形式 有限元法基础一一.高斯循序消去法高斯循序消去法对于对于n n阶线性方程组阶线性方程组1.1.消元消元 1112111212222212nnnnnnnnkkkaPkkkaPkkkaP *(1)1111211*(1)22222*(1)00
5、0nnnnnnnnnakkkPakkPakP 7.1 高斯消去法及其变化形式 有限元法基础对于对于n阶线性方程组，共需进行阶线性方程组，共需进行n-1次消元：次消元：l第第m次消元：次消元：以第以第m-1次消元结果为基础次消元结果为基础第第m行元素为消元行，行元素为消元行，为主元为主元仅对仅对m+1 n行元素进行元素进 行元，并将行元，并将m列元素中列元素中 m+1 n消为消为0 0 1mmmK7.1 高斯消去法及其变化形式 有限元法基础l对对i行行m列（列（im）消元）消元,将将m列从列从m+1列的元素列的元素消为消为0 0 称为高斯消去因子称为高斯消去因子7.1 高斯消去法及其变化形式 有
6、限元法基础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.1 高斯消去法及其变化形式
7、有限元法基础l特点特点原系数矩阵是对称的，则每次消元后矩阵依然是对称原系数矩阵是对称的，则每次消元后矩阵依然是对称的，只需存储一半的矩阵的，只需存储一半的矩阵消元结果中，消元结果中，和和 中的第中的第i i行就是（行就是（i-1)i-1)次次消元的结果消元的结果 (1)nP(1)nK(1)(1)(1)(1),1,2,1,niniijijiiKKPPinji in7.1 高斯消去法及其变化形式 有限元法基础载荷列阵载荷列阵 消元用到的元素都是矩阵消元用到的元素都是矩阵 中的元素，中的元素，因此，因此，的消元过程随时可进行，对于多载荷工况，可的消元过程随时可进行，对于多载荷工况，可以利用消元后的以
8、利用消元后的 矩阵进行消元和回代求解。这样可矩阵进行消元和回代求解。这样可大量节省求解所需的计算时间。大量节省求解所需的计算时间。这是直接法相对迭代法的一个优点。这是直接法相对迭代法的一个优点。P(1)nKPK7.1 高斯消去法及其变化形式 有限元法基础2.2.回代求解回代求解回代公式回代公式 (1)(1)(1)(1)(1)11,2,2,1()/nnnnnnnnnniiijjiij iPaKinnaPKaK 7.1 高斯消去法及其变化形式 有限元法基础l例例:用高斯消元法求下列矩阵的解用高斯消元法求下列矩阵的解7.1 高斯消去法及其变化形式 有限元法基础44320/7(20/7)5/32,45
9、/615/7aaa 34232112/5(16/5)2(4)3,214/55aaaaaa 回代求解得：回代求解得：7.1 高斯消去法及其变化形式 有限元法基础二二.三角分解法三角分解法由高斯消去法能得到对由高斯消去法能得到对 的三角分解的三角分解设设KTKLDL下三角阵下三角阵对角阵对角阵TDLS上三角阵上三角阵7.1 高斯消去法及其变化形式 有限元法基础由代数方程由代数方程可分解为可分解为高斯消元法相当于高斯消元法相当于令令KaPK=LS单位下三角阵单位下三角阵上三角阵上三角阵LSaP11L LSaL P1L PVP P在消元后的结果在消元后的结果7.1 高斯消去法及其变化形式 有限元法基础
10、l三角分解后的代数方程求解步骤三角分解后的代数方程求解步骤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 高斯消去法及其变化形式 有限元法基础 按行分解存储情况按行分解存储情况7.1 高斯消去法及其变化形式 有限元法基础按列分解按列分解 j=1j=1j j=2,=2,3
11、,n3,n 7.1 高斯消去法及其变化形式 有限元法基础 按列分解存储情况按列分解存储情况7.1 高斯消去法及其变化形式 有限元法基础 l关于三角分解法关于三角分解法称为改进称为改进CholeskiCholeski法，经典方法法，经典方法比高斯消去法效率更高比高斯消去法效率更高只是改变了高斯消去法的循环循序和存储只是改变了高斯消去法的循环循序和存储 TKUU 按行三角分解按行三角分解 Do 15 i=1,n Do 15 j=1,n Do 15 m=1,i-1 K(i,j)=K(i,j)K(m,i)*K(m,j)/K(m,m)15 continue高斯循环消去法高斯循环消去法 Do 15 m=1
12、,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 带状系数矩阵的直接法 有限元法基础 一.一.系数矩阵在计算机中的存储方法系数矩阵在计算机中的存储方法l 等带宽存储等带宽存储K的特点：的特点：对称、带状、稀疏对称、带状、稀疏7.2 带状系数矩阵的直接法 有限元法基础 二维等带宽存储二维等带宽存储(nND)7.2 带状系数矩阵的直接法 有限元法基础 l相关节点：相关节点：所有与节点所有与节点i i共单元的节点共单元的节点称为节点称为节点i i的相关节点的相关节点如果节点如果节点j j是节点
13、是节点i i的相关节点的相关节点则则如果不是相关节点如果不是相关节点则则0ijK0ijK 7.2 带状系数矩阵的直接法 有限元法基础 l一维变列高存储主对角线位置主对角线位置MM：1 1，2 2，4 4，6 6，1010，1212，1616，1818，2222j j列上最上面非零元素行号列上最上面非零元素行号在一维存储中得位置在一维存储中得位置()M jji7.2 带状系数矩阵的直接法 有限元法基础 l两种存储方式比较二维等带宽存储二维等带宽存储一维变带宽存储一维变带宽存储占内存较多占内存较多乘除法计算量相对较多乘除法计算量相对较多编程简单编程简单寻址时间较少寻址时间较少占内存较少占内存较少乘
14、除法计算量相对较少乘除法计算量相对较少程序编制复杂程序编制复杂寻址时间较多寻址时间较多变列高变列高找元素找元素7.2 带状系数矩阵的直接法 有限元法基础 二.二维等带宽的高斯消去法l工作三角形工作三角形 由于系数矩阵呈带状由于系数矩阵呈带状 每次消元只涉及包括每次消元只涉及包括 主元在内的一个三角主元在内的一个三角 形内的元素，称为工形内的元素，称为工作三角形作三角形。7.2 带状系数矩阵的直接法 有限元法基础 l二维等带宽高斯消去法公式7.2 带状系数矩阵的直接法 二维等带宽存储二维等带宽存储(n(nND)ND)采用按行分解采用按行分解I=i,J=j-i+1I=i,J=j-i+1新的循环界：
15、新的循环界：r=max(j-ND+1,i-1)r=max(j-ND+1,i-1)有限元法基础 l二维等带宽三角分解7.2 带状系数矩阵的直接法 有限元法基础 三.一维变列高存储高斯消去法l采用按列分解7.3 利用计算机外存的直接法 有限元法基础 主要解决计算机内存不足的问题，充分利用外存保存分解后的系数矩阵与未分解的系数矩阵，以达到小内存算大问题的目的。7.3 利用计算机外存的直接法 有限元法基础 一.高斯消去法的特点1）第m次消元过程中，所涉及的元素仅在三角形次消元过程中，所涉及的元素仅在三角形工作区内工作区内消元行元素消元行元素7.3 利用计算机外存的直接法 有限元法基础 2）在）在m次消
16、元过程中，前面的元素不再参加消元，次消元过程中，前面的元素不再参加消元，后面的元素尚未参加消元。后面的元素尚未参加消元。3）在整个消元过程中，工作区自上向下运动。）在整个消元过程中，工作区自上向下运动。为分块解法奠定基础为分块解法奠定基础7.3 利用计算机外存的直接法 有限元法基础 二二.分块解法分块解法 设允许使用内存为设允许使用内存为NANA=NQ ND行数NAND要求 在每一分块，在每一分块，NQNQNDND行集成完毕，可进行消元修行集成完毕，可进行消元修正，最后的正，最后的NDND行进入到下一块系数矩阵一起集成，行进入到下一块系数矩阵一起集成，消元修正。消元修正。7.3 利用计算机外存
17、的直接法 有限元法基础 l分块解法的特点分块解法的特点1)在每一分块中，系数矩阵的元素是先集成后消在每一分块中，系数矩阵的元素是先集成后消元修正元修正2)从求解的全过程看，系数矩阵的集成和消元修从求解的全过程看，系数矩阵的集成和消元修正是交替进行正是交替进行7.3 利用计算机外存的直接法 有限元法基础 分块解法简单框图分块解法简单框图7.3 利用计算机外存的直接法 有限元法基础 三三.波前法（波前法（Front MethodFront Method）高斯循序消去法和三角分解法的求解规模与高斯循序消去法和三角分解法的求解规模与带宽带宽NDND有关。有些情况下带宽会很大，占用内存有关。有些情况下带
18、宽会很大，占用内存很大，限制了计算机的求解能力。很大，限制了计算机的求解能力。波前法和分块解法基本思想都是基于对高斯波前法和分块解法基本思想都是基于对高斯消去法的再分析，由先集成后消元修正，发展到消去法的再分析，由先集成后消元修正，发展到集成和消元修正交替进行。集成和消元修正交替进行。7.3 利用计算机外存的直接法 有限元法基础 三三.波前法（波前法（Front MethodFront Method）高斯循序消去法和三角分解法的求解规模与高斯循序消去法和三角分解法的求解规模与带宽带宽NDND有关。有些情况下带宽会很大，占用内存有关。有些情况下带宽会很大，占用内存很大，限制了计算机的求解能力。很
19、大，限制了计算机的求解能力。波前法和分块解法基本思想都是基于对高斯波前法和分块解法基本思想都是基于对高斯消去法的再分析，由先集成后消元修正，发展到消去法的再分析，由先集成后消元修正，发展到集成和消元修正交替进行。集成和消元修正交替进行。7.3 利用计算机外存的直接法 有限元法基础 l波前法的特点波前法的特点1 1）刚度矩阵）刚度矩阵K K和载荷矩阵和载荷矩阵P P不按自然编号进入内存，不按自然编号进入内存，而是按计算时参加运算的顺序排列而是按计算时参加运算的顺序排列2 2）在内存中保留尽可能少的一部分）在内存中保留尽可能少的一部分K K和和P P中的元素中的元素3 3）完成消元修正的行保存在外
20、存）完成消元修正的行保存在外存 求解的自然编号是节点顺序求解的自然编号是节点顺序 参加运算的顺序是单元顺序参加运算的顺序是单元顺序7.3 利用计算机外存的直接法 有限元法基础 内存内存最大工作三角块（波前区）最大工作三角块（波前区）三角形直角边为波前数三角形直角边为波前数W最大波前数最大波前数W ND需要存储大量信息，以用于恢复完成集成，消元需要存储大量信息，以用于恢复完成集成，消元的自由度号的自由度号 7.3 利用计算机外存的直接法 有限元法基础 l波前法的计算过程波前法的计算过程1 1）按单元顺序扫描计算单元刚度矩阵，并送入内）按单元顺序扫描计算单元刚度矩阵，并送入内存进行集成存进行集成2
21、 2）检查那些）检查那些DOFDOF已完成集成，将集成完毕的已完成集成，将集成完毕的DOFDOF作为主元，对其他行、列进行消元修正作为主元，对其他行、列进行消元修正3 3）完成消元修正后，将主）完成消元修正后，将主DOFDOF行有关的行有关的K K和和P P中的中的元素移到外存元素移到外存4 4）重复）重复1 13 3步，将全部单元扫描完毕步，将全部单元扫描完毕5 5）按消元顺序，由后向前依次回代求解）按消元顺序，由后向前依次回代求解 7.3 利用计算机外存的直接法 有限元法基础 波前法一度是有限元研究者广泛采用的方法波前法一度是有限元研究者广泛采用的方法与分块解法相比，波前法利用内存更少与分
22、块解法相比，波前法利用内存更少由于频繁使用内外存交换求解效率低由于频繁使用内外存交换求解效率低编程复杂，以效率换取求解规模编程复杂，以效率换取求解规模随着计算机硬件的发展，目前已较少应用随着计算机硬件的发展，目前已较少应用 7.4 迭代法 有限元法基础 一一.雅克比迭代法雅克比迭代法方程组为方程组为方程组非奇异，且方程组非奇异，且 0(1,2,)iiain7.4 迭代法 有限元法基础 方程组可改写为方程组可改写为雅克比迭代法雅克比迭代法 设初始解设初始解 迭代方程迭代方程 0(1,2,)ixin7.4 迭代法 有限元法基础 为了便于编程，方程组可改写为为了便于编程，方程组可改写为精度检查准则精
23、度检查准则 为允许误差为允许误差当系数矩阵为严格对角优势矩阵时，方法收敛当系数矩阵为严格对角优势矩阵时，方法收敛 1(0,1,2,)kkkxxxk1(1,2,)niiijjj iaain7.4 迭代法 有限元法基础 例例1 1：用雅克比迭代法求解方程组：用雅克比迭代法求解方程组 Ax=b,Ax=b,其中其中相对误差控制为相对误差控制为初始解取为初始解取为精确解为精确解为经过经过3030次迭代可达到精度要求。次迭代可达到精度要求。61000,0,0,0Tx 1,2,3,4Tx 7.4 迭代法 有限元法基础 例例2 2：用雅克比迭代法求解方程组：用雅克比迭代法求解方程组 Ax=b,Ax=b,其中其
24、中相对误差控制为相对误差控制为初始解取为初始解取为精确解为精确解为雅克比迭代法不收敛。原因：不满足严格对角优势。雅克比迭代法不收敛。原因：不满足严格对角优势。61000,0,0,0Tx 1,2,3,4Tx 7.4 迭代法 有限元法基础 二二.高斯赛德尔高斯赛德尔(Gauss(Gauss-Seidel)迭代法迭代法雅克比迭代式雅克比迭代式在计算在计算 时，时，为已知量，将迭为已知量，将迭代式修改为代式修改为 1kix1(1,2,1)kjxji7.4 迭代法 有限元法基础 三三.超松弛超松弛迭代法迭代法 逐次超松弛迭代是逐次超松弛迭代是G-C迭代的一种加速收敛算法，迭代的一种加速收敛算法，引入松弛
25、因子引入松弛因子 111G-C迭代法迭代法低松弛迭代法低松弛迭代法超松弛迭代法超松弛迭代法7.4 迭代法 有限元法基础 三三.超松弛超松弛迭代法迭代法 逐次超松弛迭代是逐次超松弛迭代是G-C迭代的一种加速收敛算法，迭代的一种加速收敛算法，引入松弛因子引入松弛因子 111G-C迭代法迭代法低松弛迭代法低松弛迭代法超松弛迭代法超松弛迭代法7.4 迭代法 有限元法基础 例例3 3：松弛因子：松弛因子迭代法求解例迭代法求解例1 其余条件与例其余条件与例1 相同，考察松弛因子对收敛相同，考察松弛因子对收敛性的影响。性的影响。7.4 迭代法 有限元法基础 例例4 4：松弛因子：松弛因子迭代法求解例迭代法求
26、解例2 其余条件与例其余条件与例2 相同，考察松弛因子对收敛相同，考察松弛因子对收敛性的影响。性的影响。可见，当系数矩阵不是严格的对角优势，算法也可见，当系数矩阵不是严格的对角优势，算法也收敛。通常取收敛。通常取 。1.27.4 迭代法 有限元法基础 四四.共轭梯度法共轭梯度法 共轭梯度法由共轭梯度法由Hestenes和和Stiefel提出，可以提出，可以看作是称为最速下降法的梯度法发展而来，看作是称为最速下降法的梯度法发展而来，简称简称CG法（法（Conjugate Gradient Method）。）。7.4 迭代法 有限元法基础 1.1.梯度法梯度法 对与很多数学物理问题，如果方程是自伴
27、随的，对与很多数学物理问题，如果方程是自伴随的，可等效为求解对应的二次泛函的极值问题。可等效为求解对应的二次泛函的极值问题。线性方程组线性方程组 对应的二次泛函对应的二次泛函在在xk处的梯度为处的梯度为记负梯度方向记负梯度方向 1()2TTfxx Ax-b xAx=b()kkfx xxAx-bxkkrb-Ax7.4 迭代法 有限元法基础 设设 xk 为为 的一个近似解，利用最速下降法的一个近似解，利用最速下降法改善近似值：改善近似值：在在 的方向移动的方向移动xk ，即，即再选择再选择 使使 取得极小值，即令取得极小值，即令 Ax=b-grad f()x1kkk=-grad f()xxxk1(
28、)kfx()/0kk kkdf+dxrTkkkTkkr rr Ar7.4 迭代法 有限元法基础 l迭代公式迭代公式初始解初始解 x0 0 可任意选取。可任意选取。实际计算中发现收敛速度不高。实际计算中发现收敛速度不高。10,1,2,kkk kkkTkkkTkk=+kxxrrb-Axr rr Ar7.4 迭代法 有限元法基础 2.2.共轭梯度法共轭梯度法 定义：向量定义：向量 pi i 和和 pj j ,若若称为关于称为关于A正交或共轭。正交或共轭。使用共轭梯度方向搜索近似解使用共轭梯度方向搜索近似解xk1 ，可加快收，可加快收敛速度。敛速度。0()Tijijp Ap近似解为近似解为共轭梯度方向
29、共轭梯度方向选择选择 求二次泛函极值求二次泛函极值7.4 迭代法 有限元法基础 1kkkk=+xxp110kkkkTkkprpp Apk7.4 迭代法 有限元法基础 l共轭梯度法迭代公式共轭梯度法迭代公式111111kkkkkkkkTkkkTkkTkkkTkkkkkk=+xxpprppArpApp rp AprrAp00000任意取xrbAxpr0,1,2,k 7.4 迭代法 有限元法基础 l共轭梯度法特点共轭梯度法特点1）在迭代过程中，）在迭代过程中，rk 是相互正交的是相互正交的2）对）对n阶线性代数方程组通过阶线性代数方程组通过n次迭代得到的次迭代得到的rn一定一定是零是零3）xn必定是
30、原方程的解，由于有数值截断误差，计必定是原方程的解，由于有数值截断误差，计算时不一定得到精确解算时不一定得到精确解 另一加速算法为另一加速算法为PCG（Preconditioned conjugate gradient method）收敛更快。）收敛更快。7.4 迭代法 有限元法基础 3.预条件共轭梯度法预条件共轭梯度法(Preconditioned conjugate gradient method)该法简称该法简称 PCG法，是一种加速共轭梯度法。法，是一种加速共轭梯度法。PCG法通过引入预条件矩阵法通过引入预条件矩阵M，使方程系数矩阵的条件数，使方程系数矩阵的条件数降低，以达到提高收敛速
31、度的目的。降低，以达到提高收敛速度的目的。7.4 迭代法 有限元法基础 线性方程组线性方程组引入对称正定矩阵引入对称正定矩阵原方程转换为原方程转换为改写为改写为其中其中M称为预条件矩阵。当称为预条件矩阵。当M为为A的近似时，的近似时，接近单位接近单位矩阵，它的条件数近似为矩阵，它的条件数近似为1，然后用，然后用CG法求解新方程法求解新方程组。组。TMW W1TTWAW WxWbAx=bAx=b1,TTA=WAWx=Wx,b=WbA7.4 迭代法 有限元法基础 lPCG法迭代公式法迭代公式0000000任意取求解xrbAxMhrph0,1,2,k 11111111TkkkTkkkkkkkkkkkkTkkkTkkkkkk=+h rp ApxxprrApMhrhrh rphp7.4 迭代法 有限元法基础 lPCG法对法对M的要求的要求1）应是原系数矩阵）应是原系数矩阵A的近似的近似2）应容易求逆，并尽可能少占存储）应容易求逆，并尽可能少占存储例如例如如何选择如何选择M成为成为PCG法的关键问题法的关键问题。1122M=D=D D
