数值分析实验报告.docx

1、数值分析实验报告 数值分析上机实验报告课题一 线性方程组的直接算法实验源程序#include stdio.h#include math.hdouble A11010=4,2,-3,-1,2,1,0,0,0,0, 8,6,-5,-3,6,5,0,1,0,0, 4,2,-2,-1,3,2,-1,0,3,1, 0,-2,1,5,-1,3,-1,1,9,4, -4,2,6,-1,6,7,-3,3,2,3, 8,6,-8,5,7,17,2,6,-3,5, 0,2,-1,3,-4,2,5,3,0,1, 16,10,-11,-9,17,34,2,-1,2,2, 4,6,2,-7,13,9,2,0,12,4,

2、 0,0,-1,8,-3,-24,-8,6,3,-1,;double A288= 4,2,-4,0,2,4,0,0, 2,2,-1,-2,1,3,2,0, -4,-1,14,1,-8,-3,5,6, 0,-2,1,6,-1,-4,-3,3, 2,1,-8,-1,22,4,-10,-3, 4,3,-3,-4,4,11,1,-4, 0,2,5,-3,-10,1,14,2, 0,0,6,3,-3,-4,2,19,;double A31010=4,-1,0,0,0,0,0,0,0,0, -1,4,-1,0,0,0,0,0,0,0, 0,-1,4,-1,0,0,0,0,0,0, 0,0,-1,4,-1,

3、0,0,0,0,0, 0,0,0,-1,4,-1,0,0,0,0, 0,0,0,0,-1,4,-1,0,0,0, 0,0,0,0,0,-1,4,-1,0,0, 0,0,0,0,0,0,-1,4,-1,0, 0,0,0,0,0,0,0,-1,4,-1, 0,0,0,0,0,0,0,0,-1,4,;double B110=5,12,3,2,3,46,13,38,19,-21,;double B28=0,-6,20,23,9,-22,-15,45,;double B310=7,5,-13,2,6,-12,14,-4,5,-5,;void lzy(double a1010,double b10,int

4、 n) int i,j,k; double l,x10,temp; for(k=0;kn-1;k+) for(j=k,i=k;jn;j+) if(j=k) temp=fabs(ajk); else if(tempfabs(ajk) temp=fabs(ajk); i=j; if(temp=0) printf(1无解!n); return; else for(j=k;jn;j+) temp=akj; akj=aij; aij=temp; temp=bk; bk=bi; bi=temp; for(i=k+1;in;i+) l=aik/akk; for(j=k;j=0;i-) temp=0; for

5、(j=i+1;jn;j+) temp=temp+aij*xj; xi=(bi-temp)/aii; for(i=0;in;i+) printf(x%d=%lft,i+1,xi); printf(n); void pfg(double a88,double b8,int n) int i,k,m; double x8,y8,temp; for(k=0;kn;k+) temp=0; for(m=0;mk;m+) temp=temp+pow(akm,2); if(akktemp) return; akk=pow(akk-temp),1.0/2.0); for(i=k+1;in;i+) temp=0;

6、 for(m=0;mk;m+) temp=temp+aim*akm; aik=(aik-temp)/akk; temp=0; for(m=0;m=0;k-) temp=0; for(m=k+1;mn;m+) temp=temp+amk*xm; xk=(yk-temp)/akk; for(i=0;in;i+) printf(x%d=%lft,i+1,xi); printf(n); void zgf(double a1010,double b10,int n) int i; double a010,c10,d10,a110,b110,x10,y10; for(i=0;in;i+) a0i=aii;

7、 if(i0) di-1=aii-1; a10=a00; for(i=0;in-1;i+) b1i=ci/a1i; a1i+1=a0i+1-di+1*b1i; y0=b0/a10; for(i=1;i=0;i-) xi=yi-b1i*xi+1; for(i=0;in;i+) printf(x%d=%lft,i+1,xi); printf(n); void main() Printf(线性方程组的解：(Gauss列主元法)n); lzy(A1,B1,10); printf(正定线性方程组的解：(平方根法)n); pfg(A2,B2,8); printf(三对角方程组的解：(追赶法)n); zgf

8、(A3,B3,10);实验结果：结果分析 从运行结果可以看出，对方程组一和方程组三解得的结果和题中提供的理论值基本相符，解误差可以忽略，但方程组二采用平方根方法的到的结果与理论值相差甚远，原因是方程组的系数矩阵A的条件数Cond（A）非常大，导致系数矩阵A是一个病态矩阵，方程组是一个病态方程组。课题三 线性方程组的迭代法实验源程序#include#includeusing namespace std;#define maxsize 2000/ Jacobi迭代法所对应的函数void Jacobi(double e,int n,double A50,double b) int k,j,i; do

9、uble sum,y; double bestf=0; double x210=0,0,0,0,0,0,0,0,0,0; for(k=0;kmaxsize;) bestf=0; for(i=0;in;i+) sum=0; for(j=0;jn;j+) sum+=Aij*x0j; x1i=x0i+(bi-sum)/Aii; y=x0i-x1i; if(bestffabs(y) bestf=fabs(y); for(i=0;in;i+) x0i=x1i; if(bestf=e) break; k+; for(i=0;in;i+) coutx0i ; coutendl; cout迭代了k次。endl

10、;/ GaussSeidol所对应的函数void GaussSeidol(double e,int n,double A50,double b) int k,j,i; double sum,y; double bestf=0; double x210=0,0,0,0,0,0,0,0,0,0; for(k=0;kmaxsize;) bestf=0; for(i=0;in;i+) sum=0; for(j=0;jn;j+) sum+=Aij*x0j; x1i=x0i+(bi-sum)/Aii; y=x0i-x1i; if(bestffabs(y) bestf=fabs(y); x0i=x1i; k

11、+; if(bestf=e) break; for(i=0;in;i+) coutx0i ; coutendl; cout迭代了k次。endl; / sor迭代法所对应的函数void sor(double w,double e,int n,double A50,double b) int k,j,i; double sum,y; double bestf=0; double x210=0,0,0,0,0,0,0,0,0,0; for(k=0;kmaxsize;) bestf=0; for(i=0;in;i+) sum=0; for(j=0;jn;j+) sum+=Aij*x0j; x1i=x0

12、i+w*(bi-sum)/Aii; y=x0i-x1i; if(bestffabs(y) bestf=fabs(y); x0i=x1i; k+; if(bestf=e) break; for(i=0;in;i+) coutx0i ; coutendl; cout迭代了k次。endl;/主函数void main() int n,i,j; double A5050; double b50; cout请输入数组维数：n; cout请输入数组数据：endl; for(i=0;in;i+) for(j=0;jAij; cout请输入右端项：endl; for(i=0;ibi;double w,e; co

13、ut请输入精度：e; cout请输入sor的松弛因子：w; coutJacobi结果：endl; Jacobi(e,n,A,b); coutGauss-Seidol结果：endl; GaussSeidol(e,n,A,b); coutSOR结果：endl; sor(w,e,n,A,b); 运行结果： ：精度同为0.001，松弛因子为1.4的结果迭代次数比松弛因子为1.5时增加了16次。松弛因子同为1.5，精度为0.0001的结果的迭代次数比精度为0.001时多了144 ： 结果分析1.体会迭代法求解线性方程组，并能与消去法做以比较消去法是一种求解的直接方法，不考虑计算过程中的舍入误差，运用此类

14、方法经过有限次算术运算就能求出线性方程组的精确解，但由于舍入误差的存在，一般求得的也是近似解。而迭代法与直接方法不同，即使无舍入误差，迭代法也难获得精确解，是一种逐次近似的方法，适于解大型稀疏线性方程组。2. 分别对不同精度要求，如由迭代次数体会该迭代法的收敛快慢 由于在程序运行结果中采用的是精度e=0.001，在下面主要是计算了e=0.0001和e=0.00001并主要针对第二个方程来说明当e=0.001时，三种不同迭代方法的迭代次数依次为1069，133，14当e=0.0001时，三种不同迭代方法的迭代次数依次为1069，570，545当e=0.00001时，三种不同迭代方法的迭代次数依次

15、为1069，1006，1079针对上述数据，可以得出随着精度要求的提高，迭代次数依次增大，收敛的越慢。3.对方程组2，3使用SOR方法时，选取松弛因子=0.8，0.9，1，1.1，1.2等，试看对算法收敛性的影响，并能找出你所选用的松弛因子的最佳者 由于=0.8在程序运行结果中，已经给出了其他几组数的运行结果 对于方程组2， =0.7时，迭代次数为18=0.8时，迭代次数为13 =0.9时，迭代次数为14=1时，迭代次数为133 =1.1时，迭代次数为197=1.2时，迭代次数为228发现w=0.8时，迭代次数最小，为13次。对于方程组3， =0.8时，迭代次数为9 =0.9时，迭代次数为8=

16、1时，迭代次数为9 =1.1时，迭代次数为9=1.2时，迭代次数为12发现=0.9时，迭代次数最小，为8次。课题六 函数插值方法5次Lagrange插值多项式#include#includeusing namespace std;#define N 5float x = 0.4, 0.55, 0.65, 0.80, 0.95, 1.05;float y = 0.41075, 0.57815, 0.69675, 0.90, 1.00, 1.25382;float p(float xx)int i,k;float pp = 0, m1, m2;for( i=0; i=N; i+ )m1 = 1;

17、m2 = 1;for( k=0; k=N; k+ )if( k != i ) m1 *= xx - xk;m2 *= xi - xk; pp += yi * m1/m2;return pp;void main() float x,y; cout请输入x,y的值：x; ciny;cout结果为：endl;coutf(x)=p(x)endl;coutf(y)=p(y)endl; 运行结果课题九 数值积分源程序#include#includeusing namespace std;/问题1所对应的梯形公式double T1(float a,float b,int n ) double x,sum;

18、double h; int k=0; h=(b-a)/n; sum=sqrt(4-sin(a)*sin(a); for(k=1;kn;k+) x=a+k*h; sum+=2*sqrt(4-sin(x)*sin(x); sum+=sqrt(4-sin(b)*sin(b); sum=(h/2)*sum; return sum;/问题1所对应的Simpson公式double S1(float a,float b,int n ) double x,sum=0; double h; int k=0; h=(b-a)/n;sum=sqrt(4-sin(a)*sin(a); for(k=1;ke) sum=

19、0; for(i=1;i=pow(2,k-1);i+) sum+=sqrt(4-sin(a+(2*i-1)*c/pow(2,k)*sin(a+(2*i-1)*c/pow(2,k); t10=0.5*t00+c*sum/pow(2,k); for(j=1;j=k;j+) t1j=(pow(4,j)*t1j-1-t0j-1)/(pow(4,j)-1); bestf=t1j-1-t0j-1; bestf=fabs(bestf); for(j=0;j=k;j+) t0j=t1j; k+; sum=t1k-2; return sum;/问题2所对应的梯形公式double T2(float a,float

20、 b,int n ) double x,sum; double h; int k=0; h=(b-a)/n; sum=1; for(k=1;kn;k+) x=a+k*h; sum+=2*sin(x)/x; sum+=sin(b)/b; sum=(h/2)*sum; return sum;/问题2所对应的Simpson公式double S2(float a,float b,int n ) double x,sum=0; double h; int k=0; h=(b-a)/n;sum=1; for(k=1;ke) sum=0; for(i=1;i=pow(2,k-1);i+) sum+=sin(

21、a+(2*i-1)*c/pow(2,k)/(a+(2*i-1)*c/pow(2,k); t10=0.5*t00+c*sum/pow(2,k); for(j=1;j=k;j+) t1j=(pow(4,j)*t1j-1-t0j-1)/(pow(4,j)-1); bestf=t1j-1-t0j-1; bestf=fabs(bestf); for(j=0;j=k;j+) t0j=t1j; k+; sum=t1k-2; return sum; /问题3所对应的梯形公式double T3(float a,float b,int n ) double x,sum; double h; int k=0; h=(b-a)/n; sum=exp(a)/(4+a*a); for(k=1;kn;k+) x=a+k*h; sum+=2*exp(x)/(4+x*x); sum+=exp(x)/(4+x*x); sum=(h/2)*sum; return sum;/问题3所对应的Simpson公式do