非线性薛定谔方程求解代码：分步傅里叶法实现孤子演化仿真

简介这份资源提供非线性薛定谔方程的数值求解代码面向从事光纤通信、等离子体物理、流体力学及非线性光学等方向的研究生、科研人员与工程技术人员帮助解决复杂非线性偏微分方程难以手工推导与快速验证的问题。压缩包共2个文件包含1个m脚本文件与1个fig图形文件整体约15KB其中m文件承载方程离散、迭代求解与结果输出的核心逻辑fig文件则保存了配套的图形界面或结果可视化窗口便于直接运行观察演化过程。目前已有1980人学习下载说明其在相关课程设计与科研入门中具有一定参考价值。读者可借助该代码理解分步傅里叶法或有限差分法的实现思路快速复现孤子演化、波形传播等典型现象并在此基础上修改参数、更换初始条件或扩展边界处理从而节省从零搭建求解框架的时间适合作为学习非线性薛定谔方程数值解法的实践起点。1. 非线性薛定谔方程求解代码从分步傅里叶法到孤子演化的完整复现做光纤通信或者超短脉冲激光仿真的同行大概率都遇到过同一个场景想验证一个孤子传输方案或者看一眼自聚焦效应对脉冲波形的影响结果卡在数值求解这一步。非线性薛定谔方程NLSE不像线性方程那样有现成的解析解绝大多数实际参数下只能靠数值方法推进。这份求解代码就是冲着这个痛点来的——它把分步傅里叶法SSFM的完整流程封装成可直接运行的脚本覆盖从线性色散步到非线性相位旋转的核心环节适合做光脉冲传播、玻色-爱因斯坦凝聚体动力学、以及水波包演化的从业者直接拿来改参数跑结果。新手能照着跑通第一个孤子案例熟手能顺着模块拆出自己需要的边界条件和初始条件。2. 分步傅里叶法为什么是首选从算子分裂到误差阶数2.1 算子分裂的数学骨架NLSE 的标准形式可以写成i ∂A/∂z - (β₂/2) ∂²A/∂T² γ |A|² A右边两项分别对应色散算子和非线性算子。直接对整条方程做数值离散时间和空间步长必须同时压得很小计算量会迅速膨胀。分步傅里叶法的思路是把这两个算子拆开先让色散作用一小步再让非线性作用一小步交替推进。每一步都在频域或时域里做对角化运算避开了耦合项带来的矩阵求逆。具体来说色散步在频域里就是一个相位因子乘法Ã(ω, zh) Ã(ω, z) · exp(i β₂ ω² h / 2)非线性步在时域里同样是一个相位旋转A(T, zh) A(T, z) · exp(i γ |A|² h)两个步骤各自精确合起来引入的误差来自算子不对易。对称分步傅里叶法把非线性步放在两个半步色散之间局部误差降到 O(h³)全局误差 O(h²)。这就是为什么大多数开源实现默认用对称格式而不是简单的交替推进。2.2 为什么不用有限差分直接硬解有限差分法当然能解 NLSE但有两个现实问题。第一色散项的二阶导数在差分格式下要求网格足够密否则数值色散会污染真实色散第二非线性项在时域里是逐点乘法差分法处理起来反而绕远路。分步傅里叶法天然适配伪谱思路空间导数在频域里变成乘法精度随网格数指数收敛。对于孤子这类需要长距离演化、波形对相位误差敏感的场景SSFM 的性价比明显更高。常见做法是如果问题里出现陡峭梯度或者强耗散项才考虑有限差分或谱方法混合纯 NLSE 的保守演化SSFM 基本是默认选项。2.3 代码包里的模块划分拿到这份求解代码后先别急着改参数。花五分钟把目录结构过一遍后面调参会省很多时间。典型的结构是这样的文件/模块作用需要动的频率main.py入口定义网格、初始条件、调用求解器每次换问题都动ssfm_solver.py对称分步傅里叶核心推进基本不动initial_conditions.py孤子、高斯、超高斯等初始波形按需扩展analysis.py演化图、频谱、能量守恒检查出图时改config.yaml物理参数与数值参数分离调参主战场这种拆法的好处是物理参数β₂、γ、脉冲宽度和数值参数时间窗口、网格点数、步长分开管理换一个仿真案例时不用翻遍所有文件。我一般会先把config.yaml里的参数抄一遍确认量纲统一再跑main.py。3. 从零跑通第一个孤子案例网格、步长与初始条件3.1 网格设置的三个硬约束时间窗口和网格点数不是随便填的。它们必须同时满足三个条件第一时间窗口要覆盖脉冲的全部能量边缘要衰减到接近零。如果脉冲拖尾被截断频域里会出现振荡演化几步后就能看到虚假的旁瓣。第二网格点数决定频域分辨率。NLSE 里非线性效应会把能量往高频搬如果频域窗口不够宽四波混频产生的分量会折叠回来污染结果。第三步长 h 要同时满足色散相位和非线性相位的采样要求。经验规则是单个步长内最大非线性相移不超过 0.05 弧度色散相位因子的变化不超过 π。下面是一个可运行的参数配置片段# config.yaml 对应的 Python 字典形式方便直接嵌入脚本 import numpy as np # 物理参数 beta2 -20e-27 # 群速度色散单位 s^2/m负值对应反常色散 gamma 1.3e-3 # 非线性系数单位 1/(W·m) P0 1.0 # 峰值功率单位 W T0 1e-12 # 脉冲半宽单位 s # 数值参数 N 2048 # 网格点数取 2 的幂次便于 FFT T_window 100e-12 # 时间窗口单位 s约为 100 倍脉冲宽度 L 10.0 # 传播距离单位 m num_steps 20000 # 步数步长 h L / num_steps # 派生网格 T np.linspace(-T_window/2, T_window/2, N, endpointFalse) dt T_window / N omega 2 * np.pi * np.fft.fftfreq(N, ddt) h L / num_steps这段代码里N取 2048 是折中再小频域分辨率不够再大单次仿真时间明显上升。T_window取 100 倍T0是为了让孤子两侧的连续波背景充分衰减。num_steps给到 20000 是保守值实际跑孤子时可以先试 5000 步看波形是否稳定再决定要不要加密。3.2 初始条件的写法与量纲检查孤子初始条件对应 NLSE 的一阶孤子解def sech_initial(T, P0, T0): 一阶孤子初始包络返回功率归一化的复振幅 return np.sqrt(P0) / np.cosh(T / T0)这里有个容易翻车的地方np.sqrt(P0)还是P0取决于代码里 |A|² 代表功率还是振幅平方。如果求解器里非线性项写成gamma * np.abs(A)**2 * A那 A 的量纲是 sqrt(W)初始条件必须开根号。如果写成gamma * np.abs(A)**2直接乘在相位上而 A 已经是功率量纲那就不开。拿到代码后先搜一遍gamma出现的位置确认量纲约定再填初始条件。量纲检查还有一个土办法跑一步之后看能量积分是否守恒。保守 NLSE 下总能量∫|A|² dT应该几乎不变。如果第一步就掉了几个百分点多半是步长太大或者量纲错了。3.3 对称分步傅里叶的核心循环求解器的主循环不长但每一步的顺序不能乱def ssfm_step(A, omega, beta2, gamma, h): 对称分步傅里叶法推进一步 # 半步色散频域乘相位因子 A_hat np.fft.fft(A) A_hat * np.exp(1j * beta2 * omega**2 * h / 4) A np.fft.ifft(A_hat) # 全步非线性时域乘相位因子 A * np.exp(1j * gamma * np.abs(A)**2 * h) # 再半步色散 A_hat np.fft.fft(A) A_hat * np.exp(1j * beta2 * omega**2 * h / 4) A np.fft.ifft(A_hat) return A注意色散相位因子里的h/4而不是h/2因为对称格式把一步拆成两个半步色散每个半步对应h/2而相位因子里的系数是β₂ ω² / 2乘起来就是β₂ ω² h / 4。这个系数写错是高频翻车点症状是孤子要么快速展宽要么直接爆炸。非线性步里的h是全步因为非线性只作用一次。如果写成h/2脉冲峰值功率的演化速度会偏慢孤子周期对不上。3.4 跑通后的第一张验证图跑完main.py后至少看三样东西第一时域波形随 z 的演化图。一阶孤子在没有高阶效应时波形和频谱都应该保持不变。如果看到明显展宽或压缩先查步长和量纲。第二能量守恒曲线。把每一步的np.sum(np.abs(A)**2) * dt存下来画出来应该是一条水平线。漂移超过 1% 就要回头查。第三频谱对称性。初始 sech 脉冲的频谱也是 sech 形演化过程中如果出现明显不对称说明频域窗口不够或者步长太大引入了数值耗散。这三张图跑出来没问题才算真正跑通了第一个案例。后面改参数、加高阶效应、换初始条件都是在这个基础上做增量。4. 避坑与排查步长、边界与量纲的五个血泪教训4.1 孤子跑着跑着就爆炸现象前几百步波形正常之后峰值功率指数上升最终出现 NaN。原因步长 h 太大非线性相位在单个步长内超过 π相位旋转出现混叠。或者频域窗口不够宽高频分量折叠回来形成正反馈。解决先把num_steps翻倍跑一遍。如果爆炸推迟但没消失检查T_window是否覆盖了脉冲展宽后的范围。必要时把N也翻倍保持dt不变的同时扩大频域窗口。4.2 能量缓慢漂移现象能量曲线每千步掉 0.1%跑完全程掉了几个百分点。原因FFT 的周期性边界条件把脉冲拖尾从一端绕到另一端形成微弱的不连续。或者步长处于误差累积的敏感区间。解决把时间窗口加大到脉冲宽度的 200 倍以上让边缘真正衰减到零。如果还漂改用对称格式并适当减小步长。注意SSFM 本身不严格守恒能量小量漂移是正常的但持续单调下降说明有问题。4.3 频谱出现镜像分量现象频谱图上在-ω0附近出现一个不该有的峰和主峰关于零频对称。原因初始条件或者非线性项里用了实数信号FFT 后正负频率共轭对称。如果代码里只保留了正频率分量逆变换时会丢失信息产生镜像。解决确认整个流程用复数包络表示np.fft.fft和ifft成对出现不要手动截断频率轴。如果确实需要单边谱在最后分析时取不要在演化中间取。4.4 高阶效应加上去就翻车现象只加 β₂ 和 γ 时一切正常加上自陡峭项或拉曼项后几步就发散。原因高阶项包含时间导数在频域里对应乘以ω。如果频域窗口边缘的ω很大高阶项的相位因子会剧烈振荡步长要求比纯 NLSE 严格得多。解决加高阶效应时步长至少减半同时确认频域窗口边缘的相位变化不超过 π。常见做法是给高频分量加一个平滑窗函数抑制边缘振荡但要注意窗函数不能太窄否则会削掉真实的高频成分。4.5 换一组参数结果完全对不上文献现象同样的孤子阶数别人论文里演化一个周期后波形不变自己的代码跑出来展宽了。原因量纲约定不一致。孤子周期z0 π T0² / (2 β₂)里T0是半宽还是全宽β₂的符号约定都会让周期差一个因子。另外有些文献用A表示功率归一化振幅有些用sqrt(P0)。解决先算一遍孤子周期和文献里的 z 轴刻度对一下。如果差一个 π 或 2就是半宽全宽的问题。如果符号反了检查 β₂ 的正负号约定。这一步没有捷径只能逐项对齐。5. 进阶用法用守恒量做步长自适应与精度验证跑通基础案例之后真正决定这份代码能不能用在正经仿真里的是步长选择有没有依据。固定步长在弱非线性区域浪费算力在强非线性区域又不够用。一个实用的做法是用守恒量做自适应判据。NLSE 有两个核心守恒量能量E ∫|A|² dT和哈密顿量H ∫(β₂/2 |∂A/∂T|² - γ/2 |A|⁴) dT。在数值演化中这两个量的漂移速率直接反映步长是否合适。我一般会这样改主循环def adaptive_ssfm(A, omega, beta2, gamma, L, tol1e-6): 基于能量漂移的自适应步长推进 h L / 1000 # 初始步长 z 0.0 E0 np.sum(np.abs(A)**2) * (2*np.pi / (omega[1]-omega[0]) / len(A)) while z L: A_new ssfm_step(A, omega, beta2, gamma, h) E_new np.sum(np.abs(A_new)**2) * (2*np.pi / (omega[1]-omega[0]) / len(A)) drift abs(E_new - E0) / E0 if drift tol: h * 0.5 # 漂移超标步长减半重试 continue elif drift tol * 0.1: h * 1.2 # 漂移很小适当放大步长 h min(h, L - z) # 不要越过终点 A A_new z h E0 E_new return A这段代码的逻辑是每推进一步检查能量相对漂移。超过容差就退回重试步长减半远低于容差就放大步长提高效率。tol取 1e-6 是保守值对大多数孤子问题够用。如果问题里非线性很强可以放宽到 1e-5但要做收敛性验证。参数说明omega[1]-omega[0]是频域分辨率用来把离散求和还原成积分。len(A)是网格点数。这个能量计算方式假设了时域和频域的 Parseval 关系成立如果代码里 FFT 没有做归一化需要相应调整系数。验证自适应步长是否可靠可以跑两组不同tol的仿真对比最终波形。如果tol1e-6和tol1e-7的结果在视觉上无法区分说明当前精度足够。如果还有可见差异继续收紧容差直到结果收敛。还有一个容易被忽略的验证手段把传播方向反过来跑。NLSE 在无耗散时是时间可逆的正向跑 L 距离再反向跑 L 距离应该回到初始波形。如果回不去说明数值格式引入了不可逆误差。这个测试比单看能量曲线更严格我一般会在正式出结果前跑一遍。从那以后我每次换新参数或者加新效应都强制走一遍反向传播验证确认可逆性没问题再出图。希望帮到你。本文还有配套的精品资源点击获取