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

数值分析设计实验实验报告课题一迭代格式的比较问题提出设方程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<stdio.h>#include<math.h>voidmain(){ floatx1,x2,x3,q,a,z,p,e=0.00001; x1=-1.0000;x2=-1.0000;x3=1.0000; inti,y=3; printf("0%f%f%f\n",x1,x2,x3); q=x1-p;a=x2-p;z=x3-p; for(i=1;i<=60;i++) { if(q<e&&q>(0-e)) gotoa; else { p=x1; x1=(3*x1+1)/(x1*x1); printf("%d1%f\t",i,x1); q=x1-p; } a: if(a<e&&a>(0-e)) gotoz; else { p=x2; x2=(x2*x2*x2-1)/3; printf("%d2%f\t",i,x2); a=x2-p; } z: if(z<e&&z>(0-e)) gotoend; else { p=x3; x3=pow((3*x3+1),1.0/y); printf("%d3%f\n",i,x3); z=x3-p; }end:; }}六。程序运行结果七.程序运行结果讨论和分析：对于迭代格式一、二、三对于初值为-1.0000，-1.0000，1.0000分别迭代了37次，8次，10次，由此可知，简单迭代法的收敛性取决于迭代函数，以及初值x的选取，并且对初值的选取要求较高，需谨慎选取。课题二线性方程组的直接算法一、问题提出给出下列几个不同类型的线性方程组，请用适当算法计算其解。设线性方程组=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(){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");}}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]);}voidMSave(Matrix*M){inti,j;FILE*eryar;eryar=fopen("Matrix.txt","a");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;; system("pause"); 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;}程序运行结果分析误差分析：列主元高斯法解决了一部分大数吃掉小数的现象，但在计算机运算时，仍会对中间数据进行舍入，影响结果的精度，这种误差同样也存在其他两种算法中；对于三对角阵，选用追赶法可以得到更准确的解。课题三线性方程组的迭代法一、问题提出1、设线性方程组=x=(1,-1,0,1,2,0,3,1,-1,2)J迭代算法：程序设计流程图：源程序代码：#include<stdlib.h>#include<stdio.h>#include<math.h>voidmain(){floata[50][51],x1[50],x2[50],temp=0,fnum=0;inti,j,m=200,n=10,e,bk=0;for(i=1;i<n+1;i++)x1[i]=0;printf("请输入系数矩阵:\n");for(i=1;i<n+1;i++){j=1;while(j<n+1){scanf("%f",&a[i][j]);j++;}}printf("输入方程组的常数项:\n");for(i=1;i<n+1;i++){scanf("%f",&a[i][n+1]);} printf("请输入精度要求:\n"); scanf("%d",&e);while(m!=0){for(i=1;i<n+1;i++){for(j=1;j<n+1;j++){if(j!=i)temp=a[i][j]*x1[j]+temp;}x2[i]=(a[i][n+1]-temp)/a[i][i];temp=0;}for(i=1;i<n+1;i++){fnum=float(fabs(x1[i]-x2[i]));if(fnum>temp)temp=fnum;}if(temp<=pow(10,-4))bk=1;for(i=1;i<n+1;i++)x1[i]=x2[i]; m--;}printf("原方程组的解为:\t");for(i=1;i<n+1;i++){if((x1[i]-x2[i])<=e||(x2[i]-x1[i])<=e) { printf("x%d=%7.4f",i,x1[i]); }} printf("\n");}运行结果：2、设对称正定阵系数阵线方程组=x=(1,-1,0,2,1,-1,0,2)GS迭代算法：#include<iostream.h>#include<math.h>#include<stdio.h>constintm=11;voidmain(){intchoice=1;while(choice==1){doublea[m][m],b[m],e,x[m],y[m],w,se,max;intn,i,j,N,k;cout<<"Gauss-Seidol迭代法"<<endl; cout<<"请输入方程的个数：";cin>>n;for(i=1;i<=n;i++){cout<<"请输入第"<<i<<"个方程的各项系数：";for(j=1;j<=n;j++)cin>>a[i][j];}cout<<"请输入各个方程等号右边的常数项:\n";for(i=1;i<=n;i++){cin>>b[i];}cout<<"请输入最大迭代次数：";cin>>N;cout<<"请输入最大偏差：";cin>>e;for(i=1;i<=n;i++){x[i]=0;y[i]=x[i];}k=0;while(k!=N){k++;for(i=1;i<=n;i++){w=0;for(j=1;j<=n;j++){if(j!=i)w=w+a[i][j]*y[j];}y[i]=(b[i]-w)/double(a[i][i]);}max=fabs(x[1]-y[1]);for(i=1;i<=n;i++){se=fabs(x[i]-y[i]);if(se>max)max=se;}if(max<e){cout<<endl;for(i=1;i<=n;i++)cout<<"x"<<i<<"="<<y[i]<<endl;break;}for(i=1;i<=n;i++){x[i]=y[i];}}if(k==N)cout<<"迭代失败！！"<<endl; choice=0;}}3、三对角形线性方程组=x=(2,1,-3,0,1,-2,3,0,1,-1)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;}精度e=0.001 松弛因子w=0.8四、程序运行结果讨论和分析：迭代法即使不考虑计算过程中的舍入误差，迭代法也很难获得精确解。采用SOR迭代法，不同的松弛因子，会影响迭代次数。经过计算可知Gauss-Seidel方法比Jacobi方法计算量小,可是精确度底。SOR是进一步改进Gauss-Seidel，并且比Jacobi,Gauss-Seidel方法收敛速度快,改变松弛因子的取值范围来可以得到Jacobi,Gauss-Seidel方法.课题四数值积分一、问题提出选用复合梯形公式，复合Simpson公式，Romberg算法，计算I=I=I=I=二、要求编制数值积分算法的程序；分别用两种算法计算同一个积分，并比较其结果；分别取不同步长，试比较计算结果（如n=10,20等）；给定精度要求，试用变步长算法，确定最佳步长。三、目的和意义深刻认识数值积分法的意义；明确数值积分精度与步长的关系；根据定积分的计算方法，可以考虑二重积分的计算问题。四.程序及流程图：计算函数1，复化simpson：#include<stdio.h>#include<math.h>voidmain(){ doublea=0,b=1; intn=10;//或20 doubleh