企业官网建设流程全解析

简介：这份DOCX文档面向非线性光纤光学与数值仿真初学者，围绕分步傅里叶法求解非线性薛定谔方程（NLS）展开，帮助读者理解色散与非线性效应共同作用下的光脉冲传输过程。压缩包仅含1个DOCX文件，约13KB，篇幅精炼，核心是一份带注释的Matlab程序及配套解析，涵盖光纤长度、色散与非线性参数、FFT点数、时间窗口和z步长设置，sech与超高斯输入脉冲生成，以及1/2N→D→1/2N对称分步流程。文档对非线性相移因子、色散相移因子和频域变换等关键语句逐段说明，并给出输入与输出脉冲的时域、频域可视化结果，便于读者对照公式理解代码实现。已有1376人学习，适合希望快速掌握NLS数值求解思路、搭建光纤传输仿真脚本或排查分步傅里叶实现细节的读者参考。

1. 把光纤切成两百段，一段一段推着走

调孤子传输的人多半撞过同一个场面：把色散参数从 -1 改成 -1.5，图上的脉冲既没有分裂也没有偏移，而是原地炸成一团噪声。回头翻手上那份六十行的 MATLAB 脚本，问题往往不在物理，而在分步傅里叶法的步长和窗口。这份代码干的就是这件事——给定色散、非线性和输入脉冲，用分步傅里叶法（Split-Step Fourier Method，SSFM）把非线性薛定谔方程从 z=0 推到 z=distance，再把时域波形和频谱画出来。它适合两类人：一类做非线性光纤光学、孤子、超连续谱仿真，需要一个能跑能改的起点；另一类想把分步傅里叶这个数值方法真正弄懂，而不是停在公式上。整份脚本不到六十行，不做函数封装，不依赖任何工具箱，纯基础语法就能跑起来。

2. 算法骨架：NLS 方程在代码里被拆成哪两块

2.1 从方程到算子分裂：一个在频域乘，一个在时域乘

脚本开头那行注释已经把方程写死了：

% idu/dz - sgn(beta2)/2 d^2u/d(tau)^2 + N^2*|u|^2*u = 0

把它移项，得到的是u沿z的演化率：等号左边是i·∂u/∂z，右边是色散项(sgn(β2)/2)·∂²u/∂τ²减去非线性项N²|u|²u，整体再除以i。这件事的关键在于，右边这两项各有各的“舒适区”：色散项对τ求二阶导，在频域里就是一串乘法；非线性项是|u|²乘u，在时域里是逐点相乘。一个天然属于频域，一个天然属于时域，中间只隔一次 FFT。

把这两项分别看成两个算子，各自在一步h内的解析解都非常干净：

% 色散算子 D：频域乘一个相位核，h 是这一步的步长 % u_tau_tau 在频域变成 -(omega^2) * u_hat D_kernel = exp(1i * 0.5 * beta2 * omega.^2 * h); % 非线性算子 N：时域乘一个纯相位，模长不变 % 方程里是 i*N^2*|u|^2*u，所以指数上带 i N_kernel = exp(1i * N^2 * abs(u).^2 * h);

这里有个容易被忽略的点：非线性算子的指数上只有相位，没有幅度，所以|u|²在这一步内是常数——先用当时的功率算出相位，再整体乘上去，功率本身不变。这正是分步傅里叶法能在时域直接处理非线性的原因，也是它比直接打靶法快几个数量级的地方。参数h就是脚本里的deltaz，omega是按 FFT 顺序排好的角频率数组，beta2在标准归一化下取 ±1。

误差从哪来？两个算子不对易。∂²/∂τ²和|u|²乘在一起，先算哪个结果不一样，所以一次“先色散后非线性”的拼接必然有残差。这个残差的大小，决定了下面要不要用对称分步。

2.2 为什么要对称：1/2N → D → 1/2N 而不是 N → D

最朴素的分步是“整步非线性 + 整步色散”，局部误差 O(h²)，全局一阶。脚本用的是对称形式：先走半步非线性，再走整步色散，最后补半步非线性，也就是常说的 Strang 分裂。局部误差降到 O(h³)，全局二阶，代价只是多算一次非线性相位，几乎白送。

% 对称分步：半步 N -> 整步 D -> 半步 N temp = uu .* exp(abs(uu).^2 .* hhz / 2); % 前半步非线性 for n = 1:step_num f_temp = ifft(temp) .* dispersion; % 整步色散，频域乘法 uu = fft(f_temp); % 回到时域 temp = uu .* exp(abs(uu).^2 .* hhz); % 整步非线性 end uu = temp .* exp(-abs(uu).^2 .* hhz / 2); % 补后半步非线性

注意hhz = 1i * N^2 * deltaz是提前算好的，循环里只做乘法和指数。前后两个半步合起来等于一个整步，位置卡在色散步的两端，这样整条链路上非线性算子对色散算子是对称包围的，误差主项相互抵消。这也是为什么改这份代码时，不要随手把首尾两个半步删掉——删掉就退回一阶精度，同样的step_num下结果会明显变糙。

理论上还可以用更高阶的分裂（比如四阶的 Yoshida 系数），但在这类问题里，步长误差通常不是主导，窗口和步数选错带来的误差要大得多，所以二阶足够。

2.3 五段式骨架：改哪一段会动到什么

整份脚本没有函数，从上到下一条直线，按职责可以切成五段。看懂这张表，后面动刀的时候就知道自己在动哪根筋。

第一段是唯一“改数字不会心虚”的地方；第二段每改一个数都要重新验算dtau和omega的间隔；第三段最容易埋雷，尤其是omega里那段[(0:nt/2-1) (-nt/2:-1)]的拼接顺序；第四段基本不动，除非要把固定步长换成自适应；第五段可以随便重写，只要把fftshift和(nt*dtau)/sqrt(2*pi)这个因子保留。

3. 参数与网格：nt、Tmax、step_num 这三个数定了仿真成败

3.1 五个物理参数：单位全藏在归一化里

脚本里的物理量都是归一化的，没有米、没有秒、没有瓦。distance的单位是色散长度 L_D，tau的单位是脉宽 τ0，N的定义是 sqrt(L_D/L_NL)，也就是孤子阶数。基孤子就是 N=1。

beta2这个名字有点误导。方程注释里写的是sgn(beta2)，但代码里直接用了beta2本身：dispersion = exp(1i*0.5*beta2*omega.^2*deltaz)。当beta2 = -1时，这个核等于exp(-1i*0.5*ω²h)，和sgn(β2)/2 = -1/2的推导完全对得上，所以默认值是自洽的。但如果你把beta2写成 -2，色散相位会翻倍，而非线性项系数没跟着变，实际上等于换了一组归一化，物理参数已经不是原来那组了。这个点后面还会再提。

mshape的分支在初值那一段：

if mshape == 0 uu = sech(tau) .* exp(-0.5i*chirp0*tau.^2); % 双曲正割，基孤子解 else uu = exp(-0.5*(1+1i*chirp0).*tau.^(2*mshape)); % 超高斯，阶数为 mshape end

mshape = 0时生成的是标准孤子形状，峰值功率为 1，和N = 1搭配正好是基孤子。mshape > 0时生成超高斯，tau.^(2*mshape)意味着阶数越高顶部越平、边沿越陡。注意这条分支里chirp0是直接乘在指数上的，而不是像 sech 那条那样乘tau.^2，所以同样一个chirp0，在两条分支下的实际啁啾量纲不一样，跨分支对比时要留意。

选sech还是超高斯，取决于你想看什么物理。要看孤子自频移、孤子分裂、色散波辐射，用sech；要看超连续谱的平坦度、看陡边沿激发的色散波，用超高斯。这两条分支的谱宽差很多，同一个Tmax未必都装得下，这是后面踩坑的一个来源。

3.2 时间轴与频率轴：omega 的顺序决定了 FFT 能不能对上

这一段是整份代码最需要看懂的地方：

dtau = (2*Tmax)/nt; % 时间步长 tau = (-nt/2:nt/2-1)*dtau; % 时间轴，从 -Tmax 到 Tmax-dtau omega = (pi/Tmax)*[(0:nt/2-1) (-nt/2:-1)]; % 频率轴，FFT 顺序，未 fftshift

nt = 1024、Tmax = 32，算出来dtau = 64/1024 = 0.0625，tau从 -32 走到 31.9375，总共 64 个时间单位。这里的Tmax是半窗宽，真正的窗口长度是2*Tmax，这点在改参数时很容易记反。

omega的间隔是pi/Tmax，为什么不是2*pi/(2*Tmax)？因为2*pi/(2*Tmax)就是pi/Tmax，两者是一回事。这个间隔来自 FFT 的频率分辨率：采样间隔dtau对应 Nyquist 角频率pi/dtau = pi*nt/(2*Tmax)，一共nt个点，所以相邻间隔就是pi/Tmax。用nt/2个正频率加上nt/2个负频率，拼出来的顺序正好是 MATLABfft输出的顺序：第一个元素是直流，接着是正频率递增，后半段是负频率。

这个顺序非常重要。因为omega是按 FFT 顺序排的，所以主循环里可以直接写ifft(temp) .* dispersion，不需要fftshift把零频搬到中间；一旦你为了画图方便把omega改成fftshift过的顺序，主循环里的乘法就会把高频成分乘到低频上，结果看起来像频谱被“搅拌”过。只有出图那一句才需要fftshift：

temp = fftshift(ifft(uu)) .* (nt*dtau) / sqrt(2*pi); plot(fftshift(omega)/(2*pi), abs(temp).^2, 'b', 'linewidth', 2);

fftshift(omega)把零频挪到中间，除以2*pi换成普通频率，横轴就落在 ±8 的范围内。脚本里画图时用axis([-.5 .5 0 inf])只截了中间一小段，因为归一化之后脉冲的谱宽通常远小于 Nyquist，全画出来两边全是零，看不出结构。

还有一点：tau的起点是-nt/2，这在 MATLAB 里是负索引，(-512:511)一共 1024 个点，没问题。但如果有人图省事写成(-nt/2+1:nt/2)，点数还是 1024，可时间轴整体平移了半个步长，和omega的相位关系就对不上了，做频域分析时会多出一个线性相位。

3.3 步数经验公式：20distanceN^2 从哪来，什么时候失效

step_num = round(20*distance*N^2); deltaz = distance/step_num;

默认参数下step_num = round(20*10*1) = 200，deltaz = 0.05。这个 20 不是随便拍的：它约束的是单步内的最大非线性相移。基孤子峰值功率为 1，N = 1，单步非线性相位大约是N²·|u|²·deltaz ≈ 0.05弧度，远小于 1，对称分步的误差主项就很小。经验上单步非线性相移控制在 0.1 弧度以内比较稳，20 这个系数大概留了两倍余量。

这个公式在几种情况下会失灵。第一种是脉冲峰值功率远大于 1，比如你把初值改成5*sech(tau)，功率变成 25 倍，单步相位直接到 1.25 弧度，必须把系数从 20 提到 200 以上。第二种是加了损耗或者增益，功率沿 z 变化，固定步长在功率高的前段精度不够、后段又浪费。第三种是distance很小的时候，round会把它取整到 0，deltaz = distance/0 = Inf，下一句exp(1i*...*Inf)直接出 NaN，整张图空白，新手最容易在这里卡住。稳妥的写法是加一句step_num = max(step_num, 1)，或者max(step_num, 100)。

4. 主循环拆解：1/2N → D → 1/2N 对称分步的实现细节

4.1 预计算相移因子：把 exp 提到循环外

主循环之前有两行看起来不起眼、实际上是性能关键的代码：

dispersion = exp(1i*0.5*beta2*omega.^2*deltaz); % 色散相位核，长度 nt hhz = 1i*N^2*deltaz; % 非线性系数，标量

dispersion是一个长度 1024 的复数向量，只和omega、beta2、deltaz有关，和当前场u无关，所以整场仿真只需要算一次。如果把它写进循环，200 步就是 200 次长度 1024 的exp，在 MATLAB 里虽然向量化之后也不算慢，但完全没必要。这是一个很典型的“把不变量提出循环”的写法。

hhz更简单，就是个标量，后面所有非线性步都用它乘|u|²。注意它是纯虚数，所以exp(|u|²·hhz)出来的一定是单位模的复数，只改相位不改幅度。这一点可以用来做自检：任何一步之后，abs(uu)的能量应该和上一步几乎一样（色散步也不改总能量，只是重新分布），如果能量明显漂移，说明步长或者窗口出问题了。

4.2 主循环逐行：ifft 在前还是 fft 在前

for n = 1:step_num f_temp = ifft(temp) .* dispersion; % 时域 -> 频域，乘色散核 uu = fft(f_temp); % 频域 -> 时域，回到时间轴 temp = uu .* exp(abs(uu).^2 .* hhz); % 整步非线性 end

为什么是ifft在前、fft在后，而不是反过来？因为 MATLAB 的ifft是带1/N归一化的正变换（指数上是正号），fft是负号。这里用的傅里叶约定是U(ω) = ∫u(τ)e^{+iωτ}dτ，和方程里i·∂u/∂z的符号配套：色散项在频域变成-ω²，再乘上方程前面的-i·sgn(β2)/2和这里的i，最后落到exp(i·0.5·β2·ω²·h)这个核上。如果换成fft在前，符号就反了，等价于把beta2取了相反数——图上看就是脉冲往反方向色散，异常色散变正常色散。

ifft和fft互为逆变换，MATLAB 里ifft(fft(u))严格等于u，所以这一来一回不会引入额外的幅度因子，dispersion的模恒为 1，能量守恒靠的就是这一点。这也是为什么主循环里不需要任何归一化系数。

循环里temp和uu交替使用，读的时候容易晕。捋一下：进入循环时temp是“已经走过半步非线性的场”，ifft(temp)把它送进频域，乘dispersion得到f_temp，再fft回时域存进uu，这时uu是“走完一个整步色散的场”，然后temp = uu .* exp(...)得到“又走了一个整步非线性的场”，进入下一轮。所以循环体的语义是“D 然后 N”，而首尾两个半步 N 把它补成了对称形式。

4.3 收尾半步与那行有点绕的赋值

uu = temp .* exp(-abs(uu).^2 .* hhz/2); % 补上后半步非线性

这行里出现了两个变量：temp是循环最后一次赋值的“走完整步非线性”的场，uu是它之前那个“走完整步色散”的场。指数里用的是abs(uu).^2，而不是abs(temp).^2。两者模长其实一样，因为非线性步只改相位，所以数值上没问题。但从可读性讲，这里写abs(temp).^2更直白，写uu会让读代码的人停下来想一秒。

更值得说的是这个收尾半步本身。前面算过：整个流程的非线性步数是0.5 + step_num + 0.5 = step_num + 1，而色散步数只有step_num。严格对称的写法应该是step_num个色散步配step_num个非线性步，也就是首尾各半步、中间step_num - 1个整步。现在多出的这半步非线性，量级是0.05/2弧度，对结果影响很小，但在distance很大或者N很大的时候会累积成一个可观测的相位偏差。后面避坑那一节会给出改法。

4.4 出图与频谱归一化：(ntdtau)/sqrt(2pi) 是怎么来的

temp = fftshift(ifft(uu)) .* (nt*dtau) / sqrt(2*pi);

这一句把时域场转成频谱，三个动作叠在一起。ifft(uu)给出的是离散傅里叶系数，MATLAB 的定义里带1/N；乘上nt就还原成不带归一化的求和Σ u_n e^{iωτ_n}；再乘dtau，这个求和就近似成了积分∫u(τ)e^{iωτ}dτ，也就是连续傅里叶变换。最后除sqrt(2*pi)是一个约定俗成的归一化，让正反变换的系数对称。三步合起来正好是(nt*dtau)/sqrt(2*pi)，和代码里写的一模一样。

如果只是想看谱的形状和相对宽度，这个因子其实可以省掉，因为它是常数，不影响峰位和宽度。但如果要比较不同参数下的绝对谱功率，或者要算频谱积分验证 Parseval 定理，这个因子就必须留着，否则量级会差好几个数量级。

出图部分用subplot(2,1,1)和subplot(2,1,2)分别画时域功率和频域功率，axis([-5 5 0 inf])把时间轴截在 ±5，因为基孤子的能量基本都在这个范围内，画满 ±32 什么都看不清。这个截断范围要跟着脉宽走：如果换成超高斯、或者把chirp0调大导致脉冲展宽，±5 就不够用了，得改成 ±10 甚至更大。

5. 避坑与排查：五类跑崩现场的现象、原因和解法

5.1 图是空的、变量全是 NaN

现象：脚本跑完没有任何报错，但figure(2)里一片空白，工作区里uu全是NaN。

原因：step_num被round取整成 0 了。step_num = round(20*distance*N^2)，当distance小于 0.025 且N = 1时，这个值就落到 0；紧接着deltaz = distance/step_num得到Inf，dispersion = exp(1i*...*Inf)直接变NaN，后面一路污染。

解决：把步数定死一个下限，step_num = max(round(20*distance*N^2), 100)。短距离仿真本来就该用固定步数，靠经验公式反而不稳。改完顺手在deltaz后面加一句assert(isfinite(deltaz))，下次再出这种事会当场停住，而不是跑完了才发现图是空的。

5.2 脉冲峰值炸到几百，或者直接发散

现象：传输到一半，时域峰值从 1 冲到几十甚至上百，再往后变成NaN，频谱铺满整个窗口。

原因：单步非线性相位太大。deltaz由20*distance*N^2反推，只对峰值功率约为 1 的输入成立。一旦你把初值改成3*sech(tau)，功率变 9 倍，或者把N提到 3，单步相位就从 0.05 涨到 0.45 甚至更大，对称分步的误差主项压不住，数值解会往高频堆能量，最后爆掉。

解决：按峰值功率重新估步数。经验做法是先算P0 = max(abs(uu).^2)，再取step_num = max(round(200*distance*N^2*P0), 500)。如果还是不稳，就得上自适应步长：每步用局部误差估计（比如比较一个整步和两个半步的结果）来决定下一步的deltaz。这份脚本是固定步长，改造的时候把主循环里的deltaz换成变量即可，dispersion和hhz都要跟着每步重算，性能会掉一些，但换来的稳定性是值得的。

5.3 频谱两端出现规则的锯齿或者镜像

现象：频谱图在靠近 Nyquist 的边缘出现不该有的凸起，形状像是从另一边折过来的，改step_num也不消失。

原因：时间窗口不够。脉冲在传输过程中会展宽（正常色散下尤其明显），或者孤子自频移会把能量推到时间轴的一端，一旦越过±Tmax的边界，FFT 的周期性假设就把它从另一端接了回来，形成循环卷积混叠。这时候加步数没用，因为问题出在窗口不在步长。

解决：加大Tmax，同时保持dtau不变，也就是按比例加大nt。Tmax翻倍、nt翻倍，dtau才不变，频率分辨率也才不变。只改Tmax不改nt，dtau会变大，Nyquist 频率降低，反而更容易混叠。判断窗口够不够的土办法是看传输结束时的时域图，如果|u|²在边界处还没降到 1e-6 以下，窗口就是不够。

5.4 改了 beta2 的数值，结果和预期完全对不上

现象：把beta2从 -1 改成 -2，以为只是色散变强，结果脉冲形状和展宽速度都不是两倍关系。

原因：这份脚本里beta2既当符号又当大小。方程注释写的是sgn(beta2)，但色散核里直接用了beta2本身。取 -1 时sgn和原值恰好一致，所以默认参数下没问题；取 -2 时色散相位翻倍，而非线性项N²没动，等于偷偷改了L_D和L_NL的比值，物理问题已经换成另一个了。

解决：如果只是想要正常色散，用beta2 = +1；如果确实要改色散大小，应该同时按新色散长度重新归一化distance和N，而不是单独改一个数。更规范的做法是把色散核改成exp(1i*0.5*sign(beta2)*omega.^2*deltaz)，让beta2只贡献符号，大小交给归一化吸收。改完之后记得把N的定义也核一遍，N² = L_D/L_NL里的L_D是跟着|beta2|走的。

5.5 能量不守恒，或者对称性差半步

现象：跑完之后总算能量sum(abs(uu).^2)*dtau和初值差了千分之几，或者把step_num加密一倍，结果收敛得很慢。

原因：两个。一是前面说的首尾半步问题，非线性步比色散步多了整整一步，破坏了严格对称；二是最后一行用了abs(uu).^2而不是abs(temp).^2，虽然数值等价，但掩盖了“这一步应该基于哪一时刻的功率”这个语义。

解决：把循环改成for n = 1:step_num-1，让中间的整步非线性只走step_num-1次，首尾各半步，总数正好对齐step_num个色散步：

temp = uu .* exp(abs(uu).^2 .* hhz/2); for n = 1:step_num-1 uu = fft(ifft(temp) .* dispersion); temp = uu .* exp(abs(uu).^2 .* hhz); end uu = fft(ifft(temp) .* dispersion); uu = uu .* exp(abs(uu).^2 .* hhz/2);

改完之后再算能量，基孤子传输 10 个色散长度，相对误差应该能压到 1e-6 量级。如果还是差很多，那就不是对称性的问题，得回去查窗口和步长。

6. 二次开发：把脚本改成能扫参、能自检的工具函数

固定参数的脚本跑一次很爽，跑一百次就是灾难。我一般第一步先把整段逻辑包成函数，让参数从外面进来，把tau、omega、uu一起返回：

function [uu, tau, omega, E] = ssfm_nls(distance, beta2, N, mshape, chirp0, nt, Tmax) % 输入：物理参数 + 网格参数；输出：末态场、时间轴、频率轴、总能量 dtau = 2*Tmax/nt; tau = (-nt/2:nt/2-1)*dtau; omega = (pi/Tmax)*[(0:nt/2-1) (-nt/2:-1)]; step_num = max(round(200*distance*N^2), 500); deltaz = distance/step_num; if mshape == 0 uu = sech(tau).*exp(-0.5i*chirp0*tau.^2); else uu = exp(-0.5*(1+1i*chirp0).*tau.^(2*mshape)); end E0 = sum(abs(uu).^2)*dtau; % 入射能量，用于守恒校验 dispersion = exp(1i*0.5*beta2*omega.^2*deltaz); hhz = 1i*N^2*deltaz; temp = uu .* exp(abs(uu).^2 .* hhz/2); for n = 1:step_num-1 uu = fft(ifft(temp) .* dispersion); temp = uu .* exp(abs(uu).^2 .* hhz); end uu = fft(ifft(temp) .* dispersion); uu = uu .* exp(abs(uu).^2 .* hhz/2); E = sum(abs(uu).^2)*dtau; assert(abs(E-E0)/E0 < 1e-6, '能量漂移 %.2e，检查步长与窗口', abs(E-E0)/E0); end

这段改造里有三个点值得说。第一，步数系数从 20 提到 200、下限设到 500，是拿稳定性换时间，单次跑慢一点，但扫参的时候不会中途炸掉。第二，加了能量断言，这是最省事的自检手段：只要assert没触发，就说明窗口和步长至少没坏到离谱，比肉眼判断图对不对靠谱得多。第三，把tau和omega一起返回，画图和后续分析就不用在外面重算一遍，避免网格对不齐。

有了函数，扫参就是一层循环：

验证方法也很直接：把distance设成 1、N设成 1、beta2设成 -1、mshape设成 0，跑完看末态波形和初值的差。基孤子在L_D量级的距离上形状几乎不变，如果这时候就差很多，说明步长或窗口有问题，别急着往上加距离。等这个自检过了，再跑distance = 10看周期性呼吸，心里就有底了。这套流程我吃过亏：早些年直接上长距离大N，图出来花里胡哨还挺好看，后来才发现是数值噪声被非线性放大出来的假结构，白高兴一场。

从那以后，我每次动beta2或者改初值之前，都会先跑一遍distance = 1的基孤子自检，再顺手看一眼能量断言有没有触发，两步都过了才敢往下做。希望帮到你。

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