有限元法基础-7线性代数方程组的解法

1、第七章 线性代数方程组的解法 7.1 高斯消去法及其变化 7.2 带状系数矩阵的直接法 7.3 利用外存的直接法 7.4 迭代解法 7. 线性代数方程组的解法 本章要点本章要点 lGaussGauss消去法和三角分解法的原理和算法步骤消去法和三角分解法的原理和算法步骤 l二维等带宽存储和一维变带宽存储的特点二维等带宽存储和一维变带宽存储的特点 l分块解法的原理和实施方案分块解法的原理和实施方案 l几种迭代解法几种迭代解法 有限元法基础 7. 线性代数方程组的解法 关键概念关键概念 高斯循环消去法高斯循环消去法 三角分解法三角分解法 二维等带宽存储二维等带宽存储 一维变带宽存储一维变带宽存储 分

2、块解法分块解法 迭代解法迭代解法 超松弛迭代法超松弛迭代法 梯度法梯度法 共轭梯度法共轭梯度法 预条件共轭梯度法预条件共轭梯度法 有限元法基础 7. 线性代数方程组的解法 弹性力学的有限元方程为弹性力学的有限元方程为 对于弹性（本构关系线性）小变形（几何方程线性）问对于弹性（本构关系线性）小变形（几何方程线性）问 题题K与与q无关，无关，为常数矩阵，方程组为线性代数方程组。为常数矩阵，方程组为线性代数方程组。 求解是有限元方程分析中费时最多的步骤。求解是有限元方程分析中费时最多的步骤。 有限元法基础 KqQ 7. 线性代数方程组的解法 线性代数方程组的解法分为两大类，即线性代数方程组的解法分为

3、两大类，即直接解法直接解法和和迭迭 代法代法。 直接法的特点是，事先可按规定的算法步骤计算出它直接法的特点是，事先可按规定的算法步骤计算出它 所需要的算术运算操作数，直接给出最后的结果。所需要的算术运算操作数，直接给出最后的结果。 迭代法的特点是，首先假定初始解，然后按一定的算迭代法的特点是，首先假定初始解，然后按一定的算 法进行迭代，在每次的迭代过程检查解的误差，通过多法进行迭代，在每次的迭代过程检查解的误差，通过多 次迭代直至满足解的精度要求。次迭代直至满足解的精度要求。 有限元法基础 7. 线性代数方程组的解法 有限元法基础 直接解法直接解法 以高斯消去法为基础，以高斯消去法为基础， 求

4、解效率高，适用于小于求解效率高，适用于小于 一定阶数的方程组，根据计算机和软件的不同有所一定阶数的方程组，根据计算机和软件的不同有所 不同，比如不同，比如1 1万万1010万阶方程组。万阶方程组。 迭代解法迭代解法 当方程组阶数过高当方程组阶数过高 时，由于计算机有效位数的限时，由于计算机有效位数的限 制，直接解法中舍入制，直接解法中舍入 误差的积累影响精度，误差的积累影响精度， 采用采用 迭代迭代 解法。解法。 7.1 高斯消去法及其变化形式 有限元法基础 一一. .高斯循序消去法高斯循序消去法 对于对于n n阶线性方程组阶线性方程组 1.1.消元消元 1112111 2122222 12

5、n n nnnnnn kkkaP kkkaP kkkaP *( 1) 1 111211 *( 1) 2 2222 *( 1) 0 00 n n n n n n nnn akkkP akkP akP 7.1 高斯消去法及其变化形式 有限元法基础 对于对于n阶线性方程组，共需进行阶线性方程组，共需进行n-1次消元：次消元： l第第m次消元：次消元： 以第以第m-1次消元结果为基础次消元结果为基础 第第m行元素为消元行，行元素为消元行， 为主元为主元 仅对仅对m+1 n行元素进行元素进 行元，并将行元，并将m列元素中列元素中 m+1 n消为消为0 0 1m mm K 7.1 高斯消去法及其变化形式

6、有限元法基础 l对对i行行m列（列（im）消元）消元, ,将将m列从列从m+1列的元素列的元素消为消为0 0 称为高斯消去因子称为高斯消去因子 7.1 高斯消去法及其变化形式 有限元法基础 l因此消元过程可以写为因此消元过程可以写为 最终的最终的 为上三角阵。其中为上三角阵。其中 (0)(1)(0)(1) , nn KLKPLP (1)n K 21 3132121 12 100 10 , 1 n nn l ll ll LLL LL (1)(1)(1) (1)(1)(1) ,1,2, /,1,2,1, nnn iiii nnn iiijii dKin KKKinji in 为对角阵 单位上三角阵

7、 KDKD 7.1 高斯消去法及其变化形式 有限元法基础 l因此因此 因为因为K K(0) (0)为对称矩阵，所以 为对称矩阵，所以 (0)(1)n KLDK (1)nT KL (0)T KLDL 三角分解法的基础三角分解法的基础 7.1 高斯消去法及其变化形式 有限元法基础 l特点特点 原系数矩阵是对称的，则每次消元后矩阵依然是对称原系数矩阵是对称的，则每次消元后矩阵依然是对称 的，只需存储一半的矩阵的，只需存储一半的矩阵 消元结果中，消元结果中， 和和 中的第中的第i i行就是（行就是（i-1)i-1)次次 消元的结果消元的结果 (1)n P (1)n K (1)(1)(1)(1) , 1

8、,2,1, nini ijijii KKPP inji in 7.1 高斯消去法及其变化形式 有限元法基础 载荷列阵载荷列阵 消元用到的元素都是矩阵消元用到的元素都是矩阵 中的元素，中的元素， 因此，因此， 的消元过程随时可进行，对于多载荷工况，可的消元过程随时可进行，对于多载荷工况，可 以利用消元后的以利用消元后的 矩阵进行消元和回代求解。这样可矩阵进行消元和回代求解。这样可 大量节省求解所需的计算时间。大量节省求解所需的计算时间。 这是直接法相对迭代法的一个优点。这是直接法相对迭代法的一个优点。 P (1)n K P K 7.1 高斯消去法及其变化形式 有限元法基础 2.2.回代求解回代求

9、解 回代公式回代公式 (1) (1) (1)(1)(1) 1 1,2,2,1 ()/ n n n n nn n nnn iiijjii j i P a K inn aPKaK 7.1 高斯消去法及其变化形式 有限元法基础 l例例: :用高斯消元法求下列矩阵的解用高斯消元法求下列矩阵的解 7.1 高斯消去法及其变化形式 有限元法基础 4 43 20/7( 20/7)5/3 2,4 5/615/7 a aa 3423 21 12/5( 16/5)2( 4) 3,2 14/55 aaaa aa 回代求解得：回代求解得： 7.1 高斯消去法及其变化形式 有限元法基础 二二. . 三角分解法三角分解法

10、由高斯消去法能得到对由高斯消去法能得到对 的三角分解的三角分解 设设 K T KLDL 下三角阵下三角阵对角阵对角阵 T DLS 上三角阵上三角阵 7.1 高斯消去法及其变化形式 有限元法基础 由代数方程由代数方程 可分解为可分解为 高斯消元法相当于高斯消元法相当于 令令 KaP K = LS 单位下三角阵单位下三角阵上三角阵上三角阵 LSaP 11 L LSaL P 1 L PVP P在消元后的结果在消元后的结果 7.1 高斯消去法及其变化形式 有限元法基础 l三角分解后的代数方程求解步骤三角分解后的代数方程求解步骤 KaP T LDL a = P 1 V = L P 1T L a = D

11、V 1T a = L D V 7.1 高斯消去法及其变化形式 有限元法基础 l三角分解的递推公式三角分解的递推公式 K K 中任意元素中任意元素 T KLDL 7.1 高斯消去法及其变化形式 有限元法基础 按行分解按行分解 i i = =1 1 i i = =2 2 7.1 高斯消去法及其变化形式 有限元法基础 i i = =3,4,n3,4,n 7.1 高斯消去法及其变化形式 有限元法基础 按行分解存储情况按行分解存储情况 7.1 高斯消去法及其变化形式 有限元法基础 按列分解按列分解 j=1j=1 j j =2,=2,3,n3,n 7.1 高斯消去法及其变化形式 有限元法基础 按列分解存储

12、情况按列分解存储情况 7.1 高斯消去法及其变化形式 有限元法基础 l关于三角分解法关于三角分解法 称为改进称为改进CholeskiCholeski法，经典方法法，经典方法 比高斯消去法效率更高比高斯消去法效率更高 只是改变了高斯消去法的循环循序和存储只是改变了高斯消去法的循环循序和存储 T KUU 按行三角分解按行三角分解 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, n-1 Do 15

13、 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 continue 7.2 带状系数矩阵的直接法 有限元法基础 一.一. 系数矩阵在计算机中的存储方法系数矩阵在计算机中的存储方法 l 等带宽存储等带宽存储 K的特点：的特点： 对称、带状、稀疏对称、带状、稀疏 7.2 带状系数矩阵的直接法 有限元法基础 二维等带宽存储二维等带宽存储 (nND) 7.2 带状系数矩阵的直接法 有限元法基础 l相关节点：相关节点： 所有与节点所有与节点i i共单元的节点共单元的节点 称为节点称为节点i i的相关节点的相关节点

14、如果节点如果节点j j是节点是节点i i的相关节点的相关节点 则则 如果不是相关节点如果不是相关节点 则则 0 ij K 0 ij K 7.2 带状系数矩阵的直接法 有限元法基础 l一维变列高存储 主对角线位置主对角线位置MM：1 1，2 2，4 4，6 6，1010，1212，1616，1818，2222 j j列上最上面非零元素行号列上最上面非零元素行号 在一维存储中得位置在一维存储中得位置 ( )M jji 7.2 带状系数矩阵的直接法 有限元法基础 l两种存储方式比较 二维等带宽存储二维等带宽存储一维变带宽存储一维变带宽存储 占内存较多占内存较多 乘除法计算量相对较多乘除法计算量相对较

15、多 编程简单编程简单 寻址时间较少寻址时间较少 占内存较少占内存较少 乘除法计算量相对较少乘除法计算量相对较少 程序编制复杂程序编制复杂 寻址时间较多寻址时间较多 变列高变列高 找元素找元素 7.2 带状系数矩阵的直接法 有限元法基础 二.二维等带宽的高斯消去法 l工作三角形工作三角形 由于系数矩阵呈带状由于系数矩阵呈带状 每次消元只涉及包括每次消元只涉及包括 主元在内的一个三角主元在内的一个三角 形内的元素，称为工形内的元素，称为工 作三角形作三角形 。 7.2 带状系数矩阵的直接法 有限元法基础 l二维等带宽高斯消去法公式 7.2 带状系数矩阵的直接法 二维等带宽存储二维等带宽存储 (n(

16、nND)ND) 采用按行分解采用按行分解 I=i , J=j-i+1I=i , J=j-i+1 新的循环界：新的循环界： r=max(j-ND+1, i-1)r=max(j-ND+1, i-1) 有限元法基础 l二维等带宽三角分解 7.2 带状系数矩阵的直接法 有限元法基础 三.一维变列高存储高斯消去法 l采用按列分解 7.3 利用计算机外存的直接法 有限元法基础 主要解决计算机内存不足的问题，充分 利用外存保存分解后的系数矩阵与未分解 的系数矩阵，以达到小内存算大问题的目 的。 7.3 利用计算机外存的直接法 有限元法基础 一.高斯消去法的特点 1）第m次消元过程中，所涉及的元素仅在三角形次

17、消元过程中，所涉及的元素仅在三角形 工作区内工作区内 消元行元素消元行元素 7.3 利用计算机外存的直接法 有限元法基础 2）在）在m次消元过程中，前面的元素不再参加消元，次消元过程中，前面的元素不再参加消元， 后面的元素尚未参加消元。后面的元素尚未参加消元。 3）在整个消元过程中，工作区自上向下运动。）在整个消元过程中，工作区自上向下运动。 为分块解法奠定基础为分块解法奠定基础 7.3 利用计算机外存的直接法 有限元法基础 二二.分块解法分块解法 设允许使用内存为设允许使用内存为NA NA =NQ ND 行数 NAND要求 在每一分块，在每一分块，NQNQNDND行集成完毕，可进行消元修行集

18、成完毕，可进行消元修 正，最后的正，最后的NDND行进入到下一块系数矩阵一起集成，行进入到下一块系数矩阵一起集成， 消元修正。消元修正。 7.3 利用计算机外存的直接法 有限元法基础 l分块解法的特点分块解法的特点 1) 在每一分块中，系数矩阵的元素是先集成后消在每一分块中，系数矩阵的元素是先集成后消 元修正元修正 2) 从求解的全过程看，系数矩阵的集成和消元修从求解的全过程看，系数矩阵的集成和消元修 正是交替进行正是交替进行 7.3 利用计算机外存的直接法 有限元法基础 分块解法简单框图分块解法简单框图 7.3 利用计算机外存的直接法 有限元法基础 三三. .波前法（波前法（Front Me

19、thodFront Method） 高斯循序消去法和三角分解法的求解规模与高斯循序消去法和三角分解法的求解规模与 带宽带宽NDND有关。有些情况下带宽会很大，占用内存有关。有些情况下带宽会很大，占用内存 很大，限制了计算机的求解能力。很大，限制了计算机的求解能力。 波前法和分块解法基本思想都是基于对高斯波前法和分块解法基本思想都是基于对高斯 消去法的再分析，由先集成后消元修正，发展到消去法的再分析，由先集成后消元修正，发展到 集成和消元修正交替进行。集成和消元修正交替进行。 7.3 利用计算机外存的直接法 有限元法基础 三三. .波前法（波前法（Front MethodFront Method

20、） 高斯循序消去法和三角分解法的求解规模与高斯循序消去法和三角分解法的求解规模与 带宽带宽NDND有关。有些情况下带宽会很大，占用内存有关。有些情况下带宽会很大，占用内存 很大，限制了计算机的求解能力。很大，限制了计算机的求解能力。 波前法和分块解法基本思想都是基于对高斯波前法和分块解法基本思想都是基于对高斯 消去法的再分析，由先集成后消元修正，发展到消去法的再分析，由先集成后消元修正，发展到 集成和消元修正交替进行。集成和消元修正交替进行。 7.3 利用计算机外存的直接法 有限元法基础 l波前法的特点波前法的特点 1 1）刚度矩阵）刚度矩阵K K和载荷矩阵和载荷矩阵P P不按自然编号进入内存

21、，不按自然编号进入内存， 而是按计算时参加运算的顺序排列而是按计算时参加运算的顺序排列 2 2）在内存中保留尽可能少的一部分）在内存中保留尽可能少的一部分K K和和P P中的元素中的元素 3 3）完成消元修正的行保存在外存）完成消元修正的行保存在外存 求解的自然编号是节点顺序求解的自然编号是节点顺序 参加运算的顺序是单元顺序参加运算的顺序是单元顺序 7.3 利用计算机外存的直接法 有限元法基础 内存内存最大工作三角块（波前区）最大工作三角块（波前区） 三角形直角边为波前数三角形直角边为波前数W 最大波前数最大波前数W ND 需要存储大量信息，以用于恢复完成集成，消元需要存储大量信息，以用于恢复

22、完成集成，消元 的自由度号的自由度号 7.3 利用计算机外存的直接法 有限元法基础 l波前法的计算过程波前法的计算过程 1 1）按单元顺序扫描计算单元刚度矩阵，并送入内）按单元顺序扫描计算单元刚度矩阵，并送入内 存进行集成存进行集成 2 2）检查那些）检查那些DOFDOF已完成集成，将集成完毕的已完成集成，将集成完毕的DOFDOF 作为主元，对其他行、列进行消元修正作为主元，对其他行、列进行消元修正 3 3）完成消元修正后，将主）完成消元修正后，将主DOFDOF行有关的行有关的K K和和P P中的中的 元素移到外存元素移到外存 4 4）重复）重复1 13 3步，将全部单元扫描完毕步，将全部单元

23、扫描完毕 5 5）按消元顺序，由后向前依次回代求解）按消元顺序，由后向前依次回代求解 7.3 利用计算机外存的直接法 有限元法基础 波前法一度是有限元研究者广泛采用的方法波前法一度是有限元研究者广泛采用的方法 与分块解法相比，波前法利用内存更少与分块解法相比，波前法利用内存更少 由于频繁使用内外存交换求解效率低由于频繁使用内外存交换求解效率低 编程复杂，以效率换取求解规模编程复杂，以效率换取求解规模 随着计算机硬件的发展，目前已较少应用随着计算机硬件的发展，目前已较少应用 7.4 迭代法 有限元法基础 一一. .雅克比迭代法雅克比迭代法 方程组为方程组为 方程组非奇异，且方程组非奇异，且 0(

24、1,2, ) ii ain 7.4 迭代法 有限元法基础 方程组可改写为方程组可改写为 雅克比迭代法雅克比迭代法 设初始解设初始解 迭代方程迭代方程 0 (1,2, ) i xin 7.4 迭代法 有限元法基础 为了便于编程，方程组可改写为为了便于编程，方程组可改写为 精度检查准则精度检查准则 为允许误差为允许误差 当系数矩阵为严格对角优势矩阵时，方法收敛当系数矩阵为严格对角优势矩阵时，方法收敛 1 (0,1,2,) kkk xxxk 1 (1,2, ) n iiij j j i aain 7.4 迭代法 有限元法基础 例例1 1：用雅克比迭代法求解方程组：用雅克比迭代法求解方程组 Ax=b,

25、Ax=b,其中其中 相对误差控制为相对误差控制为 初始解取为初始解取为 精确解为精确解为 经过经过3030次迭代可达到精度要求。次迭代可达到精度要求。 6 10 0 0,0,0,0 T x 1,2,3,4 T x 7.4 迭代法 有限元法基础 例例2 2：用雅克比迭代法求解方程组：用雅克比迭代法求解方程组 Ax=b,Ax=b,其中其中 相对误差控制为相对误差控制为 初始解取为初始解取为 精确解为精确解为 雅克比迭代法不收敛。原因：不满足严格对角优势。雅克比迭代法不收敛。原因：不满足严格对角优势。 6 10 0 0,0,0,0 T x 1,2,3,4 T x 7.4 迭代法 有限元法基础 二二.

26、 .高斯赛德尔高斯赛德尔(Gauss(Gauss-Seidel)迭代法迭代法 雅克比迭代式雅克比迭代式 在计算在计算 时，时， 为已知量，将迭为已知量，将迭 代式修改为代式修改为 1k i x 1 (1,2,1) k j xji 7.4 迭代法 有限元法基础 三三. .超松弛超松弛迭代法迭代法 逐次超松弛迭代是逐次超松弛迭代是G-C迭代的一种加速收敛算法，迭代的一种加速收敛算法， 引入松弛因子引入松弛因子 1 1 1 G-C迭代法迭代法 低松弛迭代法低松弛迭代法 超松弛迭代法超松弛迭代法 7.4 迭代法 有限元法基础 三三. .超松弛超松弛迭代法迭代法 逐次超松弛迭代是逐次超松弛迭代是G-C迭

27、代的一种加速收敛算法，迭代的一种加速收敛算法， 引入松弛因子引入松弛因子 1 1 1 G-C迭代法迭代法 低松弛迭代法低松弛迭代法 超松弛迭代法超松弛迭代法 7.4 迭代法 有限元法基础 例例3 3：松弛因子：松弛因子迭代法求解例迭代法求解例1 其余条件与例其余条件与例1 相同，考察松弛因子对收敛相同，考察松弛因子对收敛 性的影响。性的影响。 7.4 迭代法 有限元法基础 例例4 4：松弛因子：松弛因子迭代法求解例迭代法求解例2 其余条件与例其余条件与例2 相同，考察松弛因子对收敛相同，考察松弛因子对收敛 性的影响。性的影响。 可见，当系数矩阵不是严格的对角优势，算法也可见，当系数矩阵不是严格

28、的对角优势，算法也 收敛。通常取收敛。通常取 。 1.2 7.4 迭代法 有限元法基础 四四. .共轭梯度法共轭梯度法 共轭梯度法由共轭梯度法由Hestenes和和Stiefel提出，可以提出，可以 看作是称为最速下降法的梯度法发展而来，看作是称为最速下降法的梯度法发展而来， 简称简称CG法（法（Conjugate Gradient Method）。）。 7.4 迭代法 有限元法基础 1.1.梯度法梯度法 对与很多数学物理问题，如果方程是自伴随的，对与很多数学物理问题，如果方程是自伴随的， 可等效为求解对应的二次泛函的极值问题。可等效为求解对应的二次泛函的极值问题。 线性方程组线性方程组 对应

29、的二次泛函对应的二次泛函 在在xk处的梯度为处的梯度为 记负梯度方向记负梯度方向 1 ( ) 2 TT fxx Ax -b x Ax = b ( ) k k f x x x Ax -b x kk rb- Ax 7.4 迭代法 有限元法基础 设设 xk 为为 的一个近似解，利用最速下降法的一个近似解，利用最速下降法 改善近似值：改善近似值： 在在 的方向移动的方向移动xk ，即，即 再选择再选择 使使 取得极小值，即令取得极小值，即令 Ax = b -grad f()x 1 kkk =-grad f() xxx k 1 () k f x ()/0 kk kk df+dxr T kk k T kk

30、 r r r Ar 7.4 迭代法 有限元法基础 l迭代公式迭代公式 初始解初始解 x0 0 可任意选取。可任意选取。 实际计算中发现收敛速度不高。实际计算中发现收敛速度不高。 1 0,1,2, kkk k kk T kk k T kk =+k xxr rb- Ax r r r Ar 7.4 迭代法 有限元法基础 2.2.共轭梯度法共轭梯度法 定义：向量定义：向量 pi i 和 和 pj j , ,若若 称为关于称为关于A正交或共轭。正交或共轭。 使用共轭梯度方向搜索近似解使用共轭梯度方向搜索近似解xk 1 ， ，可加快收可加快收 敛速度。敛速度。 0() T ij ijp Ap 近似解为近似解为 共轭梯度方向共轭梯度方向 选择选择 求二次泛函极值求二次泛函极值 7.4 迭代法 有限元法基础 1kkkk =+ xxp 1 1 0 kkkk T kk prp p Ap k 7.4 迭代法 有限元法基础 l共轭梯度法迭代公式共轭梯度法迭代公式 1 1 1 11 1 kkkk kkkk T kk k T kk T kk k T kk kkkk =+ xxp prp pAr