《[高等教育]东北大学数值分析实验报告.doc》由会员分享，可在线阅读，更多相关《[高等教育]东北大学数值分析实验报告.doc（33页珍藏版）》请在三一文库上搜索。
1、数值分析实验报告课题一 迭代格式的比较一、 问题提出设方程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、明确迭代收敛性与初值选取的关
2、系。程序代码：#include#includeint 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);whil
3、e(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)
4、;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，只能求得方程的一个根；由以上结果很容易发现，初值的选取对迭代敛散性有很大影响。以第一种迭代格式为例，当初值大
5、于等于-0.3时，迭代格式发散；当初值小于等于-8.3时，迭代格式也发散；只有初值在-0.3和-8.3之间时，迭代格式才收敛于1.53209。其他迭代格式也有这样的性质，即收敛于某个数值区间，超出这个区间迭代格式就是发散的，这就是所谓迭代格式的收敛性。对于不同迭代格式在不同区间具有不同的敛散性的原因，我认为可以从一下两方面理解：1、迭代法是一种逐次逼近法，其基本思想是将隐式方程归结为一组显式的计算公式，就是说，迭代过程实质上是个逐步显式化的过程。2、我们可以用几何图像来更好地理解迭代过程。由图可知，在某些区间选取的初始值随着迭代次数的增加会越来越逼近精确值，即收敛于精确值，而在另外一些区间选取
6、的初始值随着迭代次数的增加却离精确值越来越远，即不会收敛于一个确定值。课题二 线性方程组的直接算法一、问题提出给出下列几个不同类型的线性方程组，请用适当算法计算其解。 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列主元消去法；平方根法与改进平方根法；追赶法求解（选择其一）；
7、2、 应用结构程序设计编出通用程序；3、 比较计算结果，分析数值解误差的原因；4、 尽可能利用相应模块输出系数矩阵的三角分解式。三、目的和意义1、通过该课题的实验，体会模块化结构程序设计方法的优点；2、运用所学的计算方法，解决各类线性方程组的直接算法；3、提高分析和解决问题的能力，做到学以致用；4、通过三对角形线性方程组的解法，体会稀疏线性方程组解法的特点。程序代码：1.Gsuss列主元消去法#include#include#includeusing namespace std;int main()int n, i,j,k,m;coutn;float *A=new double*(n+1);f
8、or(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;jAij;cout输入向量b:;for(i=1;ibi;coutendl; for(k=1;kn;k+) temp1=abs(Akk);m=k;for(i=k;i=n;i+)/找最大值的列主元 if(temp1abs(Aik) temp1=abs(Aik);m=i;/m是确定的最列主元的行标if(te
9、mp1=0) coutno 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)coutno unique solution!=1;i-)float sum=0; for(j=i+1;j=10;j+) sum=sum+Aij*xj;xi=(bi-sum)/Aii;cout输
10、出结果向量x:endl;for(i=1;i=10;i+) coutxiendl;system(pause);return 0;2.平方根法#include#include#includeusing namespace std;int main()int n,i,j,k,m;coutn;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
11、=1;i=n;i+)for(j=1;jAij;cout输入向量b:;for(i=1;ibi;coutendl;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=1;k-)double temp3=0;for(m=k+1
12、;m=n;m+)temp3=temp3+Amk*xm;xk=(yk-temp3)/Akk;cout输出结果向量x:endl;for(i=1;i=n;i+) coutxiendl;system(pause);return 0;3.追赶法#include#include#includeusing namespace std;int main()int n,i;coutn;double *a=new doublen+1;double *c=new doublen+1;double *d=new doublen+1;double *b=new doublen+1;double *x=new doubl
13、en+1;double *y=new doublen+1;cout输入系数矩阵A数据：endl;for(i=1;iai;for(i=1;ici;for(i=1;idi;cout输入b :endl;for(i=1;ibi;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+)coutaiendl;cout输出解向量c:endl;for(i=1;i=n;i+)coutciendl;y1=b1/a1;for(i=2;i=n;i+)yi=(bi-di*yi-1)/ai;cout输出解向量y:endl;for
14、(i=1;i=n;i+)coutyi=1;i-)xi=yi-ci*xi+1;cout输出解向量x:endl;for(i=1;i=n;i+)coutxiendl;system(pause);return 0;四、 程序运行结果分析 在方法的选择上存在一定的误差，可以选一些更准确的方法求解；程序中对变量的类型设定，若设成double型，结果可以更精确；计算机在做运算时，会根据需要对中间结果进行舍入，这也会对最终结果有影响； 课题三 线性方程组的迭代法一、问题提出1、设线性方程组=x= ( 1, -1, 0, 1, 2, 0, 3, 1, -1, 2 ) J迭代算法：程序设计流程图：源程序代码：#i
15、nclude#include #include 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;in+1;i+) x1i=0; printf(输入方程组的系数矩阵:n); for(i=1;in+1;i+) j=1; while(jn+1) scanf(%f,&aij); j+; printf(输入方程组的常数项:n); for(i=1;in+1;i+) scanf(
16、%f,&ain+1); printf(n);printf(请输入迭代次数:n);scanf(%d,&m);printf(请输入迭代精度:n);scanf(%d,&e); while(m!=0) for(i=1;in+1;i+) for(j=1;jn+1;j+) if (j!=i) temp=aij*x1j+temp; x2i=(ain+1-temp)/aii; temp=0; for(i=1;itemp) temp=fnum; if(temp=pow(10,-4) bk=1; for(i=1;in+1;i+) x1i=x2i;m-; printf(原方程组的解为:n); for(i=1;in+
17、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#include#includeconst 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; coutGauss-Seidol迭代法endl;coutn; for(i=1;i=n;i+) cout请输
18、入第i个方程的各项系数：; for(j=1;jaij; cout请输入各个方程等号右边的常数项:n; for(i=1;ibi; coutN; coute; 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;imax) max=se; if(maxe) coutendl; for(i=1;i=n;i+) coutxi=yiendl; br
19、eak; 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 # include #include/*定义全局变量*/float *a; /*存放A矩阵*/float *b; /*存放b矩阵*/float *x; /*存放x矩阵*/float p; /*精确度*/float w; /*松弛因子*/int n; /*未知数个数*/int c; /*最大迭代次数*/int k=1; /*实际迭代次数*/*SO
20、R迭代法*/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;in+1;i+) t=0.0; tt=0.0; for(j=1;ji;j+) t=t+aij*xlj; for(j=i;jn+1;j+) tt=tt+aij*xkj; xli=xki+w*(bi-t-tt)/aii; t=0.0; for(i=1;in+1;i+) tt=fabs(xli-xki); tt=tt*tt; t+=tt; t=sqrt(t); f
21、or(i=1;ic) if(tp) 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;in+1;i+) ai=(float*)malloc(sizeof(float)*(n+1);printf(请输入三对角矩阵:n);for(i=
22、1;in+1;i+) for(j=1;jn+1;j+) scanf(%f,&aij);for(i=1;in+1;i+) for(j=1;jn;j+)b=(float *)malloc(sizeof(float)*(n+1);printf(请输入等号右边的值:n);for(i=1;in+1;i+) scanf(%f,&bi);x=(float *)malloc(sizeof(float)*(n+1);printf(请输入初始的x:);for(i=1;in+1;i+) scanf(%f,&xi);printf(请输入精确度:);scanf(%f,&p);printf(请输入迭代次数:);scanf
23、(%d,&c);printf(请输入w(0w2):n);scanf(%f,&w);SOR(x);printf(方程的结果为:n);for(i=1;in+1;i+) printf(x%d=%fn,i,xi);四、程序运行结果讨论和分析：迭代法具有需要计算机的存贮单元较少,程序设计简单,原始系数矩阵在计算过程中始终不变等优点.迭代法在收敛性及收敛速度等方面存在问题.注:A必须满足一定的条件下才能运用以下三种迭代法之一.在Jacobi中不用产生的新数据信息,每次都要计算一次矩阵与向量的乘法,而在Gauss利用新产生的信息数据来计算矩阵与向量的乘法.在SOR中必须选择一个最佳的松弛因子,才能使收敛加速
24、.经过计算可知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
25、） I = （3） I = （4）（5） I = 二、要求1、 编制数值积分算法的程序；2、 分别用两种算法计算同一个积分，并比较其结果；3、 分别取不同步长，试比较计算结果（如n = 10, 20等）；4、 给定精度要求，试用变步长算法，确定最佳步长。 三、目的和意义1、 深刻认识数值积分法的意义；2、 明确数值积分精度与步长的关系；根据定积分的计算方法，可以考虑二重积分的计算问题。流程图：复化梯形：复化simpson：Romberg求积法：程序代码：复化梯形：#include#includeint N;/声明被积分的函数f(x)double f(double x) return(sqrt(
26、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;iN;i+) /以h的大小为步长递增 x+=h; sum+=f(x); sum*=2;sum+=(f(x)+f(y);sum*=(h/2);return sum;void main()coutN;double a
27、 ,b;couta;coutb;/输出运行结果cout经用梯形公式计算知原函数积分近似值为：T(a, b)endl; 复化simpson公式#include#includeint 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-
28、x)/N;a=x;for(int i=1;iN;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()coutN;N*=2;/hh表示公式中h的一半double a ,b;couta;coutb;/输出运行结果cout经用Simpson公式计算知原函数的积分近似值为：S(a, b)endl;Romberg求积法1.#include #include #include #include #include doub
29、le 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
30、)/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;jTi double k1=4.0/3.0,k2=1.0/3.0; for(j=0;ji;+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
31、 main() float a,b,e;coutn请输入求积节点 a：a;coutn请输入求积节点 b：b;coutn请输入求积精度 e：e;coutRomberg(my_f,a,b,e)endl;运行结果：1、Simpsonn=10n=20Romberg2、Simpsonn=10n=20Romberg3、Simpsonn=10n=20Romberg5、Simpsonn=10n=20Romberg确定最佳步长n=2n=4n=5实验心得：通过本次试验，深刻认识数值积分法的意义；明确数值积分精度与步长的关系；根据定积分的计算方法，可以考虑二重积分的计算问题。对不同步长不仅影响到运算时间运算量还影响着运算精度.,在实际操作中应根据不同的要求选取适当的步长。