蒙特卡罗方法在热传导方程中的应用
第22卷 第11期Vol.22 No.11重庆工学院学报(自然科学)
JournalofChongqingInstituteofTechnology(NaturalScience)
2008年11月Nov.2008
蒙特卡罗方法在热传导方程中的应用
史慧会,曲小钢
(西安建筑科技大学理学院,西安 710055)
Ξ
摘要:介绍了蒙特卡罗方法的基本原理,建立了热传导方程边值问题的概率模型,并对采用蒙特卡罗方法所产生的误差进行分析,最后给出了在计算机上实现的方法.关 键 词:热传导方程;蒙特卡罗方法;随机函数;数学期望中图分类号:O21 文献标识码:A
文章编号:1671-0924(2008)11-0101-03
ApplicationofCiao2gang
Science,Xi’anUniversityofArchitectureandTechnology,Xi’an710055,China)
Abstract:ThispaperintroducesthebasicprincipleofMonte2Carlomethod,establishestheprobabilitymodelofboundaryvaluefortheheatinductionequation,analyzestheerrorcausedbyMonte2Carlomethod,andfinally,givesthemethodto
implementitincomputer.
Keywords:heatconductionequation;Monte2Carlomethod;stochasticfunction;mathematicalexpectation
由于温度不均匀,热量从温度高的地方向温度低的地方转移,这种现象称为热传导.一维与二维的热传导方程的形式为
2一维:=a2
9t9x
22
二维:=a(2+2)
9t9x9y
对于热传导边值问题的求解,可以采用有限差分法、有限元法或区域分解算法等数值方法,运算量都比较大,而且运算复杂.本文中介绍一种比较简便的方法———蒙特卡罗方法来解热传导方程的边值问题.
Ξ 收稿日期:2008-08-25
基金项目:陕西省教育厅专项研究基金资助项目(05JK239);西安建筑科技大学基础研究基金资助项目(03BR02).
),女,陕西西安人,硕士研究生,主要从事数值计算方面的研究.作者简介:史慧会(1982—
1 蒙特卡罗方法
蒙特卡罗方法也称为统计实验法或统计模拟法,是20世纪40年代中期随着科学技术的发展和电子计算机的发明而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法.蒙特卡罗方法的基本原理及思想如下:当所要求解的问题是某种事件出现的概率或者是某个随机变量的期望值时,可以通过某种“试验”的方法,得到这种事件出现的频率,或者这个随机变数的平均值,
102 重庆工学院学报
令r=a
:
2,整理得h
uk,j-uk,j-1=ruk+1,j-2ruk,j+ruk-1,j(5)
并用它们作为问题的近似解[1-2].
蒙特卡罗方法解题有3个主要步骤
[3]
1)构造或描述概率过程.对于本身就具有随
机性质的问题,主要是正确描述和模拟这个概率过程,对于不是随机性质的确定性问题,比如计算定积分,就必须事先构造一个人为的概率过程,它的某些参量正好是所要求问题的解.
2)从已知概率分布进行抽样.构造了概率模
当1-r≤0时,上式可以改写成
uk,j=
uk,j-1+u+1+2r1+2rk+1,j
1+2r
uk-1,j=α1uk-1,j+
(6)
α2uk+1,j+α3uk-1,j
α其中α0,i=1,2,3.∑i≥i=1.
i=1
型以后,由于各种概率模型都可以看作是由各种各样的概率分布构成的,因此产生已知概率分布的随机变量(或随机向量)就成为实现蒙特卡罗方法模拟实验的基本手段,这也是蒙特卡罗方法被称为随机抽样法的原因.最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布(或称矩形分布).
3)建立各种估计量.一般说来,构造了概率模
3
设想一种随机游动过程,位于点p(k,j)的粒
子分别以概率α1,α2,α3向邻点游至j=0或边界点Qk1,此时赋值以f(Qk1).如果想要计算
u(xk,tj)的值,只需从点p(k,j)开始作这样的随
机游动,取得一个u(xk,tj)f(Qk1).将此过M,得到k)k2…,f(QkM),于(,tj)M∧
型并能从中抽样后,即实现模拟实验后,一个随机变量,,偏估计..
Mi=1
∑f(Qki).
M
3 误差分析
蒙特卡罗方法是用随机变量的样本研究随机
2 热传导方程边值问题的数值解
设求解区域Ω中的问题为:在区域
R:{0≤x≤l,0≤t≤T}内求函数,满足方程
变量的分布,只不过样本不是由实验观测取得,而是用数学方法模拟产生而已.
在区域内任意一点满足方程(1)的近似解为
s
2=a,0
及初始条件
ut=0=φ(x),0≤x≤l
w(p)=E(ξ(p))=
(1)
i=1
∑f(Q)u(p,Q)
j
j
(7)
其中
u1(t),x=0
f(Q)=u(x,t)=
u2(t),x=l
(2)
边界条件
u(0,t)=μ1(t)u(l,t)=μ2(t)
φ(x),t=0
0≤t≤T
(3)
f(Qj)为f(Q)在边界网格点Qj上的值,u(p,Qj)
表示质点从p出发,游动到边界网格点Qj的概率.当每个质点游动到边界网格点时,一次模拟实验就结束了.
定义随机变量
0,当第i个质点不到达边界网格点Qj.
则M个质点随机游动到边界网格点Qj(j=1,2,
1,当第i
个质点到达边界网格点Qj,
并且φ(0)=μ1(0),φ(l)=μ2(0),其中μ1(t),μ2(t)是给定的函数.
假设h=Δx,τ=Δt分别是自变量x,t的改变量,令xk=kh,tj=τj,k=0,1,2,…;j=0,1,2,…,两组平行线构成的长方形网格覆盖了整个x-t平面,h称为沿x方向的步长,τ称为沿t方向的步长.用隐式三点差分格式代替上面的微分方程,得
h
2
ξij=
…,s)的质点数
M
s
τ
=a(4)mj=
i=1
ξ,∑m∑
ij
j=1
j
=M
史慧会,等:蒙特卡罗方法在热传导方程中的应用
E(ξij)=0・p(ξij=0)+1・p(ξij=1)=u(p,Q)
103
^
以此来判断u(k,j)取值,共取M次,则u(k,j)=Mi=1
注意对固定的j,ξij(j=1,2,…,M)是独立同分布的,故有
E(ξij)=E(
∑f(xki).
M
M
M
j=1
ξ)∑
ij
=
M
M
j=1
∑E(ξ)
ij
5 结束语
由契比晓夫不等式得
j=1
ξ∑ij-u(p,Qj)
M
≥
(ξ)
2
εM
=
2
4εM
,
研究了MC方法在解热传导边值问题的应用,主要以一维热传导方程为例来研究.首先要建立适合该问题的概率统计模型,其次对该模型进行随机抽样,最后用其算术平均值作为所求方程的近似解.另外MC方法求解程序简单且误差与空间的维数无关,这是其它方法所没有的优点.
j=1,2,…,s.
故误差为
w(p)-w(p)
^
^
≤∑≥
j=1
s
sf(Q)
2
4εM
22
其中M为实验次数,f(Qi)为f(Q)在边界网格点
Qi上的值,w(p)为满足方程(1)的解w(p)的估计
^
值.
在用蒙特卡罗方法计算时,其值具有概率的性质,从所得的误差表达式可以看出,的误差,.
参考文献:
[1].系统科学与
)].常微分方程的数值解
法[M].北京:科学出版社,1979.
4 利用Matlab语言编写了上述热传导方程边值问题数值解的程序[4-7].程序中利用随机函数rand产生(0,1)区间上的均匀分布,然后作如下处
[3] 靖新.概率论与数理统计[M].大连:大连理工大学出
版社,2005.
[4] 萧骁.QuickBasic标准编程指南[M].北京:希望电脑
公司,1991.
[5] 宋兆基,徐流美.
Matlab在科学计算中的应用[M].北
理:
ifξ0,
3
,,
u(k,j)=u(k,j-1)u(k,j)=u(k+1,j)u(k,j)=u(k-1,j)
京:清华大学出版社,2005.
[6] 王家生,刘嘉 .随机过程基础[M].天津:天津大学
出版社,2004.
[7] 李维.基于蒙特卡罗统计试验法的电算模型[J].重
ifξ,
33ifξ,1,3
庆工学院学报:自然科学,2007,21(4):24-27.
(责任编辑 刘 舸)