浙大计算方法上机报告.docVIP

学号：3100300038 姓名： 专业： 作业1：用列主元高斯消去法和列主元三角分解法解P227页第3题 1. 列主元高斯消去法 目的：用高斯消去法解Ax=b时，其中设A为非奇异矩阵，可能出现a=0 情况，这时必须进行带行交换的高斯消去法。 但在实际计算中即使 但其绝对值很小时，用作除数，会导致中间结果矩阵A(k)元素数量级严重增长和舍入误差的扩散，使得最后的计算结果不可靠。 列主元高斯消去法可以难过一般高斯法的这些缺点。 一、列主元高斯消去法解方程的Matlab程序如下： function a=columneli(a) %对矩阵a进行列主元消去 [m n]=size(a); %求取a的行数m和列数n for i=1:m-1 [maxEle,pos]=max(abs(a(i:end,i))); maxRow=pos+i-1; %在每次变换前寻找绝对值最大主所在列maxRow if a(maxRow,i)==0 disp(矩阵为奇异矩阵) return %对于非奇异矩阵，在程序中给予提示,结束程序 end if(maxRow~=i) temp=a(maxRow,:);a(maxRow,:)=a(i,:);a(i,:)=temp; end %与列主元绝对值最大的行进行行交换 for j=i+1:m a(j,i)=a(j,i)/a(i,i); %求取第j列主元 for k=i+1:n a(j,k)=a(j,k)-a(j,i)*a(i,k); %对第j列主元进行行变换 end end end function x=elisolve(a) %利用列主元消去的结果求方程的解,a为方程组的增广矩阵 a=columneli(a); [m,n]=size(a); x=zeros(m,1); for i=m:-1:1 x(i)=(a(i,n)-a(i,i:m)*x(i:m))/a(i,i); end 二、列主元高斯消去法解P227页第3题： 答案：x=[7/6 -1/3 1/2]T 三、程序流程图如下： 条件图框里面没有条件式，因为i是从m到1，所以i=1后运行下面的命令，直接输出结果。 这里也是一样的，直接运行下面的命令。 四、Matlab中的命令截图如下： 运行结果：x=[7/6 -1/3 1/2]T 上机结果和正确答案完全相同。 2. 列主元三角分解法 一，列主元三角分解法解方程的Matlab程序如下： function [l,u,p]=tridecom(a) %对矩阵a进行三角分解 [m n]=size(a); %求取a的行数m和列数n if(m~=n) disp(矩阵的行数与列数不相等，无法进行三角分解) return end l=zeros(n);u=zeros(n);p=eye(n); %p*a=l*u [maxEle,goodRow]=max(abs(a(:,1))); temp=a(1,:);a(1,:)=a(goodRow,:);a(goodRow,:)=temp; temp=p(1,:);p(1,:)=p(goodRow,:);p(goodRow,:)=temp; %行交换 u(1,:)=a(1,:);l(1:end,1)=a(1:end,1)/u(1,1); for r=2:n tempUrr=zeros(n,1); for k=r:n tempUrr(k)=a(k,r)-l(k,1:r-1)*u(1:r-1,r); end [maxEle,goodRow]=max(abs(tempUrr)); temp=a(r,:);a(r,:)=a(goodRow,:);a(goodRow,:)=temp; temp=l(r,:);l(r,:)=l(goodRow,:);l(goodRow,:)=temp; temp=p(r,:);p(r,:)=