第1页 / 共6页
第2页 / 共6页
第3页 / 共6页
第4页 / 共6页
第5页 / 共6页
第6页 / 共6页
求解数值积分的蒙特卡罗方法.pdf
中国科技论文在线
http://www.paper.edu.cn
求解数值积分的蒙特卡罗方法
魏海燕
辽宁工程技术大学理学院，辽宁 阜新 (123000)
Email: sdweihy001@163.com
摘 要：积分的近似计算在许多实际工程的应用十分广泛,可积系统的一大类问题事实上就
是不规则区域上高维积分问题.这类问题绝大多数按数学的研究处理模式出发,都是难以解决
的,而蒙特卡罗方法对计算高维空间中维数不很高且积分区域是规则的积分、以及复杂区间
甚至不连通的区域上的积分都却是很方便的.本文主要介绍了蒙特卡罗法计算数值积分的两
种方法：随机投点法和平均值法,并且具体给出了这两种方法的计算原理以及方法步骤,然后
结合 Matlab 数学软件,用具体实例给出了定积分和多重积分的计算求解过程.
关键词：蒙特卡罗方法；数值积分；随机投点法；平均值法.
中图分类号：O241.4
1 引言
蒙特卡罗方法又称随机抽样法或统计试验方法,是一种重要的利用计算机模拟的近似计
算方法,它是 20 世纪 40 年代中期由于科学技术的发展,特别是核武器研制和电子计算机发明
而提出的,它在原子能技术中有广泛应用,也可以解决许多典型的数学问题.对于计算定积分、
重积分(特别是高维积分),蒙特卡罗方法是一种较为有效的方法,并且它特别适合于并行计算,
在计算机功能越来越强大的今天,这种数值计算方法虽然计算量较大,却很容易实现,所以计
算重积分成为蒙特卡罗方法最成功和典型的应用.
2 蒙特卡罗方法求解数值积分
2.1 随机投点法求积分的近似值
用 Monte-Carlo 方法求解一个问题时,基本工具就是一组数据,称为随机数.可以用一组随
机数来代替某些随机过程中的变数.计算机中的随机数由于是用确定法则获得的,因而称为伪
随机数.由“伪”逼“真”,理论上非常深入,同时也是评价方法优劣的权重.
考虑如下形式的 n 重积分：
= ∫
其中 D 表示 n 维单位立方体{0
在 D 上有 0
f p
(
≤
I
D
f p dp
(
)
,
p
=
p x x
(
,
2
1
,
x
)n
≤
ix
≤
1,
i
= .
n
1,2,
}
= ∫
I
f p dp
(
)
f P
(
)
D
dzdx
1
.
dx
n
) 1
≤ .很明显,此时积分
可以看作为求 1n + 维空间长方
体 :
V D
×
[0,
f p
(
)]
的体积.即
I
= ∫
∫ ∫
0
D
对于这种较为一般形式的多重积分问题,采用随机投点法时,首先产生 1n + 个随机数
.然后检验ξ是否落在V 中,那
)
ξ ξ ξ ξ η
=
(
,
,
,
,
1
2
n
ξ
i
(1,2,
么取
η 和 ,构造 1n + 维随机向量：
n
)
mI
= 作为近似值,即
N
mI
≅
N
[1].
2.2 平均值法求积分的近似值
- 1 -
中国科技论文在线
http://www.paper.edu.cn
Monte-Carlo 方法可构造数学期望来求积分值,特别是积分区域不规则的高维积分.而许
多问题的求解可从数字特征方面演绎为一种随机变量的数学期望.
对于 n 重积分
I
= ∫
D
f p dp
(
)
,其中
p
=
p x x
(
,
2
1
,
, D 表示 n 维空间中的积分区域.
x
)n
在 积 分 区 域 D 上 构 造 一 个 n 维 随 即 变 量 p, 其 密 度 函 数 为 (
p D
∈
令
p D
∉
) 0
=
时，
时，
(
ρ
(
ρ
0,
≠
p
p
.
)
)pρ , 满 足
F p
(
)
则
I
=
p
) /
(
ρ
p
f p
) , (
)
0
(
≠
ρ
⎧
= ⎨
p
) 0
0 , (
=
ρ
⎩
F p
p dp E F p
)
))
(
(
(
(
ρ
=
)
⋅
∫
D
,即此时 I 是随机变量 X 的函数 (
F X 的数学期望[2].
)
按 (
ip i
)pρ 抽取 D 中的随机样本 (
ˆ
I F
n
≈
n
, )
= ,由大数定律得：
1,2,
1
N
= ∑
N =
i
1
F p
(
i
)
(1)
这就是 Monte-Carlo 方法计算重积分的表达式.很容易在计算机上编程计算.选取 (
)pρ
的最简单方法是选取区域 D 上的均匀分布,设 D 的体积为 DV ，则
(
ρ
p
)
当
DV
1/
,
⎧
= ⎨
0,
当
⎩
p D
∈
p D
∉
这时
F p
(
)
所以
ˆ
I F
n
≈
p
f p V
, (
(
)
ρ
⎧
D
= ⎨
p
0 , (
ρ
⎩
V
1
D
n
n
F p
(
i
∑
=
=
)
n
i
1
=
0
)
≠
) 0
=
n
∑
i
1
=
f p
(
i
)
当积分是定积分
I
为：
b
= ∫
a
f x dx
( )
U a b
时,选取均匀分布的随机变量 ~ ( , )
ξ
,则(1)式转化
I
≈
b a
−
n
n
∑
i
1
=
f x
(
i
)
(2)
这就是求多重积分的蒙特卡罗方法,利用中心极限定理可以证明当 n → ∞ 时,抽样算术
nF 以概率 1 收敛到 I ,即
平均 ˆ
P
(lim
n
→∞
ˆ
F
n
=
I
) 1
=
蒙特卡罗方法的误差估计可以用中心极限定理推出.若方差
2
σ
=
∫
D
(
f p
(
)
−
I
)
(
ρ
p dp
)
< ∞
则在α置信水平下,不等式
ˆ
nF
ε− ≤ =
I
αλσ
n
成立.从而计算高维积分的误差估计为：
- 2 -
中国科技论文在线
http://www.paper.edu.cn
VI
D
−
n
n
∑
i
1
=
f p
(
i
)
≤
αλσ
n
增加 n 或减少方差 2σ 都可以提高计算精度.虽然误差阶仅为
,但是对于现代计
算机而言是可以接受的.为提高计算效率,仍要考虑减少σ[3],常用的方法有:方差缩减法、分层
抽样法等.
3 数值算例
o n−
(
1/2
)
3.1 蒙特卡罗方法求解定积分
假设函数 ( )
f x 在[ , ]a b 内有界连续且 ( ) 0
,为计算出它的
值,可以构造出一个概率模型：取一个边长分别为b a− 和 M 的矩形 D,使的曲边梯形包含在
矩形域之内,然后在矩形域内随机投点,假设随机点均匀的落在整个矩形域内,则落在曲边梯
形的随机点数 k 与投点总数 N 之比 /k N ,就近似的等于曲边梯形的面积与矩形面积之比，
f x ≥ ,对定积分
f x dx
( )
= ∫
I
a
b
从而得出定积分
b a M
(
−
)
.
I
≈
k
N
x dx
2
1
2
.
例：计算
∫
0
x= 用蒙特卡罗方法求解的步骤如下：
f x
( )
对于
P x y ；
(1) 产生矩阵域 D 内的 N 个均匀随机点 (
i
(2) 统计满足条件
=
(3) 取
k
N
y
≤
i
b
0,
=
k
b a M
(
f x
(
i
n
1,
10000
M
1,
≈
=
−
=
a
.
I
)
)
,
i
i
)
的落在曲边梯形内的随机点数 k;
10000
=
,则定积分的近似值为：
Matlab 程序如下：
xi=unifrnd(0,1,100000,1);
yi=rand(100000,1);
k=0;
y=xi.^2;
for i=1:100000
if yi(i)<=y(i)
k=k+1;
end
end
I=k/100000
用 Matlab 运行五次，所的计算结果如下：
0.3338 0.3339 0.3328 0.3329 0.3335 0.3320
由于我们得到的都是近似值,每次的计算结果都在准确值 1/3 附近作很小的摆动.
3.2 蒙特卡罗方法求解二重积分
f x y dxdy
( ,
对二重积分
I
)
= ∫∫
A
,设 ( ,
f x y 为区域 A 上的有界函数且 ( ,
f x y ≥ ,据其
) 0
)
- 3 -
中国科技论文在线
http://www.paper.edu.cn
f x y 为曲面顶、A 为底的曲顶柱体C 的体积.据此,用均匀随机数计算二
几何意义,它是以 ( ,
重积分的蒙特卡罗方法基本思路为：假设曲顶柱体C 包含在已知体积为 DV 的几何体 D 的内
)
部,在 D 内产生 N 个随机点,统计出在C 内部的随机点数目 k,则
I
.
≈
V
D
k
N
x y
{( ,
例：计算
I
=
∫∫
A
(1
+
1
−
2
x
−
2
y
−
2
x
+
y dxdy
2
)
,其中
A
=
) |
2
x
+
2
y
1}
≤ .
可见该积分可以看作以
+ 为顶,以 A 为底的曲顶柱体的体
积.并且它是包含在一个边长为 2 的立方体之内,很容易计算出它的精确值为π.下面我们用蒙
特卡罗方法来求解.
−
−
−
z
1
= +
2
y
2
y
2
x
1
2
x
Matlab 程序如下：
x=unifrnd(-1,1,100000,1);
y=unifrnd(-1,1,100000,1);
zi=unifrnd(0,2,100000,1);
z=1+sqrt(1-x.^2-y.^2)-sqrt(x.^2+y.^2);
k=0;
for i=1:100000
if x(i)^2+y(i)^2<=1&zi(i)<=z(i)
k=k+1;
end
end
I=8*k/100000
用该程序运行五次,其计算结果如下：
3.1522 3.1483 3.1418 3.1472 3.1437
由上述数据我们看到,每次的计算结果都在准确值π附近作很小的摆动.
3.3 蒙特卡罗方法求解三重积分
= 内部区域的面积.
1
2
2
2
0
0
2
y
z
2
y
−
+
2
x
+
π
/4
=
∫
∫
2cos
1)
x
+ 上方和球面 2
我们用实例来介绍蒙特卡罗方法求解三重积分[4].
z
(
例：计算锥面
如果我们使用球面坐标,该区域可以表示为如下的积分：
ϕ πρ ϕ ϕ ρ θ
∫
sin d d d
0
我们很容易得到该面积的精确值为π.
下面我们用蒙特卡罗方法计算该积分.
Matlab 程序如下:
n=100000;
pop=[unifrnd(0,pi/4,n,1),unifrnd(0,2,n,1),unifrnd(0,2*pi,n,1)];
k=0;
for i=1:n
if pop(i,2)<=2*cos(pop(i,1))
k=k+1;
ff(k)=pop(i,2)^2*sin((pop(i,1)));
end
- 4 -
中国科技论文在线
http://www.paper.edu.cn
end
ff;
shuzhijifen=k/n*pi^2*(sum(ff)/k)
用该程序运行五次,其计算结果如下：
3.1502 3.1464 3.1418 3.1472 3.1437
另一方面,显然该区域可以装在如下立方体的盒子里
1
− ≤ ≤ − ≤ ≤
而该立方体的体积为 8.只要生成这个盒子里均匀分布的随机点,落入该区域点的个数与
2
≤ ≤
1,0
z
x
1, 1
y
总点数之比再乘以 8 就是其体积.
Matlab 程序如下:
clear all
n=100000;
pop=[unifrnd(-1,1,n,1),unifrnd(-1,1,n,1),unifrnd(0,2,n,1)];
k=0;
for i=1:n
if sqrt(pop(i,1)^2+pop(i,2)^2)<=pop(i,3)&&pop(i,3)<=1+
sqrt(pop(i,1)^2+pop(i,2)^2)
k=k+1;
end
end
shuzhijifen=k/n*8
用该程序运行五次,其计算结果如下：
3.8887 3.8923 3.8794 3.9114 3.8904
4 小结
可积系统的一大类问题事实上是不规则区域上高维积分问题.这类问题绝大多数按数学
的研究处理模式出发,都是难以解决的,但利用 Monte-Carlo 方法却易于解决.既使是积分中被
积函数难以用解析函数表达,也可以通过抽样点的函数值及其统计实验的正交设计来逼近真
值.本文介绍了求解数值积分的蒙特卡罗方法,并结合了 Matlab 数学软件,给出写具体积分的
计算过程.而且由分析我们也可以得到 Monte-Carlo 方法的收敛速度与问题维数无关，要达到
同一精度，用蒙特卡罗方法选取的点数与维数无关，计算时间仅与维数成比例，而其他的数
值方法，达到同样的误差，点数与维数的幂次方成比例.同时蒙特卡洛法计算多重积分存在
着收敛速度慢且具有概率性质的特点,增加抽取样本次数 n,计算量增加,其精度提高很慢，但
可以通过选择适当的随机数,如有利随机数等,可使其精度提高.
随着计算机技术和数值分析技术的完善，Monte-Carlo 方法将以其简洁性和较强适应性
在优化[5]等各个领域展现它特有功效.
- 5 -
中国科技论文在线
http://www.paper.edu.cn
参考文献
[1] 鲁珊珊,李秦. Monte-Carlo 算法求解多重积分及改进[J].重庆工学院学报,2008,22(1):39-41.
[2] 宫野.计算多重积分的蒙特卡罗方法与数论网格法[J].大连理工大学学报,2001,41(1):20-23.
[3] 柴中林,银俊成.蒙特卡罗方法计算定积分的进一步讨论[J].应用数学与计算数学学报,2007,21(2):
125-128.
[4] 李庆阳,关治,白峰杉.数值计算原理[M].北京:清华大学出版社,2000:118-139.
[5] 高雷阜. Monte-Carlo 理论与优化方法的研究[J].辽宁工程技术大学学报（自然科学版）,2002，21(3):
392-394.
Monte-Carlo Method for Calculating Multiple Integral
Liaoning Technical University ,Liaoning fuxin (123000)
Wei Haiyan
Abstract
The approximate calculation of integrals has a wide range of applications in many practical engineering,
in fact a big problem of integrable system is the issue of high-dimensional points in the irregular region.
These issues are mostly difficult to solve according to the general study of mathematical models, but
the Monte-carlo method is very convenient for calculating Multiple Integral of low dimension in
regular domain of integration or complex range or even non-connected region. This paper mainly
introduces the two methods of Monte-Carlo：Random Vote Method and Average Method, and
concretely gives their principles and method steps. Besides, it combines Mathematical software Matlab
to give the solving process of Definite integral and multiple points.
Keywords: Monte-Carlo method; Multiple Integral; Random Vote Method; Average Method.
- 6 -
相关推荐
2023年江西萍乡中考道德与法治真题及答案.doc
2012年重庆南川中考生物真题及答案.doc
2013年江西师范大学地理学综合及文艺理论基础考研真题.doc
2020年四川甘孜小升初语文真题及答案I卷.doc
2020年注册岩土工程师专业基础考试真题及答案.doc
2023-2024学年福建省厦门市九年级上学期数学月考试题及答案.doc
2021-2022学年辽宁省沈阳市大东区九年级上学期语文期末试题及答案.doc
2022-2023学年北京东城区初三第一学期物理期末试卷及答案.doc
2018上半年江西教师资格初中地理学科知识与教学能力真题及答案.doc
2012年河北国家公务员申论考试真题及答案-省级.doc
2020-2021学年江苏省扬州市江都区邵樊片九年级上学期数学第一次质量检测试题及答案.doc
2022下半年黑龙江教师资格证中学综合素质真题及答案.doc
