福建農(nóng)林大學(xué)2008數(shù)值分析報告考卷及問題詳解

1、2008-2009 年（1）福建農(nóng)林大學(xué)機電工程學(xué)院研究生數(shù)值分析課程試卷姓名：鄭金鳳專業(yè)年級：08級機械設(shè)計及理論學(xué)號：1081302 成績: 注：以下試題編程語言任選。要求上機調(diào)試通過，提供原程序清單、計算結(jié)果。1、分別用列主元法和 Gauss-Seidel迭代法(要求 XnXn 110 4)求解方程組。10x! 2x2 x332x110x2 3x315x1 2x2 5x310(1 )列主元法#in clude<stdio.h>#in clude<math.h>Gaussdieda i(n, a,x)int n;double *a,*x;int i,j,k,r;do

2、uble m;for(k=0;k< n-1;k+)r=k;m=fabs(ak*( n+1)+k);for(i=k+1;i <n ;i+)if(fabs(ai*( n+1)+k)>m)r=i;m=fabs(ai*( n+1)+k);if(r>k)for(j=k;j <n +1;j+)m=ak*( n+1)+j;ak*( n+1)+j=ar*( n+1)+j;ar*(n +1)+j=m;for(i=k+1;i <n ;i+)m=ai*( n+1)+k/ak*( n+1)+k;for(j=k+1;j <n+1;j+)ai*( n+1)+j-=m*ak*( n

3、+1)+j;xn-1=a( n-1)*( n+1)+n/a( n-1)*( n+1)+( n-1);for(i=n-2;i>=0;i-)m=ai*( n+1)+n;for( j=i+1;j<n;j+)m-=ai*( n+1)+j*xj;xi=m/ai*( n+1)+i;main ()double a34=10,-2,-1,3,-2,10,-3,15, -1,-2,5,10,x3;int i;Gaussdiedai(3,a,x);prin tf("x1=%fn",x0);prin tf("x2=%fn",x1);prin tf("x3

4、=%fn",x2);-1B0000-735000x3 =3-330000Pre sskey 右卩 cont; xnu.>(2) Gauss-Seidel 迭代法#i nclude"stdio.h"#i nclude"math.h"#i nclude"malloc.h"GaussSeidel( n, a,b,x)int n;double *a,*b,*x;int i,j;double t,u,m=0,eps;while(1)eps=0;for(i=0;i< n;i+)t=xi;for(j=0;j< n;j+)

5、if(j!=i)m+=ai* n+j *xj;xi=(bi-m)/ai* n+i;m=0;u=xi;if(fabs(u-t)>eps)eps=fabs(u_t);if(eps<1e-4) retur n 1; mai n() int i;double a33=10,-2,-1,-2,10,-3,-1,-2,5,b3=3,15,10,x3=0,0,0;GaussSeidel(3,a,b,x);printf("%sn","解得");prin tf("x1=%fn",x0);prin tf("x2=%fn",

6、x1);prin tf("x3=%fn",x2);xl=l.179988x2=2呼甘Vx3-3_329993Pi'ess： an屮 kep ho con± Inue2、用簡單迭代法求解方程x cosx 0，要求 Xn xn 110#in elude "math.h"#i nclude "stdio.h"main ()double x1=0.0,x0;x0=x1;x1=cos(x0);while(fabs(x0-x1)>0.000001)x0=x1;x1=cos(x0);printf("%sn"

7、;,"解得");prin tf("x=%fn",x1);prin tf("n");X-0 押陰£Press any key to Gontinue3、已知函數(shù)表如下，用二次插值求y(0.85)的近似值。Xi0.500.700.800.90Yi0.100.120.130.14#in clude "math.h"double n Ig3(double x,double y,i nt n double t) i nt i,j,k,m;double z,s;z=0.0;if(*1) return(z);if(n=

8、1) z=yO;return(z);if(n=2)z=(y0*(t-x1)-y1*(t-x0)/(x0-x1); return(z);if (t<=x1) k=0;m=2;else if (t>=x n-2) k= n-3;m=n-1;elsek=1;m=n;while (m-k !=1)i=(k+m)/2;if (t<xi-1) m=i;else k=i;k=k-1;m=m-1;if(fabs(t-xk)<fabs(t-xm) k=k-1;else m=m+1;z=0.0;for(i=k;i<=m;i+)s=1.0;for(j=k；j<=m;j+)if(j

9、!=i) s=s*(t-xj)/(xi-xj);z=z+s*yi;return(z); main () double t,z;static double x4=0.50,0.70,0.80,0.90;static double y4=0.10,0.12,0.13,0.14;prin tf("n");t=0.85;z=nl g3(x,y,5,t);prin tf("x=%f,f(x)=%f n",t,z);.850J-f<x>=0.135000Press any key to continue4、已知數(shù)據(jù)表如下，試求其最小二乘二次擬合多項式。X

10、i2345678910Yi654211234#in elude "math.h"#i nclude "stdio.h"void spir(double x,double y,i nt n, double a,i nt m,double dt)int i,j,k;double z,p,c,g,q,d1,d2,s20,t20,b20;for(i=0;i<=m-1;i+) ai=0.0;if(m>n) m=n;if(m>20) m=20;z=0.0;for(i=0;i<=n-1;i+) z=z+xi/(1.0*n);bO=1.O;d1=

11、1.O* n;p=0.0;c=0.0; for(i=0;i<=n _1;i+)p=p+(xi-z);c=c+yi;c=c/d1;p=p/d1;aO=c*bO;if(m>1)t1=1.0;t0=-p; d2=0.0;c=0.0;g=0.0;for(i=0;i<=n _1;i+)q=xi-z-p;d2=d2+q*q;c=c+yi*q;g=g+(xi-z)*q*q;C=c/d2;p=g/d2;q=d2/d1;d1=d2;a1=c*t1;a0=c*t0+a0;for( j=2;j<=m_1;j+)sj=tj-1;sj-1=-p*tj-1+tj-2;if( j>=3)for

12、(k=j-2;k>=1;k-)sk=-p*tk+tk-1-q*bk;sO=-p*tO-q*bO;d2=0.0;c=0.0;g=0.0;for(i=0;i<=n _1;i+)q=s j；for(k=j-1;k>=0;k-)q=q*(xi-z)+sk; d2=d2+q*q;c=c+yi*q;g=g+(xi-z)*q*q;C=c/d2;p=g/d2;q=d2/d1;d1=d2;aj=c*sj;tj=sj;for(k=j-1;k>=0;k-)ak=c*sk+ak;bk=tk;tk=sk;main ()double x9=2,3,4,5,6,7,8,9,10,y9=6,5,4,2

13、,1,1,2,3,4,a3,dt3;int i;spir(x,y,9,a,3,dt);for (i=0; i<=2; i+)prin tf("a(%2d)=%en",i,ai);8>-1.545455e*0001>3.16666?e-0012>=2.34S4S5e-001Press any key to continue5.用Romberg方法計算積分I e dx的近似值要求誤差不超過10 -4。0#i nclude "stdio.h"#in clude "math.h" double romb(double

14、a,double b,double eps,double(*f)()int m,n,i,k;double y10,h,ep,p,x,s,q;h=b-a;y0=h*(*f) (a)+(*f) (b)/2.0; m=1; n=1; ep=eps+1.0; while (ep>=eps )&&(m <=9) p=0.0;for (i=0;i<=n-1; i+)x=a+(i+0.5)*h;p=p+(*f) (x);p=(y0+h*p)/2.0;s=1.0;for (k=1; k<=m; k+)s=4.0*s;q=(s*p-yk-1)/(s-1.0);yk-1=p

15、；p=q;ep=fabs(q-ym-1);m=m+1;ym-1=q;n=n+n;h=h/2.0;return(q);main ()double a,b,eps,t,rombf(double);a=0.0; b=1.0; eps=0.0001;t=romb(a,b,eps,rombf);prin tf("%sn","解得");prin tf("I=%f n ",t);prin tf("n");double rombf(double x)double y;y=exp(_x*x);return(y);FFWI=0.74&

16、amp;834mPresskey to continn6、用RK方法求初值問題的數(shù)值解。取步長 h=0.1 ,計算結(jié)果取8位小數(shù)。2xy1 x20x2y(0) 0#in clude <stdio.h>#in clude <math.h> double f(double x,double y) double rus;rus=x-(2*x*y)/(1+x*x);return(rus);void mai n()double x0=0; /范圍上限double x1=2;/ 范圍下限double h=0.1;/ 步長double n;/計算出的點的個數(shù)double k1,k2,

17、k3,k4;double y0=0; /用于存放計算出的常微分方程數(shù)值解int i=0;int j;n=(x1-x0)/h);printf("%sn","結(jié)果如下");for( j=O;j<n +1;j+)k仁f(xO,yO); /求 K1k2=f(x0+h/2),(y0+h*k1/2); /求 K2k3=f(x0+h/2),(y0+h*k2/2); /求 K3k4=f(x0+h),(y0+h*k3); / 求 K4 prin tf("x0=%2.1f,y=%9.8fn",x0,y0);y0=y0+(h*(k1+2*k2+2*k3+k4)/6); /求 y xO=xO+h;寶吉來交EF >c0 =H _0干x:0=0 ZL , x0 0 _ 2,>c0 =0 - 9 .V=0044039X0 =0 _ 5 -1,-2,-3,-4-5,-y=0 _11250029 _ ±S&±?6