蒙特卡罗(Monte Carlo)模拟

MC模拟
假设有n个人在一起，各自的生日为365 天之一，根据概率理论，与很多人的直 觉相反，只需23个人便有大于50％的几 率人群中至少有2个人生日相同。
365 364 363 365 ( n 1) 1 ...... 365 365 365 365
n 理论几率
模拟几率
10
20 23 30 45 50
0.117
0.411 0.507 0.706 0.941 0.986
0.110
0.412 0.520 0.692 0.936 0.987
蒲丰的投针问题 Buffon’s Needle problem
Real function prob(n,npts) Integer i, npts Logical birth Real sum Sum=0 For I=1 to npts do If birth(n) then sum sum+1 Endfor Prob sum/real(npts) End function
[2] Zheng Zheng & Jordi Miralda-Escude. 2002, ApJ, 578, p33. [3]Watson, A.M. & Henney, W.J. 2001, RMxAA, 37, p221. [4]Lorne W. Avery & Lewis L. H. , 1968, ApJ, 152, p493. [5]Lawrence J. C., Peter D. N. & Jeffery D.S. 1972, ApJ, 176, 439. [6]David A. Neufeld, 1990, ApJ, 350, p216.
(averageof f over n randompointsin A)
例如： 求


sin ln(x y 1)dxdy
其中 {( x, y) : ( x 0.5)2 ( y 0.5)2 0.25}
在该区域中取5000个点，可得该积分约0.57.
0 x 1 0 y 1 0 z 1 x 2 sin y z
算法 （常微分）
蛙 跳 如Mori et al. 1998, ApJ, 494, 430 如Liu, W. J. et al. 2007, Adv. Space Res., in press
龙格库塔 CUEBIT
Hamilton et al. 2003, SP, 214, 339
保证了E=0时满足能量守恒。
17. 蒙特卡罗(Monte Carlo)模拟
蒙特卡罗方法，或称计算机随机模拟方法，是一种基于“随机数” 的计算方法。该方法源于美国在WWI后研制原子弹的“曼哈顿计 划”。 S. M. Ulam (1908-1984)和该计划的主持人之一、数学家 冯· 诺伊曼（J.von Neumann ）用驰名世界的赌城—摩纳哥的 Monte Carlo—来命名这种方法。 其基本思想很早以前就被人们所发现和利用。17世纪，人们就知道 用事件发生的“频率”来决定事件的“概率”。19世纪人们用投针 试验的方法来决定π。高速计算机的出现，使得用数学方法在计算机 上大量模拟这样的试验成为可能。
科技计算及社会中的问题比计算π要复杂得多。比如 金融衍生产品（期权、期货、掉期等）的定价及交 易风险估算，问题的维数（即变量的个数）可能高 达数百甚至数千。对这类问题，难度随维数的增加 呈指数增长，这就是所谓的“维数的灾难”(Curse of Dimensionality)，传统的数值方法难以对付（即使使 用速度最快的计算机）。Monte Carlo方法能很好地 用来对付维数的灾难，因为该方法的计算复杂性不 再依赖于维数。以前那些本来是无法计算的问题现 在也能够计算量。为提高方法的效率，科学家们提 出了许多所谓的“方差缩减”技巧。
1 b 1 n f ( x )dx f ( xi ) a ba n i 1
d b 1 1 n f ( x, y)dxdy f ( xi , yi ) c a (d c)(b a) n i 1
一般地，

A
f ( x )d ( measure of A)
蒲丰证明在相距d 的一系列平行线上 投长为 l 针。针与线的交点数除以投 的次数等于2 l / (πd)。 不妨取d ＝ l。产生随机数u和，其中 u为针中心距最近的线的距离，故u小 于等于0.5， 为夹角，根据对称性， sin u即为相交 可取0到π/2。 1 2
算时需要用上是否不自洽（姜冰）
作业
u

1 sin 2
中子屏蔽问题 Neutron Shielding problem 铅墙（长为5）
入 口
假设铅墙长为5，中子在铅中的平 均自由程为1，中子与铅原子碰撞 后各向同性散射。令碰撞8次后中 子能量耗尽，试求穿透铅墙的中 子的比例。 暂不考虑垂直纸面的运动，则中 子的水平位移是 1 cos cos ... cos 。
模拟结果
1
误差约
1 n
，它并不能和一些高级的数值积分算法比拟，
但对多维情况，MC方法却很有吸引力。

1 1 1
0 0 0
1 n f ( x, y, z )dxdydz f ( xi , yi , zi ) n i 1
我们可产生一系列随机数 1, 2 , 3 , 4 , 5 , 6 ,....... 可简单取3个随机数构成一个随机点，即 (1, 2 , 3 ), (4 , 5 , 6 ),....... 相应地，
辐射转移问题
对三维辐射转移问题,蒙特卡罗模拟几乎是必须的
计算方法
1. 将空间划分为Nx*Ny*Nz个网格 2. 选取空间中的随机点产生光子包(photon package), 并此向随机方向传播
3. 考虑光子与粒子的随机碰撞
产生具有某种分布的随机数
y f (x )
F ( x ) f ( x )dx
多随机点，当位于椭圆内的点数为1000时停止。
x
有时随机数产生函数通不过严格的随机性测试。而在一些计算（如MC 积分）中随机性非常重要。故使用更长字节的计算机更好。
应用之一
Numerical integration
估计面积和体积
选取(0, 1)中随机数序列x1， x2， x3， …… xn。则
1 n f ( xi ) 0 f ( x)dx n i 1
http://www.bjkp.gov.cn/gkjqy/rkx/rd2.htm
随机数 许多计算机系统都有随机数生成函数 F90: call random_seed call random_number(a) Matlab: x=rand(N) 产生元素在(0, 1)间随机分布的N*N矩阵 s= rand(‘state’,0) 重设该生成函数到初始状态 计算程序产生的随机数不是真正的随机数，它们是确定的，但看 上去是随机的，且能通过一些随机性的检验，故常称为伪随机数。
可由随机数序列 r (0,1)构造一系列随机数序列
如 (b-a)r+a integer((n+1)r) x (a,b) x {0,1,2,……n}
y
试产生1000个在椭圆x2+4y2=4内均匀分布的随机点：
1 y 1 中均匀产生足够 方法：在 2 x 2，
注意：上述随机数序列均具周期性，如上页random子程序的周期 约230。
Maple has a collection of random number generators. 如：
With(stats)
x :=uniform(0..1): seq(x(),i=1..10);
Matlab:
x=rand(10,1)产生10个随机数。
计算体积 Pseudo-code
x z ey 1
Program volume_region Integer parameter n 5000, iprt 1000 Real array (r ij )n3 integer i, m Real vol,x,y,z Call random((rij)) For I=1 to n do x ri1 y ri2 z ri3 If x 2 sin y z x z e y 1 then m m+1 If mod(I,iprt)=0 then Vol real(m)/real(I) Output I, vol Endif Endfor End program
xi ＝li/ (231-1) endfor
其它算法
1. 取初值x0 (0,1)，let xi be the fractional part of (π+ xi-1)5
2. Let u0 =1. For i=1,2,3,…, let ui be the remainder of (8t-3) ui-1 divided by 28, and xi = ui/2i. Here t can be any large integer
dF F 'dx 即分布密度 1 正比于 F '即f (x ) dx
试验粒子数值模拟
全粒子 研究粒子在电 磁场中的运动
d m 0 v q E v B dt dx v dt
研究微观不稳定性、 反常电阻等等
试验粒子
研究粒子加速等 忽略粒子间相互作用 及对电磁场的反馈
也可考虑库仑碰撞。
1 2 7
1
2
Hale Waihona Puke Baidu
21点问题 Blackjack 21 problem
参数拟合问题 Parameter fitting problem
[1] Ahn,S.-H., Lee, H.-W. & Lee, H.-M. 2000, JKAS, 33, p29.
辐射转移问题
Radiative transfer problem
计算体积之作业：冰淇淋的体积
z
x 2 y 2 ( z 1) 2 1 z2 x2 y2
y
x
应用之二
生日问题
Real function prob(n,npts) Integer i, npts Logical birth Real sum Sum=0 For I=1 to npts do If birth(n) then sum sum+1 Endfor Prob sum/real(npts) End function
另一类形式与Monte Carlo方法相似，但理论基础 不同的方法—“拟蒙特卡罗方法”(Quasi-Monte Carlo方法)—近年来也获得迅速发展。我国数学家 华罗庚、王元提出的“华—王”方法即是其中的 一例。这种方法的基本思想是“用确定性的超均 匀分布序列(数学上称为Low Discrepancy Sequences) 代替Monte Carlo方法中的随机数序列。对某些问 题该方法的实际速度一般可比Monte Carlo方法提 出高数百倍，并可计算精确度。
如 对32位字长的计算机
real procedure random((xi)) integer array (li)n real array (xi)n l0 any integer that 1< l0 <231-1 for i=1 to n do
li =(231-1)除以16807 li-1的余数