Octave / MATLAB 第97题：蒙特卡洛随机投点法估算圆周率 pi（向量化几何模拟）

题目

编写程序，使用蒙特卡洛方法（Monte Carlo Simulation）估算圆周率 $\pi$：在单位正方形 $[0, 1] \times [0, 1]$ 内全向量化生成 $N = 100,000$ 个独立均匀随机散点，统计落在单位圆四分之一圆弧内的散点比例，估算 $\pi$ 并评估收敛误差。

输入格式：无需外部输入。

输出格式：输出总投点数、落在圆内点数、估算的 $\pi$ 值以及与真值的绝对偏差。

数据范围与约定：蒙特卡洛随机模拟。

思路与算法

1. 几何概率原理
正方形面积为 1，四分之一圆弧面积为 $\frac{\pi}{4}$。投点在圆内的概率理论上等于 $\frac{\pi}{4}$。因此 $\pi \approx 4 \times \frac{N_{\text{inside}}}{N_{\text{total}}}$。

2. 全向量化实现
一次性生成 $N$ 维随机坐标向量 x = rand(N, 1) 与 y = rand(N, 1)；通过判定 dist_sq = x.^2 + y.^2 <= 1.0 获取逻辑掩码，消除任何慢速 for 循环。

3. 大数定律收敛
估算精度按 $O(1/\sqrt{N})$ 收敛。$N = 100,000$ 级别通常可稳定达到 2~3 位有效精度。

完整程序

% Octave / MATLAB 第97题：蒙特卡洛向量化估算圆周率 pi
% 固定随机种子保证结果可复现
rand("seed", 42);

N = 100000; % 投点总数
x = rand(N, 1);
y = rand(N, 1);

% 判定是否落在四分之一圆内: x^2 + y^2 <= 1
is_inside = (x.^2 + y.^2) <= 1.0;
inside_count = sum(is_inside);

% 估算 pi
pi_estimated = 4 * (inside_count / N);
abs_error = abs(pi_estimated - pi);

fprintf("蒙特卡洛投点总数: %d\n", N);
fprintf("落入圆内散点数:   %d\n", inside_count);
fprintf("估算的 pi 结果:   %.5f\n", pi_estimated);
fprintf("理论真值 pi:      %.5f\n", pi);
fprintf("绝对误差:         %.5f\n", abs_error);

运行示例

输出：

蒙特卡洛投点总数: 100000
落入圆内散点数:   78380
估算的 pi 结果:   3.13520
理论真值 pi:      3.14159
绝对误差:         0.00639