东北大学数值分析实验报告

实验一迭代格式的比较问题提出设方程f(x)=x-3x–1=0有三个实根x=1.8793,x=-0.34727,x=-1.53209现采用下面三种不同计算格式，求f(x)=0的根x或xx=x=x=二、要求1、编制一个程序进行运算，最后打印出每种迭代格式的敛散情况；2、用事后误差估计来控制迭代次数，并且打印出迭代的次数；3、初始值的选取对迭代收敛有何影响；4、分析迭代收敛和发散的原因。三、目的和意义1、通过实验进一步了解方程求根的算法；2、认识选择计算格式的重要性；3、掌握迭代算法和精度控制；4、明确迭代收敛性与初值选取的关系。四、程序设计流程图五、源程序代码#include<iostream>#include<math.h>floatone(floatx0){x1=(3*x0+1)/(x0*x0)；return(x1)；}floattwo(floatx0){x2=(pow(x0,3)-1)/3；return(x2)；}floatthree(floatx0){x3=pow(3*x0+1,0.33333)；return(x3)；}voidmain(){floatx,x0;floatx1,x2,x3；intk1,k2,k3;printf(pleaseinputx0=")；scanf("%f"，&x）；x0=x;x1=one(x0);printf("第一个公式迭代结果:\n");while（fabs(x0-x1)>1e-5）{printf("x1=%6.5f\n",x1);x0=x1;x1=one(x0);k1++;}printf("x1=%6.5f\n",x1);printf("k1=%i",k1);x0=x;x2=two(x0);printf("第二个公式迭代结果为:\n");while（fabs(x0-x2)>1e-5）{printf("x2=%6.5f\n",x2);x0=x;x1=two(x0);k2++;}printf("x2=%6.5f\n",x2);printf("k2=%i",k2);x0=x;x3=three(x0);printf("第三个公式迭代结果为:\n");while（fabs(x0-x3)>1e-5）{printf("x3=%6.5f\n",x3);x0=x;x1=three(x0);k3++;}printf("x3=%6.5f\n",x3);printf("k3=%i",k3);}运行结果如下：结果分析:初值对迭代结果影响解析：实验二线性方程组的直接算法一、问题提出给出下列几个不同类型的线性方程组，请用适当算法计算其解。设线性方程组=x=(1,-1,0,1,2,0,3,1,-1,2)Gsuss列主元消去法#include<math.h>#include<conio.h>#include<stdio.h>#defineMAX100typedefstruct{introw,col;floatMAT[MAX][MAX];floatSolution[MAX];}Matrix;voidGauss(Matrix*M);voidMBack(Matrix*M);voidMSave(Matrix*M);voidMInput(Matrix*M);voidMOutput(Matrix*M);voidSolution(Matrix*M);voidMSort(Matrix*M,intn);voidmain(){ printf("列主元方法如下：\n");MatrixMat;MInput(&Mat);MSave(&Mat);Gauss(&Mat);MSave(&Mat);if(Mat.row==Mat.col-1){MBack(&Mat);Solution(&Mat);}printf("Pressanykeytohalt...");getch();}voidMInput(Matrix*M){inti,j;printf("输入行数:");scanf("%d",&M->row);printf("输入列数:");scanf("%d",&M->col);for(i=0;i<M->row;i++){printf("第%d行:",i+1);for(j=0;j<M->col;j++){scanf("%f",&M->MAT[i][j]);}}for(i=0;i<M->row;i++)M->Solution[i]=0;}voidMOutput(Matrix*M){inti,j;printf("MATRIX:\n");for(i=0;i<M->row;i++){for(j=0;j<M->col;j++)printf("%10.3f",M->MAT[i][j]);printf("\n");}printf("---END----\n");}voidGauss(Matrix*M){inti,j,k;floattemp;for(i=0;i<M->row-1;i++){MSort(M,i);MOutput(M);for(j=i+1;j<M->row;j++){temp=M->MAT[j][i];for(k=0;k<M->col;k++)if(temp!=0){M->MAT[j][k]/=temp;M->MAT[j][k]*=M->MAT[i][i];M->MAT[j][k]-=M->MAT[i][k];}}}MOutput(M);}voidMSort(Matrix*M,intn){inti,j,k;floattemp[MAX];for(i=n;i<M->row-1;i++){for(j=n;j<M->row-i-1;j++){if(fabs(M->MAT[j][n])<fabs(M->MAT[j+1][n])){for(k=0;k<M->col;k++){temp[k]=M->MAT[j+1][k];M->MAT[j+1][k]=M->MAT[j][k];M->MAT[j][k]=temp[k];}}}}}voidMBack(Matrix*M){inti,j;floatsum;M->Solution[M->row-1]=M->MAT[M->row-1][M->col-1]M->MAT[M->row-1][M->row-1];for(i=M->row-2;i>=0;i--){sum=M->MAT[i][M->col-1];for(j=i+1;j<M->row;j++)sum-=M->MAT[i][j]*M->Solution[j];M->Solution[i]=sum/M->MAT[i][i];}}voidSolution(Matrix*M){inti;printf("Solution:\n");for(i=0;i<M->row;i++)printf("X[%d]=%f\n",i+1,M->Solution[i]);printf("\n---END---\n");}voidMSave(Matrix*M){inti,j;FILE*eryar;eryar=fopen("Matrix.txt","a");fprintf(eryar,"--------BEGIN--------\n");for(i=0;i<M->row;i++){for(j=0;j<M->col;j++)fprintf(eryar,"%10.3f",M->MAT[i][j]);fprintf(eryar,"\n");}fclose(eryar);}2、设对称正定阵系数阵线方程组=x=(1,-1,0,2,1,-1,0,2)平方根法:源程序:#include<iostream>#include<cmath>#include<cstdlib>usingnamespacestd;intmain() { intn,i,j,k,m; cout<<"输入维数:"; cin>>n; double**A=newdouble*[(n+1)]; for(i=1;i<=n;i++) A[i]=newdouble[n+1]; double*b=newdouble[n+1]; double*x=newdouble[n+1]; double*y=newdouble[n+1]; cout<<"输入系数对称正定矩阵A[][]:"<<endl; for(i=1;i<=n;i++) for(j=1;j<=n;j++) cin>>A[i][j]; cout<<"输入向量b[]:"; for(i=1;i<=n;i++) cin>>b[i]; cout<<endl; for(k=1;k<=n;k++) { doublesum=0; for(m=1;m<=k-1;m++) { sum=sum+pow(A[k][m],2.0); } sum=A[k][k]-sum; A[k][k]=sqrt(sum); for(i=k+1;i<=n;i++) { doubletemp1=0; for(m=1;m<=k-1;m++) { temp1=temp1+A[i][m]*A[k][m]; } temp1=A[i][k]-temp1; A[i][k]=temp1/A[k][k]; } doubletemp2=0; for(m=1;m<=k-1;m++) { temp2=temp2+A[k][m]*y[m]; } y[k]=(b[k]-temp2)/A[k][k]; } x[8]=y[8]/A[8][8]; for(k=n-1;k>=1;k--) { doubletemp3=0; for(m=k+1;m<=n;m++) { temp3=temp3+A[m][k]*x[m]; } x[k]=(y[k]-temp3)/A[k][k]; } cout<<"输出结果向量x[]:"<<endl; for(i=1;i<=n;i++)cout<<x[i]<<endl; return0;}3、三对角形线性方程组=x=(2,1,-3,0,1,-2,3,0,1,-1)追赶法:

源程序#include<iostream>#include<cmath>#include<cstdlib>usingnamespacestd;intmain(){ intn,i; cout<<"输入系数矩阵的维数:"; cin>>n; double*a=newdouble[n+1]; double*c=newdouble[n+1]; double*d=newdouble[n+1]; double*b=newdouble[n+1]; double*x=newdouble[n+1]; double*y=newdouble[n+1]; cout<<"输入系数矩阵A[]数据："<<endl; for(i=1;i<=n;i++)cin>>a[i]; for(i=1;i<=n;i++)cin>>c[i]; for(i=1;i<=n;i++)cin>>d[i]; cout<<"输入b[]:"<<endl; for(i=1;i<=n;i++)cin>>b[i]; for(i=1;i<=n-1;i++) { c[i]=c[i]/a[i]; a[i+1]=a[i+1]-d[i+1]*c[i]; } cout<<"输出解向量a[]:"<<endl; for(i=1;i<=n;i++)cout<<a[i]<<endl; cout<<"输出解向量c[]:"<<endl; for(i=1;i<=n;i++)cout<<c[i]<<endl; y[1]=b[1]/a[1]; for(i=2;i<=n;i++) { y[i]=(b[i]-d[i]*y[i-1])/a[i]; } cout<<"输出解向量y[]:"<<endl; for(i=1;i<=n;i++)cout<<y[i]<<endl; x[n]=y[n]; for(i=n-1;i>=1;i--) { x[i]=y[i]-c[i]*x[i+1]; } cout<<"输出解向量x[]:"<<endl; for(i=1;i<=n;i++)cout<<x[i]<<endl; system("pause"); return0;}运行结果课题三线性方程组的迭代法一、问题提出对课题二所列目的和意义的线性方程组，试分别选用Jacobi迭代法，Gauss-Seidol迭代法和SOR方法计算其解。二、要求1、体会迭代法求解线性方程组，并能与消去法做以比较；2、分别对不同精度要求，如由迭代次数体会该迭代法的收敛快慢；3、对方程组2，3使用SOR方法时，选取松弛因子=0.8，0.9，1，1.1，1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者；4、给出各种算法的设计程序和计算结果。三、目的和意义1、通过上机计算体会迭代法求解线性方程组的特点，并能和消去法比较；2、运用所学的迭代法算法，解决各类线性方程组，编出算法程序；3、体会上机计算时，终止步骤<或k>（予给的迭代次数），对迭代法敛散性的意义；4、体会初始解x，松弛因子的选取，对计算结果的影响。1．Jacobi迭代法源程序:#include<iostream>#include<cmath>#include<cstdlib>usingnamespacestd;intmain(){ intn,N; inti,j,k; doublee; cout<<"输入维数n:"; cin>>n; cout<<"输入最大迭代次数N:"; cin>>N; cout<<"输入精度e:"; cin>>e; double**A=newdouble*[(n+1)]; for(i=1;i<=n;i++) A[i]=newdouble[n+1]; double*b=newdouble[n+1]; double*x0=newdouble[n+1]; double*x=newdouble[n+1]; cout<<"输入系数矩阵A[][]:"<<endl; for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { cin>>A[i][j]; } } cout<<"输入右端项b[]:"<<endl; for(i=1;i<=n;i++)cin>>b[i]; cout<<"初始化解向量x0[]:"<<endl; for(i=1;i<=n;i++)cin>>x0[i]; for(k=1;k<=N;k++) { for(i=1;i<=n;i++) { doublesum=0; for(j=1;j<=n;j++) { sum=sum+A[i][j]*x0[j]; } x[i]=x0[i]+(b[i]-sum)/A[i][i]; } doublemax=0; for(i=1;i<=n;i++) { doublem; m=fabs(x[i]-x0[i]); if(m>max)max=m; } if(max<=e) { cout<<"输出迭代次k:"<<k<<endl; cout<<"迭代后的解x[]:"<<endl; for(i=1;i<=n;i++)cout<<x[i]<<endl; system("pause"); return0; } else { for(i=1;i<=n;i++) { x0[i]=x[i]; } } } cout<<"已达到最大迭代次数:"<<N<<endl; cout<<"迭代后的解x[]:"<<endl; for(i=1;i<=n;i++)cout<<x[i]<<endl; system("pause");return0;}2．Gauss-Seidol迭代法源程序#include<iostream>#include<cmath>#include<cstdlib>usingnamespacestd;intmain(){ intn,N; inti,j,k; doublee; cout<<"输入维数n:"; cin>>n; cout<<"输入最大迭代次数N:"; cin>>N; cout<<"输入精度e:"; cin>>e; double**A=newdouble*[(n+1)]; for(i=1;i<=n;i++) A[i]=newdouble[n+1]; double*b=newdouble[n+1]; double*x0=newdouble[n+1]; double*x=newdouble[n+1]; cout<<"输入系数矩阵A[][]:"<<endl; for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { cin>>A[i][j]; } } cout<<"输入右端项b[]:"<<endl; for(i=1;i<=n;i++)cin>>b[i]; cout<<"初始化解向量x0[]:"<<endl; for(i=1;i<=n;i++)cin>>x0[i]; for(k=1;k<=N;k++) { for(i=1;i<=n;i++) { doublesum=0; for(j=1;j<=n;j++) { if(j<i)sum=sum+A[i][j]*x[j]; else sum=sum+A[i][j]*x0[j]; } x[i]=x0[i]+(b[i]-sum)/A[i][i]; } doublemax=0; for(i=1;i<=n;i++) { doublem; m=fabs(x[i]-x0[i]); if(m>max)max=m; } if(max<=e) { cout<<"输出迭代次k:"<<k<<endl; cout<<"迭代后的解x[]:"<<endl; for(i=1;i<=n;i++)cout<<x[i]<<endl;system("pause"); return0; } else { for(i=1;i<=n;i++) { x0[i]=x[i]; } } } cout<<"已达到最大迭代次数:"<<N<<endl; cout<<"迭代后的解x[]:"<<endl; for(i=1;i<=n;i++)cout<<x[i]<<endl;system("pause");return0;}3．SOR方法源程序#include<iostream>#include<vector>#include<cstdlib>#include<cmath>usingnamespacestd;intmain(){intn;cout<<"输入维数:";cin>>n; vector<vector<double>>va; vector<double>vd; cout<<"输入系数矩阵:"<<endl; for(inti=0;i<n;++i){ for(intj=0;j<=n;++j){ doubledtemp; cin>>dtemp; vd.push_back(dtemp); } va.push_back(vd); vd.clear(); } cout<<"输入初始解向量:"; vector<double>vtemp; for(inti=0;i<n;++i){ doubledtemp; cin>>dtemp; vd.push_back(dtemp); } cout<<"输入精度:"; doubled; cin>>d; cout<<"输入最大迭代次数:"; intmcount; cin>>mcount; cout<<"输入松弛因子:"; doublexd; cin>>xd; vector<vector<double>>::iteratoria,iaend; ia=va.begin(); boolflag=true; intcount=0; while(flag){ doubledmax=0.0; vtemp.assign(vd.begin(),vd.end()); for(inti=0;i<n;++i){ doubledtemp=0.0; for(intj=0;j<n;++j){ if(j<i){ dtemp+=(*(ia+i))[j]*vd[j]; } elseif(j>=i){ dtemp+=(*(ia+i))[j]*vtemp[j]; } } dtemp=(*(ia+i))[n]-dtemp; vd[i]=xd*dtemp/(*(ia+i))[i]; vd[i]+=vtemp[i]; } ++count; for(inti=0;i<n;++i){ if(dmax<fabs(vd[i]-vtemp[i])) dmax=fabs(vd[i]-vtemp[i]); } if(dmax<d){ cout<<"输出迭代次数:"<<count<<endl; cout<<d-dmax<<endl; cout<<"输出迭代解:"<<"("; for(inti=0;i<n;++i){ if(i==n-1) cout<<vd[i]<<")"; else cout<<vd[i]<<","; } flag=false; } elseif(count>mcount){ cout<<"已达到最大迭代次数!"<<endl; cout<<"输出迭代解:"<<"("; for(inti=0;i<n;++i){ if(i==n-1) cout<<vd[i]<<")"; else cout<<vd[i]<<","; } flag=false; }} system("pause");return0;}方程组一1．Jacobi迭代法2．Gauss-Seidol迭代法方程组二 1．Jacobi迭代法精度e=0.0012．Gauss-Seidol迭代法3．SOR方法松弛因子w=0.8方程组三 1．Jacobi迭代法2．Gauss-Seidol迭代法3．SOR方法精度e=0.001松弛因子w=0.8松弛因子w=1.0松弛因子w=1.2结果分析 问题1.体会迭代法求解线性方程组，并能与消去法做以比较；消去法如果不考虑计算过程中的舍入误差，能求出方程组的精确解。迭代法即使不考虑计算过程中的舍入误差，迭代法也很难获得精确解。 问题2.对于不同的精度，迭代次数也不同。对第三个方程组用GS迭代法采用不同精度，迭代次数不同。精度越精确，需要迭代次数越多。 问题3．从题中W=0.81.01.2时，可以看出w=1.0时迭代的最快。课题四数值积分一、问题提出选用复合梯形公式，复合Simpson公式，Romberg算法，计算I=I=I=I=二、要求编制数值积分算法的程序；分别用两种算法计算同一个积分，并比较其结果；分别取不同步长，试比较计算结果（如n=10,20等）；给定精度要求，试用变步长算法，确定最佳步长。三、目的和意义深刻认识数值积分法的意义；明确数值积分精度与步长的关系；根据定积分的计算方法，可以考虑二重积分的计算问题。1复化梯形：复化simpson：Romberg求积法：源程序：#include<stdio.h>#include<conio.h>#include<math.h>floatf(floatx){returnsqrt(4-sin(x)*sin(x));}floatRomberg(floata,floatb,float(*f)(float),floatepsilon){