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

1、第七章 线性代数方程组的解法,7.1 高斯消去法及其变化 7.2 带状系数矩阵的直接法 7.3 利用外存的直接法 7.4 迭代解法,7. 线性代数方程组的解法,本章要点 Gauss消去法和三角分解法的原理和算法步骤 二维等带宽存储和一维变带宽存储的特点 分块解法的原理和实施方案 几种迭代解法,有限元法基础,7. 线性代数方程组的解法,关键概念 高斯循环消去法 三角分解法 二维等带宽存储 一维变带宽存储 分块解法 迭代解法 超松弛迭代法 梯度法 共轭梯度法 预条件共轭梯度法,有限元法基础,7. 线性代数方程组的解法,弹性力学的有限元方程为 对于弹性（本构关系线性）小变形（几何方程线性）问题K与q

2、无关，为常数矩阵，方程组为线性代数方程组。求解是有限元方程分析中费时最多的步骤。,有限元法基础,7. 线性代数方程组的解法,线性代数方程组的解法分为两大类，即直接解法和迭代法。 直接法的特点是，事先可按规定的算法步骤计算出它所需要的算术运算操作数，直接给出最后的结果。 迭代法的特点是，首先假定初始解，然后按一定的算法进行迭代，在每次的迭代过程检查解的误差，通过多次迭代直至满足解的精度要求。,有限元法基础,7. 线性代数方程组的解法,有限元法基础,直接解法 以高斯消去法为基础， 求解效率高，适用于小于 一定阶数的方程组，根据计算机和软件的不同有所不同，比如1万10万阶方程组。 迭代解法 当方程组

3、阶数过高 时，由于计算机有效位数的限制，直接解法中舍入 误差的积累影响精度， 采用迭代 解法。,7.1 高斯消去法及其变化形式,有限元法基础,一.高斯循序消去法 对于n阶线性方程组 1.消元,7.1 高斯消去法及其变化形式,有限元法基础,对于n阶线性方程组，共需进行n-1次消元： 第m次消元： 以第m-1次消元结果为基础 第m行元素为消元行， 为主元 仅对m+1 n行元素进 行元，并将m列元素中 m+1 n消为0,7.1 高斯消去法及其变化形式,有限元法基础,对i行m列（im）消元,将m列从m+1列的元素消为0,称为高斯消去因子,7.1 高斯消去法及其变化形式,有限元法基础,因此消元过程可以写

4、为 最终的 为上三角阵。其中,7.1 高斯消去法及其变化形式,有限元法基础,因此 因为K(0)为对称矩阵，所以,三角分解法的基础,7.1 高斯消去法及其变化形式,有限元法基础,特点 原系数矩阵是对称的，则每次消元后矩阵依然是对称的，只需存储一半的矩阵 消元结果中， 和 中的第i行就是（i-1)次消元的结果,7.1 高斯消去法及其变化形式,有限元法基础,载荷列阵 消元用到的元素都是矩阵 中的元素，因此， 的消元过程随时可进行，对于多载荷工况，可以利用消元后的 矩阵进行消元和回代求解。这样可大量节省求解所需的计算时间。 这是直接法相对迭代法的一个优点。,7.1 高斯消去法及其变化形式,有限元法基础

5、,2.回代求解 回代公式,7.1 高斯消去法及其变化形式,有限元法基础,例:用高斯消元法求下列矩阵的解,7.1 高斯消去法及其变化形式,有限元法基础,回代求解得：,7.1 高斯消去法及其变化形式,有限元法基础,二. 三角分解法 由高斯消去法能得到对 的三角分解 设,下三角阵,对角阵,上三角阵,7.1 高斯消去法及其变化形式,有限元法基础,由代数方程 可分解为 高斯消元法相当于 令,单位下三角阵,上三角阵,P在消元后的结果,7.1 高斯消去法及其变化形式,有限元法基础,三角分解后的代数方程求解步骤,7.1 高斯消去法及其变化形式,有限元法基础,三角分解的递推公式 K 中任意元素,7.1 高斯消去

6、法及其变化形式,有限元法基础,按行分解 i =1 i =2,7.1 高斯消去法及其变化形式,有限元法基础,i =3,4,n,7.1 高斯消去法及其变化形式,有限元法基础,按行分解存储情况,7.1 高斯消去法及其变化形式,有限元法基础,按列分解 j=1 j =2,3,n,7.1 高斯消去法及其变化形式,有限元法基础,按列分解存储情况,7.1 高斯消去法及其变化形式,有限元法基础,关于三角分解法 称为改进Choleski法，经典方法 比高斯消去法效率更高 只是改变了高斯消去法的循环循序和存储,按行三角分解 Do 15 i = 1, n Do 15 j = 1, n Do 15 m = 1, i-1

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

8、i的相关节点 如果节点j是节点i的相关节点 则 如果不是相关节点 则,7.2 带状系数矩阵的直接法,有限元法基础,一维变列高存储,主对角线位置M：1，2，4，6，10，12，16，18，22,j列上最上面非零元素行号,在一维存储中得位置,7.2 带状系数矩阵的直接法,有限元法基础,两种存储方式比较,变列高 找元素,7.2 带状系数矩阵的直接法,有限元法基础,二.二维等带宽的高斯消去法,工作三角形 由于系数矩阵呈带状 每次消元只涉及包括 主元在内的一个三角 形内的元素，称为工 作三角形 。,7.2 带状系数矩阵的直接法,有限元法基础,二维等带宽高斯消去法公式,7.2 带状系数矩阵的直接法,二维等

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

10、的元素不再参加消元，后面的元素尚未参加消元。 3）在整个消元过程中，工作区自上向下运动。 为分块解法奠定基础,7.3 利用计算机外存的直接法,有限元法基础,二.分块解法 设允许使用内存为NA,在每一分块，NQND行集成完毕，可进行消元修正，最后的ND行进入到下一块系数矩阵一起集成，消元修正。,7.3 利用计算机外存的直接法,有限元法基础,分块解法的特点 在每一分块中，系数矩阵的元素是先集成后消元修正 从求解的全过程看，系数矩阵的集成和消元修正是交替进行,7.3 利用计算机外存的直接法,有限元法基础,分块解法简单框图,7.3 利用计算机外存的直接法,有限元法基础,三.波前法（Front Meth

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

12、度矩阵K和载荷矩阵P不按自然编号进入内存，而是按计算时参加运算的顺序排列 2）在内存中保留尽可能少的一部分K和P中的元素 3）完成消元修正的行保存在外存,求解的自然编号是节点顺序 参加运算的顺序是单元顺序,7.3 利用计算机外存的直接法,有限元法基础,内存最大工作三角块（波前区） 三角形直角边为波前数W 最大波前数W ND 需要存储大量信息，以用于恢复完成集成，消元的自由度号,7.3 利用计算机外存的直接法,有限元法基础,波前法的计算过程 1）按单元顺序扫描计算单元刚度矩阵，并送入内存进行集成 2）检查那些DOF已完成集成，将集成完毕的DOF作为主元，对其他行、列进行消元修正 3）完成消元修正

13、后，将主DOF行有关的K和P中的元素移到外存 4）重复13步，将全部单元扫描完毕 5）按消元顺序，由后向前依次回代求解,7.3 利用计算机外存的直接法,有限元法基础,波前法一度是有限元研究者广泛采用的方法 与分块解法相比，波前法利用内存更少 由于频繁使用内外存交换求解效率低 编程复杂，以效率换取求解规模 随着计算机硬件的发展，目前已较少应用,7.4 迭代法,有限元法基础,一.雅克比迭代法 方程组为 方程组非奇异，且,7.4 迭代法,有限元法基础,方程组可改写为 雅克比迭代法 设初始解 迭代方程,7.4 迭代法,有限元法基础,为了便于编程，方程组可改写为 精度检查准则 为允许误差 当系数矩阵为严

14、格对角优势矩阵时，方法收敛,7.4 迭代法,有限元法基础,例1：用雅克比迭代法求解方程组 Ax=b,其中 相对误差控制为 初始解取为 精确解为 经过30次迭代可达到精度要求。,7.4 迭代法,有限元法基础,例2：用雅克比迭代法求解方程组 Ax=b,其中 相对误差控制为 初始解取为 精确解为 雅克比迭代法不收敛。原因：不满足严格对角优势。,7.4 迭代法,有限元法基础,二.高斯赛德尔(Gauss(Gauss-Seidel)迭代法 雅克比迭代式 在计算 时， 为已知量，将迭代式修改为,7.4 迭代法,有限元法基础,三.超松弛迭代法 逐次超松弛迭代是G-C迭代的一种加速收敛算法，引入松弛因子,G-C

15、迭代法 低松弛迭代法 超松弛迭代法,7.4 迭代法,有限元法基础,三.超松弛迭代法 逐次超松弛迭代是G-C迭代的一种加速收敛算法，引入松弛因子,G-C迭代法 低松弛迭代法 超松弛迭代法,7.4 迭代法,有限元法基础,例3：松弛因子迭代法求解例1 其余条件与例1 相同，考察松弛因子对收敛性的影响。,7.4 迭代法,有限元法基础,例4：松弛因子迭代法求解例2 其余条件与例2 相同，考察松弛因子对收敛性的影响。 可见，当系数矩阵不是严格的对角优势，算法也收敛。通常取 。,7.4 迭代法,有限元法基础,四.共轭梯度法 共轭梯度法由Hestenes和Stiefel提出，可以看作是称为最速下降法的梯度法发

16、展而来，简称CG法（Conjugate Gradient Method）。,7.4 迭代法,有限元法基础,1.梯度法 对与很多数学物理问题，如果方程是自伴随的，可等效为求解对应的二次泛函的极值问题。 线性方程组 对应的二次泛函 在xk处的梯度为 记负梯度方向,7.4 迭代法,有限元法基础,设 xk 为 的一个近似解，利用最速下降法改善近似值： 在 的方向移动xk ，即 再选择 使 取得极小值，即令,7.4 迭代法,有限元法基础,迭代公式 初始解 x0 可任意选取。 实际计算中发现收敛速度不高。,7.4 迭代法,有限元法基础,2.共轭梯度法 定义：向量 pi 和 pj ,若 称为关于A正交或共轭

17、。 使用共轭梯度方向搜索近似解xk1 ，可加快收敛速度。,近似解为 共轭梯度方向 选择 求二次泛函极值,7.4 迭代法,有限元法基础,7.4 迭代法,有限元法基础,共轭梯度法迭代公式,7.4 迭代法,有限元法基础,共轭梯度法特点 1）在迭代过程中，rk 是相互正交的 2）对n阶线性代数方程组通过n次迭代得到的rn一定是零 3）xn必定是原方程的解，由于有数值截断误差，计算时不一定得到精确解 另一加速算法为PCG（Preconditioned conjugate gradient method）收敛更快。,7.4 迭代法,有限元法基础,3.预条件共轭梯度法 (Preconditioned conjugate gradien