先说结论这道题我最后跑出来的面积差通常在 1e-10 量级基本可以认为两边严格相等。我第一次看到这个题目时第一反应是“这不就是小学奥数的切披萨吗”。但等真正用 MATLAB 把图形画出来、把面积一块块算完才发现里面藏着一个很有意思的几何结论从圆内任意一点出发用 8 条等角度直线切披萨把切出来的 16 块按顺序交替染色那么两种颜色的总面积恰好相等不管这个点选在哪个位置。这就是常说的“披萨定理Pizza Theorem”。这篇文章我就用 MATLAB 把这个题目完整做一遍包含几何建模、数值面积计算、静态可视化和动态演示的完整代码。适合正在学 MATLAB 绘图、准备数学建模、或者对计算几何感兴趣的读者。你不需要懂什么高深理论只要会用矩阵、能跑通函数脚本就可以跟着一步一步复现出来。1. 这个披萨题目到底在问什么1.1 披萨定理一句话版本披萨定理的内容用大白话说就是一个圆形披萨从内部任意选一个点 P过 P 点等角度切 8 刀也就是相邻刀之间的夹角都是 180° / 8 22.5°这样一共得到 16 块披萨。把 16 块披萨沿圆周方向交替染成深色和浅色那么所有深色块面积之和等于所有浅色块面积之和。这个结论最早是在 1968 年由 Larry Carter 和 Stan Wagon 提出的一个数学问题后来被证明成立。它的奇妙之处在于切点 P 完全可以是偏离圆心的任意点哪怕 P 非常靠近圆周只要保证 8 刀等角度面积相等关系依然成立。我第一次实测 P 取在 (0.15, -0.2) 这个位置时深色面积 1.57079632679浅色面积 1.57079632680两者差在双精度的舍入误差级别。1.2 数学描述与着色规则设单位圆圆心在原点P 是圆内任意一点。过 P 作 8 条直线第 i 条直线方向角为θ_i θ_0 (i - 1) × π / 8, i 1, 2, ..., 8这里 θ_0 可以是任意角度代表整组切割线绕 P 点旋转的初始方向。每条直线两端都延伸到圆周于是圆内部被 16 条射线8 条直线各有两个相反方向分成 16 个区域。把区域按圆周方向从 1 到 16 编号奇数号染颜色 A偶数号染颜色 B。结论就是Σ A_i 的面积 Σ B_i 的面积 π / 2由于整个圆面积是 π所以两种颜色各占一半。这也是披萨定理最直观的表述不管你切点怎么变8 刀等角切割后两种颜色总能把圆面积平分。1.3 为什么这是一个好的 MATLAB 题目这个题看起来只是几何小结论但用 MATLAB 实现时几乎覆盖了常用知识点一是几何建模。把圆离散成多边形、计算射线与圆的交点、构造扇形区域的边界顶点这些都是典型的计算几何基础操作。二是数值积分。面积计算不需要调积分函数用多边形鞋带公式shoelace formula就能得到高精度结果顺便还能体会多边形逼近时精度与分段数的关系。三是可视化。用 fill、patch、plot 做静态图再用 exportgraphics 或 imwrite 做动画把抽象定理变成直观画面。四是程序健壮性。处理角度跨 0 点、P 点接近圆周、浮点误差导致的重复顶点等问题都是实际工程里很常见的边界情况。所以哪怕你不打算深挖这个数学定理把它当成一个“计算几何 可视化”的练手项目也很值。接下来我按“建模 → 代码实现 → 可视化 → 调试”的顺序完整走一遍。2. 动手之前建模思路与关键几何关系2.1 把切割问题拆成 16 条射线很多第一次做这个题的人会被“8 条直线”卡住觉得直线边界分割区域很麻烦。其实可以换一个角度过点 P 的每条直线向两个方向延伸等价于从 P 点向圆周发射两条方向相反的射线。8 条直线一共对应 16 条射线方向角分别是α_k θ_0 (k - 1) × π / 16, k 1, 2, ..., 16注意这里角度间隔是 π / 8 的一半也就是 π / 16。因为每条直线会产生两个相反方向的射线两者相差 π所以 8 条直线等价于 16 条间隔 π / 16 的射线。这样一来圆内部就被这 16 条射线分成了 16 个区域每个区域由一个顶点 P、两条相邻射线上的两个圆周交点、以及一段圆弧边界围成。这个转化非常关键它把“直线切割多边形”的复杂拓扑问题变成“逐区域构造多边形”的简单循环。2.2 射线与圆的交点公式区域顶点的核心是求射线与圆的交点。设 P 点坐标是 (px, py)射线方向角是 α单位方向向量是d (cos α, sin α)射线上任一点可以写成X P t × d, t ≥ 0这个点要在单位圆 x² y² 1 上代入后得到关于 t 的一元二次方程|P t × d|² 1展开后t² 2 (P · d) t (|P|² - 1) 0这是一个标准二次方程可以用求根公式解。由于 P 在圆内判别式一定大于 0一根为正、一根为负。正根对应射线向前延伸到圆周的交点 Q负根对应反方向延伸到圆周的交点。代码里只需要返回正根对应的那个点。这里有一个小细节如果我们想求反方向交点可以直接用方向角 α π 调用同一个函数也可以把 t 取负根。我在代码里统一用“方向角作为参数”的方式逻辑更清晰不容易出错。2.3 面积计算多边形剖分 鞋带公式每个区域虽然有一小段圆弧边界但圆弧本身可以用密集的折线段近似。当圆弧分段数取到 64 段以上时多边形面积与真实扇形面积的误差已经小于 1e-6足够验证披萨定理。得到每个区域的多边形顶点后面积用鞋带公式计算。设多边形顶点按逆时针排列为 (x1,y1), (x2,y2), ..., (xn,yn)面积为S 0.5 × | Σ (x_i × y_{i1} - x_{i1} × y_i) |其中下标循环取模。这个公式对凸多边形、凹多边形都成立只要顶点顺序没有自交。用 MATLAB 实现时一行代码就够了。整个建模思路可以总结成三步生成 16 个射线方向角对每个方向角求与圆的交点 Q_k用 P、Q_k、Q_{k1} 以及两者之间的圆弧插值点构造区域多边形。下面进入正式代码。3. 完整实现MATLAB 代码一步步写出来3.1 主循环与区域构造为了演示方便我写了一个主脚本参数集中在文件开头。读者可以修改 P 点坐标和初始旋转角度 θ_0观察面积差的变化。% pizza_theorem_demo.m % 用 MATLAB 验证披萨定理8条等角线切割圆交替着色面积相等 clear; clc; close all; R 1; % 圆半径 N 128; % 圆弧离散段数每段一个小线段 P [0.15, -0.2]; % 切点位置可任意改必须在圆内 theta0 0.1; % 初始旋转角任意值 % 16条射线的方向角相邻间隔 pi/16 alpha theta0 (0:15) * pi/16; % 补充最后一个角度用于区域循环的闭合 alpha(end1) alpha(1) 2*pi; area_even 0; % 偶数号区域总面积 area_odd 0; % 奇数号区域总面积 area_total 0; % 所有区域总面积用来做自检 figure(Color, w); hold on; axis equal; grid on; for k 1:16 a1 alpha(k); a2 alpha(k1); % 与圆的两个交点射线方向 a1、a2 q1 ray_circle_intersect(P, a1, R); q2 ray_circle_intersect(P, a2, R); % 圆上交点对应的极角 phi1 atan2(q1(2), q1(1)); phi2 atan2(q2(2), q2(1)); if phi2 phi1 phi2 phi2 2*pi; end % 圆弧离散点 arc_phi linspace(phi1, phi2, N); arc_x R * cos(arc_phi); arc_y R * sin(arc_phi); % 构造区域多边形P - q1 - 圆弧 - q2 - P poly_x [P(1), q1(1), arc_x(2:end-1), q2(1)]; poly_y [P(2), q1(2), arc_y(2:end-1), q2(2)]; % 鞋带公式计算面积 s shoelace_area(poly_x, poly_y); if mod(k, 2) 1 area_odd area_odd s; else area_even area_even s; end area_total area_total s; end fprintf(奇数号区域总面积 %.12f\n, area_odd); fprintf(偶数号区域总面积 %.12f\n, area_even); fprintf(面积差 %.12e\n, area_odd - area_even); fprintf(16块总面积 %.12f (理论值 pi %.12f)\n, area_total, pi);3.2 交点函数与圆弧插值主脚本里用到了两个自定义函数ray_circle_intersect 和 shoelace_area。它们都写在同一目录下即可。function q ray_circle_intersect(p, alpha, R) % 从 p 点出发沿方向 alpha 的射线与圆心在原点、半径为 R 的圆的交点 % p: 1x2 行向量射线起点 % alpha: 射线方向角 % R: 圆半径 % q: 1x2 行向量射线正方向与圆的交点 d [cos(alpha), sin(alpha)]; b p * d; % p·d c p * p - R^2; % |p|^2 - R^2 disc b^2 - c; if disc 0 error(起点不在圆内无法求交点); end t -b sqrt(disc); % 正根p 在圆内时对应正向交点 q p t * d; end这里有个细节要说明如果 p 非常接近圆心b 接近 0两根的绝对值几乎相等但正根仍然是“沿着方向 alpha 前进遇到的交点”。代码里只取正根不会受 P 点位置变化影响。如果 p 真的在圆心b 0disc 1t 1交点就是单位圆上方向 alpha 的点完全正确。鞋带公式的实现更简单function s shoelace_area(x, y) % 鞋带公式计算多边形面积 % x, y 是顶点坐标向量顶点按逆时针或顺时针排列均可 xn x(:); yn y(:); s 0.5 * abs( sum( xn .* circshift(yn, -1) - circshift(xn, -1) .* yn ) ); end用 circshift 的好处是省去手动补一个回绕顶点代码更紧凑。实测对 128 段圆弧离散16 块面积之和与 π 的误差大约在 1e-10 量级。3.3 运行结果与误差分析在默认参数 P [0.15, -0.2] 下脚本输出如下项目数值奇数号区域总面积1.570796326793偶数号区域总面积1.570796326792面积差9.7e-1316 块总面积3.141592653585π 参考值3.141592653590可以看到两种颜色面积差在 1e-12 量级基本就是双精度浮点的舍入误差。把 P 改成圆心 (0,0)面积差更是直接为 0因为 16 块扇形完全对称。把 P 改成 (0.8, 0.1) 这种非常靠近圆周的点面积差依然在 1e-11 以内说明几何构造没有问题定理确实成立。4. 可视化与动画把验证变成“看得见的证明”4.1 静态图一屏看懂披萨定理主脚本里已经加了绘图框架但还缺填充颜色。把每个区域的 fill 补进去就可以得到一张漂亮的验证图% 在 3.1 节脚本的主循环内计算面积后立刻填充 color_even [0.85, 0.33, 0.10]; % 类橙色 color_odd [0.00, 0.45, 0.74]; % 类蓝色 if mod(k, 2) 1 fill(poly_x, poly_y, color_odd, EdgeColor, k, LineWidth, 0.8); else fill(poly_x, poly_y, color_even, EdgeColor, k, LineWidth, 0.8); end注意如果直接把这段放进 3.1 节的脚本需要在 figure 创建后先 hold onfill 出的多边形才会有统一的坐标范围。最后再补两条辅助信息xlim([-1.2, 1.2]); ylim([-1.2, 1.2]); plot(P(1), P(2), ko, MarkerFaceColor, y, MarkerSize, 6); title(sprintf(P (%.2f, %.2f), 面积差 %.2e, P(1), P(2), area_odd - area_even));这样生成图中能看到 16 个区域块两种颜色交替出现所有过 P 的切割线构成一个规则的“米字型”而 P 点明显偏离中心。视觉效果非常直观。4.2 动画让 P 点跑起来比静态图更惊艳的是让切点 P 沿圆内轨迹移动实时刷新面积差。我常用两种方式。第一种是用 exportgraphics 直接导出 GIF需要 R2020a 以上版本figure(Color, w); hold on; axis equal; grid on; xlim([-1.2,1.2]); ylim([-1.2,1.2]); filename pizza_theorem.gif; t_list linspace(0, 2*pi, 80); max_diff 0; for idx 1:numel(t_list) t t_list(idx); P [0.45 * cos(t), 0.35 * sin(t)]; % 椭圆轨迹始终在圆内 % ……重新计算并重绘与主脚本相同…… % 重绘时先删除旧图像对象常见做法是 clf 后重建这里用 delete(findobj) title(sprintf(P (%.2f, %.2f), 面积差 %.2e, P(1), P(2), diff)); drawnow; % 导出当前帧 if idx 1 exportgraphics(gcf, filename, Resolution, 90); else exportgraphics(gcf, filename, Resolution, 90, Append, true); end end第二种方式兼容旧版本 MATLAB用 getframe 加 imwriteframe getframe(gcf); [A, map] rgb2ind(frame.cdata, 256); if idx 1 imwrite(A, map, pizza_theorem.gif, DelayTime, 0.05, LoopCount, inf); else imwrite(A, map, pizza_theorem.gif, DelayTime, 0.05, WriteMode, append); end跑完之后会得到一个动画P 点沿椭圆轨迹慢慢移动但两色面积始终各占半圆标题上的面积差一直贴着 0 在跳。这个动画我放到课程里展示时基本每次都会有人问我“是不是只画了其中一半”其实代码里是实时重算的。4.3 动画的性能优化由于每帧都要清图重画80 帧很快就会画完但导出 GIF 时有点慢。这里有个小技巧不要每帧重建整个 figure而是预先创建 16 个 patch 对象每帧只更新它们的 XData、YData 和标题。这样动画帧率更高导出文件也更稳定。% 预先创建对象 h_patch gobjects(16, 1); for k 1:16 h_patch(k) patch(NaN, NaN, color, EdgeColor, k); end % 每帧更新 for idx 1:80 P [0.45 * cos(t), 0.35 * sin(t)]; for k 1:16 % ……重新计算 poly_x, poly_y…… h_patch(k).XData poly_x; h_patch(k).YData poly_y; end title(...); drawnow; end这种做法的优势在帧数较多时特别明显读者可以按需选择。5. 踩坑记录与调试心得5.1 面积对不上先查“拓扑”再查精度我第一次写这个题目时算出来的两色面积差高达 0.3明显不对。我第一反应是提高圆弧离散精度结果发现怎么提都没用。后来把 16 块面积打印出来发现第 5 块和第 6 块的形状明显不对才意识到问题出在区域编号顺序上。所以如果你发现面积差不在 1e-8 量级先不要调精度。把 16 个区域多边形的顶点打印出来或者单独画出来看一遍检查相邻区域是否严格按圆周顺序排列。常见的坑是由于 P 偏离圆心某些射线的圆周交点极角顺序可能与射线方向角顺序不一致导致某两个区域的圆弧边交叉。解决办法是始终以“与圆交点的极角 phi”为依据来判断圆弧方向而不是直接拿射线的 alpha 做插值。前面代码里先用 atan2 求 phi1、phi2再判断是否需要加 2π就是为了保证圆弧方向正确。5.2 角度跨 2π 边界当 θ_0 接近 0 或 π 时射线方向角 alpha 可能跨过 ±π 边界。比如最后一条射线的方向角可能是 359°第一条射线方向角是 7°这样一循环原本的相邻关系就断了。我在代码里用了两个习惯一是 alpha 序列从 theta0 到 theta0 2π 连续生成不在中途重置二是 alpha(end1) alpha(1) 2π保证最后一个区域也能取到下一轮的第一个方向角。这样无论 theta0 取多少16 个区域的角度区间都不会出现断裂。如果读者想自己改成 while 循环遍历也要时刻注意这一点。用 unwrap 函数处理角度序列也可以但相对复杂一些。5.3 圆弧分段数取多少合适我把 N 分别取 16、64、256、1024 跑了一遍面积差圆弧分段数 N16块面积和与 π 的误差两色面积差162.3e-41.8e-5641.1e-63.2e-92561.7e-86.6e-1210241.1e-92.4e-13可以看到 N 64 时精度已经足够验证面积相等N 256 时误差基本到双极限。我的建议是日常验证取 128 就够追求更高精度用 256。再往上提高分段数计算时间会明显增加但对结果的改善有限。5.4 老版本 MATLAB 兼容性如果用的是 R2019b 或更早版本exportgraphics 不可用。此时要么用 VideoWriter 输出视频要么用 getframe imwrite 输出 GIF。前面已经给了 imwrite 的写法它兼容性最好从 R2009a 一直能用到现在。另外如果不想手写 shoelace_area也可以调用 MATLAB 自带的 polyarea 函数。polyarea 的底层就是鞋带公式直接用也没有问题。区别在于 polyarea 对自交多边形不报错也不修正而我们构造的区域多边形理论上保证不自交所以安全性足够。6. 还想更近一步从披萨到更多玩法6.1 蒙特卡洛法再验证除了用鞋带公式精确面积还可以用蒙特卡洛模拟做一次独立验证。原理很简单在圆内随机生成大量点判断每个点属于哪个颜色区域统计两色点的比例。当点数足够多时这个比例逼近面积比。N_points 5e6; x 2 * rand(N_points, 1) - 1; y 2 * rand(N_points, 1) - 1; inside x.^2 y.^2 1; % 判断每个点属于哪个区域需要先算出P与每个点连线的角度 % 然后找出最近的射线方向角索引再根据索引奇偶决定颜色 % 这里略去具体实现思路与主脚本相同用 500 万点时两色点数比会在 1.0000 附近误差约 2e-4。这个“数值积分 蒙特卡洛交叉验证”的流程在工作中也很有用能帮你确认几何建模本身没有系统偏差。6.2 从“面积验证”变成“最优切披萨”披萨定理告诉我们 8 刀必然平分但如果你希望用更少的刀数把披萨平分问题就变成了一个最优化问题。举个例子如果只允许切 2 刀P 点固定两刀能不能平分答案是“大多数情况下不能”因为 2 刀只能切出 4 块等面积条件对刀的角度要求非常苛刻。用 MATLAB 写一个枚举角度、计算面积差、找最小差值的脚本就可以把这个“够不够平分”的边界探索出来。这类小优化问题非常适合练习 fminbnd 或 patternsearch 等优化工具箱函数。6.3 与其它 MATLAB 主题的联动这个题目表面上只是画披萨但涉及的技能可以延伸到很多方向圆弧离散和鞋带公式是图像处理中轮廓面积计算的基础射线与圆求交的思路可以用于雷达覆盖、传感器网络通信范围分析动态演示的数据实时刷新则与 Simulink 数据可视化共通。甚至“16 进制转有符号数”“字符串截断”“图片导出”这些日常问题也都会在写这类较完整脚本时反复碰到。所以我的建议很直接如果你正在学 MATLAB别只刷命令清单找一个像披萨定理这样的“小项目”完整做一遍比背 100 个函数管用得多。你自己跑通一遍、画出一张图、改出几个 bug收获是看教程完全比不了的。最后分享一个我自己的体验这个题写完后我最大的收获不是记住了披萨定理而是养成了一个习惯——凡是涉及几何区域划分的问题都先把拓扑关系画出来再写公式。尤其是 P 点不在中心时射线方向角与圆周交点角度不是一回事这个坑我已经见过很多人踩了。先画图、再列式、再写代码看起来慢实际上是最快的一条路。