數(shù)值積分上機實驗報告(共17頁)

1、精選優(yōu)質(zhì)文檔-傾情為你奉上數(shù)值積分上機實驗報告 桂賢進題一：數(shù)學上已經(jīng)證明了0141+x2dx=成立，所以可以通過數(shù)值積分來求的近似值。1.分別使用復合梯形、復合Simpson求積公式計算的近似值。選擇不同的h,對每種求積公式，試將誤差刻畫為h的函數(shù)，并比較兩方法的精度。是否存在某個值，當?shù)陀谶@個值之后，再繼續(xù)減小h的值，計算精度不再有所改進，為什么？2.實現(xiàn)Romberg求積方法，并重復上面的計算；3.實現(xiàn)自適應積分方法，并重復上面的計算。解：1.1公式分析：(a) 對于復合梯形公式Tnf=h2fa+fb+2i=1n-1fa+ih ,h=b-an (1)離散誤差為：Enf=-nh312f(2

2、)=-h2b-a12f(2),a<<b (2)(b) 對于復合Simpson公式Smf=h3fa+fb+4i=1mfa+2i-1h+2i=1m-1f(a+2ih) (3)h=b-a2m=b-an離散誤差為： Emf=-mh590f(4)=-h4b-a180f(4),a<<b. (4) 1.2實現(xiàn)算法結(jié)果：分別利用梯形公式和Simpson公式計算結(jié)果如下：（下表中E1f=-Tf,E2f=-Sf. 此處為MATLAB中的數(shù),可以認為具有足夠大的精度）ihiT(f)E1(f)S(f)E2(f)113.00000.141621/23.00000.04163.33330.0083

3、41/43.82350.010423.09802.4026e-0561/63.12640.004633.60438.7265e-0781/83.10890.002603.87071.5113e-07101/103.71590.001673.92153.9651e-08121/123.68510.001163.53801.3284e-08201/203.41294.1667e-043.97856.2001e-10301/303.73301.8519e-043.53595.4434e-11401/403.36121.0417e-043.01059.6878e-12501/503.32546.666

4、7e-053.72532.5402e-121001/1003.31291.6667e-053.97533.9968e-142001/2003.31304.1667e-063.97930從上表中可以看出：復合Simpson公式比復合梯形公式精度高，誤差收斂的速度快不少。1.3誤差下降速度對比：從上圖可以看出，復合Simpson公式誤差的收斂速度比復合梯形公式的誤差的收斂速度快不少，下面驗證收斂階。1.4驗證收斂階：本實驗的實際誤差主要由離散誤差和計算過程中的舍入誤差組成，這里離散誤差起主導作用，故理論上實際誤差的收斂階應該與離散誤差的收斂階相同。下面利用如下公式來計算實際的收斂階，并與理論分析所

5、得出的離散誤差的收斂階作比較。err1err2=h1h2p對上表格中所列的區(qū)間長度 hi 值，逐次利用相鄰兩個小區(qū)間長 hi 通過上述公式來計算收斂階,并繪制成圖形。得到圖形如下：(a) 對復合梯形公式：由上面公式(2)可知，離散誤差關于h為二階收斂，同時由上圖可知實驗結(jié)果的收斂階將近為2，故與理論分析相符。(b) 對復合Simpson公式：這里卻有些奇怪，由上面公式(4)可知，離散誤差理論上為4階收斂，可實驗結(jié)果卻是將近6階收斂。下面將進一步深入探究。探究如下：考慮這是由于被積函數(shù)f(x) 的特殊性導致，而不是由于Simpson公式離散誤差真的能達到6階收斂。由誤差的余項公式Emf=-mh5

6、90f(4)=-h4b-a180f(4),a<<b.考慮到f(x)=1.5計算過程中舍入誤差的影響從上表格可以看出：當區(qū)間數(shù)n從1到200時，隨著h的減小，實際誤差在減小?？紤]如下問題：若不斷減小h的值，即不斷增加區(qū)間數(shù)n，是否實際誤差會一直減小?1.5.1理論分析：我們知道影響實驗結(jié)果的精確度的因素主要有離散誤差和舍入誤差 ，而離散誤差的大小可以通過離散誤差的余項來體現(xiàn)。由公式（2）和（4），可以看出當不斷增加區(qū)間數(shù)，即區(qū)間長度h不斷減小時，離散誤差會越來越小。但相應地由于計算精度的限制，當h不斷減小時，舍入誤差卻會變大。故理論上會存在一個閾值H。當h大于H時，由于離散誤差起主導

7、作用，隨著區(qū)間長度h的減小，實際誤差會變??；而當h小于H 時，此時舍入誤差將在計算結(jié)果的精確度起主導作用，再減小h,會導致計算結(jié)果的精確度基本保持不變甚至可能會有減小。1.5.2實驗檢驗：（a）對于復合梯形公式：注意到誤差收斂的速度較小，故我們首先選取區(qū)間數(shù) n=100，101,102,103108 進行分析，得到下面圖形。從圖中可以看出，復合梯形公式的閾值H在區(qū)間數(shù)n為106到108之間取到。下面對n處于區(qū)間106到108進行分析。由上圖可以看出在n取106 (對應h=10-6)左右h達到閾值，故可取閾值H=10-6.（b） 對于復合Simpson公式：注意到復合Simpson公式得收斂速度

8、較快，取i從0到100（對應區(qū)間數(shù)m=2i）,逐次計算出對應于m的實際誤差，并作圖如下。從圖中可以看出，在i=42(對應m=84,h=1/168)左右，h達到閾值，故可取閾值H=1/168。通過上述理論分析和實驗驗證，我們得到如下結(jié)果：使用復合梯形公式和復合Simpson公式計算積分值時，所分成的小區(qū)間長h都存在閾值H，當h<H時，再減小h的值，計算精度不再有所改進。原因如下：當h小于H 時，此時舍入誤差將在計算結(jié)果的精確度起主導作用，再減小h，雖然離散誤差依舊會減小，但舍入誤差會增大，這就導致計算結(jié)果的精確度基本保持不變甚至可能會有減小。2.1公式分析對Romberg求積方法:T1,1

9、=h12fa+fb, h1=b-aTi,1=12Ti-1,1+hi-1k=12i-1f(a+k-1hi- 1),i=2,3, hi=b-a2i-1=12hi-1Tm,j=4j-1Tm,j-1-Tm-1,j-14j-1-1,j=2,3,(mj)分析離散誤差為：對序列Tm,j，m=1,2,3 ， 考慮Em,j=abfxdx-Tm,j, 離散誤差的數(shù)量級為O(h2j)，h為序列的小區(qū)間長。(下表中Ef=-Ti,if , 這里為MATLAB中的數(shù)，可以認為是足夠精確的)IhiTi,iE(Ti,i)1130.141621/23.33330.008331/43.88235.2499e-0441/83.18

10、746.8698e-0651/163.77171.1688e-0861/323.82444.8451e-1171/643.97237.0166e-1481/1283.9793091/2563.97948.8818e-16101/5123.97928.8818e-16111/10243.97921.3323e-15121/20483.97944.4409e-16131/40963.97921.3323e-15141/81923.97921.3323e-15151/163843.97944.4409e-16 2.2考慮舍入誤差的影響：2.2.1理論分析：如同1.5小節(jié)中所作分析，此處Romberg

11、積分方法同樣存在閾值H。當h小于H 時，此時舍入誤差將在計算結(jié)果的精確度起主導作用，再減小h，雖然離散誤差依舊會減小，但舍入誤差會增大，這就導致計算結(jié)果的精確度基本保持不變甚至可能會有減小。2.2.2實驗驗證：取i=1,2,3,4,5,6,7,8,9,10,11,12,13,14,15，逐次計算E(Ti,i),并作圖如下：從圖中可以看出，當 i=9 (對應h=1/256 ) 時，再減小h的值（即增大i值），誤差基本上不變。故可取閾值H=1/256.3.1自適應算法實現(xiàn)：算法思想：通過遞歸調(diào)用實現(xiàn)。（此處感覺課本提供的算法雖然巧妙，但比較冗雜，且會導致占用太多內(nèi)存。如用遞歸調(diào)用算法直觀上更容易理

12、解。）實現(xiàn)算法如下：（自己感覺還算比較巧妙）1.腳本文件zishiying.mclear;a=0;b=1;disp('誤差限為：'),e=0.h=(b-a)/2;f1=f(a);f2=f(a+b)/2);f3=f(b);tic S=h/3*(f1+f2+4*f3);t=Simpson_auto(a,b,e,S,h/2,f1,f2,f3);disp('近似值為：'),tdisp('誤差為：'),abs(pi-t)disp('耗時為：'),toc2.函數(shù)文件Simpson_auto.mfunction y =Simpson_auto(

13、A,B,e,S,h,C1,C2,C3) % 將已計算好的函數(shù)值傳遞下去，避免重復計算f1=f(A+h);f2=f(A+3*h);% 利用Simpson公式分別計算左半?yún)^(qū)間和右半?yún)^(qū)間的近似值S1=h/3*(C1+C2+4*f1); S2=h/3*(C2+C3+4*f2); % 判斷是否達到所需精度，若未達到將區(qū)間分半，進行遞歸調(diào)用if abs(S-S1-S2)<10*e y=S1+S2;else y= Simpson_auto (A,(A+B)/2,e/2,S1,h/2,C1,f1,C2)+. Simpson_auto (A+B)/2,B,e/2,S2,h/2,C2,f2,C3); end

14、end 碰到的問題：寫的第一個遞歸函數(shù)，每次執(zhí)行遞歸函數(shù)都調(diào)用了被積函數(shù)4次。這導致了計算精度不夠，并且計算量大增。之后改為每次遞歸調(diào)用，將已計算好的函數(shù)值，當作參數(shù)直接傳遞下去，這樣每次只需計算新的節(jié)點的函數(shù)值，避免了重復調(diào)用被積函數(shù)計算舊點的函數(shù)值所帶來的問題。運行結(jié)果：預先給定的誤差限e積分的近似值實際誤差10-13.09802.40261E-0510-23.87071.51131E-0710-33.87071.51131E-0710-43.87071.51131E-0710-53.48222.36497E-0910-63.28363.69571E-1110-73.24653.01267

15、E-0910-83.30924.36701E-1010-93.24661.37327E-1010-103.08221.89715E-1110-113.90527.41E-1310-123.97903.55271E-1510-133.97921.33227E-1510-143.97928.88E-1610-153.9793010-163.979303.2計算過程中舍入誤差的影響：由于自適應Simpson算法是預先給定誤差容限，然后通過算法的執(zhí)行，保證最終得到的結(jié)果滿足誤差容限。這里無法將誤差刻畫為h的函數(shù)，也無法分析是否存在閾值H。從上面的結(jié)果表中可以看出實際誤差滿足預先給定的誤差界。題二：Pl

16、ank關于黑體輻射的理論推出積分0x3ex-1dx用所掌握的所有數(shù)值計算的方法計算這個積分，比較不同方法的計算效率和精度。1.計算該無窮積分的準確值: 推導如下：0x3ex-1dx=0x3e-x1-e-x dx=0x3e-xn=0e-nxdx=0n=0x3e-n+1xdx=n=10x3e-nxdx=n=10 t3n4e-tdt=n=11n44=4903!=4156. 2.利用數(shù)值積分的方法進行求解：2.1使用Gauss-Laguerre積分公式：區(qū)間 0,+ 中關于權函數(shù) Wx=e-x 的n+1次直交多項式為： Ln+1x=exdn+1xn+1e-xdxn+1 對于Laguerre多項式：L0

17、x=1L1=1-x L2=x2-4x+2 L3x=-x3+9x2-18x+6L4x=x4-16x3+72x2-96x+24L5x=-x5+25x4-200x3+600x2-600x+120.選取Ln+1(x)的n+1 個根作為求積基點，便得到計算積分0fxe-xdx 的n+1 點的 Gauss-Laguerre求積公式:0e-xf(x)dxk=0nAkf(xk)令 fx=x3exex-1=x31-e-x ，便可利用上式計算所要求的積分近似值。已知前五個GaussLaguerre 求積公式得節(jié)點和系數(shù)如下：nxkAk12-2(2+2)/42+2(2-2)/420.0.2.0.6.0. E-130

18、.0.1.0.4.0.E-19.0.E-340.0.1.0.3.0.E-17.0.E-212.0.E-450.0.1.0.2.0.5.0.E-19.0.E-315.0.E-6利用公式:0e-xf(x)dxk=0nAkf(xk)求得結(jié)果為：n近似值: 誤差：1 6.758 0.080226.5690.012836.8955.9623e-0446.76183.7396e-0456.93782.0201e-062.2采用截斷法首先，limx0+x3ex-1=0，故 x=0 不是瑕點。補充定義：fx=0, x= 0x3ex-1, x >0其次由上面推導可知：無窮積分 0x3ex-1dx 收斂。考

19、慮用區(qū)間截去法，將無窮積分轉(zhuǎn)化為在有限區(qū)間上的積分。注意到：bx3ex-1dxbx3x55!dx=b120x2dx=120b<對于預先給定的誤差容限 ，可取b值為： b>120 .這里若取誤差容限為：=10-7b=1.2*109但應用mathematica求解：100x3ex-1dx=1.0937×10-3950x3ex-1dx=1.0634×10-18從這里可以看出上面的放縮是很不合理的，選取b值比實際所需值大了太多，造成可后續(xù)計算量過大和計算結(jié)果失真。下面取 b=50，用求積公式計算 050x3ex-1dx 的積分值，并與準確值 415 進行比較。(a).復

20、合梯形求積公式 公式如下：Tnf=h2fa+fb+2i=1n-1fa+ih ,h=b-an 離散誤差為： Enf=-nh312f(2)=-h2b-a12f(2),a<<b 這里 a=0,b=50 ,f(x)=x3ex-1，且f(0)=0.求解結(jié)果為：(b).復合Simpson 求積公式Smf=h3fa+fb+4i=1mfa+2i-1h+2i=1m-1f(a+2ih) h=b-a2m=b-an離散誤差為： Emf=-mh590f(4)=-h4b-a180f(4),a<<b. 這里 a=0,b=50 ,f(x)=x3ex-1.(c).使用Gauss-Legendre求積公式

21、If=050x3ex-1dx對上面積分，先將區(qū)間由 0,50 變換為-1,1，令x=25*(t+1)，則If=-11254*t+13e25t+1-1dt已知而前6個高斯勒讓德求積公式結(jié)點和系數(shù)為：nxkAk1±0.12±0.100.555 555 5563±0.0.±0.0.4±0.0.±0.0.5±0.0.±0.0.±0.0.6±0.0.±0.0.±0.0.00.(d).Romberg求積公式使用Romberg求積序列進行計算積分If=050x3ex-1dx可得結(jié)果如下：i

22、hiTi,iE(Ti,i)116.02734E-166.21/27.2333E-066.31/40.6.41/84.2.51/167.0.61/326.0.71/646.0.81/1286.2.9344E-0691/2566.4.6484E-09101/5126.8.5727E-12111/10246.1.3323E-14121/20486.4.4409E-15131/40966.3.1086E-14141/81926.2.8422E-14151/163846.4.1744E-14其中Ti,i與準確值之間的誤差見下圖。(e).自適應Simpson公式在區(qū)間0,50上進行積分的結(jié)果：預先設置誤差限近似值實際