简介这份资源提供非线性薛定谔方程的数值求解代码面向从事光纤通信、非线性光学、等离子体物理等方向的研究生与科研人员帮助解决方程难以获得解析解、需要数值模拟演化过程的问题。压缩包共2个文件包含1个m脚本文件和1个fig图形文件整体约15KBm文件用于实现分步傅里叶等数值求解流程fig文件则保存了对应的仿真结果图形便于直接查看波形与频谱演化。目前已有1980人学习下载说明该代码在相关领域具有一定参考价值。读者可借助这份代码快速搭建非线性薛定谔方程的数值实验环境理解孤子演化、色散与非线性效应相互作用等典型现象并在此基础上修改参数开展自己的研究也可作为课程作业或论文复现的起点。1. 非线性薛定谔方程求解代码从光纤脉冲到 Bose-Einstein 凝聚一套分步傅里叶方案能走多远如果你正在做超短光脉冲在光纤中的演化、玻色-爱因斯坦凝聚体的动力学或者等离子体中的包络孤子那你大概率绕不开非线性薛定谔方程NLSE。它最让人头疼的地方在于非线性项和色散项耦合在一起绝大多数初值条件没有解析解。我第一次接触这个方程时试图用最朴素的显式差分去推结果步长稍大就发散步长小了又慢到没法用——这就是典型的“玄学调参”阶段。后来才明白NLSE 求解的核心不是把方程离散化就完事而是要针对它的哈密顿结构选对分裂格式。这篇笔记就围绕“非线性薛定谔方程求解代码”这个主题把分步傅里叶法SSFM和四阶龙格-库塔法RK4两条主流路线拆开讲给出可直接复现的 Python 代码、参数设置逻辑和踩坑记录。适合已经写过一维波动方程、想进一步处理非线性色散系统的从业者也适合做光学仿真但被数值不稳定折磨过的同学。2. 先分清你手里的 NLSE 是哪种形式三种常见写法与离散策略2.1 光纤型 NLSE 与凝聚体型 NLSE 的符号差异动手写代码之前必须先把方程形式确认清楚。NLSE 在不同领域里的写法差异很大符号一换代码里的非线性项和色散项位置就全变了。最常见的是光纤通信里的形式i ∂A/∂z (β₂/2) ∂²A/∂T² − γ|A|²A这里 A 是脉冲包络z 是传播距离T 是延迟时间β₂ 是群速度色散γ 是非线性系数。注意等号右边第二项前面的负号——它决定了孤子能否形成。而在玻色-爱因斯坦凝聚或非线性光学腔里常见的是i ∂ψ/∂t −(1/2) ∂²ψ/∂x² g|ψ|²ψ V(x)ψ这里 ψ 是宏观波函数g 是相互作用强度V(x) 是外势。两种形式在数学上可以通过变量替换互相转化但代码里如果符号搞反孤子会指数发散而不是稳定传播。我一般会在代码开头用注释把方程形式写死然后所有系数都按这个形式赋值避免中途混淆。另一个容易翻车的地方是有些文献把非线性项写成 |A|²A有些写成 −|A|²A这取决于他们定义的是自聚焦还是自散焦。自聚焦对应孤子自散焦对应暗孤子或波坍缩。如果你跑出来的脉冲幅度越跑越大先检查这个符号而不是急着调步长。2.2 分步傅里叶法为什么适合 NLSE算子分裂的物理直觉NLSE 可以写成 i ∂A/∂z (L N)A其中 L 是线性算子色散项N 是非线性算子非线性项。分步傅里叶法的核心思想是在一个小步长 h 内先只考虑线性项演化再只考虑非线性项演化两者交替进行。线性项在频域里是对角化的——色散项在傅里叶空间就是一个乘子 exp(−i β₂ ω² h / 2)所以用 FFT 一步就能算完。非线性项在实空间里也是对角化的——|A|²A 逐点相乘即可。这就是 SSFM 高效的原因它把两个在不同空间里对角化的算子分开处理每步只需要两次 FFT。但这里有个精度陷阱一阶 SSFM 的局部误差是 O(h²)全局误差 O(h)。如果你做长距离传播比如光纤里几十个色散长度一阶格式会积累明显误差。改进办法是采用对称分步傅里叶法也叫 Strang 分裂半步线性 → 一步非线性 → 半步线性。这样局部误差降到 O(h³)全局误差 O(h²)。代码里实现起来只是把线性算子拆成两半计算量几乎不变但精度提升非常明显。我一般默认用 Strang 分裂除非有特殊理由。2.3 什么时候该换 RK4强非线性与耗散项的处理SSFM 虽然快但它有个前提非线性项和线性项可以近似对易。当非线性强度很大或者方程里还有耗散项比如光纤损耗、增益饱和SSFM 的精度会下降。这时候可以考虑四阶龙格-库塔法RK4直接对时间或空间步进。RK4 的优点是通用性强能处理任意形式的非线性项和外部势缺点是每步要算四次斜率计算量是 SSFM 的四倍左右而且色散项在实空间里是二阶导数需要用差分近似边界处理也更麻烦。我的经验是如果非线性长度 L_NL 远大于色散长度 L_D或者脉冲演化以线性色散为主SSFM 是首选如果 L_NL 和 L_D 相当甚至更小或者方程里有显式的空间依赖势 V(x)RK4 更稳。另外如果要做实时交互或者嵌入到优化循环里SSFM 的 FFT 开销在中等网格下N1024 到 4096通常比 RK4 的多次函数求值更快。3. 用 Python 写一版可复现的 SSFM 求解器从网格到孤子演化3.1 网格与初始条件高斯脉冲的参数怎么设先确定计算域和网格。以光纤型 NLSE 为例假设我们模拟一个初始高斯脉冲在光纤中的演化。时间窗口 T 要覆盖脉冲宽度和可能的展宽一般取脉冲半宽 T0 的 20 到 50 倍。网格点数 N 取 2 的幂次方便 FFT常用 1024 或 2048。频率轴用 fftfreq 生成注意角频率 ω 2πf 的换算。初始条件用高斯脉冲A(0,T) sqrt(P0) exp(−T²/(2T0²))。这里 P0 是峰值功率T0 是脉冲半宽。如果要做孤子需要满足 N² γ P0 T0² / |β₂| 1即基阶孤子条件。我一般先算一下孤子阶数 N_soliton然后调整 P0 让 N_soliton 接近整数这样能观察到稳定的孤子传播。import numpy as np import matplotlib.pyplot as plt # 物理参数归一化单位避免量纲混乱 beta2 -1.0 # 群速度色散负值对应反常色散 gamma 1.0 # 非线性系数 T0 1.0 # 脉冲半宽 P0 1.0 # 峰值功率 N_soliton np.sqrt(gamma * P0 * T0**2 / abs(beta2)) print(f孤子阶数 N {N_soliton:.3f}) # 网格设置 N 2048 # 网格点数 T_window 40.0 # 时间窗口 T np.linspace(-T_window/2, T_window/2, N, endpointFalse) dT T[1] - T[0] omega 2 * np.pi * np.fft.fftfreq(N, dT) # 角频率轴 # 初始高斯脉冲 A np.sqrt(P0) * np.exp(-T**2 / (2 * T0**2)).astype(complex)这段代码里beta2取负值是为了让色散和非线性平衡形成亮孤子。如果你取正值需要把非线性项符号也改掉否则脉冲会展宽而不是压缩。N_soliton的计算用来判断初始条件是否接近孤子如果算出来是 1.0 左右那基本能看到稳定传播如果是 1.5会看到脉冲周期性压缩和分裂这是高阶孤子的特征。T_window取 40 是为了让脉冲边缘衰减到零避免 FFT 周期性边界带来的混叠。如果窗口太小脉冲会绕回来干扰自身。3.2 Strang 分裂的代码实现半步线性、一步非线性、半步线性下面实现 Strang 分裂的主循环。线性步在频域里做非线性步在实空间里做。注意非线性步的相位因子是 exp(i γ |A|² h)符号要和方程形式一致。这里 h 是步长我一般取 0.01 到 0.001 个色散长度具体看非线性强度。def ssfm_strang(A, omega, beta2, gamma, h, n_steps): Strang 分裂分步傅里叶法求解 NLSE i dA/dz (beta2/2) d²A/dT² - gamma |A|² A A A.copy() # 线性算子的一半频域乘子 linear_half np.exp(-1j * beta2 * omega**2 * h / 4) # 非线性算子的全步实空间相位 for _ in range(n_steps): # 半步线性 A np.fft.ifft(linear_half * np.fft.fft(A)) # 一步非线性 A A * np.exp(1j * gamma * np.abs(A)**2 * h) # 半步线性 A np.fft.ifft(linear_half * np.fft.fft(A)) return A # 演化参数 h 0.01 # 步长 n_steps 1000 # 总步数对应传播距离 z 10 A_final ssfm_strang(A, omega, beta2, gamma, h, n_steps) # 绘图 plt.plot(T, np.abs(A)**2, labelinitial) plt.plot(T, np.abs(A_final)**2, labelfinal) plt.xlabel(T) plt.ylabel(|A|²) plt.legend() plt.show()这里linear_half是半步线性乘子指数上的 h/4 是因为做了两次半步总共 h/2 的线性演化再加上非线性步的 h整体精度是二阶的。非线性步里np.exp(1j * gamma * np.abs(A)**2 * h)的符号是正号对应方程里的 −γ|A|²A 项——因为方程左边是 i ∂A/∂z移到右边变成 −i再乘上非线性项最终相位因子是 exp(i γ |A|² h)。如果你把符号写反孤子会立刻散开。n_steps和h的乘积是总传播距离我一般先跑 10 个色散长度看看波形是否稳定再决定要不要加长。3.3 步长与网格的收敛性验证怎么确认结果不是数值假象数值求解最怕的是把数值误差当成物理效应。验证收敛性的标准做法是把步长减半看结果变化是否在可接受范围内。如果步长减半后波形几乎不变说明已经收敛如果变化明显说明步长还不够小。另一个要检查的是网格点数把 N 从 1024 加到 2048看频谱是否出现新的高频成分。如果高频成分突然增多可能是非线性步产生的混叠需要增加网格或减小步长。我一般会跑三组参数h0.01, h0.005, h0.0025然后计算最终波形的 L2 范数差异。如果 h0.005 和 h0.0025 的差异小于 1%就认为 h0.005 足够。另外能量守恒也是一个重要指标NLSE 在无耗散时总能量 ∫|A|² dT 应该守恒。如果能量在演化中单调下降说明格式有数值耗散如果单调上升说明不稳定。Strang 分裂在无耗散 NLSE 下能量守恒得很好误差在 10⁻⁶ 量级。4. RK4 路线怎么搭强非线性场景下的另一种选择4.1 把 NLSE 写成实部虚部耦合的 ODE 系统当非线性很强或者方程里有显式势场时RK4 更合适。但 NLSE 是复方程直接对复数用 RK4 需要把实部和虚部分开。设 A u i v则方程 i ∂A/∂z L A N(A) 可以拆成∂u/∂z Im(L A N(A)) ∂v/∂z −Re(L A N(A))其中 L 是色散算子在实空间里用二阶中心差分近似∂²A/∂T² ≈ (A[i1] − 2A[i] A[i−1]) / dT²。非线性项 N(A) −γ|A|²A 直接逐点计算。这样我们就得到一个 2N 维的实 ODE 系统可以用标准 RK4 步进。def nlse_rhs(u, v, dT, beta2, gamma): 计算 du/dz 和 dv/dz A u 1j * v # 色散项二阶中心差分 d2A np.zeros_like(A) d2A[1:-1] (A[2:] - 2*A[1:-1] A[:-2]) / dT**2 # 边界周期性边界条件 d2A[0] (A[1] - 2*A[0] A[-1]) / dT**2 d2A[-1] (A[0] - 2*A[-1] A[-2]) / dT**2 L_A 0.5 * beta2 * d2A N_A -gamma * np.abs(A)**2 * A rhs L_A N_A return np.imag(rhs), -np.real(rhs) def rk4_step(u, v, h, dT, beta2, gamma): k1u, k1v nlse_rhs(u, v, dT, beta2, gamma) k2u, k2v nlse_rhs(u 0.5*h*k1u, v 0.5*h*k1v, dT, beta2, gamma) k3u, k3v nlse_rhs(u 0.5*h*k2u, v 0.5*h*k2v, dT, beta2, gamma) k4u, k4v nlse_rhs(u h*k3u, v h*k3v, dT, beta2, gamma) u_new u h/6 * (k1u 2*k2u 2*k3u k4u) v_new v h/6 * (k1v 2*k2v 2*k3v k4v) return u_new, v_new这里d2A的边界用周期性条件填充对应 FFT 的周期性假设。如果你要模拟非周期边界需要改用吸收边界或完美匹配层否则边界反射会污染结果。nlse_rhs返回的是 du/dz 和 dv/dz注意虚部实部的符号关系i ∂A/∂z RHS所以 ∂u/∂z Im(RHS)∂v/∂z −Re(RHS)。这个符号如果搞反波形会朝反方向演化。4.2 RK4 的步长选择与稳定性边界RK4 的稳定域比显式欧拉大得多但也不是无条件稳定。对于色散项二阶中心差分的特征值在虚轴上RK4 的稳定域在虚轴上的截距大约是 2.8。这意味着步长 h 需要满足 h * |β₂| / dT² 2.8 / 2 ≈ 1.4。如果 dT 很小h 就必须非常小计算量会爆炸。这就是为什么 RK4 在精细网格下反而不如 SSFM 划算。我一般先用 SSFM 跑一个粗解然后用 RK4 在相同参数下验证。如果两者差异在 5% 以内说明 SSFM 的步长已经足够如果差异很大要么是 SSFM 步长太大要么是 RK4 的边界处理有问题。另外RK4 对非线性项的刚性比较敏感如果 γ|A|² 很大显式 RK4 也会不稳定这时候需要考虑隐式格式或算子分裂。4.3 两种方法的性能对比什么时候 SSFM 更快在 N2048、总步数 1000 的设置下SSFM 每步两次 FFT每次 FFT 大约 O(N log N) ≈ 2048×11 ≈ 22528 次浮点操作两步 FFT 加逐点乘法总操作量约 10⁵ 量级。RK4 每步四次 RHS 计算每次 RHS 包含一次二阶差分约 3N 次操作和一次非线性项约 5N 次操作四次合计约 32N ≈ 65536 次操作看起来比 SSFM 少。但 RK4 的步长受稳定性限制通常比 SSFM 小一个量级所以总步数更多。实际跑下来SSFM 在中等非线性下通常比 RK4 快 3 到 5 倍。不过如果非线性项非常复杂比如包含卷积或非局域相互作用SSFM 的非线性步就不能逐点相乘必须用其他方法这时候 RK4 的通用性优势就体现出来了。选型时先看非线性项是否局域再看色散是否主导最后看精度要求。5. 避坑与排查NLSE 求解代码里最常见的五个翻车现场5.1 脉冲越跑越大非线性符号与色散符号不匹配现象初始高斯脉冲在演化中幅度指数增长最终溢出浮点数范围。原因非线性项符号和色散项符号没有匹配。对于亮孤子需要反常色散β₂ 0和自聚焦非线性−γ|A|²A 中的 γ 0。如果 β₂ 0 且 γ 0方程描述的是自散焦脉冲会展宽而不是压缩但不会发散。真正导致发散的是非线性项符号写反比如把 exp(i γ |A|² h) 写成 exp(−i γ |A|² h)这相当于把自聚焦变成了自散焦但色散符号没变两者不平衡数值上表现为高频成分被不断放大。解决检查方程形式和代码里的符号是否一致先用小步长跑几步看 |A|² 的最大值是否守恒。5.2 频谱出现镜像FFT 频率轴没做 fftshift现象演化后的频谱在零频两侧不对称或者出现明显的镜像分量。原因np.fft.fftfreq生成的频率轴是 [0, 正频率, 负频率] 的顺序而np.fft.fft的输出也是这个顺序。如果你直接画图而不做fftshift视觉上会以为频谱错了。更严重的是如果在频域里做滤波或加窗没有对齐频率轴会把正负频率搞混。解决画频谱时用np.fft.fftshift把零频移到中心同时用np.fft.fftshift(omega)对齐频率轴。计算时不需要 shift因为 FFT 本身是对称的。5.3 能量不守恒步长太大导致数值耗散现象总能量 ∫|A|² dT 在演化中缓慢下降下降幅度随步长增大而增大。原因SSFM 的 Strang 分裂虽然精度高但当步长 h 太大时线性步和非线性步的对易误差会表现为数值耗散。特别是当非线性相位变化超过 π 时分裂误差急剧增大。解决限制每步的非线性相位变化即 γ|A|² h 0.1。如果 |A|² 峰值是 1γ1那 h 应小于 0.1。我一般取 h0.01 或更小。另外可以用自适应步长根据当前 |A|² 最大值动态调整 h。5.4 边界反射周期性边界条件带来的假回波现象脉冲传播到时间窗口边缘后突然从另一侧出现干扰主脉冲。原因FFT 隐含周期性边界条件脉冲从左边出去会从右边进来。如果时间窗口不够大脉冲尾部还没衰减到零就碰到了边界。解决增大 T_window让脉冲边缘至少衰减到 10⁻⁶ 以下。对于高斯脉冲T_window 取 20T0 到 40T0 通常够用。如果脉冲在演化中展宽很多需要动态扩大窗口或者用吸收边界。吸收边界可以通过在窗口边缘加一个逐渐衰减的掩模来实现但会引入少量反射需要仔细调参。5.5 高阶孤子分裂物理还是数值现象初始孤子阶数 N2 或 3 时脉冲在传播中周期性分裂和恢复。原因这是高阶孤子的真实物理行为不是数值误差。高阶孤子由于非线性项和色散项的平衡不是严格的一对一脉冲会经历压缩、分裂、再恢复的周期。周期长度 z_period π/2 * L_D。如果你不确定是物理还是数值可以把步长减半再跑一次如果分裂图案不变就是物理的。解决不需要解决这是正确结果。但如果你只想看基阶孤子把 P0 调到让 N_soliton ≈ 1 即可。6. 进阶技巧用谱方法加自适应步长把精度再提一档如果你已经跑通了上面的 SSFM 和 RK4但发现长距离演化后波形还是有可见误差可以考虑两个进阶方向谱方法和自适应步长。谱方法的核心思想是用切比雪夫多项式或傅里叶级数展开空间导数而不是用有限差分。对于周期性边界傅里叶谱方法就是 SSFM 的线性步精度是指数级的。对于非周期边界切比雪夫谱方法可以在少量网格点上达到高阶精度但需要处理边界条件。自适应步长的实现思路是每一步计算两个不同步长的解比较它们的差异如果差异大于容差就减小步长小于容差就增大步长。常用的算法是嵌入式 RK 方法比如 Dormand-Prince 5(4) 对。对于 SSFM可以用局部误差估计比较一步 Strang 分裂和两步半步分裂的结果差异作为误差指标。def adaptive_ssfm(A, omega, beta2, gamma, h_init, tol1e-6, max_steps10000): 自适应步长的 Strang 分裂 A A.copy() h h_init z 0.0 steps 0 while steps max_steps: # 一步全步 A_full ssfm_strang(A, omega, beta2, gamma, h, 1) # 两步半步 A_half ssfm_strang(A, omega, beta2, gamma, h/2, 2) # 误差估计 err np.linalg.norm(A_full - A_half) / (np.linalg.norm(A_full) 1e-16) if err tol: A A_full z h h h * min(2.0, max(0.5, (tol / (err 1e-16))**0.5)) else: h h * max(0.3, (tol / (err 1e-16))**0.5) steps 1 if z 10.0: # 总传播距离 break return A, z, steps这段代码里err是相对误差用全步和两个半步的解的范数差来估计。如果误差小于容差tol就接受这一步并尝试增大步长否则拒绝并减小步长。步长调整因子用(tol/err)**0.5这是基于二阶精度的经验公式。max_steps防止死循环。实际跑的时候tol取 10⁻⁶ 到 10⁻⁸初始步长取 0.01通常能在几百步内完成 10 个色散长度的演化而且精度比固定步长高一个量级。我自己的习惯是先用固定步长跑一遍看能量守恒和波形是否合理再用自适应步长跑一遍对比两者差异。如果差异很小说明固定步长已经够用如果差异明显就改用自适应步长。另外自适应步长在非线性强度变化剧烈时特别有用比如脉冲压缩阶段步长会自动减小避免误差积累。希望帮到你。本文还有配套的精品资源点击获取