计算方法实验报告习题2浙大版.docVIP

计算方法实验报告 实验名称： 实验2 列主元素消去法解方程组 引言 实验目的和要求算法原理与流程图 原理 流程图 程序代码及注释 %运用列主元素消去法求解方程组并实现PA=LU function [x,L,U,P,IP]=lzylu(A,b) %定义函数 [m,n]=size(A); %计算系数矩阵A的行列数 if m~=n error(系数矩阵不是方阵); return; end if det(A)==0 %计算矩阵A的行列式，若为零则A是奇异矩阵 error(矩阵不能被三角形分解) end u=A;v=b; q=eye(m);s=q; %定义单位矩阵 for j=1:m-1 %按照列循环 t=abs(u(j,j)); mj=j; for i=j:m if abs(u(i,j))t mj=i; end end %比较大小，选列主元素 if u(mj,j)==0 error(系数矩阵奇异); return end B=u(j,:);c=v(j); u(j,:)=u(mj,:);v(j)=v(mj); u(mj,:)=B;v(mj)=c; %交换两行，将列主元素所在行移到第j行 IP(j)=mj; %记录互换，表示第j行与第IP(j)行互换 ts=s(j,:); s(j,:)=s(mj,:); s(mj,:)=ts; %置换矩阵进行相应行交换 u(j+1:m,j)=u(j+1:m,j)/u(j,j); %求出L矩阵第j列第j+1到第m行元素 u(j+1:m,j+1:m)=u(j+1:m,j+1:m)-u(j+1:m,j)*u(j,j+1:m); %解出u矩阵做完第j次消元后的矩阵 v(j+1:m)=v(j+1:m)-v(j)*u(j+1:m,j); %v也要对应系数进行运算，相比杜利特尔（Doolittle）分解法减少了l*y=b，y=u*x先要求解y再求x的过程 end P=s; L=tril(u,-1)+q; %对u取下三角，其余部分全部为0，得到L矩阵 U=triu(u); %对u取上三角，其余部分全部为0，得到U矩阵 x(m)=v(m)/u(m,m); for i=m-1:-1:1 sum=0; for k=i+1:m sum=sum+U(i,k)*x(k); end x(i)=(v(i)-sum)/u(i,i); end %利用U*x=v，求解x得到行向量 x=x; %转置为列向量 算例分析 讨论与结论 序号 时间(s) 函数 [x,L,U,P,IP]=lzylu(A,b) x=PLU(A,b,eps) x=PLU2(A,b) x=PLU3(A,b) 1 0.05898 0.096369 0.12302 0.095162 2 0.051965 0.10792 0.11877 0.09065 3 0.052932 0.090374 0.12241 0.098791 4 0.054789 0.097687 0.12402 0.098823 0.0546665 0.0980875 0.122055 0.0958565 ※其中PLU3是指将PLU2中第2行的调用PLU1函数改为调用lu函数后的文件，其代码并未在本报告内写出。 a.比较一 function x=PLU(A,b,eps) %定义函数列主元三角分解法函数 if (length(A)~=length(b)) %判断输入的方程组是否有误 disp(输入方程有误！) return; end disp(原方程为AX=b：) %显示方程组 A b disp(------------------------) 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