MATLAB波动光学仿真：干涉、衍射与极化模拟全解析

简介本资源是一套基于MATLAB实现的电磁波基础物理现象仿真工具包面向计算机、电子信息工程及应用数学等专业的本科生服务于课程设计、期末大作业与毕业设计等实践教学场景帮助学习者直观理解干涉、衍射与极化三大核心概念。压缩包共14个文件6个主功能M脚本、4幅结果示意图JPG、2个参数配置TXT、1个Git属性文件及1个说明文档MD总大小仅138KB轻量易部署。所有代码采用参数化设计波长、相位差、孔径尺寸、偏振角等关键物理量均可一键修改配合详尽中文注释与清晰模块划分如Interference.m、Pattern1.m、Shifted_Pattern.m等显著降低复现门槛与调试难度。附赠可直接运行的案例数据与完整实验流程说明支持MATLAB 2014a至2024b多版本兼容兼顾教学稳定性与技术前瞻性。1. 整体设计与物理思路拆解1.1 为什么拿 MATLAB 做波动光学模拟干涉、衍射、极化这三个概念是物理光学里最基础也是最重要的内容。大学课堂上一讲就懂但真到做实验、写论文的时候很多人对“相位差是怎么变成明暗条纹的”“偏振片转90度为什么就全黑了”这些问题还是一头雾水。原因很简单教科书上的公式是静态的而波动的本质是动态的演化过程脑子转不过来。MATLAB 做这件事有两个天然优势一是矩阵运算是它吃饭的本领而光场的传播本质就是对二维复数矩阵做变换语言特性和物理模型完全对得上二是可视化工具成熟imagesc、surf、animatedline、quiver都能直接画不用像 C 那样还得折腾第三方绘图库。我教学生的时候常说MATLAB 就是物理人的“显微镜”——不是因为它画图好看而是因为它能把公式变成眼睛能看见的图。1.2 三个现象的内在逻辑线这个模拟项目最忌讳的做法是三个模型彼此孤立地写三段代码。表面看是三个实验但底层是贯通的光是一种电磁波可以用空间中的电场矢量场描述。干涉是两个同频率的波在空间叠加相位差决定振幅的大小分布衍射本质上是无数个次波源发出的球面波相干叠加只是把“两个源”换成了“连续分布的源”极化则是描述每个空间点上电场矢量随时间振荡的轨迹形态。所以在设计代码时我用了一条统一的思路主线先建立空间离散网格用二维矩阵存光场。干涉和衍射都只算标量振幅或者准确说标量电场复振幅核心是相位。极化单独走一条线因为它处理的是矢量振荡需要三维空间加时间轴。这样安排的好处是读者跑完干涉代码之后再去看衍射代码会立刻发现结构高度相似学第二遍的成本极低。而且代码可以做成一个 GUI 面板的三个页签切换实验时只需换输入参数和物理模型主框架完全复用。1.3 方案选型为什么不用 COMSOL 或者 Ansys有学生问过我既然要做电磁波模拟为什么不用 COMSOL 或者 HFSS这些商用有限元软件确实精度更高能处理任意复杂结构和边界条件甚至能算出天线辐射方向图这种工程结果。但它们的学习曲线陡峭而且对“波动光学的基本原理演示”这个需求来说完全是杀鸡用牛刀。商用软件的建模思路是“划分网格 求解偏微分方程”它输出的是数值解但你很难直观看到“相位差从0变到π光强从亮变暗”这个连续过程——中间过程全被求解器封装了。而 MATLAB 做模拟是“从公式出发直接计算物理量的离散近似值”每一步都知道自己在算什么非常适合教学过程和物理概念的验证。尤其值得强调的是MATLAB 的 FFT 是官方高度优化的对 4096×4096 的矩阵做二维傅里叶变换在普通台式机上也就一两秒。用fft2模拟夫琅禾费衍射速度快到可以实时调节缝宽和缝间距这是有限元软件根本做不到的交互体验。2. 干涉模拟从两列波叠加到条纹分布2.1 物理模型与关键公式杨氏双缝干涉是所有干涉现象的原型。两列相干光波分别从缝 S1 和 S2 出发传播到观察屏上某一点 P 时走过的路程不一样产生了光程差进而转化为相位差。数学上相干叠加用复数是最简洁的假设两束光在 P 点的复振幅分别是E1 A * exp(1i * k * r1)和E2 A * exp(1i * k * r2)合振幅就是两者之和光强正比于合振幅模的平方。如果观察屏足够远远场条件r1 和 r2 的差值可以做近似r2 - r1 ≈ d * x / L其中 d 是双缝间距x 是屏上点的横向坐标L 是缝到屏的距离。这样的话相位差 δ (2π/λ) * (r2 - r1)光强分布直接写成I(x) 4 * A² * cos²(δ/2)2.2 核心代码实现与参数选择在实际写代码时我不建议用上面的远场近似式子一步到位因为那样虽然能出条纹图但掩盖了物理过程。更好的做法是直接算每束波到达屏上每个点的真实距离再叠加。% 模拟参数 lambda 632.8e-9; % 氦氖激光波长单位米 d 1e-3; % 双缝间距单位米 L 1.0; % 缝到观察屏的距离单位米 A 1; % 振幅归一化 % 屏幕坐标观察屏上的横向范围取条纹能清晰可见的区间 N 2000; % 采样点数 x linspace(-0.02, 0.02, N); % 屏幕宽度范围单位米 [X, ~] meshgrid(x, 1); % 单行即可后面画二维图再扩展 % 两缝位置缝距中心对称 y1 d/2; y2 -d/2; % 从各缝到屏上每点的距离这里把坐标系简化为一维 r1 sqrt(L^2 (x - y1).^2); r2 sqrt(L^2 (x - y2).^2); % 复振幅叠加 E1 A * exp(1i * 2*pi/lambda * r1); E2 A * exp(1i * 2*pi/lambda * r2); E_total E1 E2; % 光强归一化 I abs(E_total).^2; I I / max(I); % 绘图 figure(Color, white); plot(x*1000, I, b-, LineWidth, 1.5); xlabel(屏幕位置 x (mm)); ylabel(归一化光强 I/I_0); title(杨氏双缝干涉条纹光强分布); grid on;这段代码里exp(1i * k * r)这种写法是关键电磁波的传播相位是kr - ωt我们关心的是空间稳态分布所以时间因子exp(-iωt)可以省掉。复数表示的好处是叠加时幅度和相位同时处理不用手动算三角函数和差化积。2.3 条纹间距的验证与误差分析跑完代码之后建议做一件事量一量图上相邻亮纹的间距跟理论值对照。由双缝干涉公式条纹间距为Δx λL / d代入我们的参数632.8e-9 * 1.0 / 1e-3 ≈ 0.6328 mm。如果在图上数十个亮纹间隔然后取平均结果和 0.6328 mm 的偏差应该小于 1%。如果偏差大先查采样点 N 是否足够——x 方向范围 0.04 米如果 N2000每个像素约 0.02 mm能分辨 0.63 mm 的条纹没问题但如果你把 x 的范围扩到 ±0.1 米而 N 不变每个像素变成 0.1 mm条纹就糊了。这就是后面第四部分要细说的采样率问题。还有一个容易踩的坑如果直接用I 4 * A² * cos(δ/2).^2这种解析式画出来的图非常完美但它只能画二维剖面。而用上面的数值叠加法你可以把x改成二维网格[X, Y]一次性得到整个屏上的干涉图样效果和实验照片一样。3. 衍射模拟从单缝到二维口径与 fft23.1 单缝衍射连续次波源的叠加单缝衍射可以这样理解缝的宽度 a 内每一点都发出一个次波这些次波传播到同一块屏幕上某一点时相位各不相同叠加结果不是简单相加决定了光强。这跟双缝的区别是双缝只有两个源单缝有无穷多个源。但是不能真的写一个 for 循环把无穷个源全部遍历一遍——那样代码慢到没法用。正确做法是用解析公式Fraunhofer 单缝衍射的光强分布为I(θ) I₀ * (sin β / β)² β π a sinθ / λ代码写起来很短lambda 632.8e-9; a 0.1e-3; % 缝宽 0.1 mm L 1.0; x linspace(-0.05, 0.05, 4000); theta atan(x / L); % 衍射角 beta pi * a * sin(theta) / lambda; I_sinc (sin(beta) ./ beta).^2; I I_sinc / max(I_sinc); figure(Color, white); plot(x*1000, I, r-, LineWidth, 1.5); xlabel(屏幕位置 x (mm)); ylabel(归一化光强); title(单缝夫琅禾费衍射光强分布); grid on;这里sin(beta)./beta要注意 beta 趋近于 0 的位置即屏幕中心MATLAB 会算出 NaN。解决办法是这行的分母加一个 eps或者干脆用sinc(beta/pi)——MATLAB 自带的sinc函数已经处理过奇点稳稳的。3.2 二维圆孔衍射用 besselj 与 fft2 两种方案对照单缝是一维的问题但实际实验里更常见的是圆孔衍射比如激光打在小孔上产生的艾里斑。圆孔衍射的理论光强是I(θ) I₀ * (2*J₁(ka) / (ka))²其中 k 2π/λa 为圆孔半径MATLAB 提供了第一类贝塞尔函数besselj直接就能算lambda 632.8e-9; a 0.05e-3; % 圆孔半径 0.05 mm L 1.0; N 1000; rho_max 0.02; rho linspace(0, rho_max, N); theta rho / L; ka 2*pi/lambda * a * theta; % 这是无量纲量 % 注意这里的 ka 其实是 ka sinθ 的简写更严格地写是 k*a*sinθ I_airy (2 * besselj(1, ka) ./ ka).^2; I_airy(ka 0) 1; % 处理奇点 I_airy I_airy / max(I_airy); % 转成二维图样将一维径向分布旋转成二维 [X, Y] meshgrid(linspace(-rho_max, rho_max, N)); R sqrt(X.^2 Y.^2); theta2 R / L; ka2 2*pi/lambda * a * theta2; I2D (2 * besselj(1, ka2) ./ ka2).^2; I2D(ka2 0) 1; I2D I2D / max(I2D); figure(Color, white); imagesc(linspace(-rho_max, rho_max, N), linspace(-rho_max, rho_max, N), I2D); axis image; colormap(gray); colorbar; title(圆孔夫琅禾费衍射艾里斑); xlabel(x (m)); ylabel(y (m));这个方法的优点是快理论精确。但它的局限也很明显只能算圆形口径。如果口径是矩形、三角形、双缝、甚至一个五角星形状解析公式就不存在了。这时候要用通用数值方法。通用的方法就是傅里叶变换法。夫琅禾费衍射的物理本质是衍射屏后方的远场复振幅分布等于孔径平面上复振幅分布也叫透射函数的二维傅里叶变换。这句话翻译成 MATLAB 就是经典三行mask double(sqrt(X.^2 Y.^2) a); % 圆形孔径掩膜内部为 1外部为 0 E_far fftshift(fft2(mask)); % 二维 FFT 且把零频移到中心 I_fft abs(E_far).^2;用这种fft2方法算出来的艾里斑和用besselj解析式算出来的光强分布第一暗环的半径位置应该完全一致差别只在强度曲线的细微起伏上数值离散误差。我给学生演示时会同时画出两张图并排对比让他们意识到数值方法不是“不能用”而是要在理解了原理之后才敢放心用。3.3 衍射模拟的网格分辨率与补零技巧很多人直接用fft2(mask)跑出来一团糊抱怨说“和理论一点不像”十有八九是网格分辨率没调好。这里有一个关键经验网格的离散间距 Δx 决定了可分辨的空间频率范围。孔径平面尺寸 D比如缝宽 1 mm和波长 λ 决定衍射发散的角度尺度 θ ≈ λ / D。用fft2模拟时输出矩阵的大小和输入一样所以采样点不足时就只有中心几个像素是亮的周围全是黑的。解决办法是补零zero padding。不要把 mask 直接扔进fft2先在周围补大片零值等效于在物理空间上把计算区域扩大但孔径尺寸不变这样可以显著提升衍射图样的解析度。[Ny, Nx] size(mask); mask_padded zeros(4*Ny, 4*Nx); mask_padded(Ny:3*Ny-1, Nx:3*Nx-1) mask; E_far fftshift(fft2(mask_padded));补零之后衍射图样会细腻很多。原理很简单FFT 是 DFT 的快速算法DFT 的频率采样间隔是 2π/NN 变大了频率采样就更密插值出来的衍射图自然更平滑。4. 极化模拟从琼斯矢量到动态电矢量轨迹4.1 极化的数学描述千万不要用实数硬算干涉和衍射处理的是标量场极化完全不同。极化描述的是电磁波电场矢量在垂直于传播方向的平面内的振荡模式。一列沿 z 方向传播的平面波电场有两个正交分量Ex E0x * cos(ωt - kz φx) Ey E0y * cos(ωt - kz φy)二者的振幅比和相位差决定了极化状态φx φy线极化电场矢量在一条直线上来回振荡。φx - φy ±π/2 且 E0x E0y圆极化电场矢量末端画圆。其他情况椭圆极化。用 MATLAB 模拟极化最经典的做法是画“电矢量轨迹图”取一个固定位置z0改变时间 t电场矢量的末端在一个周期内在 xy 平面上画出一条曲线。或者做动态图让矢量随时间旋转起来。我强烈建议用复数琼斯矢量来算而不是直接对 cos 函数求值再画图。原因如果直接用cos(omega*t)得到的是实数结果相位差、初始相位处理起来非常繁琐代码里到处是三角函数的加减容易出错。而复数法写成矩阵乘法加入波片延迟器时就是一个旋转矩阵乘琼斯矢量的形式概念清晰也方便以后扩展到任意偏振光学系统。4.2 三类极化的 MATLAB 绘制实现下面的代码一次画出三种极化态的轨迹对比% 参数设置 omega 2*pi; % 角频率归一化 T 1; % 一个周期 t linspace(0, T, 200); phi 0; % 初相 % 三种极化态: [E0x, E0y, 相位差] pols struct(name, {线极化, 圆极化, 椭圆极化}, ... Ex0, {1, 1, 1.2}, ... Ey0, {1, 1, 0.8}, ... dphi, {0, pi/2, pi/4}); figure(Color, white, Position, [100 100 1400 450]); for i 1:3 Ex pols(i).Ex0 * cos(omega*t phi); Ey pols(i).Ey0 * cos(omega*t phi pols(i).dphi); subplot(1, 3, i); plot(Ex, Ey, b-, LineWidth, 2); hold on; % 画一个单位圆做参照 th linspace(0, 2*pi, 100); plot(cos(th), sin(th), k--, LineWidth, 0.5); % 画坐标轴 line([-1.5 1.5], [0 0], Color, k, LineWidth, 0.5); line([0 0], [-1.5 1.5], Color, k, LineWidth, 0.5); axis equal; grid on; xlim([-1.5 1.5]); ylim([-1.5 1.5]); xlabel(E_x); ylabel(E_y); title(pols(i).name); end这段代码画出的静态轨迹图已经足够说明极化的分类。但如果你要做的项目是个可交互演示动态效果更震撼。动态实现用animatedline记录矢量端点轨迹同时用quiver画当前时刻的电场矢量箭头figure(Color, white, Position, [200 200 600 500]); ax axes; hold on; axis equal; grid on; xlim([-1.5 1.5]); ylim([-1.5 1.5]); xlabel(E_x); ylabel(E_y); title(动态圆极化电矢量旋转); % 轨迹线和箭头 trail animatedline(Color, b, LineWidth, 2); for n 1:length(t) Ex 1 * cos(omega*t(n)); Ey 1 * cos(omega*t(n) pi/2); addpoints(trail, Ex, Ey); % 删除之前的箭头重新画 if n 1, delete(qh); end qh quiver(0, 0, Ex, Ey, r, LineWidth, 2, MaxHeadSize, 0.5); drawnow limitrate; pause(0.03); end4.3 从二维轨迹到三维螺旋波极化的全貌二维轨迹图只能看到“投影”看不到传播方向上的振荡。如果想把电磁波沿 z 轴传播时电场矢量在三维空间的真实形态也表示出来可以将不同 z 位置处的电场矢量画成三维螺旋。比如圆极化波沿 z 传播时电场矢量末端在空间划出一条螺旋线。实现方式并不复杂固定一个时刻 t₀改变 z 坐标计算该位置的电场分量再把三维曲线用plot3画出来z linspace(0, 2*lambda, 500); t0 0.2 * T; Ex_z E0 * cos(omega*t0 - 2*pi/lambda * z); Ey_z E0 * cos(omega*t0 - 2*pi/lambda * z pi/2); figure(Color, white); plot3(z/lambda, Ex_z, Ey_z, b-, LineWidth, 1.8); xlabel(z / \lambda); ylabel(E_x); zlabel(E_y); title(圆极化波电场矢量三维轨迹固定时刻); grid on; view(30, 25);这张图非常直观地展示了“旋转”到底是怎么发生的它不是振动方向在转而是相位沿传播方向线性变化导致电场矢量的方向在空间中旋转前进。看到这张图的人基本不会再对“圆极化为什么角动量是 ±ℏ”这种问题困惑。5. 常见问题与调试技巧实录5.1 条纹不清晰或位置不对的排查流程这是做干涉模拟最高频的问题。我总结了一套排查顺序建议新手照着走现象排查方向条纹根本没有检查两列波是否用了同一个频率相位差是否为 0 导致完全相干但零相位差网格范围是否太小条纹有但模糊采样点数 N 不够间距和波长比值过大导致条纹细于一个像素条纹位置和公式算的不一致检查是否用了distance sqrt(L^2 (x-y)^2)而不是近似公式如果用近似公式检查远场条件 L d²/λ 是否成立中心不出现在 x0检查布尔运算是否有 bug或坐标轴方向搞反了图形边缘出现不需要的振荡这是有限孔径效应不是因为物理模型错误把范围缩小即可我记得有个学生用近似公式做L 只取了 0.2 md 是 5 mm套Δx λL/d算出 0.025 mm但图上量出来 0.028 mm。误差来自光程差的一次项近似不成立——这种情况必须回退到精确距离计算不能再近似。5.2 fft2 结果中心不亮或四角亮的问题很多人在用fft2模拟衍射时得到的光强分布中心是黑的四角反而是亮的。这是坐标习惯问题MATLAB 的fft2默认把零频放在矩阵的四个角上而不是中心。解决方法是加一行fftshiftE_far fftshift(fft2(mask_padded));但要注意fftshift只是把矩阵平移重组不改变数值。加完之后光强矩阵中心的点对应零频也就是衍射零级中心最亮向外逐渐变暗一圈圈暗环对应艾里斑的暗环。如果没加fftshift你会在四角看到四个亮点中心一团黑完全没法用。另外一个细节fft2输出的矩阵尺寸很大如果直接imagesc(abs(E_far).^2)光强动态范围极大中心峰值可能比边缘高 6 个数量级图上看就是中心一个白点周围全黑。一定要做归一化或者用log(1 I)取对数再显示I_log log(1 abs(E_far).^2); imagesc(I_log); colormap(gray); axis image;用对数显示后暗环才可见。这是所有光学衍射仿真里最重要的一条显示技巧比调任何算法参数都管用。5.3 极化模拟中“矢量长度不变”的常见疑问刚接触极化的人容易陷入一个误区认为圆极化波的电场矢量长度必须恒定所以三维轨迹应该是一个完美的圆。对但也不全对——这取决于观察时刻。如果固定一个时刻 t₀ 看 z 方向圆极化波的电场分量分别是 cos 和 sin平方和恒为 1确实是一条半径恒定的螺旋线。但如果你不是在“同一时刻”看而是在“同一 z 位置”随时间看那么电场矢量的末端在 xy 平面上画的是一个圆矢量长度当然也是恒定的。很多模拟代码画出“变半径”的螺旋原因是采样 z 范围和相位kz对应的2π/lambda * z没有对齐导致看起来半径在跳动。解决方法是确保 z 轴的范围覆盖整数个波长比如z linspace(0, 2*lambda, 500)这样在 z 从 0 到 2λ 的画面上螺旋线恰好转了整数圈半径平滑。5.4 性能优化矩阵运算代替循环最后说一个 MATLAB 老生常谈的问题循环慢。衍射模拟中如果要扫描参数比如扫缝宽从 0.05 mm 到 0.5 mm 共 100 个值用 for 循环里套fft2也不是不行但算起来明显卡。优化手段优先级如下尽量向量化避免对每个像素单独计算。能对矩阵整体操作时不写 for。比如计算距离矩阵用meshgrid生成坐标网格后直接sqrt(X.^2 Y.^2)而不是两层循环逐个点算。多组参数扫描时用parfor并行。如果电脑有 8 个物理核默认用parfor能把 100 组扫描时间从 100 秒压到 20 秒左右。前提是内存够而且注意parfor内不能依赖循环顺序。动态图每帧都刷新全图会卡用drawnow limitrate或者set更新已有图形对象的XData/YData不要无限plot新句柄。5.5 一套统一参数配置与界面封装建议如果这个项目要交给学生用或者放入课程演示建议不要让他们改代码里分散的参数单独做一个配置函数把波长、缝宽、透镜焦距、屏幕距离全部集中管理function params getSimParams() params.lambda 632.8e-9; % 单位米 params.d 1e-3; % 双缝间距 params.a 0.1e-3; % 单缝缝宽 params.L 1.0; % 传播距离 params.N 2048; % 网格点数 params.screenSize 0.04; % 屏幕范围半宽 params.polType circular;% 极化类型 end封装后只需改参数表三个实验共用一套网格生成和绘图工具函数。结语这个项目还能怎么扩展做完基础版之后如果你觉得不过瘾以下几个方向我亲测值得尝试把双缝改成多缝光栅看看主极大、次极大和缺级现象把单缝口径换成矩形或异形孔径用fw和fft2对比衍射图样在极化部分加入 1/4 波片和 1/2 波片的琼斯矩阵模拟从线极化到圆极化的转换过程甚至可以把三个部分组合成一个完整的“偏振光双缝干涉”仿真——用 two orthogonal polarization 分量分别干涉再合成输出偏振态这是许多光学实验教材里会做但现实中很难调出来的实验在 MATLAB 里却只需十几行代码。我在做这个项目时最大的体会是真正把你的理解拉开差距的不是会背公式而是能在代码里把公式“翻译”成人话。如果你做完这个模拟能对着图给别人讲清楚“为什么光栅的条纹比双缝细那么多”“为什么圆偏振光经过偏振片后强度减半”那说明你已经不是停留在公式层面的初学者了。把这份代码自己重写一遍跑通再改参数看看现象变化这个项目的价值才算真正发挥出来。本文还有配套的精品资源点击获取