列主元高斯消去法和列主元三角分解法解线性方程

列主元高斯消去法和列主元三角分解法解线性方程计算方法实验报告1【课题名称】用列主元高斯消去法和列主元三角分解法解线性方程【目的和意义】高斯消去法是一个古老的求解线性方程组的方法，但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。

用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A 约化为具有简单形式的矩阵（上三角矩阵、单位矩阵等），而三角形方程组则可以直接回带求解用高斯消去法解线性方程组b Ax =（其中A ∈Rn ×n ）的计算量为：乘除法运算步骤为32(1)(1)(21)(1)(1)262233n n n n n n n n n n nMD n ----+=+++=+-，加减运算步骤为(1)(21)(1)(1)(1)(25)6226n n n n n n n n n n AS -----+=++=。

相比之下，传统的克莱姆法则则较为繁琐，如求解20阶线性方程组，克莱姆法则大约要19510⨯次乘法，而用高斯消去法只需要3060次乘除法。

在高斯消去法运算的过程中，如果出现abs(A(i,i))等于零或过小的情况，则会导致矩阵元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠，所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法，从而使计算结果更加精确。

2、列主元三角分解法高斯消去法的消去过程，实质上是将A 分解为两个三角矩阵的乘积A=LU ，并求解Ly=b 的过程。

回带过程就是求解上三角方程组Ux=y 。

所以在实际的运算中，矩阵L 和U 可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法采用选主元的方式与列主元高斯消去法一样，也是为了避免除数过小，从而保证了计算的精确度【计算公式】1、 列主元高斯消去法设有线性方程组Ax=b ，其中设A 为非奇异矩阵。

方程组的增广矩阵为1112112122221[,]nn l a a a a b aa ab ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥a b LL M L LM MM M M M M M M第1步（k=1）：首先在A 的第一列中选取绝对值最大的元素1l a ，作为第一步的主元素:然后交换（A ，b ）的第1行与第l 行元素，再进行消元计算。

设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组 A(k)x=b(k)111max 0l i i na a ≤≤=≠(1)(1)(1)(1)(1)1112111(2)(2)(2)(2)22222()(()1)()()()()()1,1()(,)()[,][,] k k kk nkknkn k k k k k kk kn k k k k n k k k n nn a a a a b a a a b a a b a b b a a a +++⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥→=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦A b A b L M LM OM M M LM M M M MMM L第k 步计算如下： 对于k=1，2，…，n-1 （1）按列选主元：即确定t 使（2）如果t ≠k ，则交换[A ，b]第t 行与第k 行元素。

（3）消元计算消元乘数mik 满足:（4）回代求解()()max 0k k tkikk i na a≤≤=≠,(1,,)ikik ik kka a m i k n a ←=-=+L , (,1,,), (1,,)ijij ik kj i i ik k a a m a i j k n b b m b i k n ←+=+⎧⎨←+=+⎩L L ⎪⎪⎩⎪⎪⎨⎧--=-←←∑+=)1,,2,1(,)(1Λn n i a x a b x a b x ii n i j j ij i i nnn n2、 列主元三角分解法 对方程组的增广矩阵 经过k-1步分解后，可变成如下形式：第k 步分解，为了避免用绝对值很小的数kku 作除数，引进量1111 (,1,,;1,2,,) ()/ (1,2,,;1,2,,)k kj kj km mj m k ik ik im mk kkm u a l u j k k n k n l a l u u i k k n k n -=-=⎧=-=+=⎪⎪⎨⎪=-=++=⎪⎩∑∑L L L L11(,1,,)k i ik im mk m s a l u i k k n -==-=+∑L ，于是有kku =ks 。

如果 ，则将矩阵的第t 行与第k 行元素互换，将（i ，j ）位置的新元素仍记为jjl 或jja ，然后再做第k 步分解，这时[,]A A b =11121,11111222,122221,11,1,1,211,11,2121,112,112,1k k k k k k kj n k k j n k k k i i i k n n kkkj kn k ik ij in i nknjk k k j k n n nnk k n a a a b A a u u u u u u y l l l l l l l l l l l u u u u u y u u u u y a a b a a b l a -------------⎡→⎣L L L L L L M M M LM M M L L M M L MM M L L M M M M M L L MM L LL M M M L⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎦max tik i ns s ≤≤=【列主元高斯消去法程序流程图】【列主元高斯消去法Matlab 主程序】function()/ 1,2,,)1 (1,2,,),kkk k t ik i k iku s s s l s s i k k n l i k k n ===++≤=++L L 即交换前的，（且x=gauss1(A,b,c) %列主元法高斯消去法解线性方程Ax=bif(length(A)~=length(b)) %判断输入的方程组是否有误disp('输入方程有误！')return;enddisp('原方程为AX=b：') %显示方程组Abdisp('------------------------')n=length(A);fork=1:n-1 %找列主元[p,q]=max(abs(A(k:n,k))); %找出第k列中的最大值，其下标为[p,q]q=q+k-1; %q在A(k:n,k)中的行号转换为在A中的行号if abs(p)<cdisp('列元素太小，det(A)≈0');break;elseif q>ktemp1=A(k,:);%列主元所在行不是当前行，将当前行与列主A(k,:)=A(q,:);元所在行交换(包括b) A(q,:)=temp1;temp2=b(k,:);b(k,:)=b(q,:);b(q,:)=temp2;end%消元for i=k+1:nm(i,k)=A(i,k)/A(k,k); %A(k, k)将A(i,k)消为0所乘系数A(i,k:n)=A(i,k:n)-m(i,k)*A(k,k:n); %第i行消元处理b(i)=b(i)-m(i,k)*b(k); %b消元处理endenddisp('消元后所得到的上三角阵是')A %显示消元后的系数矩阵b(n)=b(n)/A(n,n);%回代求解for i=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n)*b(i+1:n)))/A(i, i);endclear x;disp('AX=b的解x是')x=b;【调用函数解题】【列主元三角分解法程序流程图】【列主元三角分解法Matlab主程序】①自己编的程序：function x=PLU(A,b,eps) %定义函数列主元三角分解法函数if (length(A)~=length(b)) %判断输入的方程组是否有误 disp('输入方程有误！')return;enddisp('原方程为AX=b：') %显示方程组Abdisp('------------------------')n=length(A);A=[A b]; %将A与b合并，得到增广矩阵for r=1:nif r==1for i=1:n[c d]=max(abs(A(:,1))); %选取最大列向量，并做行交换 if c<=eps %最大值小于e，主元太小，程序结束break;elseendd=d+1-1;p=A(1,:);A(1,:)=A(d,:);A(d,:)=p;A(1,i)=A(1,i);endA(1,2:n)=A(1,2:n);A(2:n,1)=A(2:n,1)/A(1,1); %求u(1,i)elseu(r,r)=A(r,r)-A(r,1:r-1)*A(1:r-1,r); %按照方程求取u(r,i) if abs(u(r,r))<=eps %如果u（r,r）小于e，则交换行p=A(r,:);A(r,:)=A(r+1,:);A(r+1,:)=p;elseendfor i=r:nA(r,i)=A(r,i)-A(r,1:r-1)*A(1:r-1,i); %根据公式求解，并把结果存在矩阵A中endfor i=r+1:nA(i,r)=(A(i,r)-A(i,1:r-1)*A(1:r-1,r))/A(r,r); %根据公式求解，并把结果存在矩阵A中endendendy(1)=A(1,n+1);for i=2:nh=0;for k=1:i-1h=h+A(i,k)*y(k);endy(i)=A(i,n+1)-h; %根据公式求解y(i)endx(n)=y(n)/A(n,n);for i=n-1:-1:1h=0;for k=i+1:nh=h+A(i,k)*x(k);endx(i)=(y(i)-h)/A(i,i); %根据公式求解x(i)endAdisp('AX=b的解x是')x=x'; %输出方程的解②可直接得到P,L,U并解出方程解的的程序（查阅资料得子函数PLU1，其作用是将矩阵A分解成L乘以U的形式。

PLU2为调用PLU1解题的程序，是自己编的）(Ⅰ).function [l,u,p]=PLU1(A) %定义子函数，其功能为列主元三角分解系数矩阵A[m,n]=size(A); %判断系数矩阵是否为方阵if m~=nerror('矩阵不是方阵')returnendif det(A)==0 %判断系数矩阵能否被三角分解error('矩阵不能被三角分解')endu=A;p=eye(m);l=eye(m); %将系数矩阵三角分解，分别求出P,L,Ufor i=1:mfor j=i:mt(j)=u(j,i);for k=1:i-1t(j)=t(j)-u(j,k)*u(k,i);endenda=i;b=abs(t(i));for j=i+1:mif b<abs(t(j))b=abs(t(j));a=j;endendif a~=ifor j=1:mc=u(i,j);u(i,j)=u(a,j);u(a,j)=c;endfor j=1:mc=p(i,j);p(i,j)=p(a,j);p(a,j)=c;endc=t(a);t(a)=t(i);t(i)=c;endu(i,i)=t(i);for j=i+1:mu(j,i)=t(j)/t(i);endfor j=i+1:mfor k=1:i-1u(i,j)=u(i,j)-u(i,k)*u(k,j);endendendl=tril(u,-1)+eye(m);u=triu(u,0)(Ⅱ).function x=PLU2(A,b) %定义列主元三角分解法的函数[l,u,p]=PLU1(A) %调用PLU分解系数矩阵Am=length(A); %由于A左乘p，故b也要左乘pv=b;for q=1:mb(q)=sum(p(q,1:m)*v(1:m,1));endb(1)=b(1) %求解方程Ly=bfor i=2:1:mb(i)=(b(i)-sum(l(i,1:i-1)*b(1:i-1)));endb(m)=b(m)/u(m,m); %求解方程Ux=yfor i=m-1:-1:1b(i)=(b(i)-sum(u(i,i+1:m)*b(i+1:m)))/u(i,i);endclear x;disp('AX=b的解x是')x=b;【调用函数解题】①②【编程疑难】这是第一次用matlab编程，对matlab的语句还不是非常熟悉，因此在编程过程中，出现了许多错误提示。