微环谐振腔光频梳LLE方程MATLAB仿真入门与实战

微环谐振腔的光频梳仿真，这几年是集成光子学里相当热的方向。从片上光互连、微波光子学到精密测量，微梳都算得上一个绕不开的基础工具。但刚接触这个方向时，最大的门槛不是物理概念，而是怎么把一堆偏微分方程变成能跑的代码。LLE方程（Lugiato-Lefever Equation）是目前做微腔光频梳仿真最主流的模型，能把微环里的泵浦、损耗、色散、非线性这几件事统一到一个方程里，而且用MATLAB就能在普通笔记本上跑出漂亮的梳状谱。这篇博文就围绕LLE方程的MATLAB仿真展开，从方程理解、参数归一化、代码实现到结果解读和排坑技巧，完整过一遍。

如果你正准备做微腔光频梳方向的课题，或者手头有微环器件但想先通过仿真理解非线性动力学过程，这篇文章可以直接当一份入门操作手册用。我不只写能跑的代码，还会解释每一步为什么这么做、参数怎么选取、哪些坑容易踩，以及如何判断仿真结果是否物理可信。

1. 光频梳与微环谐振腔：先想清楚我们在仿什么

1.1 微环里怎么长出一把“梳子”

光频梳的本质是一系列频率等间隔、相位相干的光谱线。锁模激光器能产生频梳，微环谐振腔也能，而且微环的结构更紧凑、功耗更低。微环产生频梳的物理过程可以这样理解：一束连续波激光从直波导通过倏逝波耦合进入微环，当激光频率与微环某个谐振模式对齐时，腔内的光场被增强，功率密度远高于输入波导中的平均值。

增强的光场通过克尔非线性效应改变材料折射率，激发四波混频过程，产生两个新的频率分量（信号光和闲频光）。这两个新分量分布在泵浦两侧，与泵浦的频率间隔等于自由光谱范围。随着腔内功率继续提升，新产生的频率分量又会作为新的“泵浦”参与更高阶的四波混频，不断有频率位于泵浦两侧、间隔为FSR整数倍的梳齿被激发，最终形成跨越一定光谱范围的等间隔频率线。整个过程实质上是一个级联的参量过程，在频域上看到的是一把逐渐铺开的梳子。

需要特别注意的是，微环产生频梳并不是必然的。腔的色散、损耗、泵浦失谐、泵浦功率都直接影响梳齿能否建立、梳齿包络是什么形状。仿真的核心就是在数值上重现这些条件的相互作用，预测一个具体的微环参数体系会在什么状态下运行。

1.2 为什么仿真选择LLE方程而不是其他模型

做微环仿真，可选的方法不止一种。用FDTD或FEM直接做全波仿真，理论上很精确，但计算资源需求非常夸张。一个半径100微米的微环，光学波长在1550纳米附近，空间网格要覆盖的尺度跨越好几个数量级，三维FDTD单次仿真就可能要跑几周甚至更久，基本不实用。耦合模理论（CMT）思路直观，把每个谐振模式当作一个微分方程来描述，但模式多了方程数量呈指数增长，而且模式之间的非线性串扰处理起来很繁琐。

LLE方程走的是另一条路：它建立在慢变包络近似之上，把光场在腔内往返一周的快变化平均掉，用一个在慢时间尺度上演化的方程描述整个腔场分布。这个思路相当于从“每个模式的独立方程”变成“全场分布的一个方程”，复杂度大幅下降，同时保留了色散、损耗、非线性、泵浦等核心物理。这也是为什么LLE已经成为微腔频梳仿真的事实标准，大量顶刊文献中的微梳仿真结果都是基于LLE方程得到的。

当然，LLE方程有适用条件：腔内的非线性效应不能太强以至于高次谐波和模式耦合超出慢变近似范围；微环最好是单模传输或者至少我们要关注的谱段内没有显著的模式交叉；泵浦线宽足够窄，可以近似为单频。对于绝大多数硅基、氮化硅微环频梳研究来说，这些条件都是满足的。还有一点：LLE是非线性偏微分方程，不存在通用的解析解，必须依靠数值方法求解，而这正是MATLAB擅长的事情。

2. 仿真整体设计与归一化参数体系

2.1 归一化到底在归一什么

做微环LLE仿真，第一件绕不开的事情就是参数归一化。物理参数一列出来会很吓人：微环半径、有效折射率、群折射率、FSR、损耗率、耦合率、色散系数、非线性折射率、有效模场面积、泵浦功率、失谐频率……直接把这些参数塞进偏微分方程里，数值范围差异巨大，不仅迭代步长不好选，数值稳定性也容易出问题。归一化的意义是把这些物理参数压缩成几个无量纲的数，让方程变得简洁，避免量纲带来的麻烦。

归一化形式的LLE方程通常写成这样：

$$\frac{\partial \psi}{\partial \tau} = \left[ -(1 + i\alpha) + i|\psi|^2 + \frac{i d_2}{2}\frac{\partial^2}{\partial \theta^2} \right] \psi + f$$

这里的$\psi$是归一化腔场，$\tau$是归一化慢时间，$\theta$是绕环一周的角坐标，$\alpha$是归一化泵浦失谐，$d_2$是归一化二阶色散，$f$是归一化泵浦强度。物理参数到归一化参数的映射关系，不同文献写法略有差异，但主流约定大致如下：

• 时间以总损耗率的倒数归一化：$\tau = \kappa t$，$\kappa$是微环总损耗率（包括耦合损耗和固有损耗）。
• 场幅以腔内非线性系数和损耗归一化，使得非线性项正好变成$i|\psi|^2\psi$。
• 失谐$\alpha$是泵浦频率相对谐振频率的偏移量，单位是损耗率。
• 泵浦项$f$的平方正比于输入泵浦功率和耦合效率，反比于损耗率的三次方。
• 色散系数$d_2$是二阶色散相对损耗率做归一化后的结果。

举个具体换算例子。一个氮化硅微环，半径100微米，自由光谱范围约200GHz，总损耗率$\kappa = 2\pi \times 1\text{ GHz}$，泵浦失谐$\delta = 2\pi \times 2.5\text{ GHz}$，那么归一化失谐$\alpha = \delta / \kappa = 2.5$。耦合效率$\eta=0.5$，泵浦功率20毫瓦，非线性系数$\gamma$约1 W⁻¹m⁻¹，腔长$L=628$微米，算出来的$f$大致在4到6这个量级。这个量级感对后面调参数特别有用：如果仿真时发现$f$值只有0.1，那意味着泵浦功率太低，梳齿大概率建不起来；如果$f$动辄几十上百，那要小心数值发散，因为非线性太强。

我建议开始仿真时不要直接换算具体器件参数，先用一组已知能出梳齿的归一化参数（比如$\alpha=2.5$，$f=4.5$，$d_2=-0.02$）实现整条链路，然后再换成自己的目标参数。这样能快速验证代码是否正确，也能逐步积累参数区间的感觉。

2.2 数值方法选型：为什么用分步傅里叶法

LLE方程是一个含一阶时间导数、二阶空间（角）导数和非线性项的偏微分方程。MATLAB求解偏微分方程有很多方法，最简单的是有限差分法，但要求空间网格足够密，时间步长足够小，否则数值色散和数值耗散会污染结果。更好的选择是分步傅里叶法（Split-Step Fourier Method），这也是非线性光纤光学里一种非常成熟的做法。

分步傅里叶法的核心思想是：把整个微分算子拆成线性部分和非线性部分，每一步让光场分别作用。线性部分包括损耗、失谐、色散，可以在频域精确计算；非线性部分只包含$i|\psi|^2\psi$，在时域直接做乘法。傅里叶变换把空间微分变成频域代数乘法，这让线性步从微分方程变成了一个指数算子相乘，既快速又稳定。非线性步在时域做乘法，物理意义就是克尔效应引起的自相位调制。

具体实现时，每步迭代分为三步：先走半个线性步（频域），再走一个完整非线性步（时域），最后再走半个线性步（频域）。这种对称结构叫对称分步傅里叶法，比简单的一步线性一步非线性精度要高一个阶次，在同样的时间步长下误差更小。泵浦项$f$是一个常数驱动项，不依赖于当前场值，按常规做法在每个小步累加$f \cdot \Delta\tau$即可。这里有个细节值得说：有些文献会把泵浦项和损耗失谐项一起放进线性算子做“推进”，但因为这涉及常数项，处理不仔细会导致泵浦驱动次数加倍，所以最好还是显式累加，逻辑清楚也不容易出错。

网格和时间步长的选择是仿真成败的关键。角坐标$\theta$通常取$[-\pi, \pi)$或$[0, 2\pi)$范围，离散点数$N$取256到1024之间。频域中最高波矢对应的频率必须覆盖我们关心的谱宽，否则高频分量会被截断或混叠。时间步长$\Delta\tau$一般取$10^{-3}$到$10^{-5}$量级，具体取决于泵浦强度和非线性强度；跑通之前可以先用小步长验证，再逐步放大到可接受的计算速度。

2.3 把泵浦、损耗、色散各项准确放进迭代

写代码前，先把迭代逻辑整清楚。归一化LLE方程右边有几类项：线性项中的$- \psi$对应损耗，$-i\alpha\psi$对应失谐，$\frac{i d_2}{2} \frac{\partial^2 \psi}{\partial \theta^2}$对应色散；非线性项$i|\psi|^2\psi$对应克尔效应；常数项$f$对应泵浦。分步傅里叶法的每次迭代，相当于把这几个物理过程按顺序作用一次。

频域线性算子的形式是：

$$L(k) = -1 - i\alpha + \frac{i d_2}{2} k^2$$

其中$k$是角坐标对应的波矢序号。注意这个算子不依赖于当前场值，所以可以在循环之前一次性构建好，每次迭代直接用指数算子的形式相乘即可。这一步相当于把损耗、失谐、色散这三个线性过程“一步到位”地处理了。

非线性步的写法是直接乘一个相位因子：$\psi \leftarrow \psi \cdot \exp(i|\psi|^2 \Delta\tau)$。这对应的是在一个小时间步内，克尔非线性引起的瞬时相位改变。如果$|\psi|^2$比较大，这个相位变化也会比较大，所以时间步长不能取得太大，否则相当于用矩形近似积分非线性效应，精度不够。

泵浦项$f$是常数，直接在每一步加$f \cdot \Delta\tau$。迭代的初值一般取极小幅度的复高斯随机场，用来模拟腔内噪声和真空涨落，这些微扰在增益条件下会被放大，最终演化出梳齿。如果初值取零，泵浦项会不断注入能量，本质上也是一种驱动方式，但随机初值更接近实验中的实际启动过程。

3. 从零搭建MATLAB仿真代码

3.1 参数初始化与频域算子构建

直接上代码。下面这段MATLAB脚本实现了归一化LLE方程的分步傅里叶求解。参数变量命名尽量贴近方程符号，方便对照文中公式检查。

%% 归一化LLE方程频梳仿真 clear; close all; clc; %% 1. 归一化参数 alpha = 2.5; % 泵浦失谐（单位：损耗率） f_pump = 4.5; % 归一化泵浦强度 d2 = -0.02; % 归一化二阶色散，负值为反常色散 N = 512; % 角坐标离散点数 dt = 5e-4; % 归一化时间步长 t_max = 400; % 总仿真时间 n_steps = round(t_max / dt); %% 2. 角坐标网格与频域波矢 theta = linspace(-pi, pi - 2*pi/N, N).'; dk = 1; % 角波矢间隔，周期长度2*pi下每阶间隔为1 k = [0:N/2, -N/2+1:-1].'; k2 = k.^2; %% 3. 频域线性算子 linear_op = -1 - 1i*alpha + 1i * d2/2 * k2; %% 4. 初始场：复高斯微扰 rng(42); psi = 1e-5 * (randn(N,1) + 1i*randn(N,1)); %% 5. 预分配记录数组 power_record = zeros(n_steps, 1); %% 6. 主迭代：对称分步傅里叶 tic; for step = 1:n_steps % 半个线性步（频域） psi_hat = fft(psi); psi_hat = psi_hat .* exp(linear_op * dt/2); psi = ifft(psi_hat); % 非线性步（时域，自相位调制） psi = psi .* exp(1i * abs(psi).^2 * dt); % 半个线性步（频域） psi_hat = fft(psi); psi_hat = psi_hat .* exp(linear_op * dt/2); psi = ifft(psi_hat); % 泵浦项驱动 psi = psi + f_pump * dt; % 记录腔内平均功率，用于观察演化 power_record(step) = mean(abs(psi).^2); end toc;

这段代码跑在普通笔记本上，N=512、t_max=400、dt=5e-4，总共80万步，大约需要几十秒到一两分钟。如果嫌太慢，可以先把t_max降到100验证链路，再逐步增加。

细节上有个地方容易踩坑：fft和ifft的归一化问题。MATLAB的fft不做归一化，ifft会自动除以N。分步傅里叶中，频域乘一个指数算子然后ifft回来，整体幅值关系是对的，不需要额外处理。但如果自己手动实现FFT，一定要记得归一化，否则场幅会差一个因子N，非线性项的强度全变了。

3.2 主迭代循环里的工程优化

上面的代码能跑，但有几个工程优化点值得说。

第一点是数据布局。MATLAB是列优先存储，fft对列向量的处理效率高于行向量。上面代码中theta定义时用了转置，确保psi是列向量，这是有意的。如果数据是行向量，循环里反复转置会额外开销不少。

第二点是预分配。power_record在循环前就分配好了，避免了每次迭代动态增长数组，这是MATLAB性能优化里最基本的一条。凡是需要在循环中逐步记录的数据，都先zeros一条。

第三点是频域算子的预计算。linear_op在循环外一次性构建，避免了循环内反复计算k2和指数。如果要做失谐扫描，线性算子会随$\alpha$变化，但这时没必要每步重建，而是在每个扫描台阶更新一次即可。

第四点是关于对称分步格式的精度。对称分步傅里叶法的局部误差是三阶，全局误差是二阶。一般情况下足够。如果发现结果对dt敏感，可以把dt减半对比一次，若结果差异明显说明dt还没收敛，需要继续减小。

我对上面代码做过一个简单测试：同样的参数，dt取1e-3和5e-4，最终频梳形状几乎一致；dt取5e-3时就会出现明显的高频数值振荡，这说明该参数下dt的收敛上限在1e-3附近。跑自己参数时，建议用这种“减半步长对比结果”的方式确认收敛性。

3.3 结果提取与画图

仿真输出的psi是最后的腔场分布，还需要转换成可读的频梳信息。提取频梳功率谱的代码如下：

%% 7. 频域结果提取 psi_f = fftshift(fft(psi)); power_comb = abs(psi_f).^2 / N; % 归一化功率 freq_index = k; % 对应梳齿阶数 %% 8. 画图：时域与频域 figure('Position', [100 100 900 700]); % 时域强度分布 subplot(2,2,1); plot(theta/pi, abs(psi).^2, 'b-', 'LineWidth', 1); xlabel('角度 \theta / \pi'); ylabel('腔内功率 |\psi|^2'); title('时域腔内强度分布'); grid on; % 频域梳齿功率线性坐标 subplot(2,2,2); plot(freq_index, power_comb, 'r-', 'LineWidth', 1); xlabel('梳齿阶数 k'); ylabel('梳齿功率'); title('频率梳频谱（线性坐标）'); xlim([-60 60]); grid on; % 频域梳齿功率对数坐标 subplot(2,2,3); semilogy(freq_index, power_comb, 'k-', 'LineWidth', 1); xlabel('梳齿阶数 k'); ylabel('梳齿功率（对数）'); title('频率梳频谱（对数坐标）'); xlim([-100 100]); grid on; % 腔内平均功率演化 subplot(2,2,4); plot((1:n_steps)*dt, power_record, 'g-', 'LineWidth', 1); xlabel('归一化时间 \tau'); ylabel('平均腔内功率'); title('腔内功率演化'); xlim([0 t_max]); grid on;

看结果时，先看时域图：如果腔内功率呈现多峰结构，说明产生了调制不稳定或孤子晶格；如果分布平坦，说明处于连续波状态。再看频域图：线性坐标看主梳齿的分布，对数坐标看低功率梳齿能否分辨出来。把梳齿功率导出为数组后，如果想看某个梳齿的功率随参数的变化，可以在循环里按参数扫描并记录该阶数的功率值。

需要提醒的是，fftshift把零频移到中间后，横坐标freq_index的正负对应泵浦两侧的高低频方向。画图时xlim的范围要根据梳齿实际展宽适当调整，不要一开始就限制太窄漏掉重要结构。

4. 仿真结果如何看懂：从调制不稳定到孤子

4.1 从噪声出发，固定参数下会得到什么

固定失谐和泵浦功率，从随机噪声出发，系统演化到稳态后会出现几种典型状态。第一种是连续波状态：腔内场强力分布平坦，频域只有泵浦峰和几个微弱边带。这对应泵浦功率太低或失谐太大，四波混频增益不足以克服腔损耗。第二种是调制不稳定态：时域出现周期性强度调制，频域泵浦两侧出现对称的调制边带，这些边带继续级联产生更多梳齿。第三种是混沌态：时域波形完全不规则，频域梳齿数量多但包络不平滑，功率分布抖动明显。第四种是孤子态：时域出现一个个孤立的窄脉冲，频域梳齿的功率包络在双对数坐标下接近双曲正割形态，形状平滑。

作为一个经验法则，在反常色散（$d_2<0$）条件下，要激发调制不稳定性，归一化失谐$\alpha$和泵浦$f$都要达到一定阈值。粗略地说，在$f$大约3到6、$\alpha$大约1到5的范围内，比较容易观察到明显的梳齿建立过程。我自己跑过的经验是$f=4.5$、$\alpha=2.5$、$d_2=-0.02$这组参数，从噪声出发后约跑200到300个归一化时间单位，频域会出现清晰的梳齿结构，时域呈现周期脉冲序列。如果想把某个状态归类，可以结合时域、频域和功率演化的特征一起判断，单看一个域容易误判。

4.2 失谐扫描法找孤子：为什么需要扫描

实验里产生孤子微梳，通常靠慢速扫描泵浦激光频率，从蓝失谐扫到红失谐，让系统经历复杂的多稳态动力学，最终捕获在孤子态上。仿真里也可以复现这个过程，做法是将$\alpha$从一个较大正值逐渐减小，让系统在每个失谐值下演化一段时间，记录最终状态。

下面这段代码实现失谐扫描：

%% 失谐扫描找孤子 alpha_start = 5; alpha_end = 1; alpha_steps = 800; alpha_list = linspace(alpha_start, alpha_end, alpha_steps); psi = 1e-5 * (randn(N,1) + 1i*randn(N,1)); power_at_alpha = zeros(alpha_steps, 1); % 每个失谐台阶的演化时间 inner_steps = 200; dt_inner = 1e-3; for ai = 1:alpha_steps alpha = alpha_list(ai); linear_op = -1 - 1i*alpha + 1i * d2/2 * k2; for step = 1:inner_steps psi_hat = fft(psi); psi_hat = psi_hat .* exp(linear_op * dt_inner/2); psi = ifft(psi_hat); psi = psi .* exp(1i * abs(psi).^2 * dt_inner); psi_hat = fft(psi); psi_hat = psi_hat .* exp(linear_op * dt_inner/2); psi = ifft(psi_hat); psi = psi + f_pump * dt_inner; end power_at_alpha(ai) = mean(abs(psi).^2); end %% 画失谐扫描迹线 figure; plot(alpha_list, power_at_alpha, 'b.-'); xlabel('归一化失谐 \alpha'); ylabel('腔内平均功率'); title('失谐扫描迹线'); grid on;

扫描方向要特别注意：实验中孤子通常在失谐减小（从蓝失谐到红失谐）过程中被捕获，仿真也沿着同样的方向扫描。如果反向扫描，系统通常会停在调制不稳定或混沌态上，很难进入孤子分支。扫描速度也很关键：每步失谐变化量要足够小，相当于绝热演化，太快了系统来不及响应。上面的参数每步$\Delta\alpha \approx 0.005$，配合200步内部演化，总体扫描耗时约16万个归一化时间单位，计算量不小。

如果扫描完发现腔内平均功率迹线出现明显台阶，说明系统在失谐扫描中发生了状态跳变，其中一个台阶很可能就是孤子分支。这时保存当前$\psi$，画时域和频域图确认一下。孤子的时域特征是单个或少数几个窄峰，频域包络平滑。如果看到一堆密集而不规则的窄峰，那多半是混沌态的局部片段，不是孤子。

4.3 典型仿真结果解读与物理判断

以$f=4.5$、$d_2=-0.02$、扫描$\alpha$从5到1这组条件为例，扫描迹线通常会是这样的走势：大失谐处功率接近线性增长段，对应连续波解；随着失谐减小，功率出现折叠和双稳结构；在某个临界失谐附近，功率发生跳变，从高功率连续波态跳到低功率孤子态。孤子态下腔内平均功率往往反而比连续波状态低，因为大部分能量转化为色散波辐射和离散梳齿了。

判断频梳质量时，我习惯先看对数频谱包络是否平滑。孤子微梳的包络在双对数坐标下接近双曲正割函数，也就是中心高、两侧单调下降，没有明显的凹陷。如果包络出现明显振荡或凹陷，可能是调制不稳定、色散波或高阶色散效应引起的。再看梳齿间距，理论上应该正好等于FSR，仿真中用离散网格表示时，梳齿阶数$k$对应频率偏移$k \times FSR$，所以频域图上相邻峰的横坐标间隔应该恒为1。

另外一个容易忽略的经验是：仿真中偶然一个随机种子会让系统落在不同分支上。也就是说，即使参数完全一样，换一个随机初值，最终可能停在连续波态、调制不稳定态或孤子态的不同位置。这不是代码错误，而是非线性系统的多稳态本性。想系统扫描参数找规律时，最好每个参数点多跑几个随机种子，做统计而非单次判断。

5. 常见问题与调试经验

5.1 仿真直接发散成NaN怎么处理

跑微环LLE仿真的头号问题是发散。现象就是跑着跑着psi出现NaN或Inf，或者时域功率突然暴涨到天文数字。排查思路按概率排序：

第一，时间步长太大。归一化LLE中的非线性和色散项在大场强下会产生很大的相位变化，dt必须足够小。保守起见，把dt从1e-3降到1e-4再跑一次，如果不再发散，就是步长问题。第二，泵浦强度过大导致初始非线性相位超过可处理范围。可以用较小的f（比如1左右）先验证程序能稳定运行，再逐步增大。第三，网格数$N$不足导致高频分量混叠。增大$N$到1024或2048试试，如果结果变化很大，说明网格不够。第四，初始场幅值设得过大。我习惯把初始噪声幅值设成$10^{-5}$量级，如果设成0.1甚至更大，一开始非线性就极强，非常容易发散。

如果问题出在失谐扫描过程中，注意扫描到某些中间失谐值时系统可能经历混沌瞬态，这个阶段场强峰值可能非常大，如果dt不够小，容易在混沌区发散。一个实用办法是给场强加一个瞬时限制：如果某个网格点的abs(psi)^2超过预设上限（比如1000），就把这一步的dt临时减半重算一遍。虽然这会让代码稍复杂，但能显著提升扫描的稳定性。

5.2 结果长期不收敛或“看起来不对”

有几种情况容易让人困惑。一种是跑了几百个时间单位，频域已经出现梳齿，但时域波形还在缓慢变化。这种情况要区分是物理混沌还是尚未收敛。判断方法很简单：把演化时间延长一倍，如果结果改变很大，说明还没收敛；如果时域形态和频谱包络基本不变，可以认为进入了稳态。

另一种是频域梳齿出现了，但梳齿间隔不是1，而是别的分数值。这通常意味着调制的空间周期不是$2\pi$的整数分之一，说明系统正在形成周期倍化或结构，而不是标准的FSR梳齿。这并不一定是错误，微环中确实存在周期倍化不稳定性。但如果想得到标准的等间隔频梳，需要检查参数是否落在标准孤子区域内。

还有一种是频域梳齿不对称，泵浦两侧功率差很大。这种不对称在实验中很常见，来自高阶色散或模式交叉。但归一化LLE里如果只保留二阶色散，理论上是左右对称的。如果仿真出现明显不对称，先检查代码有没有bug，比如频域波矢定义、fftshift使用是否正确。排查无误的话，可以在方程中添加高阶色散项模拟真实器件的模式交叉效应，这属于进阶玩法了。

5.3 参数调优经验速查表

把常用经验整理成一张速查表，方便对照调试。

这表只是起点，不同物理体系对参数的敏感程度差异很大。真正做课题时，建议做一次粗扫描找到大致区间，再在目标区间内细化参数网格。

5.4 我的几个调试小习惯

最后分享几个我自己实际用着顺手的习惯，供参考。

每次改参数前，我会先记一下当前参数下的腔内稳态平均功率和梳齿数量，作为“基线”。调参后对比这个基线，能快速判断参数变化是朝着哪个方向影响系统的。这比每次都重新看时域和频域图要高效得多。

跑失谐扫描时，保存最终psi而不是只保存功率迹线。很多时候扫完发现某个失谐处功率有奇特阶梯，但已经没保存对应时刻的场分布，只能重新跑一遍，浪费时间。

如果要做参数扫描研究趋势，千万不要把多组实验的顺序参数写死在循环里然后一次跑完。我在初期就犯过这个错误：三维参数扫描跑了整整一晚上，结果结果文件太大了，占了几十个GB磁盘，而实际有用的只有很小一部分。后来改成边算边存、只保存关键指标（稳态平均功率、梳齿数量、包络宽度），需要时再针对特定参数点单独重跑获取详细信息，效率高很多。

再补一个关于随机种子的小技巧。为了在参数变化时能区分“物理变化”和“随机种子引起的偶然差异”，我习惯固定随机种子（比如rng(42)），先在同一随机初值下做第一轮参数扫描，等判断出大致物理趋势后，再用几个不同随机种子在关键参数点做重复验证。这样做既能保证趋势的可比性，又能排除单次采样的随机性影响。

微环光频梳的LLE仿真，代码本身并不复杂，真正花时间的地方在参数理解、状态判断和调参迭代上。希望这份记录能帮刚开始接触这个方向的朋友少走一些弯路。按上面这套流程走一遍，应该能顺利跑出第一组梳齿，也对接下来的孤子态仿真有个明确的操作方向。