高斯消元法和三角分解法

《高斯消元法和三角分解法》由会员分享，可在线阅读，更多相关《高斯消元法和三角分解法（8页珍藏版）》请在装配图网上搜索。

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=leng

2、th(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)max1max1=abs(A(i,k);r=i;endendif max1v1e-10flag=failure;return;end%交换两行if rkfor 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

3、=-det;end%消元过程for i=k+l:nm=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)v1e-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三角分解法：function L,U,flag=LU_Decom(A

4、)%求矩阵A的LU分解%A为被分解的矩阵%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)veps flag=failure;return;endfor i

5、=k+1:nz=0;for q=l:k-lz=z+L(i,q)*U(q,k);endL(i,k)=(A(i,k)-z)/U(k,k);endendP108 15(a)解：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(hilb(4),b)A=1.000000000000000.500000000000000.333333333333330.250000000000000.500000000000

6、000.333333333333330.250000000000000.200000000000000.333333333333330.250000000000000.200000000000000.166666666666670.250000000000000.200000000000000.166666666666670.14285714285714b =10 0 01.0e+002 *0.15999999999999-1.199999999999922.39999999999980-1.39999999999987det =1.653439153439300e-007 flag =OK2

7、、利用三角消元，在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=1.000000000000000.500000000000000.333333333333330.250000000000000.500000000000000.333333333333330.250000000000000.200000000000000.333333333333330.250000000000000.200000000000

8、000.166666666666670.250000000000000.200000000000000.166666666666670.14285714285714b =1000L =1.000000000000000000.500000000000001.000000000000001.0000000000000000.333333333333331.00000000000000000.250000000000000.90000000000000-0.600000000000001.00000000000000U =1.000000000000000.500000000000000.3333

9、33333333330.2500000000000000.083333333333330.088888888888890.0833333333333300-0.00555555555556-0.0083333333333300.00035714285714x =1.0e+002 *0.15999999999999-1.199999999999922.39999999999980-1.39999999999987由上面的输出结果可以看出：通过高斯消元和三角分解所的出的值二者之间的相差很小，甚至从结果无法判断二者的相差多少。(b)1、利用高斯消元，在matlab的命令窗口中输入命令，并得出结果：

10、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.16670.25000.20000.16670.1429b =10 0 0 x =16.0000-120.0000240.0000-140.0000 det =1.6534e-007flag =OK2、利用三角消元，在matlab的命令窗口中输入命令，并得出结果: fo

11、rmat 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 l/7,b=l 0 0 O,L,U=lu(A),x=U(Lb)A=1.00000.50000.33330.25000.50000.33330.25000.20000.33330.25000.20000.16670.25000.20000.16670.1429b =10001.00000000.50000.0833000.33330.0889-0.005600.25000.0833-0.00830.00041.00000000.50001.00001.000000.3

12、3331.0000000.25000.9000-0.60001.000016.0000-120.0000240.0000-140.0000分析：通过三角分解求解线性方程组，当矩阵为四位有效数字表示时和用分数表示时，其 输出的解结果明显不同，当二者的差值在10-4到10-5时输出结果相差值有1点多(a、b中 的三角分解的输出结果)。因此可以看出当矩阵系数发生明显变化其输出结果变化也会明 显。其稳定性不好。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,f

13、lag=Gauss(A,b)13572-1350025-2-6-31b =123x =1.34290.6857-3.00001.8000 det =35.0000 flag =OK 2、利用三角消元，在matlab的命令窗口中输入命令，并得出结果： A=1 3 5 7;2 -1 3 5;0 0 2 5;-2 -6 -3 l,b=l 2 3 4,L,U=lu(A),x=U(Lb)A=13572-1350025-2-6-31b =12340.5000-0.50001.00001.000000000.5714-1.00001.000001.00000U =2.0000-1.00003.00005.00000-7.000006.0000003.50007.50000000.71431.34290.6857-3.00001.8000