列主元消去法解方程组实验报告

实验名称：列主元消去法解方程组1引言我们知道，高斯消去法是一个古老的解线性方程组的方法。而在用高斯消去法解Ax=b时，其中设A为非奇异矩阵，可能出现()0kkka的情况，这时必须进行带行交换的高斯消去法。但在实际计算中即使()0kkka但其绝对值很小时，用()kkka作除数，会导致中间结果矩阵()kA元素数量级严重增长和舍入误差的扩散，使得最后的结果不可靠。因此，小主元可能导致计算的失败，我们应该避免采用绝对值很小的主元素。为此，我们在高斯消去法的每一步应该在系数矩阵或消元后的低阶矩阵中选取绝对值最大的元素作为主元素，保持乘数1ikm，以便减少计算过程中舍入误差对计算解的影响。一种方式是完全主元消去法，这种消去法是在每次选主元时，选择()()max0kkkkijijkinkjnaa为主元素。这种方法是解低阶稠密矩阵方程组的有效方法，但这种方法在选取主元时要花费一定的计算机时间。实际计算中我们常采用部分选主元的的消去法。列主元消去法即在每次选主元时，仅依次按列选取绝对值最大的元素作为主元素，且仅交换两行，再进行消元计算。2实验目的和要求运用matlab编写一个.m文件，要求用列主元消去法求解方程组（实现PA=LU）：12345671111111721111118321111110432111113543211117654321122765432128xxxxxxx要求输出以下内容：（1）计算解x；（2）L,U；（3）整形数组IP（i）（i=1,2，…，n-1）（记录主行信息）3算法原理与流程图（1）算法原理设有线性方程组Ax=b，其中设A为非奇异矩阵。方程组的增广矩阵为111211212222112[,]nnnnlnnnaaabaaabaaaAabb第1步（k=1）：首先在A的第一列中选取绝对值最大的元素1ia，作为第一步的主元素:1,111max0iiinaa，然后交换（A，b）的第1行与第i1行元素，再进行消元计算。设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组()()kkAxb(1)(1)(1)(1)(1)1112111(2)(2)(2)(2)22222()()()()()()()1,1()1,()()()[,][,]kkkknknkkkkkkkkknkkkknkkknnnnkaaaabaaabaabaabaabAbAb第k步计算如下：对于k=1，2，…，n-1（1）按列选主元：即确定ik使,max0ikkikkinaa（2）如果,0ikka，则A为非奇异矩阵，停止计算。（3）如果ik≠k，则交换[A，b]第ik行与第k行元素。（4）消元计算,(1,,)ikikikkkaamikna,(,1,,),(1,,)ijijikkjiiikkaamaijknbbmbikn消元乘数ikm满足:1,(1,,)ikikkkamikna（5）回代求解1/(0)()/,1,,2,1n当nnnnnniiijjiijibbaabbabain计算解1(,,)nxx在常数项b(n)内得到。（2）流程图见图1（3）4程序代码及注释%列主元消去法解方程组Ax=b，实现PA=LUfunction[x,L,U,IP,P]=gauss(A,b)%x为方程组的解，IP用来记录行信息%每次选列主元时，将A的第k行与第IP(k)行进行交换n=length(b);[p,q]=size(A);%当输入的系数矩阵不为方阵，或方阵维数与b不符时，报错ifp~=q||p~=nfprintf('Error!Pleaseinputagain!');end%为提高运行速度，给IP，P，L，U赋初值IP=zeros(1,n-1);L=zeros(n,n);U=zeros(n,n);P=eye(n);x=zeros(1,n);det=1.0;%按列选主元，并进行行交换，记录行信息fork=1:n-1IP(k)=k;form=k+1:nifabs(A(m,k))abs(A(k,k))IP(k)=m;endendI=eye(n);ifIP(k)~=kfori=1:np(i)=I(k,i);I(k,i)=I(IP(k),i);I(IP(k),i)=p(i);endA=I*A;b=I*b';b=b';end%进行消元计算fori=k+1:nA(i,k)=A(i,k)/A(k,k);b(i)=b(i)-A(i,k)*b(k);forj=k+1:nA(i,j)=A(i,j)-A(i,k)*A(k,j);enddet=det*A(k,k);P=I*P;end%回代求解x(n)=b(n)/A(n,n);fori=n-1:-1:1sum=0.0;forj=i+1:nsum=sum+A(i,j)*x(j);endx(i)=(b(i)-sum)/A(i,i);enddet=det*A(n,n);ifdet==0fprintf('Theequationshavenouniquesolution!');end%输出PA=LU中的L,U的信息fori=1:nforj=1:nifijU(i,j)=A(i,j);elseifi==jL(i,j)=1;U(i,j)=A(i,j);elseL(i,j)=A(i,j);endendendend附图1输入n，A，b，εdet=1.0k=1,2，…，n-1按列选主元,max0ikkikkinaac0=,ikka|c0|εik=k换行,(,,),detdetkjikjkikaajknbb消元计算iki(k+1,,）,a/*(k+1,,）,b*ikkkijijikkjiikkinaaaamajnbmbdetdetkka输出det()0A停机回代求解1/(当0)()/,1,2,det*detnnnnnnniiijjiijinnbbaabbabaina输出计算解及行列式值b（i）（i=1,2，…，n）及det停机YYNN5算例分析1、测试示例A=[123456];b=[3711];[xLUIPP]=gauss(A,b)Error!Pleaseinputagain!A=[123000456];b=[3711];[xLUIPP]=gauss(A,b)Theequationshavenouniquesolution!x=NaN-InfInfL=1.0000000.25001.00000001.0000U=4.00005.00006.000000.75001.5000000IP=33P=1000010102、计算过程（1）首先输入系数矩阵A和矩阵bA=[1111111211111132111114321111543211165432117654321];b=[781013172228];（2）输出结果[xLUIPP]=gauss(A,b)x=1.00001.00001.00001.00001.00001.00001.0000L=1.00000000000.28571.0000000000.42860.80001.000000000.57140.60000.75001.00000000.71430.40000.50000.66671.0000000.1429-0.2000-0.2500-0.3333-0.50001.000000.85710.20000.25000.33330.5000-1.00001.0000U=7.00006.00005.00004.00003.00002.00001.00000-0.7143-0.4286-0.14290.14290.42860.714300-0.8000-0.6000-0.4000-0.20000.0000000-0.7500-0.5000-0.25000.00000000-0.6667-0.3333-0.0000000000.50001.00000000001.0000IP=723457P=10000000100000001000000010000000100000000100000106讨论与结论1、时间复杂度：tic;[xLUIPP]=gauss(A,b);tocElapsedtimeis0.000856seconds.2、程序优化初次编程时，没有考虑到给一个变量赋初值的情况。虽然在MATLAB中变量不赋初值是完全允许的，但是由于一个变量中含有多个元素时，每次改变该数组的长度，便会增加计算机时间。另外，给程序加上一定的判断条件及报错信息，一定程度上有程序优化的作用。因此，本程序中的以下程序段都起到了程序优化的作用。%当输入的系数矩阵不为方阵，或方阵维数与b不符时，报错ifp~=q||p~=nfprintf('Error!Pleaseinputagain!');end%为提高运行速度，给IP，P，L，U赋初值IP=zeros(1,n-1);L=zeros(n,n);U=zeros(n,n);P=eye(n);x=zeros(1,n);det=1.0;参考文献[1]易大义，沈云宝，李有法.计算方法(第2版)，浙江大学出版社.p.29-53.[2]张琨高思超毕靖编著MATLAB2010从入门到精通电子工业出版社