[数值算法与编程]高斯消去法

 在线性有限元计算中，通过单元刚度矩阵推导，整体刚度矩阵组装等步骤后，最终得到的是一个线性方程组：   

  其中，K是整体刚度矩阵，U是未知的节点位移向量，F是荷载向量。通常情况下，得到的该方程组中的K不是满秩矩阵，其矩阵行列式为0，因此此时U的解不唯一，在采用拉格朗日乘子法或者罚函数法等方式引入了边界条件进行修正后，K为满秩矩阵，从而U的解唯一。

  之后的任务就是求解这个方程组。求解该方程组的方式通常有两种：直接法和迭代法。高斯消去法，就是一种传统的直接求解法，其基本步骤如下：

例：求解  

（1）

1.通过行加减，将第2-4行的第一列变为0（即line2=line2+line1*4/5...，右侧的向量需要同样地对应操作）：

          

（2）

2. 继续行加减，将第3-4行的第二列变为0：

（3）

3. 继续行加减，将第4行的第3列变为0:

（4）

在得到（4）式以后，根据最后一行，可以直接求出U4的值，接着，再将U4的值带入第3行，求出U3，再代入第二行...，从而最终求出所有U。

由此可见，高斯消去法实际上是通过行变换，将系数矩阵转换为上三角矩阵，接着先求出最后一个值，然后依次往上回代，求出所有待求的U。

对上述过程进行算法总结：

（1）

 

上述过程完成的是对左侧的系数矩阵和右侧的向量进行变换，n是系数矩阵a的阶数，c是用于行变换的系数，通过上述变换后，得到的a为上三角矩阵，b为对应的右侧向量。

（2）

上述过程完成的是依次将求得的值回代入方程中，最终获得的b，即为方程ax=b的解。

 

编程：

 

上述例子是针对的满秩系数矩阵的求解，200阶的系数矩阵求解时间约为0.01s的级别。在实际的有限元计算中，形成的有限元刚度矩阵通常为具有对称性的带状稀疏矩阵（即非0元素主要集中在矩阵对角线附近），因此针对带状稀疏矩阵，还可进一步对高斯消去法进行优化，使得其求解效率更高。

 

 

 

参考文献：《有限元法（第二版）.下.理论、格式与求解方法》.Klaus-Jurgen Bathe著 轩建平 译

《Introduction to finite elements in engineering-Pearson (2012)》Belegundu, Ashok D._ Chandrupatla, Tirupathi R.

欢迎关注微信公众号 有限元术  分享有限元理论与实践知识