(完整word版)高斯消元法和三角分解法

1、实验四 线性方程组的直接数值解法信息与计算科学金融崔振威201002034031一、实验目的：编程实现高斯消元法和三角分解法二、实验内容：p108.15、p120.1三、实验要求：1、分别对所给算例利用高斯消元法和三角分解法进行数值求解2、分析系数矩阵p108.15的算法稳定性主程序高斯消元法：function x,det,flag=Gauss(A,b)%求线性方程组的列主元Gauss消去法%A为方程组的系数矩阵%b为方程组的右端项%x为方程组的解%det为系数矩阵A的行列式的值%flag为指标向量，flag为failure表示计算失败，为 ok时表示计算成功n ,m=size(A); nb=

2、le ngth(b);%当方程组与列的维数不相等时停止计算，并输出出错信息if n=merror(方程组与右端项列的维数不相等，错误。 );return;end%赋初始值，计算flag=OK;det=1;x=zeros (n ,1);for k=1: n-1%选主元max1=0;for i=k: nif abs(A(i,k)max1 max1=abs(A(i,k);r=i;endendif max1kfor j=k: nz=A(k,j);A(k,j)=A(r,j);A(r,j)=z;endz=b(k);b(k)=b(r);b(r)=z;det=_det;end%消元过程for i=k+1: n

3、m=A(i,k)/A(k,k);for j=k+1: nA(i,j)=A(i,j)-m*A(k,j);endb(i)=b(i)-m*b(k);enddet=det*A(k,k);enddet=det*A( n,n);%回代过程if abs(A( n,n) )1e-10flag=failure;return;endfor k=n :-1:1for j=k+1: nb(k)=b(k)-A(k,j)*x(j);endx(k)=b(k)/A(k,k);endx(k)=b(k)/A(k,k);end三角分解法：fun ctio n L,U,flag=LU_Decom(A)%求矩阵A的LU分解%A为被分解

4、的矩阵%L为单位下三角阵%U为单位上三角阵%flag为指标向量，flag为failure表示计算失败，为ok时表示计算成功n, m=size(A);%要求所分解的矩阵是方阵，否则停止计算，并输出错误信息if n=merror(所分解的矩阵不是一个方阵);return;endL=eye( n);U=zeros (n );flag=OK;for k=1: nfor j=k: nz=0;for q=1:k-1z=z+L(k,q)*U(q,j);endU(k,j)=A(k,j)-z;endif abs(U(k,k) A=1 1/2 1/3 1/4;1/2 1/3 1/4 1/5;1/3 1/4 1/5

5、 1/6;1/4 1/5 1/6 1/7,b=1 0 0 0,x,det,flag=Gauss(hilb(4),b)A =1.000000000000000.500000000000000.333333333333330.250000000000000.500000000000000.333333333333330.250000000000000.200000000000000.333333333333330.250000000000000.200000000000000.166666666666670.250000000000000.200000000000000.1666666666666

6、70.14285714285714x =b =10 0 01.0e+002 *0.15999999999999-1.199999999999922.39999999999980-1.39999999999987det =1.653439153439300e-007 flag =OK2、利用三角消元，在matlab的命令窗口中输入命令，并得出结果： A=1 1/2 1/3 1/4;1/2 1/3 1/4 1/5;1/3 1/4 1/5 1/6;1/4 1/5 1/6 1/7,b=1 0 0 O,L,U=lu(A),x=U(Lb)A = b =1.000000000000000.500000000

7、000000.333333333333330.250000000000000.500000000000000.333333333333330.250000000000000.200000000000000.333333333333330.250000000000000.200000000000000.166666666666670.250000000000000.200000000000000.166666666666670.1428571428571410001.000000000000000000.500000000000001.000000000000001.00000000000000

8、00.333333333333331.00000000000000000.250000000000000.90000000000000-0.600000000000001.000000000000001.000000000000000.500000000000000.333333333333330.2500000000000000.083333333333330.088888888888890.0833333333333300-0.00555555555556-0.0083333333333300.000357142857141.0e+002 *0.15999999999999-1.19999

9、9999999922.39999999999980-1.39999999999987由上面的输出结果可以看出：通过高斯消元和三角分解所的出的值二者之间的相差很小，甚至从结果无法判断二者的相差多少。(b)1、利用高斯消元，在matlab的命令窗口中输入命令，并得出结果： A=1 1/2 1/3 1/4;1/2 1/3 1/4 1/5;1/3 1/4 1/5 1/6;1/4 1/5 1/6 1/7,b=1 0 0 0,x,det,flag=Gauss(A,b)A =1.00000.50000.33330.25000.50000.33330.25000.20000.33330.25000.20000

10、.16670.25000.20000.16670.1429b =10 0 0 x =16.0000-120.0000240.0000-140.0000 det =1.6534e-007flag =OK2、利用三角消元，在 matlab的命令窗口中输入命令，并得出结果: format A=1 1/2 1/3 1/4;1/2 1/3 1/4 1/5;1/3 1/4 1/5 1/6;1/4 1/5 1/6 1/7,b=1 0 0 O,L,U=lu(A),x=U(Lb)1.00000.50000.33330.25000.50000.33330.25000.20000.33330.25000.20000

11、.16670.25000.20000.16670.1429b =10001.0000 00.50001.00000.33331.00000.25000.90000 01.0000 00 0-0.6000 1.00001.00000000.50000.0833000.33330.0889-0.005600.25000.0833-0.00830.000416.0000-120.0000240.0000-140.0000分析：通过三角分解求解线性方程组，当矩阵为四位有效数字表示时和用分数表示时，其输出的解结果明显不同，当二者的差值在10-4到10-5时输出结果相差值有 1点多（a、b中的三角分解的输

12、出结果）。因此可以看出当矩阵系数发生明显变化其输出结果变化也会明 显。其稳定性不好。P120.1解：1、利用高斯消元，在 matlab的命令窗口中输入命令，并得出结果: A=1 3 5 7;2 -1 3 5;0 0 2 5;-2 -6 -3 1,b=1 2 3 4,x,det,flag=Gauss(A,b)A =13572-1350025-2-6-31 b =1234 x =1.34290.6857-3.00001.8000 det =35.0000 flag =OK2、利用三角消元，在matlab的命令窗口中输入命令，并得出结果： A=1 3 5 7;2 -1 3 5;0 0 2 5;-2 -6 -3 1,b=1 2 3 4,L,U