实验五矩阵的LU分解法雅可比迭代.docx
- 文档编号:3553407
- 上传时间:2022-11-23
- 格式:DOCX
- 页数:16
- 大小:357.47KB
实验五矩阵的LU分解法雅可比迭代.docx
《实验五矩阵的LU分解法雅可比迭代.docx》由会员分享,可在线阅读,更多相关《实验五矩阵的LU分解法雅可比迭代.docx(16页珍藏版)》请在冰豆网上搜索。
实验五矩阵的LU分解法雅可比迭代
实验五矩阵的LU分解法,雅可比迭代
实
验
报
告
学院:
计算机科学与软件学院
班级:
116班
姓名:
薛捷星
学号:
112547
一、目的与要求:
Ø熟悉求解线性方程组的有关理论和方法;
Ø会编制列主元消去法、LU分解法、雅可比及高斯—塞德尔迭代法德程序;
Ø通过实际计算,进一步了解各种方法的优缺点,选择合适的数值方法。
二、实验内容:
Ø会编制列主元消去法、LU分解法、雅可比及高斯—塞德尔迭代法德程序,进一步了解各种方法的优缺点。
三、程序与实例
Ø列主元高斯消去法
算法:
将方程用增广矩阵[A∣b]=(
表示
1)消元过程
对k=1,2,…,n-1
①选主元,找
使得
=
②如果
,则矩阵A奇异,程序结束;否则执行③。
③如果
,则交换第k行与第
行对应元素位置,
j=k,┅,n+1
④消元,对i=k+1,┅,n计算
对j=l+1,┅,n+1计算
2)回代过程
①若
,则矩阵A奇异,程序结束;否则执行②。
②
;对i=n-1,┅,2,1,计算
程序与实例
例1解方程组
输出结果如下:
X[0]=-0.398234
X[1]=0.013795
X[2]=0.335144
程序如下:
#include
#include
main()
{
inti,j,p,o,l,q;
doublea[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}};
doublex[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解方程组 计算结果如下 B=-1.161954 C=1.458125 D=-6.004824 E=-2.209018 F=14.719421 程序如下: #include #include voidmain(void) { inti,j,p,o,l,q; doublea[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}}; doublex[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]); } } Ø矩阵直接三角分解法 算法: 将方程组Ax=b中的A分解为A=LU,其中L为单位下三角矩阵,U为上三角矩阵,则方程组Ax=b化为解2个方程组Ly=b,Ux=y,具体算法如下: ①对j=1,2,3,…,n计算 对i=2,3,…,n计算 ②对k=1,2,3,…,n: a.对j=k,k+1,…,n计算 b.对i=k+1,k+2,…,n计算 ③ ,对k=2,3,…,n计算 ④ 对k=n-1,n-2,…,2,1计算 注: 由于计算u的公式于计算y的公式形式上一样,故可直接对增广矩阵 [A∣b]= 施行算法②,③,此时U的第n+1列元素即为y。 程序与实例 例3求解方程组Ax=b A= b= 结果为 X[0]=3.000001 X[1]=-2.000001 X[2]=1.000000 X[3]=5.000000 程序如下: #include voidmain(void) { inti,j; doublea[4][4]={{1,2,-12,8},{5,4,7,-2},{-3,7,9,5},{6,-12,-8,3}}; doublel[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系数矩阵的对角线元素 M为迭代次数容许的最大值,ε为容许误差。 ①取初始向量x= ,令k=0。 ②对i=1,2,…,n计算 ③如果 ,则输出 ,结束;否则执行④。 ④如果k≥M,则不收敛,终止程序;否则 转②。 程序与实例 例4用雅可比迭代法解方程组 结果为 迭代次数为20 X[0]=1.000000 X[1]=2.000000 X[2]=-1.000000 程序如下: #include #include #definee0.000001 voidmain(void) { floata,b,c,x[3]; inti; 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) break; else a=x[0];b=x[1];c=x[2]; if(i>200) { printf("发散\n"); break; } } printf("迭代%d次\n",i+3); printf("方程组的解为: \n"); for(i=0;i<3;i++) { printf("x[%d]=%f\n",i,x[i]); } } 结果如下: Ø高斯-塞尔德迭代法 算法: 设方程组Ax=b的系数矩阵的对角线元素, M为迭代次数容许的最大值,ε为容许误差 ①取初始向量 令k=0。 ②对i=1,2,…,n计算 ③如果 ,则输出 结束;否则执行④。 ④如果 则不收敛,终止程序;否则 转②。 程序与实例 例5用高斯-塞尔德迭代法解方程组 结果为 X[0]=3.000000 X[1]=2.000000 X[2]=1.000000 程序如下: #include #include #definee0.000001 voidmain(void) { inti; floata[3][4]={{8,-3,2,20},{4,11,-1,33},{6,3,12,36}},x[3]; floatt,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) break; 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用雅可比迭代法解方程组 迭代4次得解 若用高斯-塞尔德迭代法则发散。 结果如下: 用高斯-塞尔德迭代法解方程组 迭代84次得解 ,若用雅克比迭代法则发散。 结果如下: _
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 实验 矩阵 LU 解法 可比
![提示](https://static.bdocx.com/images/bang_tan.gif)