數(shù)值逼近課程設(shè)計(jì)

1、<p><b>　　數(shù)</b></p><p><b>　　值</b></p><p><b>　　逼</b></p><p><b>　　近</b></p><p><b>　　課</b></p><p

2、><b>　　程</b></p><p><b>　　設(shè)</b></p><p><b>　　計(jì)</b></p><p><b>　　報(bào)</b></p><p><b>　　告</b></p><p>&

3、lt;b>　　作業(yè)一</b></p><p>　　多項(xiàng)式插值的Runge現(xiàn)象</p><p>　　對(duì)于Runge函數(shù)f(x)= ，在[-1,1]上作等距節(jié)點(diǎn)插值，分別取n=4,n=8,n=12,編出程序，畫出此插值的圖像。</p><p>　　程序代碼(matlab實(shí)現(xiàn))：</p><p>　　lagrange.m<

4、/p><p>　　function y=lagrange(x0,y0,x)</p><p>　　ii=1:length(x0);y=zeros(size(x));</p><p><b>　　for i=ii</b></p><p>　　ij=find(ii~=i);y1=1;</p><p>　　f

5、or j=1:length(ij),y1=y1.*(x-x0(ij(j)));</p><p><b>　　end</b></p><p>　　y=y+y1*y0(i)/prod(x0(i)-x0(ij));</p><p><b>　　end</b></p><p><b>　　rung

6、e.m</b></p><p>　　function runge(m1,m2,m3)</p><p>　　x1=-1+2*[0:m1]/m1;</p><p>　　y1=1./(1+25*x1.^2);</p><p>　　x=-1:0.01:1;</p><p>　　y4=lagrange(x1,y1,x

7、);</p><p>　　x2=-1+2*[0:m2]/m2;</p><p>　　y2=1./(1+25*x2.^2);</p><p>　　x=-1:0.01:1;</p><p>　　y5=lagrange(x2,y2,x);</p><p>　　x3=-1+2*[0:m3]/m3;</p><

8、;p>　　y3=1./(1+25*x3.^2);</p><p>　　x=-1:0.01:1;</p><p>　　y6=lagrange(x3,y3,x);</p><p>　　x=-1:0.01:1;</p><p>　　y=1./(1+25*x.^2);</p><p>　　plot(x,y,'k-

9、',x,y4,'r--',x,y5,'b-.',x,y6,'m:')</p><p>　　legend('原函數(shù)','n=4','n=8','n=12')</p><p>　　三、運(yùn)行結(jié)果（運(yùn)行過程）：</p><p>　　輸入runge(4,8,

10、12)，可得插值圖像：</p><p><b>　　作業(yè)二</b></p><p><b>　　Remez算法</b></p><p>　　求函數(shù)f(x)= 在[-1,1]上的二次多項(xiàng)式逼近。</p><p>　　參考文獻(xiàn)：數(shù)值分析算法描述與習(xí)題解答（徐士良著）</p><p&g

11、t;　　二、程序代碼(C語言實(shí)現(xiàn))：</p><p>　　#include<stdio.h></p><p>　　#include<math.h></p><p>　　double remezf(double x)</p><p><b>　　{</b></p><p>&

12、lt;b>　　double y;</b></p><p><b>　　y=exp(x);</b></p><p>　　return(y);</p><p><b>　　}</b></p><p>　　void remez(double a,double b,double p[],i

13、nt n,double eps)</p><p>　　{ int i,j,k,m;</p><p>　　double x[21],g[21],d,t,u,s,xx,x0,h,yy;</p><p>　　if(n>20) n=20;</p><p>　　m=n+1; d=1.0e+35;</p><p>　　f

14、or(k=0;k<=n;k++)</p><p><b>　　{ </b></p><p>　　t=cos((n-k)*3.1415926/(1.0*n));</p><p>　　x[k]=(b+a+(b-a)*t)/2.0;</p><p><b>　　}</b></p>&l

15、t;p>　　while(1==1){</p><p><b>　　u=1.0;</b></p><p>　　for(i=0;i<=m-1;i++){</p><p>　　p[i]=remezf(x[i]);</p><p><b>　　g[i]=-u; </b></p>&

16、lt;p><b>　　u=-u;</b></p><p><b>　　}</b></p><p>　　for(j=0;j<=n-1;j++){</p><p><b>　　k=m; </b></p><p><b>　　s=p[k-1];</b>

17、;</p><p>　　xx=g[k-1];</p><p>　　for(i=j;i<=n-1;i++){</p><p>　　t=p[n-i+j-1]; </p><p>　　x0=g[n-i+j-1];</p><p>　　p[k-1]=(s-t)/(x[k-1]-x[m-i-2]);</p>

18、<p>　　g[k-1]=(xx-x0)/(x[k-1]-x[m-i-2]);</p><p><b>　　k=n-i+j; </b></p><p><b>　　s=t; </b></p><p><b>　　xx=x0;</b></p><p><b>

19、　　}</b></p><p><b>　　}</b></p><p>　　u=-p[m-1]/g[m-1];</p><p>　　for(i=0;i<=m-1;i++)</p><p>　　p[i]=p[i]+g[i]*u;</p><p>　　for(j=1;j<=n-

20、1;j++){</p><p><b>　　k=n-j; </b></p><p><b>　　h=x[k-1];</b></p><p><b>　　s=p[k-1];</b></p><p>　　for(i=m-j;i<=n;i++){</p><

21、p>　　t=p[i-1]; </p><p>　　p[k-1]=s-h*t;</p><p><b>　　s=t; </b></p><p><b>　　k=i;</b></p><p><b>　　}</b></p><p><b>　

22、　}</b></p><p>　　p[m-1]=fabs(u); </p><p><b>　　u=p[m-1];</b></p><p>　　if(fabs(u-d)<=eps) return;</p><p><b>　　d=u; </b></p><p&g

23、t;　　h=0.1*(b-a)/(1.0*n);</p><p><b>　　xx=a; </b></p><p><b>　　x0=a;</b></p><p>　　while(x0<=b){</p><p>　　s=remezf(x0); </p><p><

24、b>　　t=p[n-1];</b></p><p>　　for(i=n-2;i>=0;i--)</p><p>　　t=t*x0+p[i];</p><p>　　s=fabs(s-t);</p><p>　　if(s>u) { u=s; xx=x0;}</p><p><b>　

25、　x0=x0+h;</b></p><p><b>　　}</b></p><p>　　s=remezf(xx);</p><p><b>　　t=p[n-1];</b></p><p>　　for(i=n-2;i>=0;i--)</p><p>　　t=t

26、*xx+p[i];</p><p>　　yy=s-t; i=1; j=n+1;</p><p>　　while((j-i)!=1){</p><p>　　k=(i+j)/2;</p><p>　　if(xx<x[k-1])j=k;</p><p><b>　　else i=k;</b><

27、;/p><p><b>　　}</b></p><p>　　if(xx<x[0]){</p><p>　　s=remezf(x[0]); </p><p><b>　　t=p[n-1];</b></p><p>　　for(k=n-2;k>=0;k--)</p&

28、gt;<p>　　t=t*x[0]+p[k];</p><p><b>　　s=s-t;</b></p><p>　　if(s*yy>0.0)x[0]=xx;</p><p><b>　　else{ </b></p><p>　　for(k=n-1; k>=0; k--)&

29、lt;/p><p>　　x[k+1]=x[k];</p><p><b>　　x[0]=xx; </b></p><p><b>　　}</b></p><p><b>　　}</b></p><p><b>　　else{</b>&l

30、t;/p><p>　　if(xx>x[n])</p><p><b>　　{</b></p><p>　　s=remezf(x[n]); </p><p><b>　　t=p[n-1];</b></p><p>　　for(k=n-2;k>=0;k--)</p&

31、gt;<p>　　t=t*x[n]+p[k];</p><p><b>　　s=s-t;</b></p><p>　　if(s*yy>0.0)x[n]=xx;</p><p><b>　　else{</b></p><p>　　for(k=0;k<=n-1;k++)<

32、/p><p>　　x[k]=x[k+1];</p><p><b>　　x[n]=xx;</b></p><p><b>　　}</b></p><p><b>　　}</b></p><p><b>　　else{</b></p

33、><p>　　i=i-1; j=j-1;</p><p>　　s=remezf(x[i]); </p><p><b>　　t=p[n-1];</b></p><p>　　for(k=n-2;k>=0;k--)</p><p>　　t=t*x[i]+p[k];</p><p&

34、gt;<b>　　s=s-t;</b></p><p>　　if(s*yy>0.0) x[i]=xx;</p><p>　　else x[j]=xx;</p><p><b>　　}</b></p><p><b>　　}</b></p><p>

35、<b>　　}</b></p><p><b>　　}</b></p><p>　　void main()</p><p>　　{ double a,b,eps,p[4];</p><p>　　a=-1.0; b=1.0; eps=1.0e-10;</p><p>　　r

36、emez(a,b,p,3,eps);</p><p>　　printf("a0=%e\n",p[0]);</p><p>　　printf("a1=%e\n",p[1]);</p><p>　　printf("a2=%e\n",p[2]);</p><p>　　printf(&quo

37、t;\n");</p><p>　　printf("MAX(p-f)=%e\n",p[3]);</p><p>　　printf("\n");</p><p><b>　　}</b></p><p>　　三、運(yùn)行結(jié)果及說明：</p><p>　　其

38、中a0,a1,a2分別為二次多項(xiàng)式零次，一次，二次的系數(shù)。</p><p><b>　　作業(yè)三</b></p><p><b>　　復(fù)化公式</b></p><p>　　用復(fù)化公式計(jì)算積分。</p><p>　　程序代碼(C++實(shí)現(xiàn))：</p><p>　　#include&

39、lt;iostream.h></p><p>　　#include<math.h></p><p>　　#include<stdlib.h></p><p>　　#include<iomanip.h></p><p>　　double f(double x) //被積分函數(shù)f(x)</p>

40、;<p><b>　　{</b></p><p>　　return exp(x);</p><p><b>　　}</b></p><p>　　int fun(int n)</p><p><b>　　{</b></p><p><b

41、>　　if(n>=0)</b></p><p><b>　　{</b></p><p>　　if(n==0)return 1;</p><p>　　else return 2*fun(n-1);</p><p><b>　　}</b></p><p>

42、　　else exit(0);</p><p><b>　　}</b></p><p>　　double T(int n) //復(fù)化梯形公式</p><p><b>　　{</b></p><p>　　double a=0,b=1,sum=0,h=(b-a)/n;</p><p

43、><b>　　long i;</b></p><p>　　for(i=1;i<n;i++)</p><p><b>　　{</b></p><p>　　sum+=f(a+i*h);</p><p><b>　　}</b></p><p>　　

44、sum+=(f(a)+f(b))/2;</p><p>　　return (h*sum);</p><p><b>　　}</b></p><p>　　double S(int n) //復(fù)化Simpson公式</p><p><b>　　{</b></p><p>　　do

45、uble Sn;</p><p>　　Sn=(4*T(2*n)-T(n))/3;</p><p>　　return (Sn);</p><p><b>　　}</b></p><p>　　double C(int n) //復(fù)化Cotes公式</p><p><b>　　{</b&

46、gt;</p><p>　　double Cn;</p><p>　　Cn=(16*S(2*n)-S(n))/15;</p><p>　　return (Cn);</p><p><b>　　}</b></p><p>　　void main()</p><p><b

47、>　　{</b></p><p><b>　　int i;</b></p><p>　　cout<<"1復(fù)化梯形公式,2復(fù)化Simpson公式,3復(fù)化Cotes公式:\n ";</p><p>　　for(i=1;i<=10;i++)</p><p><b&g

48、t;　　{</b></p><p>　　cout<<fun(i)<<" ";</p><p>　　cout<<setiosflags(ios::showpoint)<<setprecision(16);</p><p>　　cout<<T(fun(i))<<&q