使用蒙特卡罗方法和粒子追踪估算圆周率

2024年 5月 24日

如果说有一个数字可以统治所有的数学和科学,那就是圆周率。π 这个小小的符号有着悠久的历史,可以追溯到数千年以前。从对圆进行粗略近似的古代文明,到计算万亿位数的现代超级计算机,圆周率一直吸引着数学家和好奇者的想象力。这篇博客,我们将通过COMSOL Multiphysics®仿真软件提供的功能,以一种有趣和流行的方法来计算圆周率。

历史上对圆周率的近似计算

已知最早的圆周率近似值出现在古代文明中。巴比伦的数学家将圆周率近似为 3,这个数值在当时的建筑工程中是合理的,后来又被精确为 3.1251。埃及的数学家和印度的学者分别通过比较圆形和八边形的面积2 ,和通过巨量的计算3得出了近似的数值。包括阿基米德在内的一些希腊学者利用几何方法将圆周率的近似值精确到 3 个数量级以内4,使圆周率的计算取得了重大突破。

一幅 Archimedes的绘画。
 
Fibonacci的肖像画。

左图: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。

半径为 r 的圆嵌于边长为 2r 的正方形中。
估计圆周率的基础。

在 COMSOL Multiphysics® 中使用蒙特卡罗方法估算圆周率

为了进行这个简单的蒙特卡罗模拟,我们将使用 数学粒子追踪 接口。在 COMSOL Multiphysics® 软件平台中添加粒子追踪模块就可以使用这个接口。虽然该模块的用户通常不会使用它来随机生成点,但出于可视化和美观的目的,我们决定在这个有趣的示例中使用它。

现在,我们来举例说明。一些粒子被随机释放到一个正方形区域并保持静止。对位于正方形内切圆区域内的粒子数量进行追踪,来获取圆周率的实时估计值。可以看到,随着点数的增加,估计值(蓝色实线)逐渐接近真实值(绿色虚线)。值得注意的是,估计值的精确度并不随点的数量呈线性变化。蒙特卡罗近似的统计误差通常与 1/sqrt(n) 成正比。这意味着,要将误差减少 10 倍,通常需要将点数增加 100 倍。

图中 Y 轴为估计值,X 轴为粒子数,蓝色实线上下波动,绿色虚线在所有粒子值中都保持在略高于 3.14 的水平。
在随机放置的点数不断增加的情况下,圆周率的实时估计值(蓝色实线)与真实值(绿色虚线)的比较。

接下来,我们使用 COMSOL Multiphysics® 中的 App 开发器创建了一个基于多物理场仿真模型的仿真 App。在这个 App 中,我们可以使用一个滑块改变点的数量,并获得圆周率在不同点数的估计值以及与真实值的误差。该 App 还将随机放置的点可视化,并通过颜色协调来识别位于圆内的点。

使用仿真App根据不同的点数估算圆周率,并获得对结果的可视化解读。

下一步

欢迎从 COMSOL 案例库下载包含 App 设计和相关文件的 MPH 文件:使用蒙特卡洛法估算圆周率值


评论 (0)

正在加载...
浏览 COMSOL 博客