使用蒙特卡罗方法和粒子追踪估算圆周率
使用蒙特卡罗方法和粒子追踪估算圆周率
圆周率π是数学和科学领域中最具统治力的数字之一，其历史可以追溯到数千年以前。从古代文明对圆的粗略近似，到现代超级计算机计算万亿位数，π一直吸引着数学家和好奇者的想象力。本文将通过COMSOL Multiphysics®仿真软件提供的功能，以一种有趣和流行的方法来计算圆周率。
历史上对圆周率的近似计算
已知最早的圆周率近似值出现在古代文明中。巴比伦的数学家将圆周率近似为 3，这个数值在当时的建筑工程中是合理的，后来又被精确为 3.1251。埃及的数学家和印度的学者分别通过比较圆形和八边形的面积2，和通过巨量的计算3得出了近似的数值。包括阿基米德在内的一些希腊学者利用几何方法将圆周率的近似值精确到 3 个数量级以内4，使圆周率的计算取得了重大突破。
左图：Domenico Fetti 于 1620 年创作的 Archimedes Thoughtful（又名：Portrait of a Scholar）。图片属于公有领域，通过Wikimedia Commons共享。右图：Leonardo Fibonacci 的肖像。图片属于公有领域，通过Wikimedia Commons共享。
如今我们普遍使用的 3.14 近似值来自中国数学家刘徽，他提出这个近似值的目的是为了实用5。圆周率计算的后续发展涉及无穷级数估算和三角关系的利用。数学家们利用微积分推导出了无穷级数，可用于计算高精度的圆周率，其中Fibonacci和 Al-Khwarizmi 做出了重要贡献。
这些发展为我们使用的现代方法奠定了基础，包括计算机中使用的算法，即通过先进的数学工具和计算能力计算万亿位数的算法，以极高的精度计算圆周率。一些著名的计算方法包括 Chudnovsky 算法、Gauss–Legendre 算法、Machin 公式，以及 Monte Carlo（蒙特卡罗）方法。
通俗易懂的蒙特卡罗方法
蒙特卡罗方法是一种依靠随机抽样来估计数值结果的计算技术，特别适用于包含大量变量的问题。对于这种情况，可以利用内在的随机性来解决确定性问题。想象这样一个场景，你正在为一场聚会计算需要订购多少个披萨。这里的确定性问题是计算每个人要吃多少片披萨。与其询问每个人要吃多少片披萨，然后求和得出结果（这对一个大型聚会来说可能相当麻烦），不如随机挑选几个朋友，询问他们要吃多少片披萨，然后求平均值来解决问题。这有点像蒙特卡罗方法，即使用随机样本来估计一个值。蒙特卡洛法被广泛用于模拟复杂现象，如流体、统计力学、生物化学、密码学、社会学和心理学。
这种思维可以扩展到现在流行的一种有趣的估计圆周率的方法。这种方法是在一个正方形内随机放置一些点，然后计算有多少点位于正方形内切圆内。圆内的点数与总点数之比可以用来近似计算圆周率。由于内嵌在边长为 2r 的正方形中的圆的面积为 πr²，而正方形的面积为 (2r)²= 4r²，因此它们的面积之比为 π/4。也就是说一个点落在圆内的概率是 π/4。因此，如果我们将圆内点数与总点数之比乘以 4，就可以得到 π 的估计值。这是因为随着点数的增加，比率会趋近于实际值 π/4。
估计圆周率的基础。
在 COMSOL Multiphysics®中使用蒙特卡罗方法估算圆周率
为了进行这个简单的蒙特卡罗模拟，我们将使用数学粒子追踪接口。在 COMSOL Multiphysics®软件平台中添加粒子追踪模块就可以使用这个接口。虽然该模块的用户通常不会使用它来随机生成点，但出于可视化和美观的目的，我们决定在这个有趣的示例中使用它。
现在，我们来举例说明。一些粒子被随机释放到一个正方形区域并保持静止。对位于正方形内切圆区域内的粒子数量进行追踪，来获取圆周率的实时估计值。可以看到，随着点数的增加，估计值（蓝色实线）逐渐接近真实值（绿色虚线）。值得注意的是，估计值的精确度并不随点的数量呈线性变化。蒙特卡罗近似的统计误差通常与 1/sqrt(n) 成正比。这意味着，要将误差减少 10 倍，通常需要将点数增加 100 倍。
在随机放置的点数不断增加的情况下，圆周率的实时估计值（蓝色实线）与真实值（绿色虚线）的比较。
接下来，我们使用 COMSOL Multiphysics®中的 App 开发器创建了一个基于多物理场仿真模型的仿真 App。在这个 App 中，我们可以使用一个滑块改变点的数量，并获得圆周率在不同点数的估计值以及与真实值的误差。该 App 还将随机放置的点可视化，并通过颜色协调来识别位于圆内的点。
使用仿真App根据不同的点数估算圆周率，并获得对结果的可视化解读。
下一步
欢迎从 COMSOL 案例库下载包含 App 设计和相关文件的 MPH 文件：使用蒙特卡洛法估算圆周率值。
更多使用蒙特卡罗方法的模型:
演示如何直接对 Ishigami 函数进行蒙特卡罗模拟的示例
使用蒙特卡洛法评估自由分子流动体系中简单涡轮分子泵性能的模型
使用蒙特卡洛法计算通过 S 形弯曲几何结构的传输概率的模型
有关 App 开发器的信息:
参考文献
P. Beckmann, A History of π. New York:St. Martin’s Press, 1971
C. Rossi, Architecture and Mathematics in Ancient Egypt.Cambridge University Press, 2004
C. Krishna, A profile of Indian culture.Indian Book Company, 1975
D.B. Damini & A. Dhar, How Archimedes showed that pi is approximately equal to 22/7.arXiv e-prints, 2020
Y. Lam & T.S. Ang, Circle measurements in ancient China,Historia Mathematica, 1986
本文原文来自COMSOL官方博客