简介面向非线性光纤光学与数值仿真学习者的分步傅里叶法求解NLS方程源代码文档适合已具备Matlab基础、希望理解光脉冲传输建模与数值实现的中高级读者。内容围绕非线性薛定谔方程展开包含输入参数、FFT点数与时间窗口、步长设置sech与超高斯脉冲生成以及非线性半步与色散步骤交替推进的核心写法并给出时域和频域结果可视化思路便于对照公式理解代码流程与调试要点。压缩包仅含1个docx文件约13KB篇幅精炼但覆盖从参数配置到主循环实现的关键环节。已有1376人学习可作为课程实验、光纤传输仿真入门及代码复现的参考材料。1. 1024 点 FFT 替代偏微分方程分步傅里叶法解 NLS 的结论分步傅里叶法解非线性薛定谔方程NLS 方程的源代码里最反直觉的一点是没有迭代求解偏微分方程也没有组装矩阵核心循环只有两次 FFT 和两次逐点乘法。光纤里的色散和非线性被拆成两个算子分别在自己的表象里精确求解再用对称半步拼起来。这套做法的误差不是来自“算不动”而是来自两个算子不对易理解这一点才能明白为什么step_num要按N^2放大、为什么Tmax和nt一改频谱就变形。这份 MATLAB 源代码把分裂格式、频域相移、时域非线性相移和频谱归一化都压在一个不到百行的主循环里读透它等于把分步傅里叶法的底层账本翻了一遍。2. 算子分裂色散和非线性拆开算凭什么不出错2.1 Strang 分裂的 1/2N - D - 1/2N 顺序NLS 方程可以写成算子形式∂u/∂z (D N)u其中D是线性色散算子N是非线性算子。精确走一步h的解是exp(h(DN))u但D和N不对易这个指数不能直接拆成两个指数相乘。分步傅里叶法的做法是用乘积近似exp(h(DN)) ≈ exp(hD/2) exp(hN) exp(hD/2)这就是 Strang 分裂。它把一步分成“半个色散、整个非线性、半个色散”或者对称地写成“半个非线性、整个色散、半个非线性”。两种写法误差同阶代码选后者因为非线性在时域是逐点乘法色散在频域是逐点乘法每步只需要两次 FFT。% 对称分裂半步非线性 - 整步色散 - 半步非线性 temp uu .* exp(abs(uu).^2 .* hhz/2); % 1/2 N for n 1:step_num f_temp ifft(temp) .* dispersion; % D频域逐点乘 uu fft(f_temp); % 回到时域 temp uu .* exp(abs(uu).^2 .* hhz); % 整步 N end uu temp .* exp(-abs(uu).^2 .* hhz/2); % 扣掉多做的半步 Nhhz i*N^2*deltaz是非线性相移系数abs(uu).^2是瞬时功率。初始temp先做半步非线性循环里每次先ifft到频域乘色散因子再fft回时域做整步非线性循环结束后再乘exp(-abs(uu).^2.*hhz/2)把最后一次多做的半步扣掉。这样总分裂正好是N(h/2) D(h) N(h) ... D(h) N(h/2)每一步的色散只做一次非线性首尾各半。注意exp(-abs(uu).^2.*hhz/2)里的负号不是“反向传播”而是把循环内最后一次整步非线性修正回半步。漏掉这一项总非线性相移会多出半步长距离传播后相位明显偏移。2.2 局部误差与全局误差deltaz 的三次方去哪了Strang 分裂的局部误差是O(h^3)全局误差是O(h^2)。也就是说步长减半误差理论上降到四分之一。这个结论来自 Baker-Campbell-Hausdorff 公式两个算子乘积的对数里一阶项是h(DN)二阶项包含交换子[D,N]对称分裂让二阶项抵消剩下的三阶项才是误差主项。误差系数和交换子大小有关。色散算子D在频域是i β2 ω^2 / 2非线性算子N在时域是i N^2 |u|^2。交换子会同时带上ω^2和|u|^2所以频率越高、功率越大、N越大同样步长下的误差越大。代码里step_num round(20*distance*N^2)是一个经验公式每单位长度至少20N^2步。对N1、distance10步数是 200deltaz0.05峰值非线性相移约N^2|u|^2h 0.05弧度对N3步数变成 1800步长缩到约0.0056峰值非线性相移仍控制在 0.05 弧度附近。这里容易误解的一点是色散相移可以很大步长却不必按色散相移来缩。频域乘exp(i β2 ω^2 h/2)是精确解色散相移再大也不产生分裂误差误差来自“色散和非线性不能同时对角化”。所以判断步长够不够看的是非线性相移和交换子而不是频域相位转了多少圈。血泪经验是N增大时只把step_num按N^2放大还不够如果脉冲峰值功率远大于 1还要看实际max(abs(uu).^2)必要时再缩半步。2.3 为什么色散在频域是逐点相乘傅里叶变换把时间导数变成频率乘法∂²u/∂τ² - -ω² ũNLS 方程里的色散项是-sgn(beta2)/2 ∂²u/∂τ²频域就变成sgn(beta2)/2 ω² ũ。走一步h的色散解是ũ(zh, ω) ũ(z, ω) exp(i β2 ω² h / 2)代码里预计算dispersion exp(i*0.5*beta2*omega.^2*deltaz)主循环里只做一次逐点乘法。这里用的是beta2本身不是sgn(beta2)因为这份代码的归一化默认|beta2|1两者等价。如果把物理光纤的beta2-20直接填进去色散因子会变成exp(-i*10*ω²h)和方程里的sgn(beta2)不一致波形会明显跑偏。常见的实现思路是要么保持|beta2|1的归一化要么把代码改成sgn(beta2)并同步调整非线性系数。另一个细节是 FFT 方向。MATLAB 的ifft指数是ifft指数是-i。代码用ifft做正变换、fft做逆变换所以omega数组的符号约定要和ifft对齐。色散因子是ω²偶函数正负频率对称看起来“怎么换都行”但一旦加入线性项比如自频移或拉曼项符号错了就会翻车。读这段代码时最好把ifft和fft的方向当成固定约定不要随手互换。3. 网格与频谱FFT 点数、窗口和频率轴的对齐关系3.1 tau 网格与 omega 网格的构造代码里时间网格和频率网格是这样来的nt 1024; Tmax 32; dtau (2*Tmax)/nt; tau (-nt/2:nt/2-1)*dtau; omega (pi/Tmax)*[(0:nt/2-1) (-nt/2:-1)];tau从-Tmax走到Tmax-dtau共nt个点窗口长度是2*Tmax64。omega不是自然顺序而是 FFT 顺序先是0, dω, 2dω, ...的正频率再是-nt/2*dω, ..., -dω的负频率。dω 2π/(nt*dtau) π/Tmax。这种顺序和fft/ifft的输出顺序一致传播时可以直接逐点乘绘图时才用fftshift把零频移到中间。参数默认值含义影响nt1024FFT 点数同时决定频率分辨率和ω_maxTmax32半窗口决定时域截断和dωdtau0.0625时间步决定ω_max π/dtaudω约 0.0982频率步决定频率分辨率step_num200z 方向步数决定分裂误差deltaz0.05z 步长越小越准越慢如果tau用(0:nt-1)*dtau构造时域原点会跑到窗口左边频谱会多出一个线性相位因子幅度谱看起来不变相位和线性项却会错。代码用(-nt/2:nt/2-1)把零时刻放在窗口中心fftshift后频谱也以零频为中心读图更直观。3.2 fft/ifft 的归一化与频谱幅度MATLAB 的ifft默认带1/nt归一化fft不带。主循环里f_temp ifft(temp) .* dispersion; uu fft(f_temp);ifft的1/nt和fft的无归一化正好抵消所以色散步不改变幅度。这是主循环能逐点乘而不担心幅度的原因。频谱计算则不一样temp fftshift(ifft(uu)) .* (nt*dtau)/sqrt(2*pi);ifft(uu)含1/nt乘nt*dtau后变成dtau相当于连续傅里叶变换的离散近似。sqrt(2*pi)是傅里叶变换约定用不用它只影响频谱纵轴刻度不影响传播。如果自己写频谱图忘了这个因子峰值高度会和理论值差一个常数如果忘了fftshift零频会出现在图的两端看起来像“频谱折叠”。能量检查可以用 Parseval 对照E_time sum(abs(uu).^2) * dtau; E_freq sum(abs(temp).^2) * (pi/Tmax) / (2*pi);E_time和E_freq应该近似相等。temp已经fftshift但绝对值平方不受移位影响。两者差得多通常不是归一化写错而是时域窗口截断或频域混叠。3.3 采样窗口与混叠Tmax 怎么选时域窗口2*Tmax对应频域采样间隔dωπ/Tmax时间步dtau对应频域最大角频率ω_maxπ/dtau。这两个量通过 本文还有配套的精品资源点击获取