应用蒙特卡罗方法求解一类随机微分方程

第36卷　第4期2003年7月 天　津　大　学　学　报Journal of Tianjin U niversity Vol.36　No.4J ul.2003
应用蒙特卡罗方法求解一类随机微分方程Ξ
张　华1,练继建1,刘嘉　2
(1.天津大学建筑工程学院,天津300072;2.天津大学理学院,天津300072)
摘　要:建立一类随机微分方程初值的概率模型,应用蒙特卡罗(Monte 2Carlo )法对其抽样产生一组伪随机数,应用四阶龙格2库塔(Runge 2Kutta )法求解随机微分方程.给出了一个实例,求得其解析解和数值解,在计算次数大于50
和小于100的条件下,数值解的最大相对误差为3.60
0.关键词:蒙特卡罗法;随机微分方程;数值解
中图分类号:O211.63 文献标识码:A 文章编号:049322137(2003)0420430204
Monte 2C arlo Method for C alculating a Class of
Stochastic Differential Equation
ZHAN G Hua 1,L IAN Ji 2jian 1,Liu Jia 2kun 2
(1.School of Civil Engineering ,Tianjin University ,Tianjin 300072,China ;
2.School of Sciences ,Tianjin University ,Tianjin 300072,China )
Abstract :A probability model was first established for the starting value of a class of stochastic differential equa 2tion.Then a series of pseudo 2random number was generated using Monte 2Carlo method.Finally Runge 2Kutta method was used to solve the problem.To verify this method ,an example was presented.After the analytic solution and the numerical solution were found ,the maximum relative error of the numerical solution to the
analytic solution was 3.60
0under the condition of more than 50and less than 100times calculations.K eyw ords :Monte 2Carlo method ;stochastic differential equation ;numerical solution
在自然科学和社会科学中,建立了许多用确定性微分方程描述的数学模型.人们认为只要输入所有的初始条件和边界条件,就能确切地预测未来的系统状态.但初始条件和边界条件中的误差是不可避免的,况且有些问题本身就是随机的,使得确定性数学模型的实际价值大大减低.这时,须用随机微分方程代替确定性微分方程,建立用随机微分方程描述的随机数学模型[1～3].
仅含随机初始条件或随机边界条件的随机微分方程称为一类随机微分方程,文中就应用蒙特卡罗方法求解一类随机微分方程的问题作一探讨.
1　一类随机微分方程的数值解
1.1　基本方程
研究随机微分方程组[4～6]
d X i (t )
d (t )
=f i (t ,X 1(t ),X 2(t ),…,X n (t ))
X i (t 0)=X i 0
T ∈[t 0,t f ]
i =1,2,…,n
(1)
式中:X i 0、X i (t )都是H 中的二阶矩随机变量,H 是二阶矩随机变量空间.所有的运算都是在均方意义下的,采用向量形式,式(1)成为
d X (t )
d t =f (t ,X (t ))X (t 0)=X 0
(2)式中:X (t )=(X 1(t ),X 2(t ),…,X n (t ));
f =(f 1,f 2,…,f n );
X 0=(X 10,X 20,…,X n 0),X 0为随机变量.
X (t )=(X 1(t ),X 2(t ),…,X n (t ))是n 个二阶矩随机变量组成的n 维向量,它们组成一个线性空间,
Ξ收稿日期:2002206224;修回日期:2002210222.
作者简介:张　华(1962—
),男,博士研究生,副教授.
这个空间的范数定义为
‖X (t )‖n =max 1≤j ≤n
‖X j (t )‖
这个空间就是H n .由于H 是完备的,H n 也是完备的.
因此H n 是一个Banach 空间.
在式(2)中,若f :T ×H n →H n 满足Lipschitz 条件
‖f (t ,X (t ))-f (t ,Y (t ))‖n ≤
k (t )‖X (t )- Y (t )‖n 其中k (t )满足
∫t
f
t 0
k (t )d t <
∞
则对任意X 0∈H n ,式(2)等价于 X (t )=X 0+
∫
t
t 0f (t ,X (t ))d t (3)
在一般情况下,由f (t ,X (t ))的任意性,式(3)无法得到解析解,现寻求其数值解.
1.2　随机微分方程数值解的基本步骤 蒙特卡罗(Monte 2Carlo )方法是以概率统计理论
为基础,以随机抽样为其主要手段,通过统计试验来达
到计算某些参量的目的. 随机微分方程数值解的基本步骤如下: 1)建立随机微分方程初值X 0的概率模型F (X 0;θ),其中F 为函数形式,θ为参数向量.
2)按照概率模型F (X 0;θ)的特点,进行随机抽样产生一组伪随机数[7],即X 0.
3)将这些初值X 0分别代入随机微分方程组,应用四阶龙格2库塔(Runge 2Kutta )法[8]
求得随机微分方程的一组解. 事实上,对于式(2),当t =t i 时,X (t i )为已知,令t i +1=t i +Δt ,应用四阶龙格2库塔法求得
X (t i +1)=X (t i )+Δt
6
(K 1+2K 2+2K 3+K 4) K 1=f (t i ,X (t i ))
K 2=f (t i +12Δt ,X (t i )+1
2K 1Δt )
K 3=f (t i +12Δt ,X (t i )+1
2K 2
Δt )
K 4=f (t i +Δt ,X (t i )+K 3Δt ) 4)重复以上步骤m 次,计算X (t i +1)的数学期望,它近似等于X (t i +1)的算数平均值.
2　算　例
例1　假设从远离地面的某处,以v 0的初速度连
续不断地垂直抛出小球,v 0是服从N (0,1)的正态分布,抛球流量为每秒n 个,n =1000.设开始抛球的时
刻为t =0,垂直向下为z 轴的正方向,抛球点位于z 轴的零点.令z i =i ,记Δz i =(z i -1,z i ],i =1,2,…,10.计算t =14s 时,在Δz i 区间(i =1,2,…,10)小球数的数学期望.2.1　例1的解析解
小球运动的随机微分方程为
d 2z (t )d t
2
=g ,d z (0)d t =v 0,z (0)=0(4)
令X 1(t )=z (t ),X 2(t )=d z (t )
d t
式改写为
d X 1(t )
d t =X 2(t )d X 2(t )
d t =g
X 1(0)=0X 2(0)=v 0
(5)
因为抛球流量为每秒n 个,所以1
n
秒抛一个小
球,t 秒内抛球nt 个.先求t 时刻第一个小球处于Δz i
区间的概率.将式(5)代入式(3)得
d z (t )
d t
=v 0+gt
z (t )=v 0t +1
2
gt 2
根据文献[4],z (t )的分布密度为 f (t ;z (t ))=
1
2πt
exp [-(z -
12
gt 2)22t 2
]
第一个小球在t 时,处于Δz i 区间的概率,即处于Δz i 区间小球数的数学期望为
∫
z
i
z i -1
f (t ;z (t ))d z =
∫
z
i
z i -1
f (
nt n ;z (nt
n
))d z 故t 时刻nt 个抛出的小球处于Δz i 区间小球数的数学
期望为
∫
z i z i -1
f (1n ;z (1n
))d z +∫
z i z i -1
f (2n ;z (2
n
))d z +
…+
∫
z
i
z i -1
f (nt n ;z (
nt
n ))d z =
6nt
j =1
∫
z
i
z i -1
f (
j
n
;z (
j
n
))d z ≈ n
∫t 0∫z
i
z i -1
f (t ;z (t ))d z d t
代入参数,求得其解析解,如表1所示.
・
134・　2003年7月 张　华等:应用蒙特卡罗方法求解一类随机微分方程
表1　例1的解析解与数值解比较
T ab.1　Comparison betw een numerical solution and analytic solution in exemple 1
m 数值解
5102030405060708090100最大相对误差/0 0m <50m ≥50
解析解
Δz 1380.6383.9377.9382.2384.3381.8380.1380.4380
.7380.4382.10.90.3380.
9Δz 2179.4178.1185.1181.9182.5182.4184.3183.8184.0185.8182.83.21.0183.9Δz 3140.4143.4142.9145.0140.2142.3143.0143.4142.4143.1142.42.00.9142.1Δz 4128.0121.0122.2119.9119.7120.3119.0119.9119.5119.0120.16.60.9120.1Δz 5102.2106.7102.6107.7106.9105.6105.0105.1104.9106.0105.53.61.0106.0Δz 697.898.496.892.396.997.899.498.495.794.997.43.83.695.9Δz 789.287.990.889.089.888.086.886.588.888.488.52.82.088.3Δz 880.283.983.182.183.083.882.182.182.483.181.52.41.982.2Δz 976.877.174.479.377.476.277.476.176.676.977.23.81.677.3Δz 1074.0
68.7
76.7
69.9
72.8
73.5
73.6
72.7
75.1
72.7
73.0
6.0
2.7
73.1
2.2　例1的数值解
由于t 秒内抛球nt 个,使得例1的问题比一般问题复杂.首先需要nt 次求解式(5),得到nt 组解z (t )
和d z (t )/d t ,然后根据拉格朗日质点跟踪法,跟踪每一个小球的位置,统计在区间的小球数,这样重复m
次可得在Δz i 区间小球数的数学期望.计算结果如表1和图1所示.可见,当计算次数小于50时,数值解的
最大相对误差为6.60
0;当计算次数大于50和小于100时,数值解的最大相对误差为3.60
0.(a )　Δz 1小球数的数学期望(c )　Δz 6～Δz 8小球数的数学期望(b )　Δz 2～Δz 5小球数的数学期望
(d )　Δz 9～Δz 10小球数的数学期望
图1　例1的解析解与数值解比较
Fig.1　Comparison betw een numerical solution and analytic solution in exemple 1
・234・天 津 大 学 学 报 第36卷　第4期　
3　结　语
应用蒙特卡罗法(Monte2Carlo)和四阶龙格2库塔(Runge2Kutta)法,对于一类随机微分方程建立了—个求解模式.利用Matlab软件,通过求解小球运动的随机微分方程,得到其解析解和数值解,在计算次数大于50和小于100的条件下,数值解的最大相对误差为3.60 0.可见一类随机微分方程的求解模式是可行的.
一类随机微分方程的求解模式不仅适宜于一维问题,而且特别适宜于多维问题,计算时间仅与维数成比例,问题的维数增加不会增加求解的难度.随着计算次数的增加,数值解将收敛于解析解.本求解方法的弱点是收敛速度慢,且存在初值随机抽样随机误差和四阶龙格2库塔计算方法的误差.
在高坝挑流消能雾化课题中,水滴随机喷溅问题是至关重要的,应用解析的方法很难解决这个问题.作者应用一类随机微分方程的求解模式,在这方面做了一些尝试,成功地解决了挑流水舌入水喷溅引起的降雨问题.
参考文献:
[1]　陈建兴,戎红玉,王惠民.用随机微分方程模型求解明渠
渐变流[J].河海大学学报,1994,22(5):43—47.
[2]　金　明.一维稳态河流水质的随机微分方程模型[J].水
利学报,1991(2):19—25.
[3]　潘宏雨,王国成.地下水运动随机微分方程浅论[J].勘
察科学技术,1998(4):7—10.
[4]　刘嘉　.应用随机过程[M].北京:科学出版社,2000.
[5]　Friedman A.随机微分方程及其应用[M].吴让泉译.北
京:科学出版社,1983.
[6]　S oong T T.Random Dif f erential Equations in Science and
Engineering[M].New Y ork:Academic Press,1973. [7]　官　野.计算物理[M].大连:大连理工大学出版社,
1987.
[8]　刘则毅.科学计算技术与Matlab[M].北京:科学出版
社,2001.
・
3
3
4
・
　2003年7月 张　华等:应用蒙特卡罗方法求解一类随机微分方程。