多重積分的蒙特卡羅算法編程

1、第23卷 第4期 湖 南 文 理 學(xué) 院 學(xué) 報(bào)(自 然 科 學(xué) 版 Vol. 23 No. 4 2011年12月 Journal of Hunan University of Arts and Science(Natural Science Edition Dec. 2011多重積分的蒙特卡羅算法編程任明慧, 劉麗芳, 梅漢飛(湖南文理學(xué)院 數(shù)學(xué)與計(jì)算科學(xué)學(xué)院, 湖南 常德, 415000摘 要: 基于積分重?cái)?shù)的增加, MATLAB 中計(jì)算積分的基本庫(kù)函數(shù)算法的復(fù)雜程度不斷增加, 運(yùn)用MATLAB 的向量化編程和蒙特卡羅(Monte Carlo算法, 實(shí)現(xiàn)了多重積分的快速計(jì)算. 關(guān)鍵詞: 向

2、量化編程; 蒙特卡羅算法; 多重積分.中圖分類號(hào): O 241 文章編號(hào): 16726146(201104000102Editing program with Monte-Carlo method for calculating multiple integralREN Ming hui, LIU Li fang, MEI Han fei(College of Mathematics and Computer Science, Hunan University of Art and Science, Changde 415000, ChinaAbstract: Integral approxi

3、mate calculation is widely used in many practical projects. This paper achieved computing quickly multiple integral by applying vectoring program in MATLAB and with Monte Carlo method. Key words: vectoring program; Monte-Carlo method; multiple integral收稿日期: 20111012基金項(xiàng)目: 國(guó)家自然科學(xué)基金資助(60872129; F010401

4、蒙特卡羅算法12是計(jì)算多重積分的一種重要且有效的方法, 它可以計(jì)算被積函數(shù)非常復(fù)雜的定積分、多重積分, 并且維數(shù)沒有限制, 雖然計(jì)算量比較大, 給出的計(jì)算精度不很高, 但隨著計(jì)算機(jī)的發(fā)展, 快速得到一個(gè)有價(jià)值的結(jié)果是可行的. 本文運(yùn)用數(shù)學(xué)軟件MATLAB, 通過向量化的編程, 讓蒙特卡羅算法能快速地計(jì)算出一般區(qū)域上的多重積分.1 MATLAB 中計(jì)算積分的基本庫(kù)函數(shù)1.1 單變量積分?jǐn)?shù)值計(jì)算中人們常用矩形公式、梯形公式及辛 普森公式來計(jì)算積分, 如單變量積分(d ba f x x , 將, a b 分為n 個(gè)小區(qū)間, 在1, i i x x +上的積分i f 的 近似值分別為3:梯形公式: 1

5、(2i i i i f x f x f h +; (1辛普森公式: (4(4ii i i h f f x f x + 32(4(2412i i i i i i i h h hf x f x f x h +. (2式中1i i i h x x +=. MATLAB 基于上算法, 采用自適應(yīng)步長(zhǎng)的方法, 分別給出了庫(kù)函數(shù)trapz(, x y 和quad(fun, , a b , 這里, x y 分別是自變量的列向量和相應(yīng)被積函數(shù)列向量, fun 是被積函數(shù)名, 或用inline(函數(shù)直接定義. 雖然上述近似計(jì)算能得到比較滿意的結(jié)果, 但隨著積分重?cái)?shù)的增加, 該算法的計(jì)算量顯著增加, 以至于計(jì)算機(jī)

6、也很難完成, 另外算法的復(fù)雜程度要比蒙特卡羅更困難. 1.2 多重積分對(duì)二重積分, MATALB 有辛普森算法的庫(kù)函數(shù)dblquad(fun,a b c d , 但只能對(duì)矩形區(qū)域進(jìn)行積分;對(duì)一般區(qū)域的積分, 最近美國(guó)學(xué)者Howard Wilson 與Bryce Gardner 開發(fā)的數(shù)值積分工具箱中, 庫(kù)函數(shù)quad2dggen 解決了這個(gè)問題, 但該函數(shù)在最近的MATLAB 版本中還不能使用; 符號(hào)積分的庫(kù)函數(shù)2 湖 南 文 理 學(xué) 院 學(xué) 報(bào)(自 然 科 學(xué) 版 2011年 int(fun,a b (單變量積分, 可通過化二重積分(,d d Df x y x y 為2個(gè)一次積分21(d (

7、,d bx ax x f x y y ,運(yùn)用1int(int(, , (,f x y y x 2(, , x a b 來實(shí)現(xiàn), 但該被積函數(shù)的內(nèi)積分必須是可積的.2 蒙特卡羅算法蒙特卡羅方法是指任何利用隨機(jī)數(shù)序列來做隨機(jī)模擬的數(shù)值方法. 蒙特卡羅算法是利用服從某一種概率分布的隨機(jī)數(shù)計(jì)算數(shù)值積分, 有2種方法: 隨機(jī)投點(diǎn)法和均值估計(jì)法, 均值估計(jì)法比隨機(jī)投點(diǎn)法更好操作, 精度(若選擇的密度函數(shù)合適更高, 本文討論均值估計(jì)法.考慮k 維積分1, 45:1212(,d d d k k I f x x x x x x =""",其中為k 維積分區(qū)域.選取V 上的一個(gè)概率密

8、度函數(shù)12(,g x x "k x , (: , 1,2,i i i V V a x b i k =", 且:1212(,0, (,k k g x x x x x x "",故:1211212(,(,d d d (,k k k k f x x x I x x x x xg x x x ="""""121212121(,(,1E(,(,m k i i ik k i i ik i f x x x f x x x g x x x m g x x x ="""". (3 這

9、里m 是隨機(jī)點(diǎn)落在區(qū)域中的投點(diǎn)數(shù).特別對(duì)二重積分:12(,d ,:, (f x y a x b c x y x d , 若選取g 是上的均勻分布, 記,V a b c d =×, 由幾何概率得V S mS n, ,V S S 分別是,V 的測(cè)度,于是:11(,(,mmVi i iii i SSI f x y f x y mn=. (43 多重積分的蒙特卡羅算法MATL AB 編程實(shí)現(xiàn)以二重積分為例, 運(yùn)用算法(4編程, 一個(gè)是通常的for 循環(huán)編程, 文件名為mtc.m, 另一個(gè)是向量化編程, 文件名為mtcx.m, 后者運(yùn)行速度要快很多.function s=mtc(f , fai

10、1, fai2, a , b , c , d , n if nargin<8 n =10 000; endx =unifrnd(a , b , 1, n ; y =unifrnd(c , d , 1, n ; s =0;for k =1: nif and(y (k >=feval(fai1, x(k, y (k <=feval (fai2, x (k s =s +feval(f , x (k , y (k ;end ends =s/n*(b a *(d c ; returnfunction s =mtcx(f , fai1, fai2, a , b , c , d , n i

11、f nargin<8 n =10 000; endx =unifrnd(a , b , 1, n ; y =unifrnd(c , d , 1, n ;s =sum(feval(f , x , y .*and(y >=feval(fai1, x , y <= feval(fai2, x s=s/n*(b a *(d c ; return三重或三重以上積分只須增加若干個(gè)輸入?yún)?shù), 即新增積分變量的上、下限, 方法不變, 對(duì)程序稍作修改即可.例計(jì)算1212sin(d d x y y x +.積分的精確值為0.411 929 546 176 30, 在命令窗口調(diào)用函數(shù)M 文件xtc

12、.m 和xtcx.m 情況如下:>>tic; s 1=mtc(inline(exp(x .2/2.*sin(x .2+y ,inline(sqrt(1x .2/2, inline(sqrt(1x .2/ 2, 0.5, 1, 1, 1, 100 000, toc s 1 =0.409 010 343 471 83Elapsed time is 40.094 000 seconds.>>tic; s 2=mtcx(inline(exp(x .2/2.*sin(x .2+ y , inline(sqrt(1x .2/2, inline(sqrt(1 x .2/ 2, 0.5

13、, 1, 1, 1, 100 000, toc s 2 =0.411 839 625 630 18Elapsed time is 0.047 000 seconds.在運(yùn)用了向量化編程后運(yùn)行時(shí)間從40.094 s 縮短到0.047 s, 是原來的1/853, 效果好; 誤差為 0.000 089 920 546 120 003 66, 也比前者好.若將上面程序稍作如下修改, 即輸入?yún)?shù)減少2個(gè), 那么當(dāng)n 為1百萬時(shí), 運(yùn)行時(shí)間大約增加0.45 s, 精度不變.function s =mtcxg(f , fai1, fai2, a , b , n if nargin<6 n =10 00

14、0; end (下轉(zhuǎn)第6頁(yè)6 湖 南 文 理 學(xué) 院 學(xué) 報(bào)(自 然 科 學(xué) 版 2011年b. 7n G S ×包含一個(gè)子圖同胚于n H , 所以 7(5,2n n n cr G S cr H Z n ×=+(. 綜上所述有7(5,2n n cr G S Z n ×=+(.參考文獻(xiàn):1 Bondy J A, Murty M S R. Graph Theory withApplicationsM. New York: Am Elsvier, 1976.2 Garary M R, Johnson D S. Crossing number is NP complete

15、J. SIAM J Algebric Discrete Methods, 1993, 4: 312316.3 Kle M. The crossing numbers of 2,3n K P × and2,3n K S ×J. Tatra Mountain Math Publ, 1996, 9: 5156.4 Kle M. The crossing number of Cartesian products ofpaths with 5vertex graphsJ. Discrete Math, 2001, 233: 353359.5 Kle M. The crossing n

16、umbers of paths and stars with4-vertex graphsJ. J Graph Theory, 1994, 18: 605614. 6 Kle M. The crossing numbers of 5nK P ×J. Tatra Mountain Math Publ, 1999, 18: 6368.7 Beineke L W, Ringeisen R D. On the crossing number ofcycles and graphs of order fourJ. J Graph indent Theory, 1980, 4: 145155.8

17、 Kle M. The crossing number of certain CartesianproductsJ. Discuss Math Graph Theory, 1995, 15: 510. 9 Dean A M, Richter R B. The crossing number of 44C C ×J. J Graph Theory, 1995, 19: 125129.10 Kle M, Richter R B, Stobert Ian. The crossing numbersof 4n C C ×J. J Graph Theory, 1996, 22: 23

18、9243. 11 肖文兵, 黃元秋. 一類笛卡兒積圖的交叉數(shù)J. 湖南師范大學(xué)學(xué)報(bào): 自然科學(xué)版, 2003, 26(4: 37, 17. 12 張莉茜, 李波, 黃元秋. 關(guān)于六階圖與星的笛卡兒積圖交叉數(shù)J. 湖南文理學(xué)院學(xué)報(bào): 自然科學(xué)版, 2008, 20(1: 1619.13 周智勇, 肖文兵, 黃元秋. 星圖n S 及5個(gè)六階圖與路的笛卡兒積圖的交叉數(shù)J. 湖南文理學(xué)院學(xué)報(bào): 自然科學(xué)版, 2007, 19(1: 14.14 Woodall D R. Cyclic graphs and Zarankiewicz´s crossingnumber conjectureJ. J

19、 Combin Theory, 1993, 17: 657 671.15 Kleitman D J. The crossing number of 5,n K J. J Comb inatorial Theory, 1970, 9: 315323.16 Kouhei Asano. The crossing number of 1,3,n K and2,3,n K J. J Graph Theory, 1980, 10: 18.(責(zé)任編校: 劉曉霞(上接第2頁(yè)x =unifrnd(a , b , 1, n ;c =min(feval(fai1, x ;d =max(feval(fai2, x ; y =unifrnd(c