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

,.计算方法实验报告1【课题名称】用列主元高斯消去法和列主元三角分解法解线性方程【目的和意义】高斯消去法是一个古老的求解线性方程组的方法，但由它改进得到的选主元的高斯消去精品文档放心下载法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。精品文档放心下载用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A约化为具精品文档放心下载有简单形式的矩阵（上三角矩阵、单位矩阵等），而三角形方程组则可以直接回带求解谢谢阅读用高斯消去法解线性方程组Axb（其中A∈Rn×n）的计算量为：乘除法运算步骤为精品文档放心下载n(n1)(n1)n(2n1)n(n1)n(n1)n3nMD26223n23，加减运算步骤为AS(n1)n(2n1)n(n1)n(n1)(n1)n(2n5)6226。相比之下，传统的克莱姆法则则较为繁琐，如求解20阶线性方程组，克莱姆法则大约要51019次乘法，而用高斯消去法只需要3060次乘除法。谢谢阅读在高斯消去法运算的过程中，如果出现abs(A(i,i))等于零或过小的情况，则会导致矩阵谢谢阅读元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠，所以目前计算机上常精品文档放心下载用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法，从而使计算结果更加精确。精品文档放心下载2、列主元三角分解法高斯消去法的消去过程，实质上是将A分解为两个三角矩阵的乘积A=LU，并求解Ly=b谢谢阅读的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中，矩阵L和U可以直精品文档放心下载接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的谢谢阅读,.简略，大大提高了运算速度，这就是三角分解法采用选主元的方式与列主元高斯消去法一样，也是为了避免除数过小，从而保证了计算的精谢谢阅读确度【计算公式】1、列主元高斯消去法设有线性方程组Ax=b，其中设A为非奇异矩阵。方程组的增广矩阵为感谢阅读aaLLaMb11121n1aaLLaMb21222n2MM[a,b]al1MMMMMMaaLLaMbn1n2nnn第1步（k=1）：首先在A的第一列中选取绝对值最大的元素al1，作为第一步的主元素:精品文档放心下载amaxa0l11ini1然后交换（A，b）的第1行与第l行元素，再进行消元计算。精品文档放心下载设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与感谢阅读原方程组等价的方程组 A(k)x=b(k)a(1) 11[A,b][A(k),b(k)]精品文档放心下载

a(1)La(1)a(1)Mb(1)121k1n1a(2)La(2)a(2)Mb(2)222k2n2OMMMa(k)La(k)Mb(k)kkknka(k)a(k)Mb(k)k1,kk1,nk1MMMM,.第k步计算如下：对于k=1，2，…，n-1a(k)maxa(k)0tkkinik（1）按列选主元：即确定t使（2）如果t≠k，则交换[A，b]第t行与第k行元素。感谢阅读（3）消元计算ama,(ik1,L,n)ikikikakkaama, ij ij ik kj bbmb,i i ik k

(i,jk1,L,n)(ik1,L,n)消元乘数mik满足:（4）回代求解xbnnann(bnax)ijjixji1,(in1,n2,,1)iaii,.2、列主元三角分解法A[A,b]对方程组的增广矩阵经过k-1步分解后，可变成如下形式：uuLuuLuLuy11121,k11k1j1n1l21u22Lu2,k1u2kLu2jLu2ny2MMMMMMMlk1,1lLuuLuLuyk1,2k1,k1k1,kk1,jk1,nk1AllLlaLaLabk1k2k,k1kkkjknkMMMMMMMllLlaLaLabi1i2i,k1ikijiniMMMMMMMllLlaLaLabn1n2n,k1nknjnnn第k步分解，为了避免用绝对值很小的数ukk作除数，引进量谢谢阅读ualu(jk,k1,L,n;k1,2,L,n)k1kjkjkmmjm1l(ak1u)/u(ik1,k2,L,n;k1,2,L,n)likikimmkkkm1k1lu(ik,k1,L,n)saiikimmk，于是有ukk=sk。如果smaxs，则将矩阵的第m1tkini行与第k行元素互换，将（i，j）位置的新元素仍记为ljj或ajj，然后再做第k步分解，这感谢阅读时,.u s (s即交换前的s)，kk k k ts/s（ik1,k2,L,n)感谢阅读ik i kl1(ik1,k2,L,n),感谢阅读ik【列主元高斯消去法程序流程图】【列主元高斯消去法Matlab主程序】functionx=gauss1(A,b,c) %列主元法高斯消去法解线性方程Ax=b感谢阅读,.if(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中的行感谢阅读号ifabs(p)<cdisp('列元素太小，det(A)≈0');break;elseifq>ktemp1=A(k,:); %列主元所在行不是当前行，将当前行与列谢谢阅读主,.A(k,:)=A(q,:); 元所在行交换(包括b)精品文档放心下载A(q,:)=temp1;temp2=b(k,:);b(k,:)=b(q,:);b(q,:)=temp2;end%消元fori=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); %回代求解fori=n-1:-1:1b(i)=(b(i)-sum(A(i,i+1:n)*b(i+1:n)))/A(i,i);精品文档放心下载endclearx;,.disp('AX=b的解x是')x=b;【调用函数解题】,.【列主元三角分解法程序流程图】,.【列主元三角分解法Matlab主程序】①自己编的程序：functionx=PLU(A,b,eps) %定义函数列主元三角分解法函谢谢阅读数if(length(A)~=length(b)) %判断输入的方程组是否有误精品文档放心下载,.disp('输入方程有误！')return;enddisp('原方程为AX=b：') %显示方程组Abdisp('------------------------')感谢阅读n=length(A);A=[Ab]; %将A与b合并，得到增广矩阵谢谢阅读forr=1:nifr==1fori=1:n[cd]=max(abs(A(:,1))); %选取最大列向量，并做行交换精品文档放心下载ifc<=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)谢谢阅读ifabs(u(r,r))<=eps %如果u（r,r）小于e，则交换行感谢阅读p=A(r,:);A(r,:)=A(r+1,:);A(r+1,:)=p;elseendfori=r:nA(r,i)=A(r,i)-A(r,1:r-1)*A(1:r-1,i); %根据公式求解，并把结果存在矩阵A谢谢阅读中endfori=r+1:nA(i,r)=(A(i,r)-A(i,1:r-1)*A(1:r-1,r))/A(r,r);%根据公式求解，并把结果存在矩阵A中精品文档放心下载endendend,.y(1)=A(1,n+1);fori=2:nh=0;fork=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);fori=n-1:-1:1h=0;fork=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); %判断系数矩阵是否为方感谢阅读阵ifm~=nerror('矩阵不是方阵')returnendifdet(A)==0%判断系数矩阵能否被三角分解谢谢阅读error('矩阵不能被三角分解')endu=A;p=eye(m);l=eye(m); %将系数矩阵三角分解，分别求出谢谢阅读P,L,Ufori=1:mforj=i:mt(j)=u(j,i);fork=1:i-1,.t(j)=t(j)-u(j,k)*u(k,i);endenda=i;b=abs(t(i));forj=i+1:mifb<abs(t(j))b=abs(t(j));a=j;endendifa~=iforj=1:mc=u(i,j);u(i,j)=u(a,j);u(a,j)=c;endforj=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);forj=i+1:mu(j,i)=t(j)/t(i);endforj=i+1:mfork=1:i-1u(i,j)=u(i,j)-u(i,k)*u(k,j);谢谢阅读endendendl=tril(u,-1)+eye(m);u=triu(u,0)(Ⅱ).functionx=PLU2(A,b) %定义列主元三角分解法的精品文档放心下载函数[l,u,p]=PLU1(A) %调用PLU分解系数矩谢谢阅读阵A,.m=length(A); %由于A左乘p，故b也要左感谢阅读乘pv=b;forq=1:mb(q)=sum(p(q,1:m)*v(1:m,1));精品文档放心下载endb(1)=b(1) %求解方程Ly=bfori=2:1:mb(i)=(b(i)-sum(l(i,1:i-1)*b(1:i-1)));精品文档放心下载endb(m)=b(m)/u(m,m); %求解方程Ux=yfori=m-1:-1:1b(i)=(b(i)-sum(u(i,i+1:m)*b(i+1:m)))/u(i,i);谢谢阅读endclearx;disp('AX=b的解x