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

计算方法实验报告1【课题名称】用列主元高斯消去法和列主元三角分解法解线性方程【目的和意义】高斯消去法是一个古老的求解线性方程组的方法，但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A约化为具有简单形式的矩阵(上三角矩阵、单位矩阵等)，而三角形方程组则可以直接回带求解用高斯消去法解线性方程组Axb(其中AERnxn)的计算量为：乘除法运算步骤为n(n1)(n1)n(2n1)n(n1)n(n1)n3nMD26223n23，加减运算步骤为(n1)n(2n1)n(n1)n(n1)(n1)n(2n5)AS6226。相比之下，传统的克莱姆法则则较为繁琐，如求解20阶线性方程组，克莱姆法则大约要51019次乘法，而用高斯消去法只需要3060次乘除法。在高斯消去法运算的过程中，如果出现abs(A(i,|)于零或过小的情况，则会导致矩阵元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠，所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法，从而使计算结果更加精确。2、列主元三角分解法高斯消去法的消去过程，实质上是将A分解为两个三角矩阵的乘积A=LU，并求解Ly=b的过程。回带过程就是求解上三角方程组U*=y。所以在实际的运算中，矩阵L和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法采用选主元的方式与列主元高斯消去法一样，也是为了避免除数过小，从而保证了计算的精确度【计算公式】1、列主元高斯消去法设有线性方程组A*=b，其中设A为非奇异矩阵。方程组的增广矩阵为第1步(k=1)：首先在A的第一列中选取绝对值最大的元素气1，作为第一步的主元素:然后交换(A，b)的第1行与第l行元素，再进行消元计算。设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组A(k)*=b(k)第k步计算如下：对于k=1，2，…，n-1按列选主元：即确定t使|气"|max|")|0tk,.ik一，，k.in、如果t尹k，则交换［A，b］第t仃与第k仃元素。(4)回代求解2、列主元三角分解法_对方程组的增广矩阵经过k-A,步分解后，可变成如下形式：x，一―一yau第k步分解，为了避免用绝对值很小的数kk作除数，引进量(bnbax)k1i.ijjsax,-I—tt~，火...nn)l,n2,,1)1*'、"Ji，于是有"=I。如豪J则将矩留临第t行与第k行kin工__，-,、一,1a,,一，一、，元素互换，将(i,j)位置的新元素仍记为或打，然后再做第k步分解，这时【列主元高斯消去法程序流程图】【列主元高斯消去法Matlab主程序】function*=gauss1(A,b,列主元法高斯消去法解线性方程A*=bif(1ength(A)~=1ength(判断输入的方程组是否有误disp输入方程有误！')return;enddisp原方程为A*=b:')显示方程组Abdisp('')n=1ength(A);fork=1:n-1找列主元[p,q]=ma*(abs(A(k:n,k)我；%第k列中的最大值，其下标为[p,q]q=q+k-1;%q在A(k:n,中的行号转换为在A中的行号ifabs(p)<cdisp(列元素太小，det(AKO');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将Mk),l消为0所乘系数A(i,k:n)=A(i,k:n)-m(i,k)*A(]第ki行)消元处理b⑴=b(i)-m(i,k)*b(k);j%b^理endenddisp消元后所得到的上三角阵是')A%显示消元后的系数矩阵b(n)=b(n)/A(n,n)回代求解fori=n-1:-1:1b⑴=(b⑴-sum(A(i,i+1:n)*b(i+1:n)))/A(i,i);endclear*;dispCA*=bft解*是')*=b;【调用函数解题】【列主元三角分解法程序流程图】【列主元三角分解法Matlab主程序】①自己编的程序：function*=PLU(A,b,eps)

function*=PLU(A,b,eps)

if(length(A)=length(b))disp输入方程有误！')%定义函数列主元三角分解法函数%判断输入的方程组是否有误return;enddisp原方程为A*=b：')%显示方程组bdisp('')n=length(A);A=[Ab];将A与b合并，得到增广矩阵forr=1:nifr==1fori=1:n[cd]=ma*(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:i%按膈方程求取u(r,i)ifabs(u(r,r))<=ePWu(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:根据公)式求解，并把结果存在矩阵A中endfori=r+1:nA(i,r)=(A(i,r)-A(i,1:r-1)*A(1:r-1,根据公式求解，并把结果存在矩阵A中endendendy(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)end*(n)=y(n)/A(n,n);fori=n-1:-1:1h=0;fork=i+1:nh=h+A(i,k)**(k);end*(i)=(y(i)-h)/A(i,蹄根据公式求解*(i)enddisp('A*=b的解*是')*=*'；%输出方程的解②可直接得到P,L,U并解出方程解的的程序(查阅资料得子函数PLU1,其作用是将矩阵A分解成L乘以U的形式。PLU2为调用PLU1解题的程序，是自己编的)(I).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);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)(II).function*=PLU2(A,b)%定义列主元三角分解法的函数[l,u,p]=PLU1(A)%调用PLU分解系数矩阵Am=length(A);%由于A左乘p,故b也要左乘pv=b;forq=1:mb(q)=sum(p(q,1:m)*v(1:m,1));end81)=81)%求解方程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);%求解方程U*=yfori=m-1:-1:1b(i)=(b(i)-sum(u(i,i+1:m)*b(i+1:m)))/u(i,i);endclear*;disp('A*=b的解*是')*=b;【调用函数解题】①②【编程疑难】这是第一次用matlab编程，对