雅可比迭代法与矩阵特征值.docx
- 文档编号:12162695
- 上传时间:2023-06-04
- 格式:DOCX
- 页数:33
- 大小:284.19KB
雅可比迭代法与矩阵特征值.docx
《雅可比迭代法与矩阵特征值.docx》由会员分享,可在线阅读,更多相关《雅可比迭代法与矩阵特征值.docx(33页珍藏版)》请在冰点文库上搜索。
雅可比迭代法与矩阵特征值
实验五
矩阵的lu分解法,雅可比迭代法
班级:
学号:
姓名:
实验五矩阵的LU分解法,雅可比迭代
一、目的与要求:
Ø熟悉求解线性方程组的有关理论和方法;
Ø会编制列主元消去法、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,计算
程序与实例
程序设计如下:
#include
#include
usingnamespacestd;
voiddisp(double**p,introw,intcol){
for(inti=0;i for(intj=0;j cout< cout< } } voiddisp(double*q,intn){ cout<<"====================================="< for(inti=0;i cout<<"X["< cout<<"====================================="< } voidinput(double**p,introw,intcol){ for(inti=0;i cout<<"输入第"< "; for(intj=0;j cin>>p[i][j]; } } intfindMax(double**p,intstart,intend){ intmax=start; for(inti=start;i if(abs(p[i][start])>abs(p[max][start])) max=i; } returnmax; } voidswapRow(double**p,intone,intother,intcol){ doubletemp=0; for(inti=0;i temp=p[one][i]; p[one][i]=p[other][i]; p[other][i]=temp; } } booldispel(double**p,introw,intcol){ for(inti=0;i intflag=findMax(p,i,row);//找列主元行号 if(p[flag][i]==0)returnfalse; swapRow(p,i,flag,col);//交换行 for(intj=i+1;j doubleelem=p[j][i]/p[i][i];//消元因子 p[j][i]=0; for(intk=i+1;k p[j][k]-=(elem*p[i][k]); } } } returntrue; } doublesumRow(double**p,double*q,introw,intcol){ doublesum=0; for(inti=0;i if(i==row) continue; sum+=(q[i]*p[row][i]); } returnsum; } voidback(double**p,introw,intcol,double*q){ for(inti=row-1;i>=0;i--){ q[i]=(p[i][col-1]-sumRow(p,q,i,col))/p[i][i]; } } intmain() { cout<<"Inputn: "; intn;//方阵的大小 cin>>n; double**p=newdouble*[n]; for(inti=0;i p[i]=newdouble[n+1]; } input(p,n,n+1); if(! dispel(p,n,n+1)){ cout<<"奇异"< return0; } double*q=newdouble[n]; for(inti=0;i q[i]=0; back(p,n,n+1,q); disp(q,n); delete[]q; for(inti=0;i delete[]p[i]; delete[]p; } 1.用列主元消去法解方程 例2解方程组 计算结果如下 B=-1.461954 C=1.458125 D=-6.004824 E=-2.209018 F=14.719421 Ø矩阵直接三角分解法 算法: 将方程组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 usingnamespacestd; double**newMatrix(introw,intcol){ double**p=newdouble*[row];//行 for(inti=0;i p[i]=newdouble[col]; returnp; } voiddelMatrix(double**p,introw){ for(inti=0;i delete[]p[i]; delete[]p; } voidinputMatrix(double**p,introw,intcol){ for(inti=0;i cout<<"输入第"< "; for(intj=0;j cin>>p[i][j]; } } voiddispMatrix(double**p,introw,intcol){ for(inti=0;i for(intj=0;j cout< cout< } } voiddispVector(double*q,intn){ cout<<"====================================="< for(inti=0;i cout<<"X["< cout<<"====================================="< } voidinitMatrix(double**p,introw,intcol){ for(inti=0;i for(intj=0;j p[i][j]=0; } doublesum(double**L,double**U,introw,intcol){ intmin=(row>=col? col: row); doublesum=0; for(inti=0;i sum+=(U[i][col]*L[row][i]); returnsum; } voidresolve(double**A,double**L,double**U,introw,intcol){ for(inti=0;i U[0][i]=A[0][i];//把A的第一行给U的第一行 L[0][0]=A[0][0]; for(inti=1;i L[i][0]=A[i][0]/A[0][0]; for(inti=1;i for(intj=1;j if(i<=j) U[i][j]=A[i][j]-sum(L,U,i,j); else L[i][j]=(A[i][j]-sum(L,U,i,j))/U[j][j]; } for(inti=0;i L[i][i]=1; } doublesumRowY(double**L,double*y,introw){ doublesum=0; for(inti=0;i sum+=(L[row][i]*y[i]); returnsum; } voidbackY(double**L,double*b,double*y,introw){ for(inti=0;i y[i]=0;//初始化y向量全为零 for(inti=0;i y[i]=b[i]-sumRowY(L,y,i); } doublesumRowX(double**U,double*x,introw,intcol){ doublesum=0; for(inti=row+1;i sum+=(U[row][i]*x[i]); returnsum; } voidbackX(double**U,double*y,double*x,introw){ for(inti=0;i x[i]=0;//初始化x向量全为零 for(inti=row-1;i>=0;i--) x[i]=(y[i]-sumRowX(U,x,i,row))/U[i][i]; } intmain() { cout<<"输入方阵大小n: "; intn=0; cin>>n; cout<<"====================================="< cout<<"开始录入方阵数据....."< double**A=newMatrix(n,n);//开辟矩阵A inputMatrix(A,n,n);//录入数据到A中 cout<<"录入方阵数据完毕......"< cout<<"====================================="< cout<<"开始录入列阵....."< cout<<"输入列阵: "; double*b=newdouble[n];//开辟向量b for(inti=0;i cin>>b[i];//录入数据到b中 cout<<"录入列阵数据完毕....."< double**L=newMatrix(n,n);//开辟方阵L initMatrix(L,n,n);//初始化L全为零 double**U=newMatrix(n,n);//开辟方阵U initMatrix(U,n,n);//初始化U resolve(A,L,U,n,n);//分解A为L与U double*y=newdouble[n];//开辟向量y backY(L,b,y,n);//回代求出y double*x=newdouble[n];//开辟向量X backX(U,y,x,n);//回代求出X dispVector(x,n); //释放空间 delMatrix(A,n); delMatrix(L,n); delMatrix(U,n); delete[]b; delete[]y; delete[]x; } Ø迭代法 雅可比迭代法 算法: 设方程组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 usingnamespacestd; double**newMatrix(introw,intcol){ double**p=newdouble*[row];//行 for(inti=0;i p[i]=newdouble[col]; returnp; } voiddelMatrix(double**p,introw){ for(inti=0;i delete[]p[i]; delete[]p; } voidinputMatrix(double**p,introw,intcol){ for(inti=0;i cout<<"输入第"< "; for(intj=0;j cin>>p[i][j]; } } voiddispMatrix(double**p,introw,intcol){ for(inti=0;i for(intj=0;j cout< cout< } } voiddispVector(double*q,intn){ for(inti=0;i cout<<"X["< cout< } doublesumRow(double**A,double*x1,introw,intcol){ doublesum=0; for(inti=0;i if(row==i) continue; sum+=(A[row][i]*x1[i]); } returnsum; } voiditerat(double**A,double*b,double*x1,double*x2,introw){ for(inti=0;i x2[i]=(b[i]-sumRow(A,x1,i,row))/A[i][i]; } boolblow_error(double*x1,double*x2,introw,doublee){ doublesum=0; for(inti=0;i sum+=abs(x2[i]-x1[i]); } if(sum elsereturnfalse; } intmain() { cout<<"输入方阵大小n: "; intn=0; cin>>n; cout<<"====================================="< cout<<"开始录入方阵数据....."< double**A=newMatrix(n,n);//开辟矩阵A inputMatrix(A,n,n);//录入数据到A中 cout<<"录入方阵数据完毕......"< cout<<"====================================="< cout<<"开始录入列阵....."< cout<<"输入列阵: "; double*b=newdouble[n];//开辟向量b for(inti=0;i cin>>b[i];//录入数据到b中 cout<<"录入列阵数据完毕....."< cout<<"====================================="< cout<<"输入迭代次数容许的最大值: "; intM=0; cin>>M; cout<<"输入容许误差: "; doublee=0; cin>>e; cout<<"====================================="< double*x1=newdouble[n];//开辟x1向量 double*x2=newdouble[n];//开辟x2向量 cout<<"输入初始向量: "; for(inti=0;i cin>>x1[i]; cout<<"====================================="< intk=0;//迭代计数器 for(;;){ iterat(A,b,x1,x2,n);//迭代一次 if(blow_error(x1,x2,n,e)){ dispVector(x2,n); cout<<"迭代次数为: "< return0; }else{ if(k>=M){ cout<<"不收敛"< return0; }else{ k++; for(inti=0;i x1[i]=x2[i]; } } } //释放空间 delMatrix(A,n); delete[]b; delete[]x1; delete[]x2; } Ø高斯-塞尔德迭代法 算法: 设方程组Ax=b的系数矩阵的对角线元素, M为迭代次数容许的最大值,ε为容许误差 ①取初始向量 令k=0。 ②对i=1,2,…,n计算 ③如果 ,则输出 结束;否则执行④。 ④如果 则不收敛,终止程序;否则 转②。 程序与实例 例5用高斯-塞尔德迭代法解方程组 结果为 X[0]=3.000000 X[1]=2.000000 X[2]=1.000000 #include #include #include #include #defineN100 main() { inti; float*x; floatc[12]={8.0,-3.0,2.0,20.0, 4.0,11.0,-1.0,33.0, 6.0,3.0,12.0,36.0}; float*GauseSeidel(float*,int); x=GauseSeidel(c,3); for(i=0;i<=2;i++)printf("x[%d]=%f\n",i,x[i]); getch(); } float*GauseSeidel(float*a,intn) { inti,j,nu=0; float*x,dx,d; x=(float*)malloc(n*sizeof(float)); for(i=0;i<=n-1;i++)x[i]=0.0; do { for(i=0;i<=n-1;i++) { d=0.0; for(j=0;j<=n-1;j++) d+=*(a+i*(n+1)+j)*x[j]; dx=(*(a+i*(n+1)+n)-d)/(*(a+i*(n+1)+i)); x[i]+=dx; } if(nu>=N) { printf("foldnumber\n"); } nu++; } while(fabs(dx)>1e-6); returnx; } 例6用雅可比迭代法解方程组 迭代4次得解 若用高斯-塞尔德迭代法则发散。 #include voidmain(void) { intk,n; doublex[3]={7,2,5}; for(k=0;k<5;k++) { doublea,b; a=x[0]; b=x[1]; x[0]=(7-2.0*x[1]+2*x[2])/1; x[1]=(2-a-x[2])/1; x[2]=(5-2*a-2*b)/1; } for(n=0;n<3;n++) { printf("x[%d]=%8.6f\n",n,x[n]); } } 用高斯-塞尔德迭代法解方程组 迭代84次得解 ,若用雅克比迭代法则发散。 #include #include voidLOOP(floata[10][10],floatb[10],floatx[10],int); voidmain(void) {
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 可比 迭代法 矩阵 特征值