高效蒙特卡羅方法在偏微分方程初邊值問題中的應(yīng)用*

游 皎，李萬愛

(中山大學(xué)中法核工程與技術(shù)學(xué)院，廣東 珠海 519082)

?

高效蒙特卡羅方法在偏微分方程初邊值問題中的應(yīng)用*

游 皎，李萬愛

(中山大學(xué)中法核工程與技術(shù)學(xué)院，廣東 珠海 519082)

介紹了蒙特卡羅方法求解偏微分方程初邊值問題的基本原理，并針對該方法在多節(jié)點(diǎn)、多游動次數(shù)計(jì)算中速度慢的問題，提出了一種改進(jìn)方案。通過理論分析和算例驗(yàn)證，在相同計(jì)算條件下與傳統(tǒng)方法相比，該改進(jìn)方案不僅大大減少了計(jì)算時間，而且降低了誤差，這將使蒙特卡羅方法得到更廣泛的應(yīng)用。

初邊值問題；蒙特卡羅方法；隨機(jī)游動

蒙特卡羅方法是一種利用計(jì)算機(jī)進(jìn)行統(tǒng)計(jì)模擬的數(shù)值計(jì)算方法，已經(jīng)被應(yīng)用到包括物理、生物、化學(xué)和金融等眾多領(lǐng)域[1-2]。而在科學(xué)計(jì)算中常常需要求解偏微分方程初邊值問題，而該類問題傳統(tǒng)上均使用確定性方法如有限元方法。蒙特卡羅與有限元方法的區(qū)別在于：①后者需要使用方程的弱解形式進(jìn)行離散，而前者可以使用各種有效的游走規(guī)則；②后者求解大型線性方程組獲得求解點(diǎn)值，而前者利用概率統(tǒng)計(jì)理論建立求解值與已知值的關(guān)系；③后者在非確定性問題上的求解需要建立相關(guān)隨機(jī)模型，而前者自然地具有概率性求解的特點(diǎn)。故一些文獻(xiàn)也開始探索如何有效地利用蒙特卡羅方法來求解微分方程[3-8]，為了使結(jié)果更加精確，利用蒙特卡羅法時往往需要對求解域進(jìn)行細(xì)密的網(wǎng)格劃分。而節(jié)點(diǎn)數(shù)增大，游走次數(shù)增多后，計(jì)算量變得很大，需要較長計(jì)算時間，這嚴(yán)重限制了蒙特卡羅方法的實(shí)際應(yīng)用。文獻(xiàn)[4]提出了一種改進(jìn)方法，該方法已經(jīng)比已有方法(文獻(xiàn)[5-6])在計(jì)算效率上有較大改進(jìn)。本文繼續(xù)對文獻(xiàn)[4]方法加以改進(jìn)，方程邊界條件處理采用文獻(xiàn)[7]方法，并通過理論分析和算例結(jié)果說明，本文的改進(jìn)方法在相同的條件下，不僅可以大大減少計(jì)算時間，而且降低了誤差。這些改進(jìn)將使蒙特卡羅方法能應(yīng)用到更為占用計(jì)算資源的初邊值偏微分方程計(jì)算中。

1 基本原理及改進(jìn)方案

1.1 傳統(tǒng)方案

用蒙特卡羅方法求解微分方程，首先要構(gòu)造方程對應(yīng)的隨機(jī)概率模型。為具體起見，本文首先以二維拉普拉斯方程為例，對利用蒙特卡羅方法解決邊值問題進(jìn)行描述。設(shè)問題的控制方程為

(1)

其中u為要求解物理量，Ω為求解域。它滿足Dirichlet邊界條件

其中?Ω為Ω的邊界，up(x,y)為邊界值函數(shù)。首先對求解域Ω進(jìn)行網(wǎng)格劃分，然后相應(yīng)的離散格式作為蒙特卡羅方法的游動規(guī)則。這里以正規(guī)結(jié)構(gòu)化網(wǎng)格上的差分格式為例，采用五點(diǎn)差分格式離散上面的偏微分方程，可得

(2)

其中(i,j)為網(wǎng)格點(diǎn)編號，在均勻網(wǎng)格中有

可以將這種思路推廣到微分方程初邊值問題中，以擴(kuò)散方程為例，

(3)

其中方程的邊值和初值如下

對該方程采用時間后差、空間中心差的差分格式(BTCS)并整理可得

(4)

其中

p0=λ,p1=p0σx,p2=p0σy

1.2 改進(jìn)方案

上面使用蒙特卡羅方法求解偏微分方程時需要使用差分格式來對方程進(jìn)行離散。在工程計(jì)算里常常需要采用較密網(wǎng)格來得到可靠的結(jié)果，同時由于隨機(jī)游動過程具有隨機(jī)性，通常m要取得很大時才能得到較為準(zhǔn)確的結(jié)果。因此這種方法的計(jì)算量比較大，需要一些提高計(jì)算速度的方法來加速這個計(jì)算過程[1,4]。而文獻(xiàn)[4]的方法在加快計(jì)算速度及提高精度上都更有優(yōu)勢，它的改進(jìn)步驟為下面的隨機(jī)游動3：

在這個誤差減少的分析基礎(chǔ)上，本文提出一種先計(jì)算“大網(wǎng)格”的方法，通過改變網(wǎng)格點(diǎn)的計(jì)算順序以達(dá)到更好的優(yōu)化。具體來說，若求解域劃分為50×50的網(wǎng)格，可以首先計(jì)算i=10, 20, 30, 40和j=10, 20, 30, 40這4行4列的函數(shù)值，相當(dāng)于把求解域分成了25個小求解域，其邊界由外邊界和已知近似值的點(diǎn)構(gòu)成，這樣大大縮短了其它點(diǎn)的馬爾可夫鏈長度，因而縮短了計(jì)算時間。此外，在計(jì)算“大網(wǎng)格”時可以適當(dāng)增大游動次數(shù)以獲得更準(zhǔn)確的值。本文中采用的改進(jìn)版均是首先計(jì)算了間隔為5的網(wǎng)格點(diǎn)的值，再進(jìn)行其他點(diǎn)的計(jì)算。下面將這種方法稱為reorder方法。

1.3 誤差分析

隨機(jī)游動4其實(shí)是隨機(jī)游動3在初邊值問題的推廣。雖然文獻(xiàn)[4]已經(jīng)驗(yàn)證過隨機(jī)游動3的精度比隨機(jī)游動1更為提高，但是沒有進(jìn)行數(shù)學(xué)原理上的分析證明。這里將給出隨機(jī)游動3比1、隨機(jī)游動4比2能夠計(jì)算更準(zhǔn)確的分析。

這里有兩個誤差的問題。首先是游動規(guī)則帶來的誤差，本文的游動規(guī)則是通過有限差分方法給出的，其空間誤差為截?cái)嗾`差

時間誤差為O(Δt)。

第二個誤差為游動過程的概率統(tǒng)計(jì)誤差。首先需要研究對于傳統(tǒng)方法，游動次數(shù)與誤差之間的關(guān)系。利用文獻(xiàn)[6]的類似分析，以拉普拉斯方程為例，對于在求解域Ω內(nèi)任意一個網(wǎng)格點(diǎn)A, 由公式可以推得滿足方程的近似解為

(5)

式中s為邊界點(diǎn)個數(shù)，up(Qj)為邊界點(diǎn)Qj對應(yīng)的函數(shù)值up(x,y)；P(A,Qj)表示粒子從A點(diǎn)出發(fā)后游動到邊界Qj的概率。另外，定義隨機(jī)變量

ξi,j的期望E(ξi,j)為

方差V(ξi,j)為

用nj表示m個游動粒子中到達(dá)邊界點(diǎn)Qj的粒子數(shù)目，有

由契比雪夫不等式，結(jié)合式和式，對任意一給定的正數(shù)ε1>0，有

結(jié)合式和式，對任意小的正數(shù)ε2，有

由于ξi,j(i=1,2,…,m)具有相同的概率分布，并且相互獨(dú)立，因此有

所以

這個結(jié)論說明一點(diǎn)的游動次數(shù)m越大，則該點(diǎn)誤差越小。在離散公式(2)在均勻網(wǎng)格為例

一點(diǎn)的函數(shù)值取決于相鄰四點(diǎn)的值。如果以該點(diǎn)(i,j)為起點(diǎn)需要游動m次，則對于傳統(tǒng)方案，可以假定每個相鄰點(diǎn)的值都經(jīng)由m/4次游動后得到。若有一個鄰點(diǎn)的值是已知的，即它已經(jīng)經(jīng)歷過至少m次隨機(jī)游動，直接使用鄰點(diǎn)的已知函數(shù)值相當(dāng)于過該鄰點(diǎn)的游動次數(shù)至少為m2/4(乘法原理)。平均來說有兩個鄰點(diǎn)函數(shù)值已知則總次數(shù)為m/2+m2/2，這將大大增加隨機(jī)游動次數(shù)從而減少誤差。

2 數(shù)值例子

2.1 已知精確解的二維泊松方程

考慮一個溫度分布T(x,y)在域Ω=[0,1]×[0,1]滿足方程如下的泊松方程，

且邊界條件為T(x,y)=1+x。則此方程有解析解

分別采用傳統(tǒng)的(隨機(jī)游動1，對應(yīng)表 1方案0)和改進(jìn)后(隨機(jī)游動3+reorder方法，對應(yīng)表 1方案1)的方案對上述問題進(jìn)行求解，并改變參數(shù)條件，進(jìn)行橫向和縱向的對比。由于蒙特卡羅方法具有隨機(jī)性，表 1中每個方案都進(jìn)行了3次實(shí)驗(yàn)，結(jié)果取平均值。在表 1中，N為網(wǎng)格節(jié)點(diǎn)數(shù)，m為每個節(jié)點(diǎn)的游動次數(shù)，時間單位為秒，時間比為方案1的時間與方案0的時間之比，誤差采用L1范數(shù)，相對誤差為誤差與精確解之比?？梢钥闯?，在相同的條件下，改進(jìn)版不僅在時間上較傳統(tǒng)版有巨大優(yōu)勢，而且在降低誤差方面也表現(xiàn)得更好。

表1 泊松方程的蒙特卡羅求解方法的誤差對比

2.2 二維非定常熱傳導(dǎo)問題

考慮如下二維非定常熱傳導(dǎo)方程

其邊界條件及初始條件分別為

u(t,0,y)=u(t,1,y)=u(t,x,0)=u(t,x,1)=0,

u(x,y,0)=exp(-20(x-1/2)2-20(y-1/2)2)

由于非定常問題傳統(tǒng)方法(即隨機(jī)游動2)在二維中計(jì)算量過大，這里的傳統(tǒng)方法采用隨機(jī)游動2加上隨機(jī)游動4中的時間改進(jìn)方案。使用傳統(tǒng)方法(對應(yīng)表 2方案0)和改進(jìn)方案(隨機(jī)游動4+reorder方法，對應(yīng)表 2方案1)求解t=0.05時的分布。用二階Crank-Nicolson差分法取精細(xì)條件(N=1 000×1 000，dt=0.001)作為參考解，以此計(jì)算兩種方案的誤差。其計(jì)算結(jié)果及時間、誤差的比較如表 2所示?？梢钥闯龈倪M(jìn)方案不僅大程度地壓縮了運(yùn)行時間，而且在一定程度上減少了誤差。

為了進(jìn)一步分析計(jì)算結(jié)果，本文還作了溫度等值線進(jìn)行對比，如圖 1所示。改進(jìn)方案的等高線更加圓滑，說明該方案的隨機(jī)誤差更小，結(jié)果更穩(wěn)定。

表2 二維非定常熱傳導(dǎo)問題的蒙特卡羅法誤差比較(N=512,m=1 000)

圖1 常熱傳導(dǎo)方程蒙特卡羅法計(jì)算的溫度等值線對比Fig.1 Thermal isolines comparison using MC methods on unsteady heat conduction equation

3 總 結(jié)

本文研究了蒙特卡羅方法在微分方程初邊值問題中的應(yīng)用，針對傳統(tǒng)方法計(jì)算速度較慢的問題，提出了一種適合多節(jié)點(diǎn)、多游動次數(shù)的改進(jìn)方法并給出精度提高的分析證明過程。與較傳統(tǒng)方法相比，改進(jìn)方案無論在計(jì)算時間和結(jié)果精度上都有明顯優(yōu)勢，適合蒙特卡羅方法在更加復(fù)雜工程問題中的應(yīng)用。下一步工作將結(jié)合曲面處理技術(shù)[9],將該方法推廣到一般的幾何計(jì)算區(qū)域中。

[1] 劉軍. 科學(xué)計(jì)算中的蒙特卡羅策略[M]. 北京：高等教育出版社，2009.

[2] 洪志敏. 基于Monte-Carlo技術(shù)的積分(微分)方程數(shù)值求解方法研究[D]. 內(nèi)蒙古：內(nèi)蒙古工業(yè)大學(xué)博士學(xué)位論文.

[3] 左應(yīng)紅，王建國. 蒙特卡羅方法在解微分方程邊值問題中的應(yīng)用[J]. 強(qiáng)激光與粒子束，2012, 24(12): 3023-3026.

[4] 周泉，陳植武，詹杰民. 蒙特卡羅法的改進(jìn)：一種新型的快速算法[C]∥第二十一屆全國水動力學(xué)研討會暨第八屆全國水動力學(xué)學(xué)術(shù)會議暨兩岸船舶與海洋工程水動力學(xué)研討會文集, 2008.

[5] 林建國，周俊陶. 蒙特卡羅法的一種快速計(jì)算：在勢流問題中的應(yīng)用[J]. 水動力學(xué)研究與進(jìn)展A輯，2005, 20(3): 405-410.

[6] 史慧會，曲小鋼. 蒙特卡羅方法在熱傳導(dǎo)方程中的應(yīng)用[J]. 重慶工學(xué)院學(xué)報(bào)：自然科學(xué)版，2008, 22(11): 101-103.

[7]VAJARGAHBF,VAJARGAHKF.MonteCarlomethodforfindingthesolutionofDirichletpartialdifferentialequations[J].AppliedMathematicalSciences, 2007, 1(10): 453-462.

[8]LAUDUDP,BINDERK.AguildtoMonteCarlosimulationsinstatisticalphysics[M].CambridgeUniversityPress, 2009.

[9]JUNS,SUNGHWANY,NAMZC.AMonteCarlomethodforsolvingheatconductionproblemswithcomplicatedgeometry[J].NuclearEngineeringandTechnology, 2007, 29(3): 207-214.

The Efficient Monte-Carlo Method in Solving the Initial-Boundary Value Problem of Partial Differential Equations

YOUJiao,LIWanai

(Sino-French Institute of Nuclear Engineering & Technology, Sun Yat-sen university, Zhuhai 519082, China)

The primary principle of Monte Carlo method in solving the initial-boundary value problem for partial differential equations is introduced. And an improved strategy to solve the slow speed problem in computations with large grid points and random walks is proposed. Through the theoretical analysis and numerical validation, the improved strategy not only takes much less computational time, but also gives more accurate results comparing to the traditional methods under the same conditions.

initial-boundary value problem; Monte Carlo method; random walk

10.13471/j.cnki.acta.snus.2015.06.007

2015-01-17 基金項(xiàng)目：國家自然科學(xué)青年基金資助項(xiàng)目(11402313)；超級計(jì)算應(yīng)用領(lǐng)域資助項(xiàng)目(141gjc10)

游皎(1992年生)，男；研究方向：核工程與技術(shù)；通訊作者：李萬愛；E-mail:liwai@mail.sysu.edu.cn

O21

A

0529-6579(2015)06-0037-05