常系數(shù)三維橢圓邊值問題的蒙特卡羅算法

柯李瑞

(溫州大學(xué)數(shù)理與電子信息工程學(xué)院，浙江溫州 325035)

在巖土工程、機(jī)械設(shè)計(jì)、彈性力學(xué)、材料力學(xué)等工程和物理問題中，常會(huì)遇到只需求解某幾個(gè)特殊點(diǎn)處的函數(shù)值的集中載荷問題.這種情況下，繼續(xù)使用有限元方法，會(huì)計(jì)算很多不必要的值，造成浪費(fèi).因此，文[1]提出了一種新的求解橢圓型方程的方法——有限元的概率算法，文[2]確立了有限元概率算法的基本理論，文[3]在此基礎(chǔ)上針對(duì)一類特殊的二維橢圓邊值問題提出了一種高效蒙特卡羅算法，使工作量大大減少.本文給出了一種求解常系數(shù)三維橢圓邊值問題的蒙特卡羅算法，并通過算例說明了該方法的可行性.

1 一類橢圓邊值問題的蒙特卡羅算法

問題的提出：

本文思路.首先，將方程(1)轉(zhuǎn)化為如下方程：

然后，運(yùn)用有限元的蒙特卡羅算法求解上述方程即可.

所以有：

即BBT=A，其中

由A的對(duì)稱正定性，可以求解得到.從而，方程(1)可以轉(zhuǎn)化為如下形式：

求解方程(1)在任意一點(diǎn)數(shù)值解的問題轉(zhuǎn)化為求解方程(2)在任意一點(diǎn)數(shù)值解的問題.

方程(2)是一類經(jīng)典的方程，求解這類方程的關(guān)鍵在于Ω1和f(y)，還有c.當(dāng)Ω1，f(y)以及c中有一個(gè)改變時(shí)，方程的解也隨之改變.

當(dāng)Ω1為一般區(qū)域時(shí)，求解方程(2)的數(shù)值解比較困難.所以，本文先考慮簡單情形：Ω1為以y0為球心，以R為半徑的球.

設(shè)存在游粒β，當(dāng)它位于點(diǎn)y時(shí)，質(zhì)量為一單元的游粒擁有能量v(y).假設(shè)在初始狀態(tài)下，處有一單位質(zhì)量的游粒β，每次它都直接游動(dòng)到Γ1，且游動(dòng)到Γ1上的每一點(diǎn)的概率相同.記游到Γ1上的點(diǎn)為，則為隨機(jī)變量[4]，從而也為隨機(jī)變量，關(guān)于有以下結(jié)論.

定理1[1]隨機(jī)變量的期望與方差均存在，且

由以上分析可知，游粒游動(dòng)到Γ1上的每一點(diǎn)的概率相同，所以可以構(gòu)造如下求解的概率模型[5]：假設(shè)存在常數(shù)h＞0.對(duì)邊界Γ1進(jìn)行三角形網(wǎng)格劃分，以Ω1的內(nèi)接正八面體為基礎(chǔ)，將邊界Γ1分為8個(gè)等球面三角形，然后應(yīng)用“經(jīng)緯度平分法”進(jìn)一步細(xì)分球面三角形，直到Γ1上球面三角形的每兩個(gè)頂點(diǎn)之間的距離不超過h為止.這樣就得到邊界Γ1的近似均勻劃分.將網(wǎng)格頂點(diǎn)記為邊界Γ1上的結(jié)點(diǎn)，且依次記為，令每次游動(dòng)到的概率為：，這樣構(gòu)造的概率滿足以下條件：

上面給出了當(dāng)Ω1為以y0為球心的球時(shí)求解的方法.當(dāng)Ω1為一般區(qū)域[6]時(shí)，由于計(jì)算困難，所以不使用游粒從y0直接到達(dá)Γ1上的方法來求，但可采用下述方法求：

取定一個(gè)足夠小的數(shù)ε，在Ω1內(nèi)以y0為球心，以y0到邊界Γ1的最小距離為半徑作球，并將其邊界記為，半徑記為R0.游粒從y0出發(fā)，直接隨機(jī)游到，并且到達(dá)上的每一點(diǎn)的概率相同.記y0到達(dá)上的點(diǎn)為，有.由于未知，所以游粒需要繼續(xù)游動(dòng).設(shè)y0到達(dá)上的點(diǎn)為y1，在Ω1內(nèi)以y1為球心，以y1到邊界Γ1的最小距離為半徑作球，并將其邊界記為，半徑記為R1.游粒直接從y1出發(fā)，隨機(jī)游到，并且到達(dá)上的每一點(diǎn)的概率相同.記y1到達(dá)上的點(diǎn)為，有.讓游粒這樣一直游動(dòng)，若該游粒在游動(dòng)ni次后到達(dá)Γ1，記該次實(shí)驗(yàn)終點(diǎn)為.如果在游動(dòng)ni次后仍未到達(dá)邊界Γ1，但到邊界Γ1的距離小于或等于ε，不妨把它設(shè)為，讓游粒直接游到邊界Γ1上離它最近的點(diǎn)，記這次實(shí)驗(yàn)終點(diǎn)為.

設(shè)此次游動(dòng)共進(jìn)行了N次，由上述分析可得：

設(shè)實(shí)驗(yàn)次數(shù)為N，由以上分析，可得本文求解方程(2)的任意一點(diǎn)y0的數(shù)值解的蒙特卡羅算法的步驟如下：

1)取定N和ε；

2)取S=0；

3)取C=0；

4)取y=y0；

5)確定點(diǎn)y是否在Γ1上，如果點(diǎn)y在Γ1上，則轉(zhuǎn)12)，如果點(diǎn)y不在Γ1上，則轉(zhuǎn)6)；

6)求出y到Γ1的最小距離

8)以點(diǎn)y為球心，以Ri為半徑作球，其中；

9)令y隨機(jī)游動(dòng)到Γi上任一點(diǎn)(游動(dòng)到Γi上任一點(diǎn)的概率相同)，記游動(dòng)終點(diǎn)為yi，轉(zhuǎn)5)；

10)找到y(tǒng)距離Γ1最近的點(diǎn)，求出兩點(diǎn)間距離，并讓游粒直接到達(dá)該點(diǎn)；

13)判斷C＜N是否成立，若成立，轉(zhuǎn)(4)；

以上討論了三維情況下求方程(2)在任意一點(diǎn)y0處的數(shù)值解的方法.由于本文要求的是方程(1)在指定點(diǎn)的數(shù)值解，所以首先應(yīng)該將所求點(diǎn)x0轉(zhuǎn)化為y0，然后用上述方法，通過求解y0處的數(shù)值解來得到x0處的數(shù)值解.

2 數(shù)值算例

考慮方程

表1 點(diǎn)游動(dòng)結(jié)果與方程精確解的比較Table 1 Comparison with the Result of the Equation at point and the Exact Solution of the Equation

通過上述算例，說明本文提出的對(duì)于常系數(shù)三維橢圓邊值問題的蒙特卡羅算法是可行的.

[1] 朱起定.有限元概率算法[J].湘潭大學(xué)自然科學(xué)學(xué)報(bào)，1989，11(3)：1-5.

[2] 朱起定.有限元概率算法的基本理論[J].湘潭大學(xué)自然科學(xué)學(xué)報(bào)，2001，23(2)：121-133.

[3] 何文明，崔俊芝.一類特殊的橢圓型問題的高效蒙特卡羅算法[J].系統(tǒng)科學(xué)與數(shù)學(xué)，2004(2)：210-217.

[4] 朱起定.橢圓邊值問題的概率算法[J].系統(tǒng)科學(xué)與數(shù)學(xué)，2002(2)：168-179.

[5] 孫文彬，趙學(xué)勝，高彥麗，等.球面似均勻格網(wǎng)的剖分方法及特征分析[J].地理與地理信息科學(xué)，2009(1)：53-60.

[6] 朱起定.調(diào)和方程第一邊值問題高效概率算法[J].計(jì)算數(shù)學(xué)，2000(1)：121-128.

Monte-Carlo Method for Three-dimensional Elliptic Boundary Value Problem with Constant Coefficient