实验五矩阵的LU分解法，雅可比迭代
实
验
报
告
学院：计算机科学与软件学院班级：116班
姓名：薛捷星
学号：112547
一、目的与要求：
熟悉求解线性方程组的有关理论和方法；
会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序；
通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法。 二、 实验内容：
会编制列主元消去法、LU 分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解各种方法的优缺点。
三、 程序与实例
列主元高斯消去法
算法：将方程用增广矩阵[A ∣b ]=(ij a )1n (n )+?表示 1) 消元过程 对k=1,2,…,n-1
①选主元，找{}n ,,1k ,k i k +∈使得
k
,i k a =ik a n
i k max ≤≤ ②如果0a k
,i
k =，则矩阵A 奇异，程序结束；否则执行③。
③如果k i k ≠，则交换第k 行与第k i 行对应元素位置，
j i kj k a a ? j=k,┅,n+1
④消元，对i=k+1, ┅,n 计算
kk ik ik a a l /=
对j=l+1, ┅,n+1计算
kj ik ij ij a l a a -=
2) 回代过程
①若0=nn a ，则矩阵A 奇异，程序结束；否则执行②。 ②nn n n n a a x /1,+=；对i=n-1, ┅,2,1,计算
ii n
i j j ij n i i a x a a x /11,???
? ?
?-
=∑+=+
程序与实例 例1 解方程组
??
?
??=++-=++-=++035
.3643x .5072x .1835x .2137.2623x .43712x 347x .1 1.183
3.555x 2.304x 0.101x 321321321 输出结果如下： X[0]=-0.398234 X[1]= 0.013795
X[2]= 0.335144
程序如下：
#include
int i,j,p,o,l,q; double
a[3][4]={{0.101,2.304,3.555,1.183},{-1.347,3.712,4.623,2.137},{-2.835,1.072,5.643,3.035}};
double x[3],z[4];
printf("列主元消去法\n"); for(j=0;j<2;j++) { for(i=j+1;i<3;i++) { if(fabs(a[j][j]) { for(p=0;p<4;p++) { z[p]=a[j][p]; a[j][p]=a[i][p]; a[i][p]=z[p]; }/*交换得最大主元*/ } } for(l=j+1;l<3;l++) { for(q=3;q>=j;q--) { a[l][q]=(a[l][q]-(a[l][j]/a[j][j])*a[j][q]); } } printf("进行消去：\n"); for(o=0;o<3;o++) { for(p=0;p<4;p++) { printf("%12.6f",a[o][p]); } printf("\n"); } } x[2]=a[2][3]/a[2][2]; x[1]=(a[1][3]-x[2]*a[1][2])/a[1][1]; x[0]=(a[0][3]-x[2]*a[0][2]-x[1]*a[0][1])/a[0][0]; printf("最后的解:\n"); for(i=0;i<3;i++) { printf("x[%d]=%12.6f\n",i,x[i]); } } 结果如下： 例2 解方程组 ????? ????=++++=++++=++++=++++=++++-12.04 1.0F 1.02E 3.47D 1.04C 3.54B -6.301.0F 2.01E 2.51D 4.04C 5.05B -8.531.0F 1.21E 2.92D 1.46C 3.53B -20.071.0F 1.10E 4.48D 1.21C 4.93B -32.041.0F 1.55E 5.66D 2.40C 8.77B 计算结果如下 B=-1.161954 C= 1.458125 D=- 6.004824 E=-2.209018 F= 14.719421 程序如下： #include int i,j,p,o,l,q; double a[5][6]={{8.77,2.40,5.66,1.55,1.0,-32.04},{4.93,1.21,4.48,1.10,1.0,-20.07},{3.53,1.46,2.92,1.21,1.0,-8.53},{5.05,4.04,2.51,2.01,1.0,-6.30},{3.54,1.04,3.47,1.02,1.0,-12.04}}; double x[5],z[6]; printf("列主元消去法求五元一次方程组：\n"); for(j=0;j<4;j++) { for(i=j+1;i<5;i++) { if(fabs(a[j][j]) { for(p=0;p<6;p++) { z[p]=a[j][p]; a[j][p]=a[i][p]; a[i][p]=z[p]; }/*交换得最大主元*/ } } for(l=j+1;l<5;l++) { for(q=5;q>=j;q--) { a[l][q]=(a[l][q]-(a[l][j]/a[j][j])*a[j][q]); } } printf("消去一列：\n"); for(o=0;o<5;o++) { for(p=0;p<6;p++) { printf("%12.6f",a[o][p]); } printf("\n"); } } x[4]=a[4][5]/a[4][4]; x[3]=(a[3][5]-x[4]*a[3][4])/a[3][3]; x[2]=(a[2][5]-x[4]*a[2][4]-x[3]*a[2][3])/a[2][2]; x[1]=(a[1][5]-x[4]*a[1][4]-x[3]*a[1][3]-x[2]*a[1][2])/a[1][1]; x[0]=(a[0][5]-x[4]*a[0][4]-x[3]*a[0][3]-x[2]*a[0][2]-x[1]*a[0][1])/a[0][0]; printf("方程组的解为：\n"); for(i=0;i<5;i++) { printf("x[%d]=%12.6f\n",i,x[i]); } } 矩阵直接三角分解法 算法：将方程组A x=b 中的A 分解为A =LU ，其中L 为单位下三角矩阵，U 为上三角矩阵，则方程组A x=b 化为解2个方程组Ly =b ，Ux =y ，具体算法如下： ①对j=1,2,3,…,n 计算 j j a u 11= 对i=2,3,…,n 计算 1111/a a l i i = ②对k=1,2,3,…,n: a. 对j=k,k+1,…,n 计算 ∑-=-=1 1 k q qj kq kj kj u l a u b. 对i=k+1,k+2,…,n 计算 kk k q qk iq ik ik u u l a l /)(1 1∑-=-= ③11b y =，对k=2,3,…,n 计算 ∑-=-=1 1k q q kq k k y l b y ④nn n n u y x /=,对k=n-1,n-2,…,2,1计算 kk n k q q kq k k u x u y x /)(1 ∑+=- = 注：由于计算u 的公式于计算y 的公式形式上一样，故可直接对增广矩阵 [A ∣b ]=? ????? ? ?????+++1,2 11,2222211,111211n n nn n n n n n n a a a a a a a a a a a a 施行算法②，③，此时U 的第n+1列元素即为y 。 程序与实例 例3 求解方程组A x=b A=????????? ???-----381265973274581221, b=???? ????????4911427 结果为 X[0]= 3.000001 X[1]=-2.000001 X[2]= 1.000000 X[3]= 5.000000 程序如下： #include void main(void) { int i,j; double a[4][4]={{1,2,-12,8},{5,4,7,-2},{-3,7,9,5},{6,-12,-8,3}}; double l[4][4],b[4]={27,4,11,49},y[4],x[4]; printf("直接三角分解法求方程组的解：\n"); for(i=0;i<4;i++) { l[i][0]=a[i][0]; l[0][i]=a[0][i]; } l[1][1]=a[1][1]-l[1][0]*a[0][1];l[1][2]=a[1][2]-l[1][0]*a[0][2];l[1][3]=a[1][3] -l[1][0]*a[0][3]; l[2][1]=(a[2][1]-l[2][0]*l[0][1])/l[1][1];l[2][2]=a[2][2]-l[2][0]*l[0][2]-l[2][1] *l[1][2];l[2][3]=a[2][3]-l[2][0]*l[0][3]-l[2][1]*l[1][3]; l[3][1]=(a[3][1]-l[3][0]*l[0][1])/l[1][1];l[3][2]=(a[3][2]-l[3][0]*l[0][2]-l[3][1 ]*l[1][2])/l[2][2];l[3][3]=a[3][3]-l[3][0]*l[0][3]-l[3][1]*l[1][3]-l[3][2]*l[2][3]; printf("LU合并矩阵：\n"); for(i=0;i<4;i++) { for(j=0;j<4;j++) { printf("%12.6f",l[i][j]); } printf("\n"); } y[0]=b[0]; y[1]=b[1]-y[0]*l[1][0]; y[2]=b[2]-y[0]*l[2][0]-y[1]*l[2][1]; y[3]=b[3]-y[0]*l[3][0]-y[1]*l[3][1]-y[2]*l[3][2]; printf("Y矩阵：\n"); for(i=0;i<4;i++) printf("y[%d]=%12.6f\n",i,y[i]); x[3]=y[3]/l[3][3]; x[2]=(y[2]-x[3]*l[2][3])/l[2][2]; x[1]=(y[1]-x[3]*l[1][3]-x[2]*l[1][2])/l[1][1]; x[0]=(y[0]-x[3]*l[0][3]-x[2]*l[0][2]-x[1]*l[0][1])/l[0][0]; printf("方程组的解：\n"); for(i=0;i<4;i++) printf("x[%d]=%12.6f\n",i,x[i]); } 结果如下： 迭代法 雅可比迭代法 算法：设方程组Ax=b 系数矩阵的对角线元素)n ,,2,1i (0a ii =≠,M 为迭代次数容许的最大值,ε为容许误差。 ①取初始向量x =T (0)n )0(2)0(1)x ,,x ,x ( ，令k=0。 ②对i=1,2,…,n 计算 )x a b (a 1x n i j 1 j (k)j ij i ii )1k (i ∑≠=+-= ③如果ε<-∑=+n 1 i )k (i )1k (i x x ，则输出)1k (x +，结束；否则执行④。 ④如果k ≥M,则不收敛，终止程序；否则1+←k k ,转②。 程序与实例 例4 用雅可比迭代法解方程组 ??? ??=--=-+=++1 6321382825321 321321x x x x x x x x x 结果为 迭代次数为20 X[0]= 1.000000 X[1]= 2.000000 X[2]=-1.000000 程序如下： #include float a,b,c,x[3]; int i; printf("Jacobi 迭代法求方程组：\n"); printf("输入X1,X2,X3的初始值,以“，”间隔：\n"); scanf("%f,%f,%f",&a,&b,&c); for(i=0;;i++) { x[0]=(8-2*b-c)/5; x[1]=(21-2*a+3*c)/8; x[2]=(1-a+3*b)/(-6); if(fabs(x[0]-a) printf("迭代%d 次\n",i+3); printf("方程组的解为：\n"); for(i=0;i<3;i++) { printf("x[%d]=%f\n",i,x[i]); } } 结果如下： 高斯-塞尔德迭代法 算法:设方程组A x=b 的系数矩阵的对角线元素,),,2,1(0n i a ii =≠,M 为迭代次数容许的最大值,ε为容许误差 ①取初始向量T n x x x x ),,,()0()0(2)0(1 =,令k=0。 ②对i=1,2,…,n 计算 )a b (a 1x 1 )(11)1(ij i ii )1k (i ∑∑+=-=++--=n i j k j ij i j k j x a x ③如果ε<-∑=+n i k i k i x x 1 )()1(，则输出)1(+k x ,结束；否则执行④。 ④如果,M k ≥则不收敛，终止程序；否则1+←k k ,转②。 程序与实例 例5 用高斯-塞尔德迭代法解方程组 ??? ??=++=-+=+-36 12363311420 238321 321321x x x x x x x x x 结果为 X[0]=3.000000 X[1]=2.000000 X[2]=1.000000 程序如下： #include void main(void) { int i; float a[3][4]={{8,-3,2,20},{4,11,-1,33},{6,3,12,36}},x[3]; float t,b,c; printf("高斯--赛德尔法解方程组\n"); printf("输入X1，X2，X3的初始值，以逗号间隔：\n"); scanf("%f,%f,%f",&x[0],&x[1],&x[2]); for(i=0;;i++) { t=x[0];b=x[1];c=x[2]; x[0]=(a[0][3]-a[0][1]*x[1]-a[0][2]*x[2])/a[0][0]; x[1]=(a[1][3]-a[1][0]*x[0]-a[1][2]*x[2])/a[1][1]; x[2]=(a[2][3]-a[2][0]*x[0]-a[2][1]*x[1])/a[2][2]; if(fabs(x[0]-t) else continue; if(i>200) { printf("发散\n"); break; } } printf("迭代%d次\n",i+1); for(i=0;i<3;i++) printf("x[%d]=%12.6f\n",i,x[i]); } 结果如下： 例6 用雅可比迭代法解方程组 ??? ??=++=++=-+5 222722321 321321x x x x x x x x x 迭代4次得解T )1,2,1(-,若用高斯-塞尔德迭代法则发散。 结果如下： 用高斯-塞尔德迭代法解方程组 ?? ? ??=++=++=++7 .1x 9x .09x .00.29x .0x 9x .09 .19x .09x .0x 321321321 迭代84次得解()T 3,2,1，若用雅克比迭代法则发散。 结果如下：
相关主题
文本预览