线性代数方程组数值解法及MATLAB实现综述

1、线性代数方程组数值解法及MATLAB实现综述廖淑芳 20122090 数计学院 12计算机科学与技术1班（职教本科）一、分析课题随着科学技术的发展，提出了大量复杂的数值计算问题，在建立电子计算机成为数值计算的主要工具以后，它以数字计算机求解数学问题的理论和方法为研究对象。其数值计算中线性代数方程的求解问题就广泛应用于各种工程技术方面。因此在各种数据处理中，线性代数方程组的求解是最常见的问题之一。关于线性代数方程组的数值解法一般分为两大类：直接法和迭代法。直接法就是经过有限步算术运算，可求的线性方程组精确解的方法（若计算过程没有舍入误差），但实际犹如舍入误差的存在和影响，这种方法也只能求得近似解

2、，这类方法是解低阶稠密矩阵方程组级某些大型稀疏矩阵方程组的有效方法。直接法包括高斯消元法，矩阵三角分解法、追赶法、平方根法。迭代法就是利用某种极限过程去逐步逼近线性方程组精确解的方法。迭代法具有需要计算机的存储单元少，程序设计简单，原始系数矩阵在计算过程始终不变等优点，但存在收敛性级收敛速度问题。迭代法是解大型稀疏矩阵方程组（尤其是微分方程离散后得到的大型方程组）的重要方法。迭代法包括Jacobi法SOR法、SSOR法等多种方法。二、研究课题-线性代数方程组数值解法一、 直接法1、 Gauss消元法通过一系列的加减消元运算，也就是代数中的加减消去法，以使A对角线以下的元素化为零，将方程组化为上

3、三角矩阵；然后，再逐一回代求解出x向量。1.1消元过程1. 高斯消元法(加减消元)：首先将A化为上三角阵，再回代求解。 步骤如下：第一步： 第二步： 类似的做下去，我们有：第k步：。n1步以后，我们可以得到变换后的矩阵为：注意到，计算过程中处在被除的位置，因此整个计算过程要保证它不为0。所以，Gauss消元法的可行条件为：。就是要求A的所有顺序主子式均不为0，即因此，有些有解的问题，不能用Gauss消元求解。另外，如果某个很小的话，会引入大的误差。例 用Gauss消去法解方程组：（1）（1）对增广矩阵进行初等变换得等价方程组 回代得，。第一步：将(1)/3使x1的系数化为1,再将(2)、(3)

4、式中x1的系数都化为零，即由(2)-2×(1)(1)得 由(3)-4×(1)(1)得 第二步：将(2)(1)除以2/3，使x2系数化为1，得再将(3)(1)式中x2系数化为零，由(3)(1)-(-14/3)*(2)(2) ，得第三步：将(3)(2)除以18/3，使x3系数化为1，得经消元后，得到如下三角代数方程组：1.2 回代过程由(3)(3)得 x3=1,将x3代入(2)(2)得x2=-2,将x2 、x3代入(1)(1)得x2=1,所以，本题解为x=1,2,-1T 1.3高斯消元的公式综合以上讨论，不难看出，高斯消元法解方程组的公式为第一步，消元（1） 令aij(1) =

5、 aij , （i,j=1,2,3,n）bi(1) =bi ,（i=1,2,3,n）（2） 对k=1到n-1，若akk(k)0，进行lik = aik (k) / akk(k) ,（i=k+1,k+2,n）aij(k+1) = aij(k) - lik * akj(k), （i,j= k+1,k+2,n）bi(k+1) = bi(k) - lik * bk(k), （i= k+1,k+2,n）第二步，回代若ann(n) 0xn = bn(n) / ann(n)xi = (bi(i) sgm(aij(i) * xj )/- aii(i) ,（i = n-1,n-2,1）,( j = i+1,i+

6、2,n )2 、LU分解法求解线性代数方程组除了高斯消元法外，还常用LU分解法（三角形分解法）。LU分解法的优点是当方程组左端系数矩阵不变，仅仅是方程组右端列向量改变，即外加激励信号变化时，能够方便地求解方程组。矩阵的三角分解法可分为直接三角分解法,列主元三角分解法,平方根法,三对角方程组的追赶法。下面讨论直接三角分解法。设n阶线性方程组Ax=b 。假设能将方程组左端系数矩阵A，分解成两个三角阵的乘积，即A=LU ，式中，L为主对角线以上的元素均为零的下三角矩阵, 且主对角线元素均为1的上三角矩阵；U为主对角线以下的元素均为零所以有，LUx=b 令 Ux=y， 则 Ly=b由A=LU,由矩阵的

7、乘法公式： a1j = u1j , j=1,2,nai1 = li1u11 , i=1,2,n 推出 u1j = a1j, j=1,2,nli1 = ai1/u11, i=1,2,n 这样就定出了U的第一行元素和L的第一列元素。设已定出了U的前k-1行和L的前k-1列，现在确定U的第k行和L的第k列。由矩阵乘法：当r>k时，lkr=0, 且lkk=1，因为所以，同理可推出计算L的第k列的公式：因此得到如下算法杜利特（Doolittle）算法:(1) 将矩阵分解为A=LU,对k=1,2,n；j=k,k+1,n； i=k,k+1,n；公式1 (2) 解Ly=b (3) 解Ux=y对大规模稀疏

8、问题，如果能够通过调整方程及未知量的顺序使得方程组的系数矩阵成带状结构，则对系数矩阵使用通常的LU分解，可以保障单位下三角矩阵L及上三角矩阵U仍为带状结构.3、直接三角分解法Gauss消去法还有许多变形，有些变形是为了利用特殊技巧减少误差，把Gauss消去法改写为更紧凑的形式，还有一些变形时根据某类矩阵的特性作一些修正和简化，这些方法可统称为直接三角分解法。矩阵的三角分解设的顺序主子式，则可建立线性方程组的Gauss消去法与矩阵分解的关系，即矩阵的LU分解。这个问题前面已经讲的比较详细了，此处不再赘述。Doolittle分解法首先假设的顺序主子式都不为零，则可作Doolittle分解，即，其中

9、是单位下三角阵，有，时；是上三角阵，时。仔细写出为（2.11）在前面逐步推导和的元素公式都要借助于有关的来表示。现在强调指出，只要从给定的通过比较（2.11）式的两边就可能逐步地把和构造出来，而不必利用Gauss消去法的中间结果，这种方法称为Gauss消去法的紧凑格式。根据矩阵的乘法规则，比较（2.11）式的两边可得， （2.12）先算出的第1行，在（2.12）令，由，可得，接着在（2.12）中令，由，从而算出的第1列，若的第1至行和的第1至列已经算出，则由（2.12）可计算出的第行元素， （2.13）接着再计算出的第列元素 ， （2.14）解方程组就化为求解，分两步求解。第一步解，其中为下三

10、角阵，只要用逐次向前代入的方法；第二步解，其中为上三角阵，用逐次向后代入方法即可解除。例 用Doolittle方法求解：解：第1步算出和：，第2步解得：第3步解得：二、 迭代法迭代法具有的特点是速度快。与非线性方程的迭代方法一样，需要我们构造一个等价的方程，从而构造一个收敛序列，序列的极限值就是方程组的根。对方程组：做等价变换：如：令，则则，我们可以构造序列若则，我们可以构造序列若同时：所以，序列收敛 与初值的选取无关1. Jacobi迭代格式很简单：迭代矩阵记 易知，Jacobi迭代有，是简单迭代的迭代矩阵。看上述公式和过程很抽象，我们来举个简单例子： （）得上述变换的方程后，我们任取一向量

11、作初始近似，代入， 假设上述方程的准确解是 那么如果2. GaussSeidel迭代在GaussSeidel迭代中，使用最新计算出的分量值易知，GaussSeidel迭代有，是简单迭代的迭代矩阵。二种方法都存在收敛性问题。收敛条件迭代格式收敛的充要条件是的谱半径<1。有例子表明：Gauss-Seidel法收敛时，Jacobi法可能不收敛；而Jacobi法收敛时， Gauss-Seidel法也可能不收敛。若为严格对角占优的话则是收敛的，如假如方程组不满足收敛，有时候对其进行变换，可以改变收敛性。如求下述方程组的解：格式结果，也满足收敛。如 Jacobi，的特征方程为 解得，所以用Jacob

12、i迭代法收敛。 Gauss-Seidel，的特征方程为 解得，, 所以Gauss-Seidel迭代法是不收敛的。3. 松弛迭代记则可以看作在前一步上加一个修正量。若在修正量前乘以一个因子，有对GaussSeidel迭代格式写成分量形式，有关于直接法和迭代法的例题：1用选主元三角分解法求解下列方程组 列主这里再通过，得x的解。2用Jacobi迭代法求解下列方程组，并且使精度为。以4，3，4分别除3方程两边 即有Jacobi迭代式可证明此迭代格式是收敛的，极限值即为解。解得由此可知，用Jacobi迭代法是收敛的。取为初始值，迭代得，其实即在内，可以停止计算了。所以即为近似解。分析：此题题目可以直接

13、说是解方程组，然后求解过程我们可以先验算Jacobi和GaussSeidel迭代的收敛情况，再选择算法，一般来说GaussSeidel迭代要比Jacobi迭代快。三、用Matlab 软件求解线性方程组一。高斯消去法1.顺序高斯消去法直接编写命令文件a=d='n,n=size(a);c=n+1a(:,c)=d; for k=1:n-1a(k+1:n, k:c)=a(k+1:n, k:c)-(a(k+1:n,k)/ a(k,k)*a(k, k:c); 消去endx=0 0 0 0' 回带x(n)=a(n,c)/a(n,n);for g=n-1:-1:1x(g)=(a(g,c)-a(

14、g,g+1:n)*x(g+1:n)/a(g,g)end2.列主高斯消去法*由于“r,m=max(abs(a(k:n,k)”返回的行是“k:n,k”内的第几行，所以要通过修正来把m 改成真正的行的值。该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。直接编写命令文件a=d= 'n,n=size(a);c=n+1a(:,c)=d; (增广)for k=1:n-1r,m=max(abs(a(k:n,k); 选主m=m+k-1; (修正操作行的值) if(a(m,k)=0) if(m=k)a(k m,:)=a(m k,:); 换行enda(k+1:n, k:c)=a(k+

15、1:n, k:c)-(a(k+1:n,k)/ a(k,k)*a(k, k:c); %消去endendx=0 0 0 0' 回带x(n)=a(n,c)/a(n,n);for g=n-1:-1:1x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n)/a(g,g)end二 迭代法J迭代公式a=5 2 1;-1 4 2;2 -3 10d=-12;20;3x=0;0;0; 初始向量stop=1.0e-4 迭代的精度L=-tril(a,-1)U=-triu(a,1)D=inv(diag(diag(a)X=D*(L+U)*x+D*d; J迭代公式n=1;while norm(X-x,in

16、f)>=stop x=X;X=D*(L+U)*x+D*d;n=n+1;endXnG-S迭代公式a=5 2 1;-1 4 2;2 -3 10x=0;0;0; 初始向量d=-12;20;3 stop=1.0e-4L=-tril(a,-1)U=-triu(a,1)D=(diag(diag(a)X=inv(D-L)*U*x+inv(D-L)*d; G-S迭代公式n=1;while norm(X-x,inf)>= stop 时迭代中止否则继续x=X;X=inv(D-L)*U*x+inv(D-L)*d;n=n+1;endXnSOR迭代公式a=5 2 1;-1 4 2;2 -3 10x=0;0;

17、0; 初始向量d=-12;20;3stop=1.0e-4w=1.44; 松弛因子L=-tril(a,-1)U=-triu(a,1)D=(diag(diag(a)X=inv(D-w*L)*(1-w)*D+w*U)*x+w*inv(D-w*L)*d SOR迭代公式n=1;while norm(X-x,inf)>=stop 时迭代中止否则继续x=X;X=inv(D-w*L)*(1-w)*D+w*U)*x+w*inv(D-w*L)*d;n=n+1;endXn3. a*x=d a=5 2 1;-1 4 2;2 -3 10 d=-12;20;3(1)考察用J、G-S及SOR解此方程组的收敛性.(2)

18、分别用雅可比迭代法、高斯-赛德尔迭代法及逐次超松弛迭代法解此方程组。要求 时迭代中止，观察收敛速度。(1)    J迭代矩阵的谱半径：max(abs(eig(D*(L+U)= 0.506 G-S迭代矩阵的谱半径：max(abs(eig(inv(D-L)*U)= 0.2000SOR迭代矩阵的谱半径：max(abs(eig(inv(D-w*L)*(1-w)*D+w*U)=0.9756所以用J、G-S及SOR均收敛（且有 ）。(2)J： X =-4.0000 3.0000 2.0000 n =18G-S： X =-4.0000 3.0000 2.0000 n =8SOR( )： X =-4.0000 3.00