計算方法實驗報告

1、計算方法實驗報告（三）（一）數(shù)值微分法一、實驗問題已知函數(shù)y=f(x)的數(shù)值如下，利用兩點和三點數(shù)值微分公式計算f'1和f”(1)。x0.60.80.911.11.11.4y0.73600.83650.909511.11051.24461.6017二、問題的分析（描述算法的步驟等）基本的數(shù)值微分公式：fkxLnkx=i=0nlikxfxi其中節(jié)點等距（xi=x0+ih）且x為節(jié)點的特殊情形兩點數(shù)值微分公式：f'(x0)=(y1-y0)/hf'(x1)=(y1-y0)/h三點數(shù)值微分公式：f'(x1)=(y2-y0)/2h二階數(shù)值微分公式：f''(

2、x1)=(y2-2y1+y0)/h2三、程序設(shè)計兩點數(shù)值微分公式f0.m,f1.mf0.m:function y0 = f0( x1,x2,h )y0=(x1-x2)/h;endf1.m:function y1 = f1( x1,x2,h )y1=(x1-x2)/h;end三點數(shù)值微分公式f11.mfunction y11 = f11( x2,x0,h )y11=(x2-x0)/(2*h);end二階數(shù)值微分公式f2.mfunction yy = f2( x2,x1,x0,h )yy=(x2-2*x1+x0)/(h*h);end得出f'1和f”(1)結(jié)果的程序weifen.m：x=0.

3、6 0.8 0.9 1 1.1 1.2 1.4 ;y=0.7360 0.8365 0.9095 1 1.1105 1.2446 1.6017;for i=1:3 h=0.1*2(i-1); d(1)=f0(y(4+i),y(4),h); d(2)=f1(y(4),y(4-i),h); d(3)=f11(y(4+i),y(4-i),h); d(4)=f2(y(4+i),y(4),y(4-i),h); fprintf('h=%.1f時：',h); d=vpa(d,5)end四、計算結(jié)果五、結(jié)果分析隨著h的增大，數(shù)值微分公式的截斷誤差越來越大。同時，在兩點誤差公式中很明顯計算f

4、9;(x0)的誤差與f'(x1)的誤差的絕對值基本相等，同時在三點數(shù)值微分公式中求f'x1的誤差在三個求一階導數(shù)的公式中最小。六、實驗的總結(jié)與體會數(shù)值微分公式的截斷誤差是隨步長h的減小而減小的，但由于近似式本身均用h作除數(shù)，所以函數(shù)值計算的舍入誤差，將隨h的減小而給導數(shù)的計算帶來越來越大的誤差。因此h的選取不能過小，原則上不能使舍入誤差超過截斷誤差。在求解微分方程時，很明顯三點數(shù)值微分方程結(jié)果更加精確。（二）變步長積分法一、實驗問題用變步長辛浦生積分公式和龍貝格積分法計算下列積分。（1）02sinxxdx (2)01dx1+x二、問題的分析（描述算法的步驟等）（1）變步長辛浦生

5、積分法步驟，要求誤差小于令h=b-a2,s0fa+fb,s1fa+h s20,ss0+4s1h3,n1L:令hh2,n2n,s2s1+s2,s10 對i=1n 令s1s1+fa-h+2ih 令s*s0+4s1+2s2h3,Ds*-s15,ss*若D> 則轉(zhuǎn)L，否則輸出s并停止計算（2）龍貝格積分法步驟，要求誤差小于：1令hb-a,T0fa+fbh2,n122對i=110做令hh2,n2n,Ti0對k=1n 令TiTi+fa+2ih-h令TiTi-12+Tih,c4對j=1i做令D(Ti+1-j-Ti-j)/(c-1)令Ti-jTi+1-j+D,c4c若則輸出i-j并停止計算3打印失敗信息

6、，輸出T0并停止計算三、程序設(shè)計變步長辛浦生公式的函數(shù)Simpson.mfunction Svalue = Simpson( a,b,m )h=(b-a)/2;ss=f(a)+f(b);sss=f(a+h);ssss =0;s=(ss+4*sss)*(h/3);n=1;while m>=nh=h/2;n=2*n;ssss =sss+ssss; sss=0;for i= 1:n; sss=sss+f(a-h+2*i*h);ends=(ss+4*sss+2*ssss)*(h/3);endSvalue = vpa(s,7);end龍貝格積分法Romberg.mfunction Rvalue =

7、 Romberg (a,b,m)h=b-a;R(1,1)=(ff(a)+ff(b)*(h/2);n=1;j=0;while m>=(n/2) j=j+1; h=h/2; T=0; for i=1:n x=a+h*(2*i-1); T=T+ff(x); end R(j+1,1)=R(j,1)/2+h*T; for k=1:j R(j+1,k+1)=R(j+1,k)+(R(j+1,k)-R(j,k)/(4k-1); end n=2*n;endR=vpa(R,6)Rvaule=vpa(R(j+1,j+1),6)nend 第（1）題中的函數(shù)f.mfunction y =f(x )y=sin(x)/x;if x = 0y=1;end end 第（2）題中的函數(shù)ff.mfunction yy = ff( x )yy=1/(1+x);end調(diào)用時的調(diào)用語句：（1）Svalue=Simpson(0,pi/2,3)Rvalue=Romberg(0,pi/2,3)（2）Svalue=Simpson(0,1,3)Rvalue=Romberg(0,1,3)四、計算結(jié)果（1）辛浦生公式的結(jié)果如下：龍貝格求積公式結(jié)果如下：（2）辛浦生公式的結(jié)果如下：龍貝格求積公式結(jié)果如下：五、結(jié)果分析由運行結(jié)果可知，采用龍貝格