MATLAB牛顿环动态仿真：从干涉原理到GIF动图生成全解析

简介：这份MATLAB光学仿真资源面向光电信息、物理及相关专业学习者，生动还原了牛顿环这一等厚干涉现象的动态形成过程。通过精确计算光程差并控制厚度变化，项目以动图形式呈现干涉环级数与明暗条纹随膜厚增加或减薄而连续演变，帮助理解曲率半径、折射率测量等光学测量原理。资源包内共4个文件，包含可直接运行的MATLAB脚本、动态演示GIF、等倾干涉与牛顿环对比分析的Word文档，以及含使用说明的TXT文件，压缩包整体约362MB。目前已有1794人学习下载。脚本可直接运行，支持调整初始厚度、变化速率等参数以对比不同条件下的干涉图样，适合课堂演示、实验教学与自学复盘；对比文档则清晰梳理了两类干涉在原理和应用上的差异，是光学数值模拟与应用研究的实用案例。

1. 光学仿真的第一道坎：牛顿环不是几个同心圆那么简单

光学仿真的第一个作业，多半是从牛顿环开始的。用 MATLAB 做动态仿真，把牛顿环的干涉图像做成一张动图，表面上是为了交作业时比静态图好看，实际上它逼你把干涉条件从公式翻译成像素级的相位计算——这是后面做波前仿真、干涉测量、薄膜设计都要过的坎。我见过不少人第一版就把中心那个黑斑做成白斑，不是代码写错，是半波损失压根没进公式。这篇文章从光程差推导开始，给出可复现的 MATLAB 代码，把曲率半径、波长这些参数逐个变成动画帧，最后导出 GIF 动图。适合正在做物理光学课程设计的人、想验证干涉理论的工程师，以及所有被「环为什么会动」困扰的初学者。

2. 先把干涉条件翻译成像素相位：牛顿环的离散化建模

2.1 光程差与半波损失：为什么中心一定是暗斑

牛顿环属于等厚干涉。平凸透镜放在平板玻璃上，中间夹着一层空气膜。以接触点为中心，半径 r 处的空气膜厚度 t 可以近似为 r² / (2R)，其中 R 是透镜曲率半径。这个近似只在 r 远小于 R 时成立——实验中 R 通常是 0.5 m 到几米，r 只有几毫米，误差在 1e-8 量级，完全够用。

反射光的干涉必须算上半波损失。光从空气射向玻璃，在玻璃表面反射时相位突变 π，相当于多走了半个波长。总光程差写成 Δ = 2t + λ/2，空气折射率近似为 1。暗环条件 Δ = mλ，代入得 2t = (m - 1/2)λ；亮环条件 Δ = (m + 1/2)λ，代入得 2t = mλ。把 t = r²/(2R) 代进去，就得到暗环半径和亮环半径的解析式：

• 暗环半径：r_k = sqrt((k - 1/2)λR)，k = 1, 2, 3...
• 亮环半径：r_k = sqrt(kλR)

当 r = 0 时，光程差只剩 λ/2，反射光干涉相消，所以牛顿环中心永远是暗的。这是判断仿真结果对不对的第一个锚点。后面写代码时会用 sin²(π r²/(Rλ)) 表示反射光强度，这个函数在 r=0 处取 0，天然满足中心暗斑。你要是用了 cos²，中心就变成白斑，典型翻车现场。

2.2 meshgrid 构建二维场并计算相位

连续公式翻译成 MATLAB 代码，核心就两步：用 meshgrid 生成二维坐标网格，然后对整个矩阵批量计算相位。不要用双重 for 循环逐点算，MATLAB 矩阵运算一次搞定。

lambda = 632.8e-9; % 波长，单位米，氦氖激光常见波长 R = 1.0; % 透镜曲率半径，单位米 L = 8e-3; % 模拟区域边长，单位米 N = 2000; % 每边采样点数 x = linspace(-L/2, L/2, N); [X, Y] = meshgrid(x); % 生成 N×N 二维坐标网格 r2 = X.^2 + Y.^2; % 半径平方矩阵 phase = pi * r2 / (R * lambda); % 反射光相位差的一半 I = sin(phase).^2; % 反射光干涉强度，中心为 0

meshgrid(x) 相当于 meshgrid(x, x)，生成两个 N×N 矩阵 X 和 Y。r2 是每个像素到中心接触点的距离平方。这里最容易被忽略的是坐标单位：linspace 生成的 x 是米，lambda 也是米，如果混用毫米和纳米，条纹密度会差几千倍。习惯上我在变量名后面直接注释单位，避免算完发现图像是一张纯色板。

2.3 采样率与像素标定：让理论环半径可以直接比对

生成图像之前，先算一组采样约束。相邻条纹间距在半径 r 处大约等于 Rλ/(2r)，越靠外沿条纹越密。模拟区域边长为 L，边缘处半径 r ≈ L/2，最小条纹间距 d_min ≈ Rλ/L。MATLAB 图像每个像素对应物理尺寸 dp = L/N，要让条纹看得清，至少满足 dp ≤ d_min / 2，也就是 N ≥ 2L²/(Rλ)。

拿上面参数算：L = 8e-3，R = 1，λ = 632.8e-9，得到 N ≥ 202。我取了 N=2000，采样余量很足。如果 R 小到 0.2 m，L 又开到 20 mm，同样公式要求 N 接近 10000，2000 就会看到摩尔纹。这个检查必须写在参数前面，而不是等图像出问题再回头调。

像素标定也很简单：每像素对应的物理尺寸 dp = L/N。暗环第 k 个的理论像素半径 = sqrt((k - 1/2)λR) / dp。用图像工具读出环的像素位置，跟这个值比对，误差在 1~2 像素内基本正常，超过 5 像素就该查坐标单位。

3. 静态牛顿环最小代码：20 行生成第一幅干涉图

3.1 单帧生成代码与运行顺序

前面已经把相位矩阵算好了，静态图只是加一个显示层。完整脚本如下，直接复制到脚本文件或实时编辑器都能跑：

close all; lambda = 632.8e-9; % 氦氖激光波长，单位米 R = 1.0; % 曲率半径，单位米 L = 8e-3; % 模拟区域边长，单位米 N = 2000; % 每边采样点数 x = linspace(-L/2, L/2, N); [X, Y] = meshgrid(x); r2 = X.^2 + Y.^2; I = sin(pi * r2 / (R * lambda)).^2; imagesc(x*1e3, x*1e3, I); % 坐标轴转成毫米显示 axis image; colormap(gray); title('牛顿环反射光干涉'); xlabel('x (mm)'); ylabel('y (mm)');

运行顺序上，先算 r2，再算 I，最后 imagesc。这里有个细节：imagesc 的 x 轴和 y 轴标签我用毫米显示，但 I 矩阵本身对应的物理坐标仍然是米，只是显示层做了缩放。axis image 保证像素是正方形，不会把圆环拉成椭圆。

3.2 imagesc 与 imshow 的方向、坐标与色阶差异

静态图出来了，接下来要留意显示函数的方向问题。imagesc(x, y, C) 显示时 C 的每一列对应 x 坐标，每一行对应 y 坐标。对牛顿环这种中心对称图像，转不转置肉眼看不出来，但如果你想叠加一条半径方向的剖面线，坐标对不上就会很痛苦。我一般会显式转置一次：imagesc(x1e3, x1e3, I.')，确保后续用 xline 或 plot 叠加时方向一致。

imshow 也可以显示，但它不做坐标映射，像素位置直接对应矩阵下标。如果代码后面要读取暗环像素位置，imshow 反而方便，因为它省掉了一次坐标换算。我的习惯是：出图给别人看用 imagesc，自己做验证用 imshow，然后直接在矩阵下标层面操作。

色阶方面，gray colormap 是 64 级灰度，干涉条纹包络很光滑。不要用 jet 或 parula，会把干涉图染成伪彩色，掩盖细节。如果想让暗环更锐利，可以把 colormap 改为 gray(256)，并配合 caxis 把强度范围压缩到 [0, 1]。

3.3 反射光与透射光互补：一份代码出两幅图

反射光和透射光强度是互补的。反射光用 sin²，透射光用 cos²，两者之和恒等于 1。这可以做一个很直观的对照实验：同一组 R、lambda，把 cos² 的结果画在旁边，你会看到反射光的暗环位置正是透射光的亮环位置。

I_ref = sin(pi * r2 / (R * lambda)).^2; I_tra = cos(pi * r2 / (R * lambda)).^2; subplot(1, 2, 1); imagesc(x*1e3, x*1e3, I_ref); axis image; colormap(gray); title('反射光'); subplot(1, 2, 2); imagesc(x*1e3, x*1e3, I_tra); axis image; colormap(gray); title('透射光');

注意中心位置：反射光中心是暗的，透射光中心是亮的。这个互补关系本身就是检查代码有没有写反的试金石。如果两幅图中心同为暗斑，那多半是半波损失处理出了问题。

4. 动态仿真牛顿环：扫描曲率半径并导出 GIF 动图

4.1 扫描参数怎么选：R、λ、膜厚三个方向的视觉差异

动态仿真的本质是让干涉条件里的某个物理量随时间变化。常见做法是扫描曲率半径 R、扫描波长 λ、或者给相位加一个随时间变化的倾斜项。三种扫描的视觉效果完全不同：

我最常用的是扫描曲率半径 R。原因有两个：一是 R 的物理意义直观，平凸透镜换一枚曲率半径更大的，环就变稀；二是 R 出现在分母上，数值变化对条纹密度的影响比 λ 更明显，动画效果更易观察。实际演示时，我一般让 R 从 0.5 m 匀速变化到 2.0 m，整个动画时长控制在 2 秒左右。

4.2 帧循环与 GIF 写入的完整脚本

动态仿真的脚本结构是：预先算好不变量（x、r2），在循环里更新变量参数，重新计算强度矩阵，抓取当前画面，追加写入 GIF 文件。下面这段是我常用的模板：

gifFile = 'newton_rings_R_scan.gif'; RStart = 0.5; % 起始曲率半径 REnd = 2.0; % 结束曲率半径 nFrames = 100; % 总帧数 fig = figure('Color', 'w'); for i = 1:nFrames alpha = (i - 1) / (nFrames - 1); % 对数扫描，避免前几帧相位变化太剧烈 R = RStart * (REnd / RStart)^alpha; I = sin(pi * r2 / (R * lambda)).^2; imagesc(x*1e3, x*1e3, I); axis image; colormap(gray); title(sprintf('R = %.2f m', R)); drawnow; frame = getframe(fig); [ind, cmap] = rgb2ind(frame2im(frame), 256); if i == 1 imwrite(ind, cmap, gifFile, 'gif', 'LoopCount', Inf, 'DelayTime', 0.02); else imwrite(ind, cmap, gifFile, 'gif', 'WriteMode', 'append', 'DelayTime', 0.02); end end

这段代码有几个参数需要说明。nFrames 取 100，配合 DelayTime=0.02，播放时长约 2 秒，看起来流畅又不卡顿。LoopCount=Inf 让 GIF 无限循环。R 用对数扫描而不是线性扫描，这点很关键——线性扫描时 R 从 0.5 到 2.0，前几帧 R 变化比例很大，条纹移动速度不均匀，看起来像「先快后慢」；对数扫描保证每帧相位变化比例恒定，动画节奏均匀。

4.3 让环「平滑扩散」而不闪烁：相位步进的三个约束

动态仿真的坑不在代码能不能跑，而在画面会不会闪烁。闪烁的根源是相邻帧之间条纹移动距离超过半个条纹周期，视觉上会看到明暗环在抖动而不是平滑扩散。要避免这个，必须控制相位步进。

以扫描 R 为例，边缘处像素的相位 φ ≈ π (L/2)² / (Rλ)，R 变化时相位变化量近似为 Δφ = φ · |ΔR| / R。要让 Δφ 小于 π，帧数得满足 nFrames ≥ φ_max · ln(REnd/RStart) / π。拿前面的参数算一下：φ_max 约 159 rad，ln(4) / π 约 0.44，算出来至少需要 70 帧。这解释了为什么 30 帧的 GIF 会闪，100 帧就基本平滑。

除了帧数，还有两个小习惯能减少闪烁。第一，GIF 写入前用 rgb2ind 统一减色到 256 色，并保持 colormap 不变，不然每一帧颜色表都变，视频会像霓虹灯一样跳。第二，不要每帧重新调用 colormap(gray)，该函数会重映射颜色表，在 getframe 前加一次 drawnow 即可，颜色表只在循环开始前设置一次。

5. 牛顿环动态仿真避坑指南：现象、原因、三行解决

5.1 中心黑斑变白斑：半波损失多写或少写一个 π

现象：仿真出的牛顿环中心是亮斑，跟教材里的暗斑完全相反。即使动态演示中环的扩散方向正确，这个错误依然很显眼。 原因：反射光强度公式用了 cos²(π r2/(Rλ))，而正确的应该是 sin²。cos² 在 r=0 处取 1，中心自然变亮。半波损失的本质是反射相位突变 π，公式里少处理这半个周期。 解决：把强度公式改回 sin(pi * r2 / (R * lambda)).^2。如果想用 cos² 写，就手动给相位加 π：cos(pi * r2 / (R * lambda) + pi/2).^2，但这属于自找麻烦，直接用 sin² 最省心。

5.2 满屏摩尔纹，几条环都数不出来：采样约束没算

现象：画面上全是细密的水波纹，或者画出来是几条歪歪扭扭的弧形，跟牛顿环完全没关系。图像分辨率明明很高，却看不出环形条纹。 原因：网格密度 N 低于 2L²/(Rλ) 的约束，条纹最密处的空间频率超过采样率，发生混叠。常见于把 L 从毫米改成厘米量级、或把 R 调小到 0.2 以下，却忘了同步调整 N。 解决：在参数区加一行守卫代码，采样不足时直接报错：

dp = L / N; % 像素对应的物理尺寸 dMin = R * lambda / L; % 模拟区域边缘的条纹间距 if dp > dMin / 2 error('采样不足，请将 N 至少提高到 %d', ceil(2 * L * L / (R * lambda))); end

5.3 GIF 卡顿或文件巨大：颜色表和帧尺寸是元凶

现象：动图文件动辄几十 MB，播放器打开就卡；或者 GIF 文件很小但画面全是色块，条纹边缘像锯齿。 原因：每帧都用 rgb2ind 生成新的颜色表，颜色表一变，GIF 需要为每一帧保存完整映射，压缩率断崖下跌。另一个原因是 getframe 抓的是整个 figure，如果窗口中有标题栏、坐标轴标签等多余元素，帧的尺寸会远超图像本身。 解决：循环外先设定固定的 256 色灰度表，循环内只做 append 写入；写 GIF 前用 figure('Position') 固定窗口尺寸，删掉不必要的装饰元素。实测同样参数下，固定颜色表能让文件体积缩到原来的三分之一左右。

5.4 环没向内收反而向外爆：扫描方向与相位符号

现象：R 从 0.5 扫到 2.0，按理环应该向外扩展、条纹变稀，但画面上环却在向中心收缩，而且越来越密。动态方向跟理论完全相反。 原因：相位表达式里 R 放错了位置。r_k = sqrt(kλR)，R 增大环半径必然增大。如果代码里误写成 phase = pi * r2 * R / lambda，效果就是 R 增大时相位增大，同一位置的条纹序数变大，看起来环向中心挤。 解决：检查 phase 中 R 是否在分母。同时用解析环位置做静态验证：先固定 R=0.5 和 R=2.0 各生成一张静态图，用 2.3 节的像素标定法数三条暗环的像素半径，确认它们分别对应 sqrt((k-1/2)λR) 的预测值，再去做动态扫描。

6. 进阶：用滑动条实时调参，并用像素行核对暗环半径

6.1 基于普通 figure 的交互式动态仿真

静态 GIF 适合展示，但自己调试时，每改一次 R 都要重跑脚本太浪费时间。我一般用一个加滑动条的交互窗口，滑一下 R 图像立刻更新。普通 figure 的 uicontrol 写法最稳，兼容性也最好：

fig = figure('MenuBar', 'none'); ax = axes('Units', 'normalized', 'Position', [0.1 0.15 0.8 0.8]); ph = imagesc(ax, x*1e3, x*1e3, I); axis(ax, 'image'); colormap(ax, 'gray'); sld = uicontrol('Style', 'slider', ... 'Min', 0.2, 'Max', 3.0, 'Value', 1.0, ... 'Units', 'normalized', 'Position', [0.1 0.05 0.8 0.04], ... 'Callback', @(src, ~) updateRing(src.Value, ph)); function updateRing(R, ph) I = sin(pi * r2 / (R * lambda)).^2; set(ph, 'CData', I); end

关键是提前拿到 imagesc 的句柄 ph，回调里只更新 CData，不重建坐标轴，这样滑动时响应无延迟。R 的滑动范围设在 0.2 到 3.0，足以覆盖大多数平凸透镜的典型曲率。这个交互窗口比 GIF 更适合课堂演示，学生拖一下滑块就能直观看到环疏密随 R 的变化。

6.2 自检方法：取一行灰度，三点拟合暗环半径

最后一个自检技巧：从生成的强度矩阵中取中心行，用 findpeaks 定位暗环的像素位置，再跟理论公式比对。这个方法每次跑动态仿真前我都会先跑一遍，确保公式没写歪。

row = I(round(N/2), :); % 取中心行的强度分布 [pks, locs] = findpeaks(-row); % 找极小值，即暗环位置 rPix = abs(x(locs)); % 暗环的物理半径 rTheo = sqrt((0.5:numel(rPix))' * lambda * R); % 理论暗环半径 table(rPix(:), rTheo, 'VariableNames', {'measured', 'theoretical'})

findpeaks 对边界假峰会误报，所以我在 rPix 里加一个条件：只保留 rPix 大于 3*dp 且小于 L/2 的位置。实测中心行能数出十几个环，跟理论值误差基本在 1 个像素以内。这个习惯帮我避开了很多次玄学问题——图像看着正常、但数字对不上，最后发现是坐标某个地方少乘了 1e3。养成任何动态仿真都带一个解析公式校验线的习惯后，这类问题基本能在五分钟内定位。希望帮到你。

本文还有配套的精品资源，点击获取