非线性薛定谔方程数值求解：分步傅里叶法从原理到MATLAB实践

简介：这是一份非线性薛定谔方程（NLSE）的MATLAB数值模拟代码，面向光学、凝聚态物理领域的研究生与科研人员，适用于光孤子传输、超快光学现象以及Bose-Einstein凝聚体动力学研究。包内仅含1个m文件，压缩包大小约1KB，代码精简，内含势能项与非线性系数的可调参数，便于拓展至不同物理场景。NLSE在光纤通信与超冷原子领域有广泛应用，该程序通常基于分步傅里叶法或有限差分法求解，可复现光脉冲自相位调制、孤子稳定传播与碰撞、BEC中的振荡和分裂等典型现象，并能输出波函数演化过程，便于可视化分析。已有399人学习下载，适合希望快速上手非线性偏微分方程数值模拟的初学者，也适合作为相关课题的验证工具。通过阅读该代码，可系统理解方程离散化、边界条件设置及迭代求解的关键步骤，为后续设计光开关、光存储器件或多分量耦合模型打下基础。

1. 非线性薛定谔方程：这份 NLSE.m 到底在算什么

你手上这份NLSE.zip，核心就一个NLSE.m。很多人看到“薛定谔三个字，第一反应是量子力学，但实际拿到手里你会发现，它解决的是光学里非常常见的一类问题：脉冲在光纤里怎么传、怎么变成孤子、两个孤子撞在一起会发生什么。非线性薛定谔方程（NLSE）是这些现象的统一模型，而NLSE.m这种 MATLAB 程序，绝大多数情况下是在用分步傅里叶法做数值演化，把解析解推不出来的场景直接算给你看。

我拆过不少类似的光通信和冷原子项目，这类代码的典型画像很一致：一个主脚本管初值，一个函数管分步演化，最后画强度图和频谱图。它能帮你在一个小时之内把“孤子长什么样、色散和非线性怎么平衡”从抽象概念变成能动的曲线。适合刚接触 NLSE 的研究生，也适合要从实验数据反推参数的光通信工程师。下面我按实际拆代码的顺序，把它讲透。

2. 从方程到代码：分步傅里叶法的实现骨架

2.1 先做归一化：为什么代码里看不到 hbar 和 m

教科书里 NLSE 长这样：

[ i\hbar \frac{\partial \psi}{\partial t} = -\frac{\hbar^2}{2m} \nabla^2 \psi + V(\mathbf{r},t)\psi + g|\psi|^2\psi ]

但你打开NLSE.m大概率不会看到hbar、m这些符号。这不是写代码的人偷懒，而是大家都习惯先做无量纲化。光纤场景下，方程会写成：

[ i\frac{\partial u}{\partial z} - \frac{\beta_2}{2}\frac{\partial^2 u}{\partial T^2} + \gamma |u|^2 u = 0 ]

这里z是传播距离，T是随脉冲移动的本地时间，β2是群速度色散系数，γ是非线性系数。如果再把u、z、T分别用峰值功率、色散长度和非线性长度归一化，方程还能进一步变成：

[ i\frac{\partial \psi}{\partial \xi} - \frac{\text{sgn}(\beta_2)}{2}\frac{\partial^2 \psi}{\partial \tau^2} + \sigma|\psi|^2\psi = 0 ]

σ = ±1决定自聚焦还是自散焦。绝大多数孤子模拟取σ = +1、sgn(β2) = -1，对应光纤的反常色散区。你在动手改代码之前，先确认它用的是什么归一化约定，这一步错了后面全错。

2.2 分步傅里叶法：非线性是“增味”，色散才是“主菜”

NLSE 没有通用的解析解，数值上最主流的做法就是分步傅里叶法（Split-Step Fourier Method, SSFM）。核心思路听起来像玄学，其实很简单：在一个足够小的步长h内，把色散和非线性分开处理，先在时域只算非线性，再进频域只算色散，最后再回到时域。

对称分步格式是实际项目里最常用的，一个步长内走“半非线性步 → 全线性步 → 半非线性步”：

function psi = ssft_nlse(psi0, dz, Nsteps, opts) % 对称分步傅里叶法主循环 % psi0 : 初始波函数（复数数组，长度 Nx） % dz : 空间步长（传播距离步长） % Nsteps : 总步数 % opts : 结构体，含 beta2、gamma、dt_ 时间轴等信息 % 频率轴构造：这是最常见翻车点，后面详细讲 Nx = length(psi0); df = 1 / (Nx * opts.dt_); % 频率分辨率 f = (-Nx/2 : Nx/2-1) * df; % 以 0 为中心的频率轴 omega = 2 * pi * f; % 线性算子的频域表示：色散项 L = -1i * 0.5 * opts.beta2 * omega.^2; % 非线性算子：时域的相位旋转，|psi|^2 是光强 NL = @(x) 1i * opts.gamma * abs(x).^2; psi = psi0; for k = 1:Nsteps psi = psi .* exp(NL(psi) * dz/2); % 前半非线性步 psi = ifft(fft(psi) .* exp(L * dz)); % 频域内走完线性步 psi = psi .* exp(NL(psi) * dz/2); % 后半非线性步 end end

这段代码里有三个点你需要盯住：

第一，omega = 2*pi*f必须配合fftshift使用。MATLAB 的fft输出从 0 频率开始，如果直接用(-Nx/2:Nx/2-1)这种负正对称的频率轴，就必须先fftshift(psi)再fft，否则高频和低频是反的，波形会直接裂开。

第二，NL函数里的exp(1i * gamma * abs(psi).^2 * dz/2)对应的是自相位调制。这一步在时域做局部相位旋转，物理含义是光强改变了介质折射率，折射率变化又反馈给相位。

第三，L是纯含i的指数项，所以线性步在频域做乘法不改变能量，数值上天然稳定。分步法的稳定性主要被非线性步的局部误差限制。

2.3 初始条件和边界条件：高斯脉冲是最低配

调用上面的函数之前，你得先造一个初始波函数。最常用的初值是双曲正割孤子解，因为它对应理想孤子的精确形状：

% 时间轴 T0 = 1; % 脉冲宽度归一化 nt = 1024; dt_ = 0.02; t = (-nt/2 : nt/2-1) * dt_; % 一阶孤子初始场：sech 形，N=1 对应基孤子 N = 1; psi0 = N * sech(t / T0); % 调用分步傅里叶 opts.beta2 = -1; % 反常色散，聚焦型 opts.gamma = 1; % 自相位调制系数 opts.dt_ = dt_; psi_final = ssft_nlse(psi0, 0.001, 2000, opts);

双曲正割是孤子研究的默认起手式，因为它的波形在演化过程中能保持形状不变。N在这里不是点数，而是孤子阶数：N=1是基孤子，N=2及以上的高阶孤子会经历周期性压缩和分裂。如果你拿高斯脉冲当初始条件，除非功率正好卡在某个阈值，否则它会发散成色散波，而不是孤子。

边界条件方面，NLSE.m这类程序默认用周期边界，因为 FFT 天然假设信号在两端是循环的。波包一旦演化到边界，会从另一头钻进来，模拟结果就废了。实际工程里你会看到两种处理：一是把时间窗开得足够大，让波包在仿真结束前碰不到边界；二是在边界附近加吸收层。先检查代码里有没有dt_ * nt这个总窗口，如果窗口小于脉冲宽度的 20 倍，建议你直接把它改大。

3. 光学与 BEC 场景：参数映射和初始场该怎么设

3.1 光纤场景：从物理单位到归一化单位

很多人卡在参数换算。你拿到手的物理参数是β2、γ、峰值功率P0、脉冲半宽T0，但代码跑的是归一化坐标。换算关系是一组固定公式：

这个表值得你贴在工位旁边。**N² 是 L_D 和 L_NL 的比值，它决定了孤子阶数。**N=1 意味着色散和非线性恰好抵消，脉冲形状不变；N>1 时脉冲先压缩再分裂，重现周期是 π/2 倍的色散长度。

举个例子：标准单模光纤在 1550 nm 处β2 = -20 ps²/km，γ = 1.3 W⁻¹km⁻¹。如果你发一个T0 = 5 ps、峰值功率P0 = 10 mW的脉冲，那么：

[ N^2 = \frac{1.3 \times 0.01 \times 25}{20} = 0.01625 ]

N 只有 0.13，远小于 1，这只是一个弱的非线性修正，跑不出孤子。要让 N=1，P0 需要大约 30 W 量级。这件事让我在实验室吃过亏：用 10 mW 的皮秒脉冲跑 NLSE 仿真，怎么跑都稳定，还以为代码写错了，其实是物理参数根本没进入孤子区。

3.2 BEC 场景：初始波函数换成 Thomas-Fermi

NLSE.m如果用sech初始化，那基本是按光学场景写的。如果你要把它改到玻色-爱因斯坦凝聚（BEC）里研究涡旋、振荡这类问题，得换初值。BEC 的均值场模型是 Gross-Pitaevskii 方程，形式和 NLSE 相同，但势阱项V(r)不是零，而是一个谐振子势：

% 二维 BEC 初始场：Thomas-Fermi 近似 [Nx, Ny] = deal(256, 256); dx = 0.1; x = (-Nx/2 : Nx/2-1) * dx; [X, Y] = meshgrid(x, x); mu = 10; % 化学势 g = 100; % 原子间相互作用强度 % 谐振子势 V = 0.5 * (x^2 + y^2) V = 0.5 * (X.^2 + Y.^2); dens = max((mu - V) / g, 0); % Thomas-Fermi 密度分布 psi0 = sqrt(dens);

这里max((mu - V)/g, 0)是关键一步：密度在势阱边缘必须截断为 0，不能有负密度。然后你要把NLSE.m的主循环里加上势场项，线性算子变成：

[ L = -i\left( -\frac{1}{2}\nabla^2 + V \right) ]

常见做法是把∇²在频域乘上-k²，V在时域相乘，线性步变成两步：先频域色散，再时域势场。这已经超出原始代码的范围，但改造方向很固定。我的经验是先在 1D 上调通了再升 2D，直接上二维你会被 FFT 的坐标轴顺序折磨到怀疑人生。

3.3 网格和步长：先看这三件事，再跑循环

我拿到任何一份 NLSE 代码，第一件事不是看主循环，而是检查网格设置。三件事：

时间窗口T_total = nt * dt_。它必须远大于脉冲宽度的 10 倍以上，否则周期边界会在演化中引入假信号。读代码时找dt_和nt，乘法算一下就知道有没有问题。

空间步长dz必须满足dz << min(L_D, L_NL)。一般取dz < L_D / 100，不然非线性步的误差累积会让孤子形状在长距离演化后漂移。你可以做一个简单测试：把dz减半重跑，看结果变化是否在容忍范围内。如果变化明显，说明原步长不够密。

能量守恒测试：每跑完一轮，计算sum(abs(psi).^2) * dt_，这个量在理想 NLSE 下是常数。如果能量漂移超过 1%，先怀疑dz，再怀疑边界。

% 能量守恒检查：理想 NLSE 的 L2 范数恒定 Nt = size(psi_store, 2); power = zeros(Nt, 1); for k = 1:Nt power(k) = sum(abs(psi_store(:,k)).^2, 1) * opts.dt_; end plot(power)

功率曲线如果是平的，数值格式基本健康；如果缓慢上升或末尾跳变，大概率是边界问题或者fftshift没配对。

4. NLSE.m 常见翻车现场：避坑与排查手册

4.1 波形高频毛刺：fftshift 配对错误

现象：初始是光滑的sech脉冲，演化几步后波形上出现密密麻麻的高频振荡，频率轴对称的位置出现镜像分量。

原因：构造频域算子时用了对称频率轴，但送入fft前没有fftshift，或者用了ifftshift混用导致频谱错位半格。MATLAB 里fft的输出索引从 1 对应频率 0，正频率在左，负频率在右；而(-N/2:N/2-1)的数组是从负频率开始的。两者不配对，指数项exp(L*dz)就把相位加到了错误的频率分量上。

解决：统一用下面这个模板，走过的路再也不会错：

f = (-Nx/2 : Nx/2-1) * df; % 1D 频率轴 psi_hat = fftshift(fft(psi)); % 先 fft 再 shift psi_out = ifft(ifftshift(psi_hat .* lin_op));

关键记忆点：时域数组是正常顺序进fft，频域结果一定先fftshift再乘算子，最后ifftshift回正常顺序再ifft。这四步少一步都会出怪图。

4.2 长期演化波形分裂：步长 dz 太大

现象：短距离演化正常，传了 20 个色散长度之后，孤子中心出现塌陷，两翼散开，能量守恒仍然好，但波形不对称。

原因：对称分步法的局部误差是O(dz³)，步长不够密的时候，非线性步和高阶色散之间的耦合误差会累积，表现为孤子频谱的非对称畸变。能量守不守恒不能说明格式精确，只能说明总能量没丢。

解决：把dz至少减小一个数量级重跑。另外注意，有些代码里的dz是归一化坐标，不是物理距离，改的时候看清楚变量名。我的习惯是记录下dz * Nsteps算出来的总距离，保证它等于目标色散长度的整数倍。

4.3 波包消失或出现额外脉冲：边界反射

现象：演化一定距离后，预期单个脉冲却在两侧多出两个小包，随时间往中间走，像边界在“吐东西”。

原因：周期边界下，波包主体还没碰到边界，但波包的非线性尾翼已经越过边界，从另一端衍射回来，产生干涉条纹。这种问题在窄时间窗、低功率尾翼明显的场景下特别隐蔽，因为总能量变化很小，只有演化到中后期才显现。

解决：有两个土办法。第一，扩大时间窗，T_total至少到脉冲宽度的 50 倍；第二，加吸收层，在时域数组两端各 10% 的采样点乘一个余弦衰减窗，每步或每几步乘一次。我一般用第二个方案，因为它能同时压掉高频数值噪声：

% 吸收窗：两端各 10% 采样点余弦衰减 absorb = ones(nt, 1); w = round(nt * 0.1); absorb(1:w) = cos(linspace(0, pi/2, w)).^2; absorb(end-w+1:end) = cos(linspace(0, pi/2, w)).^2; % 每 10 步乘一次，避免吸收窗自身引入反射 if mod(step, 10) == 0 psi = psi .* absorb; end

4.4 用 N=1 的初值去验证 N=2 的行为

现象：跑出来的波形在第一个压缩点没有达到预期峰值，反而逐渐展宽。

原因：孤子的初始sech振幅要乘上阶数 N。很多人记着 N=1 的归一化方程是u = sech(t)，就忘了 N=2 时应该是u = 2*sech(t)，振幅缺了 2 倍，结果跑出的是 0.5 阶的非孤子。

解决：写一行断言在初始化后面：

% 验证孤子阶数是否满足 N^2 = gamma * P0 * T0^2 / |beta2| N_order = sqrt(opts.gamma * abs(psi0(1))^2 * T0^2 / abs(opts.beta2)); assert(abs(N_order - N) < 1e-6, '初始振幅与孤子阶数不匹配');

注意取值位置，psi0(1)对应时间中心，取它的模就是峰值振幅。

5. 让 NLSE.m 真正可复现：守恒量、判据和脚本习惯

5.1 三个守恒量：功率、动量、哈密顿量

数值代码是否正确，判断依据不是图好不好看，而是守恒量是否保持。NLSE 在均匀系统里有三个基本守恒量，你的代码必须具备至少前两个的监测能力：

% 计算三个守恒量 P = sum(abs(psi).^2) * opts.dt_; % 功率/粒子数 Q = 1i * sum(conj(psi) .* gradient(psi)) * opts.dt_; % 动量 H = 0.5 * sum(abs(gradient(psi)).^2) * opts.dt_ ... - 0.5 * opts.gamma * sum(abs(psi).^4) * opts.dt_; % 哈密顿量

功率守恒是分步傅里叶法自动满足的，因为每一步都是乘上纯相位因子。动量和哈密顿量不是自动守恒的，需要你在步长足够小时才能近似保持。如果你看到 H 单调上升或下降超过 0.1%，那就是步长不够。这个检查做完，你的结果才有资格写进论文或交付给项目。

5.2 孤子判据：形状因子和面积定理

除了守恒量，你还要验证“它确实是个孤子，而不是恰好没发散”。最直观的判据是形状因子，定义为主瓣宽度和峰值强度的乘积。基孤子的这个积在演化中应该不变：

% 计算脉冲宽度（强度二阶矩）和峰值 I = abs(psi).^2; w = sqrt(sum((t').^2 .* I) / sum(I) / 2); % 脉宽 p = sqrt(2) * max(abs(psi)) * sqrt(w); % 面积因子 % 理想基孤子的面积因子约 1.414

这个 1.414 的来历是sech函数的面积除以峰值再乘半高宽，实际中 1.4 到 1.45 之间都算合格。如果面积因子随演化明显偏离，你可以先怀疑初值阶数不对，再看步长。

5.3 代码习惯：把 NLSE.m 拆成三个文件

独立工作的习惯比代码本身重要。NLSE.m如果只有一个大脚本，你会陷入“改初值要翻半天”的麻烦。我建议拿到手先重构：

• nlse_params.m：集中放所有物理参数和归一化换算
• nlse_split_step.m：只做分步傅里叶演化
• nlse_diagnostics.m：画图、算守恒量、算面积因子

这样改完之后，换场景只是改params一个文件，主循环永远不动。曾经有个项目，我因为乱改主循环导致 2 周的数据全部作废，从那以后我立了规矩：仿真代码的主循环只允许优化，不允许往里塞业务逻辑。

6. 进阶验证：双孤子碰撞的调试技巧

分步傅里叶法的代码练手到一定程度，就该跑双孤子碰撞了。这个实验能一次性检验你的代码改对没有。

% 两个反向运动的孤子：一个带正频移，一个带负频移 psi0 = sech((t - 5) / T0) .* exp(1i * 2 * vin * t) ... + sech((t + 5) / T0) .* exp(-1i * 2 * vin * t);

vin是速度参数，决定两个孤子相向运动的速度。跑完整个碰撞过程，你观察两件事：碰撞后有没有产生色散辐射，碰撞后孤子振幅和形状是否复原。理想弹性碰撞下，两个 N=1 孤子碰撞后应该保持原有形状和速度，只会出现相位移。

碰撞实验里最容易发现的隐藏 bug 是时间窗不够大导致的相位误差。两个孤子反向运动，总距离是他们速度之和乘以演化时间，如果窗口只留 20 倍脉宽，孤子还没碰撞就先被周期边界“分身”了。我自己的习惯是把窗口开到 100 倍脉宽，跑完之后检查两翼是否平坦，如果有干涉条纹，直接加大窗口重跑。

验证的收尾动作是扫参数。把vin从 0.5 扫到 3，记录碰撞后的振幅

变化，画一条曲线。这条曲线应该在不同速度下都趋向水平，说明你的代码在参数空间里是稳定的。如果曲线在某个速度附近跳变，先回去检查有没有吸收窗乘的位置不对。

我还遇到过更隐蔽的问题：非线性系数gamma太大时，分步法在高强度碰撞瞬间的局部误差暴增，碰撞后孤子带上高频尾翼。这时候再减小dz效果有限，更有效的方案是换成四阶 Runge-Kutta 在非线性子步上做积分，但这属于代码重构了，不建议新手一上来就碰。

如果你手上这份NLSE.m是带注释的版本，直接先找注释里写的归一化约定；如果是裸代码，先把第 2 小节的频率轴配对检查做完，再跑第 5 小节的守恒量监测脚本。两轮检查都通过，就可以放心改参数跑自己的场景了。从那以后，我每次拿到来历不明的 NLSE 代码，都强制先跑一遍 N=1 基孤子作为冒烟测试：形状不变、能量守恒、面积因子约 1.414，三个条件全过才接着用。这个习惯帮我过滤掉了至少三个有隐藏 bug 的源码，希望帮到你。

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