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

计算方法实验报告 实验名称： 实验2 列主元素消去法解方程组 引言 工程实际问题中，线型方程的系数矩阵一般为低阶稠密矩阵和大型稀疏矩阵。用高斯消去法解Ax=b时，可能出现很小，用作除数会导致中间结果矩阵元素数量级严重增长和舍入误差的扩散，使结果不可靠；采用选主元素的三角分解法可以避免此类问题。高斯消去法的消去过程，实质上是将A分解为两个三角矩阵的乘积A=LU，并求解Ly=b的过程。回带过程就是求解上三角方程组Ux=y。所以在实际的运算中，矩阵L和U可以直接计算出，而不需要任何中间步骤，从而在计算过程中将高斯消去法的步骤进行了进一步的简略，大大提高了运算速度，这就是三角分解法。采用选主元的方式与列主元高斯消去法一样，也是为了避免除数过小，从而保证了计算的精确度。 实验目的和要求 通过列主元素消去法求解线性方程组，实现PA=LU。要求计算解x，L，U，整形数组IP(i)，(i=1,2,…,)（记录主行信息）。 算法原理与流程图 原理 将A分解为两个三角矩阵的乘积A=LU。 对方程组的增广矩阵经过k-1步分解后，可变成如下形式： 第k步分解，为了避免用绝对值很小的数作除数，引进量 ，于是有=。如果，则将矩阵的第t行与第k行元素互换，将（i,j）位置的新元素仍记为或，然后再做第k步分解，这时 流程图 程序代码及注释 %运用列主元素消去法求解方程组并实现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; %转置为列向量 算例分析 解方程组 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.0000 L = 1.0000 0 0 0 0