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

1、线性代数方程组数值解法及MATLAB实现妹述9* 20122090数廿学院12 it算机科学与技术1 8 （SH8本科）一、分林课题I®着科学技术的发展，提岀了大量复杂的数值廿算问题，在建立电子it算机成 为Uffiit算的主要工具以后，它以数字计算机求解数学间题的理论和方法为研究 对象。其数值廿算中线性代数方程的求解间题就广泛应用于各种工程技术方面。 因此在各种数据处理中，线11代数方程组的求解是最常见的问题之一。关干线性 代数方程组的数値解法一般分为两大类：直接法和迭代法。直接法就是经11有限步算术运算，可求的线性方程组猜确解的方法（若it算11 程没有舍人误差），但实际犹如舍人

2、误差的存在和影响，这种方法也只能求得近 似解，这类方法是解低阶囲密葩阵方程组级某些大里棉磴葩阵方程组的有效方 法。直接法色括高斯消元法，犯阵三角分解法、追赴法、平方根法。迭代法就是利用某种板限11程去逐步ilifiMtt方程组精确解的方法。迭代法具 有需要廿算机的存惆单元少，程序设廿简单，原始系数矩阵在it算11程始终不变 等优点，但存在收敛性级收敛速度冋題。迭代法是解夫塑稀磴矩阵方程组（尤其 是锁分方程离散后得到的大型方程组）的重要方法。选代法0 IS Jacobi法SOR法、 SSOR法等多种方法。二、研究课题-集性代数方程组数值解法一、直接法1、Gauss消元法通11一系列的加减消元运算

3、，也就是代数巾的加减消去法，以使A对角线以 下的元素化为零，将方程组化为上三角炭阵；然后，再逐一回代求解岀x向量。可修编.1.高斯消元法（加减消元）:首先将A化为上三角阵,再回代求解。步骤如下:第一步：第1行二+第布J = 2,丿,/2）第二步:第2行X越+筋行类似的做下去,我们有:% *心曙第 k 步：第k行x- +第布J = k + l,.,n。Clkkn-1步以后,我m可以得到变換后的犯阵为:如如如 知0迅？斎00窃 、000 a(n)注恿到，廿算il程中泌）处在被除的位置，因此整个廿算11程要保込它不为0。 所以，Gauss消元法的可行条件为："：'工0o就是要求A的

4、所有顺序主子式均不为0,即%5、det :工0=1,“I5丿因此，有些有解的间题，不能用Gauss消元求齡。另外，如果某个":很小的话，会引人大的炭差。例用Gauss消去法解方程组:(1）对増广矩阵进行初等变换12一33415-183-1-1-151111631-112 )“12-3O032700003715-521511292T 611735743612一3_3254743723474<121515T19T7盲丿一3_32515152926119103340372HTr 12一33/X'15、1)-183-1-1-151I116< 31-11 >H丿<

5、; 2丿得等价方程组12xJ - 3x2 + 3x3 + 4x4 = 1537 c 15_产+产+5兀=于1129“丁人3 +无=1 13 o91330 ft 棉兀4 = 0 ,兀3 = 3 x2 = 2 9 X = 1 o第_步:将(1)/3使人的系数化为1,再将、(3)式中人的系数都化为零，即由(2)-2 X(1f 得2 |Xl +§“2 +严=2 (1)由(3)-4 X得第二步：将(2尸除以2/3,使X2系数化为得吃+2皆。再将(3严式中&系数化为零，由(3)-(-14/3广巴得 器=-6第三步：将(3)除以18/3，使怡系数化为1,得2 1经消元后，得鸟如下三角代数方

6、程组:1+-+- = 2 (1) j2 + 2x5 = 0 忑二-10)1.2 0代过程由(3)得X3=1.«X3代入得Z,務X2、X3代入得x=l,所以,本題解为x=1,2,-1T1.3高斯消元的公式妹合以上讨论，不难看岀，高斯消元法解方程组的公式为第一步，消元(1 )令心, (i,j=1,2,3,n )b“,(i=1,2,3,-,n)(2)对 k=1 到 n-1,若宀0, »IIlik = ak<k) / akk<k, (i=k+1,k+2,n )y -lik*a*j<k), (i,j= k+1,k+2,n )b叫厲 (kk+1,k+2,n)第二步，H

7、l代若 ann<n)H 0Xn= bF/aJ"Xi= (bi(i) - s°m(ai严为)/-a? , (i = n-1,n-2,1 ) ,(j = i+1,i+2,n )2、LU分解法求解线11代数方程组除了高斯消元法外，还常用LU分解法(三角形分辭法)。LU分解法的优点是当方程组左常系数犯阵不变，仅仅是方程组右常列向量改变,即外加激励信号变化时,能够方便地求解方程组。矩阵的三角分解法可分为直接 三角分解法，列主元三轴分解法，平方根法，三对幷方程组的追tl法。下面过论直接 三角分解法。设n阶线性方程组Ax=b o假设能将方程组左常系数矩阵A,分解成两个三角阵 的乘枳

8、，K|1A=LU ,式中，L为主对轴线以上的元素均为零的下三轴矩阵，目主对 角线元素均为1的上三角矩眸；U为主对角线以下的元素血为零所以有,LUx=b 令 Ux=y,Ly=b由A=LU,由矩阵的乘法公式:3ii= Uij,j=12川aI1 = liiU1l,i=1,21,n推岀 Uij = aibj=1,2f,nlii = an/unf i=1,2,n这样就定岀了U的第一行元素和L的第一列元素。设已定出了 U的前kJ行和L的前k-1列，现在确定U的第k行和L的第k列。由葩阵乘当r>k时，1庐0,且I心1,因为nr-1n切=4灯_工5r-1j = k、k + 、nhk =(險-工厶吗&qu

9、ot;同理可推出廿算L的第k列的公式:因此棉到如下算法一杜利特(Doolittle )算法:(1)将矩阵分解为 A=LU,对心1,2,n； j=k,k+1,-n； i=k,k+1,n；解Ly=bJb-1公式2 : yk = bk 一工/好儿 k = 12,nr-l解Ux=y公式3：Itxk =(儿一 XX几)/%«=儿一1,1/上+1对大規模用聽问SL如果能舉通il调整方程及未知量的顺序使得方程组的系数矩阵成带状结构，则对系数矩阵使用通常的LU分解，可以保障单位下三角矩眸L及上三角矩阵U的为带状结飢3、直接三角分解法Gauss消去法还有许多变形，有些变形是为了利用特殊技巧减少误差，把

10、 Gauss消去法改写为更紧凑的形式，还有一些变形时根据某类矩阵的特性作一些 修正和简化，这些方法可统称为直接三角分解法。矩阵的三角分解设A的顺序主子式严0伙=1,2,，町，則可建立线牲方程组Ax=b的Gauss 消去法与矩萍分解的关系，即葩晖A的LU分解。这个间题前面已经讲的比较详 细了，此处不再贽述。Doolittle分解法首先假设A的顺序主子式乂都不为零，则A可作Doolittle分解，即A = LU, 其中厶是单位下三舟阵，有厶=1 , / v )时如=0 ; U是上三角阵，/ > j H Uy = 0 o仔细写岀为如(«11 他 2A =«21«2

11、2%=仏1.“22 4% CIJm爲1丿!.% >=LU ( 2.11 )在前面逐步推导LU的元素公式都要借助于有关的q来表示。现在強调 指出,只要从给定的A通ii比较(2.11 )式的两边就可能逐步地把LWU沟造岀来，而不必利用Gauss消去法的中同给果,逆种方法称为Gauss消去法的紧凑格 氏O根据矩阵的乘法規!M,比较(2.11)衣的两边可得5=工 *5，i、j = '2(2.12 )先算岀U的第1行，在(2.12)令心1,由/H=l, /iy=0(l<y</2)可得 接着在(2.12 )中令j = ,由创fs，从而算出E的第1列厶=ailuw =ai,aw &

12、#187; ，= 2,”若U的第1至k-lfj和厶的第1至k l列已经算岀，则由(2.12)可计算岀 U的第k行元素jfc-l1%=%一工5叫j = k，k +匕(2.13 )r-l接着再it算岀厶的第斤列元素Jt-lllk =,丿=1,/(2.14)%解方程组A.x = b就化为求解LUx = b,分两步求解。第一步解Ly = h,其中厶 为下三角阵，只要用逐次向前代人的方法；第二步解Ur = y,其中"为上三轴阵, 用逐次向后代入方法即可解除X。例用Doolittle方法求解:6 2 1 -1/ 24102-1114-1X35<-1 0-13,<-5>解：第1步

13、算出厶和(/： 13L= 11651 1<"6 To1937<6,U =210T12337To一3丄"1019174 )i可修编.第2步解Ly = b：23 191 yFy =第3步解Ux = >* fl:兀=（1,一1 丄一 1）7二、迭代袪迭代法具有的特点是速度快。与非线性方程的选代方法一样，需要我们构造一 个等价的方程，从而沟造一个收敛序列，序列的就是方程组的根。对方程组：Ax = b&等价变换：x = Gx + g如:令 A = M-N , 1Ax = b> (M - N)x = b=> Mx = h + Nx=> x =

14、 MlNx+Mlb«,我们可以梅造序列x=Gx+g若 f >a * => x* = G x*+g => Av* = b«,我们可以沟造序列严=GH + g若 x"' >a * => x* = G x*+g => Ax* = b同时:x(*+l)-x* = Gxik - Gx* = G(?ft)-x*)= . = Gfc+,(x(0)-x*)所以，序列收敛o/tO 与初值的选取无关1. Jacobi 迭代&內+.+气百=勺g+%©=»斗=二1（如心+气聲一勺）=>勺=2內 + "

15、23“ + + %© -b2） 。22兀=（色+ an -內-_仇） ann严二（如屮+%/）（如1）_一 1/ 丫丄 、"）丄 丄“ X “） A X2（°2內+。23“+ 401 一 2）< Cl22V （D _ 1 Q,丄 丄v "） U Xn（勺內 + + ®“g-乞）%格式很简单：屮讥二（壬佔+土勺屮T）aii >1jz迭代葩阵HA = D-L-U%0 'D =< oann丿'00、了0a2-Cln '_。2100L = u =0° an-n-Cln 一10><00易知,

16、Jacobi选代有（D-L-U）x = bDx = （L+U）x+bx = D（L + U）x + D'lbB0=D-|（L+（7） = Z-D-,A , /0 = Db , B°是简单迭代的迭代矩阵。看上述公式和过程很抽象,我们来举个简单例子:&序+坷2勺+坷3%3=勺< a2X + a22X2 + "23® = b2( aii 工 0 )4內+吗2兀2+吗3冯=“3西=丄(也一"2吃一。13小)得< x2 = (b2-a2lxA-a23x3)得上if变换的方程后，我们任取一向量/o)=(xr,xrr)r作初始近做，代人,

17、卍=(屮，閒北D)7，严=(輕,娜，習丁假设上述方程的准确解是"=(%勺43)册么如果limB=d,d我们就说上述迭代是收敛的。加f 3C2. Gauss - Seidel 迭代在Gauss - Seidel迭代中，使用最新it算出的分量值十一二(如屮+仏屮-勺)寸小=!2內3 +"23X3"' + +-“2)< a22(£+"_( (Jt+D(*+1) b、Xn 一 anX+ + 6 宀-i一 久)ann严“=二(土佔心+土伽屮7)Clii2+1易知,Gauss - Seidel迭代有(D-L-U)x = bDx = (L+U)

18、x+b x = D'(L + U)x + D'bGo=D'(L + U) = I-D"A , g = D"b , G()是简单迭代的迭代葩阵。二种方法都存在收敛tiiOo收敛条件迭代恪式收敛的充要条件是%的谱半g<1o有例子表明：Gauss-Seidel法收敛时,Jacobi法可能不收敛;而Jacobi法收敛时,Gauss-Seidel法也可能不收敛。1若为严格对角占优的话!M是收敛的,#0 1<20.020.03050.10.5、0.34丿假如方程组不满足收敛,有时候对其进行变换,可以改变收敛性。如求下述方程组的解:-180 A =-1

19、09,b =8,x(0) =(0,0,0/、9-1-b7/ 9-1-1>7、A =-180,b =7<-109丿6格式x：z = l/9(x+xf+7)<x；*+,)=l/8(x；*+,+7)彳z =1/9(彳z)+8)结果尤=(0.777& 0.9722,0.9753卩兀=(0.9942,0.9993,0.9994)7尤=(0.9999,0.9999,0.9999/x(4> = (1.0000,1.0000,1.0000)rmxi, mo1也满足收敛。A =-11-11-221丿2-22Jacobi, J?。的特征方程为-1A-1=0-22A解得& =

20、人=心=0 ,所以用Jacobi选代法收敛。2-2 2Gauss-Seidel, Go 的特征方程为一/I2-1 =0-22 -22 2解得人=0,/?2 =2+人=2 , /?(3() = _2_> 1,所以Gauss-Seidel选代法是不收敛的。3. 松弛选代记 Ar(=x(*+,)-xa>则严一尹+3可以看作在前一步上加一个修正量。若在修正量前乘以一个因子有严十+处严x(z=严+°(x(z 一严)对Gauss - Seidel选代招式严。=£>炯曲)+肋)+b)£如)=严+秋矿厶如)+矿也")+矿_严)?*+,) =(1-e)x

21、 + co(Dl Lxk+l) + D'Ux(k + D'lb)写成分量形式，有X"" =(1_0)召2 +0 (a”x/i + + "”£ 人 _bj 如x,a*" = (1 -e)x»"1 - (d，z"' + 仙不"'+ + %x,' -by) < °22X""-" =(l-<y)xu, +co (。”内" + + «_! a;,_|U, -bn) %关于直接法和迭代法的例題:1用选主元三

22、箱分解法求解下列方程组3兀+ 3x2 + 2x3 + 4x4 = 62x + 4x2 + 2x3 = 10<2x + 2x2 + 4x3 + 3x4 = 174%j + 2x2 + 2x3 + 6x4 = 8可修编.主 列 o 711 112% o%2268、3113-30613%Vi-x0丿48 6 H-3=y)7'33)4、/ 、"6 2420兀2102243兀317<42/>6>亠<8>42>68、/ 4/ 40、224224031017T、33>/416>0062244233422o> 3242u9 7o o

23、 1o 1 o1 %= 厶 nE再通il吩八得X的解。z 、 斗(24、X2-23<X4><-3>2用Jacobi选代法求解下列方程组，并II使精度为10。40.24-008)州'8、0.093-0.15X2=9以4, 3, 4分别除3方程两边0.04-0.084丿b丿、20)10.06-0.02)(0.031-0.05兀2=3,0.01-0.021 )( 、0-0.060.02/ 、召（2）即=-0.0300.05X2+30一0010.020 ?Z351+】）、X（0-0.060.02、八1有Jacobi迭代式甥=-0.0300.05严）2+3皆丿,-0.01

24、0.020 >3可证明此迭代榕式是收敛的，檢限值fill为解。'4A0.24-0.08''A0.06-0.02'0.0932-0.150.032-0.05,0.04-0.084久>,0.01-0.022 0.06兄2-00018A-0.020.05Z +0.0006A(2 + 0.06)(322 + 0.0 U - 0.0009)3(22-0.0018)解 II 入=一0.06,人=-0.01 + V0.0109 . _ 0.01 + 70109=由此可知，用Jacobi S代法是收敛的。=（2,3,5/为初始值，选代得，（0-0.060.02、1.

25、92、=-0.0300.053+33.19x；,) 3,-0.010.02° >込<5.04,（ 、 岸），0-0.060.02、5.92''1.909'?2)=-0.0300.053.19+3=3.194?2) 37,-0.010.020Z.5.04 /,5.045 丿e 0其实兀=-0.03< 3丿-0.01-0.060.02(1.909、，2、"1.909、00.053.194+3=3.1940.020 J5.045,<5.045,即在10=兀*,可以停止廿算了。所以兀即为近ifl解。分析：1H题題目可以直接说是解方程组

26、，然后求解过程我I可以先验算Jacobi 和Gauss - Seidel迭代的收敛怕况，再选择算法，一般来说Gauss - Seidel迭代要比Jacobi迭代快。三、用Matlab飲件求解线性方程组一.高斯消去法1噸序高斯消去法直接编写金令文件a=d=*n,n=size(a);c=n+1 a(：,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 01回带x(n)=a(n,c)/a(n,n);for g=n-1>1:1 x(g)=(a(g,c)-a(g,g+1:n

27、)*x(g+1:n)/a(g,g)end 2 列主高斯消去法"由于“,m=max(abs(a(k:n,k)”逆回的行是“k:n,k”的第几行,所以要通过修 正来把m改成真正的行的值。该程序只是演示程序，頁正机器it算不需要算主 元素所在列以下各折应为零的值。直接编写命令文件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,:);换 fienda(k+1:n, k:c)=a(k+1:n, k:c)-

28、(a(k+1 :n,k)/ a(k,k)*a(k, k:c);%消去end回带end x=0 0 0 01x(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二迭代法Jii代公式a=5 2 1;-1 4 2;2 -310初始向量选代的精度J迭代公式d=-12;20;3x=0;0;0;stop=1.0e-4L=-tril(a,-1)U=-triu(a,1)D=inv(diag(diag(a)X=D*(L+U)*x+D*d;n=1；while norm(X-x,inf)>=stopx=X;X=D*(L+U)*x+D*d;n=n+1;endG-S选代公式a=5 2 1;-1 4 2;2 -310x=0;0;0;d=-12;20;3stop=1.0e-4L=-tril(a-1)U=-triu(a,1)D=(diag(diag (a)X=inv(D-L)*U*x+inv(D-L)*d;n=1；while norm(X-x,inf)>= stopx=X;X=inv(D-L)*U*x+inv(D-L)*d;n=n+1;endXffl始向量G-S