当前位置:文档之家› 东北大学数值分析实验报告

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

数值分析设计实验实验报告课题一 迭代格式的比较一、问题提出设方程f(x)=x3- 3x –1=0 有三个实根 x*1=1.8793 ,x *2=-0.34727 ,x *3=-1.53209现采用下面三种不同计算格式,求 f(x)=0的根 x *1 或x *21、 x =213x x + 2、 x = 313-x3、 x = 313+x二、要求1、编制一个程序进行运算,最后打印出每种迭代格式的敛散情况;2、用事后误差估计k k x x -+1〈ε来控制迭代次数,并且打印出迭代的次数;3、初始值的选取对迭代收敛有何影响;4、分析迭代收敛和发散的原因。

三、目的和意义1、通过实验进一步了解方程求根的算法;2、认识选择计算格式的重要性;3、掌握迭代算法和精度控制;4、明确迭代收敛性与初值选取的关系。

四、程序设计流程图五、源程序代码#include<stdio.h>#include<math.h>void main(){float x1,x2,x3,q,a,z,p,e=0.00001;x1=-1.0000;x2=-1.0000;x3=1.0000;int i,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))goto a;else{ p=x1;x1=(3*x1+1)/(x1*x1);printf("%d 1 %f\t",i,x1);q=x1-p;}a: if(a<e&&a>(0-e))goto z;else{ p=x2;x2=(x2*x2*x2-1)/3;printf("%d 2 %f\t",i,x2);a=x2-p;}z: if(z<e&&z>(0-e))goto end;else{p=x3;x3=pow((3*x3+1),1.0/y);printf("%d 3 %f\n",i,x3);z=x3-p;}end:;}}六。

程序运行结果七.程序运行结果讨论和分析:对于迭代格式一、二、三对于初值为-1.0000,-1.0000,1.0000分别迭代了37次,8次,10次,由此可知,简单迭代法的收敛性取决于迭代函数,以及初值x 的选取,并且对初值的选取要求较高,需谨慎选取。

课题二 线性方程组的直接算法一、问题提出给出下列几个不同类型的线性方程组,请用适当算法计算其解。

1、 设线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------------------1368243810041202913726422123417911101610352431205362177586832337616244911315120130123122400105635680000121324⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x =⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-2119381346323125x *= ( 1, -1, 0, 1, 2, 0, 3, 1, -1, 2 )TGsuss 列主元消去法#include <math.h>#include <conio.h>#include <stdio.h>#define MAX 100typedef struct{int row,col;float MAT[MAX][MAX];float Solution[MAX];}Matrix;void Gauss(Matrix *M);void MBack(Matrix *M);void MSave(Matrix *M);void MInput(Matrix *M);void MOutput(Matrix *M); void Solution(Matrix *M); void MSort(Matrix *M,int n);void main(){Matrix Mat;MInput(&Mat);MSave(&Mat);Gauss(&Mat);MSave(&Mat);if(Mat.row==Mat.col-1) {MBack(&Mat);Solution(&Mat);}printf("Press any key to halt...");getch();}void MInput(Matrix *M){int i,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;}void MOutput(Matrix *M){int i,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");}}void Gauss(Matrix *M){int i,j,k;float temp;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);}void MSort(Matrix *M,int n){int i,j,k;float temp[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];}}}}}void MBack(Matrix *M){int i,j;float sum;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];}}void Solution(Matrix *M){int i;printf("Solution:\n");for(i=0;i<M->row;i++)printf("X[%d] = %f\n",i+1,M->Solution[i]); }void MSave(Matrix *M){int i,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、设对称正定阵系数阵线方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------------19243360021411035204111443343104221812334161206538114140231212200420424⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡87654321x x x x x x x x = ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---4515229232060x *= ( 1, -1, 0, 2, 1, -1, 0, 2 )T平方根法#include<iostream>#include<cmath>#include<cstdlib>using namespace std;int main(){int n,i,j,k,m;cout<<"输入维数:";cin>>n;double **A=new double*[(n+1)];for(i=1;i<=n;i++)A[i]=new double[n+1];double *b=new double[n+1];double *x=new double[n+1];double *y=new double[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++){double sum=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++){double temp1=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];}double temp2=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--){double temp3=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");return 0;}3、三对角形线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡------------------4100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x = ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----5541412621357x *= ( 2, 1, -3, 0, 1, -2, 3, 0, 1, -1 )T追赶法#include<iostream>#include<cmath>#include<cstdlib>using namespace std;int main(){int n,i;cout<<"输入系数矩阵的维数:";cin>>n;double *a=new double[n+1];double *c=new double[n+1];double *d=new double[n+1];double *b=new double[n+1];double *x=new double[n+1];double *y=new double[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");return 0;}四、 程序运行结果分析误差分析:列主元高斯法解决了一部分大数吃掉小数的现象,但在计算机运算时,仍会对中间数据进行舍入,影响结果的精度,这种误差同样也存在其他两种算法中;对于三对角阵,选用追赶法可以得到更准确的解。

相关主题