MATLAB激光器谐振腔模拟：ABCD矩阵、稳区与Fox-Li迭代

简介：面向激光物理、光电信息专业学生及科研人员，这份MATLAB资源围绕激光器谐振腔的模拟分析展开，覆盖从物理建模、参数设置到数值求解的完整流程，适合课程设计、科研入门或项目预研时快速上手。压缩包共5个文件，包含4个MATLAB脚本和1个TXT说明文档，整体仅9KB，轻量而聚焦，便于逐行阅读和调试。脚本分别实现谐振腔模式函数计算、高功率激光器稳定区分析、腔内多模传输求解以及模式尺寸随腔长变化等功能；说明文档则对平行平面腔自再现模式的模拟要点进行了梳理，可与脚本配合理解理论到代码的映射。已有354人学习下载。通过运行代码，读者可直观获得光场分布、稳定区边界等关键结果，并掌握基于传输矩阵法的数值迭代思路，为后续设计高功率激光器、优化谐振腔结构提供实用参考。

1. 激光器谐振腔模拟，为什么MATLAB比光学追迹软件更适合练手

拿到一套谐振腔仿真需求，多数人第一反应是打开现成光学软件把腔型摆进去等结果。这套流程确实高效，但整套分析里最关键的稳定性边界、束腰位置和模式损耗全都藏在操作界面背后，一旦腔型超出常见范围，计算结果对不上实测，你根本不知道是该改参数还是改软件设置。用MATLAB从矩阵开始写谐振腔模拟，等于把黑匣子拆开重装一遍：g参数稳区、往返矩阵、Fox-Li迭代、q参数传输，每一步都落到数组和公式上。这篇笔记面向激光器设计和光学实验方向的从业者，用一套完整的平凹腔模拟案例讲清楚谐振腔仿真怎么做、参数怎么设、哪些位置最容易翻车。

2. 把谐振腔写成ABCD矩阵：g参数稳定判据与可执行代码

2.1 谐振腔的元器件矩阵分解

谐振腔仿真的起点不是渲染腔体结构，而是把每一个光学元件表达成2×2的ABCD矩阵。光线在腔内往返一周，依次经过自由空间传播、镜面反射、再次自由空间传播，整个过程就是一个矩阵链。对平凹腔这种最常见腔型，三个基本矩阵就够了。

自由空间传播距离L，矩阵形式为：

function M = fresnel_prop(L) M = [1, L; 0, 1]; end

曲率半径为R的反射镜，近轴矩阵为：

function M = mirror_reflect(R) M = [1, 0; -2/R, 1]; end

注意两个细节。第一个是反射镜矩阵的写法对凹面镜和凸面镜完全不一样，R为正代表凹面镜，R为负代表凸面镜，符号写反，后续所有稳定性和光斑结果全部颠倒。第二个是平凹腔的平面镜R取Inf，在MATLAB里直接传Inf计算即可，-2/Inf得到0，矩阵退化为单位阵。

把三个矩阵按光线传播顺序相乘，就得到从平面镜出发、经凹面镜反射、再回到平面镜的往返矩阵：

R2 = 5.0; % 凹面镜曲率半径，单位m L = 1.0; % 腔长，单位m M_round = fresnel_prop(L) * mirror_reflect(R2) * fresnel_prop(L);

矩阵乘法顺序代表光线作用顺序：第一个作用在光线上的矩阵放在最右边。很多初学者习惯从左往右读，结果把传播方向搞反了，算出来的光斑尺寸与实际差出几倍。我一般会在脚本开头用注释把顺序写清楚：“从平面镜出射→传播L→凹面镜反射→再传播L”。

这个往返矩阵是整个谐振腔数值分析的骨架。后续无论是算g参数、用q参数求束腰，还是做Fox-Li迭代，都要从它出发。矩阵写对，后边所有计算才有意义。

2.2 g参数稳定性条件与稳区图

有了往返矩阵，稳定性判据可以直接从矩阵元素读出，也可以回到经典g参数定义。对腔长L、两镜曲率半径R1和R2，g参数定义为：

g1 = 1 - L/R1 g2 = 1 - L/R2

稳定腔条件为0 < g1×g2 < 1，等号对应临界腔——共焦腔或平行平面腔正落在边界上，实际工程中极少把工作点设计在边界附近，因为微小装配误差就会让g乘积越界，模式损耗急剧增大。

在MATLAB里画稳区图，直接生成g1-g2平面的等值线：

g1 = linspace(-2, 2, 401); g2 = linspace(-2, 2, 401); [G1, G2] = meshgrid(g1, g2); prod = G1 .* G2; stable = (prod > 0) & (prod < 1); imagesc(g1, g2, stable); axis xy; colormap(gray); xlabel('g1 = 1 - L/R1'); ylabel('g2 = 1 - L/R2'); title('稳定腔区域图');

linspace(-2, 2, 401)把g参数范围设得比实际腔型域更大，方便观察稳区边界；meshgrid生成二维网格后，prod > 0 & prod < 1是一次性做双重条件判断，得到的就是稳定腔区域掩模。这里用imagesc而不是contour，是因为稳定/非稳定是二值判断，用图像显示边界最直观。矩形硬边界看起来不光滑，把网格点加密到401×401后，边界位置的视觉误差已经小于1%。

画完图，再把具体腔型的g工作点叠加进去：

R1 = Inf; % 平面镜 g1_val = 1 - L/R1; % g1 = 1 g2_val = 1 - L/R2; % g2 = 1 - 1/5 = 0.8 hold on; plot(g1_val, g2_val, 'ro', 'MarkerSize', 8, 'LineWidth', 2);

这组参数计算出来g1×g2=0.8，落在稳定区域内部且离边界有足够余量，是典型的工作点设计。平凹腔里g1恒为1，所以工作点是一条竖线，腔长越接近曲率半径R2，g2越接近0，工作点越靠近稳区边界。

参数选型上需要解释一下：为什么选择R2=5m、L=1m而不是其他值？因为这个腔型g1g2离边界还有20%余量，镜片加工容差和热畸变导致的曲率变化不至于让腔直接进入非稳区。实际设计时我建议至少留15%的稳区余量，否则激光器装调时会反复出现不出光的情况，每次都要怀疑镜子装歪了，其实根源是参数选得太贴边界。

2.3 非稳腔为何值得留意

稳定腔条件解决的是低阶模低损耗问题，但很多高功率激光器用的是非稳腔。非稳腔的g1g2乘积落在0到1区间之外，基模损耗大，却有很好的横模鉴别能力，输出耦合可以做得非常高，适合大增益体积、大模体积的振荡器设计。

正向非稳腔的几何放大率由往返矩阵的迹决定：

m = (abs(M_round(1,1) + M_round(2,2)) + sqrt((M_round(1,1) + M_round(2,2))^2 - 4)) / 2;

放大率m可以粗略理解为每一往返后光斑横向扩大的倍数。这个参数对非稳腔设计很关键——输出耦合效率和镜面尺寸都要按放大率推算。不过在入门阶段，先集中把稳定腔的模拟做扎实，非稳腔的衍射迭代计算需要更谨慎的网格采样，直接套用稳定腔的代码会出很多意想不到的数值问题。

3. Fox-Li迭代求解腔模：衍射积分离散化与收敛判据

3.1 自再现模原理与角谱传播

g参数只能回答“这个腔稳不稳定”，回答不了“腔内光斑长什么样”。要得到基模横向分布和衍射损耗，必须做本征模式求解。经典Fox-Li迭代的核心是自再现思想：光在腔内往返足够多次后，横向场分布趋于稳定，每一往返只改变一个复常数增益因子，分布形状不再变化。

实现这一思想，通常借助角谱衍射传播。把镜面M1上的场分布做二维傅里叶变换，在频域乘以近轴传播传递函数，再逆变换，就得到传播L后在镜面M2上的场分布。单次往返过程包含两次自由空间传播和两次镜面反射。

3.2 谐振腔单程传播的MATLAB实现

仿真参数和网格先定义好：

lambda = 1064e-9; % Nd:YAG激光波长，单位m L = 1.0; % 腔长，单位m R2 = 5.0; % 凹面镜曲率半径，单位m D = 6e-3; % 凹面镜有效孔径直径，单位m N = 256; % 采样网格点数，方形网格 dx = 0.2e-3; % 空间采样步长，单位m x = (-N/2 : N/2-1) * dx; [X, Y] = meshgrid(x, x); r2 = X.^2 + Y.^2; % 径向坐标平方，单位m^2 % 频域坐标 fx = (-N/2 : N/2-1) / (N * dx); [FX, FY] = meshgrid(fx, fx); f2 = FX.^2 + FY.^2; % 凹面镜反射相位因子（薄透镜近似） phase_mirror = exp(-1i * pi / (lambda * R2) * r2); % 凹面镜孔径遮挡 aperture = zeros(N, N); aperture(r2 <= (D/2)^2) = 1;

这里有几个参数必须提前想清楚。N=256是采样点数，点数太少频谱分辨率低，相位调制细节丢失；dx=0.2mm决定空间网格步长，窗口总宽度是256×0.2mm=51.2mm，大约是孔径直径6mm的8.5倍，既能完整保留边缘衍射，又不会因为窗口过大浪费采样点。频域坐标fx的范围由1/(N*dx)决定，最大空间频率约1/0.2mm=5000m⁻¹，对应衍射角约0.005弧度，对5m曲率半径镜面产生的相位变化来说完全够用。

初始场选择高斯分布：

w0 = 1.5e-3; % 初始光斑半径猜测值，单位m u = exp(-r2 / w0^2); u = u / sqrt(sum(abs(u(:)).^2) * dx^2); % 能量归一化

初始光斑半径的猜测值允许有误差，迭代会逐渐收敛到自再现模。需要说明的是：如果初始光斑给得和真实基模偏差太大，比如给了平面波，迭代次数会明显增加，甚至需要几千次才收敛；如果给得太窄，高频分量太多，可能激发高阶模。用高斯分布做初值是最稳妥的常见做法。

接下来是迭代主循环：

max_iter = 2000; tol = 1e-6; last_field = u; loss_history = zeros(max_iter, 1); for iter = 1:max_iter % 从平面镜传播到凹面镜：角谱法 U = fftshift(fft2(u)); k = 2 * pi / lambda; H = exp(1i * k * L) .* exp(-1i * pi * lambda * L * f2); u2 = ifft2(ifftshift(U .* H)); % 凹面镜反射：施加曲率相位与孔径遮挡 u2 = u2 .* phase_mirror; u2 = u2 .* aperture; % 从凹面镜传播回平面镜，同样用角谱法 U2 = fftshift(fft2(u2)); u1 = ifft2(ifftshift(U2 .* H)); % 归一化，记录能量损耗 energy_before = sum(abs(u1(:)).^2) * dx^2; u1 = u1 / sqrt(energy_before); loss_per_round = 1 - energy_before; % 收敛判断：归一化场分布的变化量 diff = sum(abs(u1(:) - last_field(:)).^2) * dx^2; last_field = u1; loss_history(iter) = loss_per_round; u = u1; if diff < tol fprintf('迭代收敛于第%d次\n', iter); break; end end

这段代码的逻辑分四步：正向传播、凹面镜作用、反向传播、归一化与收敛判断。正向传播和反向传播用的是同一个传递函数H，因为两次传播距离都等于腔长L，腔结构对称时可以直接复用。凹面镜作用分两步实施，先乘相位因子反映曲率对波前的弯折，再乘孔径掩膜反映镜面有限尺寸造成的衍射损耗，两者物理意义不同，不能合并成一个矩阵。

fftshift和ifftshift的配对容易写错——傅里叶变换前要把场分布从左下角排列挪到中心排列，逆变换后再挪回来。如果只挪一次不配对，结果会多出一个整体空间偏移，而且这个偏移很难靠肉眼从强度图上发现，通常到和解析光斑半径做对比时才暴露。

3.3 收敛后的模式损耗分析与物理解读

迭代结束后，我们从两个维度看结果：横向强度分布与衍射损耗。

intensity = abs(u).^2; intensity = intensity / max(intensity(:)); % 沿x轴截取分布，拟合光斑半径 profile = abs(u(N/2+1, :)).^2; w_fitted = dx * sqrt(2 * sum(profile .* (x - 0).^2) / sum(profile)); fprintf('Fox-Li迭代基模光斑半径：%.2f mm\n', w_fitted*1e3); % 稳态损耗 steady_loss = mean(loss_history(end-100:end)); fprintf('稳态单程损耗：%.4f%%\n', steady_loss*100);

光斑半径用二阶矩定义计算而不是直接取1/e²点，因为数值解分布不严格是高斯形，二阶矩能更客观地反映能量扩展。

这里有个很反直觉的现象：归一化后的场分布形状不再变化，但每一往返能量仍然有小幅损耗，这个损耗来自镜面孔径的硬边衍射，只要孔径有限，损耗就存在。若把孔径去掉，迭代能量基本守恒，场分布也会慢慢展宽且收敛不到稳定形态。硬边衍射是谐振腔模拟中少数“必须有”的耗散机制，工程上用来抑制高阶模，数值上则保证自再现方程有非零解。

4. 高斯光束q参数传输：束腰位置与光斑尺寸的快速计算

4.1 q参数、光斑尺寸与波前曲率的关系

Fox-Li迭代给出数值解，但实际工程里更常用复光束参数q做解析计算。q参数和光斑半径w、波前曲率半径R的关系写成：

1/q = 1/R - i·λ/(π·w²)

实部对应波前曲率，虚部对应光斑尺寸。这个表达方式的工程价值在于：q参数经ABCD矩阵传输后，用双线性变换更新：

q' = (A·q + B) / (C·q + D)

只要知道入射光束在某个参考面的q参数，任意位置的光斑尺寸、曲率半径都能用复数运算算出来，不需要逐点做衍射积分。数值上它比Fox-Li快几个量级，适合做参数扫描和实时反馈。

4.2 平凹腔自再现q参数求解

稳定腔内往返一周后q参数必须自再现，即q经过往返矩阵后回到原值。这给出一元二次方程，程序里直接求根：

R2 = 5.0; L = 1.0; M_round = fresnel_prop(L) * mirror_reflect(R2) * fresnel_prop(L); A = M_round(1,1); B = M_round(1,2); C = M_round(2,1); D = M_round(2,2); % 自再现方程 C*q^2 + (D-A)*q - B = 0 p = [C, D-A, -B]; roots_q = roots(p); % 筛选虚部为负的物理有效根 valid = roots_q(imag(roots_q) < 0); q0 = valid(1); w0 = sqrt(-lambda / (pi * imag(1/q0))); R_curv = 1 / real(1/q0); fprintf('基模束腰半径 w0 = %.2f mm\n', w0*1e3);

这段代码的关键在roots函数两项选择：一元二次方程有两个根，只有虚部为负的根对应物理光束。若不加筛选直接取第一个根，可能拿到虚部为正的发散解，光斑半径变成复数，后边所有物理量全是错的。这是一个非常典型的翻车点。

自再现束腰半径计算出来后，束腰位置也可以推算。平凹腔中凹面镜的曲率中心附近会形成束腰，严格位置从平面镜出发找q虚部为零的传输距离：

z = linspace(0, L, 1001); wz_list = zeros(size(z)); for i = 1:length(z) Mz = fresnel_prop(z(i)); qz = (Mz(1,1)*q0 + Mz(1,2)) / (Mz(2,1)*q0 + Mz(2,2)); wz_list(i) = sqrt(-lambda / (pi * imag(1/qz))); end [~, idx_min] = min(wz_list); z_waist = z(idx_min); fprintf('束腰距平面镜距离：%.3f m\n', z_waist);

linspace(0, L, 1001)把平面镜到凹面镜的整个腔长范围细分成1001个截面，逐点计算q参数再提取光斑半径。[~, idx_min] = min(wz_list)这行返回光斑最细处的位置。需要注意这个束腰位置搜索依赖于先前解出的q0，如果第2.2节往返矩阵顺序写反，q0本身就是镜像解，扫描出的束腰位置也会落在错误的端面。

4.3 q参数解析与Fox-Li迭代结果的交叉验证

解析公式和高阶数值迭代不应该彼此孤立。我在实际项目中习惯同时跑两组仿真，对比结果判断参数设得对不对：

q参数法得到的光斑半径是“无限大镜面”下的理想基模，Fox-Li迭代在加了6mm孔径后得到的半径通常会小一点，因为硬边衍射限制了光场扩展。如果两者差距在3%以内，说明网格和采样参数设置合理，数值结果可信；差距超过10%，优先检查网格步长和窗口宽度，再检查初始场是否收敛到了高阶模。

这个交叉验证习惯花费的时间不到五分钟，却能把Fox-Li迭代中绝大多数肉眼察觉不到的数值瑕疵暴露出来。我每次换新腔型都会强制跑一遍对比，很大程度上避免了后边实验调试时把“仿真误差”和“装配误差”搞混。

5. 谐振腔仿真避坑：采样窗口、初始场与硬边衍射的典型翻车点

5.1 采样窗口过窄导致模式能量被截断

现象：Fox-Li迭代收敛后，光斑分布边缘出现明显的矩形条纹，强度分布不光滑，拟合的半径偏小且和q参数解析结果差距越来越大。

原因：网格总宽度N×dx小于光斑实际扩展范围，场分布碰到计算窗口边缘被强制截断，等效于人为加了一个不该有的矩形光阑。

解决：先按q参数解析公式估算光斑半径，再设定窗口宽度为预估光斑半径的6到8倍。具体操作是：用当前代码跑一次，在收敛后统计光斑半径w_fitted，如果w_fitted大于窗口宽度的1/6，把dx调大或增加N重算，重新对比解析值。

5.2 网格步长过大导致频域传递函数采样不足

现象：迭代不收敛，损耗曲线振荡且幅度不下降。频率谱在高频端出现折叠特征，强度图上有精细的莫尔条纹。

原因：角谱传播中频域相位φ=πλLf²随f²增长，当dx过大，最高空间频率分量对应的相位变化超过π，频域欠采样造成相位混叠。

解决：检查频域坐标最大值的相位是否满足πλLf_max² < π/2。以λ=1064nm、L=1m为例，f_max必须小于约3.86×10⁴ m⁻¹，对应dx必须大于1/(2·f_max·N)约0.65μm——实际工程中取dx=0.2mm已经远大于这个下限，问题更多出现在dx设置过小导致f_max过大。这里容易搞反：dx过小同样有害。我踩过的坑是优先把dx往小调，结果N固定时窗口变小，反而触发5.1的问题。

5.3 初始场给随机噪声导致收敛缓慢甚至落入高阶模

现象：相同腔型参数，换成随机相位初始场后迭代两千次还没收敛，或者收敛后光斑分布带有一个明显的相位旋转，损耗值明显偏高。

原因：随机相位场包含大量高频分量和多个横模成分，Fox-Li迭代收敛的是损耗最低的横模，但随机初值不容易把能量集中到该模式上，迭代过程在模式竞争的“中间态”里停留很久。

解决：初始场用高斯分布，宽度取预估光斑半径的0.8到1.5倍。就算给宽一点，迭代前几十轮会快速收敛到基模邻域。想验证高阶模存在性时，可故意用厄米-高斯组合做初值，但要清楚此时得到的可能是多模叠加而不是纯基模，不能直接当作单模结果用。

5.4 忘记施加镜面孔径，迭代永远不收敛

现象：去掉aperture后，每轮损耗接近零，光斑半径不断展宽且没有稳定趋势，迭代差值diff始终在10⁻³量级不下来。

原因：无孔径谐振腔在近轴近似下无衍射损耗，自再现方程没有有限尺寸的非零解。能量不断向外扩展，等价于场分布永不闭合。这在物理上是“理想腔没有稳态解”，在数值上就是迭代发散。

解决：给两个镜面中至少一个设定实际通光孔径。如果真实系统里镜架、泵浦模块本身有明显限孔，可以把等效孔径直接设为机械限制中最窄的一处直径。这里给一个经验值：孔径直径取预估基模光斑直径的3到5倍，既保留明显的基模损耗优势，又不至于把损耗压得太大导致能量利用率下降。

5.5 矩阵乘法顺序写反，腔内往返方向颠倒

现象：q参数法解出的w0和Fox-Li迭代结果差得离谱，甚至算出的束腰负值或虚值。

原因：公式上往返矩阵写成M_round = fresnel_prop(L) * mirror_reflect(R2) * fresnel_prop(L)，但如果习惯性地把矩阵沿着“平面镜→凹面镜→平面镜”的物理顺序从左往右排，就变成了T(L) * R(R2) * T(L)的自右向左作用，方向反了，等价于交换了R1和R2，对非对称腔结果完全不同。

解决：给矩阵链加注释，标明靠近向量的矩阵是第一个作用在光线上的元件。验证方法很简单：把R1和R2交换位置，看两个镜面分别是无穷大平镜和5m凹镜时的结果是否交换。如果交换后结果一样，说明本来就结构对称，检验失效；对平凹腔这种不对称腔，交换后w0应当显著变化，若没变化说明矩阵顺序写错。

6. 从稳区图到参数扫描：验证仿真结果的三个进阶习惯

习惯一：做参数扫描前，先把单点自洽检查跑通。我见过不少人在参数扫描里得到一条光滑曲线，却忘了验证其中一个点是否物理正确。正确流程是先固定一组工作点，用q参数法算光斑半径和束腰位置，再让Fox-Li迭代收敛，比较两个结果差异在3%以内，然后才开始扫描。

习惯二：腔长扫描时注意稳区边界突变。以R2=5m、R1=Inf的平凹腔为例，腔长L从0.5m扫到4.9m，g2从0.9单调降到0.02，g1g2乘积从0.9降到0.02，束腰半径在L接近5m时发散。用一个循环把这段趋势画出来：

L_list = linspace(0.5, 4.9, 200); w0_list = zeros(size(L_list)); for i = 1:length(L_list) M = fresnel_prop(L_list(i)) * mirror_reflect(R2) * fresnel_prop(L_list(i)); A = M(1,1); B = M(1,2); C = M(2,1); D = M(2,2); rts = roots([C, D-A, -B]); rts = rts(imag(rts) < 0); if isempty(rts) w0_list(i) = NaN; else w0_list(i) = sqrt(-lambda / (pi * imag(1/rts(1)))); end end plot(L_list, w0_list*1e3); xlabel('腔长 L/m'); ylabel('束腰半径 w0/mm');

这段扫出来的曲线会在L接近R2时急剧抬升，曲线尽头w0发散。看到这种现象不要慌，这不是程序错误，而是谐振腔物理上在临界点失效的体现。用这个图选工作点，我会把L设在1m到4m之间的中段区域，这个范围既远离稳区边界，又有相对小的光斑，工程调试容差好。

习惯三：双参数稳区扫描，确认设计点不是孤点。把g1和g2同时扫描，画g1g2乘积的等高线，看目标工作点周围是否存在连续的稳定区域。只做单参数扫描容易漏掉二维空间里的窄带通道。双参数扫描用contourf实现，稳定区内填充绿色，不稳定区留白，目标工作点用星号标记，一眼就能看出余量方向。

从那以后，我每次搭建新的谐振腔模型都强制走一遍这套流程：g参数验稳区、q参数定束腰、Fox-Li验证横向分布、扫描确认设计余量，单点合格再做扫描。这套组合拳帮我躲过了很多数不清的翻车现场，尤其是那些仿真图上看起来挺漂亮、实际装调却根本出不了光的腔型。希望帮到你。

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