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

实验名称： 列主元消去法解方程组 引言 我们知道，高斯消去法是一个古老的解线性方程组的方法。而在用高斯消去法解Ax=b时，其中设A为非奇异矩阵，可能出现的情况，这时必须进行带行交换的高斯消去法。但在实际计算中即使但其绝对值很小时，用作除数，会导致中间结果矩阵元素数量级严重增长和舍入误差的扩散，使得最后的结果不可靠。因此，小主元可能导致计算的失败，我们应该避免采用绝对值很小的主元素。为此，我们在高斯消去法的每一步应该在系数矩阵或消元后的低阶矩阵中选取绝对值最大的元素作为主元素，保持乘数，以便减少计算过程中舍入误差对计算解的影响。 一种方式是完全主元消去法，这种消去法是在每次选主元时，选择为主元素。这种方法是解低阶稠密矩阵方程组的有效方法，但这种方法在选取主元时要花费一定的计算机时间。实际计算中我们常采用部分选主元的的消去法。列主元消去法即在每次选主元时，仅依次按列选取绝对值最大的元素作为主元素，且仅交换两行，再进行消元计算。 实验目的和要求 运用matlab编写一个.m文件，要求用列主元消去法求解方程组（实现PA=LU）： 要求输出以下内容： 计算解x； L,U； 整形数组IP（i）（i=1,2，…，n-1）（记录主行信息） 算法原理与流程图 算法原理 设有线性方程组Ax=b，A为非奇异矩阵。方程组的增广矩阵为 第1步（k=1）：A的第一列中选取绝对值最大的元素，作为第一步的主元素: ，然后交换（A，b）的第1行与第i1行元素，再进行消元计算。 设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组 第k步计算如下： 对于k=1，2，…，n-1 （2）如果，则A为非奇异矩阵，停止计算。 （3）如果ik≠k，则交换[A，b]第ik行与第k行元素。 （4）消元计算 消元乘数满足: （5）回代求解 计算解在常数项b(n)内得到。 流程图见图1 程序代码及注释 %列主元消去法解方程组Ax=b，实现PA=LU function [x,L,U,IP,P] =gauss(A,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+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 %回代求解 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); end det=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 ij