[高等教育]东北大学数值分析实验报告.doc

数值分析实验报告课题一 迭代格式的比较一、 问题提出设方程f(x)=x- 3x 1=0 有三个实根 x=1.8793 , x=-0.34727 ,x=-1.53209现采用下面三种不同计算格式，求 f(x)=0的根 x 或x 1、 x = 2、 x = 3、 x = 二、要求1、编制一个程序进行运算，最后打印出每种迭代格式的敛散情况；2、用事后误差估计来控制迭代次数，并且打印出迭代的次数；3、初始值的选取对迭代收敛有何影响；4、分析迭代收敛和发散的原因。三、目的和意义1、通过实验进一步了解方程求根的算法；2、认识选择计算格式的重要性；3、掌握迭代算法和精度控制；4、明确迭代收敛性与初值选取的关系。程序代码：#include<iostream>#include<math.h>int k=1,k1=1,k2=1,k3=1;float x1,x2,x3;float one(float x0)x1=(3*x0+1)/(x0*x0);return(x1);float two(float x0)x2=(pow(x0,3)-1)/3;return(x2);float three(float x0)x3=pow(3*x0+1,0.33333);return(x3);main()float x,x0;printf("输入x0=");scanf("%f",&x);x0=x;x1=one(x0);printf("第一个公式迭代结果: n");while(fabs(x0-x1)>1e-5)printf("x1=%6.5fn",x1);x0=x1;x1=one(x0);k1+;printf("x1=%6.5f n",x1);printf("k1=%in",k1);x0=x;x2=two(x0);printf("第二个公式迭代结果: n");while(fabs(x0-x2)>1e-5)printf("x2=%6.5fn",x2);x0=x2;x2=two(x0);k2+;printf("k2=%in",k2);x0=x;x3=three(x0);printf("第三个公式迭代结果: n");while(fabs(x0-x3)>1e-5)printf("x3=%6.5fn",x3);x0=x3;x3=three(x0);k3+;printf("x3=%6.5fn",x3);printf("k3=%in",k3);scanf("%");实验结果：四、程序运行结果讨论和分析：对于第一种迭代格式，收敛区间-8.2 -0.4，在该收敛区间内迭代收敛于-1.53209，只能求得方程的一个根；对于第二种迭代格式，收敛区间-1.5 1.8，在该收敛区间内迭代收敛于-0.34730，同样只能求得方程的一个根；对于第三种迭代格式，收敛区间-0.3 +)，在该收敛区间内迭代收敛于1.87937，只能求得方程的一个根；由以上结果很容易发现，初值的选取对迭代敛散性有很大影响。以第一种迭代格式为例，当初值大于等于-0.3时，迭代格式发散；当初值小于等于-8.3时，迭代格式也发散；只有初值在-0.3和-8.3之间时，迭代格式才收敛于1.53209。其他迭代格式也有这样的性质，即收敛于某个数值区间，超出这个区间迭代格式就是发散的，这就是所谓迭代格式的收敛性。对于不同迭代格式在不同区间具有不同的敛散性的原因，我认为可以从一下两方面理解：1、迭代法是一种逐次逼近法，其基本思想是将隐式方程归结为一组显式的计算公式，就是说，迭代过程实质上是个逐步显式化的过程。2、我们可以用几何图像来更好地理解迭代过程。由图可知，在某些区间选取的初始值随着迭代次数的增加会越来越逼近精确值，即收敛于精确值，而在另外一些区间选取的初始值随着迭代次数的增加却离精确值越来越远，即不会收敛于一个确定值。课题二 线性方程组的直接算法一、问题提出给出下列几个不同类型的线性方程组，请用适当算法计算其解。 1、设线性方程组=x= ( 1, -1, 0, 1, 2, 0, 3, 1, -1, 2 ) 2、设对称正定阵系数阵线方程组= x = ( 1, -1, 0, 2, 1, -1, 0, 2 ) 3、三对角形线性方程组 = x= ( 2, 1, -3, 0, 1, -2, 3, 0, 1, -1 ) 二、要求1、 对上述三个方程组分别利用Gauss顺序消去法与Gauss列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）；2、 应用结构程序设计编出通用程序；3、 比较计算结果，分析数值解误差的原因；4、 尽可能利用相应模块输出系数矩阵的三角分解式。三、目的和意义1、通过该课题的实验，体会模块化结构程序设计方法的优点；2、运用所学的计算方法，解决各类线性方程组的直接算法；3、提高分析和解决问题的能力，做到学以致用；4、通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。程序代码：1.Gsuss列主元消去法#include<iostream>#include<cmath>#include<cstdlib>using namespace std;int main()int n, i,j,k,m;cout<<"输入维数:" cin>>n;float *A=new double*(n+1);for(i=1;i<=n;i+)Ai=new doublen+1;float *b=new doublen+1;float *x=new doublen+1;float l;float temp1,temp2,temp3;cout<<"输入系数矩阵A:"<<endl;for(i=1;i<=n;i+)for(j=1;j<=n;j+)cin>>Aij;cout<<"输入向量b:"for(i=1;i<=n;i+)cin>>bi;cout<<endl; for(k=1;k<n;k+) temp1=abs(Akk);m=k;for(i=k;i<=n;i+)/找最大值的列主元 if(temp1<abs(Aik) temp1=abs(Aik);m=i;/m是确定的最列主元的行标if(temp1=0) cout<<"no unique solution!"<<endl; exit(0);if(m!=k)/换行for(j=1;j<=n;j+)temp2=Akj;Akj=Amj;Amj=temp2;temp3=bk;bk=bm;bm=temp3;for(i=k+1;i<=n;i+)/消元l=Aik/Akk;for(j=k+1;j<=n;j+)Aij=Aij-l*Akj;bi=bi-l*bk;if(Ann=0)cout<<"no unique solution!"<<endl; exit(0);xn=bn/Ann;/回代求解for(i=n-1;i>=1;i-)float sum=0; for(j=i+1;j<=10;j+) sum=sum+Aij*xj;xi=(bi-sum)/Aii;cout<<"输出结果向量x:"<<endl;for(i=1;i<=10;i+) cout<<xi<<endl;system("pause");return 0;2.平方根法#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+)Ai=new doublen+1;double *b=new doublen+1;double *x=new doublen+1;double *y=new doublen+1;cout<<"输入系数对称正定矩阵A:"<<endl;for(i=1;i<=n;i+)for(j=1;j<=n;j+)cin>>Aij;cout<<"输入向量b:"for(i=1;i<=n;i+)cin>>bi;cout<<endl;for(k=1;k<=n;k+)double sum=0;for(m=1;m<=k-1;m+)sum=sum+pow(Akm,2.0);sum=Akk-sum;Akk=sqrt(sum);for(i=k+1;i<=n;i+)double temp1=0;for(m=1;m<=k-1;m+)temp1=temp1+Aim*Akm;temp1=Aik-temp1;Aik=temp1/Akk;double temp2=0;for(m=1;m<=k-1;m+)temp2=temp2+Akm*ym;yk=(bk-temp2)/Akk;x8=y8/A88;for(k=n-1;k>=1;k-)double temp3=0;for(m=k+1;m<=n;m+)temp3=temp3+Amk*xm;xk=(yk-temp3)/Akk;cout<<"输出结果向量x:"<<endl;for(i=1;i<=n;i+) cout<<xi<<endl;system("pause");return 0;3.追赶法#include<iostream>#include<cmath>#include<cstdlib>using namespace std;int main()int n,i;cout<<"输入系数矩阵的维数:"cin>>n;double *a=new doublen+1;double *c=new doublen+1;double *d=new doublen+1;double *b=new doublen+1;double *x=new doublen+1;double *y=new doublen+1;cout<<"输入系数矩阵A数据："<<endl;for(i=1;i<=n;i+) cin>>ai;for(i=1;i<=n;i+) cin>>ci;for(i=1;i<=n;i+) cin>>di;cout<<"输入b :"<<endl;for(i=1;i<=n;i+) cin>>bi;for(i=1;i<=n-1;i+)ci=ci/ai;ai+1=ai+1-di+1*ci;cout<<"输出解向量a:"<<endl;for(i=1;i<=n;i+)cout<<ai<<endl;cout<<"输出解向量c:"<<endl;for(i=1;i<=n;i+)cout<<ci<<endl;y1=b1/a1;for(i=2;i<=n;i+)yi=(bi-di*yi-1)/ai;cout<<"输出解向量y:"<<endl;for(i=1;i<=n;i+)cout<<yi<<endl;xn=yn;for(i=n-1;i>=1;i-)xi=yi-ci*xi+1;cout<<"输出解向量x:"<<endl;for(i=1;i<=n;i+)cout<<xi<<endl;system("pause");return 0;四、 程序运行结果分析 在方法的选择上存在一定的误差，可以选一些更准确的方法求解；程序中对变量的类型设定，若设成double型，结果可以更精确；计算机在做运算时，会根据需要对中间结果进行舍入，这也会对最终结果有影响； 课题三 线性方程组的迭代法一、问题提出1、设线性方程组=x= ( 1, -1, 0, 1, 2, 0, 3, 1, -1, 2 ) J迭代算法：程序设计流程图：源程序代码：#include<stdlib.h>#include<stdio.h> #include<math.h> void main() float a5051,x150,x250,temp=0,fnum=0; int i,j,m,n,e,bk=0; printf("使用Jacobi迭代法求解方程组:n"); printf("输入方程组的元：nn="); scanf("%d",&n); for(i=1;i<n+1;i+) x1i=0; printf("输入方程组的系数矩阵:n"); for(i=1;i<n+1;i+) j=1; while(j<n+1) scanf("%f",&aij); j+; printf("输入方程组的常数项:n"); for(i=1;i<n+1;i+) scanf("%f",&ain+1); printf("n");printf("请输入迭代次数:n");scanf("%d",&m);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=aij*x1j+temp; x2i=(ain+1-temp)/aii; temp=0; for(i=1;i<n+1;i+) fnum=float(fabs(x1i-x2i); if(fnum>temp) temp=fnum; if(temp<=pow(10,-4) bk=1; for(i=1;i<n+1;i+) x1i=x2i;m-; printf("原方程组的解为:n"); for(i=1;i<n+1;i+) if(x1i-x2i)<=e|(x2i-x1i)<=e)printf("x%d=%7.4f ",i,x1i); 运行结果：2、设对称正定阵系数阵线方程组= x = ( 1, -1, 0, 2, 1, -1, 0, 2 )GS迭代算法：#include<iostream.h>#include<math.h>#include<stdio.h>const int m=11;void main() int choice=1; while(choice=1) double amm,bm,e,xm,ym,w,se,max; int n,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>>aij; cout<<"请输入各个方程等号右边的常数项:n" for(i=1;i<=n;i+) cin>>bi; cout<<"请输入最大迭代次数：" cin>>N; cout<<"请输入最大偏差：" cin>>e; for(i=1;i<=n;i+) xi=0; yi=xi; 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+aij*yj; yi=(bi-w)/double(aii); max=fabs(x1-y1); for(i=1;i<=n;i+) se=fabs(xi-yi); if(se>max) max=se; if(max<e) cout<<endl; for(i=1;i<=n;i+) cout<<"x"<<i<<"="<<yi<<endl; break; for(i=1;i<=n;i+) xi=yi; if(k=N) cout<<"迭代失败！"<<endl;choice=0;3、三对角形线性方程组 = x= ( 2, 1, -3, 0, 1, -2, 3, 0, 1, -1 )SOR方法：# include <stdio.h># include <math.h>#include<stdlib.h>/*定义全局变量*/float *a; /*存放A矩阵*/float *b; /*存放b矩阵*/float *x; /*存放x矩阵*/float p; /*精确度*/float w; /*松弛因子*/int n; /*未知数个数*/int c; /*最大迭代次数*/int k=1; /*实际迭代次数*/*SOR迭代法*/void SOR(float xk) int i,j; float t=0.0; float tt=0.0; float *xl; xl=(float *)malloc(sizeof(float)*(n+1); for(i=1;i<n+1;i+) t=0.0; tt=0.0; for(j=1;j<i;j+) t=t+aij*xlj; for(j=i;j<n+1;j+) tt=tt+aij*xkj; xli=xki+w*(bi-t-tt)/aii; t=0.0; for(i=1;i<n+1;i+) tt=fabs(xli-xki); tt=tt*tt; t+=tt; t=sqrt(t); for(i=1;i<n+1;i+) xki=xli; if(k+1>c) if(t<=p) printf("nReach the given precision!n"); else printf("nover the maximal count!n"); printf("nCount number is %dn",k); else if(t>p) k+; SOR(xk); else printf("nReach the given precision!n"); printf("nCount number is %dn",k); /*程序*开始*/void main() int i,j;printf("SOR方法n");printf("请输入方程个数:n");scanf("%d",&n);a=(float *)malloc(sizeof(float)*(n+1); for(i=0;i<n+1;i+) ai=(float*)malloc(sizeof(float)*(n+1);printf("请输入三对角矩阵:n");for(i=1;i<n+1;i+) for(j=1;j<n+1;j+) scanf("%f",&aij);for(i=1;i<n+1;i+) for(j=1;j<n;j+)b=(float *)malloc(sizeof(float)*(n+1);printf("请输入等号右边的值:n");for(i=1;i<n+1;i+) scanf("%f",&bi);x=(float *)malloc(sizeof(float)*(n+1);printf("请输入初始的x:");for(i=1;i<n+1;i+) scanf("%f",&xi);printf("请输入精确度:");scanf("%f",&p);printf("请输入迭代次数:");scanf("%d",&c);printf("请输入w(0<w<2):n");scanf("%f",&w);SOR(x);printf("方程的结果为:n");for(i=1;i<n+1;i+) printf("x%d=%fn",i,xi);四、程序运行结果讨论和分析：迭代法具有需要计算机的存贮单元较少,程序设计简单,原始系数矩阵在计算过程中始终不变等优点.迭代法在收敛性及收敛速度等方面存在问题.注:A必须满足一定的条件下才能运用以下三种迭代法之一.在Jacobi中不用产生的新数据信息,每次都要计算一次矩阵与向量的乘法,而在Gauss利用新产生的信息数据来计算矩阵与向量的乘法.在SOR中必须选择一个最佳的松弛因子,才能使收敛加速.经过计算可知Gauss-Seidel方法比Jacobi方法剩点计算量,也是Jacobi方法的改进.可是精确度底,计算量高,费时间,需要改进.SOR是进一步改进Gauss-Seidel而得到的比Jacobi,Gauss-Seidel方法收敛速度快,综合性强.改变松弛因子的取值范围来可以得到Jacobi,Gauss-Seidel方法. 选择一个适当的松弛因子是关键.结论：线性方程组1和2对于Jacobi 迭代法，Gauss-Seidol迭代法和SOR方法均不收敛，线性方程组3收敛。课题四 数值积分一、问题提出 选用复合梯形公式，复合Simpson公式，Romberg算法，计算（1） I = （2） I = （3） I = （4）（5） I = 二、要求1、 编制数值积分算法的程序；2、 分别用两种算法计算同一个积分，并比较其结果；3、 分别取不同步长，试比较计算结果（如n = 10, 20等）；4、 给定精度要求，试用变步长算法，确定最佳步长。 三、目的和意义1、 深刻认识数值积分法的意义；2、 明确数值积分精度与步长的关系；根据定积分的计算方法，可以考虑二重积分的计算问题。流程图：复化梯形：复化simpson：Romberg求积法：程序代码：复化梯形：#include<iostream.h>#include<math.h>int N;/声明被积分的函数f(x)double f(double x) return(sqrt(4-sin(x)*sin(x);/return(x=0? 1:sin(x)/x);/float e=2.718281828;/return(pow(e,x)/(4+x*x);/return(log(1+x)/(1+x*x);inline double T(double x ,double y)double sum,h,a;sum=0;h=(y-x)/N;a=x;for(int i=1;i<N;i+) /以h的大小为步长递增 x+=h; sum+=f(x); sum*=2;sum+=(f(x)+f(y);sum*=(h/2);return sum;void main()cout<<"输入你想执行的步长："cin>>N;double a ,b;cout<<"输入下界："cin>>a;cout<<"输入上界："cin>>b;/输出运行结果cout<<"经用梯形公式计算知原函数积分近似值为："<<T(a, b)<<endl; 复化simpson公式#include<iostream.h>#include<math.h>int N;double f(double x) return(sqrt(4-sin(x)*sin(x);/return(x=0? 1:sin(x)/x);/float e=2.718281828;/return(pow(e,x)/(4+x*x);/return(log(1+x)/(1+x*x);inline double S(double x, double y) double sum, hh,a;sum=0;hh=(y-x)/N;a=x;for(int i=1;i<N;i+)/以hh的大小为步长递增x+=hh;/对i的单双数进行不同处理if(i%2=1)sum+=4*f(x);elsesum+=2*f(x);sum+=(f(x)+f(y);sum*=hh/3;return sum; void main()cout<<"输入你想执行的步长："cin>>N;N*=2;/hh表示公式中h的一半double a ,b;cout<<"输入下界："cin>>a;cout<<"输入上界："cin>>b;/输出运行结果cout<<"经用Simpson公式计算知原函数的积分近似值为："<<S(a, b)<<endl;Romberg求积法1.#include <iostream.h>#include <stdio.h>#include <string.h>#include <malloc.h>#include <math.h>double my_f(double x) return(sqrt(4-sin(x)*sin(x);/return(x=0? 1:sin(x)/x);/float e=2.718281828;/return(pow(e,x)/(4+x*x);/return(log(1+x)/(1+x*x);double Romberg(double (*f)(double),double a,double b,double eps) double T64,h=b-a; int n=1; T0=h*(f(a)/4.0+f(b)/4.0+f(a+h/2.0)/2.0);/复化梯形公式 T1=h*(f(a)/6.0+f(b)/6.0+f(a+h/2.0)/1.5);/Simpson公式 for(int i=2;fabs(Ti-2-Ti-1)>eps;+i) /计算递推梯形值,base h/=2.0;n<<=1;/分点数 int j; double base=T0/h; double x=a+h/2.0; for(j=0;j<n;+j) base+=f(x); x+=h; base=base*h/2.0; /计算外推加速值,T0->Ti double k1=4.0/3.0,k2=1.0/3.0; for(j=0;j<i;+j) double hand=k1*base-k2*Tj; Tj=base; base=hand; k2=k2/(4.0*k1-k2); k1=k2+1.0; Ti=base; return Ti-1;void main() float a,b,e;cout<<"n请输入求积节点 a："<<endl; cin>>a;cout<<"n请输入求积节点 b："<<endl; cin>>b;cout<<"n请输入求积精度 e："<<endl; cin>>e;cout<<Romberg(my_f,a,b,e)<<endl;运行结果：1、Simpsonn=10n=20Romberg2、Simpsonn=10n=20Romberg3、Simpsonn=10n=20Romberg5、Simpsonn=10n=20Romberg确定最佳步长n=2n=4n=5实验心得：通过本次试验，深刻认识数值积分法的意义；明确数值积分精度与步长的关系；根据定积分的计算方法，可以考虑二重积分的计算问题。对不同步长不仅影响到运算时间运算量还影响着运算精度.,在实际操作中应根据不同的要求选取适当的步长。