解线性方程组的列主元素高斯消去法和LU分解法实验报告.docx
- 文档编号:6789131
- 上传时间:2023-01-10
- 格式:DOCX
- 页数:9
- 大小:32.61KB
解线性方程组的列主元素高斯消去法和LU分解法实验报告.docx
《解线性方程组的列主元素高斯消去法和LU分解法实验报告.docx》由会员分享,可在线阅读,更多相关《解线性方程组的列主元素高斯消去法和LU分解法实验报告.docx(9页珍藏版)》请在冰豆网上搜索。
解线性方程组的列主元素高斯消去法和LU分解法实验报告
解线性方程组的列主元素高斯消去法和LU分解法
一、实验目的:
通过数值实验,从中体会解线性方程组选主元的必要性和LU分解法的优点,以及方程组系数矩阵和右端向量的微小变化对解向量的影响。
二、实验内容:
解下列两个线性方程组
三、实验要求:
(1)用你熟悉的算法语言编写程序用列主元高斯消去法和LU分解求解上述两个方程组,输出Ax=b中矩阵A及向量b,A=LU分解的L及U,detA及解向量x.
(2)将方程组
中系数3.01改为3.00,0.987改为0.990,用列主元高斯消去法求解变换后的方程组,输出列主元行交换次序,解向量x及detA,并与
中结果比较。
(3)将方程组
中的2.099999改为2.1,5.900001改为5.9,用列主元高斯消去法求解变换后的方程组,输出解向量x及detA,并与
中的结果比较。
(4)用MATLAB的内部函数inv求出系数矩阵的逆矩阵,再输入命令x=inv(A)*b,即可求出上述各个方程组的解,并与列主元高斯消去法和LU分解法求出的解进行比较,体会选主元的方法具有良好的数值稳定性。
用MATLAB的内部函数det求出系数行列式的值,并与
(1)、
(2)、(3)中输出的系数行列式的值进行比较。
四、实验过程:
(1)列主元高斯消去法的主程序为
function[RA,RB,n,X]=liezhuY(A,b)
B=[Ab];n=length(b);RA=rank(A);
RB=rank(B);zhica=RB-RA;
D=det(A)
ifzhica>0,
disp('请注意:
因为RA~=RB,所以此方程组无解.')
return
end
ifRA==RB
ifRA==n
disp('请注意:
因为RA=RB=n,所以此方程组有唯一解.')
X=zeros(n,1);C=zeros(1,n+1);
forp=1:
n-1
[Y,j]=max(abs(B(p:
n,p)));C=B(p,:
);
B(p,:
)=B(j+p-1,:
);B(j+p-1,:
)=C;
fork=p+1:
n
m=B(k,p)/B(p,p);
B(k,p:
n+1)=B(k,p:
n+1)-m*B(p,p:
n+1);
end
end
b=B(1:
n,n+1);A=B(1:
n,1:
n);X(n)=b(n)/A(n,n);
forq=n-1:
-1:
1
X(q)=(b(q)-sum(A(q,q+1:
n)*X(q+1:
n)))/A(q,q);
end
else
disp('请注意:
因为RA=RB end end 解方程组 (1) 在MATLAB工作窗口输入 >>A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];b=[1;1;1];[RA,RB,n,X]=liezhuY(A,b) 运行后输出结果为 请注意: 因为RA=RB=n,所以此方程组有唯一解. D=-0.1225 RA=3RB=3n=3 X=397.8654 -157.6242 -123.1120 解方程组 (2) 在MATLAB工作窗口输入 >>A=[10-701;-32.09999962;5-15-1;2102];b=[8;5.900001;5;1];[RA,RB,n,X]=liezhu(A,b) 运行后输出结果为 请注意: 因为RA=RB=n,所以此方程组有唯一解. D=-762.0000 RA=4RB=4n=4 X=0.0000 -1.0000 1.0000 1.0000 LU分解法及MATLAB主程序为 functionhl=zhjLU(A) [nn]=size(A);RA=rank(A); D=det(A) ifRA~=n disp('请注意: 因为A的n阶行列式hl等于零,所以A不能进行LU分解.A的秩RA如下: '),RA,hl=det(A); return end ifRA==n forp=1: n h(p)=det(A(1: p,1: p)); end hl=h(1: n); fori=1: n ifh(1,i)==0 disp('请注意: 因为A的r阶主子式等于零,所以A不能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下: '),hl;RA return end end ifh(1,i)~=0 disp('请注意: 因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下: ') forj=1: n U(1,j)=A(1,j); end fork=2: n fori=2: n forj=2: n L(1,1)=1;L(i,i)=1; ifi>j L(1,1)=1;L(2,1)=A(2,1)/U(1,1);L(i,1)=A(i,1)/U(1,1); L(i,k)=(A(i,k)-L(i,1: k-1)*U(1: k-1,k))/U(k,k); else U(k,j)=A(k,j)-L(k,1: k-1)*U(1: k-1,j); end end end end hl;RA,U,L end end 解方程组 (1) 在MATLAB工作窗口输入 >>A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];h1=zhjLU(A) 运行输出结果为 请注意: 因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下: D=9.8547 RA=3 U=3.01006.03001.9990 04.1600-2.0734 005.3016 L=1.000000 0.42191.00000 0.3279-1.63161.0000 h1=3.01004.8635-0.1225 解方程组 (2) 在MATLAB工作窗口输入 >>A=[10-701;-32.09999962;5-15-1;2102];h1=zhjLU(A) 运行后输出结果为 请注意: 因为A的各阶主子式都不等于零,所以A能进行LU分解.A的秩RA和各阶顺序主子式值hl依次如下: D=-762.0000 RA=4 U=10.0000-7.000001.0000 02.10006.00002.3000 00-2.1429-4.2381 0-0.0000012.7333 L=1.0000000 -0.30001.000000 0.50001.19051.0000-0.0000 0.20001.14293.20001.0000 h1=10.0000-0.0000-150.0001-762.0001 (2)在MATLAB工作窗口输入 >>A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];b=[1;1;1];A(1,1)=3;A(1,3)=0.990;[RA,RB,n,X]=liezhu(A,b) 请注意: 因为RA=RB=n,所以此方程组有唯一解. RA=3RB=3n=3 X=-4.0264 1.9193 1.5210 hi=3.00004.82199.8547 在MATLAB工作窗口输入 x=[397.8654;-157.6242;-123.1120]';x1=[-4.0264;1.9193;1.5210]';wucha=x1-x 运行后输出结果为 wucha=-401.8918159.5435124.6330 (3)在MATLAB工作窗口输入 >>A=[10-701;-32.09999962;5-15-1;2102];A(2,2)=2.1;b(2,1)=5.9;b=[8;5.900001;5;1];[RA,RB,n,X]=liezhu(A,b) 运行后输出结果为 请注意: 因为RA=RB=n,所以此方程组有唯一解. RA=4RB=4n=4 X=0.0000 -1.0000 1.0000 1.0000 h1=10.0000-0.0000-150.0000-762.0000 在MATLAB工作窗口输入 >>x=[0;-1;1;1]';x1=[0;-1;1;1]';wucha=x1-x 运行后输出结果为 wucha=0000 (4)解方程组 (1) 在MATLAB工作窗口输入 >>A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];B=inv(A) 运行后结果为 B=-268.9293538.3418128.4529 106.7599-213.4281-50.9561 83.3992-166.8022-39.7090 在MATLAB工作窗口输入 >>b=[1;1;1];x=inv(A)*b 运行后结果为 x=397.8654 -157.6242 -123.1120 在MATLAB工作窗口输入 >>A=[3.016.031.999;1.274.16-1.23;0.987-4.819.34];A(1,1)=3;A(1,3)=0.990;B=inv(A) 运行输出结果为 B=3.3424-6.1983-1.1705 -1.32692.74420.5020 -1.03652.06820.4893 在MATLAB工作窗口输入 >>b=[1;1;1];x=inv(A)*b 运行后输出结果为 x=-4.0264 1.9193 1.5210 解方程组 (2) 在MATLAB工作窗口输入 >>A=[10-701;-32.09999962;5-15-1;2102];B=inv(A) 运行后结果为 B=-0.0223-0.09840.11810.1686 -0.1601-0.11810.14170.2690 0.01080.10630.0724-0.0755 0.10240.1575-0.18900.1969 在MATLAB工作窗口输入 >>b=[8;5.900001;5;1];x=inv(A)*b 运行后输出结果为 x=0 -1.0000 1.0000 1.0000 在MATLAB工作窗口输入 >>A=[10-701;-32.09999962;5-15-1;2102];A(2,2)=2.1;B=inv(A) 运行后输出结果为 B=-0.0223-0.09840.11810.1686 -0.1601-0.11810.14170.2690 0.01080.10630.0724-0.0755 0.10240.1575-0.18900.1969 在MATLAB工作窗口输入 >>b=[8;5.900001;5;1];b(2,1)=5.9;x=inv(A)*b 运行后输出结果为 x=-0.0000 -1.0000 1.0000 1.0000 五、实验结果分析: 实验的数学原理很容易理解,也容易上手。 把运算的结果带入原方程组,可以发现符合的还是比较好。 这说明列主元消去法计算这类方程的有效性。 当A可逆时,能够将计算进行到底,列主元法就能确保算法的稳定,而且计算量不大。 直接三角消去过程,实质上是将A分解为两个三角矩阵的乘积A=LU,并求解Ly=b的过程。 回带过程就是求解上三角方程组Ux=y。 所以在实际的运算中,矩阵L和U可以直接计算出,而不需要任何中间步骤,从而在计算过程中将高斯消去法的步骤进行了进一步的简略,大大提高了运算速度,这就是三角分解法。 通过以上的计算比较,方程组 (1)具有严重的病态性。 当系数矩阵有微小的变化时,wucha=-401.8918159.5435124.6330,所得的解与原方程组的解有很大的相对误差。 方程组 (2)中当系数矩阵A和b有微小变化时,wucha=0000,所得的解与方程组的解没有相对误差。 所以方程组 (2)是良性的。 用MATLAB内部函数inv通过求逆矩阵,然后通过x=inv(A)*b也可以求出方程组的解,但是没有列主元高斯消去法具有良好的稳定性。 det函数求方程组系数矩阵的行列式时所得结果和高斯消去法和三角法所得结果相同,具有方便快捷的优点。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 线性方程组 元素 消去 LU 解法 实验 报告