计算物理 蒙特卡罗方法 http://125.217.162.13/lesson/ComputationalPhysics

蒙特卡罗方法 蒲丰投针 收敛性、误差和优缺点 任意分布的随机数 粒子输运问题 随机过程模拟 梅氏抽样 √

蒲丰投针(1/5) 蒙特卡罗方法 蒲丰投针 又称随机抽样技巧或统计试验方法 以概率统计理论为基础的 能够比较逼真地描述事物的特点及物理实验过程 解决一般数值方法难以解决的问题 随着电子计算机的发展而发展 首先在核武器的试验与研制中得到了应用 蒲丰投针 法国数学家蒲丰的1777年出版的著作：“在平面上画有一组间距为 d 的平行线，将一根长度为 l (l<d) 的针任意掷在这个平面上，求此针与平行线中任一条相交的概率。” √

蒲丰投针(2/5) 步骤 在桌面上画出间距为 2d 的平行线 准备长度为 2l (l<d) 的针 向桌面随机投针 如果针与平行线相交，则计数器 n 加 1 计算：计数器 n 与总投针数 N 的比例(视作相交概率 P ) 概率分析 P = ? 各条平行线地位等同，仅考虑某条平行线附近的情况 平行线方向的 x 坐标对概率没影响 针的中点的 y 坐标在线之间等概率落入(均匀分布在 [0, d])，仅当 yl 才可能针-线相交 针-线的夹角 q 均匀分布在 [0, p]，q 与 y 独立 y x √

蒲丰投针(3/5) 概率 P = 2l / (p d) ，可计算圆周率 x y 实验者 时间 针长 总投数 相交数 p 值 Wolf 1850 0.8 5,000 2,532 3.1596 Smith 1855 0.6 3,204 1,218 3.1554 De Morgan,C 1860 1.0 600 382 3.137 Fox 1884 0.75 1,030 489 3.1595 Lezzerini 1901 0.83 3,408 1,808 3.1415929 Reina 1925 0.5419 2,520 859 3.1795 √

蒲丰投针(4/5) 关于蒙特卡罗方法的分析和总结 基本思想 确定所求问题的解是某事件的概率(或某随机变量的数学期望、或与概率/数学期望有关的量，如 p ) 通过试验方法，得出事件发生的频率(或该随机变量若干个具体观察值的算术平均值，如 P )，求解 数学期望与概率：当随机变量的取值仅为 1 或 0 时，它的数学期望就是事件的概率；反之亦然 数学期望与算术平均值 用随机试验的方法计算积分，将积分看作服从分布密度函数 f (r) 的随机变量 g (r) 的数学期望 通过试验，得到 N 个观察值 r1, r2, …, rN (从 f (r) 中抽取 N 个子样 r1, r2, …, rN )，求 g (r) 的算术平均值 √

蒲丰投针(5/5) 试验方法和次数 试验方法不一定可行 精确的近似解需要巨量的试验次数，但人工方法作大量的试验相当困难，甚至是不可能的 用计算机模拟随机试验过程，完成巨量的次数 抽象 x 的分布密度函数： q 的分布密度函数： 产生任意的 (x,q ) = 由 f1(x)抽样 x + 由 f2(q)抽样 q 对应的随机变量： 数学期望与算术平均值 新的问题：误差？收敛？ √

收敛性、误差和优缺点(1/4) 收敛性 求解：以随机变量 X 的简单子样 X1, X2, …, XN 的算术平均值，作为求解的近似值 近似值的收敛性 大数定理：当试验次数足够多时，事件发生的频率无穷接近于该事件发生的概率 如果 X1, …, XN 独立分布，且期望值有限( E(X)< )，那么 当随机变量 X 的简单子样数 N 充分大时，其算术平均值以概率 1 收敛于期望值 E(X) √

收敛性、误差和优缺点(2/4) 误差 概率论的中心极限定理：如果随机变量序列 X1, X2, …, XN 独立分布，且具有有限非零的方差 s 2，即 其中的 f(x) 是 X 的分布密度函数，则 当 N 足够大时(蒙特卡罗方法) 不等式 的概率约 1-a 误差定义为 ，收敛速度为 a 0.500 0.050 0.003 la 0.674 1.960 2.968 √

收敛性、误差和优缺点(3/4) 蒙特卡罗方法的误差为概率误差 均方差 s 是未知的，估计值为 减少误差的技巧(在确定的置信度 a 前提下) 误差 e 与试验次数的开根号 N1/2 成反比：精度一个数量级，次数 N 两个数量级——巨大的代价 误差 e 与均方差 s 成正比：精度一个数量级，均方差 s 一个数量级——可接受的代价 效率 降低方差增加观察子样的时间固定时间内样本数减少代价增加 蒙特卡罗方法中的效率：由均方差 s 和观察一个子样的费用(计算机时) c 衡量 = s 2 c √

收敛性、误差和优缺点(4/4) 优点 缺点 主要应用范围 较逼真地描述具有随机性质的事物的特点及物理实验过程 受几何条件限制小 收敛速度与问题的维数无关 具有同时计算多个方案与多个未知量的能力 误差容易确定 程序结构简单，易于实现 缺点 收敛速度慢 误差具有概率性 在某些问题中，计算结果与系统大小有关 主要应用范围 粒子输运，统计物理，典型数学，真空技术，激光技术，医学，生物，探矿等 √

任意分布的随机数(1/11) 随机数(蒙特卡罗方法的关键) 直接抽样方法 基本的：均匀地在 (0, 1) 内分布 任意的：非均匀地在 (a, b) 内分布 方法：直接的(离散、连续)、舍选的(简单、乘、乘加) 直接抽样方法 离散型：产生随机量 x 的抽样值 xi ，概率为 Pi (i = 1, 2, …) 方法：先计算 x 的分布函数 对于产生的 R，如果满足 Fi-1< R  Fi，则令抽样 x = xi 证明： √

任意分布的随机数(2/11) 例：二项分布 直接抽样方法： 产生 R，如果满足 ，则令 x = n 例：泊松分布 直接抽样方法： √

任意分布的随机数(3/11) 连续型：产生随机量 x 的抽样值，概率密度函数为 f (x) 方法：先计算 x 的分布函数 对于随机数 R'，解方程 R' = F(x')，得到 x' = F-1(R')，则令 x = x' 证明： √

任意分布的随机数(4/11) 例：均匀介质中，粒子运动的自由程 S 是随机量，其概率密度函数为 f (S) = Se -SS；其中的 S = ns 是宏观总截面，s 是原子截面，n 为介质中的原子数密度 分布函数： 介质中粒子的碰撞过程  随机量 R 中的抽样过程 例：散射方位角余弦分布 √

任意分布的随机数(5/11) 舍选抽样方法 直接法的特点：简单，但分布函数的反函数不一定有解析解 简单分布 定理：如果 Z 是 [a, b] 上均匀分布的随机数，那么利用条件 f(Z)/MR (M 是 f(Z) 的上界)选出的 Z 将是 [a, b] 上概率密度函数为 f(Z) 的随机数 证明：如右图所示 R Z f (Z) / M a b 1 Z 为横坐标，R 为纵坐标，实曲线为函数 f(Z)/M R 在 (0, 1) 内、Z 在 [a, b] 内均匀分布  随机点 (Z', R) 在虚框内均匀分布 随机点落入窄条的概率 = 两面积之比 舍 dZ 选 (Z', R) 已知概率密度函数抽样方法 √

任意分布的随机数(6/11) 简单分布抽样方法的流程 效率(有一部分随机数被舍弃) 例：受限的散射方位角余弦分布 效率 产生随机数 x 和 R 计算 Z = a + (b-a) x f (Z) / M  R ? Z' = Z Y N 舍 dZ (Z', R) R Z f (Z) / M a b 1 选 效率(有一部分随机数被舍弃) 例：受限的散射方位角余弦分布 效率 √

任意分布的随机数(7/11) 乘分布：f(x) 有锐锋时效率很低，需要改进的简单分布 方法：将函数写成 f(Z) = h(Z) f1(Z)，由容易抽样的 f1(Z) 抽样出 Z，代入 h(Z)，如果满足条件 h(Z)/MR (M 是 h(Z) 的上界)，那么得到概率密度函数为 f(Z) 的抽样值 证明：如右图所示 R F1(Z) h (Z) / M 1 f1(Z) 的分布函数 F1(Z) 为横坐标，R 为纵坐标，实曲线为函数 h(Z)/M R 和 F1 在 (0, 1) 内均匀分布  随机点 (F1', R) 在虚框内均匀分布 随机点落入窄条的概率 = 两面积之比 舍 dF1 选 (F1', R) 效率 √

任意分布的随机数(8/11) 乘分布抽样方法的流程 例：半正态分布，概率密度函数 产生随机数 F1 和 R 利用 F1， 由 f1(Z) 抽样 Z h (Z) / M  R ? Z' = Z Y N 例：半正态分布，概率密度函数 √

任意分布的随机数(9/11) 乘加分布 方法：如果随机量 x 具有以下乘加形式 其中 f1(x) 和 f2(x) 满足密度函数的要求，即 并且 h1(x) 和 h2(x) 满足 那么从 f(x) 的以下变形，可以得到乘加分布 其中 和 为加权因子 证明：参照直接法和舍选法(略) √

任意分布的随机数(10/11) 乘加分布抽样方法的流程 效率 产生随机数 g, F, R 利用 F 由 f1(x) 抽样 x1 g  g1 ? Y N h1(x1)/M1  R ? 由 f2(x) 抽样 x2 h2(x2)/M2  R ? x = x1 x = x2 效率 √

任意分布的随机数(11/11) 例：散射光子能量抽样。能量为 E0 的入射光子，经原子散射后的能量 E 按某种概率分布。如果令 x = E0/E，那么概率密度函数 f(x) 为(其中 K = K(E0) 为归一因子) √

粒子输运问题(1/7) 蒙特卡罗模拟 粒子输运是随机过程，运动规律是大量粒子运动的统计 蒙特卡罗模拟：模拟一定数量粒子在介质中的运动，再现粒子运动的统计规律 H 例(平板介质模型)：由单一物质组成的均匀介质，厚度为 H，能量为 E00 的平行光子束垂直射入板内，求光子对板的投射率 抽样：自由程、作用类型、散射能量、散射方向 自由程抽样(PPT15)：粒子运动的自由程 S 是随机量，其概率密度函数为 f (S) = Se -SS；其中的 S = ns 是宏观总截面，s 是原子截面，n 为介质中的原子数密度 分布函数：F(S') = 1 - e-SS' 粒子的碰撞过程  随机量 R 中的抽样过程：S' = -lnR' / S √

粒子输运问题(2/7) 作用类型抽样 光子-介质：散射(康普顿散射)和吸收(光电效应) 散射截面 吸收截面 总截面：St = Ss + Sa 光子随机地被散射或吸收的过程  随机量 R 中的抽样过程：如果 Ss / St  R，则散射，否则为吸收 √

粒子输运问题(3/7) 散射能量抽样(PPT22)：能量为 E0 的入射光子，经原子散射后的能量 E 按某种概率分布。如果令 x = E0/E，那么概率密度函数 f(x) 为(其中 K = K(E0) 为归一因子) √

粒子输运问题(4/7) x y z 散射方向抽样 坐标系：以入射方向 (由 q 和 f 确定)为参考系，散射方向 由散射角 q' 和散射方位角 f' 确定 流程图 输入 g0, E0, E 计算 m =1+1/E0-1/E 计算 g 产生随机数 R 计算 f' = 2pR √

粒子输运问题(5/7) 直接模拟方法 模拟粒子在介质中的真实物理过程 粒子在介质中的状态：空间位置，能量和运动方向 碰撞点的状态参数 (Zm, Em, gm) 表示从源出发的粒子在介质中经过 m 次碰撞后的状态 Z H O Z2 E2 q2 Z0 E0 q0 Z1 E1 q1 粒子的运动过程 = 碰撞点的状态序列 模拟粒子的运动过程 = 确定状态序列的问题 √

粒子输运问题(6/7) 记录和计算的内容 穿透率和误差估计(设 N 为入射粒子数，N1 为透射数) 透射粒子的能量和方向分布 离散化能量：Emin = E0 < E1 <  < Ei <  < EI 离散化角度：0 = q0 < q1 <  < qj <  < qJ = p / 2 能量和方向分布 粒子的轨迹 粒子运动关于 Z 轴对称，只需要记录 (Zim, qim) 流程图 √

粒子输运问题(7/7) 输入N, H, E00, Emin, g00 i>N ? N Y m>M ? N Y Z  0 ? 类型=吸收 ? E  Emin ? Y N Z  H ? N Y 光子 i 的初态m=0, Z=0, E=E00, g=g00 利用 Ss , St 抽样 作用类型 利用 E 计算截面 Ss , St mm+1 利用 Ss 抽样 自由程 S 利用 E0 E 抽样 散射能量 E 计算碰撞点位置 ZZ+Sg ii+1 利用 g0g, E0 , E 抽样散射方向 g 计数器 N1N1+1 输出 P=N1/N, s, DP √

随机过程模拟(1/3) 蒙特卡罗模拟随机过程的两种情况 已知随机过程的概率分布函数随机过程的统计特征 建立与随机过程的模型，形成随机量，并使其数字特征(概率、平均值、方差等)是问题的解 由已知的概率分布函数进行大量的抽样 统计处理抽样结果，给出问题的解及其误差 已知随机现象的观测数据概率分布函数 分析现象的特征，假设随机量服从某种分布，建立概率模型 由观测数据推断分布中的参数 按推断的分布进行大量的抽样 比较抽样值和观测值，根据误差，判断假设的准确性 √

随机过程模拟(2/3) 例：a 粒子衰变的蒙特卡罗模拟 随机现象的观测数据概率分布函数 k nk 57 6 273 1 203 7 139 2 383 8 45 3 525 9 27 4 532 10 16 5 408 例：a 粒子衰变的蒙特卡罗模拟 随机现象的观测数据概率分布函数 观测数据：每隔 Dt<<半衰期 测量一次放射的 a 粒子数，共测量 N=2608 次，测得 k 个粒子的次数为 nk 分析现象，建立模型 Dt 内的衰变数 k 是随机事件的发生次数，是随机量 k 在有限的平均值 a 上下波动 k 可看作大量( N 足够大)独立(每个原子核的衰变不受其它原子核的影响)试验的结果 事件有相同的很小的概率(原子核衰变的机会相等) k 的上述特征表明： k 近似服从泊松分布 √

随机过程模拟(3/3) 蒙特卡罗模拟 产生 2608 个泊松分布的随机数，统计其中数值等于 k 的个数并与观测值 nk 比较 √

梅氏抽样(1/1) 特别适合于多维随机量的系统 设 x 是多维随机量，f(x) 是概率密度函数，流程 梅氏1步 输入 x, dx 产生多维随机数R x'=x+dx(R-0.5) 产生随机数r f(x')/f(x)r ? x=x' 返回 x Y N 输入 Nt, Ng, Nf, dx 初始化 x 由 x 走梅氏 Nt 步 间隔 Nf 抽样： 由 x 走梅氏 Nf 步 由 x 计算 u 和 u2 AA+u, BB+u2 m>Ng ? Y N mm+1 计算平均值, 方差, 误差 主程序 结束 √

作业 用蒲丰投针在计算机上计算 p 值，取 d = 4, l = 3 分子热运动的速率概率密度函数是 取 f1(x) = 2e-2x/3/3，找出 h(x) 和 M ，用直接法对 f1(x) 抽样并写出程序