非线性薛定谔方程数值求解:用MATLAB分步傅里叶法实现光孤子模拟

非线性薛定谔方程数值求解:用MATLAB分步傅里叶法实现光孤子模拟 简介面向量子光学与凝聚态物理研究者的非线性薛定谔方程仿真脚本基于 MATLAB 编写可用于光孤子传输、玻色-爱因斯坦凝聚体动力学等场景的快速入门与复现。压缩包体积仅 1KB内含 1 个 m 文件结构精简适合直接阅读、运行或嵌入到自有研究流程中。资源定位偏重数值实验适合具备一定量子力学或光学基础的读者也可作为课程设计与科研入门参考。目前已有 398 人浏览学习说明该脚本在方程初学者与相关课题人员中有一定参考价值。脚本可能采用分步傅里叶法或有限差分法求解方程通过修改色散、非线性系数与初始波形可观察孤立子形成、碰撞及凝聚体振荡等现象也能帮助使用者更直观地把握非线性项与色散项之间的平衡机制为理解非线性量子系统提供可操作的数值实验平台。1. 从光纤里的一个波包说起做过光纤通信仿真的人都遇到过这个场景一个皮秒级光脉冲在色散作用下本该被越拉越宽但实际波形却在某个功率阈值之上保持了形状甚至在多脉冲碰撞后还能各自复原。这不是什么玄学而是非线性薛定谔方程NLSE在支配着光场的演化。NLSE 中的非线性项与色散项相互平衡时就产生了光孤子。NLSE.m 这个 MATLAB 源码的价值不在于它有多复杂而在于用分步傅里叶法将线性色散和非线性效应拆开处理几十行代码就能复现孤子传播、碰撞和凝聚体振荡等物理过程。对研究光纤通信、超快光学、BEC 动力学的工程师和研究生来说这是一份可以直接改参数、换初值、加势场的实用工具。2. 从线性到非线性NLSE 的归一化与算子拆分2.1 为什么标准薛定谔方程不够用线性薛定谔方程描述的是孤立粒子在势场中的波函数演化波包在自由空间中会持续展宽。但真实物理系统里存在粒子间的相互作用例如光纤中 Kerr 效应导致的折射率变化正比于光强BEC 中原子间的 s 波散射长度引入了平均场势。这两种场景都出现了波函数模平方对演化过程的反馈标准线性方程无法描述这种自作用效应于是方程右侧需要加上非线性项。NLSE 的一般形式为[ i\hbar \frac{\partial \psi}{\partial t} -\frac{\hbar^2}{2m}\nabla^2 \psi V(\mathbf{r},t)\psi g|\psi|^2\psi ]其中 (g) 是非线性系数(|\psi|^2) 是粒子密度或光强。在光学里时间变量 (t) 换成传播距离 (z)空间坐标换成归一化时间 (\tau)方程变成[ i\frac{\partial A}{\partial z} -\frac{\beta_2}{2}\frac{\partial^2 A}{\partial \tau^2} \gamma|A|^2A ](\beta_2) 是群速度色散参数(\gamma) 是 Kerr 非线性系数。代码包里的 NLSE.m 大概率采用的正是这种归一化形式的约定因为光学模拟几乎都用这套记号。2.2 分步傅里叶法的算子拆分逻辑NLSE 的优势是线性部分和非线性部分在数学形式上可以分离。如果将传播步长 (dz) 取得足够小可以近似认为色散和非线性在 (dz/2) 内独立作用这就是分步傅里叶法Split-Step Fourier Method, SSFM的核心。具体做法是对 (dz) 步长内的演化先走半个步长的纯色散项在频域完成再走一个完整的非线性项在时域完成最后再走半个步长的色散项。写成算子形式[ A(zdz, \tau) \approx e^{\frac{dz}{2}\hat{D}} \cdot e^{dz\cdot \hat{N}} \cdot e^{\frac{dz}{2}\hat{D}} \cdot A(z, \tau) ]其中 (\hat{D}) 是频域中的色散算子(\hat{N} i\gamma|A|^2) 是时域中的非线性算子。色散算子作用于频域因为 ( \partial^2/\partial \tau^2 ) 在频域中就是乘法因子 (-i\omega^2)无需做差分近似这也正是 SSFM 比有限差分法精度更高的原因线性部分精确解析求解仅非线性部分做局部近似。2.3 为什么用 MATLAB 而不是 C 或 PythonMATLAB 的 fft/ifft 自带 FFTW 库数值精度高矩阵化操作让频域滤波和非线性乘法都避免了显式 for 循环。对于一个 (2^{12}) 点的网格SSFM 单步演化只需一次 fft 和一次 ifft比 Python 手动写 FFT 或 C 调用 cuFFT 的工程成本低得多。NLSE.m 定位是教学验证与参数扫描MATLAB 是最合适的选择。3. NLSE.m 代码逐段拆解数据流与数值参数3.1 物理参数与离散网格设定NLSE.m 的头部通常是一段参数声明区。读取代码的第一步不是看运算逻辑而是先弄清楚物理参数到数值参数的映射关系% 物理参数 beta2 -1; % 色散系数 (ps^2/km)负值对应反常色散 gamma 1; % 非线性系数 (1/W/km) P0 1; % 峰值功率 (W) T0 1; % 脉冲宽度 (ps) % 数值参数 N 2^12; % 时间网格采样点数 Tmax 32 * T0; % 时间窗口半宽 dt 2 * Tmax / N; % 时间网格采样间隔 dz 0.002; % 传播步长 (km) Z 2 * pi / 2 * 1; % 总传播距离 Nz round(Z / dz); % 总步数参数说明beta2 -1是光纤通信里反常色散区的典型值负号保证孤子解存在。如果改成beta2 1波包会快速展宽可以对比观察正常色散区没有孤子的现象。gamma 1在放大比例后相当于普通单模光纤在 1550nm 窗口的典型非线性系数。N 2^12是为了调用 FFT 时的效率最高2 的幂次对 MATLAB 的 fft 性能最友好。dz的选取直接决定计算精度后面排错部分专门讨论。3.2 频域相位因子与分步演化色散算子在频域中是一个纯相位因子。代码中的频域数组构造方式如下omega fftshift((-N/2 : N/2 - 1) * (2 * pi / (2 * Tmax))); % 频率轴 disp_op exp(-0.5 * 1i * beta2 * omega.^2 * dz); % 半步步长的色散算子需要特别注意的是fftshift的使用构造的频率轴必须和 MATLAB fft 官方约定的一致即先fftshift频点再进入fft计算结果否则高频低频会互相混叠。核心演化循环通常写成% 初始双曲正割脉冲 A sqrt(P0) * sech(linspace(-Tmax, Tmax, N) / T0); for k 1 : Nz % 第一步半步长色散频域实现 A ifft( fft(A) .* disp_op ); % 中间步完整非线性时域实现 A A .* exp(1i * gamma * abs(A).^2 * dz); % 第二步另半步长色散 A ifft( fft(A) .* disp_op ); end这里的乘法方向体现了 SSFM 的算子顺序fft(A)将时域信号转到频域乘上disp_op完成色散作用再ifft回到时域执行非线性项。非线性项采用显式指数形式因为它本身是局域的微分方程没有导数项的耦合。3.3 为什么非线性项不放到频域去乘非线性项包含 ( |A|^2)这是时域局部量乘到频域上会变成卷积运算复杂度从 (O(N)) 变成 (O(N^2))完全得不偿失。SSFM 的设计哲学是什么域操作简单就去什么域做层间用 FFT 切换。4. 孤子模拟实战初值条件、演化规律与边界陷阱4.1 一阶孤子的判据与初值设定一阶孤子需要满足孤子面积条件( N_s^2 \gamma P_0 T_0^2 / |\beta_2|)当 (N_s 1) 时实现基阶孤子。在给定参数下当初值设为sqrt(P0) * sech(t/T0)时脉冲在传播中既不变宽也不窄化。改参数验证孤子特性P0 (1) / (gamma * T0^2) * abs(beta2); % 保证 N_s 1 N_s sqrt(gamma * P0 * T0^2 / abs(beta2)); fprintf(孤子阶数 N_s %g\n, N_s); % 期望输出N_s 1输出确认N_s 1后再跑演化观察每个位置 (z) 处的abs(A).^2峰值波动。基阶孤子应该是完美的平台起伏不超过 0.1%。4.2 高阶孤子与呼吸行为把P0放大四倍得到 (N_s 2)这就是二阶孤子。二阶孤子的演化特征是周期性的压缩—展宽—恢复振荡并伴随频谱展宽。若想观察多个孤子之间的相互作用初值可以改成两个重叠的双曲正割脉冲A sqrt(P0) * (sech((tau 3*T0)/T0) sech((tau - 3*T0)/T0)); % 两个相距 6T0 的初值脉冲这里脉冲间距要大于 (2T0) 才能清晰观察到碰撞后恢复原形的过程。间距太近会导致初始时刻两个孤子尚未分离演化行为更接近高阶孤子而不是孤子碰撞。4.3 边界反射与吸收层设置周期性边界条件是 FFT 自带的属性。如果脉冲在到达计算窗口边缘时还有显著强度边缘会折回形成虚假干扰。操作性规则时间窗口宽度至少要取脉冲宽度的 16 倍并且记录每个传播距离处窗口边缘的累计能量比例。如果边缘累计能量的对数超过 -30dB说明边界污染物已经不可忽略处理方法有两种扩展窗口Tmax 64 * T0; % 从 32*T0 扩到 64*T0或者加吸收层人为将边缘信号乘以衰减函数absorb exp(-((abs(tau) - 0.8 * Tmax) / (0.1 * Tmax)).^2); absorb(abs(tau) 0.8 * Tmax) 1; % 中央区域无衰减 A A .* absorb;吸收层本质是牺牲边缘几个格点的数据换取内部区域的干净。窗口越大可用的吸收过渡带越宽但 FFT 网格点数也随之增加。我的做法是先扩窗口再在窗口内加 10% 宽度的吸收带。4.4 SSFM 与有限差分法的实际差异有限差分法Crank-Nicolson 格式在时间步长选取上受限于数值稳定性条件而 SSFM 没有硬性稳定性约束步长只由物理精度决定通常可达到差分的 5 到 10 倍。然而 SSFM 的精度取决于 ( dz ) 取值过大时非线性项与色散项的耦合被忽略会出现能量漂移和相位误差。实际比对过两种方法在同一参数下模拟孤子碰撞SSFM 用 ( dz 0.002 ) 的 1000 步与 Crank-Nicolson 用极细 ( dz 0.0002 ) 的结果吻合良好。也就是说SSFM 在粗步长下就有相当高的精度这正是它在很多商业光子学仿真软件里成为默认算法的原因。5. 守恒量诊断法确认代码算对了拿到 NLSE.m 后不要急着换参数跑结果先做一组建模验证。NLSE 系统在无势场时存在三个守恒量粒子数或称功率、动量、能量。数值模拟是否有误、步长是否足够小在图像上不容易一眼看出来而守恒量的漂移是敏感的指标。粒子数守恒的离散计算公式为N_total sum(abs(A).^2) * dt; % 在每一传播步记录 N_total绘制相对变化曲线物理意义是sum(abs(A).^2)乘以时间格点间隔就是积分近似。如果这个值在演化 100 步后相对初始值漂移超过 1%说明dz取值过大需要缩小步长。标准检验是取三个不同dz如 0.001、0.002、0.004分别模拟相同物理长度画出脉冲峰值功率曲线看看是否有系统性偏离。如果粒子数守恒良好但初值孤子形状不稳定问题大概率出在参数不满足孤子阶数条件检查beta2与gamma符号和重力。如果粒子数随着演化单调递增多半是边界反射造成的伪能量注入此时应缩短传播距离或扩大时间窗。另一个非常实用的验证是调制不稳定性测试。在反常色散区给连续波叠加一个小的余弦扰动扰动频率选择在增益谱峰值附近演化后应观察到扰动呈现指数增长——这是孤子形成的「前奏」。实际代码如下% 连续波背景加扰动 A sqrt(P0) * (1 0.01 * cos(2 * pi * 0.1 * tau)); % 跑演化后观察时域强度包络是否出现周期结构的放大若这一现象不能被复现通常不是代码 bug而是时间窗口内扰动周期数太少一般至少要有 8 个扰动周期落在窗口里。此外注意连续波的背景功率设置要使得扰动增长率在gamma * P0量级让增长发生在传播距离 (Z) 内可见的尺度上。在 NLSE.m 中把守恒量检查封装成一个子函数在演化主循环中每 10 步调用一次即可在仿真中途发现问题不必等最终图出来再返工。这一诊断习惯适用于任何分步傅里叶法代码无论是光孤子与 BEC 模拟还是超快光纤激光器建模均可直接沿用。本文还有配套的精品资源点击获取