計算多重積分的均勻隨機數蒙特卡羅法的實現
2008-12-31 00:00:00劉輝玲葉鋒
電腦知識與技術 2008年35期
摘要：蒙特卡羅法是用一系列隨機數來近似解決問題的一種方法。采用均勻隨機數蒙特卡羅法計算多重積分是一種簡單而有效的方法，其程序結構簡單，易于編制和調試。本文給出了采用均勻隨機數蒙特卡羅法計算多重積分的步驟、算法流程，并給出了實例的具體實現過程。
關鍵詞：蒙特卡羅法；多重積分；均勻隨機數
中圖分類號：TP393文獻標識碼：A文章編號：1009-3044(2008)35-2289-03
Realization of Monte Carlo Method by Using Uniform Random Sequence for Calculating Multiple Integrals
LIU Hui-ling1，YE Feng2
(1.Department of Electronic and Electrical Engineering， Wuhan Institute of Shipbuilding Technology， Wuhan 430050，China; 2. School of Mathematics and Computer Science， Jianghan University， Wuhan 430056， China)
Abstract: Monte Carlo method is a approximate way to solve problems by using a series of random sequence. Monte Carlo method by adopting uniform random number is a simple and effective way to calculate multiple integrals， its structure is simple and easy to program and debug. This paper describes the realization steps and algorithm using Monte Carlo method of by using uniform random sequence to calculate multiple integrals， and gives the realization of concrete examples.
Key words: monte carlo method; multiple integrals; uniform random sequence
蒙特卡羅法[1-2]，也稱為統計試驗方法，是用一系列隨機數來近似解決問題的一種方法，是通過尋找一個概率統計的相似體并用實驗取樣過程來獲得該相似體的近似解的處理數學問題的一種手段。例如，f(x)在a 在計算多重積分時，如果被積函數的原函數難以求出或原函數不是初等函數獲取其解析解就變得較為困難。此時，可采用蒙特卡羅法進行計算，獲取其近似解[3]。均勻隨機數蒙特卡羅法多重積分的計算簡單，結構清晰，是一種有效的方法。本文給出了其計算步驟和流程，并用C語言予以實現。 1 多重積分的蒙特卡羅方法[4-6] 任何一個積分都可看作某個隨機變量的數學期望。因此，在利用蒙特卡洛法計算多重積分時，采用了這個隨機變量的算術平均值來作為其近似值。 以一般的S重積分為例 式(1)中：P＝P（x1，x2， …，xs）表示S維空間的點；Vs表示積分區域。 在使用蒙特卡洛方法解算問題時，首先構造或描述概率過程。由于計算定積分本身不是隨機性質的確定性問題，在使用蒙特卡洛法求解時，故要將其由不具有隨機性質的問題，轉化為隨機性質的問題。因此在求解時必須構造一個概率過程（分布密度函數），使其某些參量正好是所要求問題的解。 (1)在所求積分區域上構造一個分布密度函數，取Vs上任一概率密度函數f(P)，它滿足條件f(P)≠0。當P∈VS，G(P)≠0時，令 則式(1)可改寫為 即θ是隨機變量g(P)的數學期望。 其次，建立各種估計量，使其期望值是所要求解問題的解。 (2)用算術平均值來近似g(P)的數學期望。現從f(P)中抽取隨機變量P的N個樣本：{Pi}Ni=1，則算術平均值為 就是積分值θ的一個近似估計。 2 采用均勻隨機數的多重積分蒙特卡羅法 蒙特卡羅法計算多重積分時，采用均勻隨機數法得到隨機變量是常用方法之一。此時，選取f(P)的方法是取VS上的均勻分布，即 這里，VS表示積分區域的體積，此時g(P)=VSG(P)。 根據(1)，(3)，(4)和(5)式，采用均勻分布隨機數計算多重積分時，有 設f(x，y)為區域D上的有界函數，區域D：a≤x≤b，c≤y≤d，其面積A=(b－a)(d－c)，不妨設(xi ， yi )，i=1，…，N為落在區域D上的均勻分布隨機數列，根據(6)式，則用均勻隨機數計算二重積分 將(7)式推廣到多重積分，有 其中s為積分變量的個數，N的取值要充分大。 采用(8)式進行積分計算時，首先需要獲取隨機數列，（t1，t2， …，tk）代入被積函數，然后乘上對應積分變量上下限之差的累積，重復N次，并累加求和，當N充分大時，所得結果即為積分的近似結果。使用(8)式編寫程序計算積分時，其算法流程如圖1所示。 圖1所示流程可分解為以下8個步驟： ① 設置變量Sum，并令其初值為0，設置變量N并對其賦初值，N的取值要求充分大； ② 取0到1之間均勻分布的隨機數列，（t1，t2， …，ts），其中s為積分變量的個數； ③ 計算積分變量xj的取值uj，uj取值介于aj和bj之間，由uj＝aj+(bj-aj)*tj計算可實現，j=0，1，2， …，k，aj，bj為積分變量xj的積分上下限，重復s次，求得隨機數列（u1，u2，…，us）； ④ 將第③步求得隨機數列依次代入計算被積函數f(x0，x1，…，xs)，求得被積函數的值V； ⑤ 計算積分變量上下限差值的累積， ⑥ 計算第④步所求得的V值與第⑤步所求得的Mul的乘積，與Sum相加，并賦值給Sum； ⑦ 重復第②—⑥步N次； ⑧ Sum除以N，所得結果就是積分 3 算例及其實現 用蒙特卡羅求解 對于本例來說，其內層積分上限為一個表達式x，在采用蒙特卡羅計算積分時，其取值依賴于外層積分變量的取值，隨外層積分變量取值的變化而變化。故在實現時需對(8)式變化為(9)式。 C語言程序如下： double mtml(int s，int n，double a[]，double b[]，double (*f)()) //s為積分重數，n為抽取隨機變量樣本的次數 //a，b為數組，分別存儲各層積分的上下限 //f為指向被積分函數的指針 { int m，i; double r，Sum，*x; x=(double *)malloc(s*sizeof(double)); //x為存儲隨機變量樣本各個分量的數組 Sum=0.0; for (m=0; m { x[0]=a[0]+(b[0]-a[0])*Rnd(); //Rnd()為在0，1之間產生均勻分布隨機數的函數 b[1]=x[0]; x[1]=a[1]+(b[1]-a[1])*Rnd(); Sum＝Sum+(*f)(s，x，a，b); } Sum = Sum /d; free(x); return(Sum); } double mtmlf(int s，double x[]，double a[]，double b[]) //s為積分重數，x為存儲隨機變量樣本各個分量的數組 //a，b為數組，分別存儲各層積分的上下限 { int i; double f=1.0，Value=1.0; for (i=0; i<=s-1; i++) Value=Value*x[i]*(b[i]-a[i]); f=f+Value; return(f); } 本例的解析計算值為1.25。在抽取樣本變量的個數為10000時，其近似解為1.119648，誤差為0.005352；在抽取樣本變量的個數為100000時，其近似解為1.125504，誤差為0.000504；在抽取樣本變量的個數為1000000時，其近似解為1.124758，誤差為；由此可見，采用蒙特卡洛法計算多重積分，當n的取值充分大時0.000242，其近似解已經相當接近解析解。 4 結論 采用均勻隨機數的蒙特卡洛法計算多重積分，是一種簡單而有效的手段，其程序結構簡單，便于編制和調試；對于多維問題有普遍適用性，其收斂速度與問題維數無關，要達到同一精度，用蒙特卡羅方法選取的點數與維數無關，計算時間僅與維數成比例，而其他的數值方法，達到同樣的誤差，點數與維數的冪次方成比例。同時，均勻隨機數的蒙特卡洛法計算多重積分存在著收斂速度慢且具有概率性質的特點，增加抽取樣本次數n，計算量增加，其精度提高很慢，但可以通過選擇適當的隨機數，如有利隨機數，可使其精度提高。 參考文獻： [1] 宮野.計算物理[M].大連:大連理工大學出版社，1987. [2] 徐鐘濟.蒙特卡羅方法[M].上海:上海科學技術出版社，1985. [3] 何鳳霞，張翠蓮.蒙特卡羅方法的應用及算例[J].華北電力大學學報，2005，32(3):110-112. [4] 尹增謙，管景峰，等.蒙特卡羅方法及應用[J].物理與工程，2002，12(3):46-50. [5] 宮野.計算多重積分的蒙特卡羅方法與數論網格法[J].大連理工大學學報，2001，41(1):20-23. [6] 裴鹿成，張孝澤.蒙特卡羅方法及其在粒子輸運問題中的應用[M].北京:北京科學出版社，1980.
