1、數(shù)值積分上機實驗報告141110038桂賢進題一：數(shù)學上已經(jīng)證明了f1 4ydx = nJ。1+、2成立，所以可以通過數(shù)值積分來求71的近似值。1 .分別使用復合梯形、復合Simpson求積公式計算11的近似值。選擇不同的h,對每種求 枳公式，試將誤差刻畫為h的函數(shù)，并比較兩方法的精度。是否存在某個值，當?shù)陀谶@ 個值之后，再繼續(xù)減小h的值，計算精度不再有所改進，為什么？2 .實現(xiàn)Romberg求積方法，并重復上面的計算；3 .實現(xiàn)自適應積分方法，并重復上面的計算。解：1.1公式分析：(a)對于復合梯形公式hb - a器(D =5 Lf 3) + f 3) + 2 273 + 訪)，h=f-i=

2、l離散誤差為：,、 nh3 、 h2(b - a) ,、E(f) = -:Ga gb (2) JL4工乙(b)對于復合Simpson公式八mm-1%)= g / + f。) + 4W f(a + - l)h) + 2 fa + 2ih)(3)i=lf=l. b-a b-ah =2m n離散誤差為：Em (T)= 一尊八旬 G) = 一 /i4(i8o n)/(4) ). a b. (4)1.2 實現(xiàn)算法結(jié)果：分別利用梯形公式和Simpson公式計算結(jié)果如下：(下表中Ei(f) = n-T(f)fE2(f)=忻- S(f).此處Tl為MATLAB中的數(shù)，可以認為具有 足夠大的精度)1打T(f)E

3、i(f)S(f)E2(f)113.0000000000000000.141621/23.1000000000000000.04163.1333333333333330.008341/43.1311764705882350.010423.1415686274509802.4026e-0561/63.1369630664712640.004633.1415917809360438.7265e-0781/83.1389884944910890.002603.1415925024587071.5113e-07101/103.1399259889071590.001673.141592613939215

4、3.9651e-08121/123.1404352468468510.001163.1415926403053801.3284e-08201/203.1411759869541294.1667e-043.1415926529697856.2001e-10301/303.1414074684073301.8519e-043.1415926535353595.4434e-ll401/403.1414884869236121.0417e-043.1415926535801059.6878e-12501/503.1415259869232546.6667e-053.1415926535872532.5

5、402e-121001/1003.1415759869231291.6667e-053.1415926535897533.9968e-142001/2003.1415884869231304.1667e-063.1415926535897930從上表中可以看出：復合Simpson公式比復合梯形公式精度高，誤差收斂的速度快不少。1.3 誤差下降速度對比：誤差下降速度對比0 2 4C 4- 00 1 | 11 - 7 o o o o o o O 1 1 1 1 1 1 1-e復合梯形公式0復合Simpson公式20406080100120140160180200區(qū)間數(shù)n從上圖可以看出，復合Simp

6、son公式誤差的收斂速度比復合梯形公式的誤差的收斂速度 快不少，下面驗證收斂階。1.4 驗證收斂階：本實驗的實際誤差主要由離散誤差和計算過程中的舍入誤差組成，這里離散誤差起 主導作用，故理論上實際誤差的收斂階應該與離散誤差的收斂階相同。下面利用如下公式來il算實際的收斂階，并與理論分析所得出的離散誤差的收斂階作比較。errl /ii err 2 h2)對上表格中所列的區(qū)間長度hi值，逐次利用相鄰兩個小區(qū)間長l】i通過上述公式來計算收 斂階，并繪制成圖形。得到圖形如下：驗證收斂階i=12345E78HD1,12(a)對復合梯形公式：由上面公式(2)可知，離散誤差關(guān)于h為二階收斂，同時由上圖可知實

7、驗結(jié)果的收斂階將 近為2,故與理論分析相符。(b)對復合Simpson公式：這里卻有些奇怪，由上面公式(4)可知，離散誤差理論上為4階收斂，可實驗結(jié)果卻是將 近6階收斂。下面將進一步深入探究。探究如下：考慮這是由于被積函數(shù)f(x)的特殊性導致，而不是由于Simpson公式離散誤 差真的能達到6階收斂。由誤差的余項公式Em(/)= 一 */(9 =工 b.考慮到f(x)=1.5 計算過程中舍入誤差的影響從上表格可以看出：當區(qū)間數(shù)n從1到200時，隨著h的減小，實際誤差在減小。考慮如下問題：若不斷減小h的值，即不斷增加區(qū)間數(shù)n,是否實際誤差會一直減小？ 1.5.1理論分析：我們知道影響實驗結(jié)果的精

8、確度的因素主要有離散誤差和舍入誤差，而離散誤差的 大小可以通過離散誤差的余項來體現(xiàn)。由公式（2）和（4）,可以看出當不斷增加區(qū)間 數(shù)，即區(qū)間長度h不斷減小時，離散誤差會越來越小。但相應地由于計算精度的限制, 當h不斷減小時，舍入誤差卻會變大。故理論上會存在一個閾值H。肖h大于H時，由于離散誤差起主導作用，隨著區(qū)間長度h的減小，實際誤差會變小: 而當h小于H時，此時舍入誤差將在計算結(jié)果的精確度起主導作用，再減小h,會導致計 算結(jié)果的精確度基本保持不變甚至可能會有減小。1.5.2實驗檢驗：（a）對于復合梯形公式：注意到誤差收斂的速度較小，故我們首先選取區(qū)間數(shù)n= 10,注1,1O，1O3 .108

9、進行分析，得到下面圖形。10復合梯形公式的誤差10分段區(qū)間數(shù)nB從圖中可以看出，復合梯形公式的閾值H在區(qū)間數(shù)n為10A6到10A8之間取到。下面對 n處于區(qū)間10A6到10A8進行分析。一 - 1010復合梯形公式的誤差。一復合梯形107分段區(qū)間數(shù)n10由上圖可以看出在n取10A6 （對應h=10A_6）左右h達到閾值，故可取閾值H=10a-6.（b）對于復合Simpson公式：注意到復合Simpson公式得收斂速度較快，取i從0到100 （對應區(qū)間數(shù)m=2i）,逐次計 算出對應于m的實際誤差,并作圖如下。從圖中可以看出,在i=42（對應m=84,h=l/l68） 左右，h達到閾值，故可取閾值

10、H=l/l68。復a Simpson公式的1天差隨區(qū)間藪m的變化情況102030405060708090100i .其中分段區(qū)間數(shù)m=2i通過上述理論分析和實驗驗證，我們得到如下結(jié)果：使用復合梯形公式和復合Simpson公式計算積分值時，所分成的小區(qū)間長h都存在閾值 H,當hH時，再減小h的值，計算精度不再有所改進。原因如下：當h小于H時，此時舍入誤差將在計算結(jié)果的精確度起主導作用，再減小h, 雖然離散誤差依舊會減小，但舍入誤差會增大，這就導致計算結(jié)果的精確度基本保持不 變甚至可能會有減小。2.1公式分析對Romberg求積方法：Ti,i=y(/(a)+/W), hb-a 211b a 1Gi

11、 = 51一1，1 + %_i f(a + (k -= 2,3,.，% = -r =萬兒一】k=lTmj =上與= 2,3,.分析離散誤差為：對序列小力m= 1,2,3.,考慮1 = |，7(%)八一7/離散誤差 的數(shù)量級為0(h2J), h為序列的小區(qū)間長。(下表中E(/)=恒一 (D|，這里冗為MATLAB中的數(shù)，可以認為是足夠精確的)IruEg1130. 141621/23. 1333333333333330. 008331/43. 1421176470588235. 2499e-0441/83.1415857837618746. 8698e-0651/163. 141592665277

12、7171. 1688e-0861/323.1415926536382444. 8451e-1171/643. 1415926535897237.0166e-1481/1283. 141592653589793091/2563. 1415926535897918. 8818e-16101/5123.1415926535897928. 8818e-16111/10243. 1415926535897921. 3323e-15121/20483. 1415926535897944. 4409e-16131/40963.1415926535897921. 3323e-15141/81923. 1415

13、926535897921. 3323e-15151/163843.1415926535897914. 4409e-162.2考慮舍入誤差的影響：2.2.1 理論分析：如同1.5小節(jié)中所作分析，此處Romberg積分方法同樣存在閾值H。當h小于H時, 此時舍入誤差將在計算結(jié)果的精確度起主導作用，再減小h,雖然離散誤差依舊會減小, 但余人誤差會增大，這就導致計算結(jié)果的精確度基本保持不變甚至可能會有減小。2.2.2 實驗12:取 i=l,2,3,4,5,6,7,8,9,10,11,12,13,14,15,逐次計算EQ。,并作圖如下：從圖中可以看出，當i=9(對應h=l/256 )時，再減小h的值(即

14、增大i值)，誤差基本上 不變。故可取閾值H=l/256.3.1自適應算法實現(xiàn)：算法思想：通過遞歸調(diào)用實現(xiàn)。(此處鎮(zhèn)覺課本提供的算法雖然巧妙，但比較冗雜，且 會導致占用太多內(nèi)存。如用遞歸調(diào)用算法直觀上更容易理解。)實現(xiàn)算法如下：(自己感覺還算比較巧妙)L腳本文件zishiying.mclear;a=O;b=l;disp(,誤差限為：),e=O.OOOOO(Hh=(b-a)/2;fl=f(a);f2=f(a+b)/2);f3=f(b);ticS=h/3*(fl+f2+4*f3);t=Simpson_auto(a/b/e/Szh/2/fl/f2/f3);dispC近似值為:),tdisp(誤差為：)

15、,abs(pi-t)disp(,耗時為：力toe2.函數(shù)文件 Simpson_auto.mfunction y =Simpson_auto(A,B,e,S,h,Cl,C2,C3) %將已計算好的函數(shù)值傳遞下去，避免重復 計算fl=f(A+h);f2=flA+3*h);%利用Simpson公式分別計算左半?yún)^(qū)間和右半?yún)^(qū)間的近似值Sl=h/3*(Cl+C2+4*fl);S2=h/3*(C2+C3+4*f2);%判斷是否達到所需精度，若未達到將區(qū)間分半，進行遞歸調(diào)用if abs(S-Sl-S2)0其次由上面推導可知：無窮積分察為dx收斂。考慮用區(qū)間截去法，將無窮積分轉(zhuǎn)化為在有限區(qū)間上的積分。注意到：0

16、03I k00 x3 wdx5;00 120120=dx = .6.2899450830.103892565E- 130.3225476900.6031541041.74576110103574186924.5366202970.388879085E-19.3950709120.539294706E-340.2635603200.5217556111.4134030590.3986668113.5964257710.759424497E-17.0858100060361175868E-212.6408008440.233699724E-450.2228466040.4589646741.188

17、9321020.4170008312.9927363260.1133733825.7751435690.103991975E 19.8374674180.261017203E-315.9828739810.898547906E-6這里若取誤差容限為： = 10-7 = b = 1.2 * 109但應用mathematica求解：C0043dx = 1.138566453220937 x 10-39Jioo-Tr00 x3-dx = 1.505582131320634X 10 -18Jso ” T從這里可以看出上面的放縮是很不合理的，選取b值比實際所需值大了太多，造成可后 續(xù)計算量過大和計算結(jié)果

18、失真。下面取b=50,用求枳公式計算二言心的枳分值，并與準確值總進行比較。.復合梯形求積公式 公式如下:711hb al；(f)=2 / + / Q) + 2 2+ ih)f h = ki=l離散誤差為:En(/)= 一含/G)=l (b )/(2)(f),Q b這里 a = 0,b = 50 ,f(x) =且 f(0)=0.(b).復合Simpson求積公式hmm-1%(/)=： IfW + f(b) + 4+ (2i - 1)九)+ 2 W + 2訪)i=li=lb a b ah =2 mn離散誤差為：Em(f)= - 察-旬= 叫 Q。h-這里 a = 0,b = 50 ,f(x) =.

19、 /一1復 Simpson公式的1天差隨區(qū)間數(shù)m的變化情況1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 L其中分段反間金m=2i2O1.使用Gauss-Legendre求積公式產(chǎn)x3I=L對上面積分,先將區(qū)間由0,50變換為-1,1,令x=25*(t+l),則 f1 (f1)3已知而前6個高斯一勒讓德求積公式結(jié)點和系數(shù)為：nXkAkI0.577350269I20.774596669I00.55555555630.86113631203478548540.3399810440.65214515540.9061798460.236926885

20、0.5384693100.47862867050.9324695140.1713244920.66120938603607615730.2386191860.46791393560.9491079120.1294849660.7415311860.2797053910.405845151038183005100.417959184(d).Romberg求積公式】使用Romberg求枳序列進行計算枳分產(chǎn) x3I(0 = l可得結(jié)果如N 1hEQ。116.02734E-166.493939421/27.2333E-066.4939321731/40.1293996096.3645397941/84

21、.2871052812.2068341251/167.409694370.9157549761/326.5384733880.0445339971/646.4928298410.0011095681/1286.4939423372.9344E-0691/2566.4939394074.6484E-09101/5126.4939394028.5727E-12111/10246.4939394021.3323E-14121/20486.4939394024.4409E-15131/40966.4939394023.1086E-14141/81926.4939394022.8422E-14151/1

22、63846.4939394024.1744E44其中九與準確值之間的誤差見下圖。(e).自適應Simpson公式在區(qū)間0,50上進行積分的結(jié)果:預先設(shè)置誤差限近似值實際誤差0.10.121312585637262680.016.494986150.00104670.0016.4940241488.475E-050.00016.4939174212.198E-050.000016.4939406141.212E-061E-066.4939396612.584E-071E-076.4939396242.219E-071E-086.4939394031.145E-091E-096.4939394021337E-101E