简介面向光子学研究者的波导仿真MATLAB代码采用有限差分光束传播法FD-BPM计算光场在波导中的传播适合波导、光栅及微结构光纤等结构的数值模拟与优化设计。压缩包包含1个m格式脚本整体大小约1KB为极简的平面波导FD-BPM实现便于快速部署与修改调试。脚本完整涵盖初始化、网格划分、光场离散、传播迭代、边界条件处理及结果分析等流程并支持通过调整波导几何尺寸、光源波长与初始场分布来考察模式分布、功率传输和损耗特性。已有347人学习下载对于从事光子器件课题研究或课程设计的学生与工程师这份轻量代码可直接运行参考省去从零编写算法的过程也能加深对BPM数值方法的理解。得益于FD-BPM的数值特性脚本在常见PC上即可高效运行适合教学演示与前期方案验证。1. 平面波导模拟为何绕不开 FD-BPM 这一步做光子集成的人大概都有过这种经历结构画好了、参数算好了一进工艺线出来的器件响应和设计值差出一大截。问题常常不出在工艺而出在模拟阶段——很多人直接拿 FDTD 当万能工具却忽略了在平面波导这种弱约束、长传播场景下全波方法计算量大得离谱几毫米的传播长度足够让普通工作站跑上一整天。FD-BPM 的价值正在于此它把麦克斯韦方程降维成单向传播的抛物型方程用有限差分在横截面离散在传播方向上步进求解几秒钟就能看到整个光场演化过程。MATLAB 里用 FD-BPM 研究波导不需要昂贵的商用光学软件也不需要集群资源一台笔记本就能完成从模场分析到传输损耗估算的全部工作。这篇要讲的FD_BPM_planar_waveguide.m正是典型实现适合刚接触光子学仿真的人快速跑通也适合老手用来做结构初筛和参数扫描。2. 从亥姆霍兹方程到 Crank-Nicolson 格式FD-BPM 的理论地基2.1 缓变包络近似与单向传播假设标量亥姆霍兹方程是光场传播的起点∇²E k₀²n²(x, z)E 0直接数值求解这个椭圆型方程需要同时处理前向和后向传播的场边界条件复杂且计算量极大。BPM 的核心做法是假设光场可以写成快变相位和缓变包络的乘积E(x, z) ψ(x, z) · exp(-jk₀n₀z)其中 n₀ 是参考折射率ψ 是缓变包络。把这个代入亥姆霍兹方程并忽略 ψ 对 z 的二阶导数项这就是缓变包络近似SVEA得到∂ψ/∂z (j/2k₀n₀) · [∂²ψ/∂x² k₀²(n²(x, z) - n₀²)ψ]这个方程是抛物型方程只允许前向传播所以可以用步进法从 z0 逐层推进到 zL。这个近似的物理意义就是要求波导结构沿着传播方向变化足够缓慢折射率对比度也不宜太大——大部分平面波导和脊形波导都满足这个条件但光子晶体波导或者强限制纳米线波导就要小心了后者的折射率对比度太高FD-BPM 会高估损耗。2.2 Crank-Nicolson 离散隐式格式为什么更稳接下来把上面的偏微分方程离散到网格上。横向 x 方向用均匀网格间距为 dx传播方向 z 用 dz 作步长。对 z 方向的偏导采用 Crank-Nicolson 格式即在 z 和 zdz 两个半程处取平均(ψ^(m1) - ψ^(m)) / dz (j/2k₀n₀) · [H_op ψ^(m1) H_op ψ^(m)] / 2其中 H_op 是横向算符包含二阶差分和折射率项。把含 ψ^(m1) 的项移到左边含 ψ^(m) 的移到右边就得到三对角线性方程组A · ψ^(m1) B · ψ^(m)A 和 B 都是三对角矩阵用 Thomas 算法追赶法可以高效求解每步计算量只有 O(Nx)Nx 是横向网格数。隐式格式的好处在于不管选多大的 dz数值上都是无条件稳定的——不过稳定不代表准确dz 还是得足够小才能分辨传播方向的相位变化。2.3 代码骨架网格初始化与算符构建打开FD_BPM_planar_waveguide.m观察它的初始化部分核心就是构建横向二阶差分算子和折射率分布数组。我一般会把它拆出来单独写成函数function [D2, n_profile] build_operator(x, lambda, n_clad, n_core, wg_width) % 构建横向二阶差分算子D2和折射率分布 % x: 横向坐标数组 % lambda: 工作波长 % n_clad, n_core: 包层和芯层折射率 % wg_width: 波导芯层宽度 dx x(2) - x(1); N length(x); % 二阶差分算子中心差分 diag_main -2 * ones(N, 1) / dx^2; diag_off ones(N-1, 1) / dx^2; D2 spdiags([diag_off diag_main diag_off], -1:1, N, N); % 波导折射率分布 n_profile n_clad * ones(size(x)); n_profile(abs(x) wg_width/2) n_core; end这段代码里的spdiags是 MATLAB 稀疏矩阵构建函数三对角矩阵在 N 达到几万时依然只需要 O(N) 的存储。折射率分布用逻辑索引直接赋值比 for 循环快一个数量级。实际使用中要注意dx的取值——我一般控制在 λ/(10·n_core) 以内太粗会引入数值色散太细则三对角求解时间线性增长不划算。3. 计算窗口与边界条件透明边界条件 TBC 的实现逻辑3.1 为什么不能用零边界波导模场是指数衰减的所以窗口边缘的场理论上应该是零。但实际数值模拟中如果强制设ψ 0在边界上衰减尾迹会在这个位置被反射回来形成驻波振荡在传播场分布图上表现为干涉条纹——这跟真实的波导辐射完全是两回事。解决思路有两个吸收边界加 PML和透明边界条件TBC。PML 效果好但要额外铺吸收层增加网格数平面波导这种一维横向问题用 TBC 就够了它不需要额外网格只需要在边界处估算波的出射方向然后反向补偿。3.2 TBC 的离散公式以右边界 x x_N 为例。假设边界处的场是出射平面波ψ(x) ≈ ψ(x_N) · exp(jk_x · (x - x_N))其中 k_x 是横向波数。利用边界内侧两点的场值来估算衰减率k_x (j/dx) · ln(ψ(x_N) / ψ(x_{N-1}))然后把这个 k_x 用在外推方程里ψ(x_{N1}) ψ(x_N) · exp(jk_x · dx)。实际操作中不需要真的弄一个虚拟点而是把 Crank-Nicolson 方程组里 A 矩阵的边界行改掉直接耦合边界处的传输条件。左边界对称处理。3.3 MATLAB 中的边界条件实现function [A, B] build_cn_matrices(D2, n_profile, k0, n0, dx, dz) % 构建 Crank-Nicolson 步进所需的稀疏矩阵 A 和 B % k0: 自由空间波数 2*pi/lambda % n0: 参考折射率 N length(n_profile); I speye(N); % 横向算符 H (1/2*k0*n0) * (D2 k0^2*(n_profile.^2 - n0^2)) H (1/(2*k0*n0)) * (D2 spdiags(k0^2*(n_profile.^2 - n0^2), 0, N, N)); A I - (1j*dz/2) * H; B I (1j*dz/2) * H; % ---- 右边界透明条件 ---- % 用端点两行修正 A 和 B使边界允许出射波透过 % 计算边界处横向波数近似值用当前步的场实现在步进循环内 % 这里预留结构A(1,:), B(1,:) 和 A(N,:), B(N,:) 需要动态更新 end注意H中的折射率项是逐点相乘对应的是横向非均匀波导结构。TBC 最大的坑在于 k_x 计算依赖当前步的场值所以 A 和 B 矩阵必须在每一步都重建边界行——完全不能在循环外面一次性算好。很多初学者把[A,B] build_cn_matrices(...)放到循环外结果边界反射一直消不掉就是这个原因。性能上不用担心MATLAB 稀疏矩阵重分配的耗时远小于三对角求解本身。3.4 参数选型的实际考量平面波导模拟中窗口宽度要取到模场直径的 35 倍以上。比如芯层宽度 4 μm、折射率差 0.01 的对称波导基模模场直径大约 10 μm窗口取 4060 μm横向网格点数 5121024 就够了。窗口太小即使有 TBC模场边界处衰减不充分引入的数值误差也会偏大。初学阶段可以用高斯光束做入射场观察传播过程有没有边界反射——如果能量曲线在某个位置出现回弹基本都是窗口宽度或 TBC 参数的问题。4.FD_BPM_planar_waveguide.m主流程拆解与参数初始化4.1 脚本主流程框图与文件结构整个脚本的核心流程可以归纳为四个阶段参数定义 → 网格与折射率构建 → 传播循环含 TBC 更新→ 可视化与分析。FD_BPM_planar_waveguide.m的主体就是围绕这四个阶段展开的只是把 TBC 跟主循环写在了一起。传播循环内部每步要做的事有三件根据当前场更新边界行 → 求解三对角方程 → 记录功率或场分布用于后处理。4.2 关键参数表与推荐取值平面波导 FD-BPM 模拟的精度九成取决于参数设置。把常见的可调参数整理一下建议拿到脚本先检查这些值再运行参数符号典型取值影响工作波长lambda1.55 μm决定 k₀ 和衍射尺度芯层折射率n_core3.45硅或 1.5SiO₂决定模式约束强度包层折射率n_clad与芯层差 0.0050.1差太大时 SVEA 失效芯层宽度wg_width28 μm直接决定模式数量横向网格间距dxλ/(10·n_core)太小影响模场分辨率传播步长dz0.11 μm影响相位精度传播总长度L数百至数千 μm观察模式稳定演化关于dz有个经验法则每步的相位变化控制在 π/20 以内也就是dz λ/(20·n_eff)n_eff 是模式有效折射率。对于 1.55 μm 波长、n_eff ≈ 1.5 的情况dz 取 0.050.1 μm 足够取大了可能出现非物理的振荡取小了纯粹浪费时间。如果只想看能量分布不关心相位细节dz 可以放宽到 0.5 μm 甚至更大。4.3 入射场构造高阶模激发与高斯近似FD-BPM 的初始场决定了最终激发哪些模式。最常见的是直接输入高斯光束% 高斯光束入射场构造 w0 5e-6; % 光束腰半径单位米 x linspace(-30e-6, 30e-6, 1024); % 横向窗口 % 高斯场分布注意归一化保持总功率为1 E_in exp(-(x/w0).^2); E_in E_in / sqrt(sum(abs(E_in).^2) * dx);这里把入射场归一化到总功率为 1方便后续分析传输损耗——功率衰减的 dB 值直接就是10*log10(P_z/P_0)。值得注意的是高斯场跟波导基模并不是一回事除非你特意把 w0 调成与基模模场宽度相当否则入射后必然激发部分高阶模和辐射模表现为模场横向展宽后趋于稳定。想看纯基模演化需要先用近似解析公式或迭代法算出基模分布再作为输入场。4.4 传播主循环稀疏矩阵求解与功率记录循环体内的核心操作是求解三对角方程。MATLAB 中直接对稀疏矩阵使用左除号即可% 主传播循环 n_steps round(L / dz); P_z zeros(n_steps, 1); psi E_in.; % 列向量 for m 1:n_steps % 更新透明边界条件需要当前场的边界值 [A, B] update_tbc_boundary(A, B, psi, k0, n0, dx); % Crank-Nicolson 步进求解 rhs B * psi_new; % 注意这一步需要上一时刻的 psi psi A \ rhs; % 记录功率 P_z(m) sum(abs(psi).^2) * dx; end注意A \ rhs是利用稀疏矩阵求解器解线性方程组MATLAB 会自动选择追赶法还是 LU 分解。对于典型的一维问题N 1024、n_steps 20000这段循环运行时间在几秒到十几秒之间。如果单位换算没问题这个循环的性能瓶颈在每一次都更新稀疏矩阵如果觉得慢可以固定 A 矩阵中与折射率无关的部分只更新边界行。5. 基模与高阶模的数值验证、辐射损耗分析与收敛性判断5.1 传播常数和模场形态的提取平直波导的 FD-BPM 跑到一定距离后场分布会趋于稳定——高次模和辐射模逐渐衰减殆尽剩下的是最低阶导模。这时候可以通过互相关操作提取传播常数% 提取有效折射率 psi_ref psi(:, end); % 选取末端场分布 n_eff -angle(psi(N/2, end) / psi(N/2, end-50)) / (k0 * 50 * dz); fprintf(提取的有效折射率: %.6f\n, n_eff);取中心点相邻两传播位置的相位差除以波数和传播距离就是有效折射率。这样做比直接看模式分布更定量也更适合对照解析解。对称三层平板波导的 n_eff 解析表达式是用特征方程解的MATLAB 的fzero函数很容易求两者对比误差一般能做到 1e-4 以内——如果误差到 1e-3 量级多半是横向窗口不够宽模场被边界压缩了。5.2 功率泄漏曲线和边界反射的判别模拟完成后画出P_z随传播距离的变化曲线dB 为单位理想情况下是一条平缓下降的直线斜率对应波导的泄漏损耗。如果曲线出现抖动优先怀疑两个原因一是 z 步长过大Crank-Nicolson 离散误差累积二是边界反射——TBC 更新频率不够或者是窗口宽度不够。判断方法是把窗口宽度加倍再看曲线如果抖动明显改善就是窗口问题反之则是 dz 问题。5.3 模式重叠积分分析入射场匹配度工程上更常用的是模式重叠积分overlap integral来分析入射场激发效率% 基模场分布 psi_fundamental可用解析近似或迭代得到 % 入射场 psi_in overlap abs(sum(conj(psi_fundamental) .* psi_in) * dx)^2 / ... (sum(abs(psi_fundamental).^2) * dx * sum(abs(psi_in).^2) * dx); fprintf(基模激发效率: %.2f%%\n, overlap * 100);这个值直接告诉你入射高斯场和波导基模的匹配程度。平面波导设计中耦合效率低于 90% 意味着需要调整光束腰半径或者加装模式转换结构。FD-BPM 的好处在于你可以直接在 z0 处扫描不同的 w0每个 w0 跑一次模拟看输出端的 overlap几步就能找到最优值——这个流程不需要改动脚本结构只需要把入射场构造部分参数化。6. 模式匹配法验证 FD-BPM 精度与收敛性FD-BPM 算得准不准不能只看模拟自洽必须锚定已知解。对称三层平板波导是最理想的对象因为它的 TE 模式特征方程可以精确求解MATLAB 里用fzero十几行代码就能搞定% TE0 模式特征方程求解 function n_eff solve_TE0(lambda, n_core, n_clad, wg_width) k0 2*pi/lambda; a wg_width / 2; u (n) sqrt(k0^2 * (n_core^2 - n.^2)); w (n) sqrt(k0^2 * (n.^2 - n_clad^2)); f (n) u(n)*tan(u(n)*a) - w(n); % 对称模特征方程 n_eff fzero(f, [n_clad, n_core]); end注意fzero需要一个区间作为初始猜测这里[n_clad, n_core]就是物理上允许的有效折射率范围。拿这个值对照 FD-BPM 提取的结果可以画一条误差随 dz 变化的曲线验证收敛阶数——理论上 Crank-Nicolson 格式在 z 方向是二阶精度误差应该随 dz 平方递减。如果实测误差下降速度明显慢于二阶说明空间分辨率不够要先加密 dx 再测。这个方法同样适用于多模波导——每种模式用不同初始场激发然后对比各自的传播常数。做完这一轮验证脚本的可靠性就有了锚点之后再拿它去设计弯曲波导、锥形过渡或者 Y 分支结论才站得住。FD-BPM 的优势在于它把三维问题降成了二维扫描模式匹配法给的又是精确基准两者配合是平面波导设计前期最实用的一对组合拳。本文还有配套的精品资源点击获取 基于MATLAB的FD-BPM平面波导模拟：从理论到代码实现
简介面向光子学研究者的波导仿真MATLAB代码采用有限差分光束传播法FD-BPM计算光场在波导中的传播适合波导、光栅及微结构光纤等结构的数值模拟与优化设计。压缩包包含1个m格式脚本整体大小约1KB为极简的平面波导FD-BPM实现便于快速部署与修改调试。脚本完整涵盖初始化、网格划分、光场离散、传播迭代、边界条件处理及结果分析等流程并支持通过调整波导几何尺寸、光源波长与初始场分布来考察模式分布、功率传输和损耗特性。已有347人学习下载对于从事光子器件课题研究或课程设计的学生与工程师这份轻量代码可直接运行参考省去从零编写算法的过程也能加深对BPM数值方法的理解。得益于FD-BPM的数值特性脚本在常见PC上即可高效运行适合教学演示与前期方案验证。1. 平面波导模拟为何绕不开 FD-BPM 这一步做光子集成的人大概都有过这种经历结构画好了、参数算好了一进工艺线出来的器件响应和设计值差出一大截。问题常常不出在工艺而出在模拟阶段——很多人直接拿 FDTD 当万能工具却忽略了在平面波导这种弱约束、长传播场景下全波方法计算量大得离谱几毫米的传播长度足够让普通工作站跑上一整天。FD-BPM 的价值正在于此它把麦克斯韦方程降维成单向传播的抛物型方程用有限差分在横截面离散在传播方向上步进求解几秒钟就能看到整个光场演化过程。MATLAB 里用 FD-BPM 研究波导不需要昂贵的商用光学软件也不需要集群资源一台笔记本就能完成从模场分析到传输损耗估算的全部工作。这篇要讲的FD_BPM_planar_waveguide.m正是典型实现适合刚接触光子学仿真的人快速跑通也适合老手用来做结构初筛和参数扫描。2. 从亥姆霍兹方程到 Crank-Nicolson 格式FD-BPM 的理论地基2.1 缓变包络近似与单向传播假设标量亥姆霍兹方程是光场传播的起点∇²E k₀²n²(x, z)E 0直接数值求解这个椭圆型方程需要同时处理前向和后向传播的场边界条件复杂且计算量极大。BPM 的核心做法是假设光场可以写成快变相位和缓变包络的乘积E(x, z) ψ(x, z) · exp(-jk₀n₀z)其中 n₀ 是参考折射率ψ 是缓变包络。把这个代入亥姆霍兹方程并忽略 ψ 对 z 的二阶导数项这就是缓变包络近似SVEA得到∂ψ/∂z (j/2k₀n₀) · [∂²ψ/∂x² k₀²(n²(x, z) - n₀²)ψ]这个方程是抛物型方程只允许前向传播所以可以用步进法从 z0 逐层推进到 zL。这个近似的物理意义就是要求波导结构沿着传播方向变化足够缓慢折射率对比度也不宜太大——大部分平面波导和脊形波导都满足这个条件但光子晶体波导或者强限制纳米线波导就要小心了后者的折射率对比度太高FD-BPM 会高估损耗。2.2 Crank-Nicolson 离散隐式格式为什么更稳接下来把上面的偏微分方程离散到网格上。横向 x 方向用均匀网格间距为 dx传播方向 z 用 dz 作步长。对 z 方向的偏导采用 Crank-Nicolson 格式即在 z 和 zdz 两个半程处取平均(ψ^(m1) - ψ^(m)) / dz (j/2k₀n₀) · [H_op ψ^(m1) H_op ψ^(m)] / 2其中 H_op 是横向算符包含二阶差分和折射率项。把含 ψ^(m1) 的项移到左边含 ψ^(m) 的移到右边就得到三对角线性方程组A · ψ^(m1) B · ψ^(m)A 和 B 都是三对角矩阵用 Thomas 算法追赶法可以高效求解每步计算量只有 O(Nx)Nx 是横向网格数。隐式格式的好处在于不管选多大的 dz数值上都是无条件稳定的——不过稳定不代表准确dz 还是得足够小才能分辨传播方向的相位变化。2.3 代码骨架网格初始化与算符构建打开FD_BPM_planar_waveguide.m观察它的初始化部分核心就是构建横向二阶差分算子和折射率分布数组。我一般会把它拆出来单独写成函数function [D2, n_profile] build_operator(x, lambda, n_clad, n_core, wg_width) % 构建横向二阶差分算子D2和折射率分布 % x: 横向坐标数组 % lambda: 工作波长 % n_clad, n_core: 包层和芯层折射率 % wg_width: 波导芯层宽度 dx x(2) - x(1); N length(x); % 二阶差分算子中心差分 diag_main -2 * ones(N, 1) / dx^2; diag_off ones(N-1, 1) / dx^2; D2 spdiags([diag_off diag_main diag_off], -1:1, N, N); % 波导折射率分布 n_profile n_clad * ones(size(x)); n_profile(abs(x) wg_width/2) n_core; end这段代码里的spdiags是 MATLAB 稀疏矩阵构建函数三对角矩阵在 N 达到几万时依然只需要 O(N) 的存储。折射率分布用逻辑索引直接赋值比 for 循环快一个数量级。实际使用中要注意dx的取值——我一般控制在 λ/(10·n_core) 以内太粗会引入数值色散太细则三对角求解时间线性增长不划算。3. 计算窗口与边界条件透明边界条件 TBC 的实现逻辑3.1 为什么不能用零边界波导模场是指数衰减的所以窗口边缘的场理论上应该是零。但实际数值模拟中如果强制设ψ 0在边界上衰减尾迹会在这个位置被反射回来形成驻波振荡在传播场分布图上表现为干涉条纹——这跟真实的波导辐射完全是两回事。解决思路有两个吸收边界加 PML和透明边界条件TBC。PML 效果好但要额外铺吸收层增加网格数平面波导这种一维横向问题用 TBC 就够了它不需要额外网格只需要在边界处估算波的出射方向然后反向补偿。3.2 TBC 的离散公式以右边界 x x_N 为例。假设边界处的场是出射平面波ψ(x) ≈ ψ(x_N) · exp(jk_x · (x - x_N))其中 k_x 是横向波数。利用边界内侧两点的场值来估算衰减率k_x (j/dx) · ln(ψ(x_N) / ψ(x_{N-1}))然后把这个 k_x 用在外推方程里ψ(x_{N1}) ψ(x_N) · exp(jk_x · dx)。实际操作中不需要真的弄一个虚拟点而是把 Crank-Nicolson 方程组里 A 矩阵的边界行改掉直接耦合边界处的传输条件。左边界对称处理。3.3 MATLAB 中的边界条件实现function [A, B] build_cn_matrices(D2, n_profile, k0, n0, dx, dz) % 构建 Crank-Nicolson 步进所需的稀疏矩阵 A 和 B % k0: 自由空间波数 2*pi/lambda % n0: 参考折射率 N length(n_profile); I speye(N); % 横向算符 H (1/2*k0*n0) * (D2 k0^2*(n_profile.^2 - n0^2)) H (1/(2*k0*n0)) * (D2 spdiags(k0^2*(n_profile.^2 - n0^2), 0, N, N)); A I - (1j*dz/2) * H; B I (1j*dz/2) * H; % ---- 右边界透明条件 ---- % 用端点两行修正 A 和 B使边界允许出射波透过 % 计算边界处横向波数近似值用当前步的场实现在步进循环内 % 这里预留结构A(1,:), B(1,:) 和 A(N,:), B(N,:) 需要动态更新 end注意H中的折射率项是逐点相乘对应的是横向非均匀波导结构。TBC 最大的坑在于 k_x 计算依赖当前步的场值所以 A 和 B 矩阵必须在每一步都重建边界行——完全不能在循环外面一次性算好。很多初学者把[A,B] build_cn_matrices(...)放到循环外结果边界反射一直消不掉就是这个原因。性能上不用担心MATLAB 稀疏矩阵重分配的耗时远小于三对角求解本身。3.4 参数选型的实际考量平面波导模拟中窗口宽度要取到模场直径的 35 倍以上。比如芯层宽度 4 μm、折射率差 0.01 的对称波导基模模场直径大约 10 μm窗口取 4060 μm横向网格点数 5121024 就够了。窗口太小即使有 TBC模场边界处衰减不充分引入的数值误差也会偏大。初学阶段可以用高斯光束做入射场观察传播过程有没有边界反射——如果能量曲线在某个位置出现回弹基本都是窗口宽度或 TBC 参数的问题。4.FD_BPM_planar_waveguide.m主流程拆解与参数初始化4.1 脚本主流程框图与文件结构整个脚本的核心流程可以归纳为四个阶段参数定义 → 网格与折射率构建 → 传播循环含 TBC 更新→ 可视化与分析。FD_BPM_planar_waveguide.m的主体就是围绕这四个阶段展开的只是把 TBC 跟主循环写在了一起。传播循环内部每步要做的事有三件根据当前场更新边界行 → 求解三对角方程 → 记录功率或场分布用于后处理。4.2 关键参数表与推荐取值平面波导 FD-BPM 模拟的精度九成取决于参数设置。把常见的可调参数整理一下建议拿到脚本先检查这些值再运行参数符号典型取值影响工作波长lambda1.55 μm决定 k₀ 和衍射尺度芯层折射率n_core3.45硅或 1.5SiO₂决定模式约束强度包层折射率n_clad与芯层差 0.0050.1差太大时 SVEA 失效芯层宽度wg_width28 μm直接决定模式数量横向网格间距dxλ/(10·n_core)太小影响模场分辨率传播步长dz0.11 μm影响相位精度传播总长度L数百至数千 μm观察模式稳定演化关于dz有个经验法则每步的相位变化控制在 π/20 以内也就是dz λ/(20·n_eff)n_eff 是模式有效折射率。对于 1.55 μm 波长、n_eff ≈ 1.5 的情况dz 取 0.050.1 μm 足够取大了可能出现非物理的振荡取小了纯粹浪费时间。如果只想看能量分布不关心相位细节dz 可以放宽到 0.5 μm 甚至更大。4.3 入射场构造高阶模激发与高斯近似FD-BPM 的初始场决定了最终激发哪些模式。最常见的是直接输入高斯光束% 高斯光束入射场构造 w0 5e-6; % 光束腰半径单位米 x linspace(-30e-6, 30e-6, 1024); % 横向窗口 % 高斯场分布注意归一化保持总功率为1 E_in exp(-(x/w0).^2); E_in E_in / sqrt(sum(abs(E_in).^2) * dx);这里把入射场归一化到总功率为 1方便后续分析传输损耗——功率衰减的 dB 值直接就是10*log10(P_z/P_0)。值得注意的是高斯场跟波导基模并不是一回事除非你特意把 w0 调成与基模模场宽度相当否则入射后必然激发部分高阶模和辐射模表现为模场横向展宽后趋于稳定。想看纯基模演化需要先用近似解析公式或迭代法算出基模分布再作为输入场。4.4 传播主循环稀疏矩阵求解与功率记录循环体内的核心操作是求解三对角方程。MATLAB 中直接对稀疏矩阵使用左除号即可% 主传播循环 n_steps round(L / dz); P_z zeros(n_steps, 1); psi E_in.; % 列向量 for m 1:n_steps % 更新透明边界条件需要当前场的边界值 [A, B] update_tbc_boundary(A, B, psi, k0, n0, dx); % Crank-Nicolson 步进求解 rhs B * psi_new; % 注意这一步需要上一时刻的 psi psi A \ rhs; % 记录功率 P_z(m) sum(abs(psi).^2) * dx; end注意A \ rhs是利用稀疏矩阵求解器解线性方程组MATLAB 会自动选择追赶法还是 LU 分解。对于典型的一维问题N 1024、n_steps 20000这段循环运行时间在几秒到十几秒之间。如果单位换算没问题这个循环的性能瓶颈在每一次都更新稀疏矩阵如果觉得慢可以固定 A 矩阵中与折射率无关的部分只更新边界行。5. 基模与高阶模的数值验证、辐射损耗分析与收敛性判断5.1 传播常数和模场形态的提取平直波导的 FD-BPM 跑到一定距离后场分布会趋于稳定——高次模和辐射模逐渐衰减殆尽剩下的是最低阶导模。这时候可以通过互相关操作提取传播常数% 提取有效折射率 psi_ref psi(:, end); % 选取末端场分布 n_eff -angle(psi(N/2, end) / psi(N/2, end-50)) / (k0 * 50 * dz); fprintf(提取的有效折射率: %.6f\n, n_eff);取中心点相邻两传播位置的相位差除以波数和传播距离就是有效折射率。这样做比直接看模式分布更定量也更适合对照解析解。对称三层平板波导的 n_eff 解析表达式是用特征方程解的MATLAB 的fzero函数很容易求两者对比误差一般能做到 1e-4 以内——如果误差到 1e-3 量级多半是横向窗口不够宽模场被边界压缩了。5.2 功率泄漏曲线和边界反射的判别模拟完成后画出P_z随传播距离的变化曲线dB 为单位理想情况下是一条平缓下降的直线斜率对应波导的泄漏损耗。如果曲线出现抖动优先怀疑两个原因一是 z 步长过大Crank-Nicolson 离散误差累积二是边界反射——TBC 更新频率不够或者是窗口宽度不够。判断方法是把窗口宽度加倍再看曲线如果抖动明显改善就是窗口问题反之则是 dz 问题。5.3 模式重叠积分分析入射场匹配度工程上更常用的是模式重叠积分overlap integral来分析入射场激发效率% 基模场分布 psi_fundamental可用解析近似或迭代得到 % 入射场 psi_in overlap abs(sum(conj(psi_fundamental) .* psi_in) * dx)^2 / ... (sum(abs(psi_fundamental).^2) * dx * sum(abs(psi_in).^2) * dx); fprintf(基模激发效率: %.2f%%\n, overlap * 100);这个值直接告诉你入射高斯场和波导基模的匹配程度。平面波导设计中耦合效率低于 90% 意味着需要调整光束腰半径或者加装模式转换结构。FD-BPM 的好处在于你可以直接在 z0 处扫描不同的 w0每个 w0 跑一次模拟看输出端的 overlap几步就能找到最优值——这个流程不需要改动脚本结构只需要把入射场构造部分参数化。6. 模式匹配法验证 FD-BPM 精度与收敛性FD-BPM 算得准不准不能只看模拟自洽必须锚定已知解。对称三层平板波导是最理想的对象因为它的 TE 模式特征方程可以精确求解MATLAB 里用fzero十几行代码就能搞定% TE0 模式特征方程求解 function n_eff solve_TE0(lambda, n_core, n_clad, wg_width) k0 2*pi/lambda; a wg_width / 2; u (n) sqrt(k0^2 * (n_core^2 - n.^2)); w (n) sqrt(k0^2 * (n.^2 - n_clad^2)); f (n) u(n)*tan(u(n)*a) - w(n); % 对称模特征方程 n_eff fzero(f, [n_clad, n_core]); end注意fzero需要一个区间作为初始猜测这里[n_clad, n_core]就是物理上允许的有效折射率范围。拿这个值对照 FD-BPM 提取的结果可以画一条误差随 dz 变化的曲线验证收敛阶数——理论上 Crank-Nicolson 格式在 z 方向是二阶精度误差应该随 dz 平方递减。如果实测误差下降速度明显慢于二阶说明空间分辨率不够要先加密 dx 再测。这个方法同样适用于多模波导——每种模式用不同初始场激发然后对比各自的传播常数。做完这一轮验证脚本的可靠性就有了锚点之后再拿它去设计弯曲波导、锥形过渡或者 Y 分支结论才站得住。FD-BPM 的优势在于它把三维问题降成了二维扫描模式匹配法给的又是精确基准两者配合是平面波导设计前期最实用的一对组合拳。本文还有配套的精品资源点击获取