分步傅里叶法解NLS方程：从源码到孤子模拟全解析

简介：分步傅里叶法解非线性薛定谔方程的源代码，面向光纤通信、非线性光学方向的研究生与工程师，用于模拟光脉冲在光纤中的传输演化。包内共1个docx文档，压缩包仅13KB，文档内嵌完整Matlab源代码，包含输入参数设置、FFT网格划分、色散相移因子计算、非线性/色散交替迭代以及输入输出脉冲时域频域绘图等模块。代码结构紧凑，核心思路为将传输过程拆分为非线性步与线性色散步，通过傅里叶变换快速求解，并给出了sech脉冲与超高斯脉冲两种初始条件下的运行示例。已有1367人学习下载，适合正在学习分步傅里叶法或需要快速搭建NLS仿真程序的读者。结合代码解析，可帮助理解β2色散参数、非线性系数N、啁啾chirp等因素对脉冲演化与频谱展宽的影响，省去从零编写和调试的时间。

1. 分步傅里叶法解NLS方程：先弄清楚这段源代码在算什么

非线性光纤里，光脉冲的传输几乎都可以用非线性薛定谔方程（NLS）来描述。真正动手模拟之前，很多人第一反应是直接做分步傅里叶法，然后栽在符号约定、FFT方向、步长选取这些细节上。这段源代码给了一个可以直接跑的参考实现：用对称分步傅里叶法，把一个双曲正割或超高斯脉冲沿光纤推10个色散长度，并画出输入输出脉冲的时域波形和频谱。适合刚接触光纤非线性模拟的研究生、光通信链路工程师，以及想验证自己另写的 C++/Python 求解器的人。我建议先跑通默认参数，再动手改 chirp、孤子阶数 N，观察基阶孤子传输不变形、N=2 周期性呼吸这些经典现象。理解了这一段，split-step 的骨架基本就清楚了。

2. 从归一化NLS到split-step实现：源码结构与参数怎么对应

2.1 归一化NLS方程与符号约定

先看代码开头注释里的方程：

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

更规范的写法是：

i∂u/∂z - (sgn(β2)/2) ∂²u/∂τ² + N²|u|²u = 0

这里的 u 是归一化慢变包络，τ 是以输入脉冲宽度 T0 归一化的时间，z 是以色散长度 L_D 归一化的传输距离。β2 是群速度色散参数，N 是孤子阶数，且满足 N² = L_D / L_NL。当 β2 < 0 即反常色散区时，方程才可能支持亮孤子；β2 > 0 的正常色散区则对应暗孤子和正常色散展宽。代码默认 distance=10，β2=-1，N=1，mshape=1。这里要特别注意：方程里写的是 sgn(β2)/2，但代码的 dispersion 因子里直接用 beta2，相当于把 β2 的符号和大小都放进去了。做参数扫描时如果只把 beta2 改成 -5，色散项强度也会跟着变，距离轴含义会不一样，不要把它简单当成“换个符号”。

下面这张表把代码开头的输入参数和模拟参数一次说清。

2.2 网格、步长与频率轴

模拟网格是这段代码里最容易被抄错的点。时间步长取 dtau = 2*Tmax/nt，所以时间数组覆盖 [-Tmax, Tmax) 共 nt 个点。频率轴没有用 linspace，而是：

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

这行构造的 omega 是角频率，单位 1/τ，范围大约是 [-π/dtau, π/dtau)。关键是排列顺序：先正频率后负频率，和 Matlab 的 fft/ifft 输出顺序一致。这样后面 temp=fftshift(ifft(uu)) 画频谱时才能用 fftshift 挪到对称区间。如果自己重写代码时用[-nt/2:nt/2-1]*domega这类写法，必须同时调好 fftshift 的位置，否则色散相移exp(i*0.5*beta2*omega.^2*deltaz)作用在错误的频率分量上，脉冲会直接散掉。

纵向步长由下面两行决定：

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

默认 distance=10，N=1 时 step_num=200，deltaz=0.05。N=2 时同样是 distance=10，step_num=800，步长缩小到 0.0125，因为非线性强了以后对步长更敏感。这个“20”是经验系数，不是严格精度上限；后面讲精度检查时会说怎么判断它够不够。

2.3 可直接运行的完整源代码

我保留了原代码的逻辑，只补了少量注释，方便一行行对照。

% 分步傅里叶法解归一化NLS方程 % idu/dz - sgn(beta2)/2 d^2u/d(tau)^2 + N^2*|u|^2*u = 0 % 输入参数 distance = 10; % 光纤长度 beta2 = -1; % 色散系数，带符号 N = 1; % 非线性参数（孤子阶数） mshape = 1; % 0: sech脉冲; >0: 超高斯脉冲 chirp0 = 0; % 输入啁啾 % 模拟参数 nt = 1024; Tmax = 32; step_num = round(20*distance*N^2); deltaz = distance / step_num; dtau = (2*Tmax) / nt; % 时间与频率网格 tau = (-nt/2:nt/2-1) * dtau; omega = (pi/Tmax) * [(0:nt/2-1) (-nt/2:-1)]; % 输入脉冲 if mshape == 0 uu = sech(tau) .* exp(-0.5i * chirp0 * tau.^2); else uu = exp(-0.5 * (1 + 1i*chirp0) .* tau.^(2*mshape)); end % 初始频谱（仅用于显示） temp0 = fftshift(ifft(uu)) .* (nt*dtau) / sqrt(2*pi); % 预计算色散相移与非线性相移 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 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); % 输出频谱（用于显示） tempf = fftshift(ifft(uu)) .* (nt*dtau) / sqrt(2*pi);

代码整体分四步：构造网格、生成初始场、预计算两个相移因子、进入纵向循环。dispersion 和 hhz 放在循环外，是因为它们每个步长都相同；把复数指数提前算好，能避免循环里重复算 exp。deltaz 缩小到原来一半时，dispersion 里的 omega²*deltaz 要一起变，不能只改 step_num。这也是我把 distance 和 step_num 分开写的原因。

temp = uu .* exp(abs(uu).^2 .* hhz/2);中 hhz/2 对应半非线性步，主循环内 hhz 对应全非线性步，最后的-hhz/2用于抵消最后一次循环里多算的半步，让每个 z 步长都等价为“半非线性-色散-半非线性”。这个在第 3 章细说。

3. 对称分步傅里叶的细节：色散相移、非线性相移与FFT约定

3.1 为什么色散在频域乘相因子

线性部分把方程里的非线性项拿掉，得到：

i∂u/∂z - (sgn(β2)/2) ∂²u/∂τ² = 0

做傅里叶变换后，∂²/∂τ² 对应 -ω²，于是：

i ∂u~/∂z + (sgn(β2)/2) ω² u~ = 0

即 ∂u~/∂z = i (sgn(β2)/2) ω² u~，所以一个步长 dz 的解析解是在频域乘以：

exp(-i (sgn(β2)/2) ω² dz)

代码里 beta2=-1 时 dispersion=exp(i0.5beta2omega²deltaz)=exp(-i0.5omega²*deltaz)，完全对得上。也就是说，色散项的数值误差来源只有一个：omega 网格和 dtau 是否匹配。如果 dtau 改变而 omega 没有按2π/(nt*dtau)更新，那么这个相位因子就错了，而且错得很有隐蔽性——波形看起来还是光滑的，只是宽度、速度不对。

常见错误是把omega = (2*pi/(nt*dtau))*[0:nt/2-1 -nt/2:-1]写成了omega = 2*pi*[0:nt/2-1 -nt/2:-1]/(nt*dtau)忘记括号，导致低频分量的因子全错。建议先打印omega(2)-omega(1)，确认它等于2*pi/(nt*dtau)。这段代码用pi/Tmax作为频率分辨率，因为nt*dtau = 2*Tmax，两种写法等价。

3.2 半步非线性技巧与hhz的由来

非线性步的处理基于同一 z 位置色散不起作用的假设。此时方程退化为 du/dz = i N² |u|² u，解就是乘一个纯相位exp(i N² |u|² dz)。代码里把hhz = i * N^2 * deltaz预存，每个步长内的非线性相位就是abs(uu).^2 .* hhz。

问题在于，色散和非线性在数学上并不对易。如果每个大步长先做整段色散、再做整段非线性，误差是一阶 O(dz)。对称分步法把每个步长拆成“半非线性 → 色散 → 半非线性”，复合算子变成：

exp(h/2) · exp(D) · exp(h/2)

其中 h 表示非线性算子，D 表示色散算子，整体误差降到二阶 O(dz²)。代价是每个步长要算两次非线性指数，但对 FFT 次数没有影响。

回到代码，主循环体每轮做的是temp = uu .* exp(abs(uu).^2 .* hhz)，看起来是整段非线性。实际上由于循环外提前做了一个hhz/2，循环结束后又用-abs(uu).^2 .* hhz/2退掉最后一步的多余半相，整条 z 链路上每一段色散之前和之后都各有一半非线性相移。我一般建议读者自己加中间输出：把uu的峰值功率和脉宽在每个 step_num 节点打出来，观察 N=2 时周期性压缩，就能直观理解这个“预补偿-后补偿”设计。

3.3 fft/ifft方向：这段代码的约定可能和你习惯的相反

很多教科书上的写法是“时域→fft→频域乘色散→ifft→时域”。这段代码故意反着用：temp=ifft(uu)进频域，乘完 dispersion 后用fft回时域。从数学上讲，Matlab 的 fft 和 ifft 只是相差一个归一化和指数符号，只要配对正确结果完全一样。关键是不能在循环里混用。下面这个表可以对照：

如果你习惯fft(uu)进频域，那色散因子必须取共轭符号，回时域用ifft。否则脉冲会在每个步长被错误地反号，最终结果看起来像“噪声”，但又不是纯噪声。

另外，显示频谱时temp=fftshift(ifft(uu)).*(nt*dtau)/sqrt(2*pi)里的nt*dtau是连续傅里叶变换的数值近似系数，只影响纵轴幅度。fftshift只是为了让画图时零频在中间。千万不要把这个fftshift也搬到色散循环里。判断方法很简单：单独跑 beta2=-1、N=0（不带非线性）时，任意输入脉冲的频谱幅度应该保持不变，只改变相位；如果幅度变了，就是 FFT 方向或 fftshift 用错了。

4. 实操复现：跑通基阶孤子与N=2高阶孤子

4.1 基阶孤子验证：N=1, sech, beta2=-1

把 mshape 改成 0，N 保持 1，distance=10，输入就是双曲正割脉冲。基阶孤子的特征是形状沿 z 不变。由于数值模拟有窗口截断，你会看到时域波形几乎完全重合，频谱也基本不变；只有在时间窗边界有轻微杂散，因为 sech 尾部没完全衰减到 0。Tmax=32 对 1 个时间单位宽的脉冲足够大，杂散很小。

为了量化，可以在循环里记录峰值功率：

peak = zeros(1, step_num); for n = 1:step_num f_temp = ifft(temp) .* dispersion; uu = fft(f_temp); temp = uu .* exp(abs(uu).^2 .* hhz); peak(n) = max(abs(uu).^2); end figure; plot((1:step_num)*deltaz, peak); xlabel('z'); ylabel('Peak Power');

这段代码把每次循环后的瞬时峰值功率存下来。基阶孤子时这条线应该平稳在 1 附近。如果峰值缓慢下降，说明 Tmax 不够大或步长太大；如果峰值快速振荡，说明窗口或频率网格有问题。

4.2 参数改动与观察：N=2, 超高斯, chirp

一次只改一个参数，才能看出因果。下面是我常用的对照实验表：

比如把 mshape 改为 1、N=1，运行后输出频谱从窄高斯变成两侧陡峭的展宽谱，这正是自相位调制（SPM）的特征；如果再叠加上 chirp0=1，时域波形的不对称性会更明显。做实验时，建议把 figure(1) 的输入谱和 figure(2) 的输出谱用 hold on 画在一起，能直接看到谱宽变化方向。

对于 N=2，输入仍是sech(tau)但方程里的非线性系数变成 4，等效于标准 NLS 方程中幅度为 2 的二阶孤子。输出会看到脉冲在 z=0 附近先压缩，随后分裂再恢复，一个完整呼吸周期约 1.57。distance=10 约包含 6 个周期，能清楚地看到周期重复。

4.3 常见报错与结果异常排查

有几个高频问题，基本覆盖了大多数人跑这段源代码的踩坑点。

第一，sech函数不存在。Octave 或部分旧版 Matlab 没有内置 sech，需要自己在脚本开头定义：

sech = @(x) 1 ./ cosh(x);

第二，峰值功率发散或出现 NaN。通常是 dtau 太大导致abs(uu).^2 .* hhz的值超过浮点表示范围，或者超高斯脉冲的tau.^(2*mshape)在窗口边缘产生 inf。这时优先减小 mshape，比如从 3 降到 1，或者增大 Tmax 让窗口边缘的 tau 值变小。

第三，频谱不对称。先检查 omega 数组长度是否为 nt，再看循环内是否误用 fftshift。一个快速验证法：令 N=0、beta2=-1，输入任意脉冲跑一遍，频谱幅度应该前后完全一致。

第四，结果对 step_num 依赖明显。把 step_num 公式前的系数 20 提高到 40 对比，如果脉宽、峰值变化超过 1%，说明当前步长不够。对于强非线性场景如 N=3 以上，这个系数可能要提高到 50～100。

第五，step_num=0 报错。当 distance=0 或 N=0 时round()可能得到 0，for 循环直接不执行。如果 N=0 表示纯色散，建议单独处理：把循环里temp=uu.*exp(...)去掉，只保留色散步骤。

5. 进阶：用守恒量检查这段源代码的正确性

5.1 用能量守恒判断时间窗和步长

NLS 方程有一个重要守恒量：∫|u|²dτ，也就是脉冲总能量。分步傅里叶法的优点在于，非线性相位只在时域乘一个模长为 1 的复数指数，色散相位只在频域乘一个模长为 1 的复数指数，所以单步内能量理论上精确守恒。实际模拟中能量变化主要来自时间窗截断，而不是分步误差。

我一般会在主循环里加一个能量记录：

energy_z = zeros(1, step_num+1); energy_z(1) = sum(abs(uu).^2) * dtau; for n = 1:step_num f_temp = ifft(temp) .* dispersion; uu = fft(f_temp); temp = uu .* exp(abs(uu).^2 .* hhz); energy_z(n+1) = sum(abs(uu).^2) * dtau; end plot((0:step_num)*deltaz, energy_z/energy_z(1)-1);

如果曲线是数值噪声般的小波动，说明窗内没问题；如果呈单调下降，说明脉冲能量跑出了计算窗，优先增大 Tmax 而不是盲目加密步长。如果曲线随步长剧烈变化，就要检查 dispersion 里的 omega 网格是否和 dtau 匹配。

5.2 用半程输出代替逐帧画图

主循环里每步都画图会慢到不可接受。可以把plot移出循环，只保存若干 z 位置的场，例如：

save_step = round(step_num/20); uu_z = zeros(21, nt); k = 1; for n = 1:step_num f_temp = ifft(temp) .* dispersion; uu = fft(f_temp); temp = uu .* exp(abs(uu).^2 .* hhz); if mod(n, save_step) == 1 k = k + 1; uu_z(k,:) = uu; end end

循环结束后用imagesc(abs(uu_z).^2)画 z-τ 二维演化图，比只看输入输出两个截面信息量大得多。

如果以后要把这段源代码移植到 Python，注意 scipy.fft 的约定和 Matlab 一样，循环内同样不要用 fftshift；用 numpy 的ifft对应这里的ifft，fft对应fft。色散因子写成exp(1j*0.5*beta2*omega**2*deltaz)，hhz 同理，其余逻辑完全一致。最后再提醒一个容易忽略的点：当你看到谱边缘出现高频噪声时，最优先检查时间窗有没有截断 sech 尾巴，而不是步长。通常把 Tmax 从 32 加到 64，噪声就会消失。

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