東北大學(xué)數(shù)值分析實(shí)驗(yàn)報(bào)告

實(shí)驗(yàn)一迭代格式的比較問題提出設(shè)方程f(x)=x-3x–1=0有三個(gè)實(shí)根x=1.8793,x=-0.34727,x=-1.53209現(xiàn)采用下面三種不同計(jì)算格式，求f(x)=0的根x或xx=x=x=二、要求1、編制一個(gè)程序進(jìn)行運(yùn)算，最后打印出每種迭代格式的斂散情況；2、用事后誤差估計(jì)來控制迭代次數(shù)，并且打印出迭代的次數(shù)；3、初始值的選取對(duì)迭代收斂有何影響；4、分析迭代收斂和發(fā)散的原因。三、目的和意義1、通過實(shí)驗(yàn)進(jìn)一步了解方程求根的算法；2、認(rèn)識(shí)選擇計(jì)算格式的重要性；3、掌握迭代算法和精度控制；4、明確迭代收斂性與初值選取的關(guān)系。四、程序設(shè)計(jì)流程圖五、源程序代碼#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("第一個(gè)公式迭代結(jié)果:\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("第二個(gè)公式迭代結(jié)果為:\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("第三個(gè)公式迭代結(jié)果為:\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);}運(yùn)行結(jié)果如下：結(jié)果分析:初值對(duì)迭代結(jié)果影響解析：實(shí)驗(yàn)二線性方程組的直接算法一、問題提出給出下列幾個(gè)不同類型的線性方程組，請(qǐng)用適當(dāng)算法計(jì)算其解。設(shè)線性方程組=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("輸入行數(shù):");scanf("%d",&M->row);printf("輸入列數(shù):");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、設(shè)對(duì)稱正定陣系數(shù)陣線方程組=x=(1,-1,0,2,1,-1,0,2)平方根法:源程序:#include<iostream>#include<cmath>#include<cstdlib>usingnamespacestd;intmain() { intn,i,j,k,m; cout<<"輸入維數(shù):"; 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<<"輸入系數(shù)對(duì)稱正定矩陣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<<"輸出結(jié)果向量x[]:"<<endl; for(i=1;i<=n;i++)cout<<x[i]<<endl; return0;}3、三對(duì)角形線性方程組=x=(2,1,-3,0,1,-2,3,0,1,-1)追趕法:

源程序#include<iostream>#include<cmath>#include<cstdlib>usingnamespacestd;intmain(){ intn,i; cout<<"輸入系數(shù)矩陣的維數(shù):"; 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<<"輸入系數(shù)矩陣A[]數(shù)據(jù)："<<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;}運(yùn)行結(jié)果課題三線性方程組的迭代法一、問題提出對(duì)課題二所列目的和意義的線性方程組，試分別選用Jacobi迭代法，Gauss-Seidol迭代法和SOR方法計(jì)算其解。二、要求1、體會(huì)迭代法求解線性方程組，并能與消去法做以比較；2、分別對(duì)不同精度要求，如由迭代次數(shù)體會(huì)該迭代法的收斂快慢；3、對(duì)方程組2，3使用SOR方法時(shí)，選取松弛因子=0.8，0.9，1，1.1，1.2等，試看對(duì)算法收斂性的影響，并能找出你所選用的松弛因子的最佳者；4、給出各種算法的設(shè)計(jì)程序和計(jì)算結(jié)果。三、目的和意義1、通過上機(jī)計(jì)算體會(huì)迭代法求解線性方程組的特點(diǎn)，并能和消去法比較；2、運(yùn)用所學(xué)的迭代法算法，解決各類線性方程組，編出算法程序；3、體會(huì)上機(jī)計(jì)算時(shí)，終止步驟<或k>（予給的迭代次數(shù)），對(duì)迭代法斂散性的意義；4、體會(huì)初始解x，松弛因子的選取，對(duì)計(jì)算結(jié)果的影響。1．Jacobi迭代法源程序:#include<iostream>#include<cmath>#include<cstdlib>usingnamespacestd;intmain(){ intn,N; inti,j,k; doublee; cout<<"輸入維數(shù)n:"; cin>>n; cout<<"輸入最大迭代次數(shù)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<<"輸入系數(shù)矩陣A[][]:"<<endl; for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { cin>>A[i][j]; } } cout<<"輸入右端項(xiàng)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<<"已達(dá)到最大迭代次數(shù):"<<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<<"輸入維數(shù)n:"; cin>>n; cout<<"輸入最大迭代次數(shù)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<<"輸入系數(shù)矩陣A[][]:"<<endl; for(i=1;i<=n;i++) { for(j=1;j<=n;j++) { cin>>A[i][j]; } } cout<<"輸入右端項(xiàng)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<<"已達(dá)到最大迭代次數(shù):"<<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<<"輸入維數(shù):";cin>>n; vector<vector<double>>va; vector<double>vd; cout<<"輸入系數(shù)矩陣:"<<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<<"輸入最大迭代次數(shù):"; 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<<"輸出迭代次數(shù):"<<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<<"已達(dá)到最大迭代次數(shù)!"<<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結(jié)果分析 問題1.體會(huì)迭代法求解線性方程組，并能與消去法做以比較；消去法如果不考慮計(jì)算過程中的舍入誤差，能求出方程組的精確解。迭代法即使不考慮計(jì)算過程中的舍入誤差，迭代法也很難獲得精確解。 問題2.對(duì)于不同的精度，迭代次數(shù)也不同。對(duì)第三個(gè)方程組用GS迭代法采用不同精度，迭代次數(shù)不同。精度越精確，需要迭代次數(shù)越多。 問題3．從題中W=0.81.01.2時(shí)，可以看出w=1.0時(shí)迭代的最快。課題四數(shù)值積分一、問題提出選用復(fù)合梯形公式，復(fù)合Simpson公式，Romberg算法，計(jì)算I=I=I=I=二、要求編制數(shù)值積分算法的程序；分別用兩種算法計(jì)算同一個(gè)積分，并比較其結(jié)果；分別取不同步長(zhǎng)，試比較計(jì)算結(jié)果（如n=10,20等）；給定精度要求，試用變步長(zhǎng)算法，確定最佳步長(zhǎng)。三、目的和意義深刻認(rèn)識(shí)數(shù)值積分法的意義；明確數(shù)值積分精度與步長(zhǎng)的關(guān)系；根據(jù)定積分的計(jì)算方法，可以考慮二重積分的計(jì)算問題。1復(fù)化梯形：復(fù)化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){