用Python数值积分估算π：蒙特卡洛、梯形法与辛普森法实战对比

在数值计算这个圈子里算 π 几乎是最经典的开胃菜。不管是为了练手、验证算法还是给团队新人做入门培训用积分法估算 π 总是一个绕不开的话题。今天我把这条路从头到尾走一遍从数学原理到 Python 代码再到误差分析和踩坑记录希望能给你一份可以直接拿去用的参考。先说清楚这里讲的“积分法”指什么。严格来说围绕 π 和积分的关系有两套完全不同的玩法一套是蒙特卡洛随机投点利用概率积分去估计圆的面积另一套是数值求积把 ( \int_0^1 \frac{1}{1x^2} dx \frac{\pi}{4} ) 这个定积分用梯形法、辛普森法这类确定性算法算出来。两条路都能算出 π但收敛速度、实现的复杂度和适用场景差异巨大本文会把两种方法都拆开讲透并附上完整可复现的代码和实测对比。如果你在学数值分析、准备算法面试或者只是单纯想看看“代码怎么把数学变成数字”这篇文章都适合你。读完你不仅能跑出 π 的近似值还能理解为什么有些方法跑得再久也不够精确为什么另一些方法区间数加到几百就已经稳得离谱。1. 思路拆解为什么积分能算 π这个题目乍看有点绕π 明明是个几何常数跟积分有什么关系但只要你把 π 的定义和积分的几何意义放在一起看关系立刻变得非常直接。1.1 两个祖师爷级别的积分公式第一个思路来自圆面积。半径为 1 的圆面积是 π那四分之一圆的面积就是 π/4。如果用定积分表示这个四分之一圆就是函数 ( y \sqrt{1 - x^2} ) 在区间 [0, 1] 上与 x 轴围成的面积[ \int_0^1 \sqrt{1 - x^2} , dx \frac{\pi}{4} ]这是最直观的几何解释。只要能用数值方法求出这个定积分再乘以 4就得到 π 的近似值。第二个思路来自反正切函数的泰勒展开。( \arctan(x) ) 的导数是 ( \frac{1}{1 x^2} )而 ( \arctan(1) \frac{\pi}{4} )于是[ \int_0^1 \frac{1}{1 x^2} , dx \frac{\pi}{4} ]这个公式的好处是函数是一条光滑的、单调递减的曲线没有根号数值求积时误差更容易控制。实际操作里我更推荐用这一个式子它在数学上干净在计算上也少很多麻烦。1.2 积分数值化的两条路线有了积分等式接下来要解决的核心问题只有一个怎么在计算机上求这个定积分的近似值计算机不会解解析式只能做离散化。路线一叫做随机模拟法也就是蒙特卡洛方法。思路是在一个 1×1 的单位正方形里随机撒点统计落在四分之一圆内的点占全部点的比例。因为均匀撒点时落点落入某个区域的概率等于该区域的面积占比所以大量投点之后圆内点数除以总点数就逼近 π/4 这个面积值。这个方法的定义域很广稍后我们重点介绍。路线二叫做确定性求积法也就是用梯形、辛普森这类数值积分公式。思路是先把 [0,1] 区间切成一堆小区间在每个小区间上用简单函数直线或抛物线去近似原函数再把这些小面积加起来。区间切得越细近似就越准确。两条路线各有优劣。蒙特卡洛实现起来最无脑但收敛速度是 O(1/√N)N 是投点数想多一位精度就要增加 100 倍的计算量。数值求积里梯形法是 O(h²) 误差辛普森法是 O(h⁴) 误差h 是步长。后者在算 π 这种光滑函数时效率远超蒙特卡洛这也是为什么后面你会看到一个区间数只需要几千就能算到小数点后十几位。2. 蒙特卡洛投点法原理、实现与收敛速度蒙特卡洛方法虽然在实际工程里常被当作“最后手段”但它完美地体现了概率论与积分之间的联系也最适合拿来理解为什么随机算法会有方差、需要大量样本。2.1 核心原理用频率逼近概率先看一张“虚拟画面”你在一个边长 1 的正方形内随机撒黄豆同时画出这个正方形内嵌的四分之一圆。如果豆子落点完全均匀那么按理说落在圆内的豆子数占总豆子数的比例应该等于圆面积占正方形面积的比例。正方形面积是 1四分之一圆面积是 π/4所以有[ \frac{N_{\text{inside}}}{N_{\text{total}}} \approx \frac{\pi}{4} ]整理一下就得到 π 的估计式[ \pi \approx 4 \times \frac{N_{\text{inside}}}{N_{\text{total}}} ]这里有个容易混淆的点为什么不直接撒点统计一整个圆的面积因为计算机产生均匀随机数最简单的方式是生成 [0,1] 区间上的随机数所以通常只在第一象限投点最后乘以 4。如果你使用圆心对称的整个圆来投点原理不变但多了一个生成负数的步骤没有必要。2.2 代码实现从零写一个估算器代码不复杂核心就一个循环加一个判断。我用 Python 做了一个最简单的版本完全不用第三方库只用标准库里的 random。import random def estimate_pi_mc(num_points: int) - float: inside 0 for _ in range(num_points): x random.random() y random.random() if x * x y * y 1.0: inside 1 return 4.0 * inside / num_points这里判断x*x y*y 1.0就是判断点是否落在以原点为圆心、半径为 1 的圆内。运行一下分别投 1 万、10 万、100 万个点我这边一次典型输出是1 万点3.137610 万点3.14312100 万点3.141836你多跑几次会发现每次结果都在变这正是蒙特卡洛方法的特点它是一个随机估计量本身的方差就存在。2.3 收敛慢但思路开阔误差怎么估蒙特卡洛方法的标准差可以通过概率知识直接算出来。记每次投点是否落在圆内为随机变量 X其期望为 p π/4方差为 p(1-p)。根据中心极限定理用 N 个点估计 p 的标准误差大约是[ \sqrt{\frac{p(1-p)}{N}} ]换算成 π 的估计误差就是 4 倍的这个值。代入 p ≈ 0.7854化简后大约是[ \pi_{\text{err}} \approx \frac{1.653}{\sqrt{N}} ]这意味着每提升一位小数精度N 需要扩大约 100 倍。我从 100 万点到 1 亿点计算量翻了 100 倍精度大概只从小数点后 3 位提升到 4 位左右非常不划算。所以如果你追求高精度蒙特卡洛显然不是首选如果你需要在高维积分或者没有明确函数表达式的问题上求期望值那它几乎是唯一通用手段。这就是“收敛慢但思路开阔”的含义。3. 数值积分法梯形与辛普森求 π走完蒙特卡洛这条路接下来看看确定性数值积分方法。与随机撒点不同这类方法用规则图形去逼近函数曲线下方的面积。只要函数足够光滑误差往往能做到非常小。3.1 用 1/(1x^2) 积分算 π 的原理我们选用公式[ \int_0^1 \frac{1}{1 x^2} , dx \frac{\pi}{4} ]这个式子在数值上比根号函数友好函数没有不可导点也没有垂直切线而且在 [0,1] 区间上变化平缓从 1 单调降到 0.5。这样的函数用多项式逼近误差会非常理想。把区间 [0,1] 切成 n 等份每份宽度 h 1/n记节点 ( x_i i \times h )函数值 ( f_i \frac{1}{1 x_i^2} )。接下来要做的就是通过这些离散点构造近似面积。3.2 梯形法用直线贴曲线梯形法的几何意义非常朴素把每个小区间上的曲线段用连接两端的直线段代替形成一个下底、上底、高的梯形然后求面积。整个积分近似为[ \int_a^b f(x) , dx \approx \frac{h}{2} \left[ f(x_0) 2\sum_{i1}^{n-1} f(x_i) f(x_n) \right] ]为什么中间点的权重是 2因为每个内部点同时作为左边梯形的右底和右边梯形的左底被计算了两次所以合并后系数是 2。端点只属于一个梯形系数是 1。实现代码如下def pi_trapezoid(n: int) - float: h 1.0 / n total 0.0 for i in range(n 1): x i * h f 1.0 / (1.0 x * x) if i 0 or i n: total f else: total 2.0 * f return 4.0 * total * h / 2.0用 n 1000 跑一次结果已经能稳定到 3.1415926 附近。梯形法的误差量级是 O(h²)理论上区间数从 1000 提高 10 倍到 10000误差缩小到原来的 1/100效果非常明显。3.3 辛普森法用抛物线让精度起飞梯形法用直线逼近曲线本质上只利用了函数的一次多项式信息。辛普森法的思路升级了一层在每个小区间上不再用直线而是用一个二次抛物线去拟合函数。为了确定一条抛物线需要三个点所以辛普森法会把区间分成偶数份每两个小区间作为一个整体。公式为[ \int_a^b f(x) , dx \approx \frac{h}{3} \left[ f(x_0) 4\sum_{\text{odd}} f(x_i) 2\sum_{\text{even}} f(x_i) f(x_n) \right] ]这里奇数下标的系数是 4偶数下标不含端点的系数是 2端点系数是 1。记忆口诀很简单“端点 1奇 4偶 2”。实现时注意 n 必须为偶数否则会导致每个抛物线区间无法完整覆盖def pi_simpson(n: int) - float: if n % 2 ! 0: raise ValueError(n must be even) h 1.0 / n total 0.0 for i in range(n 1): x i * h f 1.0 / (1.0 x * x) if i 0 or i n: total f elif i % 2 1: total 4.0 * f else: total 2.0 * f return 4.0 * total * h / 3.0只用 n 100 个区间辛普森法就能算出 3.1415926536 左右比梯形法用 1000 个点还要准得多。它的误差量级是 O(h⁴)这就是为什么用抛物线拟合光滑函数的效果如此显著。注意辛普森法要求被积函数在区间上足够光滑至少要有连续的四阶导数。对于 1/(1x²) 这个函数条件完美满足。如果你换成一个带尖角的函数比如 |x|辛普森法的精度优势就会大打折扣。4. 实操全记录三套方案跑一遍理论讲完了落到实际工程里还需要考虑实现细节、运行效率、可视化验证等一堆问题。我照着自己常用的实验流程给出一份可以直接复制的完整对比测试包含代码、结果和精度分析。4.1 完整对比脚本与一次运行结果为了公平对比我把三种方法封装成同一个调用接口并统一输出绝对误差。绝对误差用 abs(estimate - math.pi) 计算这样可以直观看到每种方法距离真实 π 有多远。import math import random import time def estimate_pi_mc(num_points: int) - float: inside 0 for _ in range(num_points): x random.random() y random.random() if x * x y * y 1.0: inside 1 return 4.0 * inside / num_points def pi_trapezoid(n: int) - float: h 1.0 / n total 0.0 for i in range(n 1): x i * h f 1.0 / (1.0 x * x) if i 0 or i n: total f else: total 2.0 * f return 4.0 * total * h / 2.0 def pi_simpson(n: int) - float: if n % 2 ! 0: raise ValueError(n must be even) h 1.0 / n total 0.0 for i in range(n 1): x i * h f 1.0 / (1.0 x * x) if i 0 or i n: total f elif i % 2 1: total 4.0 * f else: total 2.0 * f return 4.0 * total * h / 3.0跑一次得到下面这张对比表我取了几组有代表性的参数方法参数投点数/区间数计算结果绝对误差耗时约蒙特卡洛1,000,0003.1418362.4e-041.2s蒙特卡洛100,000,0003.14161262.0e-05约 120s梯形法10003.14159248691.7e-070.002s梯形法1000003.14159265342.3e-100.08s辛普森法1003.14159265364.2e-120.0006s辛普森法10003.141592653589 许3.4e-140.006s可以看到同样是百万级别的计算量辛普森法的精度远超蒙特卡洛。蒙特卡洛哪怕用了 1 亿个点绝对误差还在 10⁻⁵ 量级辛普森法只切了 100 个区间误差就已经到了 10⁻¹²差距非常直观。4.2 提升精度的三个常用技巧在实际实验中有几个细节能明显改善结果稳定性。第一蒙特卡洛一定要引入固定随机种子。比如random.seed(42)。没设种子的话你每次跑出来的结果都不一样调试时很容易误判算法好坏。设了种子之后至少在同一份稳定代码下结果可复现方便对比不同参数。第二数值积分要用 Python 内置的 decimal 或者 numpy 的高精度数组来减少浮点累积误差。上面对比用的还是普通 float当 n 非常大时求和顺序会造成尾数损失。一个经验做法是不要从前往后加而是用 numpy 数组保存所有函数值再用np.sum()因为 numpy 的求和会做一部分误差补偿比逐项累加稳。第三对结果做误差事后估计。比如蒙特卡洛可以顺便计算方差数值积分可以分别用 n 和 2n 的结果做一次 Richardson 外推。外推公式对于梯形法很实用如果 T(h) 是步长 h 的结果那么改进后的值约等于 ( \frac{4T(h) - T(2h)}{3} )能额外提升一到两位精度。4.3 不同方法的收敛速度对比收敛速度是选择算法的核心依据。我把三种方法的误差随计算量变化画在了一张表里方便直观理解方法误差量级计算量加倍后的效果蒙特卡洛O(N^(-1/2))精度只提升约 40%梯形法O(n^(-2))精度提升到原来的约 1/4辛普森法O(n^(-4))精度提升到原来的约 1/16这个表解释了为什么工程上很少用蒙特卡洛去算一维定积分。一维问题有太多确定性算法又快又准蒙特卡洛真正的用武之地是高维积分比如 10 维甚至 100 维的数值积分这时确定性方法会面临维数灾难而蒙特卡洛误差与维度无关。所以项目中投点法适合拿来理解原理真到了要出数值结果的时候首选一定是梯形法或辛普森法。5. 常见问题与排查技巧实录实操中没有一个人是顺顺利利一次跑对的下面几个坑我全踩过逐个拿出来说说我的排查思路。这些问题也是初学者最容易困惑的几处。5.1 为什么我投了 100 万点结果还在 3.14 上下晃这是蒙特卡洛方法本身的统计涨落不是 bug。前面算过100 万点对应的误差标准差约 0.00165换算成 π 估计量就是标准差约 0.00165 的 4 倍除以 √N 关系实际常见偏差在 ±0.002 量级。所以看到结果在 3.139 到 3.144 之间浮动是完全正常的。如果想让结果更稳定第一是增加点数第二是使用低差异数列比如 Sobol 序列代替纯随机数这会显著降低方差。使用 numpy 生成 Sobol 序列需要安装scipy.stats.qmc官方文档有现成示例计算时间几乎不变但精度能提升一到两个量级。5.2 梯形法区间数拉满为什么反而不准了这是我调试时踩过的一个经典坑。把 n 从 1000 提到 100 万理论上误差应该继续减小但实际结果在某个区间数之后开始出现轻微抖动甚至误差变大。原因是浮点数累加的舍入误差开始累积超过 10 万次浮点累加总和的尾数误差可能与梯形法的截断误差处于同一量级导致精度不升反降。解决办法有三个方向使用math.fsum代替普通循环累加它能做高精度舍入补偿把求和改成二分段求和或者使用 Kahan 求和算法统一用 numpy 数组加np.sum它的实现内部做了配对求和。我在代码里用math.fsum跑 n 100 万结果比普通 for 循环累加干净很多。5.3 辛普森法报错说 n 必须为偶数但我明明传了偶数这种情况多半出现在你使用动态传参时。比如从配置文件读取 n不小心读成了字符串然后n % 2就报类型错误。更微妙的情况是传了浮点数比如n 100.0100.0 % 2 0.0虽然在 Python 里是 True但循环里i * h的浮点误差会让最后一步计算不稳定。我的建议是一进来就强制转成 int并且校验n 0和n % 2 0。别嫌麻烦这几行保护在真实项目里能帮你省下大量调试时间。5.4 蒙特卡洛结果与别人跑出来的完全不一样大概率是随机数种子问题。不同 Python 版本、不同操作系统随机数生成器的内部状态可能不同。如果团队协作或者写博客做演示一定要固定种子否则结果不可复现问题很难定位。另外一个隐蔽因素是随机坐标的范围。如果你生成 x 和 y 时用了不同的随机区间比如 x 在 [-1, 1]y 却在 [0, 1]面积比例公式就会彻底失效。这种错误不会让代码崩溃只会让你得到某个奇怪的数值。排查时首先检查投点范围是否一致再检查圆内判断条件。5.5 排查速查表症状可能原因解决方案蒙特卡洛结果每次不同未固定随机种子或样本量不足设置 seed提高点数至百万以上蒙特卡洛长期偏离 3.14 以上坐标范围不一致圆内判断写错检查 x、y 是否都取 [0,1]检查 x²y²1梯形法 n 增大误差反而增大浮点累加舍入误差累积使用 math.fsum、numpy sum 或 Kahan 求和辛普森法报 n 必须为偶数n 被传入浮点数或字符串强制 int 并校验 n % 2 0梯形法和辛普森结果几乎一样函数不光滑或 n 太小改用更高阶公式或检查被积函数性质计算结果稳定但总是差 4 倍忘记乘 4或面积公式比例用反核对 π 4 × 圆内点比例/总面积比例6. 我的实际操作体会最后分享一点个人经验这类项目别看表面简单真正动手跑一遍你对“误差”这件事的理解会上升一个台阶。我建议你别只复制代码而是亲手改改参数、画一下误差曲线看看不同方法的误差下降趋势有多不一样。蒙特卡洛那条误差曲线总是抖抖索索往下走梯形法和辛普森法却能快速触底这种直观感受比看任何教科书都来得牢。另外一个非常实用的扩展方向是把代码改成 numpy 向量化版本。蒙特卡洛方法用数组运算一次生成所有随机点数值积分用np.linspace生成节点、np.sum求和运行速度快几个数量级。配合 seaborn 或 matplotlib 画出投点图和函数逼近图你就得到了一份既能展示原理又能演示性能的完整小项目。如果你对数值积分进一步感兴趣还可以试试高斯-勒让德求积法。它用更少节点达到更高的代数精度是很多现代计算库背后的核心算法。算 π 只是一个入口理解了这些数值方法后面做复杂积分、微分方程求解时会有大量用武之地。