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

简介分步傅里叶法解非线性薛定谔方程的源代码面向光纤通信、非线性光学方向的研究生与工程师用于模拟光脉冲在光纤中的传输演化。包内共1个docx文档压缩包仅13KB文档内嵌完整Matlab源代码包含输入参数设置、FFT网格划分、色散相移因子计算、非线性/色散交替迭代以及输入输出脉冲时域频域绘图等模块。代码结构紧凑核心思路为将传输过程拆分为非线性步与线性色散步通过傅里叶变换快速求解并给出了sech脉冲与超高斯脉冲两种初始条件下的运行示例。已有1367人学习下载适合正在学习分步傅里叶法或需要快速搭建NLS仿真程序的读者。结合代码解析可帮助理解β2色散参数、非线性系数N、啁啾chirp等因素对脉冲演化与频谱展宽的影响省去从零编写和调试的时间。1. 分步傅里叶法解NLS方程先弄清楚这段源代码在算什么非线性光纤里光脉冲的传输几乎都可以用非线性薛定谔方程NLS来描述。真正动手模拟之前很多人第一反应是直接做分步傅里叶法然后栽在符号约定、FFT方向、步长选取这些细节上。这段源代码给了一个可以直接跑的参考实现用对称分步傅里叶法把一个双曲正割或超高斯脉冲沿光纤推10个色散长度并画出输入输出脉冲的时域波形和频谱。适合刚接触光纤非线性模拟的研究生、光通信链路工程师以及想验证自己另写的 C/Python 求解器的人。我建议先跑通默认参数再动手改 chirp、孤子阶数 N观察基阶孤子传输不变形、N2 周期性呼吸这些经典现象。理解了这一段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 的正常色散区则对应暗孤子和正常色散展宽。代码默认 distance10β2-1N1mshape1。这里要特别注意方程里写的是 sgn(β2)/2但代码的 dispersion 因子里直接用 beta2相当于把 β2 的符号和大小都放进去了。做参数扫描时如果只把 beta2 改成 -5色散项强度也会跟着变距离轴含义会不一样不要把它简单当成“换个符号”。下面这张表把代码开头的输入参数和模拟参数一次说清。参数默认值作用改参数会影响什么distance10光纤长度归一化单位总演化距离决定能否看到孤子周期beta2-1色散参数带符号色散相移大小与方向β20 为反常色散N1非线性参数/孤子阶数非线性相移强度N2 时出现周期性呼吸mshape1脉冲形状0 为 sech0 为超高斯光谱形状和初始波形的边缘陡峭程度chirp00输入脉冲线性啁啾不为 0 时脉冲先压缩或展宽nt1024FFT 点数时间窗内采样密度点数不足会混叠Tmax32时间窗半宽窗口总宽度必须远大于脉冲宽度step_numround(20distanceN²)纵向总步数步长精度N 越大需要越多步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 输出顺序一致。这样后面 tempfftshift(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;默认 distance10N1 时 step_num200deltaz0.05。N2 时同样是 distance10step_num800步长缩小到 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 时 dispersionexp(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 节点打出来观察 N2 时周期性压缩就能直观理解这个“预补偿-后补偿”设计。3.3 fft/ifft方向这段代码的约定可能和你习惯的相反很多教科书上的写法是“时域→fft→频域乘色散→ifft→时域”。这段代码故意反着用tempifft(uu)进频域乘完 dispersion 后用fft回时域。从数学上讲Matlab 的 fft 和 ifft 只是相差一个归一化和指数符号只要配对正确结果完全一样。关键是不能在循环里混用。下面这个表可以对照约定时域到频域频域回时域色散因子本代码ifftfftexp(i0.5beta2omega²deltaz)常用写法fftifftexp(-i0.5beta2omega²deltaz)如果你习惯fft(uu)进频域那色散因子必须取共轭符号回时域用ifft。否则脉冲会在每个步长被错误地反号最终结果看起来像“噪声”但又不是纯噪声。另外显示频谱时tempfftshift(ifft(uu)).*(nt*dtau)/sqrt(2*pi)里的nt*dtau是连续傅里叶变换的数值近似系数只影响纵轴幅度。fftshift只是为了让画图时零频在中间。千万不要把这个fftshift也搬到色散循环里。判断方法很简单单独跑 beta2-1、N0不带非线性时任意输入脉冲的频谱幅度应该保持不变只改变相位如果幅度变了就是 FFT 方向或 fftshift 用错了。4. 实操复现跑通基阶孤子与N2高阶孤子4.1 基阶孤子验证N1, sech, beta2-1把 mshape 改成 0N 保持 1distance10输入就是双曲正割脉冲。基阶孤子的特征是形状沿 z 不变。由于数值模拟有窗口截断你会看到时域波形几乎完全重合频谱也基本不变只有在时间窗边界有轻微杂散因为 sech 尾部没完全衰减到 0。Tmax32 对 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 参数改动与观察N2, 超高斯, chirp一次只改一个参数才能看出因果。下面是我常用的对照实验表实验编号修改项预期现象1mshape0, N1, chirp00波形、频谱不随 z 变化2mshape0, N2, chirp00周期性呼吸归一化孤子周期约 π/2≈1.573mshape1, N1, chirp00超高斯脉冲频谱快速展宽时域出现多峰结构4mshape0, N1, chirp01正啁啾使脉冲先展宽再压缩频谱宽度明显变大比如把 mshape 改为 1、N1运行后输出频谱从窄高斯变成两侧陡峭的展宽谱这正是自相位调制SPM的特征如果再叠加上 chirp01时域波形的不对称性会更明显。做实验时建议把 figure(1) 的输入谱和 figure(2) 的输出谱用 hold on 画在一起能直接看到谱宽变化方向。对于 N2输入仍是sech(tau)但方程里的非线性系数变成 4等效于标准 NLS 方程中幅度为 2 的二阶孤子。输出会看到脉冲在 z0 附近先压缩随后分裂再恢复一个完整呼吸周期约 1.57。distance10 约包含 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。一个快速验证法令 N0、beta2-1输入任意脉冲跑一遍频谱幅度应该前后完全一致。第四结果对 step_num 依赖明显。把 step_num 公式前的系数 20 提高到 40 对比如果脉宽、峰值变化超过 1%说明当前步长不够。对于强非线性场景如 N3 以上这个系数可能要提高到 50100。第五step_num0 报错。当 distance0 或 N0 时round()可能得到 0for 循环直接不执行。如果 N0 表示纯色散建议单独处理把循环里tempuu.*exp(...)去掉只保留色散步骤。5. 进阶用守恒量检查这段源代码的正确性5.1 用能量守恒判断时间窗和步长NLS 方程有一个重要守恒量∫|u|²dτ也就是脉冲总能量。分步傅里叶法的优点在于非线性相位只在时域乘一个模长为 1 的复数指数色散相位只在频域乘一个模长为 1 的复数指数所以单步内能量理论上精确守恒。实际模拟中能量变化主要来自时间窗截断而不是分步误差。我一般会在主循环里加一个能量记录energy_z zeros(1, step_num1); 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(n1) 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对应这里的ifftfft对应fft。色散因子写成exp(1j*0.5*beta2*omega**2*deltaz)hhz 同理其余逻辑完全一致。最后再提醒一个容易忽略的点当你看到谱边缘出现高频噪声时最优先检查时间窗有没有截断 sech 尾巴而不是步长。通常把 Tmax 从 32 加到 64噪声就会消失。本文还有配套的精品资源点击获取