企业官网建设流程全解析

简介：分步傅里叶法解非线性薛定谔方程的源代码，面向非线性光纤光学与光通信仿真方向的初学者和研究者，系统讲解如何用分步傅里叶法数值求解NLS方程。包内为1个docx文档，压缩包仅13KB，文档完整给出MATLAB程序及其逐段解析，内容覆盖输入参数设置、模拟参数设置、输入脉冲profile生成、分步傅里叶法实现与结果可视化五大模块，并说明参数对仿真结果的影响；其中输入参数包括光纤长度、色散参数、非线性参数、脉冲形状与chirp等，模拟参数包括FFT点数、时间窗口、z方向步长等。非线性步骤通过temp = uu.*exp(abs(uu).^2.*hhz/2)实现，色散步骤则利用ifft(temp).*dispersion完成，文档对hhz、dispersion等关键相移因子的构造逻辑进行了说明，并展示了输入输出脉冲在时域和频域上的对比可视化方法。读者可据此快速掌握分步傅里叶法的对称分步编程框架，并迁移到孤子传输、超短脉冲演化等仿真场景；已有1368人学习下载，适合需要理解NLS方程数值实现细节的科研人员与学生参考。

1. 分步傅里叶法解NLS方程：一份能直接改参数跑通的光纤仿真源码

做超短光脉冲传输仿真时，最烦人的不是方程本身，而是画出来的波形和频谱怎么看都不对劲——脉冲分裂、频谱漂移、能量乱跳，最后发现大多是分步傅里叶法里那几个看似不起眼的设置：步长取多大、FFT窗口开多宽、频域网格怎么排。分步傅里叶法，正是解非线性薛定谔方程（NLS）最主流的数值手段，说白了就是绕开“非线性项和色散项耦合在一起没法直接积分”的难题，把问题拆成两步：在频域做色散相乘，在时域做非线性相乘，交替推进。这份MATLAB源码恰好把这些细节都落在了一百多行里，适合光通信、非线性光纤光学方向的学生和工程师作为起点，改改参数就能复现孤子演化、脉冲展宽、超连续谱产生这类现象。

2. 先把输入参数读透：distance、beta2、N 和网格的连带关系

分步傅里叶法看起来是个纯数值算法，但它的第一关其实是物理参数映射到仿真参数。源码开头那段看起来只是赋值，实际上每一步都影响最终结果的物理解释，不能跳过不管。

2.1 物理参数组：光纤长度、色散、非线性阶数

代码第一段给出了五个物理量，如下：

distance = 10; % 光纤长度（归一化单位） beta2 = -1; % 群速度色散参数，负号为反常色散 N = 1; % 非线性参数，对应孤子阶数 mshape = 1; % 脉冲形状：0 为双曲正割，>0 为超高斯 chirp0 = 0; % 初始线性啁啾，0 表示无啁啾

这里要特别解释beta2的符号。NLS 方程原文写作i du/dz - sgn(beta2)/2 d^2u/dtau^2 + N^2|u|^2u = 0，代码里直接用了负的beta2，意味着默认是反常色散区，即普通光纤在近红外波段的情况。如果你仿真正常色散区，把beta2改成正值即可。

N不是单纯的非线性强度，它是孤子阶数的标志：N=1对应基阶孤子，脉冲在传输中形状基本不随距离改变；N=2以上会出现周期性压缩和分裂，也就是高阶孤子效应。初学阶段建议从N=1起步，先把程序的守恒性验证通过，再调高数值观察孤子周期。

这段参数的典型取值经验是：distance应与色散长度差不多量级。色散长度为T0^2/|beta2|，在归一化单位体系里，如果初始脉宽大约是 1 个时间单位、|beta2|=1，色散长度大约也是 1。代码里distance=10已经足够让非线性效应和色散效应充分作用几十个周期，适合用来观察孤子周期演化。

2.2 仿真网格：nt、Tmax、dtau 和 omega 的排列方式

网格是分步傅里叶法最容易翻车的地方。源码里的设置是：

nt = 1024; % FFT 点数 Tmax = 32; % 时间窗半宽 step_num = round(20*distance*N^2); % z 方向总步数 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)]; % 频率网格

先解释时间网格：tau从-Tmax到Tmax-dtau，共nt个点，时间窗口宽度是2*Tmax=64。采样率dtau越小、窗口越宽，频谱分辨率越高，但计算量也越大。对于孤子仿真，nt=1024相对富余，可以让光谱边缘的微弱分量不被窗口截断。

omega这一行值得单独讲。MATLAB 的fft输出频率顺序是直流在第一个位置，然后正频率一直到奈奎斯特频率，再是负频率从负奈奎斯特到接近零。代码里omega = (pi/Tmax)*[(0:nt/2-1) (-nt/2:-1)]就是按这个顺序构造的角频率数组，和fft的结果天然对应。

我见过不少人图省事，直接写omega = linspace(-pi/Tmax, pi/Tmax, nt)，这样虽然看着对称，但和fft的输出顺序是错位的。一旦错位，色散相位因子exp(i*0.5*beta2*omega.^2*deltaz)就会把不同频率分量挪到错误的位置上，表现出的症状是：初始脉冲形状看不出问题，但传播几步后波形在时域边缘出现锯齿状振荡，频谱向一侧偏。为了避免这个坑，我一般在生成omega后先做一个自检：对某个已知函数做fft，画fftshift(omega)对应的频谱，确认峰值落在预期的频率位置。

步数step_num的设定也隐含一个数值经验：round(20*distance*N^2)。这里放大N^2倍，是因为非线性相位随N^2增大而增大，要保证每步非线性相移足够小。若把N从 1 改到 4，而步数不增加，自相位调制造成的高频分量会迅速增长，频谱尾部容易出现非物理的振荡。

2.3 输入脉冲生成：sech 孤子与超高斯脉冲的区别

源码用mshape分支生成两种典型的入射脉冲：

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

双曲正割脉冲sech(tau)是无啁啾基阶孤子的精确解形式，所以验证程序时一定用mshape=0。超高斯脉冲的边缘更陡峭，上升沿接近矩形，适合模拟实际激光器输出的非理想脉冲。

这里的chirp0是初始线性啁啾系数，指数里的-0.5i*chirp0*tau^2表示脉冲瞬时频率随时间线性变化。啁啾对孤子的影响很直接：正啁啾会使脉冲先展宽再压缩，负啁啾则相反。如果只想复现教科书上的标准孤子演化，务必保持chirp0=0。

还有一个细节：超高斯脉冲的表达式里指数是tau.^(2*mshape)，而不是平方。mshape=1是标准高斯，mshape>1时脉冲顶部越来越平坦。传输时超高斯脉冲会由于色散和非线性的共同作用产生比双曲正割脉冲更复杂的频谱结构，适合用来观察自相位调制导致的光谱振荡。

3. 核心循环：非线性与色散如何交替推进

分步傅里叶法的核心代码量不大，难的是理解每一步在做什么、为什么这样排。把这段逻辑理顺，后面调参和排错都有依据。

3.1 循环外预计算色散相位因子与非线性相位因子

在进入主循环之前，代码先算好了两个关键的相位乘子：

dispersion = exp(i*0.5*beta2*omega.^2*deltaz); % 色散相位因子 hhz = i*N^2*deltaz; % 非线性相位因子

dispersion是线性色散在频域的传递函数。把光场变换到频域后，每传播deltaz距离，各频率分量只增加一个与omega^2成正比的相位。注意这里系数是0.5，原因是后面采用了对称分步格式——色散算子以二分之一的形式出现在一个完整步长的两侧，先走半个步长的色散，再走整步的非线性，最后再走半个步长的色散，这样整体的局部截断误差比“先非线性后色散”的非对称格式小一个数量级。

hhz则是每步非线性效应的相位贡献，正比于|u|^2。之所以把它提出去，是因为在每次迭代里都需要直接使用，预先计算能省掉重复乘法的开销。对于step_num达到几千步的仿真，这个预计算能让总耗时下降明显，尤其当nt较大时效果更明显。

3.2 对称分步循环：为什么第一个半非线性步放在循环外面

主循环的写法初看有点别扭，第一次读这份代码的人大多会问：为什么循环里只看到dispersion和整步非线性，那个hhz/2却出现在循环外？

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); % 最后一步：补上半个非线性步

仔细拆一遍：光场先被施加半个非线性相位，然后进入循环。每轮循环里，先把时域场变到频域，乘上色散相位因子，再变换回时域，接着施加整步的非线性相位。循环结束后再补上结尾的半个非线性相位。整体效果等价于：1/2 非线性 -> 色散 -> 1/2 非线性，如此反复。这种对称安排让非线性相位平均分到每步的两端，等效于在z+delta z/2处计算非线性效应，可以显著减小步进误差。

f_temp = ifft(temp).*dispersion;这一句是分步傅里叶法的“灵魂”：因为色散在时域是二阶微分算子，直接做数值微分不但慢而且精度差；但变换到频域后它就退化成乘法运算，速度和精度都理想。ifft和fft成对出现，保证计算始终在正确的域中完成。

循环内temp = uu.*exp(abs(uu).^2.*hhz);用到了当前步的光场强度abs(uu).^2，即自相位调制效应与强度成正比。这里的.*是逐元素乘法，不是矩阵乘法，一旦误写成*，MATLAB 会尝试做矩阵乘法，轻则报错，重则算出一个莫名其妙的“结果”。

最后一个半非线性步里的相位用的是-abs(uu).^2，仔细看会发现它和第一个半非线性步的相位符号一致——都是正的+号。这里看似不对称，实际上循环内最后一步已经施加了整步非线性，结束后需要再补半步行成完整的对称结构。负号本身只表示复数指数中的虚部符号，对应非线性折射率的正负号，表达式整体保持连续。

3.3 频谱输出的归一化：乘 dtau 还是不乘

程序里两次画频谱图都用了同样一行换算：

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

这行的作用是把时域波形变成物理意义明确的频谱幅度。MATLAB 的ifft默认归一化成“平均”形式，而连续傅里叶变换的离散近似需要乘采样间隔dtau。nt*dtau = 2*Tmax，正好是时间窗宽度，所以这行等价于乘上窗口总宽度再除以sqrt(2*pi)。如果不乘dtau，画出来的频谱幅度会随窗口大小变化，换一个Tmax结果就完全对不上。

fftshift的作用是把前面构造的“直流在两端”的频率顺序换成“直流在中间”的自然顺序，这样画图时频谱才能在横轴零频处对称。注意这里是对ifft的结果做fftshift，而不是对omega数组做。若弄反了，频谱会和横轴错位半个频点间距，虽然幅度形状不变，但峰值频率读出来总是差一点。

4. 避坑与常见问题：波形振荡、频谱漂移、能量不守恒的排查

这份代码我在不同版本的MATLAB上跑过，也见周围同学改参数时踩过不少坑。下面几条是最常见的，每一条都按照“现象、原因、解决”顺序写清楚。

4.1 传播几步后波形边缘出现锯齿状振荡

现象：初始脉冲很光滑，但从第几步开始时域波形两侧出现密集的细碎振荡，且振荡幅度随步数增加变大，看起来像噪声。

原因：最常见的是step_num取得太小，导致每步非线性相位过大。分步傅里叶法要求每步内|u|^2 * N^2 * deltaz远小于 1，当N=3以上还沿用step_num = 20*distance，自相位调制产生的频谱分量在单个大步长内没有被充分解析，混叠成时域振荡。

解决：把步数改成round(50*distance*N^2)或更高，观察振荡是否消失。同时对照deltaz与非线性长度1/N^2的数量级，一般取deltaz < 1/(5*N^2)比较稳妥。我通常先跑一个短距离测试，把distance临时设为 1，确认波形正常后再拉长传播距离。

4.2 频谱峰值位置总是偏移

现象：输入是原点对称脉冲，画出的频谱峰值却不在零频，而是偏左或偏右一个固定的量；改变nt时偏移比例变化。

原因：频率网格构造顺序与fft输出顺序不一致。最常见的误用是omega = linspace(-pi/Tmax,pi/Tmax,nt)，或者先定义对称网格再用fft而不是ifft，导致频域相乘时频率对不上。

解决：严格按照源码里的顺序构造omega，并把色散相位因子里的omega换成fftshift(omega)进行自检对照。如果两种方式结果相同，说明当前代码顺序正确；如果不同，就用源码方式。另外要注意，对fft的结果做fftshift与对频率数组本身做fftshift是两回事，必须配合使用。

4.3 能量随时间不守恒

现象：计算sum(abs(uu).^2)*dtau发现总能量随仿真步进单调下降或上升，且误差超过百分之一。

原因：分步傅里叶法保持能量守恒的前提是每一步都严格遵循对称格式。若把主循环写成“非线性-色散”交替而少了开头或结尾的半步行，或者循环内误把hhz/2与hhz都用上，都会引入额外相位，导致能量漂移。另一个常见原因是时间窗口太小，脉冲展宽后溢出边界，造成能量泄漏。

解决：先用基阶孤子N=1、mshape=0做测试，打印每步能量，理想情况下误差波动应在1e-10量级。若误差逐步累积，检查是否漏掉循环外的hhz/2步骤；若能量在某个特定步数突然下降，检查时域波形是否碰到窗口边缘，适当增大Tmax。

4.4 矩阵乘法与逐元素乘法混用

现象：运行时报错“矩阵维度必须一致”，或者不报错但结果全然不对。

原因：uu和exp(abs(uu).^2.*hhz)都是行向量，用*会触发矩阵乘法，维度对不上时直接报错；维度刚好能对上时（比如方阵情形），会执行矩阵乘法而不是逐元素相乘，结果面目全非。

解决：养成习惯，对数组运算一律用.*。在核心循环里，凡是涉及uu、temp、dispersion的乘法，都检查一遍是否为逐元素运算。调试时可以用whos查看变量尺寸，再断点执行看乘积维度是否保持为1 x nt。

4.5 chirp0 不为零时结果与理论不符

现象：设chirp0=1后，得到的脉冲峰值功率变化趋势和书本上的啁啾孤子演化对不上。

原因：chirp0的作用是给初试脉冲加上线性啁啾，强力改变孤子的演化行为。但注意指数里的-0.5i*chirp0*tau.^2中i的位置，如果误写成-0.5*chirp0*tau.^2，啁啾就成了实指数，相当于改变脉冲幅度分布，结果自然不会对。

解决：检查exp(-0.5i*chirp0*tau.^2)的写法，确认虚数单位确实在指数里。构造一个无啁啾基准，再对比加入 chirp 后的峰值功率曲线，这样可以直观看到啁啾对孤子演化的影响是否符合预期。

5. 进阶验证：用基阶孤子周期和能量守恒确认程序可靠

代码能够跑通只是第一关，真正让这份源码变得可信的，是拿已知解析特征去验证数值结果。我拿到这份源码之后，做的第一件事不是改参数，而是把N设为 1、mshape设为 0、chirp0设为 0，然后观察一段足够长的距离，确认基阶孤子的行为。

基阶孤子理论预期有两个特征：第一，传输过程中脉冲形状保持双曲正割不变；第二，孤子相位随距离线性增加，单个时间周期内累积相移为pi/2。实际操作中，我记录每一步的峰值功率和脉冲宽度，程序跑完后检查这两个量是否恒定。为了避免手写检查代码，我直接在源码的循环后面追加一段监测代码：

% 记录峰值功率与有效宽度，验证孤子稳定性 peak_power = max(abs(uu).^2); width_rms = sqrt(sum(abs(tau).^2 .* abs(uu).^2) / sum(abs(uu).^2)); energy = sum(abs(uu).^2) * dtau; fprintf('peak=%.6f width=%.6f energy=%.6f\n', peak_power, width_rms, energy);

如果峰值功率在1e-6精度内保持不变、能量守恒误差低于1e-10，说明色散和非线性的平衡被程序精确复现了。接下来才是调N=2观察孤子周期压缩效应、调distance观察脉冲分裂，一步步把这份源码当成自己的仿真工具箱来用。

还有一个小技巧是用N值估算步长：把deltaz和1/N^2的比值打印出来，若大于 0.2，果断把step_num翻倍。这个比值相当于每步非线性相位的“危险系数”，控制在 0.1 以下基本就不会出现前面提到的振荡。

从那以后我每次调参数，都强制走一遍峰值功率、宽度和能量守恒三项监测再画图，再也没在步长和网格上翻过车。希望帮到你。

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