張 濤,周仲禮,安素珍,張 軍
(成都理工大學(xué)數(shù)學(xué)地質(zhì)四川省重點實驗室,四川成都 610059)
關(guān)于計算多重積分的擬蒙特卡羅方法
張 濤,周仲禮,安素珍,張 軍
(成都理工大學(xué)數(shù)學(xué)地質(zhì)四川省重點實驗室,四川成都 610059)
介紹了擬蒙特卡羅方法計算多重積分的基本原理,對Halton序列和Rand函數(shù)產(chǎn)生的序列的均勻性進(jìn)行了比較,并給出了擬蒙特卡羅方法計算多重積分的步驟、實例、和Matlab語言編寫的計算程序.實例充分體現(xiàn)了采用擬蒙特卡羅方法計算多重積分的有效性、精確性和優(yōu)越性.
擬蒙特卡羅方法;Halton序列;多重積分
在數(shù)學(xué)和工程科學(xué)計算中,求解多重積分的近似值是一個至關(guān)重要的環(huán)節(jié).通常人們所采用的方法有3種,即蒙特卡羅方法、擬蒙特卡羅方法和數(shù)論網(wǎng)格方法.然而,在實際應(yīng)用中,數(shù)論網(wǎng)格法很難解決高維問題,如對20維這樣的問題,就要有一百萬多個點[1],計算量大體上隨維數(shù)的冪次增加,幾乎是不可能計算的.蒙特卡羅方法是采用單和來得到多重積分的近似值的[2],計算跟問題的維數(shù)無關(guān),克服了“維數(shù)災(zāi)難”,但其誤差卻是概率性的,得到解的精度不高.本文研究了擬蒙特卡羅方法在計算多重積分中的應(yīng)用,并給出Matlab程序的實現(xiàn).實例表明,采用擬蒙特卡羅方法計算多重積分,既能克服維數(shù)的制約,又能得到確定性的誤差估計和更高精度的解.
擬蒙特卡羅方法也稱低差異序列法,是在蒙特卡羅方法基礎(chǔ)上發(fā)展起來的一種模擬方法.?dāng)M蒙特卡羅方法與蒙特卡羅方法相似,但所用的理論基礎(chǔ)卻不同.?dāng)M蒙特卡羅方法是通過構(gòu)造所謂的低差異序列,即用的是確定性的點,然后用Koksma-Hlawka不等式來確定誤差階的,而不是根據(jù)大數(shù)定律[3].
卡羅方法模擬解積分問題的誤差界可以表示為[3-4]:
對于擬蒙特卡羅方法積分,我們有確定性的誤差估計.在計算積分中還可以明智地選擇這些確定的點來減小式(1)中的誤差,得到高精度的解.因此,就確定性和高精度性這兩點,擬蒙特卡羅積分要優(yōu)于蒙特卡羅積分.對于更高維的積分問題,擬蒙特卡羅法優(yōu)于數(shù)論網(wǎng)格法.
在過去的幾十年里,擬蒙特卡羅法得到了快速發(fā)展,人們已經(jīng)構(gòu)造出大量優(yōu)質(zhì)的低差異序列,如Halton序列、Sobol'序列、Niderreiter的(t,m,s)網(wǎng)格和(t,s)序列等[6].本文主要研究Halton序列在計算多重積分中的應(yīng)用.
下面將給出產(chǎn)生Halton序列的Matlab程序[7],并對Halton序列與Rand函數(shù)產(chǎn)生的隨機(jī)數(shù)序列的均勻性進(jìn)行比較,見圖1.
圖1 Halton序列與Rand函數(shù)產(chǎn)生的十維隨機(jī)數(shù)點的均勻性比較
由圖1可明顯看出,Rand函數(shù)產(chǎn)生的偽隨機(jī)數(shù)序列有抱團(tuán)現(xiàn)象,分布不均勻,而Halton序列則要均勻得多.有了這樣優(yōu)質(zhì)的均勻隨機(jī)序列,用擬蒙特卡羅方法求解得到的將是確定性的誤差,從而避免了蒙特卡羅方法得到概率誤差的缺陷.
下面給出兩個特殊數(shù)值算例驗證擬蒙特卡羅方法在計算多重積分中的有效性和優(yōu)越性.
例1 用Halton序列計算s重積分
取N=2 000,在計算機(jī)上分別用Halton序列和Rand函數(shù)計算的結(jié)果見表1.
的估計值.
表1 分別用Halton序列和rand函數(shù)計算的積分結(jié)果
表2 分別用Halton序列和rand函數(shù)計算的結(jié)果
通過計算可以看出,擬蒙特卡羅方法是有效的,積分重數(shù)對積分結(jié)果的誤差無明顯影響,精度較高,計算機(jī)程序?qū)崿F(xiàn)也是很簡單的.本文只是用了原始的Halton序列,由于Halton序列基數(shù)越大,相關(guān)性也越大,均勻性也會越差,這將會影響結(jié)果的精度,因此,如何選擇更優(yōu)質(zhì)的低差異序列有待進(jìn)一步學(xué)習(xí)和研究.
[1] 雷桂圓. 關(guān)于蒙特卡羅及擬蒙特卡羅方法的若干研究[D]. 杭州: 浙江大學(xué)理學(xué)院, 2003: 3-5.
[2] 杜紹洪. 高維積分的新型求積公式[D]. 成都: 四川大學(xué)數(shù)學(xué)學(xué)院, 2004: 11-15.
[3] Morokoff W J, Caflisch R E. Quasi-random Sequences and Their Discrepancies [J]. SIAM J SCI Comput, 1994, 15(6):1251-1279.
[4] Niederreiter H. Random Number Generation and Quasi Monte Carlo Methods [M]. Philadelphia: SIAM, 1992: 10-21.
[5] Morohosi H, Fushimi M. A Practical Approach to the Error Estimation of Quasi-Monte Carlo Integrations [EB/OL]. [2012-01-08]. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.37.9396&rep=rep1&type= pdf.
[6] Krykova I. Evaluting of Path-dependent Securities with Low Discrepancy Methods [D]. Massachusetts: Worcester Polytechnic Institute, 2003: 7-30.
[7] Brandimarte Paolo. Numerical Method in Finance: A MATLAB-based Introduction [M]. New York: John Wiley and Sons, 2002: 265-280.
[8] 宮野. 計算多重積分的蒙特卡羅方法與數(shù)論網(wǎng)格法[J]. 大連理工大學(xué)學(xué)報, 2001, 41(1): 21-22.
[9] Wang X Q, Fang K T. The Effective Dimension and Quasi-Monte Carlo Integration [J]. J Complexity, 2003, 19(2): 101-124.
Quasi Monte Carlo Method for Calculating Multiple Integrals
ZHANG Tao, ZHOU Zhongli, AN Suzhen, ZHANG Jun
(Key Laboratory of Mathematical Geology of Sichuan Province, Chengdu University of Technology, Chengdu, China 610059)
Basic principles of Quasi Monte-Carlo method for calculating multiple integrals were introduced. The uniformities of Halton sequences and of sequences produced by the Rand function were compared then. Later, steps, examples, and computer program written in MATLAB language of Quasi Monte- Carlo method for calculating multiple integrals were given. The effectiveness, accuracy and superiority of the proposed Quasi Monte- Carlo method could be fully demonstrated by examples.
Quasi Monte-Carlo Method; Halton Sequence; Multiple Integral
O242
A
1674-3563(2012)05-0033-06
10.3875/j.issn.1674-3563.2012.05.006 本文的PDF文件可以從xuebao.wzu.edu.cn獲得
(編輯:王一芳)
2011-12-16
張濤(1986- ),男,河南信陽人,碩士研究生,研究方向:數(shù)學(xué)計算方法