简介面向光子学与光纤通信研究者的有限差分光束传播法波导计算资源在MATLAB中模拟平面波导光场传播解决复杂结构光传输特性的快速仿真与参数分析问题。压缩包内共1个M脚本文件体量仅1KB小巧精简便于直接查看和修改代码逻辑。已有347人学习适合正在学习导波光学、需要参考数值实现思路的本科生或工程师。脚本完整涵盖了波导几何初始化、网格划分、光场离散、传播迭代及边界条件处理等核心环节通过修改波长、波导宽度、折射率分布等参数可直观观察到模式分布、功率传输与损耗特性的变化帮助读者将算法从公式落地为可运行代码并作为进一步开发三维或复杂结构仿真工具的基础模板。1. FD BPM 在波导计算里到底算什么FD BPM 这个名字直觉上会让人以为是“用有限差分计算波导模式”。实际上 FD-BPM 算的不是模式本征值而是一束光在波导里沿纵向传播的场演化过程。对于只有直波导、截面不变的结构用 FDE有限差分本征模更直接但涉及弯曲波导、Y 分支、锥形过渡、模式失配和损耗提取时FDE 没有现成的解析解这时就轮到 FD-BPM 上场。用 MATLAB 实现一个可用的二维标量 FD-BPM核心代码能压在 100 行内不需要任何工具箱。真正决定结果质量的不是代码而是三个选择离散格式、步长和边界条件。下面按一套能直接落地的顺序讲所有长度单位统一用微米。2. FD-BPM 的波动方程离散从旁轴近似到 Crank-Nicolson 格式2.1 缓变包络近似SVEA成立的条件标量亥姆霍兹方程在二维波导截面 x-z 平面内写作[ \frac{\partial^2 E}{\partial x^2} \frac{\partial^2 E}{\partial z^2} k_0^2 n^2(x,z) E 0 ]FD-BPM 的第一步是假设场可以拆成快变载波和慢变包络[ E(x,z) \psi(x,z) e^{-i k_0 n_r z} ]这里 (n_r) 是参考折射率一般取包层或衬底折射率。代入亥姆霍兹方程后忽略 (\partial^2 \psi / \partial z^2) 这一项得到傍轴波动方程[ 2 i k_0 n_r \frac{\partial \psi}{\partial z} \frac{\partial^2 \psi}{\partial x^2} k_0^2 \left(n^2(x,z) - n_r^2\right)\psi ]这个近似就是缓变包络近似也叫旁轴近似。它是否成立要看两个条件。第一折射率对比度不能太大标量 FD-BPM 对 (\Delta n / n 0.05) 的弱导波导比较可靠硅波导那种 (\Delta n \approx 1.0) 的大对比结构标量近似会出现明显误差需要矢量修正。第二后向反射必须可以忽略纵向出现强突变界面时BPM 会把反射能量错误地处理成前向传播这时要么改用 FDTD要么在界面处单独做反射修正。实际工程中判断 SVEA 是否够用的另一个经验指标是弯曲半径。当波导弯曲半径小于几十倍芯宽时波前不再垂直于传播轴旁轴假设逐渐失真。后面第 5 章会看到弯曲半径对 BPM 结果的影响可以通过等效折射率扰动来评估但前提仍然是 (R) 足够大。2.2 用 Crank-Nicolson 格式推进 z 方向把上面的傍轴方程写成算子形式[ \frac{\partial \psi}{\partial z} -i \left( D V \right) \psi ]其中 (D) 是二阶空间导数相关的扩散算子(V k_0 (n^2 - n_r^2) / (2 n_r)) 是折射率扰动算子。对 z 方向做 Crank-Nicolson 离散后每一步需要解一个复三对角线性方程组[ \left( I \frac{i \Delta z}{2}(DV) \right) \psi_{m1} \left( I - \frac{i \Delta z}{2}(DV) \right) \psi_m ]用中心差分把 x 方向的二阶导数离散成标准三对角结构后左端矩阵的主对角、下对角、上对角都可以显式写出右端则是矩阵乘向量。这就是 FD-BPM 在 MATLAB 里能高效运行的原因每一步求解一个三对角系统复杂度是 O(N)而不是 O(N²) 或更高。Crank-Nicolson 格式无条件稳定这让很多初学者误以为 (\Delta z) 可以随便取。无条件稳定只代表误差不会指数发散不代表数值色散和相位误差可以忽略。(\Delta z) 过大时每步推进的相位积累会偏离真实传播常数最终提取出的 neff 会偏大或偏小。这一点在后面第 4 章的参数规则里会专门展开。2.3 二维标量还是三维矢量按问题规模定标题这个活最常见的第一版实现是二维标量 BPM因为它能直接回答“场怎么传播”“模式相位是多少”“弯曲损耗大致多大”这类问题。二维 BPM 的网格是一个二维矩阵内存占用小调试方便适合把物理机制先摸清楚。三维条形波导如果要严格仿真有两种路线。一种是把三维场拆成两个一维问题交替推进也就是 ADI-BPM每一轮分别对 x 方向和 y 方向做隐式推进另一种是直接建立三维体网格做全矢量 BPM。后者的内存增长非常快例如 256×256×800 的复数网格用 MATLAB 默认的 complex double 存储大约需要 1 GB 内存还没算折射率数组和临时变量。所以三维场景下网格切片、数据类型转换和传播距离都要重新规划不能直接把二维代码搬过去。离散格式稳定性每步复杂度适合场景显式FDTD式推进受 CFL 条件约束O(N)时间域反射和复杂结构Crank-Nicolson BPM无条件稳定O(N) 三对角波导纵向前向传播ADI-BPM3D无条件稳定两次三对角 O(N)3D 条形波导纵向仿真表格里 ADI-BPM 的实现比二维要麻烦不少但仍是三维标量仿真里较经济的方案。做项目前先评估折射率差和内存再决定是从二维起步还是直接上三维。3. MATLAB 实现 FD-BPM 的最小可运行程序3.1 参数定义与长度单位选择所有物理长度统一用微米这一步看似简单但能省掉大量调试时间。把波长设为 1.55 μm 时(k_0 2\pi/1.55 \approx 4.05\ \mu m^{-1})配合 dx0.02 μm中间计算量在 (10^{-3}) 到 (10^2) 之间浮点误差可控。如果全部用国际单位k0 接近 (4 \times 10^6)矩阵元素数量级跨度大打印和调参都不方便。这里以一个弱导直波导为例参数表如下。参数符号本文取值说明工作波长lambda1.55 μm通信 C 波段参考折射率nr1.444取包层折射率芯层折射率差dn0.003弱导符合标量近似芯区宽度core_w4 μm矩形波导芯宽横向网格间距dx0.02 μm约为波长的 1/50纵向推进步长dz0.5 μm约 1/3 波长吸收层宽度w_abs5 μm左右各一段吸收层最大虚部alpha_max0.05 μm^-1二次型渐变3.2 主脚本折射率剖面、高斯初始场和吸收层下面这段代码是完整的可运行主脚本函数定义放在后面单独列出。代码里用复数折射率虚部模拟吸收层这是最常见的做法不需要额外写 PML 的坐标拉伸。% fd_bpm_waveguide.m % 二维标量 FD-BPM 直波导传播示例 % 运行环境MATLAB R2023b较早版本也兼容 clear; close all; clc; % 物理参数长度单位 um lambda 1.55; k0 2*pi/lambda; nr 1.444; dn 0.003; core_w 4; % 横向网格 dx 0.02; Lx 32; x (-Lx/2:dx:Lx/2).; nx numel(x); % 折射率分布矩形芯 n nr * ones(nx, 1); n(abs(x) core_w/2) nr dn; % 吸收层在左右两端加二次型复数折射率虚部 w_abs 5; alpha_max 0.05; idx_left x -Lx/2 x -Lx/2 w_abs; t_left (x(idx_left) - (-Lx/2)) / w_abs; n(idx_left) n(idx_left) - 1i * alpha_max * t_left.^2; idx_right x Lx/2 x Lx/2 - w_abs; t_right (x(idx_right) - (Lx/2 - w_abs)) / w_abs; n(idx_right) n(idx_right) - 1i * alpha_max * t_right.^2; % 高斯初始场束腰 w0 w0 3; phi exp(-(x/w0).^2); phi phi / sqrt(sum(abs(phi).^2) * dx); % 传播参数 Lz 400; dz 0.5; nz round(Lz/dz); z_list (0:nz) * dz; % 输出与监视变量 power zeros(nz1, 1); power(1) sum(abs(phi).^2) * dx; [~, imon] min(abs(x)); % x0 处监视相位 phase_mon zeros(nz1, 1); phase_mon(1) angle(phi(imon)); % 传播主循环 for step 1:nz phi bpm_2d_step(phi, dx, dz, k0, n, nr); power(step1) sum(abs(phi).^2) * dx; phase_mon(step1) angle(phi(imon)); end % 提取有效折射率跳过初始瞬态 phase_unwrap unwrap(phase_mon); idx_fit 500:nz1; pfit polyfit(z_list(idx_fit), phase_unwrap(idx_fit), 1); neff_fitted nr - pfit(1) / k0; fprintf(neff ~ %.6f\n, neff_fitted); % 输出强度演化 figure; imagesc(z_list, x, abs(phi_stack).^2); xlabel(z (um)); ylabel(x (um)); axis xy; colorbar; title(Intensity evolution); exportgraphics(gcf, fd_bpm_intensity.eps, ContentType,vector);代码里有几个容易忽略的点。折射率虚部是二次型渐变而不是常数这样吸收区与物理区交界处的介电常数突变会小一些反射也小。高斯初始场在芯区束腰为 3 μm直波导基模模场半径未必正好等于 3 μm所以这个初始场会同时激发出基模和少量高阶模。如果只关心直波导的 neff传播距离要足够长让高阶模在吸收边界作用下衰减掉或者干脆用模式源入射。neff 提取用到了 unwrap 和 polyfit。相位在传播中不断累积直接取 angle 会在 -π 到 π 之间跳变必须先用 unwrap 展开。拟合斜率 pfit(1) 是慢变包络的相位变化率总场相位变化率是 (k_0 n_r pfit(1))所以 neff nr pfit(1)/k0。由于这里的相位在传播中是递减的pfit(1) 是负值实际公式写作 neff nr - pfit(1)/k0才能得到大于 nr 的结果。3.3 传播步函数与三对角求解器主循环里的bpm_2d_step负责单步推进tridiagonal_solve是标准的 Thomas 算法。把求解器单独拆出来方便后续换成稀疏矩阵求解或 GPU 版本。function phi_new bpm_2d_step(phi, dx, dz, k0, n, nr) % 二维标量 FD-BPM 单步推进 % phi 为列向量n 为折射率列向量 phi phi(:); N numel(phi); alpha dz / (4 * k0 * nr * dx^2); beta dz * k0 * (n.^2 - nr^2) / (4 * nr); % 左端矩阵 A I - i*alpha*T - i*beta ld -1i * alpha * ones(N-1, 1); md 1 2i * alpha - 1i * beta; ud -1i * alpha * ones(N-1, 1); % 右端向量 B*phiB I i*alpha*T i*beta rhs (1 - 2i * alpha 1i * beta) .* phi; rhs(1:N-1) rhs(1:N-1) 1i * alpha .* phi(2:N); rhs(2:N) rhs(2:N) 1i * alpha .* phi(1:N-1); phi_new tridiagonal_solve(ld, md, ud, rhs); end function x tridiagonal_solve(ld, md, ud, rhs) % Thomas 算法求解复三对角方程组 N length(md); cp zeros(N-1, 1); dp zeros(N, 1); cp(1) ud(1) / md(1); dp(1) rhs(1) / md(1); for i 2:N-1 denom md(i) - ld(i-1) * cp(i-1); cp(i) ud(i) / denom; dp(i) (rhs(i) - ld(i-1) * dp(i-1)) / denom; end x zeros(N, 1); x(N) (rhs(N) - ld(N-1) * dp(N-1)) / (md(N) - ld(N-1) * cp(N-1)); for i N-1:-1:1 x(i) dp(i) - cp(i) * x(i1); end end传播函数里 beta 是向量因为它随折射率逐点变化。注意1i * alpha和1i * beta前面符号相反这是由 Crank-Nicolson 左端和右端展开决定的。写代码时如果 neff 结果异常优先检查这一组符号。Thomas 算法相比 MATLAB 内置的A\b优势在于避免了每次步进时创建稀疏矩阵和做符号解析。对于 nx 在 1000 点以上的网格手写 Thomas 的循环在 R2023b 下仍可能比稀疏矩阵直接求逆快 20% 到 50%而且内存占用更稳定。想进一步加速可以把for循环改成mex或 MATLAB Coder 生成的 MEX 文件物理不变。3.4 功率监测和可视化怎么布点功率监测是判断程序是否稳定的第一道防线。在吸收层存在时总功率会有少量下降这是正常现象如果功率异常增长说明矩阵符号错误或步长过大。实际工程中我会在 z 方向每 20 μm 记录一个截面而不是只记录最终场。这样既能看拍频也能定位在哪一段出现了异常反射。上面主脚本里用exportgraphics导出 eps 矢量图。这个函数在 R2023b 之后的版本比较稳定老版本环境可以用print -depsc2代替。注意exportgraphics导出彩色图时默认走 RGB 空间如果投稿需要灰度图记得先colormap gray再导出。4. FD-BPM 参数怎么调步长、边界窗口与折射率剖面4.1 dx 与 dz 的经验规则和收敛验证横向步长 dx 的经验值是 (\lambda / (20 \sim 40, n_{core}))。对于 1.55 μm、n≈1.444 的弱导波导dx 在 0.03 到 0.05 μm 之间即可。dx 取得过小比如 0.005 μm矩阵规模翻好几倍但对 neff 精度提升有限因为限制精度的主要因素已经变成旁轴近似本身。dz 的经验范围是 0.2 到 2 μm弱导场景下不建议超过 1 μm。网格参数经验范围对结果的影响dx(\lambda/(20 \sim 40,n_{core}))横向模场分辨率和数值色散dz(\lambda/2 \sim 2\lambda)相位累积准确度Lx模场宽度的 8 到 10 倍边界反射强度吸收层宽度2 到 6 μm反射抑制效果收敛验证的通用做法是跑三组基准网格、dx 减半、dz 减半。观察 neff 和输出模场。如果 neff 变化小于 (10^{-4})网格基本收敛。注意 dz 减半时传播步数翻倍总传播距离不变所以相位比较要在相同 z 位置取值最稳妥的是每个循环都存 z_list不要只存 step 编号。4.2 吸收边界与计算窗口从衰减窗到 PML 実部衰减窗的虚部折射率不是越大越好。alpha_max 取 0.05 μm^-1 时吸收区对入射场有适度衰减但如果取值到 0.5吸收区起始处会出现新的折射率跳变反而反射增强。常见做法是二次型渐变配合 4 到 6 μm 的吸收宽度反射可以压到 -60 dB 以下。如果还需要更低的反射就把吸收区宽度加大到 8 μm而不是一味提高 alpha_max。PML 可以看成衰减窗的升级版区别在于 PML 通过坐标拉伸让电磁波在介质中指数衰减不依赖介电常数虚部。实现上要把方程里的 (\partial/\partial x) 替换成 (\partial/\partial x \cdot (1/(1i\sigma(x)/\omega)))在 MATLAB 里会改变三对角矩阵的系数但仍然保持三对角结构。对大多数弱导波导仿真二次型衰减窗已经够用PML 更适合强反射结构或需要极高精度损耗计算时使用。4.3 折射率剖面的锯齿与平滑处理矩形芯在离散网格上会产生阶梯边界模拟结果里会出现数值散射和模式畸变。网格越粗越明显。工程上常见做法是给折射率剖面加一个边缘过渡过渡宽度取 0.1 到 0.3 μm。transition 0.15; d_core abs(x) - core_w/2; smooth_profile (1 erf(-d_core / transition)) / 2; n_smooth nr dn * smooth_profile;erf 平滑函数的好处是导数连续没有跳变。加了平滑后neff 会比直角芯的矩形波导略低这是正常现象。需要和实验对比时芯区实际折射率分布往往也不是完美矩形这个平滑参数可以当成工艺拟合量来调。4.4 程序没问题的基准测试与常见症状对照表写好的 BPM 程序先不要直接跑波导先跑均匀介质中的高斯光束传播。自由空间高斯光束的束宽解析解为[ w(z) w_0 \sqrt{1 \left(\frac{z}{k_0 w_0^2}\right)^2} ]用 BPM 算出每个截面的二阶矩宽度与上面的解析解对比偏差应该小于 1%。如果这个基准测试过不了后面的波导结果都不可信。症状可能原因检查方法功率在传播中增大Crank-Nicolson 矩阵符号错检查 alpha、beta 前的 i 符号输出场在左右边界出现亮线吸收层不足或 alpha_max 过大加宽吸收层降低 alpha_maxneff 始终偏大dz 过大导致相位累积误差减小 dz 到 0.2 μm 附近场出现高频抖动dx 太大或初始场有高频分量减小 dx检查初始场是否光滑输出模场有明显拍频初始场激发出多个模式增加传播距离或用模式源这些症状里拍频最容易被忽略。高斯束腰 3 μm 入射一个 4 μm 芯宽的弱导波导基模和高阶模的重叠积分都不为零模场会在传播方向上周期起伏周期大约等于两个模式的拍频波长。要避免拍频影响 neff 提取可以把初始场束腰调到与基模模场半径接近通常在芯宽的 1.1 到 1.2 倍之间。5. 进阶用 FD-BPM 算弯曲波导损耗与基模筛选5.1 弯曲波导的等效折射率扰动法弯曲波导的严格仿真要建立曲线坐标系但对弯曲半径远大于芯宽的常见场景可以直接用等效折射率扰动。原理是弯曲后波导外侧光程长、内侧光程短等效直波导的折射率在横截面上变成[ n_{eq}(x) n(x) \left(1 \frac{x}{R}\right) ]其中 x 的方向是弯曲平面内远离曲率中心的方向R 是弯曲半径。在 MATLAB 里这段扰动可以在每个需要模拟弯曲的 z 区间临时加到折射率数组上。R_bend 500; % 弯曲半径微米 n_bend n .* (1 x / R_bend); % 将这个 n_bend 传入 bpm_2d_step 替代原来的 n弯曲半径 R 越小等效折射率梯度越大模场会向外侧偏移并产生辐射损耗。扫描 R 从 2000 μm 逐渐减小到 200 μm观察监测点功率剩余量就能得到弯曲损耗曲线。注意这个近似在 R 小于几十微米时明显失真因为旁轴近似本身和曲线坐标的偏差已经不可忽略。5.2 基模筛选用 padded mode 代替特征求解如果不想额外写本征模求解器有一种很实用的替代方案用一个宽初始场或者随机场入射让它在几百微米长的 BPM 窗口里传播。基模损耗最低高阶模损耗高再加上吸收层传播足够距离后剩下来的场基本就是基模分布。这就是常见做法里的 padded mode 方法。实际使用时把传播距离设到 1000 μm 以上每隔 50 μm 记录一次功率。当功率衰减曲线进入单指数下降阶段后取那段距离的末态场做归一化再作为后续更精确仿真的入射源。需要注意这种方法只适合模式损耗有明显差异的波导如果两个模式损耗接近分离效果就会变差。对 1.55 μm 弱导直波导的最终推荐参数是dx0.02 μmdz0.5 μmLx32 μm吸收层 5 μmalpha_max0.05 μm^-1。先用均匀介质基准测试确认程序正确再跑 200 μm 检查功率曲线最后正式跑 400 μm 提取 neff 和模场。这样能避开绝大多数参数坑拿到稳定可复现的结果。本文还有配套的精品资源点击获取