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

1、精选优质文档-倾情为你奉上实验名称： 列主元消去法解方程组 1 引言我们知道，高斯消去法是一个古老的解线性方程组的方法。而在用高斯消去法解Ax=b时，其中设A为非奇异矩阵，可能出现的情况，这时必须进行带行交换的高斯消去法。但在实际计算中即使但其绝对值很小时，用作除数，会导致中间结果矩阵元素数量级严重增长和舍入误差的扩散，使得最后的结果不可靠。因此，小主元可能导致计算的失败，我们应该避免采用绝对值很小的主元素。为此，我们在高斯消去法的每一步应该在系数矩阵或消元后的低阶矩阵中选取绝对值最大的元素作为主元素，保持乘数，以便减少计算过程中舍入误差对计算解的影响。一种方式是完全主元消去法，这种消去法是在

2、每次选主元时，选择为主元素。这种方法是解低阶稠密矩阵方程组的有效方法，但这种方法在选取主元时要花费一定的计算机时间。实际计算中我们常采用部分选主元的的消去法。列主元消去法即在每次选主元时，仅依次按列选取绝对值最大的元素作为主元素，且仅交换两行，再进行消元计算。2 实验目的和要求运用matlab编写一个.m文件，要求用列主元消去法求解方程组（实现PA=LU）：要求输出以下内容：（1） 计算解x；（2） L,U；（3） 整形数组IP（i）（i=1,2，n-1）（记录主行信息）3 算法原理与流程图（1） 算法原理设有线性方程组Ax=b，其中设A为非奇异矩阵。方程组的增广矩阵为第1步（k=1）：首先在

3、A的第一列中选取绝对值最大的元素，作为第一步的主元素: ，然后交换（A，b）的第1行与第i1行元素，再进行消元计算。 设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组第k步计算如下：对于k=1，2，n-1（1）按列选主元：即确定ik使 （2）如果，则A为非奇异矩阵，停止计算。（3）如果ikk，则交换A，b第ik行与第k行元素。（4）消元计算消元乘数满足: （5）回代求解计算解在常数项b(n)内得到。（2） 流程图见图1（3）4 程序代码及注释%列主元消去法解方程组Ax=b，实现PA=LUfunction x,L,U,IP,P =gauss(A

4、,b)%x为方程组的解，IP用来记录行信息%每次选列主元时，将A的第k行与第IP(k)行进行交换n=length(b);p,q=size(A);%当输入的系数矩阵不为方阵，或方阵维数与b不符时，报错if p=q|p=n fprintf('Error! Please input again!');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;%按列选主元，并进行行交换，记录行信息for k=1:n-1 IP(k)=k; for m=k+

5、1:n if abs(A(m,k)>abs(A(k,k) IP(k)=m; end end I=eye(n); if IP(k)=k for i=1:n p(i)=I(k,i); I(k,i)=I(IP(k),i); I(IP(k),i)=p(i); end A=I*A; b=I*b' b=b' end%进行消元计算 for i=k+1:n A(i,k)=A(i,k)/A(k,k); b(i)=b(i)-A(i,k)*b(k); for j=k+1:n A(i,j)=A(i,j)-A(i,k)*A(k,j); end det=det*A(k,k); P=I*P;end%回

6、代求解x(n)=b(n)/A(n,n);for i=n-1:-1:1 sum=0.0; for j=i+1:n sum=sum+A(i,j)*x(j); end x(i)=(b(i)-sum)/A(i,i);enddet=det*A(n,n);if det=0 fprintf('The equations have no unique solution!');end%输出PA=LU中的L,U的信息for i=1:n for j=1:n if i<j U(i,j)=A(i,j); elseif i=j L(i,j)=1; U(i,j)=A(i,j); else L(i,j)

7、=A(i,j); end endendend输入n，A，b，det=1.0k=1,2，n-1按列选主元c0=|c0|<ik=k换行消元计算输出停机回代求解输出计算解及行列式值b（i）（i=1,2，n）及det停机YYNN附图15算例分析1、测试示例>> A=1 23 45 6;>> b=3 7 11;>> x L U IP P=gauss(A,b)Error! Please input again!>> A=1 2 30 0 04 5 6;>> b=3 7 11;>> x L U IP P=gauss(A,b)The

8、 equations have no unique solution!x = NaN -Inf InfL = 1.0000 0 0 0.2500 1.0000 0 0 0 1.0000U = 4.0000 5.0000 6.0000 0 0.7500 1.5000 0 0 0IP = 3 3P = 1 0 0 0 0 1 0 1 02、计算过程（1）首先输入系数矩阵A和矩阵b>> A=1 1 1 1 1 1 12 1 1 1 1 1 13 2 1 1 1 1 14 3 2 1 1 1 15 4 3 2 1 1 16 5 4 3 2 1 17 6 5 4 3 2 1;>>

9、 b=7 8 10 13 17 22 28;（2）输出结果>> x L U IP P=gauss(A,b)x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000L = 1.0000 0 0 0 0 0 0 0.2857 1.0000 0 0 0 0 0 0.4286 0.8000 1.0000 0 0 0 0 0.5714 0.6000 0.7500 1.0000 0 0 0 0.7143 0.4000 0.5000 0.6667 1.0000 0 0 0.1429 -0.2000 -0.2500 -0.3333 -0.5000 1

10、.0000 0 0.8571 0.2000 0.2500 0.3333 0.5000 -1.0000 1.0000U = 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 0 -0.7143 -0.4286 -0.1429 0.1429 0.4286 0.7143 0 0 -0.8000 -0.6000 -0.4000 -0.2000 0.0000 0 0 0 -0.7500 -0.5000 -0.2500 0.0000 0 0 0 0 -0.6667 -0.3333 -0.0000 0 0 0 0 0 0.5000 1.0000 0 0 0 0

11、 0 0 1.0000IP = 7 2 3 4 5 7P = 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 1 06讨论与结论1、 时间复杂度：>> tic;x L U IP P=gauss(A,b);tocElapsed time is 0. seconds.2、 程序优化初次编程时，没有考虑到给一个变量赋初值的情况。虽然在MATLAB中变量不赋初值是完全允许的，但是由于一个变量中含有多个元素时，每次改变该数组的长度，便会增加计算机时间。另外，给程序加上一定的判断条件及报错信息，一定程度上有程序优化的作用。因此，本程序中的以下程序段都起到了程序优化的作用。%当输入的系数矩阵不为方阵，或方阵维数与b不符时，报错if p=q|p=n fprintf('Error! Please input again!');end%为提高运行