分步傅里叶法解NLS方程：Python源码与孤子模拟实践

简介分步傅里叶法求解非线性薛定谔方程NLS的源代码解析文档适合学习非线性光纤光学、光脉冲传输仿真以及数值求解方法的读者使用。文档聚焦科研与教学中的经典问题如何用分步傅里叶法模拟色散与非线性共同作用下的脉冲演化。内容覆盖 NLS 方程建模、非线性步与色散步的拆分原理、傅里叶变换实现以及 sech 脉冲和超高斯脉冲的参数设置与 chirp 影响。资源包整体小巧仅 1 个 docx 文件、13KB内含完整 MATLAB 源代码与逐段注释便于对照运行和迁移修改。已有 1369 人学习下载适合本科高年级、研究生及光通信方向初级研究者快速上手。通过该文档读者可以掌握从方程离散、频域相移因子设计到结果可视化的完整流程理解分步傅里叶法在孤子传播、非线性频谱演化等场景中的实际运用是一份兼顾理论与代码实践的轻量参考资料。1. 分步傅里叶法解NLS方程的源代码核心循环不到二十行凭什么成为孤子模拟的默认方案分步傅里叶法解NLS方程的源代码核心其实不到二十行。非线性薛定谔方程NLS是光脉冲在光纤中传播、玻色-爱因斯坦凝聚、水波调制不稳定性背后共用的主方程而分步傅里叶法Split-Step Fourier Method是求解它最常用也最容易被低估的数值方案把传播拆成线性衍射/色散和非线性相位两步前者用 FFT 在频域精确对角化后者在时域逐点相乘。这个思路让一段能稳定跑出孤子演化图的代码短到让人意外也正因为短参数怎么配、边界怎么处理、步长怎么选全藏在细节里。这篇笔记从数学拆解讲到可直接运行的 Python 源代码再列出五个踩坑记录最后给到自适应步长和 GPU 加速的升级路径。适合要拿这套代码做课题仿真、写论文、复现理论结果的从业者新手能照步骤跑通熟手能在这看到边界与取舍。2. 分步傅里叶法的数学拆解线性步、非线性步与对称分裂为什么是二阶精度2.1 NLS方程的标准形式与参数物理意义β2、γ、α 分别管什么先约定本文使用的方程形式这是所有源代码的地基。考虑沿光纤或一般色散介质传播的复包络 u(z,t)用随脉冲移动的延迟时间坐标系NLS 方程写成∂u/∂z i(β2/2)·∂²u/∂t² iγ|u|²u − (α/2)u其中 z 是传播距离t 是延迟时间。这个形式的好处是线性部分在频域里是精确的对角乘子非线性部分在时域里是逐点乘子正好对应分步傅里叶法的两个半步。参数物理含义常见取值取向说明β2二阶色散系数反常色散取负值如 −1负值为聚焦型 NLS支持孤子解γ克尔非线性系数正值自相位调制强度与α线性损耗无源仿真取 0有损耗时每公里扣功率P0初始峰值功率归一化后常取 1决定非线性长度T0脉冲宽度标度归一化后常取 1决定色散长度我在给初学者讲这段时通常强调一句话分步傅里叶法只是在解一个“线性导数算子 乘法算子”的演化方程物理内容全在系数里。β2 的正负决定你看到的是光孤子还是脉冲迅速展宽γ 和 P0 的乘积决定非线性相移的积累速度α 则老老实实管损耗。很多翻车现场不是算法错了是这三个系数的符号或量纲先错了。2.2 对称分步格式推导从一阶算子分裂到消灭 h² 误差把方程写成算子形式u_z (D N)u其中 D i(β2/2)∂²/∂t² − α/2 是线性微分算子N iγ|u|² 是一个依赖 u 的乘法算子严格说 N 随传播变化这就是 D 和 N 不对易的根本原因。形式上 u(zh) exp(h(DN))u(z)但 D 和 N 不对易没法直接算指数。最常见的做法是算子分裂。一阶格式是u(zh) ≈ exp(hD)exp(hN)u(z)先走一整步线性再走一整步非线性。用 Baker-Campbell-Hausdorff 公式展开exp(hD)exp(hN) 与 exp(h(DN)) 的差首项是 [D,N]·h²/2也就是说局部误差是 O(h²)整体误差是 O(h)。对称分步格式则是把线性步拆成两个半步行进u(zh) ≈ exp(hD/2) · exp(hN) · exp(hD/2) · u(z)BCH 展开后h² 项被对称结构抵消局部误差降到 O(h³)整体误差升到二阶 O(h²)。代价只是每个传播步多一次 FFT换来的是同样步长下精度大幅提升。这个格式在量子力学、光学、流体数值计算里都叫对称分裂或 Strang 分裂本质是同一件事。实现上exp(hD/2) 作用在频域很简单D 在傅里叶域是乘子 D̂(ω) −iβ2ω²/2 − α/2所以exp(hD/2)就是频域里逐点乘exp((−iβ2ω²/2 − α/2)·h/2)。而 exp(hN) 作用在时域是逐点乘exp(iγ|u(z)|²h)。这里有一个很多人误解的点非线性步里的 |u|² 用的是当前步起点处的场强不是步长中间的值这恰恰是对称格式把线性步拆到两边的原因——尽量让非线性相移作用在“线性演化到中间”的场上减少 D 与 N 不对易带来的误差。2.3 离散化约定FFT 下标、时间窗口与周期性边界条件把连续方程离散化之前必须先想清楚 FFT 的三个约定否则后面全是玄学报错。第一时间窗口是 2·Twin采样点数为 N时间间隔 dt 2·Twin/N。第二FFT 隐含信号在窗口上周期延拓下标 0 对应 t −Twin而不是 t 0频域下标 0 对应频率 0。第三np.fft.fftfreq 返回的频率排列是 0, 1/(N·dt), … 一直到奈奎斯特频率 ±1/(2dt)再翻转到负频。我在代码里统一这样构造下标避免后续所有错位问题import numpy as np N 2**12 # 采样点数取 2 的幂 Twin 40.0 # 时间窗口半宽 dt 2.0 * Twin / N # 时间采样间隔 t (np.arange(N) - N / 2) * dt # 时域坐标中心为 0 omega 2.0 * np.pi * np.fft.fftfreq(N, dt) # 角频率坐标这里t与omega的配对关系是时域数组第 N/2 个点对应 t 0正好对应频域数组第 0 个点零频。如果用np.linspace(-Twin, Twin, N)构造时间轴实际间隔是 2·Twin/(N−1)与 fftfreq 假定的 dt 不一致色散乘子就会系统性偏大这是后面 5.4 节要展开讲的经典坑。另一个必须接受的现实是周期性边界条件FFT 等价于把信号当成周期函数脉冲展宽碰到窗口边缘后会从另一侧“钻回来”。这不是 bug是 FFT 的固有属性所以时间窗口必须留足余量。奈奎斯特频率 f_max 1/(2·dt) 决定了能表示的光谱上限脉冲在传播中因非线性产生的新的高频成分如果超过它就会混叠回低频出现锯齿状伪结构。这两条是后面第四章参数设置和第五章排查的基础。3. 用 Python 写分步傅里叶法源代码从孤子初始条件到可运行的最小实现3.1 参数定义与初始条件一阶孤子 sech 脉冲怎么构造写代码先定基准测试。一阶孤子是分步傅里叶法最好的验证用例在 β2 −1、γ 1、P0 1 的归一化条件下初始场 u(0,t) sech(t) 应当以 sech(t)·exp(i·z/2) 的形式无畸变传播峰值功率恒为 1只是相位在转。这个解析解够简单任何数值错误都能被它逼出来。import numpy as np from numpy.fft import fft, ifft # 物理参数孤子单位制 beta2 -1.0 # 反常色散聚焦型 NLS gamma 1.0 # 非线性系数 alpha 0.0 # 线性损耗验证阶段必须为 0 # 离散网格 N 2**12 Twin 40.0 dt 2.0 * Twin / N t (np.arange(N) - N / 2) * dt omega 2.0 * np.pi * np.fft.fftfreq(N, dt) # 一阶孤子初始条件 P0 1.0 u0 np.sqrt(P0) / np.cosh(t) # sech 脉冲峰值场强 1 u u0.copy() # 传播参数 Lz 10.0 # 总传播距离 dz 0.01 # 步长 nsteps int(Lz / dz)这段代码的每个参数都能直接对上 2.1 节的物理量。sech 脉冲是孤子方程的本征解形状np.sqrt(P0)/np.cosh(t)保证了峰值功率恰好是 P0。如果你做的是普通色散展宽实验把初始条件换成高斯脉冲更贴近工程场景C 0.5 # 啁啾参数正值为前沿低频 T0 1.0 u0 np.sqrt(P0) * np.exp(-(t / T0)**2 / 2) * np.exp(1j * C * (t / T0)**2 / 2)啁啾项exp(iC(t/T0)²/2)是模拟通信系统中真实光源的常见做法它让脉冲在频域出现线性频率调制传播中先压缩后展宽。验证代码用 sech工程代码用高斯加啁啾这两套初始条件建议都备着。3.2 核心迭代循环线性半步、非线性步、线性半步的完整代码核心循环就是按 2.2 节的对称分裂格式推进整套代码短到让人怀疑它在正经工作# 线性半步算子预计算频域包含损耗 lin_half np.exp(1j * beta2 * omega**2 * dz / 4 - alpha * dz / 4) # 记录演化过程 z_record [0.0] u_record [u0.copy()] for n in range(nsteps): # 前半线性步频域乘算子IFFT 回时域 u ifft(lin_half * fft(u)) # 非线性步时域逐点乘用当前 |u|² 计算相位 u u * np.exp(1j * gamma * np.abs(u)**2 * dz) # 后半线性步再走一次频域半步行进 u ifft(lin_half * fft(u)) # 每 100 步记录一次快照 if (n 1) % 100 0: z_record.append((n 1) * dz) u_record.append(u.copy())注意lin_half指数里的系数。半步行进距离是 dz/2线性算子 D̂(ω) −iβ2ω²/2 − α/2所以整段指数是(−iβ2ω²/2 − α/2)·(dz/2) −iβ2ω²·dz/4 − α·dz/4。代码里写成1j·beta2·omega²·dz/4 − alpha·dz/4负号来自算式里 −i 与 1j 的对应β2 −1 时代入正好给正指数。这个系数是我见过最容易被改错的地方少乘或错位都会让色散强度翻倍。非线性步的相位是gamma·|u|²·dz|u|² 是当前时域场强的平方每一步都要重新计算。因为 u 是复数数组np.abs(u)**2与u·np.conj(u)完全等价用哪个都行但前者可读性更好。每一步三个 FFT 操作fft、ifft、fft1000 步就是 3000 次变换N 4096 时总共几秒验证阶段完全够用。3.3 验证与可视化能量守恒、孤子相位和光谱演化代码跑完后第一件事不是画图而是验证两个理论结果。其一是能量守恒整个对称分裂格式的线性步和非线性步都是模长为 1 的乘子理论上每一步都保模能量 ∫|u|²dt 的相对偏差应该在机器精度附近。其二是孤子相位一阶孤子在 z 处的理论相位是 z/2对比数值结果必须对 2π 取模。# 能量守恒检查 energy np.array([np.sum(np.abs(v)**2) * dt for v in u_record]) rel_err np.abs(energy - energy[0]).max() / energy[0] print(能量相对偏差:, rel_err) # 预期 1e-12 # 孤子相位与峰值功率校验 center N // 2 phi np.angle(u[center]) phi_theory Lz / 2 dphi abs(((phi - phi_theory np.pi) % (2 * np.pi)) - np.pi) print(f峰值功率: {np.abs(u[center])**2:.6f} (理论 1.0)) print(f峰值相位偏差: {dphi:.2e} (理论 {phi_theory:.4f}))这套校验能一次性抓出三类错误能量偏差大说明某一步的乘子模长不是 1最常见是丢了 1j相位偏差大说明线性算子的系数或符号不对峰值功率变了说明非线性强度没有正确抵消色散N² 不是 1。可视化用 matplotlib 就够了。把每个记录点的时域包络和频域光谱画出来再用 pcolormesh 画 z–t 二维演化图能直观看到孤子的“腰”在传播中保持宽度不变这是分步傅里叶法所有后续工程验证的基准图。import matplotlib.pyplot as plt freq np.fft.fftfreq(N, dt) fig, axes plt.subplots(1, 3, figsize(14, 3.5)) for zz, uu in zip(z_record, u_record): axes[0].plot(t, np.abs(uu)**2, labelfz{zz:.1f}) axes[1].plot(freq, np.abs(fft(uu))**2, labelfz{zz:.1f}) axes[0].set_xlim(-8, 8) axes[0].set_title(时域包络) axes[1].set_xlim(-0.5, 0.5) axes[1].set_title(光谱) # z-t 演化图 field np.array([np.abs(uu)**2 for uu in u_record]) Z, Tgrid np.meshgrid(z_record, t, indexingij) axes[2].pcolormesh(Z, Tgrid, field, shadingauto) axes[2].set_title(z-t 演化) plt.tight_layout() plt.show()我一般在验证通过后会顺手把初始条件改成 2·sech(t) 跑第二组测试二阶孤子应当以 π/2 为周期“呼吸”峰值功率周期性起伏。这个测试比一阶孤子更严格能暴露步长不足引起的呼吸周期畸变也顺便确认代码能处理高阶孤子场景。4. 分步傅里叶法源代码的参数怎么设步长、窗口、归一化三个关键决策4.1 步长 dz 的选取准则用最大非线性相移兜底步长 dz 是分步傅里叶法唯一的自由参数也是最难拍脑袋定的事。经验法则是用每一步的非线性相移兜底Δφ_NL γ·|u|²·dz 应当控制在 0.01 到 0.05 弧度以内。对一阶孤子γ·P0 1所以 dz 取 0.01 能留出七八倍余量但对高阶孤子脉冲在压缩阶段峰值功率可以涨到初始值的数倍这时候必须按演化过程中的最大峰值算而不是初始值。线性部分也有约束色散引起的相位变化是 |β2|·ω²·dz/4其中 ω 要取光谱有效支撑的边缘而不是奈奎斯特边缘。奈奎斯特边缘频率极高线性相位极大但那里光谱能量为零不影响精度真正要管住的是光谱半高全宽附近那一带的相位。我在工程上通常同时看两个指标非线性相移小于 0.02色散相移在光谱支撑边缘小于 π。两个都不超再开始扫参数。一个实用的检查方法是固定其他参数把 dz 减半重跑一遍看结果变化。如果时域包络或呼吸周期的差异小于你设定精度的 1/10当前 dz 就是合格的如果差异明显继续减半直到收敛。这个“步长减半对比”是分步傅里叶法最朴素也最可靠的自检手段。4.2 时间窗口与采样点数按展宽倍数和光谱支撑反推时间窗口 Twin 的选取取决于传播距离内的最大展宽。对色散主导的脉冲宽度约按 √(1(z/LD)²) 增长其中 LD T0²/|β2| 是色散长度。例子T0 1、β2 −1 时 LD 1z 10 处宽度大约展宽到初始的 10 倍Twin 40 勉强够留到 80 更稳妥。孤子虽然自身不展宽但色散波和辐射会把能量推向窗口边缘所以即使做孤子模拟我也不会把 Twin 压到小于 20·T0。采样点数 N 反过来由光谱支撑决定。脉冲展宽越窄、峰值功率越高非线性产生的新频率越多光谱越宽需要的奈奎斯特频率越高。我用两个公式快速估算频率分辨率为 1/(2·Twin)光谱最大频率为 1/(2·dt)。需要 dt 2·Twin/N所以 N ≥ 4·Twin·f_max。实际操作中我一般从 2^12 起步先看频域图两端的值是否已经掉到峰值功率的 1e-6 以下没掉就翻倍到 2^13 或 2^14。N 每翻倍内存和计算量翻倍但换来的是光谱混叠这个最难排查的隐性错误被提前消灭。4.3 无单位化与孤子阶数 N²让参数从玄学变成公式很多源码在开头直接写 β2 −1、γ 1、P0 1看起来像魔法数字实际是做完无单位化后的约定。标准做法是引入色散长度 LD T0²/|β2| 和非线性长度 LNL 1/(γ·P0)用 ξ z/LD、τ t/T0、U A/√P0 替换原方程得到只剩下一个参数的归一化方程i·∂U/∂ξ (1/2)·∂²U/∂τ² N²·|U|²U这个 N² γ·P0·LD γ·P0·T0²/|β2| 就是孤子阶数的平方也是整套参数设计的核心判据。N² 取值物理行为数值模拟注意点N² ≪ 1色散主导脉冲展宽步长由色散相位约束N² 1一阶孤子形状不变最佳验证用例N² 1高阶孤子周期呼吸峰值功率随压缩升高dz 要相应减小理解了 N²参数设置就从“试出来的玄学”变成“按公式算出来的工程决策”。典型流程是先定 T0 和 β2 算出 LD再根据想要的现象定 N²反推 P0 N²·|β2|/(γ·T0²)。模拟完成后用 U A/√P0 把归一化场换算回真实单位。我一般会在代码注释里写清楚每一步对应哪个归一化量防止三个月后回来看不懂自己写的参数。5. 分步傅里叶法常见问题排查五个高频翻车点与解决记录5.1 模拟后期窗口边缘出现“复活”的伪峰现象传播到中后段时域图窗口边缘出现对称的伪峰有时看起来像能量从两边“流”回中央中央演化图出现干涉条纹。原因FFT 隐含周期性边界条件脉冲在窗口内展宽或产生色散波后触及边缘绕到窗口另一侧重新进入。这不是分步傅里叶法算错了是时间窗口没能容纳整个演化区域。解决加大 Twin按 4.2 节的展宽公式预留 5 倍以上余量。对于强非线性场景还可以额外在窗口边缘加吸收层常见做法是在时域乘一个超高斯窗函数但注意这会破坏能量守恒记录能量时要在滤波之前算。我通常是先加宽窗口吸收层只在窗口已经很大但仍然出现边界效应时才上。5.2 光谱高频端不归零演化图出现锯齿现象频域图上光谱在奈奎斯特频率附近不跌到数值噪声水平而是形成高台阶时域演化图出现细密的锯齿状波纹。原因dt 太大奈奎斯特频率 1/(2·dt) 小于光谱实际占有带宽。非线性自相位调制不断产生新的频率成分超出采样能力后被混叠回低频污染真实光谱。解决增大 N 或减小 Twin减小 dt。判断标准是看光谱两端的功率密度是否低于峰值 1e-6 量级。我习惯初始跑一遍就看np.abs(fft(u0))**2的边界值如果初始光谱已经触边后续必然混叠直接翻倍 N 再开始。5.3 明明给了 sech 脉冲孤子却传播几下就散开现象初始条件是标准的 sech 脉冲参数看着也对但传播几个单位距离后脉冲明显展宽完全不孤子。隔壁同学用同一套代码却能跑出漂亮的孤子。原因绝大部分是归一化没对上。你心里想的是 γ·P0 1 的聚焦 NLS实际参数组合对应的 N² 远小于 1色散主导自然展宽或者 β2 取了正号正常色散聚焦型孤子根本不存在。另一个常见原因是 dz 过大非线性相移误差累积后破坏了孤子条件。解决先用三件套 β2 −1、γ 1、P0 1 跑 sech 验证确认能出孤子后再改参数。每改一次参数用 4.3 节的公式算一次 N²确保它在 1 附近。要检查 dz 是否满足 γ·|u|²_max·dz 0.02如果不满足一起改。5.4 色散强度与理论对不上差一个微小系数现象脉冲展宽速度比理论快一点或慢一点孤子相位验证总是差一个固定比例怎么调 dz 都不消失。原因时间轴与频率轴的采样间隔不一致。典型错误是用np.linspace(-Twin, Twin, N)构造时间轴实际间隔是 2·Twin/(N−1)却给 fftfreq 传了 dt 2·Twin/N导致频率轴系统性偏大或偏小一个因子N−1)/N。这个因子在 N 较小时非常明显参数扫描时表现为色散长度始终对不上。解决统一用(np.arange(N) - N/2) * dt构造时间轴dt 严格等于 2·Twin/N。如果坚持用 linspace必须把 dt 改成 2·Twin/(N−1) 再传给 fftfreq。代码里加一行打印omega[1]与 2π/(N·dt) 是否相等能在一分钟内揪出这类错位。5.5 能量随传播距离持续上涨或下跌现象设定 α 0但记录的 ∫|u|²dt 随 z 单调变化有时涨有时跌幅度远超机器精度。原因非线性步里丢了虚数单位 1j把exp(iγ|u|²dz)写成了exp(γ|u|²dz)能量按指数涨或跌或者线性损耗的符号写反−α·dz/4写成了α·dz/4损耗变成了增益或者中间插入了模长小于 1 的滤波窗而没有补偿能量。还有人是把 u 声明成了 float32累积误差在长距离传播后放大。解决先把 α 设成 0用 3.3 节的能量检查确认相对偏差在 1e-12 量级。如果能量不再守恒逐行检查乘子的模长线性乘子abs(lin_half)必须全为 1非线性乘子模长也恒为 1。确认是浮点精度问题时整个数组改为 float64这是默认值不需要额外操作。这一步排查熟练后五分钟内就能定位新手往往在这里耗一下午。6. 让源代码更耐用的三个进阶习惯守恒校验、自适应步长与 GPU 加速6.1 守恒量校验能量与光子数必须盯住基础验证之外我建议把守恒量校验写成一个独立函数每次跑新参数前先过一遍。NLS 在 α 0 时除了能量 ∫|u|²dt还有光子数频域模平方和守恒。两者同时守恒能确认频域和时域两个方向的数值操作都没做手脚。这个习惯在加滤波窗、做分步输出、改格式时特别有用能第一时间暴露非物理的增益或损耗。def check_conservation(u, dt): energy np.sum(np.abs(u)**2) * dt photon np.sum(np.abs(fft(u))**2) / u.size return energy, photon6.2 自适应步长半步比较法自动调整 dz固定步长在参数扫描时要么保守浪费算力要么激进导致误差。常见做法是步长减半比较把大步 dz 的结果和两个 dz/2 的结果对比用最大相对误差估计局部误差误差超阈值就折半步长误差远小于阈值就翻倍推进时用半步结果更准。def take_step(u, h, beta2, gamma, alpha, omega): lin np.exp(1j * beta2 * omega**2 * h / 4 - alpha * h / 4) v ifft(lin * fft(u)) v v * np.exp(1j * gamma * np.abs(v)**2 * h) v ifft(lin * fft(v)) return v def adaptive_step(u, h, tol1e-8): u_full take_step(u, h) u_split take_step(take_step(u, h/2), h/2) err np.max(np.abs(u_full - u_split)) / np.max(np.abs(u_split)) return u_split, err # errtol 时折半步长重算errtol/10 时下次翻倍这个自适应策略的代价是每一步做 6 个 FFT 而不是 2 个贵但换来的是对整个传播距离的误差可控。我在长距离仿真z 超过几十个色散长度和高阶孤子压缩这类峰值功率剧烈变化的场景里默认开自适应固定步长只用在验证阶段。6.3 cupy 一行切到 GPU二维问题同样套路分步傅里叶法天然适合 GPUFFT 和逐点乘法都是纯数值操作没有复杂分支。把 numpy 换成 cupyFFT 换成 cupy.fft循环体一行不用改import cupy as cp from cupy.fft import fft, ifft # 所有 np. 换成 cp.其余代码不变 # 画图前用 .get() 把数据拷回 CPU实战中注意两个边界N 小于 2^14 时 CPU 已经够快GPU 的优势不明显反而要承担显存拷贝开销涉及大量参数扫描时让 GPU 常驻、把多个 z 段串行计算收益更大。这套分裂思想同样能推广到二维 NLS比如 BEC 涡旋模拟只是把 FFT 换成二维 FFT线性算子和非线性算子的结构完全不变。我自己每写一版分步傅里叶法源码都会先跑一阶孤子验证再开始加功能这个习惯帮我挡掉了大部分参数错误。另一条教训是不信任一次性写对的参数每次改完都打印能量和相位偏差确认没问题再继续。这套代码的价值不在代码本身而在你愿意为它配一套可复现的验证流程希望帮到你。本文还有配套的精品资源点击获取