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

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

2、量化編程; 蒙特卡羅算法; 多重積分.中圖分類號: O 241 文章編號: 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)目: 國家自然科學(xué)基金資助(60872129; F010401

4、蒙特卡羅算法12是計算多重積分的一種重要且有效的方法, 它可以計算被積函數(shù)非常復(fù)雜的定積分、多重積分, 并且維數(shù)沒有限制, 雖然計算量比較大, 給出的計算精度不很高, 但隨著計算機(jī)的發(fā)展, 快速得到一個有價值的結(jié)果是可行的. 本文運(yùn)用數(shù)學(xué)軟件MATLAB, 通過向量化的編程, 讓蒙特卡羅算法能快速地計算出一般區(qū)域上的多重積分.1 MATLAB 中計算積分的基本庫函數(shù)1.1 單變量積分?jǐn)?shù)值計算中人們常用矩形公式、梯形公式及辛 普森公式來計算積分, 如單變量積分(d ba f x x , 將, a b 分為n 個小區(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)步長的方法, 分別給出了庫函數(shù)trapz(, x y 和quad(fun, , a b , 這里, x y 分別是自變量的列向量和相應(yīng)被積函數(shù)列向量, fun 是被積函數(shù)名, 或用inline(函數(shù)直接定義. 雖然上述近似計算能得到比較滿意的結(jié)果, 但隨著積分重數(shù)的增加, 該算法的計算量顯著增加, 以至于計算機(jī)

6、也很難完成, 另外算法的復(fù)雜程度要比蒙特卡羅更困難. 1.2 多重積分對二重積分, MATALB 有辛普森算法的庫函數(shù)dblquad(fun,a b c d , 但只能對矩形區(qū)域進(jìn)行積分;對一般區(qū)域的積分, 最近美國學(xué)者Howard Wilson 與Bryce Gardner 開發(fā)的數(shù)值積分工具箱中, 庫函數(shù)quad2dggen 解決了這個問題, 但該函數(shù)在最近的MATLAB 版本中還不能使用; 符號積分的庫函數(shù)2 湖 南 文 理 學(xué) 院 學(xué) 報(自 然 科 學(xué) 版 2011年 int(fun,a b (單變量積分, 可通過化二重積分(,d d Df x y x y 為2個一次積分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ù)計算數(shù)值積分, 有2種方法: 隨機(jī)投點(diǎn)法和均值估計法, 均值估計法比隨機(jī)投點(diǎn)法更好操作, 精度(若選擇的密度函數(shù)合適更高, 本文討論均值估計法.考慮k 維積分1, 45:1212(,d d d k k I f x x x x x x =""",其中為k 維積分區(qū)域.選取V 上的一個概率密

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ù).特別對二重積分: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 的測度,于是:11(,(,mmVi i iii i SSI f x y f x y mn=. (43 多重積分的蒙特卡羅算法MATL AB 編程實(shí)現(xiàn)以二重積分為例, 運(yùn)用算法(4編程, 一個是通常的for 循環(huán)編程, 文件名為mtc.m, 另一個是向量化編程, 文件名為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三重或三重以上積分只須增加若干個輸入?yún)?shù), 即新增積分變量的上、下限, 方法不變, 對程序稍作修改即可.例計算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)行時間從40.094 s 縮短到0.047 s, 是原來的1/853, 效果好; 誤差為 0.000 089 920 546 120 003 66, 也比前者好.若將上面程序稍作如下修改, 即輸入?yún)?shù)減少2個, 那么當(dāng)n 為1百萬時, 運(yùn)行時間大約增加0.45 s, 精度不變.function s =mtcxg(f , fai1, fai2, a , b , n if nargin<6 n =10 00

14、0; end (下轉(zhuǎn)第6頁6 湖 南 文 理 學(xué) 院 學(xué) 報(自 然 科 學(xué) 版 2011年b. 7n G S ×包含一個子圖同胚于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é)報: 自然科學(xué)版, 2003, 26(4: 37, 17. 12 張莉茜, 李波, 黃元秋. 關(guān)于六階圖與星的笛卡兒積圖交叉數(shù)J. 湖南文理學(xué)院學(xué)報: 自然科學(xué)版, 2008, 20(1: 1619.13 周智勇, 肖文兵, 黃元秋. 星圖n S 及5個六階圖與路的笛卡兒積圖的交叉數(shù)J. 湖南文理學(xué)院學(xué)報: 自然科學(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頁x =unifrnd(a , b , 1, n ;c =min(feval(fai1, x ;d =max(feval(fai2, x ; y =unifrnd(c