简介：这套基于 MATLAB 的数值模拟脚本聚焦非线性薛定谔方程（NLSE），面向量子光学、光纤通信、光孤子与冷原子物理等方向的研究者或高年级学生，用于求解含非线性项的波函数演化问题，并可进一步分析自相位调制等典型现象。压缩包仅含 1 个 m 文件，整体大小约 1KB，轻量易读，适合作为 NLSE 仿真入门模板。脚本可能基于分步傅里叶或有限差分法实现，能够模拟光孤子的形成与稳定传播、自相位调制效应，以及玻色-爱因斯坦凝聚体中的动力学过程；结合方程形式与物理背景注释，读者可快速理解非线性系数、色散与势能的作用。该资源已有 399 人学习浏览，对需要快速验证 NLSE 基本算法或开展二次开发的用户而言，是一个简洁实用的计算起点，也便于在此基础上拓展更复杂的实际场景。

1. NLSE 是什么：为什么一个 zip 包能装下非线性薛定谔方程的“标准玩法”

NLSE.zip 这个压缩包，第一次见的人容易误会：薛定谔方程不是量子力学吗？怎么还来一个非线性版本？实际上，非线性薛定谔方程（NLSE）是很多波动系统共用的数学模型，光纤里光脉冲怎么走、深水波怎么调制、玻色-爱因斯坦凝聚体里涡旋怎么动，都可以落到它头上。NLSE.zip 里装的东西，通常不是几百页公式推导，而是一套用数值方法把偏微分方程“煮”成波形演化图的脚本。它能解决的问题很具体：给定初始波形，看它沿传播方向怎么变形、分裂、碰撞。适合那些不想从零重写求解器，又想快速验证光学仿真或物理模型的人。

真正反直觉的是：方程看着难，主流解法却短得惊人。分步傅里叶法三十行就能写出来，但很多人栽在归一化和 FFT 顺序上。这类包的核心价值不是“帮你算出结果”，而是替你把这些坑填好。拿到压缩包之后，能不能信任它，仍然取决于你有没有办法验证它。

2. 从方程到可算的代码：分步傅里叶法怎么拆解 NLSE

2.1 非线性薛定谔方程的两种常见归一化写法

NLSE 之所以叫非线性薛定谔，是因为它长得很像量子力学里带势函数的薛定谔方程，只不过势被换成了波包自己的模平方。最常碰到的无量纲形式是：

i ∂u/∂z + 0.5 ∂²u/∂τ² + |u|²u = 0

这里 u 是慢变包络，z 是传播距离，τ 是随波包运动的延迟时间。这个形式干净，适合做孤子数学推导。但很多工程人员手里拿到的是光纤传输方程：

i ∂A/∂z - (β2/2) ∂²A/∂T² + γ|A|²A = 0

β2 是群速度色散系数，γ 是非线性系数。两者描述同一个物理过程，但归一化方式不同。拿到 NLSE.zip 之后，第一件事不是看主循环，而是确认里面的方程用的是哪个符号约定：是+0.5还是+1，是|u|²u还是2|u|²u。同一个孤子解在不同约定下幅值会差 √2，后面的传播对比全都会错。

分步傅里叶法（Split-Step Fourier Method）能在几十行内搞定 NLSE，靠的是把方程拆成线性部分和非线性部分。线性部分在傅里叶空间里是一个乘法，非线性部分在实空间里是逐点乘。如果直接整体积分，非线性项让问题无法用解析传播；但拆开之后，线性步用 FFT 精确算，非线性步乘一个相位因子，交替推进即可。

2.2 线性算子、非线性算子与半步走法

标准的分步法是先走半程线性步，再走一个完整非线性步，最后再走半程线性步。这样整体误差从一阶升到二阶。具体到代码，以无量纲方程i du/dz + 0.5 d²u/dt² + |u|²u = 0为例：

import numpy as np from numpy.fft import fft, ifft, fftfreq # 假设 u0 已按离散时间窗给出 N = u0.size dt = tau[1] - tau[0] # 这个频率序列与 fft(u) 的输出顺序一致 omega = 2 * np.pi * fftfreq(N, dt) # 线性半步传播子：完整线性 dz 是 exp(-0.5j * omega^2 * dz) lin_half = np.exp(-0.25j * omega**2 * dz) u = u0.copy() for _ in range(steps): u = ifft(lin_half * fft(u)) # 前半步线性 u = u * np.exp(1j * np.abs(u)**2 * dz) # 完整非线性 u = ifft(lin_half * fft(u)) # 后半步线性

这里的关键是omega直接用fftfreq，不要额外做fftshift。fftfreq返回的频率顺序恰好像fft输出的顺序一样，从直流、正频率到负频率。如果先fftshift(fftfreq())却不处理fft(u)的排列，那么传播子会作用到错误的频率分量上，这是新手最常踩的翻车点，后面会专门讲。

非线性项写成np.abs(u)**2，不能写成u**2。u 是复数，u**2在数学上是复数的平方，不是强度；只有模平方才对应波包的强度。

2.3 拿到 NLSE.zip 后先找哪几个文件

我一般不会急着看主函数，而是先找三样东西：核心求解器、初始条件生成器、验证脚本。

核心求解器通常是一个短函数，输入u0, tau, dz, total_z，返回传播后的u。初始条件生成器负责把物理单位换算成归一化的 sech 脉冲或高斯脉冲。验证脚本里应该跑一个已知解析解，比如基态孤子，来确认能量守恒和波形不变。

常见做法是解压后先建一个干净环境，把核心求解器 import 进来，跑最小算例，不碰任何绘图依赖。如果包里有run_demo这类入口，先跑通 demo 再改自己的初始条件。如果你发现压缩包里只有孤零零一个脚本，没有验证手段，我的建议是：不要直接用，先按下一章的方法自己写一个最小求解器，把它的输出和包里的结果对比，再决定是否信任。

3. 用 NLSE 包跑通第一个算例：单孤子传播的最小配置

3.1 最小 Python 实现：从解压到跑起来

这里给一个可以直接放进去跑的最小脚本。它模拟基态亮孤子u0 = sech(tau)在反常色散条件下传播，理论上波形不应该变形，只有相位整体旋转。

import numpy as np from numpy.fft import fft, ifft, fftfreq # 离散网格 N = 512 # 采样点数 T_max = 30.0 # 时间窗半宽，soliton 衰减到 0 后才截断 tau = np.linspace(-T_max, T_max, N, endpoint=False) dt = tau[1] - tau[0] # 频率轴，与 fft 输出顺序一致 omega = 2 * np.pi * fftfreq(N, dt) # 传播参数 dz = 0.005 # 单步距离 z_max = 4.0 # 总传播距离 steps = int(round(z_max / dz)) # 初始条件：基态亮孤子 u0 = 1.0 / np.cosh(tau) def propagate(u0, omega, dz, steps): u = u0.copy() lin_half = np.exp(-0.25j * omega**2 * dz) # 半步线性传播子 for _ in range(steps): u = ifft(lin_half * fft(u)) u = u * np.exp(1j * np.abs(u)**2 * dz) u = ifft(lin_half * fft(u)) return u u_end = propagate(u0, omega, dz, steps) # 能量守恒验证 Q0 = np.sum(np.abs(u0)**2) * dt Q1 = np.sum(np.abs(u_end)**2) * dt print(f"Q0={Q0:.8f}, Q1={Q1:.8f}, diff={Q1-Q0:.2e}")

这段代码里，tau用linspace(..., endpoint=False)生成，保证时间窗首尾不相邻，适合 FFT 隐含的周期性边界。omega没有 shift，直接与fft(u)相乘。lin_half对应半个 dz 的线性传播，每次循环用两次，合起来正好是一个完整 dz 的色散相位。

如果执行后diff在1e-10量级附近，说明传播子写对了。如果diff明显漂移，比如1e-3，先从 FFT 归一化和频率轴顺序查起。

3.2 从 NLSE.zip 换成自己的物理参数：单位换算怎么做

很多 NLSE 求解器只接受无量纲参数，因为分步傅里叶法本质上是在解归一化方程。你需要把物理量换算成无量纲量。常见做法是取脉冲宽度 T0 和峰值功率 P0，定义：

N² = γ P0 T0² / |β2|

这个 N 是孤子阶数。N=1 时对应基态孤子，波形稳定传播；N=2 时会出现周期性分裂和复原，这也是验证求解器是否可靠的好观察对象。

假设你的脉冲是双曲正割型，峰值功率 P0 要满足 N=1，那么初始振幅sqrt(P0)的时间宽度和色散长度必须匹配。比如物理参数 β2=-1 ps²/km、γ=1 W⁻¹/km、T0=1 ps，想得到 N=1 的孤子，峰值功率应该取 P0 = |β2| / (γ T0²) = 1 W。如果 NLSE.zip 里要求输入的是无量纲 u0，你就直接给u0 = sech(tau/T0)，并把 z 的单位折算成色散长度LD = T0² / |β2|。

这里最容易犯的错是只换功率不换时间尺度。u 的自变量 τ 必须用 T0 归一，否则波形的宽度和相位演化都不对。

3.3 画波形演化图时最该存哪几个量

验证单孤子传播，不能只看最终波形。我通常在每个 dz 间隔处保存若干帧 u 的值，然后画成强度随时间和传播距离的二维图。对基态孤子，二维图应该是一条笔直亮线；如果出现周期性呼吸，说明初始条件不是 N=1，或者纵向步长太大。

save_step = 50 record = [] for k in range(steps): u = ifft(lin_half * fft(u)) u = u * np.exp(1j * np.abs(u)**2 * dz) u = ifft(lin_half * fft(u)) if k % save_step == 0: record.append(np.abs(u)**2)

这里save_step控制存储间隔，不要每个循环都存，否则内存会爆。记录np.abs(u)**2而不是复数 u，因为功率图更直观；但如果你想回看相位演化，必须存复数。

4. 参数怎么设才不翻车：网格、步长与非线性强度的取值边界

4.1 时间窗与频域分辨率：一对必须同时满足的约束

NLSE 模拟的网格选择不是拍脑袋。时间窗半宽 T_max 必须大到让波形在边界处基本衰减到 0，否则 FFT 的周期延拓会把边界上的残余值“折返”进窗口内，形成假的干涉条纹。

频域分辨率由时间窗总宽度决定：

domega = 2π / (2 T_max)

如果 domega 太大，孤子的频谱细节采不进去，波形会出现弥散。经验法则是：频谱半高宽内至少要采 10 个点，且时间窗半宽至少是初始波形半宽的 5 倍。对于sech(tau)，半宽约 1，T_max 取 10 以上才稳。

采样点数 N 决定最大可表示频率：

omega_max = π / dt

dt 越小，最高频率越高，计算越贵。N 常见取 512 或 1024。对于孤子模拟，N=512 已经够；如果你要模拟高阶孤子或者强压缩，N=2048 也不夸张。下面是常见取值表：

4.2 步长 dz：先看相位增量，再试能量守恒

分步傅里叶法虽然是二阶精度，但每个完整步里线性相位和非线性相位都要求变化缓慢。一个实用的准则是让单个 dz 内最大线性相位增量不超过 0.1 弧度：

max(0.5 * omega^2 * dz) < 0.1

对非线性步，要求max(|u|^2 * dz) < 0.1。你可以先粗选 dz，然后程序打印这两个量：

phase_linear_max = np.max(np.abs(0.5 * omega**2 * dz)) phase_nonlinear_max = np.max(np.abs(u0)**2 * dz) print(f"linear phase={phase_linear_max:.3e}, nonlinear phase={phase_nonlinear_max:.3e}")

如果某个值超过 0.1，就把 dz 减半。这里不是玄学，而是因为传播子展开exp(i φ)时，φ 超过 0.1 后泰勒展开项开始影响误差。注意 dz 减半会双倍增加计算量，所以不要盲目取极小值；先用能量守恒做锚点，再逐步加 dz 观察结果是否变化。

4.3 非线性强度和孤子阶数：决定了你做的是验证还是找物理

NLSE 的非线性强度通常用孤子阶数 N 描述。N=1 是基态孤子，传播不变，适合验证求解器。N=2 及以上会出现呼吸或分裂，适合研究相互作用和频谱展宽。初始条件幅值乘以 N，就得到高阶孤子的起点：

u0 = N * sech(tau)

高阶孤子不能用基态的“波形不变”标准去检查，否则一定会觉得求解器坏了。正确做法是先跑 N=1 验证能量和形状，再切到 N=2 观察演化，看是否出现预期的周期。周期长度与归一化坐标有关，通常与π/2同量级。如果你发现演化周期和文献对不上，查方程系数是0.5还是1，以及 z 是否按色散长度归一。

5. 常见问题排查：NLSE 模拟中五个“看起来对、实际错”的坑

5.1 孤子宽度周期性抖动：FFT 归一化写错

现象：初始 sech 脉冲传播几个 dz 后，宽度开始周期性鼓包，压缩、展宽、再压缩，即使参数已经是 N=1。

原因：FFT 和 IFFT 的归一化系数不匹配。NumPy 的fft和ifft是互逆的，但如果你自己写了傅里叶变换循环，或者在传播子前多乘了一个1/N，每个循环都会引入微小增益，波形就被“泵”起来了。

解决：统一用 NumPy 的fft/ifft对，不要手写。可以加一条自检：

u_test = np.random.randn(N) + 1j * np.random.randn(N) assert np.allclose(ifft(fft(u_test)), u_test)

如果这条没过，说明环境或者版本里傅里叶变换被改过。先修这个，再谈孤子。

5.2 能量不守恒：模平方写成了u**2

现象：脉冲峰值慢慢变高或变矮，总功率Q = sum(|u|^2) dt随 z 单调漂移。

原因：非线性相位因子应该是exp(i |u|² dz)，却写成了exp(i u² dz)。复数 u 的平方不等于强度，它会改变相位调制幅度，能量守恒被破坏。

解决：只要出现模平方，一律用np.abs(u)**2。在 Python 里np.abs(u)**2和u.real**2 + u.imag**2等价。如果因为性能想避免临时数组，可以写u.real*u.real + u.imag*u.imag，不要图省事写u**2。

5.3 时间窗边界反射：周期边界的“回声”

现象：孤子传播一段距离后，波形主峰旁边出现一个逐渐靠近的小包，而且传播得越远越多，像回声。

原因：FFT 假设信号在窗口内周期延拓。只要窗口边界处波形没有衰减到浮点误差量级，边界上的残余功率就会在下一周期“绕回来”，叠加在真实波形上，产生干涉条纹。

解决：把 T_max 加宽到波形初始峰值功率低于1e-6的位置。如果加宽后计算量太大，可以在窗口两端加吸收层，比如每步把边界 10% 区域乘以高斯衰减，但这类处理会引入人为损耗，最好先用足够宽的时间窗规避。

5.4 频率轴顺序错乱：fftshift 到底要不要用

现象：波形整体错位，峰值不在 t=0，或者演化后出现明显的不对称展宽。

原因：频率轴顺序和fft(u)输出顺序不匹配。fftfreq(N, dt)返回的顺序本身就是和fft输出对齐的，不需要 shift。如果用了fftshift(fftfreq(...))，那么必须相应地对频谱做fftshift，算完再ifftshift回来。最安全的方法是彻底不 shift：

omega = 2 * np.pi * fftfreq(N, dt) lin = np.exp(-0.25j * omega**2 * dz) u = ifft(lin * fft(u))

这里没有fftshift，也就没有机会弄反。如果你在别人的代码里看到 shift，先搞清他是在哪个域里操作，不要照搬。

5.5 结果和文献对不上：方程系数约定不一样

现象：同一个 N=1 孤子，你的宽度比论文窄 √2 倍，或者演化周期正好差一半。

原因：NLSE 的系数归一化有多个流派。有人用i∂u/∂z + 0.5 ∂²u/∂τ² + |u|²u = 0，有人用i∂u/∂z + ∂²u/∂τ² + 2|u|²u = 0。这个差异直接改变孤子基底形状的缩放关系。

解决：拿到任何新的 NLSE 包，先做一次“基准对拍”：用一个已知解析解，比如基态孤子，分别在你的代码和包里各跑一次，对比幅值、宽度、相位演化。只要对拍通过，形式约定就不重要了。对拍不过，就去查包里的方程系数。血泪经验：不要直接拿论文图形的坐标来比，先看论文的方程再谈对比。

6. 把 NLSE.zip 改造成自己的工具：验证技巧与迭代习惯

前面五章讲的是把求解器跑通，最后这一步做的是“接手别人的包”时最值得养成的习惯。

拿到 NLSE.zip，我通常先建一个tests/目录，把基态孤子能量守恒、宽度不变、边界功率三条写成固定测试。每次改参数或者换初始条件，先跑测试再改主代码。这样做的好处是，当你把包里的 demo 换成自己的物理参数后，如果结果开始奇怪，你能分清是物理模型变复杂了，还是求解器被改坏了。

一个简单的测试函数可以这样写：

def test_mono_soliton(): u0 = 1.0 / np.cosh(tau) u_end = propagate(u0, omega, dz, 100) Q0 = np.sum(np.abs(u0)**2) * dt Q1 = np.sum(np.abs(u_end)**2) * dt assert np.isclose(Q0, Q1, rtol=1e-8) assert np.max(np.abs(np.abs(u_end) - np.abs(u0))) < 1e-4

第二个断言对基态孤子成立，因为传播理论上只有相位变化。如果你的研究方向是双孤子碰撞或自陡峭效应，可以把测试换成双脉冲碰撞前后的光谱守恒，思路一样。

另一个习惯是：每次跑完保存三样东西——初始条件、参数表、最终波形。不要只存一张演化图。很多问题要回去对参数才能定位，没有原始记录，等于没有做过这次仿真。

我早期踩过最大的坑，就是拿到一个现成脚本后直接改参数，结果波形怎么都不对。折腾了三天，最后发现是别人代码里的“时间窗”和“空间步长”用的是不同单位制。现在我会先在笔记本上写一行字：“这个包是什么归一化约定？哪个变量是传播方向？哪个变量是时间？”想清楚这三个问题，NLSE.zip 就变成了一个真正顺手的工具，而不是黑匣子。希望帮到你。

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