计算方法实验报告习题2(浙大版).doc

专业：电气工程及其自动化姓名： 李X 计算方法实验报告实验名称： 实验2 列主元素消去法解方程组 1 引言工程实际问题中，线型方程的系数矩阵一般为低阶稠密矩阵和大型稀疏矩阵。用高斯消去法解Ax=b时，可能出现很小，用作除数会导致中间结果矩阵元素数量级严重增长和舍入误差的扩散，使结果不可靠；采用选主元素的三角分解法可以避免此类问题。高斯消去法的消去过程，实质上是将A分解为两个三角矩阵的乘积A=LU，并求解Ly=b的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中，矩阵L和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法。采用选主元的方式与列主元高斯消去法一样，也是为了避免除数过小，从而保证了计算的精确度。2 实验目的和要求通过列主元素消去法求解线性方程组，实现PA=LU。要求计算解x，L，U，整形数组IP(i)，(i=1,2,)（记录主行信息）。3 算法原理与流程图（1） 原理将A分解为两个三角矩阵的乘积A=LU。对方程组的增广矩阵经过k-1步分解后，可变成如下形式：第k步分解，为了避免用绝对值很小的数作除数，引进量，于是有=。如果，则将矩阵的第t行与第k行元素互换，将（i,j）位置的新元素仍记为或，然后再做第k步分解，这时（2） 流程图4 程序代码及注释%运用列主元素消去法求解方程组并实现PA=LUfunction x,L,U,P,IP=lzylu(A,b) %定义函数m,n=size(A); %计算系数矩阵A的行列数if m=n error(系数矩阵不是方阵); return;endif det(A)=0 %计算矩阵A的行列式，若为零则A是奇异矩阵 error(矩阵不能被三角形分解)endu=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的过程endP=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; %转置为列向量5 算例分析解方程组 A=1 1 1 1 1 1 1;2 1 1 1 1 1 1;3 2 1 1 1 1 1;4 3 2 1 1 1 1;5 4 3 2 1 1 1;6 5 4 3 2 1 1;7 6 5 4 3 2 1; b=7 8 10 13 17 22 28; x,L,U,P,IP=lzylu(A,b)x = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000L = 1.0000 0 0 0 0 0 0 0.2857 1.0000 0 0 0 0 0 0.4286 0.8000 1.0000 0 0 0 0 0.5714 0.6000 0.7500 1.0000 0 0 0 0.7143 0.4000 0.5000 0.6667 1.0000 0 0 0.1429 -0.2000 -0.2500 -0.3333 -0.5000 1.0000 0 0.8571 0.2000 0.2500 0.3333 0.5000 -1.0000 0.0000U = 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 0 -0.7143 -0.4286 -0.1429 0.1429 0.4286 0.7143 0 0 -0.8000 -0.6000 -0.4000 -0.2000 0.0000 0 0 0 -0.7500 -0.5000 -0.2500 0.0000 0 0 0 0 -0.6667 -0.3333 -0.0000 0 0 0 0 0 0.5000 0.0000 0 0 0 0 0 0 0.0000P = 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 0IP = 7 2 3 4 5 7本程序将方程组系数矩阵进行列主元素三角形分解，得到矩阵L,U,P和记录主行变换的数组IP(i)。6 讨论与结论（1）程序lzylu.m具有以下特点：1. 当系数矩阵是奇异矩阵时输出提示错误信息，停止运算2. 当系数矩阵不是方阵时，提示错误并停止运算3. 能够输出变换后的L,U,P矩阵，并记录主行变换信息（2）将自编程序与其他程序从运行耗时方面进行比较序号时间(s)函数x,L,U,P,IP=lzylu(A,b)x=PLU(A,b,eps)x=PLU2(A,b)x=PLU3(A,b)10.058980.0963690.123020.09516220.0519650.107920.118770.0906530.0529320.0903740.122410.09879140.0547890.0976870.124020.098823AVERAGE0.05466650.09808750.1220550.0958565其中PLU3是指将PLU2中第2行的调用PLU1函数改为调用lu函数后的文件，其代码并未在本报告内写出。a.比较一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: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(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:nA(r,i)=A(r,i)-A(r,1:r-1)*A(1:r-1,i); %根据公式求解，并把结果存在矩阵A中 end for i=r+1:n A(i,r)=(A(i,r)-A(i,1:r-1)*A(1:r-1,r)/A(r,r); %根据公式求解，并把结果存在矩阵A中 end endendy(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)endx(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)endAdisp(AX=b的解x是)x=x;%输出方程的解此方法的优点之一在于，将A与b进行了合并，称为增广矩阵，便于后续交换行的运算。但正因如此，导致后面没有输出L和U矩阵；此外，用max函数求第一列的最大值，比for循环要快，但其他列并没有求最大值，而是比较了对角元素和eps的大小关系，决定是否换行，所以理论上此方法未完全达到选取列主元素的功能。b.比较二().function l,u,p=PLU1(A) %定义子函数，其功能为列主元三角分解系数矩阵Am,n=size(A); %判断系数矩阵是否为方阵if m=n error(矩阵不是方阵) returnendif det(A)=0 %判断系数矩阵能否被三角分解 error(矩阵不能被三角分解)endu=A;p=eye(m);l=eye(m); %将系数矩阵三角分解，分别求出P,L,Ufor 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); end 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)=a; 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)*u(k,j); end endendl=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);endb(1)=b(1) %求解方程Ly=bfor i=2:1:m b(i)=(b(i)-sum(1:i-1)*b(1:i-1);endb(m)=b(m)/u(m,m); %求解方程Ux=yfor i=m-1:-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;PLU1文件中，用了太多的for循环语句，所以耗时较长。此方法调用了两个程序，耗时要比较长，若将第二个程序