列主元高斯消去法和列主元三角分解法解线性方程(13页).doc

《列主元高斯消去法和列主元三角分解法解线性方程(13页).doc》由会员分享，可在线阅读，更多相关《列主元高斯消去法和列主元三角分解法解线性方程(13页).doc（12页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、-计算方法实验报告1 【课题名称】用列主元高斯消去法和列主元三角分解法解线性方程【目的和意义】高斯消去法是一个古老的求解线性方程组的方法，但由它改进得到的选主元的高斯消去法则是目前计算机上常用的解低阶稠密矩阵方程组的有效方法。用高斯消去法解线性方程组的基本思想时用矩阵行的初等变换将系数矩阵A约化为具有简单形式的矩阵（上三角矩阵、单位矩阵等），而三角形方程组则可以直接回带求解用高斯消去法解线性方程组（其中ARnn）的计算量为：乘除法运算步骤为，加减运算步骤为。相比之下，传统的克莱姆法则则较为繁琐，如求解20阶线性方程组，克莱姆法则大约要次乘法，而用高斯消去法只需要3060次乘除法。在高斯消去法运

2、算的过程中，如果出现abs(A(i,i)等于零或过小的情况，则会导致矩阵元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠，所以目前计算机上常用的解低阶稠密矩阵方程的快速有效的方法时列主元高斯消去法，从而使计算结果更加精确。2、列主元三角分解法高斯消去法的消去过程，实质上是将A分解为两个三角矩阵的乘积A=LU，并求解Ly=b的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中，矩阵L和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法采用选主元的方式与列主元高斯消去法一样，也是为了避免除数过小

3、，从而保证了计算的精确度【计算公式】1、 列主元高斯消去法设有线性方程组Ax=b，其中设A为非奇异矩阵。方程组的增广矩阵为 第1步（k=1）：首先在A的第一列中选取绝对值最大的元素，作为第一步的主元素: 然后交换（A，b）的第1行与第l行元素，再进行消元计算。 设列主元素消去法已经完成第1步到第k-1步的按列选主元，交换两行，消元计算得到与原方程组等价的方程组 A(k)x=b(k) 第k步计算如下： 对于k=1，2，n-1 （1）按列选主元：即确定t使 （2）如果tk，则交换A，b第t行与第k行元素。 （3）消元计算 消元乘数mik满足: （4）回代求解2、 列主元三角分解法对方程组的增广矩阵

4、 经过k-1步分解后，可变成如下形式：第k步分解，为了避免用绝对值很小的数作除数，引进量，于是有=。如果 ，则将矩阵的第t行与第k行元素互换，将（i，j）位置的新元素仍记为或，然后再做第k步分解，这时【列主元高斯消去法程序流程图】【列主元高斯消去法Matlab主程序】function x=gauss1(A,b,c) %列主元法高斯消去法解线性方程Ax=bif (length(A)=length(b) %判断输入的方程组是否有误 disp(输入方程有误！) return;enddisp(原方程为AX=b：) %显示方程组Abdisp(-)n=length(A);for k=1:n-1 %找列主元

5、 p,q=max(abs(A(k:n,k); %找出第k列中的最大值，其下标为p,q q=q+k-1; %q在A(k:n,k)中的行号转换为在A中的行号 if abs(p)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:n m(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

6、); %第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:1 b(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)

7、 %判断输入的方程组是否有误 disp(输入方程有误！) return;enddisp(原方程为AX=b：) %显示方程组Abdisp(-)n=length(A);A=A b; %将A与b合并，得到增广矩阵for r=1:n if r=1 for i=1:n c d=max(abs(A(:,1); %选取最大列向量，并做行交换 if c=eps %最大值小于e，主元太小，程序结束 break; else end d=d+1-1; p=A(1,:); A(1,:)=A(d,:); A(d,:)=p; A(1,i)=A(1,i); end A(1,2:n)=A(1,2:n); A(2:n,1)=A

8、(2:n,1)/A(1,1); %求u(1,i) else u(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; else end for i=r:n A(r,i)=A(r,i)-A(r,1:r-1)*A(1:r-1,i); %根据公式求解，并把结果存在矩阵A中 end for i=r+1:nA(i,r)=(A(i,r)-A(i,1:r-1)*A(1:r-1,r)/A(r,r); %根据公式求解，并

9、把结果存在矩阵A中 end endend y(1)=A(1,n+1); for i=2:n h=0; for k=1:i-1 h=h+A(i,k)*y(k); end y(i)=A(i,n+1)-h; %根据公式求解y(i) end x(n)=y(n)/A(n,n); for i=n-1:-1:1 h=0; for k=i+1:n h=h+A(i,k)*x(k); end x(i)=(y(i)-h)/A(i,i); %根据公式求解x(i) end Adisp(AX=b的解x是)x=x; %输出方程的解可直接得到P,L,U并解出方程解的的程序（查阅资料得子函数PLU1，其作用是将矩阵A分解成L乘

10、以U的形式。PLU2为调用PLU1解题的程序，是自己编的）().function l,u,p=PLU1(A) %定义子函数，其功能为列主元三角分解系数矩阵A m,n=size(A); %判断系数矩阵是否为方阵 if m=n error(矩阵不是方阵) return end if det(A)=0 %判断系数矩阵能否被三角分解 error(矩阵不能被三角分解) end u=A;p=eye(m);l=eye(m); %将系数矩阵三角分解，分别求出P,L,U for i=1:m for j=i:m t(j)=u(j,i); for k=1:i-1 t(j)=t(j)-u(j,k)*u(k,i); e

11、nd end a=i;b=abs(t(i); for j=i+1:m if babs(t(j) b=abs(t(j); a=j; end end if a=i for j=1:m c=u(i,j); u(i,j)=u(a,j); u(a,j)=c; end for j=1:m c=p(i,j); p(i,j)=p(a,j); p(a,j)=c; end c=t(a); t(a)=t(i); t(i)=c; end u(i,i)=t(i); for j=i+1:m u(j,i)=t(j)/t(i); end for j=i+1:m for k=1:i-1 u(i,j)=u(i,j)-u(i,k)

12、*u(k,j); end end end l=tril(u,-1)+eye(m); u=triu(u,0)().function x=PLU2(A,b) %定义列主元三角分解法的函数l,u,p=PLU1(A) %调用PLU分解系数矩阵A m=length(A); %由于A左乘p，故b也要左乘pv=b;for q=1:m b(q)=sum(p(q,1:m)*v(1:m,1);end b(1)=b(1) %求解方程Ly=b for i=2:1:m b(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:-

13、1:1 b(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的语句还不是非常熟悉，因此在编程过程中，出现了许多错误提示。并且此次编程的两种方法对矩阵的运算也比较复杂。问题主要集中在循环控制中，循环次数多了一次或者缺少了一次，导致数据错误，一些基本的编程语句在语法上也会由于生疏而产生许多问题，但是语句的错误由于系统会提示，比较容易进行修改，数据计算过程中的一些逻辑错误，比如循环变量的控制，这些系统不会提示错误，需要我们细心去发现错误，不断修正，调试。-第 12 页-