简介围绕蒙特卡罗方法在MATLAB中的落地应用这份文档面向正在学习数值计算、数学建模与随机模拟的初、中级读者借“计算阴影区域面积”这一具体问题演示均匀网格分割与随机投点两种实现思路并对照二者在精度与稳定性上的差异。包内仅含1个doc文件约33KB篇幅精炼集中呈现问题背景、算法分析、MATLAB源码及运行结果便于直接阅读与复现。已有1626人学习下载。读者可从中获得蒙特卡罗求面积的完整推导脉络、两套可直接运行的脚本、循环与向量化写法示例以及通过增大样本量、多次取均值控制误差的实践思路适合作为课程作业、建模入门或算法对照的参考材料。1. 从一袋豆子到 7.2150这块阴影面积为什么值得用蒙特卡罗算一遍有人贴出一张图两个椭圆加一个圆三者重叠出一块阴影问面积怎么算。纯解析解要解椭圆与圆的交点、再分段积分式子长到没人愿意手推最后被采纳的答案是一段四十来行的蒙特卡罗 MATLAB 脚本跑出来7.2150。这就是命中率法hit-or-miss最典型的落地场景——不碰边界方程只在包围盒里大量撒点用落在区域内的点数比例乘以盒子面积。蒙特卡罗方法适合的正是这种边界由若干不等式定义、解析积分难以下手、又只要求两三位有效数字的场合。下面把这份实例拆开三个不等式怎么来的、均匀网格和随机撒点差在哪、参数怎么调、坑在哪。2. 阴影区域的数学描述三个不等式如何把图形锁进代码拿到图形先别写循环先把阴影翻译成能判断真假的布尔表达式。这一步错了后面撒一千万个点也是白费。2.1 三个几何体的半轴、圆心与叠加关系原始脚本里那三个a1、a2、a3不是随手写的每一个都对应一个标准二次曲线。把它们还原成几何语言变量表达式几何含义关键参数a1x²/9 y²/36 1椭圆长轴沿 y半轴 3x 向、6y 向a2x²/36 y² 1椭圆长轴沿 x半轴 6x 向、1y 向a3(x−2)² (y1)² 9圆圆心 (2, −1)半径 3阴影是两个椭圆的交集再与这个圆取交集。原脚本用了三层嵌套if逻辑上等价于三个条件同时成立。这里有个容易忽略的细节a1的 y 半轴是 6a2的 x 半轴也是 6所以整块图形被限制在 x ∈ [−6, 6]、y ∈ [−6, 6] 的方形里——这正是后面包围盒取 12×12 的来源不是随便拍的数。2.2 包围盒取 12×12 的理由与 5% 的命中率包围盒也叫参考区域必须完整盖住目标区域否则漏掉的点会被当成区域外面积结果系统性偏小。12×12 满足这个条件但远远不紧。两个椭圆的交集其实被压在 x ∈ (−3, 3)、y ∈ (−1, 1) 里面积只有 8圆的圆心在 (2, −1)右边界 x 5 也在椭圆交集之外不扩大范围。拿 144 去乘意味着命中率只有 7.2/144 ≈ 5%。换个说法每撒 100 个点95 个是废点。这就是为什么原文里 N 要开到 10000 才勉强稳定——不是蒙特卡罗本身慢是包围盒选得太松。这个数字在第 4 章算方差时会再次出现它直接决定了估计量的相对标准差。提示包围盒宁可先用解析方法收紧一次这比把 N 翻十倍划算得多。均匀步长撒点和随机撒点都吃这个亏。2.3 把嵌套 if 压成一行逻辑与三层if在 MATLAB 里意味着逐点标量判断144 万次循环就是这个写法堆出来的。判断三个区域是否同时命中本质是布尔与运算完全可以让运算符直接作用在矩阵上% 区域指示函数输入同尺寸的 x、y 矩阵返回 logical 矩阵 inShadow (x, y) ( x.^2/9 y.^2/36 1 ) ... % 椭圆一 ( x.^2/36 y.^2 1 ) ... % 椭圆二 ( (x-2).^2 (y1).^2 9 ); % 圆这里必须用而不是。是短路标量运算符只对单个 true/false 有效套在矩阵上会直接报 Operands to the || and operators must be convertible to logical scalar values。另外三处.^2的点号不能省省掉就变成矩阵幂1201×1201 的方阵自乘一次内存瞬间爆掉。这个匿名函数是后面所有代码的地基写完之后网格法和随机法都能复用同一份判定逻辑不用重复抄三遍不等式。3. 均匀网格法步长 0.01 背后的 144 万次判定方法一的思路是把豆子均匀摆在布上用固定步长的格点代替随机点。它的优点是结果可复现、没有随机波动缺点是精度完全由步长决定而且格点法在边界处有系统性偏差。3.1 原脚本逐行拆1201×1201 个格点从哪来clear x -6:0.01:6; % 每维 1201 个格点 y x; % 正方形包围盒两维取同一向量 s size(x); zs s(1,2)^2; % 总格点数 1201*1201 1442401 k 0; for i 1:s(1,2) for j 1:s(1,2) a1 (x(i)^2)/9 (y(j)^2)/36; a2 (x(i)^2)/36 y(j)^2; a3 (x(i)-2)^2 (y(j)1)^2; if a1 1 if a2 1 if a3 9 k k 1; % 三个条件同时满足才计数 end end end end end mj (12^2) * k / zs; % 命中比例 × 包围盒面积-6:0.01:6这个冒号表达式生成的是 1201 个点6−(−6) 除以 0.01 再加 1不是 1200 个。总格点 1442401双重循环跑 144 万次标量判定在普通笔记本上要几秒到十几秒。y x之所以成立是因为包围盒是正方形且两维步长一致一旦 x、y 的步长不同这个赋值就会让格点退化成一条对角线面积会错得离谱。zs用格点总数而不是格子数是格点法的一个约定每个格点代表它右上方那一小块面积边界格点会被整体算进或整体排除。3.2 meshgrid 向量化省掉双重循环后快多少同样的逻辑用meshgrid展开成矩阵判断一次成型clear; clc x -6:0.01:6; [X, Y] meshgrid(x, x); % X 沿列复制Y 沿行复制都是 1201×1201 in ( X.^2/9 Y.^2/36 1 ) ... ( X.^2/36 Y.^2 1 ) ... ( (X-2).^2 (Y1).^2 9 ); k nnz(in); % 命中格点数等价于原脚本里的 k area 12^2 * k / numel(in); % numel(in) 1442401 fprintf(k %d, area %.4f\n, k, area);meshgrid(x, x)返回的两个矩阵X的每一行都是xY的每一列都是x这样X(i,j)和Y(i,j)才构成格点坐标 (x_j, y_i)。nnz直接数非零元素个数等价于把布尔矩阵求和比sum(in(:))少一次类型转换。numel(in)就是总格点数写死 1442401 也能跑但换步长就会对不上。这段代码的关键参数只有两个步长 0.01 决定numel包围盒边界 −6 和 6 决定乘的系数 144。改步长的时候记得同步改zs或改用numel我见过有人把步长改成 0.02 却忘了改分母结果面积直接翻四倍。内存上1201×1201 的 double 矩阵是 1442401 × 8 ≈ 11.5 MBmeshgrid一次生成两个就是 23 MB加上三个中间布尔矩阵每个约 1.4 MB峰值在 30 MB 量级。这个规模在现在的机器上毫无压力但如果把步长压到 0.001单矩阵就是 1.15 GBmeshgrid会直接把内存吃满——网格法的天花板不是算力是内存。3.3 步长、内存与 O(h) 量化误差的取舍网格法本质上是用格点上的值近似每个小方格的平均值边界处必然有一圈格子被判错明明一半在区域内却因为格点落在外面而被整体排除。误差量级是周长 × 步长所以它随 h 线性下降而不是平方下降。步长 h每维点数总格点数单个 double 矩阵内存结果趋势0.10121146410.11 MB偏差在 1e-1 量级0.05241580810.46 MB偏差收敛到 5e-2 量级0.026013612012.9 MB偏差 2e-2 量级0.011201144240111.5 MB输出 7.2150从表里能看出一条经验线步长减半误差大致减半但内存翻四倍。想把精度再提高一位就得付出 100 倍的内存代价——这就是方法一在实操中很快撞墙的原因。它真正的价值在于可复现同一份代码在任何机器上跑出来都是 7.2150不像随机法每次都不一样。做对拍验证的时候我一般先跑网格法拿一个稳定参考值再去校随机法的结果。4. 随机撒点法N10000 重复 100 次均值为什么更稳方法二是标准的蒙特卡罗在包围盒里均匀掷随机点统计命中数。单次结果会飘所以原文重复了 100 次取均值。这个100不是拍脑袋它是有方差公式可算的。4.1 单次估计的波动用二项分布把 ±0.31 算出来每次撒点命中的个数 k 服从二项分布 B(N, p)p 就是真实面积占包围盒的比例约 0.05。二项分布的标准差是 √(Np(1−p))换算到面积上p 7.2 / 144; % 命中率约 0.05 N 1e4; sd_k sqrt(N * p * (1-p)); % 命中数的标准差 ≈ 21.8 sd_A 144 * sd_k / N; % 面积的标准差 ≈ 0.314也就是说N 10000 时单次跑出来的面积落在 7.25 ± 0.31一个标准差之内都算正常。有人看到两次结果差 0.4 就怀疑代码写错了其实完全在噪声范围内。把 N 提到 100000这个标准差降到 0.1 左右提到 100 万降到 0.03。代价是每提升一个数量级运行时间也翻十倍。相对标准差有个更直观的公式√((1−p)/(Np))。p 0.05 时它等于 4.36/√NN 10000 对应 4.36%这就是为什么要重复 100 次——100 次均值的标准差是单次的 1/10压到 0.44%。4.2 100 次重复的向量化写法与随机种子原脚本用外层for j 1:100套内层for i 1:N一共 100 万次标量循环。用矩阵一次铺开更省事rng(20240501); % 固定种子保证每次跑出同一结果 N 1e4; R 100; x 12*rand(N, R) - 6; % 每列是一次独立试验 y 12*rand(N, R) - 6; in ( x.^2/9 y.^2/36 1 ) ... ( x.^2/36 y.^2 1 ) ... ( (x-2).^2 (y1).^2 9 ); k sum(in, 1); % 1×100每列的命中数 est 144 * k / N; % 1×100每次试验的面积估计 mj mean(est); % 最终结果 se std(est) / sqrt(R); % 均值标准误 fprintf(area %.4f, se %.4f\n, mj, se);三个参数的含义要分清N是每次试验撒多少点决定单次精度R是重复多少次决定均值精度rng的种子决定结果能不能复现。sum(in, 1)里的1表示沿第一维行方向求和得到 1×R 的行向量如果写成sum(in)会默认也是按列求和结果一样但语义不清晰换个维度就容易翻车。12*rand(N,R) - 6把 [0,1) 的均匀随机数映到 [−6, 6)。注意rand取不到 1所以上边界 6 永远取不到——对连续分布来说这不影响面积但如果你在做离散格子模拟这个差一会让边界格子永远空着。4.3 要几位有效数字就得撒多少点把相对标准差 4.36/√N 反解一下就能按精度需求倒推采样量单次采样 N单次相对标准差配 100 次重复后适用场景1 00013.8%1.38%只做趋势演示10 0004.36%0.44%原文设置够出四位数字100 0001.38%0.14%需要和解析解对拍到小数点后两位1 000 0000.44%0.044%验证算法收敛性时用这张表的用法是反着查先想清楚结果要保留几位小数再定 N。如果只是交个作业7.25 这个量级足够如果要写进报告并和别的算法对比至少得把标准误一起报出来否则别人没法判断两个结果差 0.05 是算法差异还是随机噪声。注意std(est)/sqrt(R)才是均值的标准误直接写std(est)会把重复 100 次带来的收益全部抹掉报出来的误差大 10 倍。5. 避坑与排查五个让面积悄悄算错的细节这段代码短但出错的方式一点不少而且大部分错误不会报异常只会给你一个看起来差不多的数。5.1 逻辑与写成短路与现象改成之后要么直接报错要么k永远是个 0 或 1。原因是标量短路运算符遇到矩阵输入时 MATLAB 要求操作数能转成逻辑标量元素多于一个就抛错如果恰好是 1×1 的输入它只判断第一个元素后面的点全被忽略。解决矩阵逐元素判断一律用只有if条件里连接两个标量条件时才用。判断完可以用class(in)确认返回的是logical而不是double。5.2 meshgrid 参数顺序把圆心搬到镜像位置现象网格法算出来的面积和随机法差了 0.3 以上而且换个步长也不收敛到同一个值。原因写成[Y, X] meshgrid(x, x)或者meshgrid(y, x)坐标轴被转置。两个椭圆关于原点对称转置看不出问题但圆心在 (2, −1) 的圆转置后跑到 (−1, 2)交集形状完全变了。因为包围盒和椭圆都对称这个错误极具迷惑性结果只偏一点点。解决固定写[X, Y] meshgrid(x, x)并在判定前打印X(1,2)-X(1,1)和Y(2,1)-Y(1,1)确认前者是 x 步长、后者是 y 步长。5.3 不固定随机种子现象同一份代码连跑三次结果分别是 7.19、7.31、7.24被人质疑程序有 bug。原因rand默认按时间初始化每次启动都是新的随机流。在 N 10000、R 100 的设置下均值本身就有 ±0.03 的标准误三次结果差 0.1 属于正常波动。解决在生成随机数之前加rng(固定整数)。调试阶段固定种子确认逻辑正确正式出结果时如果想展示稳定性反而应该故意换几个不同种子把几次均值的离散程度报出来。5.4 把网格法当无偏估计现象步长从 0.02 改到 0.01面积从 7.21 变到 7.215再改到 0.005 又变回 7.21看起来在抖。原因格点法不是无偏估计。每个格点代表它右上方的一小块边界上的格子被整体算进或整体排除误差方向取决于边界穿过格子的位置所以结果会随步长振荡着趋近极限而不是单调逼近。解决别拿单一步长的网格法结果当真值。要么做步长减半的收敛序列看它是否稳定到某一位要么用integral2对同一个指示函数做数值积分对照那个结果和步长无关。5.5 N 太小就下结论现象N 取 500 时算出 6.8于是判断蒙特卡罗方法不准。原因N 500 时单次相对标准差是 4.36/√500 ≈ 19.5%6.8 距离 7.2 只差 5.5%完全在一倍标准差之内。这不是方法不准是样本量不够。解决先按第 4 章的表估一遍标准误再看偏差是否超过两倍标准误。超过才说明判定逻辑或包围盒有问题没超过就是正常噪声。这个习惯能省掉大量无效调试。6. 从这一块面积到通用积分器句柄封装、收敛验证与包围盒瘦身把判定逻辑抽成函数句柄之后同一套代码可以算任意区域不必每次重抄不等式。我一般的写法是把撒点、判定、估计、报误差封成一个函数输入区间接口和采样参数输出面积和标准误function [A, se] hitOrMiss(f, xlim, ylim, N, R) % f : 指示函数句柄输入两个同尺寸矩阵返回 logical % xlim : [xmin xmax]ylim : [ymin ymax] % N : 每次试验采样点数R : 独立重复次数 box (xlim(2)-xlim(1)) * (ylim(2)-ylim(1)); x xlim(1) (xlim(2)-xlim(1)) * rand(N, R); y ylim(1) (ylim(2)-ylim(1)) * rand(N, R); k sum(f(x, y), 1); A box * mean(k) / N; se box * std(k/N, 0, 2) / sqrt(R); end调用时把包围盒从 12×12 换成紧包围盒 [−3, 3]×[−1, 1]命中率从 5% 跳到 90%相对标准差从 4.36/√N 降到 0.33/√N同样的 N 精度直接提升 13 倍——这就是最省力的一次方差缩减不需要改算法只需要把盒子收紧f (x,y) ( x.^2/9 y.^2/36 1 ) ... ( x.^2/36 y.^2 1 ) ... ( (x-2).^2 (y1).^2 9 ); rng(7); [A1, se1] hitOrMiss(f, [-6 6], [-6 6], 1e4, 100); % 松包围盒 [A2, se2] hitOrMiss(f, [-3 3], [-1 1], 1e4, 100); % 紧包围盒验证环节我习惯再补一刀用integral2对同一个指示函数做数值积分两次结果对不上就说明撒点逻辑有问题而不是随机误差。注意被积函数是不连续的区域内外跳变默认的自适应方法会很慢加Method,tiled和AbsTol能稳不少g (x,y) double(f(x,y)); Aref integral2(g, -3, 3, -1, 1, Method, tiled, AbsTol, 1e-9);最后一步是收敛性检查把 N 按 10 的幂扫一遍画 log-log 图看斜率理论上误差应该以 N^(−1/2) 下降如果斜率明显偏离 −0.5多半是包围盒没盖全或者判定条件写漏了一项。这个检查做过一次之后我就不再靠看起来差不多来判断结果对不对了。从那以后我每次写蒙特卡罗类脚本都会先把指示函数单独抽出来、用小样本对拍一遍integral2再放开跑大 N——省下的调试时间远比多写那几行代码多。希望帮到你。本文还有配套的精品资源点击获取