第30卷 第5期 岩 土 工 程 学 报 Vol.30 No.5 2008年 5月 Chinese Journal of Geotechnical Engineering May, 2008 蒙特卡罗数值积分在自然单元法中的应用
李 武1，姚文娟1，朱合华2，蔡永昌2
（1.上海大学土木工程系，上海 200072；2.同济大学土木工程学院，上海 200092）
摘 要：自然单元法采用自然邻点插值方法在全域构造近似函数和试函数，该方法基于整个求解域内离散结点的Voronoi结构。当采用标准Galerkin法形成系统的平衡控制方程时，对弱形式的积分通常在Voronoi图的对偶图Delaunay 三角形内进行，但由于自然邻接插值形函数的特性，自然单元法数值积分存在明显误差。分析了自然单元法数值积分产生误差的各种可能的原因，并提出使用蒙特卡罗方法解决这一问题。该方法权系数直接与精度相关，确定方法简单有效。采用Delaunay三角形内布积分点，使得这种概率积分结果接近数学期望。给出最少积分点数的确定方法，尽可能提高蒙特卡罗积分的计算效率。通过分片试验和悬臂梁等算例验证蒙特卡罗方法解决这些误差的可行性和有效性。
关键词：自然单元法；自然相邻结点插值；蒙特卡罗方法；误差分析
中图分类号：O242.2 文献标识码：A 文章编号：1000–4548(2008)05–0698–07
作者简介：李 武(1978–)，男，黑龙江省哈尔滨人，博士，从事数值计算的研究。E-mail: jessy521@https://www.doczj.com/doc/4d12644316.html,。
Application of Monte Carlo numerical integration in natural element method
LI Wu1, YAO Wen-juan1, ZHU He-hua2, CAI Yong-chang2
(1. Department of Civil Engineering, Shanghai University, Shanghai 200072,China; 2. Department of Geotechnical Engineering, Tongji
University, Shanghai 200092, China)
Abstract: In natural element method, the trial and test functions were constructed with the natural neighbor interpolation method. The interpolation was based on the Voronoi tessellation of the scattered nodes in the problem domain. The integration of the weak form was performed in the Delaunay triangles which were the dual diagram of the Voronoi tessellation when the Galerkin method was used to form the discrete system equation. But there was obvious error in the numerical integration of the natural element method due to the characteristics of natural neighbor interpolation function. Every factor that produced the error possibly was analyzed, and a method using Monte Carlo integration was proposed to solve this problem. The weight coefficient was directly related to the precision, and its determination was also simple and effective. Integral point was cast in the Delaunay triangle, so the result of probability integral was close to mathematical expectation. The definite method of the least integral points, was given so as to enhance computational efficiency of the Monte Carlo integral as far as possible. Finally, the numerical example of the patch experiment and the cantilever beam confirmed the validity and feasibility of using Monte Carlo integration to solve these errors.
Key words: natural element method; natural neighbor interpolation; Monte Carlo method; error analysis
0 引 言
应用无网格方法求解固体力学问题时一般采用Galerkin法获得整体系统平衡方程，这就需要对“弱形式”在整个求解域内进行积分，如EFGM[1]，RPKM[2]等。目前流行的积分方法有两种：一种是基于背景网格的高斯积分方法[1-3]，另一种是无需背景网格的点积分方法[4-5]。高斯积分方法的实施依赖于对求解域进行某种形式的网格划分，然而背景网格的存在有悖于无网格方法的初衷，使得现行的大多数方法都不是真正的无网格方法[6]。另外，由于采用无网格方法构造的形函数是没有显式表达式的分式函数；无网格方法中形函数的局部支撑域与背景积分域往往不一致，从而导致积分误差[7]，因此一般采用高阶的高斯积分方法以保证积分精度，这无疑大大增加了计算量，而且精度提高不大。点积分方法无需背景网格，但采用这种方法求解稳定性较差；为保证数值求解的精度和稳定性，文献[4]提出将平衡方程的残差平方后作为稳定项加入到势能泛函数中，计算表明这种方法对于原本不
───────
基金项目：国家自然科学基金资助项目（50579093）；上海市青年科技启明星计划（06QA14050）；上海大学博士创新基金（A16-0118-07-002）收稿日期：2007–04–29
第5期李武，等. 蒙特卡罗数值积分在自然单元法中的应用699
稳定的点积分格式可以增加其稳定性，但对于原本稳定的点积分格式在加入稳定性项后则有损数值计算精度。另外，该稳定项中含有需人为给定的计算参数，该计算参数的合理选取一般需要经验或进行一定的数值试验，对于非均布结点合理确定结点的积分权重（即结点的代表域）较为困难，一般需要根据结点的分布对整个求解域进行某种几何划分，以确定每个结点所代表的“域”的大小，这个“域”在平面问题中体现为每个结点所代表的面积，在三维问题中则为结点所代表的体积。文献[7～9]提出了计算求解域内分布结点的V oronoi结构来确定结点的代表域，在二维情况下节点的代表域为该结点所对应的V oronoi多边形的面积，采用这种方法得到结点的代表域比文献[4]中的方法相对精确。
自然单元法（NEM）中采用基于求解域内离散结点V oronoi结构的自然相邻结点插值方法（NNI）构造形函数。采用标准Galerkin方法建立整体求解的系统平衡控制方程时，在V oronoi图的直线对偶图Delaunay 三角形内完成对弱形式的全域积分[10-12]。由于NEM 形函数非多项式形式，且背景积分域与NEM形函数的局部支撑域不一致从而导致积分误差[6]，在进行分片试验时其位移和应变的精度分别为1×10-4和1×10-3。
单位分解积分方法[13]是利用无网格方法形函数满足单位分解性质，将对弱形式的全域积分转化为在各个结点的影响域内进行。虽然采用单位分解积分方法的思想可以将在整个求解域内的积分转化为在形函数的局部支撑域内进行，从而无需整体的背景网格，但由于NEM 形函数的紧支域其形状较为复杂，无法直接应用现有的积分准则，故需要对其进行必要的分解映射，并且在弧段要映射大量高斯点以保证积分精度，从而导致计算量较大。
本文分析了自然单元法形函数其局部支撑域与背景积分域不一致、被积函数是有理式、动态插值形函数没有统一表达式等引起的误差，并比较这些误差对刚度矩阵K IJ的影响，采用蒙特卡罗积分[14-17]解决这些误差。
1 自然单元法（NEM）
自然相邻结点的概念建立在V oronoi图的基础上。图1中左图所示为平面7个结点的V oronoi划分及待插值点X的二阶V oronoi结构。自然相邻结点根据空圆规则确定——若Delaunay三角形的外接圆包含待插值点，则该三角形的3个顶点即为待插值点的自然相邻结点。标号为1，2，3，4的4个结点构成了待插值点X的自然相邻结点，待插值点与自然相邻点一起形成V oronoi图称为待插值点的二次V oronoi结构。
图1 Sibson和non-Sibsonian插值示意图 Fig. 1 Schematic diagram for construction of the Sibson and non-Sibsonian interpolation
NEM形函数的构造是通过采用自然相邻点插值思想Sibson或non-Sibsonian插值方法，如图1所示。
在构造NEM形函数的过程中，不涉及到矩阵的运算，而且形函数构造简单，导数的计算也相对容易，与EFGM 或RKPM等方法相比计算量大大减少。从形函数的构造和性质来看，NEM 兼有无网格的性质和有限元的优点。
2 形函数积分误差
根据文献[18]推导自然邻接形函数的计算公式：
1111
1111
()()()()
1
(,)[
2()()()()
i i i i i i
i
i i i i i i
x x x x y y y y
x y
x x y y x x y y
α????
????
??+??
=+
?????
1111
1111
()()()()
]
()()()()
i i i i i i
i i i i i i
x x x x y y y y
x x y y x x y y
++++
++++
??+??
?????
，(1)
()
()
()
I
J
I
S x
x
h x
α=，(2)
式中，()
I
S x是与结点I关联的V oronoi边的长度，()
I
h x 是插值点x到结点I的V oronoi边的垂直距离（如图1）。
1
()
()
()
I
I n
J
J
x
x
x
α
Φ
α
=
=
∑，(3)
利用式（1），（2）计算(,)
I
x y
Φ，
22
2
(1)()
(,)
(1)(21)(1)(1)
y x x x y y
x y
x y x x y x x y y
Φ
????
=
??++?+?+
，(4)
22
3
()
(,)
(1)(21)(1)(1)
xy x x y y
x y
x y x x y x x y y
Φ
?+?
=
??++?+?+
，(5)
700 岩 土 工 程 学 报 2008年
224(1)()
(,)(1)(21)(1)(1)
x y x x y y x y x y x x y x x y y Φ?+?++=??++?+?+，(6)
229(1)(1)()
(,)(1)(21)(1)(1)
x y x x y y x y x y x x y x x y y Φ?+?++=??++?+?+。(7)
由文献[13]知PUQ 积分在计算
T
d IJ I J K B DB ΩΩ
=∫ (8)
和
ext d d I I I f b t Ω
Ω
ΦΩΦΓ=+∫∫ (9)
时存在三方面误差：一方面，形函数是有理式，对有理式使用高斯积分或Hammer 积分时引起的误差；另一方面，形函数支持域与积分区域不一致引起的误差；第三，动态插值形函数没有统一表达式引起的误差。 2.1 高斯积分对有理式积分时引起误差
观察4邻结点的形函数发现形函数是有理函数，且分子和分母最高幂指数相等或分母比分子高一次。根据二元函数的泰勒公式，有理函数可以展成多项式，所以通过增加高斯点数可以减少由有理式积分引起的误差。下面通过有理函数?(x )，g (x ，y )反映这一趋势。
例如有理函数
1
()1f x x =+ (0 1 (,)1 y g x y x += + (0 2.2 形函数支持域与积分区域不一致引起的误差 高斯积分或Hammer 积分都是在背景网格上布积分点积分，而形函数的支持域却是由规则的圆弧组成的区域，它不能与背景网格完全重合。所以在使用高斯积分或Hammer 积分在背景网格上积分时引起误差。尽管这种误差可以使用单位分解积分消除——方法是在积弧形区域时过多的增加积分点，然而弧形上积分点布置和坐标转换较为烦琐。 2.3 插值形函数表达式不统一引起的积分误差 观察4邻结点的形函数发现形函数在同一个积分区域内没有统一表达式，如果只用少数积分点的表达式来代替积分区域内所有点的表达式，则引起误差，并且这个误差与积分点位置有关。如果积分点落在其权函数区域的平均值点上，就不会引起误差，否则误差将很大。又因为常用的数值积分方法的积分点位置都按一定的规律确定，所以该积分点的表达式不能完全代表所有点的被积函数的表达式，它与真实值之间 必然存在误差，所以也不能通过增加积分点来解决这一问题。 3 蒙特卡罗积分 对比各种积分方法后，通过研究发现蒙特卡罗积分方法能很好地解决上述原因带来的自然单元法的积分误差。利用蒙特卡罗方法进行积分，采用随机方法直接在矩形区域上投点，设矩形的X 方向坐标区域为[a ，b ]，Y 方向坐标区域为[c ，d ]。则首先采用随机函数rand()，产生两个[0，1]之间的随机数s 1，s 2，且X ，Y 方向生成的这两个随机数是不相关的，则所投点P 的坐标为 12(,)((),())P x y a b a s c d c s =+?+? 。 (12) 采用随机方法在任意有界区域上投点（见图2），该点函数值的数学期望就是函数在区域上的平均值。 图2 任意区域投点示意图 Fig. 2 Random point in a region 虽然随机投点的函数值的数学期望是函数的平均值，但是由于投一点的随机性较大，可以采取投多点的方式，然后计算其平均值，作为函数在区域的平均值。设该平均值为m ，则积分可以这样计算， (,)d f x y mA ΩΩ=∫ ， (13) 式中，A 为区域面积。 还可以把上式右边变化为每一点函数值乘以该区域面积除以投的点数n 的和，即 1(,)d (,)n i i i A f x y f x y n ΩΩ==∑∫ 。 (14) 这种积分方法是真正的完全无网格积分方法。它 能很好地解决由形函数是有理式引起的误差、形函数支持域与积分区域不一致引起的误差及动态插值形函数表达式不统一引起的积分误差。 4 方法的实施 根据蒙特卡罗积分理论，自然邻接点法的积分流程（以式（8）的计算为例）如下。 4.1 近似确定最小积分点数K （1）根据精度确定每个积分点的权系数，例如悬臂梁的一端y 方向位移精确到cm ，则可以采用1 cm 2 第5期 李 武，等. 蒙特卡罗数值积分在自然单元法中的应用 701 作为每个积分点权系数。 （2）由式（14）知，n =A /1（A 为积分区域面积）。 （3）为保证足够的积分点数K = nL （L 为扩大系数，可根据积分区域形状确定，一般矩形区域取为1.2）。 4.2 利用式（8）布置蒙特卡罗积分点 （1）取积分区域X 和Y 方向的最大值和最小值，并利用这四个数形成矩形区域。 （2）利用随机函数rand()生成两个数。 （3）把步骤（2）中生成的数经区间变换，转成步骤（1）矩形区域内的点。 （4）重复步骤（2）和（3）生成K 个点。 4.3 对该矩形区域的所有蒙特卡罗积分点x ξl （l =1， 2，…K ）循环 （1）判断蒙特卡罗积分点l x ξ的位置，如果位于域Ω外，则忽略该积分点，处理下一个积分点，如果位于边界上，用整形变量N 记录它，如果是位于域Ω内，则用整形变量N t 记录它。 （2）计算有效的蒙特卡罗积分点数，即N +N t 。 4.4 对所有有效蒙特卡罗积分点循环 （1）对蒙特卡罗积分点l x ξ（l =1，2，…N +N t ）循环，判断蒙特卡罗积分点l x ξ是否位于单元格内。如是，用计数器记录为N y ，则忽略该积分点，处理下一个积分点。 （2）结束对蒙特卡罗积分点的循环。 （3）利用所有位于单元格内的蒙特卡罗积分点，计算T ()()l l l w B x DB x ξξ，并将其组装到矩阵K 中。式中w l ＝y /A N （w l 为蒙特卡罗积分点l x ξ的权系数，A 为积分区域的面积，N y 为计数器记录的实际积分点数）。 4.5 结束对蒙特卡罗积分点的循环 5 数值算例 误差定义为 L = 。 (15) 5.1 不同积分方法的分片试验 通过图3的分片试验说明各种积分方法的误差大小。分析表1结果可知，Hammer 积分在积分点增加时，误差在同一量级上波动，未趋向真实解，积分时间随积分点数非线性增长。然而蒙特卡罗方法在积分点数达到50个前随着积分点数增加，误差出现波动，但在50个后随着积分点数增加，数值解趋近真实值，同样积分时间也随积分点数非线性增长。比较蒙特卡罗积分和Hammer 积分知，在相同积分点数时，蒙特 卡罗积分计算效率和精度都明显提高。蒙特卡罗积分能有效减少形函数其局部支撑域与背景积分域不一致、被积函数是有理式、动态插值形函数没有统一表达式等引起的误差，而Hammer 积分则不能减少这三种引起的误差。在自然单法中，计算形函数不需要求解逆矩阵，所以积分点数对计算效率的影响不显著。 图3 五结点分片试验图 Fig. 3 Five points lamination experiment 表1 各种积分方法分片试验 Table 1 The patch experiment of different integration methods Hammer 蒙特卡罗 积分点时间/s 误差/10-2 时间/s 误差/10-2 4 4.0 8.12 4.0 19.10 12 4.6 8.74 4.5 1.69 24 5.3 8.43 5.1 14.30 50 — — 6.6 8.85 52 6.9 8.50 6.8 7.59 100 9.4 8.32 9.1 0.59 200 — — 22.9 0.043 500 — — 31.1 0.00036 5.2 悬臂梁 一悬臂梁长2 m ，高0.3 m ，E ＝3.0×107 kPa/m 2， v =0.25，左端受集中力作用P ＝100 kN 。如图4。 图4 右端固定悬臂梁 Fig. 4 Right and fixed cantilever beam 分析图5可知蒙特卡罗积分比Hammer 积分更接近解析解，说明蒙特卡罗积分可以改进积分精度。由图6知对83（118个单元）和323（552个单元）结点数进行蒙特卡罗积分时，积分精度随着积分点的增加而提高，最终收敛于解析解[19]。而且机时消耗增长缓 慢。由图7知Hammer 积分在83和323结点时， 积分精度随着积分点的增加在某个值附近波动，未收敛于 702 岩 土 工 程 学 报 2008年 解析解。再由图8知两种积分在采用相同积分点时，蒙特卡罗的计算效率更高。 图5 在相同积分点数下数值解与解析解比较 Fig. 5 Comparison between analytic and numerical solutions 图6 蒙特卡罗积分数值解误差和积分时间 Fig. 6 The error of the numerical solution for Monte Carlo Integration and the integration time 图7 Hammer 数值解误差 Fig. 7 The error of the numerical solution for Hammer integration 图8 蒙特卡罗与Hammer 积分时间对比 Fig. 8 Comparison of the integration time between Monte Carlo and Hammer 5.3 尖劈 如图9所示，设受X 方向集中力P ＝100 kN 的尖劈，E ＝3.0×107 kPa/m 2，v =0.25，按平面应力分析。 取X 方向为2 m ，Y 方向为0.4 m ，作为分析模型。尖 劈在64结点离散情况下， 用Hammer 积分和蒙特卡罗积分方法计算得到的r 定值时，截面上的σr 向应力与解析解的比较。 图9 尖劈示意图 Fig. 9 Splitter diagram 根据弹性力学尖劈的解析解，半径相等的点径向力相等，因此取到尖点半径在区间（0.9，1.2）内的结点近似为半径相等的点，然后取这些点X 方向应力的平均值作为基准值，再用这些点的σX 与基准值相比得到在1附近波动的值。由图10知，按上述方法计算，解析解符合这个规律，蒙特卡罗方法比Hammer 积分更接近这个规律。由表2（r =1.05）知，蒙特卡罗积分 比 Hammer 积分计算效率高，这与图9反映情况一致。 图10 r =1.05的σX 应力比值图 Fig. 10 The comparative diagram of the stress at r =1.05 表2 计算时间与相对误差 Table 2 Computing time and comparing errors Hammer 蒙特卡罗 结点数 时间/s 时间/s 64 278 259 83 353 316 103 689 667 203 1209 1156 5.4 无限大圆孔板 设一个在左侧x 方向承受均匀拉伸的平板，其拉应力 1 Pa σ=，板中有半径为 1 m a =的小圆孔，板的宽度取为 5 m b =，按平面应力分析。考虑到对称性， 可取结构的四分之一作为分析模型（见图11） 。圆孔板在336结点离散情况下，用Hammer 积分和蒙特卡罗积分方法计算得到的x =0截面上的x 向正应力σx 与 解析解的比较如图12所示， 可见，对于这种应力集中 第5期 李 武，等. 蒙特卡罗数值积分在自然单元法中的应用 703 问题蒙特卡罗积分仍然能够得到比Hammer 积分更接近于解析解的数值解答。 图11 无限大圆孔板及其结点布置 Fig. 11 Node distributions for the infinite plate with hole 图12 无限大圆孔板x =0截面的正应力σx 比较 Fig. 12 The comparison of the stress at x =0 5.5 条形基础 如图13所示，一地基上的条形基础受均布压力使 土体产生大变形。将此问题简化为平面应变问题进行计算，材料采用弹塑性硬化本构关系，满足Mises 屈服准则，等向双线性硬化。材料参数为：弹性模量E =2.0×108 kPa ，泊松比v =0.25，不计体力，等效初始屈服应力σs =2.4×105 kPa ，切线模量E T =5.0×106 kPa 。 图13 条形基础简图 Fig. 13 The sketch of strip foundation 边界条件为底部固定，左右边界约束水平位移。假设条形基础为刚性的，在荷载的作用下均匀下沉，计算时在条形基础位置约束x 向位移，在y 方向上用给定位移的边界条件来模拟。由于几何模型及边界条件均对称，因此仅取右边一半建模进行计算。图14为采用自然单元法大变形程序计算的模型图，共布置了114个离散结点，采用增量法求解。 图14 结点布置与边界条件 Fig. 14 The sketch of node distributions and boundary condition 由图15知，分别采用蒙特卡罗积分和Hammer 积分对该问题进行计算，并把计算结果与有限元计算结果进行比较。由于网格出现畸变，使计算结果出现差异。根据数值积分的数学原理知，蒙特卡罗积分不受网格畸变影响，而Hammer 积分则因积分点布置在网格内受到影响。畸变网格不仅对有限元的积分有影响，而且对插值的影响更大，所以结果可能更不可靠。 图15 基础沉降4 m 时模型地表结点y 向变形曲线 Fig. 15 The deformation curves of surface nodes of the model with settlement of 4 m at y direction 6 结 论 分析了自然单元法数值积分误差产生的各种原因，通过研究提出使用蒙特卡罗积分消除这些误差是一种可行且十分有效的措施。该积分方法不受积分区域形状的限制，也不会增大刚度矩阵的维数；在有高精度计算要求的工程中，它是最节省时间的积分方法，且蒙特卡罗积分不需要背景网格，可以实现真正的无网格积分。从另一个角度讲，它保持计算效率几乎不 变的条件下，把精度提高10倍， 是一种可以通过增加 704 岩土工程学报2008年 积分点，明显改善自然单元法积分精度的有效方案，而且权系数与精度直接相关，确定方法也简单有效。在Delaunay三角内确定有效积分点，保证概率积分精度的稳定性。尤其在大变形计算时，该方法抗网格畸变的优点更为突出。 参考文献： [1] BELYTSCHKO T, LU Y Y, GU L. Element-free Galerkin methods[J]. Int J Numer Meth Engng, 1994, 37(2): 229– 256. [2] LIU W K, JUN S, ZHANG Y F. Reproducing kernel particle methods[J]. Int J Numer Meth Fluids, 1995, 20(11): 1081– 1106. [3] GONZALEZ D, CUETO E, MARTINEZ M A, DOBLARE. Numerical integration in natural neighbour Galerkin methods[J]. Int J Numer Meth Engng, 2004, 60: 2077–2104. [4] BEISSEL S, BELYTSCHKO T. Nodal integration of the element-free Galerkin method[M]. Computer Methods in Applied Mechanics and Engineering 1996, 139: 49–74. [5] RICARDO J, ALVESDE SOUSA1. A new one-point quadrature enhanced assumed strain (EAS) solid-shell element with multiple integration points along thickness, Part II: nonlinear applications[J]. Int J Numer Meth Engng, 2006, 67: 160–188. [6] CHEN J S，WU C T, YOON S, et a1. A stabilized conforming nodal integration for Galerkin meshfree methods[J]. Int J Numer Methods Engrg, 2001, 50(2): 435–466. [7] ZHOU J X, W EN J B, ZHAN G H Y, et a1. A nodal integration and post-processing technique based on V oronoi diagram for Galerkin meshless methods[J]. Comput Methods App1 Mech Engrg, 2003, 192(35): 383l–3843. [8] BRAUN J, SAMBRIDGE M. A numerical method for solving partial differential equations on highly irregular evolving grids[J]. Nature, l995, 376(24): 655–660. [9] SUKUMAR N, M ORAN B, BELYTSCHKO T. The natural element method in solid mechanics[J]. International Journal for Numerical Methods in Engineering, l998, 43(5): 839– 887. [10] SUKUMAR N, MORAN B, SEMENOV A Y, et a1. Natural neighbour Galerkin methods[J]. Int J Numer Meth Engng, 2001, 50(1): 1–27. [11] 蔡永昌, 朱合华, 王建华. 基于V oronoi结构的无网格局 部Petrov-Galerkin 方法[J]. 力学学报, 2003, 35(2): 187– l93. (CAI Yong-chang, ZHU He-hua, WANG Jian-hua. The meshless local Petrov-Galerkin method based on the V oronoi cells[J]. Acta Mechanica Sinica, 2003, 35(2): l87–193. (in Chinese)) [12] MARC D, HUNG N D. A truly meshless Galerkin method based on a moving least squares quadrature[J]. Commun Numer Meth Engng, 2002, 18(1): l–9. [13] CARPINTERI A, FEROG, VENTURAG. The partition of unity quadraturein meshfree method[J]. International Journal for Numerical Methods in Engineering, 2002, 54(7): 987– l006. [14] ACHILLEAS Lazopoulos. Error estimates in Monte Carlo and Quasi-Monte Carlo integration[J]. Acta Phys Pol Ser B, 2006, 35: 1089–1099. [15] LI Xiao-liang, YANG Jie, XIE Kai, ZHU Y M. High-quality rendering with depth cueing of volumetric data using Monte Carlo integration[J]. Computers in Biology and Medicine, 2006, 36: 1014–1025. [16] HEINZ Hofbauer, ANDREAS Uhl, PETER Zinterhof. Quasi- Monte Carlo integration in grid environments: Further leaping effects[J]. Parallel Processing Letters, 2006, 16: 285 –311. [17] MARC Fasnacht, SWENDSEN Robert H, ROSENBERG John M. Adaptive integration method for Monte Carlo simulations[J]. Physical Review E, 2004, 69(056704): 1–15. [18] 王兆淸, 冯伟. Delaunay多边形单元的有理函数插值格 式[J]. 力学季刊, 2004, 25(3): 375–381. (WANG Zhao-qing, FENG Wei. Rational function interpolation scheme of Delaunay ploygonal elements[J]. Chinese Quarterly of Mechanics, 2004, 25(3): 375–381. (in Chinese)) [19] TIMOSHENKO S P, GOODIER J N. Theory of Elasticity[M]. 3rd ed. New York: McGraw Hill, 1970. 文献综述 信息与计算科学 蒙特卡罗方法的应用 在解决实际问题的时候, 为了模拟某一过程, 产生各种概率分布的随机变量和对于那些由于计算过于复杂而难以得到解析解或者根本没有解析解的问题, 我们应该怎么办? 蒙特·卡罗是一种十分有效的求出数值解的方法. 蒙特卡罗法（ monte-carlo method ）简称M -C 法 通过构造概率模型并对它进行随机试验来解算数学问题的方法. 以计算函数的定积分()()1 0I f x d x =?, ()01f x ≤≤为例, 首先构造一个概率模型: 取一个边长分别为和-的矩形, 并在矩形内随机投点M , 假设随机点均匀地落在整个矩形之内, 当点的掷点数N 充分大时, 则落在图中阴影区内的随机点数与投点总数N 之比M N 就近似等于积分值I . 蒙特卡罗法历史悠久. 1773年法国G.-L.L.von 布丰曾通过随机投针试验来确定圆周率π的近似值, 这就是应用这个方法的最早例子. 蒙特卡罗是摩纳哥著名赌城, 1945年 J.von 诺伊曼等人用它来命名此法, 沿用至今. 数字计算机的发展为大规模的随机试验提供了有效工具, 遂使蒙特卡罗法得到广泛应用. 在连续系统和离散事件系统的仿真中, 通常构造一个和系统特性相近似的概率模型, 并对它进行随机试验, 因此蒙特卡罗法也是系统仿真方法之一. 蒙特卡罗法的步骤是: 构造实际问题的概率模型; ②根据概率模型的特点, 设计和使用降低方差的各类方法, 加速试验的收敛; ③给出概率模型中各种不同分布随机变量的抽样方法; ④统计试验结果, 给出问题的解和精度估计. 概率模型用概率统计的方法对实际问题或系统作出的一种数学描述. 例如对离散事件系统中临时实体的到达时间、永久实体的服务时间的描述（见离散事件系统仿真方法）就是采用概率模型. 虽然由这些模型所确定的到达时间、服务时间可能与具体某一段时间内实际到达时间、服务时间有出入, 但它是通过多次统计获得的结果, 所以从概率分布的规律来说还是相符的. 概率模型不仅可用来描述本身就具有随机特性的问题或系统, 也可用来描述一个确定型问题. 例如参数寻优中的随机搜索法（见动力学系统参数寻优）就是将参数最优化问题构造为一个概率模型, 然后用随机投点、统计分析的方法来进行搜索. 浅析蒙特卡洛方法原理及应用 于希明 （英才学院1236103班测控技术与仪器专业6120110304） 摘要：本文概述了蒙特卡洛方法产生的历史及基本原理，介绍了蒙特卡洛方法的最初应用——蒲丰投针问题求圆周率，并介绍了蒙特卡洛方法在数学及生活中的一些简单应用，最后总结了蒙特卡洛方法的特点。 关键词：蒙特卡洛方法蒲丰投针生活应用 蒙特卡洛方法（Monte Carlo method），也称统计模拟方法，是二十世纪四十年代中期由于科学技术的发展和电子计算机的发明，而被提出的一种以概率统计理论为指导的一类非常重要的数值计算方法。它是以概率统计理论为基础, 依据大数定律( 样本均值代替总体均值) , 利用电子计算机数字模拟技术, 解决一些很难直接用数学运算求解或用其他方法不能解决的复杂问题的一种近似计算法。蒙特卡洛方法在金融工程学，宏观经济学，计算物理学（如粒子输运计算、量子热力学计算、空气动力学计算）等领域应用广泛。 一、蒙特卡洛方法的产生及原理 蒙特卡洛方法于20世纪40年代美国在第二次世界大战中研制原子弹的“曼哈顿计划”计划的成员S.M.乌拉姆和J.冯·诺伊曼首先提出。数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种方法，为它蒙上了一层神秘色彩。在这之前，蒙特卡洛方法就已经存在。1777年，法国数学家蒲丰（Georges Louis Leclere de Buffon，1707—1788）提出用投针实验的方法求圆周率π。这被认为是蒙特卡洛方法的起源。 其基本原理如下：由概率定义知，某事件的概率可以用大量试验中该事件发生的频率来估算，当样本容量足够大时，可以认为该事件的发生频率即为其概率。因此，可以先对影响其可靠度的随机变量进行大量的随机抽样，然后把这些抽样值一组一组地代入功能函数式，确定结构是否失效，最后从中求得结构的失效概率。蒙特卡洛法正是基于此思路进行分析的。 设有统计独立的随机变量Xi(i=1，2，3，…，k)，其对应的概率密度函数分别为fx1，fx2，…，fxk，功能函数式为Z=g(x1，x2，…，xk)。首先根据各随机变量的相应分布，产生N组随机数x1，x2，…，xk值，计算功能函数值Zi=g(x1，x2，…，xk)(i=1，2，…，N)，若其中有L组随机数对应的功能函数值Zi≤0，则当N→∞时，根据伯努利大数定理及正态随机变量的特性有：结构失效概率，可靠指标。 二、蒲丰投针问题 作为蒙特卡洛方法的最初应用, 是解决蒲丰投针问题。1777 年, 法国数学家蒲丰提出利用投针实验求解圆周率的问题。设平面上等距离( 如为2a) 画有一些平行线, 将一根长度为2l( l< a) 的针任意投掷到平面上, 针与任一平行线相交的频率为p 。针的位置可以用针的中心坐标x 和针与平行线的夹角θ来决定。任意方向投针, 便意味着x与θ可以任意取一值, 只是0≤x ≤a, 0≤θ≤π。那么, 投针与任意平行线相交的条件为x ≤ l sinθ。相交频率p 便可用下式求 蒙特卡洛方法 1、蒙特卡洛方法的由来 蒙特卡罗分析法（Monte Carlo method），又称为统计模拟法，是一种采用随机抽样（Random Sampling）统计来估算结果的计算方法。由于计算结果的精确度很大程度上取决于抽取样本的数量，一般需要大量的样本数据，因此在没有计算机的时代并没有受到重视。 第二次世界大战时期，美国曼哈顿原子弹计划的主要科学家之一，匈牙利美藉数学家约翰·冯·诺伊曼（现代电子计算机创始人之一）在研究物质裂变时中子扩散的实验中采用了随机抽样统计的手法，因为当时随机数的想法来自掷色子及轮盘等赌博用具，因此他采用摩洛哥著名赌城蒙特卡罗来命名这种计算方法，为这种算法增加了一层神秘色彩。 蒙特卡罗方法提出的初衷是用于物理数值模拟问题, 后来随着计算机的快速发展, 这一方法很快在函数值极小化、计算几何、组合计数等方面得到应用, 于是它作为一种独立的方法被提出来, 并发展成为一门新兴的计算科学, 属于计算数学的一个分支。如今MC方法已是求解科学、工程和科学技术领域大量应用问题的常用数值方法。 2、蒙特卡洛方法的核心—随机数 蒙特卡洛方法的基本理论就是通过对大量的随机数样本进行统计分析，从而得到我们所需要的变量。因此蒙特卡洛方法的核心就是随机数，只有样本中的随机数具有随机性，所得到的变量值才具有可信性和科学性。 在连续型随机变量的分布中, 最基本的分布是[0, 1]区间上的均匀分布, 也称单位均匀分布。由该分布抽取的简单子样ξ1，ξ2ξ3……称为随机数序列, 其中每一个体称为随机数, 有时称为标准随机数或真随机数, 独立性和均匀性是其必备的两个特点。真随机数是数学上的抽象, 真随机数序列是不可预计的, 因而也不可能重复产生两个相同的真随机数序列。真随机数只能用某些随机物理过程来产生, 如放射性衰变、电子设备的热噪音、宇宙射线的触发时间等。 实际使用的随机数通常都是采用某些数学公式产生的，称为伪随机数。真随机数只是一种数学的理想化概念，实际中我们所接触到的和使用的都是伪随机数。要把伪随机数当成真随机数来使用, 必须要通过随机数的一系列的统计检验。 无论伪随机数用什么方法产生,它的局限性都在于这些随机数总是一个有限长的循环集合, 而且序列偏差的上确界达到最大值。所以若能产生低偏差的确定性序列是很有用的,产生的序列应该具有这样的性质, 即任意长的子序列都能均匀地填充函数空间。 人们已经产生了若干种满足这个要求的序列,如Halton序列、Faure序列、Sobol序列和Niederreiter序列等。称这些序列为拟随机数序列。伪随机序列是为了模拟随机性, 而拟随机序列更致力于均匀性。 3、蒙特卡洛方法的原理 当问题可以抽象为某个确定的数学问题时，应当首先建立一个恰当的概率模型，即确定某个随机事件A或随机变量X，使得待求的解等 《数学实验》报告 班级：序号：姓名： 1．问题描述 I、用蒙特卡罗方法计算以下函数在区间上的积分，并改变随机点数 目观察对结果的影响。 （1）y=1/(1+x), 0= 0.6975 >> p=shell1(0,1,10000) p = 0.6922 >> p=shell1(0,1,100) p = 0.7001 >> p=shell1(0,1,500) p = 0.6890 结果分析：改变了四次随机点数，结果都趋近于0.69，说明积分值约等于 0.69，但是点数越多，值越接近。 I、(2)使用均值估计法 程序： function p=shell2(a,b,n) z=0; x=unifrnd(a,b,1,n); fori=1:n u=(exp(3*x(i)))*sin(2*x(i)); z=z+u; end p=(b-a)*z/n; 运行结果： >> p=shell2(0,2,1000) p = -24.4911 >> p=shell2(0,2,100) p = -43.8720 >> p=shell2(0,2,10000) p = -30.8699 >> p=shell2(0,2,500) p = -23.2955 >> p=shell2(0,2,100000) p = 图1-1 蒙特卡罗方法学习总结 核工程与核技术2014级3班张振华20144530317 一、蒙特卡罗方法概述 1.1蒙特卡罗方法的基本思想 1.1.1基本思想 蒙特卡罗方的基本思想就是，当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。 1.1.2计算机模拟打靶游戏 为了能更为深刻地理解蒙特卡罗方法的基本思想，我们学习了蒲丰氏问题和打靶游戏两大经典例子。下面主要对打靶游戏进行剖析、计算机模拟（MATLAB 程序）。 设某射击运动员的弹着点分布如表1-1 所示， 首先用一维数轴刻画出已知该运动员的弹 着点的分布如图1-1所示。研究打靶游戏，我 们不用考察子弹的运动轨迹，只需研究每次“扣动扳机”后的子弹弹着点。每一环数对应唯一确定的概率，且注意到概率分布函数有单调不减和归一化的性质。首先我们产生一个在(0,1)上均匀分布的随机数（模拟扣动扳机），然后将该随机数代表的点投到P 轴上（模拟子弹射向靶上的一个确定点），得到对应的环数（即子弹的弹着点），模拟打靶完成。反复进行N 次试验，统计出试验结果的样本均值。样本均值应当等于数学期望值，但允许存在一定的偏差，即理论计算值应该约等于模拟试验结果。 clear all;clc; N=100000;s=0; for n=1:N %step 4.重复N 次打靶游戏试验 x=rand(); %step 1.产生在(0,1)上均匀分布的随机数if(x<=0.1) %step 2.若随机数落在(0.0,0.1)上，则代表弹着点在7环g=7; s=s+g; %step 3.统计总环数elseif(x<=0.2) %step 2.若随机数落在(0.1,0.2)上，则代表弹着点在8环g=8;s=s+g; elseif(x<=0.5) %step 2.若随机数落在(0.2,0.5)上，则代表弹着点在9环g=9;s=s+g; else %step 2.若随机数落在(0.5,1.0)上，则代表弹着点在10环 g=10;s=s+g; end end gn_th=7*0.1+8*0.1+9*0.3+10*0.5; %step 5.计算、输出理论值fprintf('理论值：%f\n',gn_th); gn=s/N; %step 6.计算、输出试验结果 fprintf('试验结果：%f\n',gn);1.2蒙特卡罗方法的收敛性与误差 1.2.1收敛性 由大数定律可知，应用蒙特卡罗方法求近似解，当随机变量Z 的简单子样数N 趋向于无穷大（N 充分大）时，其均值依概率收敛于它的数学期望。 1.2.2误差 由中心极限定理可知，近似值与真值的误差为N Z E Z N αλ<-)(?。式中的αλ的值可以根据给出的置信水平，查阅标准正态分布表来确定。 1.2.3收敛性与误差的关系 在一般情况下，求具有有限r 阶原点矩()∞ 蒙特卡罗方法（MC） 蒙特卡罗（Monte Carlo）方法： 蒙特卡罗（Monte Carlo）方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在本世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。 传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。这也是我们采用该方法的原因。 蒙特卡罗方法的基本原理及思想如下： 当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并 用它们作为问题的解。这就是蒙特卡罗方法的基本思想。蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。它是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。可以把蒙特卡罗解题归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。 蒙特卡罗解题三个主要步骤： 构造或描述概率过程： 对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。 实现从已知概率分布抽样： 构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。 建立各种估计量： 一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。 例如：检验产品的正品率问题，我们可以用1表示正品，0表示次品，于是对每个产品检验可以定义如下的随机变数Ti，作为正品率的估计量： 于是，在N次实验后，正品个数为： 概率论实验报告 ——蒙特卡洛方法估计积分值 姓名： 学号： 班级： 实验内容：用蒙特卡洛方法估计积分值 1用蒙特卡洛方法估计积分 20sin x xdx π ?，2-0x e dx +∞?和 22221x y x y e dxdy ++≤??的值，并将估 计值与真值进行比较。 2用蒙特卡洛方法估计积分 21 0x e dx ? 和 22x y +≤??的值， 并对误差进行估计。 要求：（1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方法； （2）利用计算机产生所选分布的随机数以估计积分值； （3）进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误差（针对第1类题）或样本方差（针对第2类题）以评价估计结果的精度。 目的：（1）能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数及其期望、方差、协方差等； （2） 熟练使用 MATLAB 对样本进行基本统计，从而获取数据的基本信息； （3） 能用 MATLAB 熟练进行样本的一元回归分析。 实验一、估计2 sin x xdx π ?的值，并将估计值与真值进行比较。 MATLAB 代码： s=0;m=0;f=0;r=0;n=50; h(1:10)=0; for j=1:10 for i=1:n a=unifrnd(0,pi/2,n,1); x=sort(a); y=pi/2*mean(x.*sin(x)); s=s+y; end b=s./n; fprintf('b=%.4f\n',b); h(j)=b; s=0; m=m+b; end p=m./10 z=1 for j=1:10 r=(h(j)-z).^2; f=f+r; end f=f./10; fprintf('f=%.6f\n',f) 运行结果： b=1.0026 b=1.0061 b=1.0037 b=1.0135 b=0.9932 b=0.9988 b=1.0213 b=1.0310 b=0.9813 b=1.0041 p = 1.0056 z = 1 f=0.000207 >> （运行截图） 结果显示f=0.000207，表明估计结果与理论值非常接近。 实验二、估计 2-0x e dx +∞ ?的值，并将估计值与真值进行比较。 I=dx e x ?+∞-02=1/2*pi dx e pi e x x *2***2/1*2/2/22-+∞ ∞--? =)(x f x 2/2**2/1x e pi - g(x)=e pi x *2*2/2- )(x f x 为标准正态分布的概率密度.分别取10个估计值h(j)，求得估计值的均值p ，对照积分的真实值求得估计均方误差f 。 蒙特卡罗方法及应用 实验讲义 东华理工大学核工系 2016.8 实验一 蒙特卡罗方法基本思想 一、实验目的 1、了解蒙特卡罗方法方法的基本思想； 2、掌握蒙特卡罗方法计算面积、体积的方法； 3、掌握由已知分布的随机抽样方法。 二、实验原理 Monte Carlo 方法，又称统计模拟方法或计算机随机模拟方法，是一种基于“随机数”进行数值模拟的方法，一种采用统计抽样理论近似求解物理或数学问题的方法。 如待求量可以表述成某些特征量的期望值、某些事件出现的概率或两者的函数形式，那么可采用蒙特卡罗方法求解。在求解某些特征量的期望值或某些事件出现的概率时，必须构建合符实际的数学模型。例如采用蒙特卡罗方法计算某函数所围面积时，构建的数学模型是构造一已知面积的可均匀抽样区域，在该区域投点，由伯努利定理大数定理可知，进入待求区域投点的频率依概率1收敛于该事件出现的概率（面积之比）。 由已知分布的随机抽样方法指的是由已知分布的总体中抽取简单子样。具体方法很多，详见教材第三章。 三、实验内容 1、安装所需计算工具（MATLAB 、fortran 、C++等）； 2、学习使用rand(m,n)、unifrnd(a,b,m,n)函数 3、求解下列问题： 3.0、蒲丰氏投针求圆周率。 3.1、给定曲线y =2 – x 2 和曲线y 3 = x 2，曲线的交点为：P 1( – 1，1 )、P 2( 1，1 )。曲线围成平面有限区域，用蒙特卡罗方法计算区域面积； 3.2 、计算1z z ?≥??≤??所围体积 其中{(,,)|11,11,02}x y z x y z Ω=-≤≤-≤≤≤≤。 4、对以下已知分布进行随机抽样： 统计计算课程设计 题目基于蒙特卡洛方法求数值积分 中文摘要 蒙特卡洛方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在上世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。传统的经验方法由于不能逼近真实的 物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问 题与实际非常符合，可以得到很圆满的结果。 利用随机投点法，平均值法，重要性采样法，分层抽样法，控制变量法，对偶变量法，运用R 软件求 1 d x e x θ- =?，42d x e x θ- =?和12 d 1 x e x x θ - = + ?数值积分。计算以上各种估计的方差，给出 精度与样本量的关系，比较各种方法的效率， 关键字蒙特卡洛随机投点法平均值法R软件 1 绪论 蒙特卡洛的基本思想是，当所求解问题是某种随机事件出现的概率，或者是某个随机变量的期望值时，通过某种“实验”的方法，以这种事件出现的频率估计这一随机事件的概率，或者得到这个随机变量的某些数字特征，并将其作为问题的解。 蒙特卡洛方法解题过程的三个主要步骤： （1）构造或描述概率过程 对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。 （2）实现从已知概率分布抽样 构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡洛方法模拟实验的基本手段，这也是蒙特卡洛方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是（0，1）上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡洛模拟的基本工具。 （3）建立各种估计量 一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。 ? 陈强,《高级计量经济学及Stata 应用》课件，第二版，2014 年，高等教育出版社。 第 19 章蒙特卡罗法与自助法 19.1 蒙特卡罗法的思想与用途 通过计算机模拟从总体抽取大量随机样本的计算方法统称为“蒙特卡罗法”(Monte Carlo Methods，简记MC)。 例（计算圆周率π）：在边长为1 的正方形中内接1单位圆。正方形面积为1，1 4圆面积为π 4。如知道1 4单位圆占 正方形面积的比例，就可计算π。 图19.1 计算圆周率 的随机实验 向这个正方形随机地射箭，落点在正方形上服从二维均匀分布。重复实验n 次，其中有m 次落在1 4圆内。 2 3 ? ? 根据大数定律，m n ?p ?→π 4，故π≈ 4m n 。 在计量中，常用 MC 来确定统计量的小样本性质。 【例】对于y i = x i 'β + εi (i = 1, , n )，对H 0 : R β = r 进行显著性水平 为 5%的大样本检验： W ≡ n (R β? - r )' ? R A var(β?)R '?-1 (R β? - r ) ?d ?→ χ 2 (m ) 其中β? 为 OLS 估计量，m 为线性约束个数。 渐近χ 2 分布只是真实分布的近似，故“5%”可能只是“名义显著 性水平”(nominal size)，而非“真实显著性水平”(true or actual size)，二者之差称为“显著性水平扭曲”(size distortion)。 可用MC 来确定“真实显著性水平”。 第一步，给定β的具体取值，以及x 与的概率分布。第二步，从x 与的分布中随机抽样，得到{x1, {ε 1, ε2 , , εn }。 x 2 , , x n }与 第三步，根据方程y i=x i'β +εi 计算{y1,y2 , , y n }。 第四步，对此样本进行OLS 估计，计算统计量W ，与χ2 (m)的5%临界值比较，确定是否拒绝原假设H0 : Rβ =r 。 第五步，大量重复第二至第四步，得到M 个随机样本(比如，M =1 000)，进行M 次检验，则拒绝原假设的比例就是真实显著性 4 第三章蒙特卡罗方法简介 3.1 Monte Carlo方法简介 Monte Carlo方法是诺斯阿拉莫斯实验室在总结其二战期间工作（曼哈顿计划）的基础上提出来的。Monte Carlo的发明，主要归功于Enrico Fermi、Von Neumann和Stanislaw Ulam等。自二战以来，Monte Carlo方法由于其在解决粒子输运问题上特有的优势而得到了迅速发展，并在核物理、辐射物理、数学、电子学等方面得到了广泛的应用。Monte Carlo的基本思想就是基于随机数选择的统计抽样，这和赌博中掷色子很类似，故取名Monte Carlo。 Monte Carlo方法非常适于解决复杂的三维问题，对于不能用确定性方法解决的问题尤其有用，可以用来模拟核子与物质的相互作用。在粒子输运中，Monte Carlo技术就是跟踪来自源的每个粒子，从粒子产生开始，直到其消亡（吸收或逃逸等）。在跟踪过程中，利用有关传输数据经随机抽样来决定粒子每一步的结果[6]。 3.2 Monte Carlo发展历程 MCNP程序全名为Monte Carlo Neutron and Photon Transport Code (蒙特卡罗中子-光子输运程序)。Monte Carlo模拟程序是在1940年美国实施“发展核武器计划”时，由洛斯阿拉莫斯实验室（LANL）提出的，为其所投入的研究、发展、程序编写及参数制作超过了500人年。1950年Monte Carlo方法的机器语言出现, 1963年通用性的Monte Carlo方法语言推出，在此基础上，20世纪70年代中期由中子程序和光子程序合并,形成了最初的MCNP程序。自那时起，每2—3年MCNP更新一次, 版本不断发展，功能不断增加，适应面也越来越广。已知的MCNP程序研制版本的更新时间表如下：MCNP-3：1983年写成，为标准的FORTRAN-77版本，截面采用ENDF /B2III。 MCNP-3A：1986年写成，加进了多种标准源，截面采用ENDF /B2I V[20]。 用蒙特卡洛方法估计积分方法及matlab编程实现 专业班级：材料43 学生姓名：王宏辉 学号：2140201060 指导教师：李耀武 完成时间：2016年6月8日 用蒙特卡洛方法估计积分 方法及matlab 编程实现 实验内容： 1用蒙特卡洛方法估计积分 2 sin x xdx π ?，2 -0 x e dx +∞ ?和 2 2 221 x y x y e dxdy ++≤?? 的值， 并将估计值与真值进行比较。 2用蒙特卡洛方法估计积分 2 1 x e dx ?和 224 4 1 11x y dxdy x y +≤++?? 的值， 并对误差进行估计。 要求： （1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方法； （2）利用计算机产生所选分布的随机数以估计积分值； （3）进行重复试验，通过计算样本均值以评价估计的无偏性；通过计算均方误差（针 对第1类题）或样本方差（针对第2类题） 以评价估计结果的精度。 目的： （1）能通过 MATLAB 或其他数学软件了解随机变量的概率密度、分布函数 及其期望、方差、协方差等； （2） 熟练使用 MATLAB 对样本进行基本统计，从而获取数据的基本信息； （3） 能用 MATLAB 熟练进行样本的一元回归分析。 实验原理： 蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从而估计相应的积分值。 具体操作如下： 一般地，积分?=b dx x g a )(S 改写成??==b b dx f h dx f g a a )(x )(x )(x f(x)) (x S 的形 式，（其中为)f(x 一随机变量X 的概率密度函数，且)f(x 的支持域） （｝｛b f ,a 0)(x |x ?>)，f(x) ) (x )(x g h =）;令Y=h(X)，则积分S=E （Y ）；利用matlab 软件，编程产生随机变量X 的随机数，在由 ?? ??∈==) b (a,,) b (a,,01I(x) ,)(x )(x y x x I h ，得到随机变量Y 的随机数，求出样本均值，以此估计积分值。 积分??=A dxdy g S )y (x,的求法与上述方法类似，在此不赘述。 概率密度函数的选取： 一重积分，由于要求)f(x 的支持域） （｝｛b f ,a 0)(x |x ?>，为使方法普遍适用，考虑到标准正态分布概率密度函数22 e 21)(x x f -=π 支持域为 R ，故选用2 2e 21)(x x f - = π 。 类似的，二重积分选用2 2 221)y (x,y x e f +-=π ，支持域为2R 。 估计评价： 用蒙特卡洛方法估计积分 方法及matlab编程实现 专业班级：材料43 学生姓名：王宏辉 学号： 指导教师：李耀武 完成时间：2016年6月8日 用蒙特卡洛方法估计积分 方法及matlab编程实现 实验内容: 31 2 乂 1用蒙特卡洛方法估计积分xsin xdx , e-x dx和II e x y dxdy的值， 0 0 x2：:;y2 d 并将估计值与真值进行比较。 2用蒙特卡洛方法估计积分e x dx和 __ dxdy的值，并对误 0 x2#g J l + x4+ y4 差进行估计。 要求： （1）针对要估计的积分选择适当的概率分布设计蒙特卡洛方 法； （2）利用计算机产生所选分布的随机数以估计积分值； （3）进行重复试验，通过计算样本均值以评价估计的无偏性； 通过计算均方误差（针对第1类题）或样本方差（针对第2类题）以评价估计结果的精度。 目的： （1）能通过MATLAB或其他数学软件了解随机变量的概率密 度、分布函数及其期望、方差、协方差等； （2）熟练使用MATLAB对样本进行基本统计，从而获取数据的 基本信息； （3）能用MATLAB熟练进行样本的一元回归分析 实验原理： 蒙特卡洛方法估计积分值，总的思想是将积分改写为某个随机变 量的数学期望，借助相应的随机数，利用样本均值估计数学期望，从 而估计相应的积分值。 具体操作如下: 式，(其中为f(x) 一随机变量 X 的概率密度函数，且f(x)的支持域 {x|f(x)〉O}=( a,b)) , h(x^-g ^x) )；令 Y=h(X)，则积分 S=E( Y ；利用 f(x) matlab 软件，编程产生随机变量X 的随机数，在由 y = h(x)l(x), l(x) =」 公匸⑻①，得到随机变量丫的随机数，求出样本均 0 X(a,b) 值，以此估计积分值。 积分S = .. g(x, y)dxdy 的求法与上述方法类似，在此不赘述。 A 概率密度函数的选取： 一重积分，由于要求f(x)的支持域｛x|f(x)?0｝二(a,b)，为使方法普 1 故选用 f(x F 彳 ￥廿 类似的，二重积分选用f(x, y) — e 2，支持域为R 2 2兀 估计评价: 进行重复试验，通过计算样本均值以评价估计的无偏性；通过计 算均方误（针对第1类题，积得出）或样本方差（针对第2类题，积 不出）以评价 般地，积分 b S = g(x)dx 改写成 a b g(x) a f(X ) b f (x)dx = j h(x) f (x)dx 的形 遍适用，考虑到标准正态分布概率密度函数 2 x _ f (x)二-^—e 2支持域为 蒙特卡罗方法及其在中子输运问题中得应用 目录 蒙特卡罗方法及其在中子输运问题中得应用 (1) 1蒙特卡罗方法简介 (3) 1.1蒙特卡罗方法的基本原理 (3) 1.2 蒙特卡罗方法的误差 (4) 2 随机变量的抽样方法 (4) 2.1 直接抽样方法 (5) 2.1.1 离散型随机变量的抽样 (5) 2.1.2 连续型随机变量的抽样 (5) 2.2 挑选抽样法 (5) 2.3 复合抽样法 (6) 3 蒙特卡罗方法模拟中子输运过程 (6) 3.1 源抽样 (6) 3.2 输运距离的抽样 (7) 3.3 碰撞核素的抽样值 (7) 3.4 反应类型的抽样值 (7) 3.5 反应后中子状态的确定 (7) 3.5.1 弹性散射 (7) 3.5.2 非弹性散射 (8) 3.5.3 裂变反应 (8) 4 蒙特卡罗方法的减方差技巧 (8) 4.1 权 (8) 4.2 统计估计法 (9) 4.3 权窗 (10) 5 蒙特卡罗方法求解通量 (10) 5.1 通量的定义 (10) 5.2 点通量的计算 (11) 5.3 面通量的计算 (11) 5.3.1 统计估计法 (11) 5.3.2 加权法 (12) 5.4 体通量的计算 (12) 5.4.1 统计估计法 (12) 5.4.2 径迹长度法 (13) 5.4.3 碰撞密度法 (13) 5.4.4 几种体通量计算方法的比较 (14) 5.5 最终结果的统计 (14) 6 蒙特卡罗方法求解k eff (15) 6.1 有效增值因子k eff的定义 (15) 6.2 蒙特卡罗方法求解k eff (15) 6.2.1 吸收估计法 (15) 6.2.2 碰撞估计法 (15) 6.2.3 径迹长度估计法 (16) 第一章：Monte Carlo方法概述 一、Monte Carlo历史渊源 Monte Carlo方法的实质是通过大量随机试验，利用概率论解决问题的一种数值方法，基本思想是基于概率和体积间的相似性。它和Simulation有细微区别。单独的Simulation只是模拟一些随机的运动，其结果是不确定的；Monte Carlo在计算的中间过程中出现的数是随机的，但是它要解决的问题的结果却是确定的。 历史上有记载的Monte Carlo试验始于十八世纪末期（约1777年），当时布丰（Buffon）为了计算圆周率，设计了一个“投针试验”。（后文会给出一个更加简单的计算圆周率的例子）。虽然方法已经存在了200多年，此方法命名为Monte Carlo则是在二十世纪四十年，美国原子弹计划的一个子项目需要使用Monte Carlo方法模拟中子对某种特殊材料的穿透作用。出于保密缘故，每个项目都要一个代号，传闻命名代号时，项目负责人之一von Neumann灵犀一点选择摩洛哥著名赌城蒙特卡洛作为该项目名称，自此这种方法也就被命名为Monte Carlo 方法广为流传。 十一、Monte Carlo方法适用用途 （一）数值积分 计算一个定积分，如，如果我们能够得到f(x)的原函数F(x)，那么直接由表达式: F(x1)-F(x0)可以得到该定积分的值。但是，很多情况下，由于f(x)太复杂，我们无法计算得到原函数F(x)的显示解，这时我们就只能用数值积分的办法。如下是一个简单的数值积分的例子。 数值积分简单示例 如图，数值积分的基本原理是在自变量x的区间上取多个离散的点，用单个点的值来代替该小段上函数f(x)值。 常规的数值积分方法是在分段之后，将所有的柱子（粉红色方块）的面积全部加起来，用这个面积来近似函数f(x)（蓝色曲线）与x轴围成的面积。这样做当然是不精确的，但是随着分段数量增加，误差将减小，近似面积将逐渐逼近真实的面积。 Monte Carlo数值积分方法和上述类似。差别在于，Monte Carlo方法中，我们不需要将所有方柱的面积相加，而只需要随机地抽取一些函数值，将他们的面积累加后计算平均值就够了。通过相关数学知识可以证明，随着抽取点增加，近似面积也将逼近真实面积。 在金融产品定价中，我们接触到的大多数求基于某个随机变量的函数的期望值。考虑一个欧式期权，假定我们已经知道在期权行权日的股票服从某种分布（理论模型中一般是正态分布），那么用期权收益在这种分布上做积分求期望即可。 （五）随机最优化 蒙特卡罗方法的解题过程可以归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。 蒙特卡罗方法解题过程的三个主要步骤： （1）构造或描述概率过程 对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。即要将不具有随机性质的问题转化为随机性质的问题。 （2）实现从已知概率分布抽样 构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。最简单、最基本、最重要的一个概率分布是（0，1）上的均匀分布（或称矩形分布）。随机数就是具有这种均匀分布的随机变量。随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。产生随机数的问题，就是从这个分布的抽样问题。在计算机上，可以用物理方法产生随机数，但价格昂贵，不能重复，使用不便。另一种方法是用数学递推公式产生。这样产生的序列，与真正的随机数序列不同，所以称为伪随机数，或伪随机数序列。不过，经过多种统计检验表明，它与真正的随机数，或随机数序列具有相近的性质，因此可把它作为真正的随机数来使用。由已知分布随机抽样有各种方法，与从(0,1)上均匀分布抽样不同，这些方法都是借助于随机序列来实现的，也就是说，都是以产生随机数为前提的。由此可见，随机数是我们实现蒙特卡罗模拟的基本工具。 （3）建立各种估计量 一般说来，构造了概率模型并能从中抽样后，即实现模拟实验后，我们就要确定一个随机变量，作为所要求的问题的解，我们称它为无偏估计。建立各种估计量，相当于对模拟实验的结果进行考察和登记，从中得到问题的解。 蒙特卡洛法模拟蒲丰（Buffon）投针实验-使用Matlab 2010年03月31日星期三8:47 蒲丰投针实验是一个著名的概率实验，其原理请参见此页： https://www.doczj.com/doc/4d12644316.html,/reese/buffon/buffon.html 现在我们利用Matlab来做模拟，顺便说一下，这种随机模拟方法便是传说中的“蒙特- 蒙特卡洛方法及其应用 1风险评估及蒙特卡洛方法概述 １.１蒙特卡洛方法。 蒙特卡洛方法，又称随机模拟方法或统计模拟方法，是在20世纪40年代随着电子计算机的发明而提出的。它是以统计抽样理论为基础，利用随机数，经过对随机变量已有数据的统计进行抽样实验或随机模拟，以求得统计量的某个数字特征并将其作为待解决问题的数值解。 蒙特卡洛模拟方法的基本原理是：假定随机变量X1、X2、X3……X n、Y，其中X1、X2、X3……X n 的概率分布已知，且X1、X2、X3……X n、Y有函数关系：Y=F（X1、X2、X3……X n），希望求得随机变量Y的近似分布情况及数字特征。通过抽取符合其概率分布的随机数列X1、X2、X3……X n带入其函数关系式计算获得Y的值。当模拟的次数足够多的时候，我们就可以得到与实际情况相近的函数Y的概率分布和数字特征。 蒙特卡洛法的特点是预测结果给出了预测值的最大值,最小值和最可能值，给出了预测值的区间范围及分布规律。 １.２风险评估概述。 风险表现为损损益的不确定性，说明风险产生的结果可能带来损失、获利或是无损失也无获利，属于广义风险。正是因为未来的不确定性使得每一个项目都存在风险。对于一个公司而言，各种投资项目通常会具有不同程度的风险，这些风险对于一个公司的影响不可小视，小到一个项目投资资本的按时回收，大到公司的总风险、公司正常运营。因此，对于风险的测量以及控制是非常重要的一个环节。 风险评估就是量化测评某一事件或事物带来的影响的可能程度。根据“经济人”假设，收益最大化是投资者的主要追求目标，面对不可避免的风险时，降低风险，防止或减少损失，以实现预期最佳是投资的目标。 当评价风险大小时，常有两种评价方式：定性分析与定量分析法。定性分析一般是根据风险度或风险大小等指标对风险因素进行优先级排序，为进一步分析或处理风险提供参考。这种方法适用于对比不同项目的风险程度，但这种方法最大的缺陷是在于，在多个项目中风险最小者也有可能亏损。而定量分析法则是将一些风险指标量化得到一系列的量化指标。通过这些简单易懂的指标，才能使公司的经营者、投资者对于项目分风险有正确的评估与判断，蒙特卡罗方法的应用【文献综述】
浅析蒙特卡洛方法原理及应用
蒙特卡洛方法
大学数学实验之蒙特卡洛方法
蒙特卡罗方法学习总结
蒙特卡罗方法(MC)
概率论实验报告蒙特卡洛方法估计积分值
蒙特卡罗方法及应用实验讲义2016资料
基于蒙特卡洛方法求数值积分与R
第19章-蒙特卡罗法与自助法
蒙特卡罗方法简介
用蒙特卡洛方法估计积分方法及matlab编程实现
用蒙特卡洛方法估计积分方法及matlab编程实现
蒙特卡罗方法及其在中子输运问题中得应用
matlab算法和蒙特卡罗计算教程
蒙特卡罗方法的解题过程可以归结为三个主要步骤
蒙特卡洛方法及其在风险评估中的应用
相关主题
文本预览