维普资讯 http://www.cqvip.com 物理与工程Ｖｏ１．１２ Ｎｏ．３ ２００２ ４５ 蒙特卡罗方法及应用① 尹增谦 管景峰？ 张晓宏 （ 曹春梅 华北电力大学物理教学部，河北保定 ０７１００３） （收稿日期：２００２ ０１ ０７） （？ 河北省保定市环境保护监测站，河北保定 ０７１０００） 摘 要 介绍蒙特卡罗方法的基本原理及其在计算物理中的应用． 关键词 蒙特卡罗方法 ＴＨＥ ＭｏＮＴＥ ＣＡＲＬｏ ＭＥＴＨｏＤ ＡＮＤ ＩＴＳ ＡＰＰＬＩＣＡＴＩｏＮ Ｙｉｎ Ｚｅｎｇｑｉａｎ （ Ｇｕａｎ Ｊｉｎｇｆｅｎｇ。 Ｚｈａｎｇ Ｘｉａｏｈｏｎｇ Ｃａｏ Ｃｈｕｎｍｅｉ Ｉ）ｅｐａｒｔｍｅｎｔ ｏｆ Ｐｈｙｓｉｃｓ，Ｎｏｒｔｈ Ｃｈｉｎａ Ｕｎｉｖｅｒｓｉｔｙ ｏｆ Ｅｌｅｃｔｒｉｃ Ｐｏｗｅｒ，Ｂａｏｄｉｎｇ，Ｈｅｂｅｉ ０７１００３） （ Ｅｎｖｉｒｏｎｍｅｎｔａｌ Ｄｅｔｅｃｔｉｎｇ Ｉｎｓｔｉｔｕｔｅ ｏｆ Ｂａｏｄｉｎｇ Ｃｉｔｙ，Ｂａｏｄｉｎｇ．Ｈｅｂｅｉ ０７１０００） Ａｂｓｔｒａｃｔ Ｔｈｅ ｐｒｉｎｃｉｐｌｅ ｏｆ Ｍｏｎｔｅ Ｃａｒｌｏ ｍｅｔｈｏｄ ａｎｄ ｉｔｓ ａｐｐｌｉｃａｔｉｏｎ ａｒｅ ｉｎｔｒｏｄｕｃｅｄ． Ｋｅｙ Ｗｏｒｄｓ Ｍｏｎｔｅ Ｃａｒｌｏ ｍｅｔｈｏｄ 可以追溯到１８世纪著名的蒲丰问题．１７７７ ｌ 引言 年，法国科学家蒲丰（Ｂｕｆｆｏｎ）提出用投针试验 计算圆周率丌值的问题．这里我们用蒲丰问 蒙特卡罗方法，又称随机抽样方法，是一 种与一般数值计算方法有本质区别的计算方 法，属于试验数学的一个分支，起源于早期的 用几率近似概率的数学思想，它利用随机数 进行统计试验，以求得的统计特征值（如均 值、概率等）作为待解问题的数值解．随着现 代计算机技术的飞速发展，蒙特卡罗方法已 经在原子弹工程的科学研究中发挥了极其重 要的作用 ］，并正在日益广泛地应用于物理 工程的各个方面，如气体放电中的粒子输运 过程等 ］．本文介绍蒙特卡罗方法的基本原 理及其在计算物理中的应用． ２ 蒙特卡罗方法的基本原理 就数学特性而言，蒙特卡罗方法的发展 题来初步说明蒙特卡罗方法的基本原理和解 决问题的基本手续． 蒲丰问题是这样一个古典概率问题：在 平面上有彼此相距为２“的平行线，向此平面 任意投一长度为２ 的针，假定 ＜“，显然，所 投的针至多可与一条直线相交，那么，此针与 任意条平行线相交的概率可以求出，由下面 的分析可知，此概率与所取针长２ 、平行线间 距２“有关，并且包含有丌值．在这里，任投一 针的概率含义有以下三点：（１）针的中点Ｍｌ 在平行线之间等概率落入，即Ｍｌ距平行线的 距离３２均匀分布在区间Ｅｏ，“］之内；（２）针与 ｒ １ 线的夹角 均匀分布在区间Ｉ 号，号ｌ之内； Ｌ 厶 厶＿Ｊ （３）．２Ｃ与 互相． 如图１所示，建立与平行线垂直且原点 ① 本文工作获华北电力大学“大学物理教学研究”教改项Ｉ；１基金资助 维普资讯 http://www.cqvip.com ４６ 物理与工程Ｖｏ１．１２ Ｎｏ．３ ２００２ 然，这是一个均匀分布，即针中心处于区间 （ｚ，ｚ＋ｄｘ）内的概率为 ｄＰ ．一一ｄｘ （２） ／， 这样，一次投掷，针中心落入（ｚ， ’＋ｄｘ）且与 ｌ ｌ线相交的概率为 ｄＰ—Ｐ１ｄＰ２一二ａｒｃｃｏｓ÷ｄｌ’ （３） 则一次投掷，针与线相交的总概率为 １ 蒲 题的概率分析 在某一条平行线上的 轴，不失一般性，假定 Ｐ—ｆＪ ｄＰ—ｆＪ ａ０ ７ｃ“ ｒｃｃ。ｓ车ｄｌ￡ ｒ．７ｃ“ ２ｌ（４） 针的中心处于图示中的ｚ轴上．由于对称性， 即： 丌一 ‘５） 我们只需分析针中心处在ｚ∈（０，“）范围的情 况即可．令探针中心的坐标值为 ，显然，只有 从（５）式可见，可利用投针试验计算丌值：设 ≤ｚ时才可能发生相交的事件．我们来分析 投针Ｎ次，其中 次针与线相交，则可用频率 在条件 ≤ｚ满足时，针与线相交的概率：只 值 ／Ｎ作为概率Ｐ的估计值，从而求得 的 估计值为 有当 ≤ 一ａｒｃｃｏｓ｝时才能相交，且相交的 丌≈型 （６） 概率为 以 Ｐ。一 ２ 这就是早期的用频率值作为概率近似值的方 ７ｃ ａｒｃｃ。ｓ手 ￡ （１） 法的应用实例，表ｌ是在历史上一些有名的 下面再来分析针中心位置在轴上的分布，显 用投针试验计算丌值的结果 ］，其中针长以“ 为单位． 表１投针试验计算 值的结果 实验者 时间（年份） 针长 投针次数 相交次数 的估值 ｗＯｌｆ ｌ８５０ ０．８ ５，０００ ２， ３２ ３．１ ５９６ Ｓｍｌｔｈ １ ８ｊ ０．６ ３．２Ｏ１ １、２１ ８ ｊ ３．１ ５ｊ１ Ｄｅ Ｍｏｒｇａｎ．Ｃ １ ８６０ １ ０ ６００ ３８２．ｊ ３．１ ３７ ＦｏＸ ｌ ８８ｄ ０．７５ １，０３０ ｌ８９ ３ ｌ ５９５ ｌ ９Ｏ１ ０．８３ ３，１０８ １、８０８ ３¨１ ５９２９ Ｒｅｉｎａ １９２ｊ ０ ５４１ ９ ２．ｊ２Ｏ ８５９ ３．１ ７９５ 需要指出的是，上述由投针试验求得丌 子的输运过程及粒子输运的总效应，若要用 的近似值的方法，是进行真正的试验，并统计 多次掷骰子的方法近似求出就是不可能的 试验结果，要使获得的频率值与概率值偏差 了．所以，在现代计算机技术出现之前，用频 小，就要进行大量的试验，这在实际中，往往 率近似概率的方法一抑或称为雏形时代的 难以做到．可以设想，对蒲丰问题这样一个简 蒙特卡罗方法——并没有得到实质上的 单的概率问题，若要进行１０万次投针试验， 应用． 以每次投针、作出是否相交判断并累加相交 若用数值模拟方法代替上述的真正的投 次数用时５秒钟计算，则需用时５Ｏ万秒，即 针试验，是利用均匀分布于（０，１）之间的随机 大约１３９个小时．那么，可以设想，对于象上 数序列，并构造出随机投针的数学模型，然后 述确定条件下的核裂变、直流气体放电中粒 进行大量的随机统计并求得丌的近似值． 维普资讯 http://www.cqvip.com 物理与工程Ｖｏ１．１２ Ｎｏ．３ ２００２ ｌ Ｉ ▲ ≈ Ｉ ｊ 『 图２ 用数值模拟方法计算蒲丰问题 如图２建立坐标系，平面上一根针的位 置可以用针中心Ｍｔ的坐标 和针与平行线 的夹角 来决定，在ｙ方向上的位置不影响 相交性质．任意投针，意味着 与 都是任意 取的．但 的范围可限于［０， ］， 的范围可 限于［０，ａ］．在这种情况下，针与平行线相交 的数学条件是 ≤ｌｓｉｎ０，０≤ ≤ｎ （７） 其次，怎样模拟投针呢？亦即如何产生任意 的［ ， ］． 在［０，ａ］任意取值，意味着 在 ［０，“］上取哪一点的概率都一样，即ｚ的概率 密度函数为 ｆ （ ）一 ）一』一｝—ａ１’当０—’ ’ ＼≤ｚ≤ｄ时（＼Ｉ、 、上＼＼“Ｕ 。 （）８） ｌ ０， 当 为其他值时 类似的， 的概率密度函数为 （ ）一ｊ １， 当０≤ ≤７【时 （９） １ ０， 当 为其他值时 由此，产生任意（ｚ， ）的过程就变为由ｆ （ｚ） 抽样 ，由ｆ。（ ）抽样 的过程．容易得到 ３－＂一ｎ １ ∞ ： 式中， ， 均为（０，１）上均匀分布的随机数． 只要随机数的均匀性和性良好，如此构 造的数值模型就很好地模拟了实际试验中的 一次投针，并用下式判断是否相交且记录统 计结果： ｓ（３７，，Ｏｉ）一 如果投针Ｎ次，那么 一 ∑ （７３ ， ） （１１） ４７ 是相交几率Ｐ的估计值．这样就实现了用数 值方法模拟真正的投针试验．用此方法计算 的 的近似值的情况如表２所示． 表２用蒙特卡罗方法计算的 的近似值 投针次数 １０．０００ ２０．００（） １００．０００ ２００．０００ 的近似值 ３ １ ６２２３３ ３ １ ３７９９３ ３ １１１１ ７９ ３．１ ｒ１１ ３５１ 表２中的计算结果表明，随着模拟投针 次数的增大。所计算的 的近似值越来越接 近于其真值．而要进行这样的数值模拟，就需 要很大的计算量，只有利用计算机才能实现． 从蒲丰问题可以看出，用蒙特卡罗方法 求解问题时，应建立一个概率模型，使待解问 题与此概率模型相联系，然后通过随机试验 求得某些统计特征值作为待解问题的近似 解．与此相似，在一些物理问题，如核裂变、直 流气体放电等过程中，粒子的输运过程及粒 子输运的总效应，也是可以与某些概率过程 联系起来，例如，电子与原子、分子、离子的碰 撞过程，实际上就是与碰撞截面有关的概率 过程，这样，从数学物理特征来说，类似于用 随机投针方法计算 的近似值，确定条件下 的核裂变、直流气体放电中粒子的输运过程 及粒子输运的总效应可以用多次掷骰子的方 法近似求出． 随着现代计算机技术的出现和飞速发 展，用计算机模拟概率过程，实现多次模拟试 验并统计计算结果，进而可获得所求问题的 近似结果．计算机的大存储量、高运算速度使 得在短时间内，获得精度极高且内容丰富的 模拟结果．在历史上，也正是原子弹工程研究 初期阶段的工作，为模拟裂变物质的中子随 机扩散，提出了运用大存储量、高运算速度计 算机的要求，这也成为当时推动计算机技术 发展的重要动力．也就是在第二次世界大战 期间，冯・诺依曼和乌拉姆两人把他们所从 事的与研制原子弹有关的秘密工作一 对裂 变物质的中子随机扩散进行直接模拟——以 摩纳哥国的世界闻名蒙特卡罗（Ｍｏｎｔｅ Ｃａｒｌｏ）作为秘密代号来称呼．用名比喻 随机模拟，风趣又贴切，很快得到广泛接受， 维普资讯 http://www.cqvip.com ４８ 此后，人们便把这种计算机随机模拟方法称为 蒙特卡罗方法．需要指出的是，正是由于广泛 领域的物理问题中存在着大量的随机过程，如 粒子间的碰撞等，使得蒙特卡罗方法在计算机 物理和物理工程中得到Ｅｔ益广泛的应用，并成 为沟通理论与实验研究的一个桥梁． 需要指出的是，蒙特卡罗方法不仅在处 理具有概率性质的问题方面获得广泛的应用， 对于具有确定性问题的计算也因其程序简单 等优点获得了广泛的应用．这里以定积分的计 算简要说明其处理确定性问题的手续． 对于定积分 ｓ一 池 通过变量替换，可以转换为下面的形式 ＝＝＝ｋ ｌ ｇ（Ｌｚ）ｄｘ ‘ 其中ｇ（Ｌｚ）∈（０，１），当Ｌｚ∈（０，１）时 （１２） ｝ 即转换为求积分Ｉ ｇ（ｘ）ｄ：ｃ亦即求边长为１的正 方形中一个曲边梯形的面积的问题，如图３所示． 图３用蒙特卡罗方法求定积分 我们可以设想这样一种随机投点求定积 分Ｉ ｇ（ｚ）ｄｘ的方法：在一个边长为１的正方形 ｔ 上并以其两边分别为坐标轴画出曲线ｇ（ ｚ）， 实际上就是图３，然后随机地正方形投掷小 球，那么，小球击中ｇ（ ｚ）曲线下部分的概率 ｝ 就等于所要求的积分Ｉｇ（ ｚ）ｄｘ，这样就将确 物理与工程Ｖｏ１．１２ Ｎｏ．３ ２００２ 定性的定积分问题转化为一个概率问题，同 样可以通过数值模拟方法——蒙特卡罗方 法求得其近似解．用此方法，我们计算了积分 ｊ’ ｒ ｄｘ，当投球数为１万次时，得到的积分近 似值为０・３３２８００，与其真值专极为接近． ３ 蒙特卡罗方法的解题手续和特点 在用蒙特卡罗方法解算问题时，一般需 要这样几个过程：构造或描述概率过程，对 于本身就具有随机性质的问题，如粒子输运 问题，主要是正确地描述和模拟这个概率过 程．对于本来不是随机性质的确定性问题， 比如计算定积分、解线性方程组、偏微分方 程边值问题等，要用求蒙特卡罗方法求解， 就必须事先构造一个人为的概率过程，它的 某些参量正好是所要求问题的解．即要将不 具有随机性质的问题，转化为随机性质的问 题．这构成了蒙特卡罗方法研究与应用上的 重要问题之一．然后建立各种估计量，使其 期望值是所要求解问题的解．最后根据所构 造的概率模型编制计算程序并进行计算，获 得计算结果． 与其他的数值计算方法相比，蒙特卡罗 方法有这样几个优点：（１）收敛速度与问题 维数无关，换句话说，要达到同一精度，用蒙 特卡罗方法选取的点数与维数无关；计算时 间仅与维数成比例．但一般数值方法，比如在 计算多重积分时，达到同样的误差，点数与维 数的幂次成比例，即计算量要随维数的幂次 方而增加．这一特性，决定了对问题的适 用性．（２）受问题的条件的影响小．（３） 程序结构简单，在电子计算机上实现蒙特卡 罗计算时，程序结构清晰简单，便于编制和 调试．（４）对于模拟象粒子输运等物理问题 具有其他数值计算方法不能替代的作用．蒙 特卡罗方法的弱点是收敛速度慢，误差大的 概率性质．这一情况在解粒子输运问题中仍 然存在．除此之外，经验证明，只有当系统的 维普资讯 http://www.cqvip.com 物理与工程Ｖｏ１．１２ Ｎｏ．３ ２００２ 中的应用．北京：科学出版社，１９８０． ４９ 大小与粒子的平均自由程可以相比较时，一 般在１０个平均自由程左右，这方法算出的结 果较为满意．而对于大系统深穿透问题，算出 －１２３ 方再根．计算机模拟和蒙特卡罗方．北京：北京工 业学院出版社，１９８８． ［３］ Ｂ．Ｅｌｉａｓｓｏｎ ｅｔ ａ１．Ｍｏｎｔｅ Ｃａｒｌｏ Ｓｉｍｕｌａｔｉｏｎ ｏｆ Ｒｕｎａｗａｙ 的结果往往偏低．对于大系统，其他数值方法 往往很适应，能算出较好的结果．因此，已有人 Ｅｌｅｃｔｒｏｎｓ ｉｎ ０２／Ｎ２ Ｍｉｘｔｕｒｅ．Ｈｅｒｂｓｔｔａｇｕｎｇ ｄｅｒ ＳＰＧ／ＳＳＰ．１９８７，Ｖｏ１．６０ ＰＰ．２４１～２４７． 将数值方法与蒙特卡罗方法联合起来使用，克 服这种局限性，取得了一定的效果． 随着计算机技术的飞速发展，蒙特卡罗方 法以其独特的优点广泛应用于计算物理和物 ［４］ 陈俊英．Ｈ ｚ／ＣＨ 系统ＥＡＣＶＤ动力学过程研究．河 北大学硕士学位论文，２０００． １－５３ 张玉红．Ｈ ｚ／Ｃ２Ｈ２系统电子助进热丝化学气相沉积 动力学过程研究．河北大学硕士学位论文，２００１． 理工程领域，对蒙特卡罗方法的推广必将使其 对物理学学科的发展发挥更大的作用． 参 考 文 献 作者简介 尹增谦，男。河北省定州市人，１９９１年毕业于哈尔滨工 业大学并获工学硕士学位，现任华北电力大学物理教学部 副教授，河北大学光学工程专业在职博士研究生，主要从事 光学工程领域及计算物理领域的科学研究工作和物理教学 工作，已发表学术论文二十余篇． ［１］裴鹿成，张孝泽．蒙特卡罗方法及其在粒子输运问题 （上接第４４页） 常用热管的工作温度范围与典型的工作介质及其相容壳体材料表 种类 低 温 热 工作介质 氮 氟里昂一２１（ＣＨＣＩ２Ｆ） 氟里昂一１１（ＣＣＩ３Ｆ） 工作温度／℃ 一６０～１００ ４０～１００ ４０～１ ２０ 相容壳体材料 铝、不锈钢、低碳钢 铝、铁 铝、不锈钢、铜 管 氟里昂一１１３（ＣＣ１２Ｆ・ＣＣＩＦ２） 己烷 —１０～１００ ０～１００ ０～１ ２０ ０～１３０ １０～１３Ｏ ０～２９０ ３０～２５０ １４ ７～３５０ １ ４７～３００ ｌ５０～３９５ １４ ７～３００ ２５０～６５０ ４００～１．０００ ４００～１１００ ５００１．２００ ．～铝、铜 黄铜、不锈钢 铝、铜、不锈钢 铜、不锈钢 铜、不锈钢、碳钢 不锈钢、常 温 热 管 丙酮 乙醇 甲醇 甲苯 水 萘 低碳钢、低合金钢 铜、碳钢（内壁经化学处理） 铝、不锈钢、碳钢 不锈钢碳钢 、中 温 热 管 联苯 导热姆Ａ 铜、不锈钢、碳钢 不锈钢、碳钢、镲 奥氏体不锈钢 不锈钢 钛、铌 导热姆一Ｅ 隶 钾 高 温 热 管 铯 钠 锂 银 不锈钢、因康镍台金 钨、钽、钼、铌 钨、钽 １．０００～１．８Ｏ０ １．８Ｏ０～２．３Ｏ０ 社，１９８８． 参 考 文 献 ［２３胡亚范等．热管技术应用于钻井余热回收的方案设 计．应用能源技术．２０００．５． ［１］ 商政宋等译．实用热管技术．北京：化学工业出版
