简介：面向光学仿真学习者与二维光子晶体研究者的COMSOL-MATLAB联合计算示例，围绕周期性结构的光子带隙求解展开，适合已具备一定电磁场与脚本基础、希望把建模流程自动化的读者。资源包共4个文件，含MATLAB脚本、FIG图形、PNG结果图与Markdown说明，约144KB，分别对应仿真控制、结果可视化、效果预览和项目说明，便于先看说明再复现带隙计算与能带图。已有1066人学习，说明该方向具有一定关注度。读者可从中了解如何用脚本设置晶格与频率范围、调用COMSOL求解并读取场分布或S参数，再据此识别带隙边界、比较结果并优化设计，获得一套可拆解、可复用的计算思路。

1. 从一张空白的能带图说起：这套脚本补的是哪一段

正方晶格、介质柱、折射率差拉满，物理课本上三行字就能讲完的 2D 光子晶体带隙，真正落到电脑上往往卡在第一张能带图上：单胞怎么切、周期条件挂在哪两条边上、k 点该给弧度还是给无量纲、特征频率求解器为什么吐出一堆负数和复数。photonic-bandgap-comsol-matlab 这套 COMSOL 与 MATLAB 的联合脚本，价值不在算法有多新，而在于把 2D 光子晶体带隙计算这条链路上的建模、扫 k、判带隙三步固定成了可重复执行的代码。适合两类人：一是刚接触 photonic crystals 仿真、被 GUI 点选绕晕的研究生，二是需要批量跑参数、画带隙地图又不想每次重头点一遍的从业者。读完你应该能自己跑出 Γ-X-M-Γ 路径上的能带图，并知道哪几个参数一改结果就变脸。

2. 计算链拆解：COMSOL 单胞建模与 MATLAB 驱动的分工

先把这条链路的职责分清楚，后面写代码才不会一边建模一边怀疑物理。COMSOL 负责的是有限元本征值求解这件重活：离散化单胞、装配 Bloch 周期边界、解广义本征值方程；MATLAB 负责的是循环、判据、出图和落盘。两者之间走 LiveLink 接口，MATLAB 手里握着的是同一个 model 对象的句柄，改一个参数不需要存盘再导入。

2.1 为什么带隙计算要坚持脚本化而不是 GUI 点选

单点求解用 GUI 五分钟能点完，问题在于带隙从来不是一个点的事。一条 Γ-X-M-Γ 路径按每段 20 个点算就是 60 次本征值求解，再叠上 r/a 从 0.10 扫到 0.45，规模立刻变成几百次。GUI 在这种量级下有兩個硬伤：一是每次改参数都要重新点一遍几何、材料、网格，出错概率随次数线性上升；二是结果散落在各个节点里，事后想复原“这张图是哪个参数跑的”只能靠记忆。

脚本化解决的正是这两点。参数写在函数签名里，一次调用就是一次完整实验，想复现某张图只要把那次调用的参数抄回来。我在实际项目里的一般做法是：建模封成一个pc_build_cell函数，求解扫点封成另一个函数，中间的模型对象不销毁、只改属性，这样一次装配多次复用，速度能差出好几倍。

提示：LiveLink 的连接是有状态的，MATLAB 里ModelUtil.tags能看到当前挂着哪些模型，跑批之前先清一遍，避免拿到上一次的残留句柄。

这里还有个选型上的判断：如果你的目标是快速看个大概，COMSOL 自带的周期性结构案例模板改一改更快；但只要涉及参数扫描、带隙宽度统计、和平面波展开法交叉验证，脚本就是唯一不折腾的路。别指望 GUI 里存个 model 文件然后用 MATLAB 打开改参数，那条路的版本兼容坑比想象中多。

2.2 单胞里到底放了什么：几何、材料与 Bloch 相位

正方晶格的单胞就是一个边长 a 的正方形，中心放一根半径 r 的圆柱。这里有个容易忽略的几何前提：只有当 r/a < 0.5 时柱体才完全不触碰单胞边界，边界落在均匀背景介质里，Floquet 周期条件挂上去才干净。如果你的参数扫到 r/a = 0.48 附近，柱体已经贴边，网格在边界处的质量和周期配对都会变得微妙，带隙位置会开始抖——这不是物理，是离散误差。

材料上只需两组数：柱体相对介电常数 εr（常见取值 11.4 左右对应近红外波段的硅，或者 12 附近对应砷化镓），背景取 1（空气）。介质柱和空气孔的带隙位置差别很大，脚本里通过交换这两个数就能切换，不需要改几何。

Bloch 相位是整件事的物理核心。周期边界上场的相位差是 exp(i k·a)，COMSOL 里对应 Floquet 周期条件，需要给的是波矢分量 kx、ky，单位 rad/m，不是无量纲的 ka。这一点是新手最容易翻车的地方：给成 ka 会让整条能带图横向压缩 1/a 倍，看着像带隙消失了，其实是横轴错了。

几何与材料的落点整理成表更直观：

把 a 取成 1 是个省事的做法，所有长度都变成相对值，归一化频率直接就是 a/λ，不用每次除来除去；代价是材料参数里的长度量纲（比如色散模型）没法直接写进去。带隙计算通常只关心无色散的介电常数，取 a = 1 完全够用。

2.3 LiveLink 的调用骨架与对象树对应关系

COMSOL 的 MATLAB API 是照抄 GUI 的对象树的，理解这一点比背函数名重要。model.component('comp1')对应模型开发器里的组件节点，下面挂geom、mesh、physics、material、study、sol六个子树，路径写法和 GUI 里点开的层级一模一样。所以你在 GUI 里能点到的设置，API 里一定有对应的feature和set。

装配顺序有讲究：几何必须先run完再建物理场，否则物理场选不到域；材料要在物理场之后赋，不然默认材料会覆盖你的设置；网格必须最后建，因为它依赖几何的最终形态。这个顺序搞反了，报的错往往和真正的病因隔着十万八千里。

% 建立空模型并挂上 2D 组件 import com.comsol.model.* import com.comsol.model.util.* model = ModelUtil.create('Model'); % 新建模型，tag 可自定义 model.modelNode.create('comp1'); % 组件节点 model.geom.create('geom1', 2); % 二维几何 % 正方形单胞（左下角在原点）+ 中心圆柱 model.geom('geom1').create('r1', 'Rectangle'); model.geom('geom1').feature('r1').set('size', {'a' 'a'}); model.geom('geom1').feature('r1').set('pos', {'0' '0'}); model.geom('geom1').create('c1', 'Circle'); model.geom('geom1').feature('c1').set('r', 'r'); model.geom('geom1').feature('c1').set('pos', {'a/2' 'a/2'}); model.geom('geom1').run; % 几何必须显式 run，否则后续选不到域

上面这段只做了一件事：把单胞的拓扑定下来。size接受字符串是为了让参数以变量形式存在，后面改a、r的值不需要重建几何。pos对矩形是左下角坐标，对圆是圆心坐标，这两个语义不一样，写错一个就会得到偏心的单胞。

参数有两种给法：一种是在model.param.set('a','1')里登记全局参数，几何里引用名字；另一种是直接把数值塞进去。前者适合做参数扫描，后者适合一次性求解。我一般用前者，因为参数扫描时只需要改model.param，几何、网格、物理场会自动跟着变，省掉大量重建。

3. 跑通第一个带隙：环境对齐、参数标定与单点求解

单点跑通是后面所有工作的地基。这一章的目标很具体：定一个 k 点（比如 X 点），把本征频率解出来，确认前几条带的量级落在合理区间。如果这一步的频率已经离谱，后面的能带图全是白费功夫。

3.1 版本与接口对齐，避免第一行就报错

LiveLink 的连接方式取决于你的 COMSOL 是客户端模式还是服务端模式。桌面版装了 LiveLink for MATLAB 之后，MATLAB 里用mphstart起一个本地服务，端口号在 COMSOL 安装目录的配置文件里能查到；也可以从 COMSOL 那一侧启动带 MATLAB 的会话，这样 MATLAB 打开时连接已经建立好了。两种方式都行，前者的好处是脚本可以独立运行。

版本错配是新手第一个坑。MATLAB 的版本要比 COMSOL 支持的版本新或者持平，反过来容易在import com.comsol.model.*这一行就抛异常。另外 COMSOL 6.x 之后部分物理场字段名有调整，比如二维电磁波频域里关于极化方向的属性，老脚本抄过来会提示“属性不存在”。遇到这种情况不要硬猜，用mphgetproperties把当前对象的属性全打出来看：

% 打印物理场属性，确认字段名与当前版本一致 props = mphgetproperties(model.physics('ewfd').prop('Components')); disp(props);

这一步看着笨，但比在论坛上翻半小时帖子快得多。我现在的习惯是每换一个 COMSOL 大版本，先把关键节点的属性名扫一遍，把差异记在自己的小本子上，后面就能少流点血。

3.2 六个决定结果的物理参数

参数不是越多越好，真正影响带隙位置的其实只有六个，其余都是数值细节。下面这张表是我在多次扫描之后固定下来的检查清单，每次跑新结构先对一遍：

shift这一项特别值得说。特征频率求解器本质上是移位反迭代，你给的中心频率离真实本征值越近，收敛越快、越不容易漏解。给 0.35c/a 是因为正方晶格介质柱结构的第一带隙通常在归一化频率 0.3 到 0.5 之间，这个值覆盖了大多数常见填充比。如果你的结构是空气孔、介电常数对比反过来，带隙会往低频挪，shift 得跟着降。

网格尺寸hmax取 a/15 是个折中，跑得快但精度一般。真要出图，我一般用 a/20 起步，关键结论再用 a/30 复核一遍，看频率有没有收敛。这个复核步骤后面第 6 章会展开。

3.3 从建模到特征频率的完整脚本

下面这段是可以直接抄的：给定 a、r、εr 和一个 k 点，返回该点前若干条带的归一化频率。

function fNorm = solve_one_k(a, r, epsr, kx, ky, neigs) % 单个 k 点的本征频率求解 % 输入：a 晶格常数, r 柱半径, epsr 柱介电常数 % kx,ky 波矢分量(rad/m), neigs 求解模态数 % 输出：fNorm 归一化频率 (a/lambda)，已排序 import com.comsol.model.* import com.comsol.model.util.* c0 = 299792458; % 真空光速 model = ModelUtil.create('Model'); model.modelNode.create('comp1'); model.geom.create('geom1', 2); % --- 几何 --- model.geom('geom1').create('r1', 'Rectangle'); model.geom('geom1').feature('r1').set('size', {a, a}); model.geom('geom1').feature('r1').set('pos', {'0' '0'}); model.geom('geom1').create('c1', 'Circle'); model.geom('geom1').feature('c1').set('r', r); model.geom('geom1').feature('c1').set('pos', {a/2, a/2}); model.geom('geom1').run; % --- 材料：柱体与背景分开赋 --- model.component('comp1').material.create('mat1', 'Common'); model.component('comp1').material('mat1').selection.set([1]); % 圆域 model.component('comp1').material('mat1').propertyGroup('def')... .set('relpermittivity', {epsr}); model.component('comp1').material.create('mat2', 'Common'); model.component('comp1').material('mat2').selection.set([2]); % 背景域 model.component('comp1').material('mat2').propertyGroup('def')... .set('relpermittivity', {1}); % --- 物理场：二维电磁波频域，出平面极化 --- model.physics.create('ewfd', 'ElectromagneticWavesFrequencyDomain', 'geom1'); model.physics('ewfd').selection.all; model.physics('ewfd').prop('Components').set('components', 'outofplane'); % --- 周期边界：左右一对、上下各一对 --- model.physics('ewfd').feature.create('pc1', 'PeriodicCondition', 2); model.physics('ewfd').feature('pc1').selection.set([1 3]); % 左右边 model.physics('ewfd').feature('pc1').set('PeriodicType', 'Floquet'); model.physics('ewfd').feature('pc1').set('kFloquet', {kx, 0}); model.physics('ewfd').feature.create('pc2', 'PeriodicCondition', 2); model.physics('ewfd').feature('pc2').selection.set([2 4]); % 上下边 model.physics('ewfd').feature('pc2').set('PeriodicType', 'Floquet'); model.physics('ewfd').feature('pc2').set('kFloquet', {0, ky}); % --- 网格 --- model.mesh.create('mesh1', 'geom1'); model.mesh('mesh1').create('ftri1', 'FreeTri'); model.mesh('mesh1').feature('size').set('custom', 'on'); model.mesh('mesh1').feature('size').set('hmax', a/15); model.mesh('mesh1').run; % --- 特征频率研究 --- model.study.create('std1'); model.study('std1').create('eig', 'Eigenfrequency'); model.study('std1').feature('eig').set('neigs', neigs); model.study('std1').feature('eig').set('shift', 2*pi*0.35*c0/a); model.study('std1').run; % --- 取解并归一化 --- freq = mphglobal(model, 'freq', 'dataset', 'dset1'); % 单位 Hz fNorm = sort(real(freq)) * a / c0; % 归一化 a/lambda end

逐段说清楚几个关键点。几何部分的selection.set([1])依赖域的编号，这个编号在不同几何构造顺序下会变，稳妥做法是用坐标箱选域（mphselectbox），而不是写死 1 和 2。周期边界那两行里的selection.set([1 3])同理，边界编号也依赖几何，第一次跑的时候务必先用mphgetadj或者 GUI 里点一下确认哪条边是 1 号。

shift里乘了 2π，因为 COMSOL 的角频率和普通频率差一个 2π，这个转换漏掉会让搜索中心偏到十万八千里，表现为求解器反复报“找不到特征值”。mphglobal取出来的freq是普通频率（Hz），归一化除以 c/a 得到 a/λ，这就是文献里最常用的横轴定义。

调用方式和预期输出：

% 在 X 点 (pi/a, 0) 求解前 8 条带 fx = solve_one_k(1, 0.2, 11.4, pi, 0, 8); disp(fx');

r/a = 0.2、εr = 11.4的正方晶格介质柱，X 点前几条归一化频率大致落在 0.1 到 0.6 的区间里。如果你拿到的是 3.14 这样的数，八成是忘了除 c/a；如果是负数或者带虚部，八成是 shift 偏了或者网格太粗。

4. 沿高对称路径扫描：从单点本征值到能带图

单点跑通只说明模型没错，带隙有没有、多宽，得看整条路径。这一章把单点循环成路径，再把路径结果折算成带隙判据，最后落到文件上。

4.1 归一化频率与带隙判据怎么写

能带图的横轴是沿路径累积的 k 距离，纵轴是归一化频率 a/λ。带隙的判定就是逐条带比较：第 n 条带在整条路径上的最大值，是否小于第 n+1 条带的最小值。这个判据朴素但容易写错——错在拿单个 k 点的值比，或者忘了带隙要在整条路径上成立。

function [gapW, gapMid, idx] = bandgap_from_bands(fNorm) % fNorm: 能带数 x k点数 的归一化频率矩阵，已按带排序 % gapW : 相对带隙宽度 (f_upper - f_lower) / f_mid % gapMid: 带隙中心归一化频率 % idx : 带隙位于第 idx 与 idx+1 条带之间 [nb, ~] = size(fNorm); gapW = 0; gapMid = NaN; idx = 0; for n = 1:nb-1 fLow = max(fNorm(n, :)); % 下带在整个路径上的最高点 fHigh = min(fNorm(n+1, :)); % 上带在整个路径上的最低点 if fHigh > fLow % 严格重叠才叫带隙 w = (fHigh - fLow) / ((fHigh + fLow)/2); if w > gapW gapW = w; gapMid = (fHigh + fLow)/2; idx = n; end end end end

这里的相对宽度用的是常见定义Δf / f_mid，有些文献用Δf / f_center里的中心换成带隙中点，两者差别不大但要前后统一。fLow取 max、fHigh取 min，这两个方向写反了会得到“处处有带隙”的荒谬结果，我见过不止一次。

还要注意一点：带隙只有在采样足够密的时候才判得准。路径上只取 20 个点时，带的极值可能落在两个采样点之间，算出来的宽度偏小。经验值是每段至少 30 个点，做正式结论时用 50 个点复核。

4.2 k 点路径生成与批量求解循环

正方晶格的高对称路径是 Γ(0,0) → X(π/a,0) → M(π/a,π/a) → Γ(0,0)，三段长度不同，采样点数按段分配更均匀。

function [kpath, kdist] = kpath_gxm(a, nSeg) % 生成 Gamma-X-M-Gamma 路径，每段 nSeg 个点 t = linspace(0, 1, nSeg)'; kx = [pi/a*t; pi/a*ones(nSeg,1); pi/a*(1-t)]; ky = [zeros(nSeg,1); pi/a*t; pi/a*(1-t)]; kpath = [kx ky]; % 累积距离作为横轴 dk = sqrt(sum(diff(kpath).^2, 2)); kdist = [0; cumsum(dk)]; end

kpath是三段拼起来的 N×2 矩阵，kdist是沿路径的累积长度，画图时用它当横轴，X 点和 M 点的位置才对得上。三段拼接处首尾点重复了一次，实际使用时可以去重，不影响结果。

批量求解有两条路。稳妥的一条是老老实实循环，每次改kFloquet再run，缺点是每次都要重新装配矩阵。快的一条是在研究节点里加参数扫描，把 kx、ky 登记成全局参数，一次提交跑完全部 k 点，COMSOL 内部会复用装配结果。参数一大，后者能快出三到五倍。

% 循环版：直观但慢，适合调试 nk = size(kpath, 1); fAll = zeros(8, nk); for i = 1:nk fAll(:, i) = solve_one_k(a, r, epsr, kpath(i,1), kpath(i,2), 8); fprintf('k 点 %d/%d done\n', i, nk); end

这里每次调用solve_one_k都会新建一个模型，跑 90 个点就是 90 次完整建模，时间全花在重复劳动上。正式跑之前把这段换成“建一次模型、循环改 kFloquet、循环 run”的写法，效率差别非常明显。

4.3 结果落盘与可复现性校验

结果不落盘等于白跑。我一般存两个文件：一个是原始矩阵（带数 × k 点数），一个是元信息（a、r、εr、网格尺寸、COMSOL 版本）。后者看着多余，但三个月后想复现某张图时，它是唯一的后悔药。

% 保存原始数据与元信息 meta = struct('a', a, 'r', r, 'epsr', epsr, ... 'hmax', a/15, 'neigs', 8, ... 'comsol', char(mphversion), 'stamp', datestr(now)); save('bands_r020.mat', 'fAll', 'kpath', 'kdist', 'meta'); writematrix([kdist, fAll'], 'bands_r020.csv');

mphversion返回当前 COMSOL 版本字符串，char转一下才能存进 struct。CSV 里第一列是 k 距离，后面每列一条带，这样用别的工具画图也不用再解析 mat 文件。

校验的做法很简单：把同一组参数跑两遍，比较矩阵是否逐元素相等。如果不等，说明模型里有随机因素（比如网格自适应），需要固定网格或者关掉自适应。这个检查只要十秒，能省掉后面一堆“这张图怎么和上次不一样”的困惑。

5. 避坑与排查：五类让带隙凭空消失的现场

下面这五条都是我在实际跑批时遇到的，按现象、原因、解决三段写，遇到类似症状可以直接对号入座。

5.1 几何与单位类的翻车

现象一：能带图上一片平坦，所有 k 点的频率几乎一样，看不到任何色散。原因：周期边界挂到了相邻的两条边上，比如把左边和下边配成一对。Floquet 条件要求配对的两条边是平移关系，配错了等于强行给场加了一个错误的相位约束，色散关系直接被抹平。另一种可能是kFloquet给的是无量纲的 ka 而不是 kx，量级差 1/a 倍，路径上所有点挤在一起。 解决：用坐标箱显式选边，别依赖编号。左边是 x=0、右边是 x=a，上下同理，选完用mphgetadj打印一下确认配对。kFloquet的单位统一用 rad/m，a 取 1 的时候 kx 最大是 π，这个数量级一眼能看出来对不对。

现象二：改了 r/a 之后带隙位置跳变，相邻两个填充比的结果对不上。原因：几何用了参数化但没有重新run，或者圆和矩形做布尔运算时因为 r/a 接近 0.5 出现了退化边。r/a = 0.45 时柱体离边界只剩 0.05a，网格在边角处质量骤降。 解决：每次改参数后强制model.geom('geom1').run，并在 r/a > 0.4 时把边界附近的网格单独加密，或者直接把这个区间排除在结论之外。

现象三：归一化频率算出来是 3 左右，和文献差了一个数量级。原因：归一化时用了角频率而不是普通频率，或者忘了乘 a。mphglobal取freq是 Hz，取omega是 rad/s，两者差 2π。 解决：统一在脚本里写一个normalize小函数，内部固定f * a / c0，所有出口都走它，杜绝各处手写。

5.2 求解器与后处理类的翻车

现象四：特征频率求解器报“未找到特征值”，或者返回一堆负数和复数。原因：shift给的搜索中心离真实本征值太远，移位反迭代收敛到了错误的谱段；或者网格过粗导致高频伪解混进来。 解决：先把 shift 设到目标频段中部，用2*pi*0.35*c/a这个量级试；如果还不行，把neigs调大，让求解器多吐几条，看看真实的本征值到底在哪，再回头设 shift。负数频率是求解器在找代数解，物理上没意义，直接丢掉就行。

现象五：一个 k 点要跑十几分钟，九十点路径跑到天亮。原因：默认求解器用的是直接求解器，每个 k 点都重新装配、重新分解一次矩阵；另外每次循环都新建模型，几何、网格、材料全部重建。 解决：建模和求解分离，模型建一次，循环里只改kFloquet和shift；更进一步用参数扫描，让 COMSOL 在一次提交里跑完所有 k 点。实测下来，同样的 90 个点，从四小时压到四十分钟是常态。

6. 进阶校验：拿平面波展开法给有限元结果上保险

有限元的结果好看，但你得知道它可不可信。收敛性检查是一种办法，交叉验证是另一种，两种都做一遍，出图的时候心里才踏实。

6.1 网格与 k 点收敛的量化做法

收敛检查不是“加密一次看看差不多就行”，要有量化指标。做法是固定 k 点路径，把hmax从 a/10 依次取到 a/30，记录同一个带隙的上下边界频率，看相对变化。

hList = [1/10 1/15 1/20 1/30]; res = zeros(numel(hList), 2); for i = 1:numel(hList) [fAll, ~] = run_bands(a, r, epsr, hList(i), 8); [~, ~, ~] = deal(0); [~, fLow, fHigh] = gap_bounds(fAll); % 返回带隙上下边界 res(i, :) = [fLow fHigh]; fprintf('hmax/a = 1/%d : %.4f ~ %.4f\n', round(1/hList(i)), fLow, fHigh); end

判据是相邻两档的带隙边界相对差小于 1%。正方晶格介质柱在 r/a = 0.2 附近，a/15 和 a/20 之间通常差 0.5% 以内，a/20 之后基本不再动。如果到 a/30 还在漂，多半是柱体边界网格没有贴合曲率，加一层边界层网格再试。

k 点收敛同理，把每段采样从 20 加到 50，看带隙宽度变化。经验上 30 点之后变化就很小了。

6.2 一段三十行的 PWE 交叉验证代码

平面波展开法（PWE）是光子晶体带隙的经典方法，思路是把场展开成倒格矢上的平面波，把微分方程变成矩阵本征值问题。它和有限元是完全独立的数值路径，两者结果接近，才说明你的有限元模型没有隐藏问题。

function [wTE, wTM] = pwe_square(a, r, epsr, kx, ky, Nmax) % 正方晶格介质柱的 PWE 本征频率（归一化 a/lambda） % Nmax 控制平面波截断，总平面波数 (2*Nmax+1)^2 n = -Nmax:Nmax; [Nx, Ny] = meshgrid(n, n); G = 2*pi/a * [Nx(:) Ny(:)]; % 倒格矢集合 NG = size(G, 1); f = pi*r^2/a^2; % 填充比 % 傅里叶系数矩阵（圆形柱的解析式，含贝塞尔函数） Cinv = zeros(NG); Ceps = zeros(NG); % 1/eps 与 eps 的系数 c0i = f/epsr + (1-f); c1i = (1/epsr - 1)*f; c0e = f*epsr + (1-f); c1e = (epsr - 1)*f; for i = 1:NG for j = 1:NG dG = norm(G(i,:) - G(j,:)); if dG < 1e-12 Cinv(i,j) = c0i; Ceps(i,j) = c0e; else J = 2*besselj(1, dG*r)/(dG*r); Cinv(i,j) = c1i*J; Ceps(i,j) = c1e*J; end end end % 动能矩阵：(k+G)·(k+G') K = G + [kx ky]; T = K * K'; % E 平行柱（TE）用 1/eps 系数，H 平行柱（TM）用 eps 系数 evTE = eig(T .* Cinv); evTM = eig(T .* Ceps); wTE = sort(sqrt(abs(real(evTE)))) * a / (2*pi); wTM = sort(sqrt(abs(real(evTM)))) * a / (2*pi); end

逐块说明。倒格矢集合用meshgrid生成二维网格再摊平，Nmax = 6时共 169 个平面波，矩阵 169×169，eig几毫秒就出结果。傅里叶系数那块是整套代码的核心：圆形柱的介电常数分布在倒空间里有解析形式，2*J1(Gr)/(Gr)这个因子在 G 趋近 0 时极限为 1，所以要单独处理dG < 1e-12的分支，否则会除零。

T .* Cinv是逐元素乘，不是矩阵乘，这一点写错会得到完全无意义的本征值。最后的归一化sqrt(ev)*a/(2*pi)，ev的量纲是 (ω/c)²，开方后乘 a 再除 2π 才得到 a/λ，和前面有限元的定义对齐。

交叉验证的做法是取同一个 k 点，比较两边的低阶带：

上面这组数只是某次运行的样例，具体数值随参数变化，看的是偏差量级。低于 1% 就说明两套方法在同一个物理图像上，有限元这边没有系统性的建模错误。如果偏差超过 3%，先查 PWE 的Nmax够不够（加到 10 再看），再查有限元的周期边界和极化方向是不是设反了。

PWE 也有它的边界：介电常数对比特别大、或者柱体形状不规则（六边形、椭圆）时，解析傅里叶系数就没有了，得换成数值积分，收敛会慢不少。所以我的定位是——PWE 当交叉验证的参照，不当主力求解器，主力还是有限元，因为有限元能处理任意几何和色散材料。

注意：两套方法在带隙边缘附近偏差会放大，这是因为带隙边缘的态密度高、对离散化误差敏感，比较时优先看带隙中心的带，别看边缘。

从那以后我每次出一张能带图，不管多赶时间，都会强制走一遍“PWE 对照 + 网格收敛”这两步，哪怕只是把hmax从 a/15 换到 a/20 看一眼。这个习惯救过我不止一次——有一次带隙宽度算出来比预期宽了将近一倍，就是靠 PWE 对照发现有限元的周期边界挂错了一条边，改成坐标箱选边之后两条曲线立刻贴上了。希望帮到你。

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