关于线性方程组的数值解法一般有两类: 直接法和 迭代法
1. 直接法:直接法就是经过有限步算术运算,可求得线性方程组精确解的方法(若计算过程中没有舍入误差)。常用于求解低阶稠密矩阵方程组及某些大型稀疏矩阵方程组(如大型带状方程组)。
2. 迭代法:迭代法就是用某种极限过程去逐步逼近线性方程组精确解的方法。 优点:存储单元少、程序设计简单、原始系数矩阵在计算过程始终不变等; 缺点:存在收敛性及收敛速度问题。常用于求解大型稀疏矩阵方程组(尤其是由微分方程离散后得到的大型方程组)。
本文主要讲直接方法。
1 高斯消去法
高斯消去法也称为逐次消去法,是一个古老的求解线性方程组的方法,特别是由其改进而得到的选主元素消去法、三角分解法仍然是目前计算机上常用的有效方法。
高斯消去法的过程写成定理即为:
设,其中。
1. 如果,则可通过高斯消去法将约化为等价的三角形线性方程组(图1),且计算公式为:
(1)消元计算
(2)回代计算
2. 如果为非奇异矩阵,则可通过高斯消去法(及交换两行的初等变换)将方程组约化为方程组(图1)。
上述消元和回代过程总的乘除法次数为
高斯消去法写成MATLAB代码为:
function
但是,由算法可知,高斯消去法对于某些简单的矩阵可能会失败,例如
2 列主元消去法
由高斯消去法可知,在消元过程中可能出现
由此来看,对于一般矩阵来说,最好每一步选取系数矩阵(或消元后的低阶矩阵)中绝对值最大的元素作为主元素,以使得高斯消去法具有较好的数值稳定性。这就是全主元素消去法的思想,不足之处在于其在选主元时要花费较多机器时间,目前主要使用的是列主元消去法。
列主元消去法总体过程与高斯消去法基本一致,唯一不同之处在于在选主元之前增加了一步按列选主元,其过程直接写成MATLAB代码为:
function
3 矩阵三角分解法
高斯消去法有很多变形,有的是高斯消去法的改进、改写,有的是用于某一类特殊性质矩阵的高斯消去法的简化。
这里矩阵的三角分解一般用的是矩阵的LU分解,即将矩阵分解为一个下三角矩阵L与一个上三角矩阵U的乘积,当然前提条件是矩阵的顺序主子式不为
function
输入矩阵
可得其LU分解得到的矩阵为:
L
一旦实现了矩阵的
(1)
(2)
这里只展示不选主元的三角分解法。其计算公式为:
(1)
计算
(2)
(3)
求解
(4)
(5)
写成MATLAB代码则为:
function
除了以上几种方法,还有针对特定矩阵的平方根法(对称正定矩阵)、追赶法(三对角矩阵),其方法也都比较简单,呃……还是来讲一下吧。
4 平方根法
应用有限元法求解结构力学问题时,最后归结为求解线性方程组,系数矩阵大多具有对称正定性质。所谓平方根法,就是利用对称正定矩阵的三角分解而得到的求解对称正定方程组的一种有效方法,目前在计算机上广泛应用平方根法解此类方程组。
先给出两个三角分解定理
对称阵的三角分解定理 设为阶对称矩阵,且的所有顺序主子式均不为零,则可唯一分解为,其中为单位下三角矩阵,为对角矩阵。对称正定矩阵的三角分解或楚列斯基(Cholesky)分解 如果为阶对称正定矩阵,则存在一个实的非奇异下三角矩阵使,当限定的对角元素为正时,这种分解是唯一的。
从而可得平方根法得计算过程为:
对于
(1)
(2)
求解
(3)
(4)
其MATLAB代码为:
function
从平方根法计算过程可以看出,在计算
其计算过程为:
对于
(1)
(2)
(3)
求解
(4)
(5)
写成MATLAB代码为:
function
5 追赶法
在一些实际问题中,例如解常微分方程边值问题,解热传导方程以及船体数学放样中建立三次样条函数等,都会要求解系数矩阵为对角占优的三对角线方程组(图4)。
图4可简记为
我们利用矩阵的三角分解,得到
求解
从而得到求解三对角线方程组的追赶法公式:
(1)计算
(2)解
(3)解
我们将计算系数
其MATLAB代码为:
function
6 总结
下面把几种方法放在一起,来看看这几种方法的运行效率,这里给的系数矩阵是一个
clear
运行得到数值结果为:
Method
图像为:
从数值结果或图形上都可以很清楚地看到各种方法的效率高低,很直观地展现了追赶法对于解三对角矩阵的优势。
参考
- ^注意这里课本上把y错写成b了!