高斯消元法

求解线性方程组的直接解法5.1Gauss消去法①三角方程组先举一个简单的例子来说明消去法的基本思想.例1.用消去法解方程组(3).122(2),54(1),632132321xxxxxxxx解第一步.将方程(1)乘上-2加到方程(3)上去,消去(3)中的未知数1x,得到(4).11432xx第二步.将方程(2)加到方程(4)上去,消去方程(4)中的未知数2x,得到与原方程组等价的三角形方程组(5).62,54,6332321xxxxxx显然,方程组(5)是容易求解的,解为.)3,2,1(Tx上述过程相当于332331(-2)6-5620014011111-56140140111156122140111)|(rrrrrrbA其中用ir表示矩阵的第i行.下面我们讨论求解一般线性方程组的高斯消去法.一般地nnnnnnnnbxabxaxabxaxaxa2222211212111当a11a22…ann≠0时,可解出xn=bn/annfork=n-1:1xk=(b1-ak,k+1xk+1-…-aknxn)/akkend注:kkbx,可用同一组单元.并可解出一个未知数即代入其它方程消去该未知数Gauss消元法的流程图为:流程图中，,(,1,2,...,)ijiabijn分别为线性方程组的系数矩阵和常数向量；k是循环次数。②顺序消去法一般地,k=1对n阶方程组消去第k个元(akk≠0):nkknnknnnkkkkkknkkkknkknnknnknkkkkkknkkkkbbbaamaamaaabbbaaaaaaaaa11,1,11,1,11,11,,11,1,11,)()(这里各行变换:i行-k行×mik,其中mik=aik/akk,i=k+1,…,n而后,k=2对n-1阶方程组消第2个元…我们有如下顺序消元算法:fork=1:n-1ifakk≠0fori=k+1:nmik=aik/akki行=i行-k行×mikendelsestopendend（每行包括右端项!）可细化,也可存储mik于aik得:fork=1:n-1ifakk≠0fori=k+1:naik=aik/akkforj=k+1:naik=aik-aik×akjendbi=bi-aik×bkendelsestopendend顺序消元过程和回代过程连起来就可得精确解.顺序消元算法也可将系数矩阵和右端项分开:fork=1:n-1ifakk≠0fori=k+1:naik=aik/akkforj=k+1:naik=aik-aik×akjendendelsestopendend(注意mik在aik)fork=1:n-1fori=k+1:nbi=bi-aik×bkendendGauss消去法运算量消去第k个元素时，对矩阵作加法和乘法运算各(n-k)×(n-k)次,除法(n-k)次.对右端作加法和乘法运算各(n-k)次.分别共12+22+…+(n-1)2=n(n-1/2)(n-1)/3和1+2+…+(n-1)=n(n-1)/2次加法乘法,消元时还有1+2+…+(n-1)=n(n-1)/2次除法.另外回代过程中加法和乘法运算各n(n-1)/2次,除法n次.运算量主要是消元的贡献,加法和乘法运算各约n3/3.定理1.设Ax=b,其中A.Rnn(1)如果),,1,2,(k0)(nakkk则可通过Gauss消去法将Ax=b约化等价的三角形方程组)()2(2)1(121)()2(2)2(22)1(1)1(12)1(11nnmnnnnnbbbxxxaaaaaa,且计算公式为:(a)消元过程设0)(kkka,对1,,2,1nk计算nkkjibmbbamaaaamkkikkikikkjikkijkijkkkkikik,,2,1,/)()()1()()()1()()((b)回代过程1,2,,1/)(/1)()()()()(niaxabxabxnijiiijiijiiinnnnnn(2)如果A为非奇异矩阵,则可通过Gauss消去法(及交换两行的初等变换)将方程组Ax=b约化为)()2(2)1(121)()2(2)2(22)1(1)1(12)1(11nnmnnnnnbbbxxxaaaaaa.行列式和逆矩阵易知顺序消元过程中，行列式不变.因此det(A)=det(U)=u11u22…unn,这里U是顺序消元过程结束时的上三角矩阵.A的顺序主子式。也是同样情况,即有det(Ak)=det(Uk)=u11u22…ukk.。顺便指出,消第k个元素时，左上角的元素(称之为主元素)非零时才能进行.所有主元素非零等价于det(Ak)≠0,k=1,2,…,n-1.求A的逆矩阵相当于解n个方程组,它们的系数矩阵是A,n个右端项是n个单位向量.用增广矩阵表示即(A∣I)，其中I是单位矩阵.消元过程可对所有右端项同时进行,消元过程结束再分别回代.也可巧妙安排就地求逆.③主元素法如果akk=0顺序消元无法进行,如果akk很小，可以进行。但将引出很大误差.例2.用顺序消去法解方程组211110001.021xx解用精确运算:999819999010001.0211110001.0解得(10000/9999,9998/9999)≈(1.00010001000100,0.99989998999900)若在十进三位尾数舍入的浮点计算机系统中运算,第二行将是(0-10000∣-10000)从而得到解1,021xx,与实解相去甚远.究其原因,小主元作分母扩大误差所致.如果两个方程(两行)交换次序再消元:1210111210001.011211110001.0得到解1,121xx,与实解很近.下面引进的列主元素法在消元前，先选该列中绝对值最大的做主元(交换两行,每行包括右端项!):fork=1:n-1找p:nkkkkkpkaaaa,,,max,1p行k行ik=pifakk≠0fori=k+1:nmik=aik/akki行=i行-k行×mikendelsestopendendik=p记录每次选主元(两行交换)的指标.跟前面一样,算法还可细化,把矩阵和右端分开算(这时就要用指标ik)。列主元素法乘数mik绝对值不大于1,不会增加误差。列主元素法用来求行列式时，要注意两行交换行列式变号。还有一种全主元素法,它是在整个右下(n-k)×(n-k)矩阵找绝对值最大的做主元(交换行及列).这对误差控制有利,但搜索太费时.通常列主元素法误差控制就已可以了.以下是列主元素法的流程图:流程图中，,(,1,2,...,)ijiabijn分别为线性方程组的系数矩阵和常熟向量；k是循环次数；d是主元素；二实验部分本章实验内容：实验题目：Gauss消元法,追赶法,范数。实验内容：①编制用Gauss消元法求解线性方程组Ax=f的程序。②编制用追赶法求解线性方程组Ax=f的程序。③编制向量和矩阵的范数程序。实验目的：①了解Gauss消元法原理及实现条件,熟练掌握Gauss消元法解方程组的算法,并能计算行列式的值。②掌握追赶法，能利用追赶法求解线性方程组。③理解向量和矩阵范数定义,性质并掌握其计算方法.编程要求：利用Gauss消元法,追赶法解线性方程组。分析误差。计算算法:①Gauss消元法:1.消元过程设0)(kkka,对1,,2,1nk计算nkkjibmbbamaaaamkkikkikikkjikkijkijkkkkikik,,2,1,/)()()1()()()1()()(⒉回代过程1,2,,1/)(/1)()()()()(niaxabxabxnijiiijiijiiinnnnnn②追赶法：1.分解Ax=f:)/(,/1111iiiiiabcbc(1,,3,2ni)2.解Lg=f,求g:)/()(,/11111iiiiiiiabgafgbfg(ni,,3,2)3．解Ux=g,求x:1,iiiinnxyxyx（1,2,,2,1nni）③范数：常用向量范数有：（令x=(x1,x2,…,xn)T）1-范数：║x║1=│x1│+│x2│+…+│xn│2-范数：║x║2=(│x1│2+│x2│2+…+│xn│2)1/2∞-范数：║x║=max(│x1│,│x2│,…,│xn│)常用的三种向量范数导出的矩阵范数是：1-范数:║A║1=max{║Ax║1/║x║1=1}=nijjnja11max2-范数:║A║2=max{║Ax║2/║x║2=1}=1,λ1是ATA的最大特征值.∞-范数:║A║=max{║Ax║/║x║=1}=njijnia11max实验例题⑴：用Gauss消元法解方程组654131321112321xxx实验例题⑵:用追赶法解三对角方程组Ax=b,其中00001,2100012100012100012100012bA实验例题⑶：设3.01.05.06.0A，计算A的行列范数.程序①:Gauss消元法functionx=Gauss(A,b)%A是线性方程组的系数矩阵,b为自由项.n=length(A)fork=1:n-1m(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m(k+1:n,k)*A(k,k+1:n);b(k+1:n)=b(k+1:n)-m(k+1:n,k)*b(k);endx=zeros(n,1);x(n)=b(n)/A(n,n);fork=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);endx=x';disp(sprintf('kx(k)'));fori=0:ndisp(sprintf('%d%f',i,x(i+1)));end数值结果：x=Gauss(A,b)n=3kx(k)01.11111110.77777822.555556程序②:追赶法function[x,y,beta]=zhuiganfa(a,b,c,f)%a,b,c是三对角阵的对角线上的元素,f是自由项.n=length(b);beta(1)=c(1)/b(1);fori=2:nbeta(i)=c(i)/(b(i)-a(i)*beta(i-1));endy(1)=f(1)/b(1);fori=2:ny(i)=(f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1));endx(n)=y(n);fori=n-1:-1:1x(i)=y(i)-beta(i)*x(i+1);enddisp(sprintf('kx(k)y(k)beta(k)'));fori=0:ndisp(sprintf('%d%f',i,x(i+1),y(i+1),beta(i+1)));end数值结果:a=[0-1-1-1-1]';b=[22222]';c=[-1-1-1-10]';f=[10000]';[x,y,beta]=zhuiganfa(a,b,c,f)kx(k)y(k)beta