雅可比迭代法及矩阵的特征值

-.z.实验五矩阵的lu分解法，雅可比迭代法班级：**：：实验五 矩阵的LU分解法，雅可比迭代一、目的与要求：熟悉求解线性方程组的有关理论和方法；会编制列主元消去法、LU分解法、雅可比及高斯—塞德尔迭代法德程序；通过实际计算，进一步了解各种方法的优缺点，选择适宜的数值方法。二、实验内容：会编制列主元消去法、LU分解法、雅可比及高斯—塞德尔迭代法德程序，进一步了解各种方法的优缺点。三、程序与实例列主元高斯消去法算法：将方程用增广矩阵[A∣b]=(表示消元过程对k=1,2,…,n-1①选主元，找使得=②如果，则矩阵A奇异，程序完毕；否则执行③。③如果，则交换第k行与第行对应元素位置，j=k,┅,n+1④消元，对i=k+1,┅,n计算对j=l+1,┅,n+1计算回代过程①假设，则矩阵A奇异，程序完毕；否则执行②。②；对i=n-1,┅,2,1,计算程序与实例程序设计如下：#include<iostream>#include<cmath>usingnamespacestd;voiddisp(double**p,introw,intcol){for(inti=0;i<row;i++){for(intj=0;j<col;j++)cout<<p[i][j]<<'';cout<<endl;}}voiddisp(double*q,intn){cout<<"====================================="<<endl;for(inti=0;i<n;i++)cout<<"*["<<i+1<<"]="<<q[i]<<endl;cout<<"====================================="<<endl;}voidinput(double**p,introw,intcol){for(inti=0;i<row;i++){cout<<"输入第"<<i+1<<"行:";for(intj=0;j<col;j++)cin>>p[i][j];}}intfindMa*(double**p,intstart,intend){intma*=start;for(inti=start;i<end;i++){if(abs(p[i][start])>abs(p[ma*][start]))ma*=i;}returnma*;}voidswapRow(double**p,intone,intother,intcol){doubletemp=0;for(inti=0;i<col;i++){temp=p[one][i];p[one][i]=p[other][i];p[other][i]=temp;}}booldispel(double**p,introw,intcol){for(inti=0;i<row;i++){intflag=findMa*(p,i,row);//找列主元行号if(p[flag][i]==0)returnfalse;swapRow(p,i,flag,col);//交换行for(intj=i+1;j<row;j++){doubleelem=p[j][i]/p[i][i];//消元因子p[j][i]=0;for(intk=i+1;k<col;k++){p[j][k]-=(elem*p[i][k]);}}}returntrue;}doublesumRow(double**p,double*q,introw,intcol){doublesum=0;for(inti=0;i<col-1;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<n;i++){p[i]=newdouble[n+1];}input(p,n,n+1);if(!dispel(p,n,n+1)){cout<<"奇异"<<endl;return0;}double*q=newdouble[n];for(inti=0;i<n;i++)q[i]=0;back(p,n,n+1,q);disp(q,n);delete[]q;for(inti=0;i<n;i++)delete[]p[i];delete[]p;}1.用列主元消去法解方程例2解方程组计算结果如下矩阵直接三角分解法算法：将方程组A*=b中的A分解为A=LU，其中L为单位下三角矩阵，U为上三角矩阵，则方程组A*=b化为解2个方程组Ly=b，U*=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求解方程组A*=bA=,b=结果为#include<iostream>usingnamespacestd;double**newMatri*(introw,intcol){double**p=newdouble*[row];//行for(inti=0;i<row;i++)//列p[i]=newdouble[col];returnp;}voiddelMatri*(double**p,introw){for(inti=0;i<row;i++)delete[]p[i];delete[]p;}voidinputMatri*(double**p,introw,intcol){for(inti=0;i<row;i++){cout<<"输入第"<<i+1<<"行:";for(intj=0;j<col;j++)cin>>p[i][j];}}voiddispMatri*(double**p,introw,intcol){for(inti=0;i<row;i++){for(intj=0;j<col;j++)cout<<p[i][j]<<'';cout<<endl;}}voiddispVector(double*q,intn){cout<<"====================================="<<endl;for(inti=0;i<n;i++)cout<<"*["<<i+1<<"]="<<q[i]<<endl;cout<<"====================================="<<endl;}voidinitMatri*(double**p,introw,intcol){for(inti=0;i<row;i++)for(intj=0;j<col;j++)p[i][j]=0;}doublesum(double**L,double**U,introw,intcol){intmin=(row>=col"col:row);doublesum=0;for(inti=0;i<min;i++)sum+=(U[i][col]*L[row][i]);returnsum;}voidresolve(double**A,double**L,double**U,introw,intcol){for(inti=0;i<col;i++)U[0][i]=A[0][i];//把A的第一行给U的第一行L[0][0]=A[0][0];for(inti=1;i<row;i++)//填充L的第一列L[i][0]=A[i][0]/A[0][0];for(inti=1;i<row;i++)for(intj=1;j<col;j++){if(i<=j)U[i][j]=A[i][j]-sum(L,U,i,j);elseL[i][j]=(A[i][j]-sum(L,U,i,j))/U[j][j];}for(inti=0;i<row;i++)L[i][i]=1;}doublesumRowY(double**L,double*y,introw){doublesum=0;for(inti=0;i<row;i++)sum+=(L[row][i]*y[i]);returnsum;}voidbackY(double**L,double*b,double*y,introw){for(inti=0;i<row;i++)y[i]=0;//初始化y向量全为零for(inti=0;i<row;i++)y[i]=b[i]-sumRowY(L,y,i);}doublesumRow*(double**U,double**,introw,intcol){doublesum=0;for(inti=row+1;i<col;i++)sum+=(U[row][i]**[i]);returnsum;}voidback*(double**U,double*y,double**,introw){for(inti=0;i<row;i++)*[i]=0;//初始化*向量全为零for(inti=row-1;i>=0;i--)*[i]=(y[i]-sumRow*(U,*,i,row))/U[i][i];}intmain(){cout<<"输入方阵大小n:";intn=0;cin>>n;cout<<"====================================="<<endl;cout<<"开场录入方阵数据....."<<endl;double**A=newMatri*(n,n);//开辟矩阵AinputMatri*(A,n,n);//录入数据到A中cout<<"录入方阵数据完毕......"<<endl;cout<<"====================================="<<endl<<endl;cout<<"开场录入列阵....."<<endl;cout<<"输入列阵:";double*b=newdouble[n];//开辟向量bfor(inti=0;i<n;i++)cin>>b[i];//录入数据到b中cout<<"录入列阵数据完毕....."<<endl;double**L=newMatri*(n,n);//开辟方阵LinitMatri*(L,n,n);//初始化L全为零double**U=newMatri*(n,n);//开辟方阵UinitMatri*(U,n,n);//初始化Uresolve(A,L,U,n,n);//分解A为L与Udouble*y=newdouble[n];//开辟向量ybackY(L,b,y,n);//回代求出ydouble**=newdouble[n];//开辟向量*back*(U,y,*,n);//回代求出*dispVector(*,n);//释放空间delMatri*(A,n);delMatri*(L,n);delMatri*(U,n);delete[]b;delete[]y;delete[]*;}迭代法雅可比迭代法算法：设方程组A*=b系数矩阵的对角线元素,M为迭代次数容许的最大值,ε为容许误差。①取初始向量*=，令k=0。②对i=1,2,…,n计算③如果，则输出，完毕；否则执行④。④如果k≥M,则不收敛，终止程序；否则,转②。程序与实例例4用雅可比迭代法解方程组结果为迭代次数为20#include<iostream>#include<cmath>usingnamespacestd;double**newMatri*(introw,intcol){double**p=newdouble*[row];//行for(inti=0;i<row;i++)//列p[i]=newdouble[col];returnp;}voiddelMatri*(double**p,introw){for(inti=0;i<row;i++)delete[]p[i];delete[]p;}voidinputMatri*(double**p,introw,intcol){for(inti=0;i<row;i++){cout<<"输入第"<<i+1<<"行:";for(intj=0;j<col;j++)cin>>p[i][j];}}voiddispMatri*(double**p,introw,intcol){for(inti=0;i<row;i++){for(intj=0;j<col;j++)cout<<p[i][j]<<'';cout<<endl;}}voiddispVector(double*q,intn){for(inti=0;i<n;i++)cout<<"*["<<i+1<<"]="<<q[i]<<'\t';cout<<endl;}doublesumRow(double**A,double**1,introw,intcol){doublesum=0;for(inti=0;i<col;i++){if(row==i)continue;sum+=(A[row][i]**1[i]);}returnsum;}voiditerat(double**A,double*b,double**1,double**2,introw){for(inti=0;i<row;i++)*2[i]=(b[i]-sumRow(A,*1,i,row))/A[i][i];}boolblow_error(double**1,double**2,introw,doublee){doublesum=0;for(inti=0;i<row;i++){sum+=abs(*2[i]-*1[i]);}if(sum<e)returntrue;elsereturnfalse;}intmain(){cout<<"输入方阵大小n:";intn=0;cin>>n;cout<<"====================================="<<endl;cout<<"开场录入方阵数据....."<<endl;double**A=newMatri*(n,n);//开辟矩阵AinputMatri*(A,n,n);//录入数据到A中cout<<"录入方阵数据完毕......"<<endl;cout<<"====================================="<<endl<<endl;cout<<"开场录入列阵....."<<endl;cout<<"输入列阵:";double*b=newdouble[n];//开辟向量bfor(inti=0;i<n;i++)cin>>b[i];//录入数据到b中cout<<"录入列阵数据完毕....."<<endl;cout<<"====================================="<<endl;cout<<"输入迭代次数容许的最大值:";intM=0;cin>>M;cout<<"输入容许误差:";doublee=0;cin>>e;cout<<"====================================="<<endl;double**1=newdouble[n];//开辟*1向量double**2=newdouble[n];//开辟*2向量cout<<"输入初始向量:";for(inti=0;i<n;i++)cin>>*1[i];cout<<"====================================="<<endl;intk=0;//迭代计数器for(;;){iterat(A,b,*1,*2,n);//迭代一次if(blow_error(*1,*2,n,e)){dispVector(*2,n);cout<<"迭代次数为:"<<k<<endl;return0;}else{if(k>=M){cout<<"不收敛"<<endl;return0;}else{k++;for(inti=0;i<n;i++)*1[i]=*2[i];}}}//释放空间delMatri*(A,n);delete[]b;delete[]*1;delete[]*2;}高斯-塞尔德迭代法算法:设方程组A*=b的系数矩阵的对角线元素,,M为迭代次数容许的最大值,ε为容许误差①取初始向量,令k=0。②对i=1,2,…,n计算③如果，则输出,完毕；否则执行④。④如果则不收敛，终止程序；否则,转②。程序与实例例5用高斯-塞尔德迭代法解方程组结果为#include<stdio.h>#include<conio.h>#include<malloc.h>#include<math.h>#defineN100main(){ inti; float**; 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); *=GauseSeidel(c,3); for(i=0;i<=2;i++)printf("*[%d]=%f\n",i,*[i]); getch();}float*GauseSeidel(float*a,intn){ inti,j,nu=0; float**,d*,d; *=(float*)malloc(n*sizeof(float)); for(i=0;i<=n-1;i++)*[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)**[j]; d*=(*(a+i*(n+1)+n)-d)/(*(a+i*(n+1)+i)); *[i]+=d*; } if(nu>=N) { printf("foldnumber\n"); } nu++; } while(fabs(d*)>1e-6); return*;}例6用雅可比迭代法解方程组迭代4次得解,假设用高斯-塞尔德迭代法则发散。#include<stdio.h>voidmain(void){ intk,n; double*[3]={7,2,5}; for(k=0;k<5;k++) { doublea,b; a=*[0]; b=*[1]; *[0]=(7-2.0**[1]+2**[2])/1; *[1]=(2-a-*[2])/1; *[2]=(5-2*a-2*b)/1; } for(n=0;n<3;n++) { printf("*[%d]=%8.6f\n",n,*[n]); }}用高斯-塞尔德迭代法解方程组迭代84次得解，假设用雅克比迭代法则发散。#include<stdio.h>#include<math.h>voidLOOP(floata[10][10],floatb[10],float*[10],int);voidmain(void){floata[10][10],b[10],*[10],A[10]; doubleS; intM,n,i,j; printf("请输入方阵阶数:"); scanf("%d",&n); printf("请输入最大允许迭代次数:"); scanf("%d",&M); printf("请按行输入各方程系数:"); for(i=0;i<=n-1;i++) { for(j=0;j<=n-1;j++) { scanf("%f",&a[i][j]); } } printf("请输入各方程值:"); for(i=0;i<=n-1;i++) { scanf("%f",&b[i]); } printf("请依次输入首次迭代*值:"); for(i=0;i<=n-1;i++) { scanf("%f",&*[i]); } do { S=0.0; for(i=0;i<=n-1;i++) { A[i]=*[i]; } LOOP(a,b,*,n); M--; for(i=0;i<=n-1;i++) { S=S+fabs(*[i]-A[i]); } }while(M>=0&&S>=0.000001); if(M>=0) { printf("迭代次数M=%d\n",M); for(i=0;i<=n-1;i++) { printf("*[%d]=%f\n",i,*[i]); } } else { printf("该迭代发散\n"); }}voidLOOP(floata[10][10],floatb[10],float*[10],intn){ floatS1,S2,A[10]; inti,j; for(i=0;i<=n-1;i++) { A[i]=*[i]; } for(i=0;i<=n-1;i++) { S1=0.0; S2=0.0; for(j=0;j<=i-1;j++) { S1=S1+a[i][j]**[j]; } for(j