简介：面向非线性光纤光学与数值仿真方向的学习者，这份资源提供用分步傅里叶法求解非线性薛定谔方程（NLS）的完整源代码解析，适合已具备Matlab基础、希望理解光脉冲在光纤中传输演化过程的高年级本科生与研究生。压缩包内共1个docx文档，约13KB，以文字讲解配合代码片段的形式，依次说明NLS方程形式、分步傅里叶法的非线性步与色散步拆分思路，以及输入参数、模拟参数、脉冲生成、主循环与结果可视化等模块的组织方式。读者可据此掌握色散相移因子与非线性相移因子的构造方法，理解1/2N→D→1/2N对称分步格式的实现细节，并对照sech与超高斯脉冲的时域、频域绘图验证仿真结果。目前已有1377人学习下载，可作为光纤光学课程设计或科研入门阶段理解分步傅里叶法数值实现的参考材料。

1. 一份能跑出孤子演化的 MATLAB 源码，到底解决了什么

做非线性光纤光学仿真的人，迟早会撞上同一个问题：想验证一段脉冲在光纤里怎么展宽、怎么被非线性压回来、什么时候裂成多个峰，手头却没有一个能直接改参数就跑的脚本。这份分步傅里叶法解 NLS 方程的 MATLAB 源码，就是干这件事的——它把非线性薛定谔方程拆成色散算符和非线性算符交替作用，用 FFT 在时域和频域之间来回切换，几十行代码就能把基孤子、高阶孤子、超高斯脉冲的传输演化跑出来。适合两类人：一类是刚进非线性光纤课题、需要快速拿到可视图的研究生；另一类是已经写过仿真但想对照一份干净参考实现的工程师。它不依赖任何工具箱，纯基础 MATLAB 语法，拿到就能跑，改参数就能换物理场景，这正是它值得拆开看的原因。

2. 方程与算法骨架：NLS 怎么被拆成两步算子

2.1 从 NLS 方程到归一化形式

源码第一行注释就把方程摆出来了：

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

这里 u 是脉冲包络，z 是传输距离，tau 是随脉冲移动的本地时间，beta2 是群速度色散，N 是非线性参数。注意源码里 beta2 直接取 -1，N 取 1，这对应反常色散区的基孤子条件 N=1。方程里 sgn(beta2) 的写法意味着色散项符号由 beta2 正负决定，改 beta2 的符号就能从反常色散切到正常色散，物理图像完全不同。N 的平方出现在非线性项前，源码里hhz = i*N^2*deltaz正是把这个系数吸收进了非线性相位因子，所以改 N 不需要动主循环，只改一个输入参数即可。

归一化的好处是长度、时间、功率都无量纲，distance=10 表示 10 个色散长度或孤子周期量级，具体物理单位由使用者自己映射。这种写法在文献里很常见，代价是初学者容易忘记自己到底在算什么尺度，后面避坑章节会专门说这一点。

2.2 对称分步傅里叶：为什么是 1/2N → D → 1/2N

分步傅里叶法的核心思想是把一个步长 deltaz 内的演化近似成两段：先让非线性单独作用半步，再让色散作用整步，最后再补半步非线性。源码注释写得很直白：

% scheme: 1/2N -> D -> 1/2N; first half step nonlinear temp = uu.*exp(abs(uu).^2.*hhz/2); % note 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); % Final field

这段是整个项目的灵魂。为什么用对称形式而不是简单的 N→D？因为对称分步傅里叶的局部误差是 O(deltaz^3)，比非对称的 O(deltaz^2) 高一阶，同样的步数下精度明显更好，这是数值方法里的经典结论。循环里每次做的是：当前场在时域乘非线性相位（整步），然后 ifft 到频域乘色散相位，再 fft 回时域。循环外的首尾两个半步非线性，加上循环内的整步非线性，合起来正好构成对称结构。

参数上，hhz/2和hhz的区别就是半步和整步，dispersion是预先算好的频域相位因子，避免在循环里重复算 exp。这种「预计算相位因子」的写法是性能优化的常见做法，尤其当 step_num 上千时，省下的 exp 调用很可观。

2.3 频域色散因子的构造

色散步在频域完成，因为 d^2/dtau^2 在频域就是乘以 -omega^2。源码这样构造：

omega = (pi/Tmax)*[(0:nt/2-1) (-nt/2:-1)]; % frequency grid dispersion = exp(i*0.5*beta2*omega.^2*deltaz); % phase factor

omega 的排列是 MATLAB fft 的标准顺序：前半是正频率 0 到 pi/Tmax，后半是负频率。这个顺序必须和 fft/ifft 的输出顺序严格对应，否则色散符号会错，脉冲会往反方向展宽——这是新手最常翻车的地方之一。0.5*beta2*omega.^2*deltaz里的 0.5 来自方程中二阶导前的 1/2 系数，deltaz 是步长，整体就是这一步的色散相位。

注意这里没有用 fftshift，因为 fft 本身输出的就是这种排列，只有在画频谱图时才需要 fftshift 把零频移到中间。源码在绘图时确实用了fftshift(omega)和fftshift(ifft(uu))，这个区分要拎清楚。

3. 参数与网格：nt、Tmax、step_num 怎么定

3.1 时间窗口与 FFT 点数的取舍

nt = 1024; Tmax = 32; % FFT points and window size dtau = (2*Tmax)/nt; % step size in tau tau = (-nt/2:nt/2-1)*dtau; % temporal grid

Tmax=32 表示时间窗口从 -32 到 +32，总宽 64，dtau 约 0.0625。nt=1024 是 2 的幂，FFT 效率最高。这两个参数不是随便定的：窗口要足够宽，让脉冲演化过程中展宽的能量不撞到边界；点数要足够密，让频谱分辨率够高。如果脉冲在传输中显著展宽，Tmax 太小会导致能量绕回（周期性边界条件），出现假的干涉条纹。

经验上，初始脉冲宽度如果是几个时间单位，Tmax 取脉冲宽度的 10 到 20 倍比较稳妥。源码里 sech(tau) 的宽度约 1，Tmax=32 给了充足余量。改脉冲形状或加 chirp 后如果发现边缘有异常振荡，第一反应就该是加大 Tmax 或 nt。

3.2 步长与步数的自动估算

step_num = round(20*distance*N^2); % No. of z steps to deltaz = distance/step_num; % step size in z

这里用20*distance*N^2估算步数，是个经验公式：N 越大非线性越强，需要的步数越多；distance 越长步数越多。20 这个系数保证每个孤子周期大约 20 步，对基孤子够用。但要注意，这个公式对强非线性或高阶孤子可能偏少，跑出来如果看到频谱不对称或能量不守恒，就该手动加大系数，比如改成 50 或 100 再对比。

deltaz 由 distance 除以 step_num 得到，保证总距离精确等于 distance。这种「先定步数再算步长」的顺序比「先定步长再算步数」更稳，因为后者可能因为取整导致总距离对不上。

3.3 输入脉冲的两种形态

if mshape==0 uu = sech(tau).*exp(-0.5i*chirp0*tau.^2); % soliton else uu = exp(-0.5*(1+i*chirp0).*tau.^(2*mshape)); end

mshape=0 走 sech 孤子分支，mshape>0 走超高斯分支。sech 是基孤子的标准形态，配 N=1、beta2=-1 就是严格的不变传输。超高斯分支里tau.^(2*mshape)的指数随 mshape 增大而变陡，mshape=2 时接近方波，常用于模拟带限脉冲或某些激光器输出。chirp0 是初始啁啾，非零时脉冲会先压缩或先展宽，取决于 chirp 符号和 beta2 符号的配合。

这里有个细节：超高斯分支没有归一化系数，幅度峰值固定为 1，而 sech 分支峰值也是 1，所以两种形态的峰值功率可比，但能量不同。做定量对比时要留意这一点。

4. 主循环逐段拆解：从半步非线性到最终频谱

4.1 首半步非线性与循环内整步

前面贴过主循环，这里把每一步的物理含义说透。temp = uu.*exp(abs(uu).^2.*hhz/2)是第一个半步非线性，abs(uu).^2是瞬时功率，乘 hhz/2 得到半步的非线性相移。注意这里用的是自相位调制（SPM）的相位，没有考虑拉曼或自陡峭等高阶效应，属于标准 NLS 范畴。

循环内f_temp = ifft(temp).*dispersion先做逆傅里叶变换到频域，乘色散相位，uu = fft(f_temp)再变回时域。这里 ifft 和 fft 的配对顺序不能反：因为 MATLAB 的 fft 定义里正变换带负号，逆变换带正号，而色散算符 exp(i*...) 的符号要和方程里的 i 对应。如果发现脉冲展宽方向反了，多半是 fft/ifft 用反了或者 omega 符号错了。

循环内temp = uu.*exp(abs(uu).^2.*hhz)是整步非线性，为下一步的色散做准备。整个循环执行 step_num 次，每次推进 deltaz。

4.2 尾半步与最终频谱计算

uu = temp.*exp(-abs(uu).^2.*hhz/2); % Final field temp = fftshift(ifft(uu)).*(nt*dtau)/sqrt(2*pi); % Final spectrum

尾半步用-abs(uu).^2.*hhz/2，负号是因为对称结构要求首尾两个半步的非线性相位方向相反，合起来抵消掉多算的部分。这一步很容易写错符号，写错的话总非线性相位会多算或漏算一个整步，孤子会明显变形。

最终频谱那行fftshift(ifft(uu)).*(nt*dtau)/sqrt(2*pi)做了三件事：ifft 变到频域、fftshift 把零频移到中间、乘nt*dtau/sqrt(2*pi)做归一化。这个归一化系数来自连续傅里叶变换和离散傅里叶变换的尺度关系，少了它频谱幅度没有物理意义，只能看形状不能看绝对值。绘图时用abs(temp).^2得到谱功率。

4.3 输入输出的可视化对照

源码画了两张图：figure(1) 是输入脉冲的时域和频域，figure(2) 是输出。时域用plot(tau, abs(uu).^2)，频域用plot(fftshift(omega)/(2*pi), abs(temp).^2)，横轴除以 2*pi 把角频率换成普通频率。axis([-5 5 0 inf])把时域视野限制在 ±5，因为脉冲主要能量集中在这里，全窗口画出来两边全是零，看不清细节。

这种「输入输出对照」的画法是调试利器：如果输出时域出现明显不对称，说明步长不够或边界反射；如果频谱出现不该有的旁瓣，可能是 Tmax 太小导致混叠。建议每次改参数后都先看这两张图再下结论。

5. 避坑与排查：五个真实会踩的坑

5.1 现象：脉冲能量随传输单调下降

原因：多半是 fft/ifft 配对或归一化系数写错，导致每次变换损失幅度。也可能是 deltaz 太大，非线性相位指数溢出。

解决：先检查ifft和fft是否成对出现且顺序正确，再确认dispersion里没有多余的归一化因子。把 step_num 临时加大 10 倍看能量是否守恒，如果守恒了就是步长问题，回到 3.2 节调系数。

5.2 现象：频谱出现对称的假边带

原因：Tmax 太小，脉冲展宽后能量绕到窗口另一端，周期性边界条件造成干涉。

解决：把 Tmax 加倍、nt 同步加倍，保持 dtau 不变，再看假边带是否消失。如果消失就是窗口问题。经验是 Tmax 至少取初始脉冲宽度的 15 倍。

5.3 现象：孤子传输几步就裂开

原因：N 和 beta2 的符号或数值不匹配。基孤子要求 N=1 且 beta2<0（反常色散），如果 beta2 取正，色散和非线性同向，脉冲会崩溃。

解决：确认 beta2 符号，反常色散区取负。再检查hhz = i*N^2*deltaz里的 i 是否保留，丢了 i 就变成纯幅度调制，物理全错。

5.4 现象：改 mshape 后超高斯脉冲形状不对

原因：tau.^(2*mshape)在 mshape 非整数时可能产生复数或负值幂，MATLAB 会给出警告或 NaN。

解决：mshape 取正整数，常用 1、2、3。如果要连续调节陡峭度，改用exp(-0.5*(tau/tau0).^(2*mshape))并显式给 tau0，避免 0 的负幂。

5.5 现象：绘图横轴范围不对，脉冲被截断

原因：axis([-5 5 0 inf])是硬编码的，脉冲展宽超过 ±5 后看不全。

解决：把 axis 范围改成[-Tmax/4 Tmax/4 0 inf]或根据实际展宽动态设置，别让绘图限制误导你对物理结果的判断。

6. 二次开发：从固定脚本到可复用仿真函数

这份源码是脚本形态，所有参数写死在开头，跑一次改一次。真要拿它做研究或工程验证，第一步就是把它包成函数，输入参数、输出场和频谱，这样才能批量扫参数。我一般会这样改：

function [uu_out, spec_out, tau, omega] = ssprop_nls(distance, beta2, N, mshape, chirp0, nt, Tmax, step_factor) % 分步傅里叶法解 NLS，返回输出场与频谱 % step_factor 控制步数系数，默认 20 if nargin < 8, step_factor = 20; end dtau = (2*Tmax)/nt; tau = (-nt/2:nt/2-1)*dtau; omega = (pi/Tmax)*[(0:nt/2-1) (-nt/2:-1)]; step_num = max(1, round(step_factor*distance*N^2)); deltaz = distance/step_num; dispersion = exp(1i*0.5*beta2*omega.^2*deltaz); hhz = 1i*N^2*deltaz; if mshape == 0 uu = sech(tau).*exp(-0.5i*chirp0*tau.^2); else uu = exp(-0.5*(1+1i*chirp0).*tau.^(2*mshape)); end temp = uu.*exp(abs(uu).^2.*hhz/2); for n = 1:step_num temp = fft(ifft(temp).*dispersion); temp = temp.*exp(abs(temp).^2.*hhz); end uu_out = temp.*exp(-abs(temp).^2.*hhz/2); spec_out = fftshift(ifft(uu_out)).*(nt*dtau)/sqrt(2*pi); end

改造要点有三个：一是把 step_factor 暴露出来，方便做步长收敛性测试；二是循环内用temp自更新，减少变量名混乱；三是输出 tau 和 omega，调用方不用重复构造网格。包成函数后，扫 N 从 1 到 3、扫 beta2 从 -2 到 -0.5，一个 for 循环就能出参数图，比改脚本快得多。

验证方法上，我习惯做两件事：一是能量守恒检查，无损耗 NLS 的总能量sum(abs(uu).^2)*dtau应该在整个传输中保持不变，偏差超过 1% 就说明步长或边界有问题；二是步长收敛测试，把 step_factor 从 20 加到 40、80，看输出峰值功率是否收敛，收敛了才敢信结果。这两步做完，再去看物理现象，心里才有底。

还有一个容易忽略的点：源码里clear all会清掉工作区所有变量，包成函数后要去掉，否则调用方环境被污染。另外figure(1)、figure(2)这种硬编码图号在批量跑时会互相覆盖，改成返回数据、由调用方决定画不画，复用性会好很多。

从那以后我每次改这类仿真脚本，都强制先跑一遍能量守恒和步长收敛，再谈物理结论——这两个检查花不了几分钟，却能挡掉大部分「看起来像新现象其实是数值伪影」的坑。希望帮到你。

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