北航数值分析第一次大作业（高斯gauss lu分解）

一、问题分析及算法描述

编写程序，分别用列主元的Gauss消去法和LU分解法求解下面线型代数方程组AX=b的解，其中A为N×N矩阵，N=50，其中第i(i≥1)行、第j(i≥1)列元素

1i+j?1

aij=

，

右端向量b的第i(i≥1)个分量为

10i+j?1

bi= Nj=1

.

列主元素Gauss消去过程中，要用到两种初等行变换。第一种，交换两行的位置；第二种，用一个数乘某一行加到另一行上。在第k次消元之前，先对增广矩阵 A(k),b(k) 作第一种行变换，使得aik中绝对值最大的元素交换到第k行的主对角线位置上，然后再使用第二种行变换进行消元。如此往复，最后得到一个上三角系数矩阵，并回代求解解向量。由于每次消元前选取了列主元素，因此与顺序Guass消元法相比，可提高数值计算的稳定性，且其计算量与顺序Guass消元法相同。列主元的Gauss消去法要求系数矩阵A非奇异。

(k)

LU分解法，即通过一系列初等行变换将系数矩阵A分解成一个下三角矩阵L与一个上三角矩阵U的乘积，进一步通过求解两个三角矩阵得出解向量。若L为单位下三角矩阵，U是上三角矩阵，则称为Doolittle分解；若L为下三角矩阵，U是单位上三角矩阵，则称为Crout分解。若系数矩阵A的前n-1阶顺序主子式不为零，则Doolittle\\Crout分解具有唯一性。若在每步行变换中选取主元，可提高数值计算稳定性。本算例中采用选主元的Doolittle分解。

通过分析可知，本算例中待求解线型方程组系数矩阵为非奇异矩阵，且其前n-1阶顺序主子式不为零。方程组的解向量为x= 10,10,?,10 T。满足列主元高斯消去法以及LU分解法的基本使用条件。为了验证上述两种方法对本算例的适用性，笔者利用Microsoft Visual C++6.0编写了该算例的列主元高斯消去法以及LU分解法的程序代码，并进行了运算求解。

二、计算结果及分析

1.高斯求解程序运行结果：

**********************************各行主元大小*********************************

第1行主元大小：1.00000e+000 第11行主元大小：-7.11093e-012 第2行主元大小：8.33333e-002 第12行主元大小：-2.06444e-013 第3行主元大小：6.48148e-003 第13行主元大小：1.41512e-014 第4行主元大小：1.13636e-003 第14行主元大小：-1.24316e-015 第5行主元大小：1.30273e-004 第15行主元大小：3.50355e-016 第6行主元大小：7.84929e-006 第16行主元大小：3.99615e-016 第7行主元大小：-5.19257e-007 第17行主元大小：3.02521e-016 第8行主元大小：3.74111e-008 第18行主元大小：6.93669e-017 第9行主元大小：-2.11897e-009 第19行主元大小：1.11654e-016 第10行主元大小：1.27171e-010

***********************************方程组的解**********************************

x1：9.999996 x11：108.556590 x2：10.000578 x12：53.247546 x3：9.977938 x13：-123.827690 x4：10.360343 x14：116.991522 x5：6.904310 x15：-164.520188 x6：25.281558 x16：291.762560 x7：-33.905604 x17：-166.575110 x8：76.093315 x18：-8.426663 x9：-10.091680 x19：73.980824 x10：-76.424166 x20：-9.385980 2.LU求解程序运行结果：

***********************************方程组的解**********************************

x1：9.999996 x11：108.556590 x2：10.000578 x12：53.247546 x3：9.977938 x13：-123.827690 x4：10.360343 x14：116.991522 x5：6.904310 x15：-164.520188 x6：25.281558 x16：291.762560 x7：-33.905604 x17：-166.575110 x8：76.093315 x18：-8.426663 x9：-10.091680 x19：73.980824 x10：-76.424166 x20：-9.385980 两种数值计算方法结果均与实际解相比产生了较大的误差，这是由于在消元过程中产生了绝对值很小的主元值，以列主元Guass消元法为例，最小主元值为6.93669e-017，这样就造成行乘数绝对值非常大，出现较小数加不到较大数中的现象，舍入误差的积累很大，运算结果完全失真。为了提高运算精度，笔者采用long double精度，但误差减小程度不显著。

三、结论

本方程组在进行高斯消元过程中出现小列主元情况，因而是病态矩阵。单纯通过提高运算精度并不能有效改善其病态程度，可采用平衡方法、残差校正法对该病态矩阵进行求解。

四、附录（VC++程序代码）

1.列主元Guass消去法：

#include #include #include

#define N 20

#define e 1e-20//主元绝对值下限

void main() { long double A[N][N]; long double b[N]; long double x[N]; long double m,S,temp,max; int i,j,k,col; //************************原始方程的计算和输出************************ //方程组原始条件的代入 for(i=0;i

printf(\高斯列主元消去法***********************\\n\printf(\//输出原方程组系数矩阵

printf(\原方程组系数矩阵A：\\n\for(i=0;i

}

}

if((j+1)%5==0) printf(\

//输出原方程组右端列向量

printf(\原方程组右端列向量b：\\n\for(i=0;i

//****************************消去过程********************************* printf(\各行主元大小***************************\\n\for(k=0;k

//查找最大元素所在的行 max=fabs(A[k][k]); col=k;

for(i=k;i

printf(\第%d行主元大小：%1.5e\\n\//判断主元大小 if(max<=e) { }

printf(\主元过小，造成舍入误差过大，高斯求解失效！\\n\exit(0);

//换行过程 for(j=k;j

temp=b[col];b[col]=b[k];b[k]=temp; //消元过程

for(i=k+1;i

}

}

for(j=0;j

//***************************消元结果********************************** printf(\消元结果*************************\\n\//高斯消元后新的系数矩阵

printf(\消元后系数矩阵A'：\\n\for(i=0;i

//高斯消元后新的右端列向量

printf(\原方程组右端列向量b'：\\n\for(i=0;i

printf(\第%d行：%1.5e\\n\ //*****************************回代过程******************************** x[N-1]=b[N-1]/A[N-1][N-1]; for(k=N-2;k>=0;k--) { S=b[k]; for(j=k+1;j

2.LU分解法： #include #include

#define N 20 void main() { long double a[N][N],l[N][N],u[N][N]; long double b[N],x[N],y[N],s[N],temp; long double max; int i,j,k,t; int M[N];

//*************************原始方程的计算和输出************************ //方程组原始条件的代入 for(i=0;i

printf(\printf(\选主元的Doolittle消去法*****************\\n\printf(\//输出原方程组系数矩阵

printf(\原方程组系数矩阵A：\\n\for(i=0;i

printf(\第%d列至第%d列：\\n\ printf(\ if((j+1)%5==0) printf(\ } }

//输出原方程组右端列向量

printf(\原方程组右端列向量b：\\n\for(i=0;i

//*****************************分解过程******************************** for(k=0;k

{ s[i]=a[i][k]; for(t=0;t<=k-1;t++) s[i]-=l[i][t]*u[t][k]; } max=0; for(i=k;i

//输出LU矩阵

printf(\矩阵：\\n\for(i=0;i

}

}

}

printf(\第%d列至第%d列：\\n\if(j=i) printf(\if((j+1)%5==0) printf(\

//求Qb

for(k=0;k

//求解y和x y[0]=b[0];

for(i=1;i

x[N-1]=y[N-1]/u[N-1][N-1]; for(i=N-2;i>=0;i--) { x[i]=y[i]; for(t=i+1;t

//输出结果

printf(\方程组的解**************************\\n\for(i=0;i

本文来源：https://www.bwwdw.com/article/3nbx.html