简介：面向非线性光纤光学与数值计算方向的学习者，这份资源提供用分步傅里叶法求解非线性薛定谔方程（NLS）的完整源代码解析，适合已具备Matlab基础、希望理解光脉冲在光纤中传输仿真实现的中高级读者。压缩包内共1个docx文档，约13KB，以文字讲解配合代码片段的方式组织，围绕NLS方程形式、分步傅里叶法的非线性步与色散步拆分、输入参数与模拟参数设置、sech与超高斯脉冲生成、色散相移因子构造以及时域频域结果可视化等模块逐段展开，并给出主循环中半步非线性、色散、再半步非线性的实现细节。已有1377人学习下载，读者可据此理清算法框架与代码对应关系，掌握步长、FFT点数等参数对仿真结果的影响，并可直接借鉴到自己的光纤传输数值实验中。

1. 一份能跑通的 NLS 方程数值解法：从光纤脉冲到分步傅里叶

做非线性光纤光学仿真的人，迟早会撞上非线性薛定谔方程（NLS）。解析解只在极少数理想条件下存在，真实场景里色散和非线性同时作用，光脉冲在传输过程中会展宽、压缩、产生自相位调制展宽频谱，甚至演化成孤子。这时候绕不开数值解法，而分步傅里叶法（Split-Step Fourier Method，SSFM）是工程上最常用、也最容易上手实现的一种。这份 MATLAB 源码做的事情很直接：给定光纤长度、色散系数、非线性系数和输入脉冲形状，用 SSFM 把 NLS 方程一步步推下去，最后画出输入输出的时域波形和频谱。适合正在做光纤通信、超连续谱产生、锁模激光器仿真，或者单纯想搞懂 SSFM 到底怎么落地的人。代码不长，但该有的都有——参数配置、网格生成、脉冲初始化、主循环、可视化，一条链路完整走通。

2. 代码结构拆解：从参数配置到主循环的完整链路

2.1 参数区：哪些量决定了仿真行为

拿到这份代码，第一眼应该看的是开头那段参数定义。它把整个仿真的自由度集中在了几个变量上，改这几个数就能覆盖大部分场景：

distance = 10; % 光纤长度（归一化单位） beta2 = -1; % 群速度色散系数，负值对应反常色散 N = 1; % 非线性参数，N=1 对应基阶孤子条件 mshape = 1; % 脉冲形状选择：0 为 sech，>0 为超高斯 chirp0 = 0; % 输入脉冲初始啁啾

beta2的符号很关键。负值意味着反常色散，和自相位调制（SPM）的相互作用可以形成孤子；正值是正常色散，脉冲会单调展宽。N的定义是 (N^2 = \gamma P_0 L_D)，当N=1且输入为基阶 sech 脉冲时，理论上应该看到孤子稳定传输——这也是验证代码正确性的一个经典测试用例。mshape控制初始脉冲包络，chirp0给脉冲加一个二次相位，模拟预啁啾或初始频率调制。

这些参数不是随便填的。我一般会先跑N=1, beta2=-1, mshape=0这组“教科书条件”，确认输出脉冲形状基本不变，再改其他参数做对比。如果一上来就调复杂参数，出了问题根本分不清是代码 bug 还是物理效应。

2.2 网格与频率轴：FFT 点数、窗口和步长的取舍

仿真参数区决定了数值精度和计算量：

nt = 1024; % FFT 点数 Tmax = 32; % 时间窗口半宽 step_num = round(20*distance*N^2); % z 方向步数 deltaz = distance/step_num; % z 步长 dtau = (2*Tmax)/nt; % 时间步长

nt和Tmax共同决定了时间分辨率dtau和频率分辨率。nt越大，频率轴越细，但每次 FFT 的开销也越大。Tmax要足够大，保证脉冲在传输过程中不会“撞墙”——也就是时域上跑出窗口边界，否则 FFT 的周期性假设会导致伪影。常见做法是让窗口至少覆盖脉冲宽度的十几倍。

step_num这里用了一个经验公式20*distance*N^2，意思是每单位长度大约走 20 步，再乘上非线性强度的平方做缩放。非线性越强，需要的步数越多，因为 SSFM 的误差主要来自把非线性算子和色散算子分开处理。步长太大，算子对易关系不成立，结果会偏。

频率轴的生成方式值得单独看一眼：

tau = (-nt/2:nt/2-1)*dtau; omega = (pi/Tmax)*[(0:nt/2-1) (-nt/2:-1)];

omega用的是 FFT 标准排列：前半段是正频率，后半段是负频率，没有做fftshift。这样后面ifft和fft出来的结果可以直接和omega对齐，省掉一次移位操作。如果自己改代码时把频率轴顺序搞乱了，色散相位就会加错方向，脉冲会往反方向展宽——这个坑我踩过不止一次。

2.3 输入脉冲生成：sech 与超高斯两条分支

脉冲初始化用了一个if-else分支：

if mshape==0 uu = sech(tau).*exp(-0.5i*chirp0*tau.^2); % 双曲正割脉冲 else uu = exp(-0.5*(1+i*chirp0).*tau.^(2*mshape)); % 超高斯脉冲 end

mshape=0走 sech 分支，这是孤子理论里的标准输入。mshape>0走超高斯分支，指数里的2*mshape控制边缘陡峭程度，mshape越大，脉冲越接近矩形。超高斯在模拟平坦脉冲或某些锁模激光器输出时有用。

注意 sech 分支的啁啾项是exp(-0.5i*chirp0*tau.^2)，超高斯分支把啁啾直接塞进了指数里(1+i*chirp0)。两种写法在chirp0=0时退化为纯实包络，没问题；但如果你要加比较大的啁啾，建议确认一下符号约定和你的物理定义是否一致。代码里chirp0默认是 0，先跑通再动它。

2.4 主循环：半非线性—色散—半非线性的对称分步

这是整份代码的核心。SSFM 的主循环采用对称分步格式（也叫 Strang 分裂），顺序是：半步非线性 → 全步色散 → 半步非线性。代码里把第一步非线性提到了循环外面：

dispersion = exp(i*0.5*beta2*omega.^2*deltaz); % 色散相位因子 hhz = i*N^2*deltaz; % 非线性相位系数 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); % 收尾半步非线性

逐段拆开看。dispersion是频域上的色散传播算子，exp(i*0.5*beta2*omega.^2*deltaz)对应 NLS 里色散项的解析解。hhz是非线性相位系数，abs(uu).^2是瞬时功率，乘上hhz就是 SPM 引入的相位。

循环体里，ifft(temp)把场变到频域，乘上dispersion完成一步色散传播，再fft回时域。然后temp = uu.*exp(abs(uu).^2.*hhz)做一整步非线性相位旋转。循环外的首尾两个半步拼起来正好是一个完整步，这种对称格式的局部误差是 (O(h^3))，比非对称的“先非线性后色散”精度高一阶。

有一个细节容易看漏：循环结束后uu = temp.*exp(-abs(uu).^2.*hhz/2)用的是-hhz/2，负号是因为前面已经多做了一个半步非线性，这里要“退回去”半个。如果你自己重写循环，把首尾半步搞反或者漏掉，结果在短距离内可能看不出大问题，但长距离传输会累积明显误差。

2.5 频谱计算与可视化：fftshift 和归一化因子

输出频谱的计算方式：

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

这里ifft之后做了fftshift，把零频移到中心，方便画图。乘上(nt*dtau)/sqrt(2*pi)是归一化，让频谱幅度和时域功率在量纲上对应。这个因子来自连续傅里叶变换和离散 FFT 之间的缩放关系，少了它频谱幅度没有物理意义，但形状不受影响。如果你只关心频谱展宽趋势，这个因子可以忽略；要做功率预算就得保留。

画图部分用了subplot(2,1,1)和subplot(2,1,2)分别画时域功率abs(uu).^2和频域功率abs(temp).^2，输入和输出各一张图。axis([-5 5 0 inf])把横轴限制在 ±5 的归一化时间范围内，纵轴从 0 开始。这个范围对 sech 脉冲够用，但如果你的脉冲展宽到窗口边缘，记得把axis范围调大，不然图上看不出展宽的全貌。

3. 跑起来：环境准备、参数修改与结果验证

3.1 最小运行路径

这份代码是纯 MATLAB 脚本，不依赖任何工具箱，基础 MATLAB 或兼容环境就能跑。步骤很简单：

1. 新建一个.m文件，把源码完整粘贴进去，保存为ssfm_nls.m。
2. 在命令行窗口cd到文件所在目录。
3. 直接输入文件名ssfm_nls回车，或者点编辑器里的运行按钮。

跑完会弹出两个图窗：Figure 1 是输入脉冲的时域和频谱，Figure 2 是输出脉冲的时域和频谱。如果 Figure 2 的时域波形和 Figure 1 几乎重合，说明孤子条件设置正确，代码工作正常。

提示：第一次跑之前把clear all保留着。脚本里用了clear all清空工作区，如果你在同一会话里反复调参数，不清空容易残留旧变量导致结果对不上。

3.2 用孤子条件做正确性自检

验证代码有没有写对，最省事的办法是跑基阶孤子：

distance = 10; beta2 = -1; N = 1; mshape = 0; % sech 脉冲 chirp0 = 0;

理论上，N=1的基阶孤子在反常色散光纤中传输时形状不变，只累积一个相位。跑完之后对比 Figure 1 和 Figure 2 的时域功率曲线，应该基本重合。如果输出脉冲明显展宽或压缩，先检查step_num是不是太小——把20改成50或100再跑一次，看结果是否收敛。收敛了说明是步长问题，不收敛就得回头查色散相位符号或者非线性系数。

再做一个对照实验：把N改成 0，也就是关掉非线性，只留色散。这时候脉冲应该单调展宽，展宽量大致符合色散长度 (L_D = T_0^2/|beta2|) 的预期。如果关掉非线性后脉冲反而压缩，那beta2的符号或者dispersion因子的写法肯定有问题。

3.3 改参数时先动哪些、后动哪些

参数不是平等的。我一般按这个顺序调：

nt和Tmax属于数值参数，不是物理参数。改它们之前先确认物理结果已经稳定，否则你分不清变化是物理效应还是数值误差。step_num的经验公式在N较大时可能不够，手动把系数从 20 提到 50 是常见操作。

3.4 输出图怎么读

Figure 1 和 Figure 2 各有上下两个子图。上图是归一化功率随归一化时间的变化，下图是频谱功率随归一化频率的变化。看几个典型特征：

• 时域波形变窄、频谱变宽：SPM 主导，非线性效应强。
• 时域波形展宽、频谱形状基本不变：色散主导。
• 时域出现多峰、频谱出现振荡结构：孤子裂变或者高阶效应开始显现。
• 频谱明显不对称：检查omega轴排列和fftshift是否用对。

如果 Figure 2 的曲线在窗口边缘突然翘起或者出现锯齿，大概率是脉冲跑出了时间窗口，FFT 周期性边界产生了混叠。把Tmax从 32 加到 64 再跑，问题通常消失。

4. 避坑与排查：SSFM 实现里最容易翻车的几个点

4.1 脉冲跑出时间窗口导致频谱出现假振荡

现象：输出频谱在高频端出现规则振荡或突然抬升，时域波形在窗口边缘不为零。

原因：Tmax不够大，脉冲展宽后尾部超出了[-Tmax, Tmax]范围。FFT 隐含周期性边界条件，跑出去的部分会从另一侧绕回来，形成混叠。

解决：把Tmax加倍，或者减小distance先看短距离行为。判断标准是输出时域波形在窗口两端应该衰减到接近零。如果加了Tmax后计算量太大，可以适当减小nt，但保证dtau仍能分辨脉冲特征。

4.2 步长过大导致能量不守恒或波形畸变

现象：传输过程中总能量sum(abs(uu).^2)*dtau明显漂移，或者孤子形状逐渐失真。

原因：SSFM 把非线性和色散分开处理，步长deltaz太大时两个算子不对易，分裂误差累积。非线性越强（N越大），对步长越敏感。

解决：把step_num公式里的系数从 20 提到 50 甚至 100，观察结果是否收敛。一个实用的自检习惯是每次跑完算一下输入输出能量比，偏差超过百分之几就加步数。对称分步格式比非对称格式宽容，但也不是无条件稳定。

4.3 频率轴顺序错位导致色散加反方向

现象：本该展宽的脉冲压缩了，或者孤子往错误方向演化。

原因：omega的排列必须和fft/ifft的输出顺序严格对应。代码里omega = (pi/Tmax)*[(0:nt/2-1) (-nt/2:-1)]是标准 FFT 顺序，如果你手动改成fftshift后的顺序，dispersion因子就会加错。

解决：不要单独改omega的生成方式。如果画图需要中心化频率轴，在画图时对频谱做fftshift，不要动omega本身。改完任何涉及 FFT 的代码后，用N=0的纯色散 case 验证脉冲展宽方向是否符合beta2符号。

4.4 首尾半步非线性写反或漏写

现象：短距离看不出问题，长距离传输后脉冲形状和理论预期有系统性偏差，且加步数也不收敛。

原因：对称分步要求“半步非线性 → 全步色散 → 半步非线性”，代码里第一个半步在循环外，最后一个半步用-hhz/2收尾。如果收尾符号写错，相当于多加了或漏减了半个非线性步。

解决：把主循环单独拎出来，用step_num=1跑一步，手动验算：输入经过半步非线性、一步色散、半步非线性后，结果应该和解析近似对得上。另一个办法是临时把非线性关掉（N=0），此时首尾半步都退化为恒等操作，如果结果和纯色散解析解一致，说明色散部分没问题，再打开非线性查半步。

4.5clear all掩盖了变量依赖

现象：单独跑脚本正常，但在已有工作区的会话里跑，结果不对或者报错。

原因：脚本开头clear all清空了所有变量，但如果你的运行环境或调用方式跳过了这一行（比如把代码嵌到函数里），残留变量可能覆盖参数。

解决：把脚本改造成函数形式，输入参数显式传入，避免依赖工作区状态。函数化之后也方便批量扫参数——用循环调函数比反复改脚本再运行干净得多。

5. 二次开发：把脚本改造成可复用的仿真函数

这份代码作为教学和快速验证够用，但要拿来做批量仿真或者嵌入更大的项目，直接改脚本会很别扭。我的做法是把它重构成一个函数，输入物理参数、输出时域和频域结果，画图交给调用方决定。

function [tau, u_out, omega, U_out] = ssfm_nls_solve(distance, beta2, N, mshape, chirp0, nt, Tmax, steps_per_unit) % SSFM_NLS_SOLVE 用对称分步傅里叶法求解 NLS 方程 % 输入： % distance - 传输距离 % beta2 - 群速度色散系数 % N - 非线性参数 % mshape - 0 为 sech，>0 为超高斯阶数 % chirp0 - 初始啁啾 % nt - FFT 点数 % Tmax - 时间窗口半宽 % steps_per_unit- 每单位长度的步数（替代硬编码的 20） % 输出： % tau, u_out - 时间轴与输出场 % omega, U_out - 频率轴与输出频谱 dtau = (2*Tmax)/nt; tau = (-nt/2:nt/2-1)*dtau; omega = (pi/Tmax)*[(0:nt/2-1) (-nt/2:-1)]; if mshape == 0 u = sech(tau).*exp(-0.5i*chirp0*tau.^2); else u = exp(-0.5*(1+i*chirp0).*tau.^(2*mshape)); end step_num = max(1, round(steps_per_unit*distance*max(N^2,1))); deltaz = distance/step_num; dispersion = exp(i*0.5*beta2*omega.^2*deltaz); hhz = i*N^2*deltaz; u = u.*exp(abs(u).^2.*hhz/2); for k = 1:step_num u = fft(ifft(u).*dispersion); u = u.*exp(abs(u).^2.*hhz); end u = u.*exp(-abs(u).^2.*hhz/2); u_out = u; U_out = fftshift(ifft(u)).*(nt*dtau)/sqrt(2*pi); end

改造的关键点有三个。第一，把steps_per_unit暴露成参数，不再硬编码20，这样调用方可以根据N的大小自己控制精度。第二，step_num里加了max(N^2,1)，避免N=0时步数变成零。第三，去掉了所有figure和plot，函数只负责算，画图由外部做——这样批量扫参数时不会弹出一堆图窗。

函数化之后，扫参数就变成几行循环：

N_list = [0.5, 1, 1.5, 2]; for idx = 1:length(N_list) [tau, u, omega, U] = ssfm_nls_solve(10, -1, N_list(idx), 0, 0, 1024, 32, 50); P(:,idx) = abs(u).^2; end plot(tau, P); legend('N=0.5','N=1','N=1.5','N=2');

验证改造是否等价，最直接的办法是拿原始脚本跑一组参数，再用函数跑同一组，对比abs(uu).^2的最大绝对误差。只要在浮点精度范围内一致，改造就没引入行为变化。我一般会固定N=1, mshape=0这组孤子条件做回归，因为它的理论预期最明确，改坏了立刻能看出来。

再往深走一步，可以加一个能量监控：每步记录sum(abs(u).^2)*dtau，如果漂移超过阈值就自动减小步长重跑。这个自适应步长的思路在非线性很强的时候特别有用，比手动试步数省事。不过要注意，自适应逻辑本身会改变计算路径，验证的时候得和固定步长版本对得上。

从那以后我每次改这类分步傅里叶代码，都强制先跑一遍N=1孤子回归和N=0纯色散对照，两个都过了才敢动物理参数。希望帮到你。

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