第3章解线性方程组的直接方法●●3.1高斯消去法设有线性方程组a11x1+a12x2+…+a1nxn=b1,a21x1+a22x2+…+a2nxn=b2,an1x1+an2x2+…+annxn=bn,(3.1)其矩阵形式为Ax=b,其中A=a11a12…a1na21a22…a2nan1an2…ann为非奇异矩阵，而x=x1x2xn，b=b1b2bn。首先举一个例子来说明消去法的基本思想。
例1用消去法解线性方程组7x1+8x2+11x3=-3,5x1+x2-3x3=-4,x1+2x2+3x3=1。(3.2)解第一步，消去线性方程组(3.2)中的第2、第3个方程中的x1。将第1个方程乘以-57加到第2个方程，第1个方程乘以-17加到第3个方程，得到等价的线性方程组7x1+8x2+11x3=-3,-337x2-767x3=-137,67x2+107x3=107。(3.3)第二步，消去方程组(3.3)第3个方程中的x2。将第2个方程乘以633加到第3个方程，得到等价的线性方程组7x1+8x2+11x3=-3,-337x2-767x3=-137,-611x3=1211。(3.4)这是一个上三角线性方程组，这种线性方程组非常容易求解。从最后一个方程解出x3，然后将x3的值代入第2个方程，解出x2，最后依次将x2、x3的值代入第1个方程，从中解出x1，从而得到线性方程组(3.4)的解为x1=-3，x2=5，x3=-2。这种先把线性方程组化为同解的上三角线性方程组(称为消元过程)，然后按由下到上的相反顺序求解上三角线性方程组(称为回代过程)，得到原线性方程组解的方法，称为高斯消去法。
由线性代数知识，上述过程就是用行的初等变换将原线性方程组的增广矩阵(Ab)化为行阶梯形矩阵，然后再回代求解，从而将求解原线性方程组(3.1)的问题转化为求解简单线性方程组(3.4)的问题。
下面我们来讨论一般的解n阶线性方程组的消去法的问题。
为了符号统一，令(Ab)def(A(1)b(1))。当a(1)11≠0时，第一步消元后得(A(2)b(2))=a(1)11a(1)12…a(1)1nb(1)1a(2)22…a(2)2nb(2)2a(2)n2…a(2)nnb(2)n，(3.5)简记为A(2)x=b(2).此步骤所需计算的是下面的量：
(1) 行乘数li1=a(1)i1a(1)11,i=2,3,…,n；
(2) a(2)ij=a(1)ij-li1a(1)1j,i,j=2,3,…,n,b(2)i=b(1)i-li1b(1)1,i=2,3,…,n。
一般地，假定已完成了k-1步消元，得到(A(k)b(k))=a(1)11a(1)12………a(1)1nb(1)1a(2)22………a(2)2nb(2)2a(k)kk…a(k)knb(k)ka(k)nk…a(k)nnb(k)n。若a(k)kk≠0，则可以进行第k步消元，从而得到(A(k+1)b(k+1))，简记为A(k+1)x=b(k+1)。此步骤需要计算的是下面的量：
(1) 行乘数lik=a(k)ika(k)kk,i=k+1,k+2,…,n； (3.6)
(2) a(k+1)ij=a(k)ij-lika(k)kj,i,j=k+1,k+2,…,n,b(k+1)i=b(k)i-likb(k)k,i=k+1,k+2,…,n。(3.7)
这样完成n-1步消元后，得到(A(n)b(n))=a(1)11a(1)12…a(1)1nb(1)1a(2)22…a(2)2nb(2)2a(n)nnb(n)n。或记为A(n)x=b(n)。(3.8)因为A为非奇异矩阵，故a(n)nn≠0，用回代过程求解上述三角形线性方程组，得 xn=b(n)na(n)nn,xi=b(i)i-∑nj=i+1a(i)ijxja(i)ii,i=n-1,n-2,…,1。如果a(1)11=0，由于A为非奇异矩阵，所以A的第1列一定有元素不等于零，例如a(1)k1≠0，于是可交换两行元素，将a(1)k1对调到第1行的位置，然后进行消元计算，这时A(2)的右下角矩阵(n-1阶)也为非奇异矩阵。继续这一过程，高斯消去法照样可以进行计算。
总结上述讨论即有下面的定理。
定理1如果A为n阶非奇异矩阵，则可通过高斯消去法(及交换两行的初等变换)将线性方程组(3.1)化为上三角线性方程组(3.8)。
我们知道，只有当a(k)kk≠0(k=1,2,3,…,n)时，才能顺利进行高斯消去法，那么矩阵A在什么条件下才能保证a(k)kk≠0(k=1,2,…,n)呢？下面的引理给出了这个问题的回答。
引理主元素a(k)kk≠0(k=1,2,…,l)的充分必要条件是矩阵A的顺序主子式Dk≠0(k=1,2,…,l)。
推论如果n阶矩阵A的前n-1个顺序主子式Dk≠0(k=1,2,…,n-1)，则a(1)11=D1；a(k)kk=DkDk-1(k=2,3,…,n)。
定理2如果n阶矩阵A的所有顺序主子式Dk≠0(k=1,2,…,n)，则可通过高斯消去法(不进行交换两行的初等变换)将线性方程组(3.1)化为三角形线性方程组(3.8)。
例2用高斯消去法求解线性方程组12x1-3x2+3x3=15,18x1-3x2+x3=15,-x1+2x2+x3=6。解利用高斯消去法容易计算12-331518-3115-1216→12-3315032-72-15207454294→12-3315032-72-1520016316。于是得到与原线性方程组同解的上三角线性方程组12x1-3x2+3x3=15,32x2-72x3=-152,163x3=16。通过回代得到原线性方程组的解为x1=1,x2=2,x3=3。现在估计高斯消去法的计算量。
由求解过程可知，第k步消元计算需做(n-k)2+2(n-k)次乘除法和(n-k)2+(n-k)次加减法，所以完成n-1步消元后，总运算量是乘除法次数=∑n-1k=1(n-k)(n-k+2)=n33+n22-56n，加减法次数=∑n-1k=1(n-k)(n-k+1)=n33-n3；回代求解所需的总运算量是乘除法次数=∑n-1i=1(n-i+1)=n(n+1)2，加减法次数=∑n-1i=1(n-i)=n(n-1)2。因此，求解线性方程组(3.1)的高斯消去法所需的总乘除法及加减法的次数分别为乘除法次数=n33+n2-n3≈n33(当n较大时)；加减法次数=n33+n22-56n≈n33(当n较大时)。由此可见，高斯消去法求解线性方程组比用克莱姆法则求解线性方程组的计算量少很多，从而优越得多。
实现顺序高斯消去法的MATLAB函数文件 gauss.m 如下。function x=gauss(a,b,flag)
%线性方程组求解的顺序高斯消去法,a为系数矩阵,b为常向量
%flag若为0,则显示中间过程,否则不显示,x为解
\[na,m\]=size(a);n=length(b);if na~=m,error('系数矩阵必须是方阵');return;end
if n~=m,error('系数矩阵a的列数必须等于b的行数');return;end
if nargin<3,flag=0; end;
a=\[a,b\];
for k=1:(n-1)
a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1))-a((k+1):n,k)/a(k,k)a(k,(k+1):
(n+1));
a((k+1):n,k)=zeros(n-k,1);
if flag==0,a,end
end
x=zeros(n,1);x(n)=a(n,n+1)/a(n,n);
for k=n-1:-1:1
x(k,:)=(a(k,n+1)-a(k,(k+1):n)x((k+1):n))/a(k,k);
end例3在MATLAB命令窗口求解例1。
解输入format long;a=\[7 8 11;5 1 -3;1 2 3\];b=\[-3; -4; 1\];
x=gauss(a,b)●3.2高斯列主元消去法3.2高斯列主元消去法
在高斯消去法中，a(k)kk(k=1,2,…,n-1)称为主元。消元过程中若出现a(k)kk=0，则消元法不能顺利进行下去。即使a(k)kk≠0，但绝对值很小，用它作为除数会使其他元素的数量级急剧增大，因而舍入误差也随之增大，致使计算结果不可靠（见下述例4）。
实用的高斯消去法应该在消元过程的每一步，都在可能范围内选择绝对值较大的元素作为主元，以免舍入误差的增长，这就是高斯主元消去法。
常用的一种主元消去法是按列选主元，即对1≤k≤n-1，在消元过程的第k步时，选ck=max k≤i≤n|a(k)ik|作为主元，并将其所在的行与第k行对换，再进行第k步消元，这样便保证每个行乘数的绝对值都不大于1。
高斯列主元消去法的算法 
1 消元过程： 消元结果冲掉A，行乘数冲掉aik(i>k), det存放行列式。
(1) 1det。
(2) 对于k=1,2,…,n-1, 做①～⑤： 
① 按列选主元ck, 即确定ik使ck=|aikk|=maxk≤i≤n|aik|。② 若ck=0，则输出错误信息；停机。
③ 若ik=k，转④，否则执行换行akjaikj,j=k,k+1,…,n，bkbik；-det det。④ 消元计算： lik=aikakkaik,i=k+1,…,n,aij-likakjaij,i,j=k+1,…,n,bi-likbkbi,i=k+1,…,n。⑤ akk·detdet。
(3) 若ann=0，则判断bn是否等于零： 
① bn≠0，则输出“无解”信息；停机。
② bn=0，则输出“非唯一解”信息，停机。
若ann≠0，则ann·detdet。
2 回代过程： 解冲掉b。xn=bnannbn，xi=bi-∑nj=i+1aijbjaiibi,i=n-1,…,2,1。3 输出x1,x2,…,xn及det。det值的大小可以作为此线性方程组是否病态的参考。
高斯列主元消去法在计算机上编程实现，其框图见图3.1。
图3.1高斯列主元消去法求解线性方程组
例4求解线性方程组0.0003x1+3.0000x2=2.0001,1.0000x1+1.0000x2=1.0000（准确解为x1=1/3,x2=2/3）。
解用四舍五入的五位有效数字浮点数计算。
(1) 不选主元的高斯消去法0.00033.00002.00011.00001.00001.0000→0.00033.00002.00013333.3－9999.0－6666.0,其中右端矩阵左下角虚线下的数是行乘数l21=1.00000.0003=3333.3。
回代求解，得x2=23≈0.66667,x1=2.0001－3.0000×(2/3)0.0003≈0.30000。我们看到，x2仅有0.00000333…的误差，但在求x1时，由“小主元”作除数得到的近似值只有一位有效数字，绝对误差为0.0333…, 相对误差达11.11%，计算结果的有效数字严重丢失，这组解是无效的。
(2) 选列主元的高斯消去法0.00033.00002.00011.00001.00001.0000→1.00001.00001.00000.00033.00002.0001→1.00001.00001.000002.99971.9998，回代求解，得x2=23≈0.66667,x1=1.0000－(2/3)1.0000≈0.33333这组解的五位数字都是有效数字。
选列主元的高斯消去法与不选主元的高斯消去法有什么关系呢？当线性方程组的系数矩阵满足一定的条件时，可以得到一些结论。为此下面先给出（按行）严格对角占优的概念，再给出相应的结论。
定义1(对角占优阵)设A=(aij)n×n，若|aii|>∑nj=1,j≠i|aij|(i=1,2,…,n)，称A为（按行）严格对角占优阵(强对角占优阵)。
如矩阵 －8433－5－11－46, 由于 |－8|>4+3,|－5|>3+|－1|,6>1+|－4|,因此该矩阵为（按行）严格对角占优矩阵。
定理3如果一个线性方程组的系数矩阵对称且（按行）严格对角占优，则不选主元的高斯消去法即选列主元的高斯消去法。 
当线性方程组的系数矩阵元素的数量级相差很大时，即使用高斯列主元消去法，求解过程中有效数字也会严重损失。为避免这一现象发生，可先将系数矩阵元素的数量级大体均衡一下，例如用每行元素的绝对值最大模的数除该行各元素，称为行标度化，然后再施行高斯列主元消去法。但要指出的是，标度化的目的是找到真正的主元，消元过程仍是对原线性方程组进行的。
实现高斯列主元消去法的MATLAB函数文件lgauss.m如下。function x=lgauss(a,b,flag)
%线性方程组求解的高斯列主元消去法,a为系数矩阵,b为常向量
%flag若为0,则显示中间过程,否则不显示,x为解
\[na,m\]=size(a);n=length(b);
if na~=m,error('系数矩阵必须是方阵');return;end
if n~=m,error('系数矩阵a的列数必须等于b的行数');return;end
if nargin<3,flag=0; end;
a=\[a,b\];
for k=1:(n-1)
%选主元
\[ap,p\]=max(abs(a(k:n,k)));p=p+k-1;
if p>k,t=a(k,:);a(k,:)=a(p,:);a(p,:)=t;end
%消元
a((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1))-a((k+1):n,k)/a(k,k)a(k,(k+1):(n+1));
a((k+1):n,k)=zeros(n-k,1);
if flag==0,a,end,end
%回代
x=zeros(n,1);x(n)=a(n,n+1)/a(n,n);
for k=n-1:-1:1
x(k,:)=(a(k,n+1)-a(k,(k+1):n)x((k+1):n))/a(k,k);
end例5在MATLAB命令窗口求解例4。
解输入format long;a=\[0.0003 3.0000;1.0000 1.0000\];b=\[2.0001;1.0000\];
x=lgauss(a,b)●3.3矩阵分解在解线性方程组中的应用3.3矩阵分解在解线性方程组中的应用
本节介绍矩阵的杜利特尔(Doolittle)分解、楚列斯基(Cholesky)分解及其在解线性方程组中的应用，解三对角线性方程组的追赶法。
1 杜利特尔分解
从上述两节我们知道，消去法是利用矩阵的行初等变换把原线性方程组化为一个等价的三角形线性方程组进行求解的方法。例如，对于线性方程组x1+2x2+x3=0,2x1+2x2+3x3=3,-x1-3x2=2,(3.9)求解的过程是，对它的增广矩阵施行行初等变换12102233-1-302→12100-2130-112→12100-213001212，得到与线性方程组(3.9)等价的三角形线性方程组x1+2x2+x3=0,-2x2+x3=3,12x3=12。(3.10)然后经回代求得解，由线性方程组(3.10)的最后一个方程得出x3=1，将x3的值代入第2个方程解得x2=-1，最后将x3和x2的值代入第1个方程得到x1=1。
上述对增广矩阵实施行初等变换，实质上是先用初等矩阵L1，再用初等矩阵L2左乘该增广矩阵，其中L1=100-210101,L2=1000100-121。若记线性方程组(3.9)的系数矩阵为A，线性方程组(3.10)的系数矩阵为R，则容易验证L2L1A=R，或者说A=(L2L1)-1R=L-11L-12R，且有L-11=100210-101,L-12=1000100121,L-11L-12=100210-1121defL。于是A=LR。(3.11)主对角线元都是1的下(上)三角矩阵称为单位下(上)三角矩阵。式(3.11)表明，方阵A可以分解成一个单位下三角矩阵与一个上三角矩阵的乘积。
定义2如果方阵A可写成式(3.11)的形式，则称A可三角分解(或LU分解或LR分解)，而式(3.11)称为A的一个三角分解或杜利特尔分解。
我们指出，当A为可逆矩阵时，R也是可逆矩阵，其主对角线元都不为零，因而R可进一步化为对角矩阵D与单位上三角矩阵R1的乘积，即A=LDR1。(3.12)其中L,R1分别为单位下三角矩阵和单位上三角矩阵，D是对角矩阵，形如式(3.12)的分解为A的LDR分解。
关于矩阵的杜利特尔分解，我们不加证明地给出下述结论。
定理4n阶方阵A有唯一的杜利特尔分解的充分必要条件是A的前n-1个顺序主子式Dk≠0(1≤k≤n-1)。
如果已知方阵A的杜利特尔分解A=LR，那么解线性方程组Ax=b(3.13)相当于解两个三角形线性方程组Ly=b,Rx=y。 (3.14)当A可逆时，由式(3.14)的第1个线性方程组依次递推地解出y1,y2,…,yn，这个过程称为前