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

学号: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]=(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,:)=p(goodRow,:);p(goodRow,:)=temp;
for i=r:n
u(r,i)=a(r,i)-l(r,1:r-1)*u(1:r-1,i);
end
for i=r:n
l(i,r)=(a(i,r)-l(i,1:r-1)*u(1:r-1,r))/u(r,r);
end
end
function x=lusolve(a) %用三角分解法解线性方程组Ax=b,a为方程组的增广矩阵[A|b]
[m,n]=size(a);
[l,u,p]=(a(:,1:m)); %对A进行三角分解A=p*l*u
x=p*a(:,n); %对a的最后一列进行与方程组系数三角分解时相同的行变