简介围绕非线性光学仿真与计算的代码资源包面向光学工程、物理专业学生及科研人员用于研究强光照射下材料的非线性光学响应如二次谐波产生SHG、参量下转换SPDC、四波混频FWM等效应并通过仿真分析光场传播、相位匹配与能量转换涵盖非线性材料与典型应用场景。资源共577个文件以MATLAB脚本、C源代码和MEX编译模块为主体辅以Visual Studio工程文件、CHM/HTML帮助文档、图形与实验数据文件整个zip压缩包仅5.23MB。内容包含可运行的仿真程序、核心算法实现、MEX接口库、示例数据及说明文档能帮助学习者从非线性极化率、波动方程求解到相位匹配计算逐步深入适合课程设计、课题预研或自主复现非线性光学仿真实验也可用于光学频率转换、量子信息等方向的前期探索。目前已有853人学习下载兼具代码参考与实践指导价值。1. 非线性光学仿真到底在仿什么一份代码包能帮你省掉哪些事第一次拿到 Nonlinear-Optics-master 这个名字的人多数是冲着“非线性光学”四个字来的。但很多人第一次跑出来的图不是漂亮的孤子演化也不是平滑的超连续谱而是一团噪声或者功率稍微调大一点就输出 NaN。这不是代码损坏而是单位、步长、采样三个环节的连环坑。这个标题背后是一套用数值方法求解非线性薛定谔方程和耦合波方程的仿真框架它能帮你算脉冲在介质里的自相位调制、自聚焦、倍频、四波混频也能帮你做光纤激光器和超连续谱光源的前期方案验证。适合刚接手激光器仿真任务的人、做光纤通信研究的在校生、想用仿真验证实验方案的工程师。这篇笔记按“先选方程、再跑通、后调参”的顺序来拆最后落在你一定会用到的验证技巧上。2. 从耦合波方程到广义非线性薛定谔方程仿真引擎的底层逻辑2.1 先选方程二阶过程用耦合波三阶过程用 GNLSE非线性光学的仿真不是一套代码吃所有场景。拿到标题里的“非线性光学”四个字第一步是判断你要算的是二阶非线性还是三阶非线性这决定了你该解哪类方程。倍频、和频、差频、光参量振荡这类二阶过程本质是三个光场通过二阶极化率 χ⁽²⁾ 交换能量对应的数值模型是耦合波方程。三个振幅在传播方向上各自演化彼此通过非线性耦合系数 κ 连在一起同时还要考虑相位失配 Δk。方程形式大致是 dA₁/dz iκ A₂ A₃ e^{-iΔkz}这个 Δk 是这类仿真里最敏感的物理量稍稍偏离相位匹配倍频效率可能掉一个数量级。三阶过程则是另一套逻辑。自相位调制、交叉相位调制、四波混频、光孤子、超连续谱产生都来源于三阶极化率 χ⁽³⁾在波导和光纤场景里被统一写进广义非线性薛定谔方程GNLSE。它把色散写成一组高阶色散系数 β₂、β₃、β₄ 的泰勒展开把 Kerr 非线性写成一个与强度有关的相移项对超短脉冲还要加入拉曼延迟响应。GNLSE 里每个项都有明确的物理意义这也是为什么它成了光纤和非线性波导仿真的事实标准。选型原则我一般这样掌握连续光或宽脉冲的倍频/参量过程优先耦合波脉冲在波导里传播、要看到频谱展宽或孤子行为直接用 GNLSE。如果你发现目标论文里出现“自相位调制”“孤子阶数”“超连续谱”这些词那它十有八九跑的是 GNLSE。在动手写代码前还有一个常被忽略的中间变量有效模面积 A_eff。非线性系数 γ 2πn₂/(λ A_eff)它把材料本征非线性 n₂ 换算成实际波导里能感受到的强度依赖相移。同样一段光纤模场直径小一点γ 就大很多。这也是很多人拿到别人的代码后发现同样的泵浦功率自己仿真结果差一截的原因之一。2.2 分步傅里叶法SSF为什么是首选算法一旦确认要解 GNLSE第一个要面对的数值问题就是怎么处理传播算符。GNLSE 里线性部分色散、损耗在频域是一个简单的乘法算符非线性部分Kerr 相移、拉曼项在时域是一个局部操作。直接用有限差分法做全时域推进步长要压得非常小才能稳定计算量让人很难受。所以工程上几乎都用分步傅里叶法把传播拆成线性步和非线性步交替执行。线性步在频域完成因为色散的微分算子在对角化后就是乘一个相位因子非线性步在时域完成因为 |A|² 是当地量直接乘相位就可以了。这样每个 z 步只需要两次 FFT 和一次逆 FFT复杂度是 O(N log N)比全时域差分快几十倍。分步法也有精度档次差异。最简单的非对称版本先走整个 Δz 的线性步再走 Δz 的非线性步误差是 O(Δz²)。对称分步法则把线性步拆成前后各 Δz/2中间夹着非线性步误差提升到 O(Δz³)实现代价只是多一次 FFT非常划算。我自己的习惯是默认写对称分步除非只是做量级估算才用非对称版。数值方法线性步处理方式单步精度适用场景非对称 SSF整个 Δz 一次处理O(Δz²)快速摸底不追求精细结构对称 SSF前后各 Δz/2O(Δz³)常规仿真首选自适应 RK4IP积分因子 四阶龙格库塔更高飞秒脉冲、强拉曼增益、光孤子裂变另外要提一句“失步”这个老概念分步法的本质是把两个不对易的算符强行拆开所以哪怕是对称版也存在算子交换带来的残留误差。只有把 Δz 压到足够小让每一步里非线性相位变化远小于 1 弧度结果才会收敛。后面第 4 章会专门讲怎么定这个步长。3. 把代码包跑起来环境、文件结构和最小复现实例3.1 代码包里通常有哪些文件它们分别干什么这类以 master 命名的光学仿真代码包组织方式大多比较套路化拿到手先别急着跑主脚本先花五分钟把目录扫一遍。常见的布局是一个入口脚本负责定义物理参数并调用传播核心一个 solver 文件夹放分步傅里叶或 RK4 的求解函数一个 materials 文件夹放折射率、色散、非线性系数的查找表examples 文件夹放示例配置docs 放说明。入口脚本通常长这样先设波长、脉宽、峰值功率、传播长度再调用 solver最后画图和存数据。核心求解函数一般不区分二阶还是三阶过程它只负责“给我初始脉冲和参数我把演化结果还给你”真正决定物理的是你传入的 β、γ、拉曼响应函数。如果你发现代码里没有单独的材料参数文件参数全写在主脚本开头也正常。特别是教学性质的代码为了降低阅读门槛会把所有物理参数摊开放在一个区块里。这种结构对初学者其实更友好因为改参数只需盯一个区域。我一般拿到这类代码后做的第一件事不是运行而是把参数区块里的单位全部写进注释里并且确认主脚本最后输出了哪些量。很多包默认只输出时域波形频率谱需要自己额外算一次 FFT这个不提前看清楚容易误以为包“只给了一半结果”。运行环境的差别也比较明显。如果是 MATLAB 写的只需要确认用的是哪个版本以及有没有用到信号处理工具箱里的 fftshift 之类的函数。如果是 Python 写的常见的依赖就是 numpy、scipy、matplotlib 三件套。我的建议是建一个干净的虚拟环境再装依赖不要直接往系统 Python 里塞否则后面你会发现装某个科学计算库时把系统环境搞坏了所有项目一起翻车。3.2 最小可复现脚本一段对称分步傅里叶法跑通脉冲传播与其先研究某个具体包里的花哨封装不如先在本地把一段最核心的对称分步傅里叶法跑通。下面这段脚本不依赖任何特定仓库是这类非线性脉冲传播仿真的公共骨架也是很多代码包 solver 部分的简化版。我把高斯脉冲放进反常色散光纤观察它在一段距离里的演化并画出末端的时域波形和频谱。import numpy as np from math import factorial import matplotlib.pyplot as plt # 物理参数单位必须和后面傅里叶频率一致 beta2 -20.0 # 群速度色散ps^2/m负号为反常色散 gamma 0.02 # 非线性系数W^-1·m^-1 alpha 0.0 # 损耗m^-1 P0 1.0 # 峰值功率W T0 1.0 # 脉冲半宽度1/e 强度点ps # 数值参数 Nt 4096 # 采样点数取 2 的幂便于 FFT Twin 100.0 # 时域窗口宽度ps L 100.0 # 传播长度m dz 0.2 # 纵向步长m t np.linspace(-Twin / 2, Twin / 2, Nt, endpointFalse) dt t[1] - t[0] # 角频率单位是 rad/ps对应上面长度单位取 ps^2/m 时相位无量纲 w 2 * np.pi * np.fft.fftfreq(Nt, dt) # 初始高斯脉冲 aw np.sqrt(P0) * np.exp(-(t / T0) ** 2 / 2) # 对称分步傅里叶法推进 z np.arange(0, L dz, dz) lin_half np.exp(-alpha / 2 * dz 1j * beta2 / 2 * w**2 * dz) for zi in z[1:]: aw np.fft.ifft(np.fft.fft(aw) * lin_half) # 前半个线性步 aw aw * np.exp(1j * gamma * np.abs(aw)**2 * dz) # 非线性步 aw np.fft.ifft(np.fft.fft(aw) * lin_half) # 后半个线性步 # 输出末端时域波形与光谱 spec np.abs(np.fft.fftshift(np.fft.fft(aw)))**2 f_axis np.fft.fftshift(w) / (2 * np.pi) # 频率轴THz 量级 fig, ax plt.subplots(1, 2, figsize(10, 4)) ax[0].plot(t, np.abs(aw)**2) ax[0].set_xlabel(t (ps)); ax[0].set_ylabel(功率 (W)) ax[1].plot(f_axis, spec) ax[1].set_xlabel(f (THz)); ax[1].set_ylabel(强度) plt.tight_layout(); plt.show()逻辑说明这段代码的核心是对称分步线性算子在频域乘一次相位因子分前后两半执行非线性算子在时域乘一个与强度相关的相位。β₂ 前的负号对应反常色散是孤子能存在的前提如果把它改成正值你会看到脉冲展宽而不是压缩。代码里的 lin_half 把一半线性步提前算好因为整个传播过程中它不变这是一个很小的性能优化习惯也可以避免每次循环里重复构建 exp 数组。参数说明dz 取 0.2 m在 100 m 传播距离上共 500 步第一步跑完如果能量不守恒再缩小。Nt 与 Twin 的配合关系到频率分辨率频谱频率间隔 Δf 1/Twin 0.01 THz可以做精细光谱时域采样间隔 Δt ≈ 0.024 ps对应奈奎斯特频率约 20 THz对 1 ps 脉冲的 SPM 展宽绰绰有余。跑通这段脚本后你就有了一份可以对照任何代码包的“基准输出”你知道一个 1 W、1 ps 高斯脉冲在 β₂-20 ps²/m、γ0.02 W⁻¹m⁻¹、100 m 光纤里应该长什么样。4. 必调参数与单位换算仿真结果不“玄学”的关键4.1 单位、前缀和材料参数差一个 1000 倍就是事故这个标题下的代码包里最常出现的“玄学”问题根本不是算法而是单位。非线性光学领域文章爱用混合单位色散用 ps²/km损耗用 dB/km非线性系数用 W⁻¹·km⁻¹功率用 mW长度用 m。你把这些数字原样敲进上面那段脚本结果一定离谱。先记住一个换算表它是我复现任何别人代码前都要核对一遍的清单。物理量常见写法 A常见写法 B换算关系二阶色散 β₂ps²/kmps²/m除以 1000非线性系数 γW⁻¹·km⁻¹W⁻¹·m⁻¹除以 1000损耗 αdB/kmm⁻¹除以 8.686脉冲宽度fsps除以 1000波长nmμm除以 1000最阴险的是 γ。同样一段光纤资料 A 写 1.5 W⁻¹·km⁻¹资料 B 写 0.0015 W⁻¹·m⁻¹数字差一千倍但物理量相同。如果你把 1.5 当成 W⁻¹·m⁻¹ 直接填进脚本Kerr 相移会被放大一千倍脉冲传播不到几米就发散。每拿到一组材料参数第一步先统一成“米制组合”再进代码这一步省掉后面所有痛苦。还有一个单位制陷阱出现在材料非线性折射率 n₂ 上。很多材料手册用的是 esu 单位制而仿真公式里要求的是 m²/W。这两者的换算没有一个简单的整数倍关系跟材料体系有关。我的规避办法是尽可能只从同一篇文献里拿同一组单位制的数据如果必须混用先从该文献的公式推导里反推出它用的是哪套单位制而不是凭感觉猜。这里最容易让人踩坑的是某些文献的 n₂ 值其实是“强度归一化”后的等效值只适用于该文献特定的有效模面积假设。β₂ 的符号也是老生常谈但反复出错的点。反常色散β₂0会让人看到脉冲压缩、孤子和调制不稳定性正常色散β₂0则只有展宽。很多工程问题里泵浦波长选在零色散波长附近此时 β₂ 很小β₃ 反而成了主导参数如果你的代码只写到了 β₂结果自然对不上。拿到一套色散数据时先看零色散波长在哪再决定你需要的色散阶数。4.2 采样点数、时域窗口与纵向步长怎么设数值参数的选择直接决定结果可信度。首先引入两个物理长度色散长度 L_D T₀²/|β₂| 和非线性长度 L_NL 1/(γP₀)。它们衡量的是“跑多远色散效应和非线性效应才会彻底改变脉冲形状”。分步傅里叶法的纵向步长必须同时远小于两者工程经验是取 min(L_D, L_NL) 的 1/50 到 1/100 作为 Δz。比如 T₀1 ps、β₂-20 ps²/m 时 L_D0.05 mGVD 效应比非线性效应先起作用那 Δz 就应该取 0.001 m 这个量级而不是 0.2 m。如果输入的功率很低又不关注色散细节倒是可以适当放大步长换速度。时域窗口 Twin 必须足够大到“装下”整个脉冲在整个传播过程中的全部包络。初始脉冲可能只占窗口的十分之一但随着色散展宽或孤子演化脉冲会变宽甚至分裂窗口太小会切掉尾部产生数值反射。我一般至少让窗口覆盖初始脉宽的 20 倍并且在结果检查时看一眼时域末端有没有异常突起。如果出现了说明窗口和步长里至少有一个不合适。采样点数 Nt 直接决定奈奎斯特频率 f_max ≈ 1/(2Δt)也就是能表征的最大频率分量。自相位调制会产生新的频谱分量频率偏移量大约等于 γP₀ 对应的相位调制速率脉宽越窄、功率越高需要的高频容量越大。我的做法是先把 Nt 设成 2¹²看频谱在奈奎斯特边界附近是否出现“撞墙”一样的噪声如果撞墙就说明高频分量溢出需要更大采样点或更大时域窗口。5. 非线性光学仿真避坑指南五个典型翻车现场5.1 总能量随时间“凭空增长”有一个很典型的翻车现象关闭损耗后传播到一半计算出的总能量 ∑|A|² 开始慢慢变大最后明显超过初始脉冲总能量。如果总能量只在一个 dz 附近波动且幅度很小那是分步法截断误差的体现如果是持续单调增长几乎可以肯定是线性步或非线性步里乘法因子写错了。常见原因是在频域线性步里把 exp(iβ₂ω²Δz/2) 写成了 exp(iβ₂ω²Δz)前半步和后半步相位因子重复累加等效于每步都多走了一半线性距离。排查方法是先设置 γ0、α0只跑线性传播此时分步法的能量应该严格守恒到机器精度。如果这一关就过不了先修代码里的算符再进入非线性调试。5.2 频域出现对称的伪峰输出频谱时主峰两侧出现两个高度对称的小峰这通常不是物理上的调制不稳定性而是数值混叠。原因要么是时域采样点不够高频成分折叠回低频要么是时域窗口太小脉冲在窗口边缘周期性复制后与自身干涉。区分方式很简单把窗口扩大一倍伪峰仍存在且位置随窗口移动就是混叠伪峰消失就是原来的窗口确实截断了频域信息。解决方法是同时提高 Nt 和 TwinNt 按 2 的幂增长Twin 至少覆盖脉冲动态范围的全程。还有一个习惯容易引发这个问题画频谱时忘记 fftshift直接把未移位的数组画出来会在图两端看到一条看似噪声的强度带那不是物理是坐标轴没归零。5.3 功率稍微一调大输出直接变 NaN脉冲峰值功率从 1 W 提到 5 W程序立刻 NaN这是老手也躲不过的筛查重灾区。直接原因是非线性相位因子 exp(iγ|A|²Δz) 在单步内增长过快数值上出现了幅度发散。本质上则是步长没有按 L_NL 缩放功率增大五倍L_NL 缩小五倍步长不变就踩进了不稳定区。另一个隐蔽原因是在飞秒脉冲仿真里加入了拉曼增益项后增益谱的最高频率部分获得指数增长而没有对应损耗来抑制就会在几百步内数值爆炸。处理办法先把拉曼响应函数去掉确认骨架稳定再把 Δz 缩到新 L_NL 的 1/50 以下如果仍不稳定换成自适应步长积分器。不要一上来就调代码的数值格式先改步长这是最低成本的后悔药。5.4 材料参数和别人对不上复现某篇文献时光谱形状大体对但横轴宽度差了几倍或者自相位调制产生的峰值间距不对这种情况九成出在材料参数。最常见的是把 β₂ 符号取反把反常色散写成了正常色散脉冲演化规律直接反转其次是 γ 计算时用的 n₂ 波长跟仿真波长不一致色散已经被软件记了一笔账。还有一类问题是有效模面积 A_eff 用的数值和原文不统一某些论文给的是模场直径你换算成面积时用的是 π(d/2)² 还是 πd²/4 都容易对不上。我的习惯是把所有输入参数写到一个配置文件里同行评审式地逐项标注来源文献、适用波长、单位制。一旦结果不对别人或者两周后的我都能快速定位哪一项可疑。5.5 复现论文时频谱“差一截”就差一步但很难受的情况末端频谱的宽度量级对但细节结构完全不同比如论文有明显的孤子峰而你这边是平滑连续谱。最常见的原因不是数值错而是初始条件不一致。论文里用的是双曲正割脉冲 sech(t/T₀)而你用了高斯脉冲或者论文里的所谓峰值功率其实是平均功率你得换算还有的是初始脉冲带了一个小 chirp这在实验里非常常见而仿真里没人刻意加。解决办法是先做输入条件消歧把高斯换成双曲正割跑一次再按论文给的时间-带宽乘积反推初始 chirp跑一次都不对再检查算法里是否漏了拉曼延迟项。到这一步通常问题就水落石出了。6. 进阶验证法用孤子数 N 判断你的仿真离超连续谱还有多远第 5 章处理完之后你应该有了一套能稳定输出的骨架。这时候别急着上复杂介质先用孤子数 N 做一通验证这是最能证明代码正确性的测试。孤子数定义为 N² L_D / L_NL γP₀T₀²/|β₂|它是非线性长度和色散长度的比值直接刻画了传播过程中谁的力量更强。当 N1 时理想的无损耗标准 NLSE 中脉冲应该保持双曲正割形状不变这是一个不需要任何复杂物理就能检验代码的基准跑完 10 个 L_D波形如果明显畸变那代码大概率有 bug。N2 时脉冲会在传播中先压缩再分裂再复原频谱上出现明显的振荡结构。N 大于 3 之后高阶孤子受到噪声或高阶色散扰动就容易发生裂变这正是超连续谱产生的最常见物理机制。我自己带人做这类仿真时有一个固定流程先固定 β₂、γ调整 P₀ 让 N 分别等于 1、2、3各跑一遍观察演化是否符合孤子周期公式。N1 孤子周期是无穷大形状不变N2 孤子周期是 (π/2)·L_DN3 则能明显看到分裂再重组。如果你的代码能复现这三个预期再进入真实材料参数才有意义。这个“先证明工具可靠、再信任结果”的习惯帮我和当时的模拟项目 X 省下了大量排错时间。做网格收敛性检查也值得坚持把 Δz 减半、Nt 加倍各跑一遍三组结果叠加画在同一张图里曲线几乎重合说明数值方案已经收敛如果三组曲线还在“打架”那当前任何一组数据都不能作为物理结论引用。输出方面除了画图我一般会把末端时域和频域数据存成 CSV 或 MAT方便和实验光谱做定量对比。这一步不需要额外依赖就把 np 数组按列写入文本文件就行。最后说一个自己的习惯不要迷信“代码能跑”就等于“结果正确”。我第一次做飞秒脉冲仿真时代码跑得很顺频谱也很平滑直到我用 N1 验证才发现脉冲没守住形状问题出在我把 β₂ 的单位写错了一个数量级整个仿真物理上都不自洽。从那以后每套新的非线性仿真代码我都会先跑通最低配置的基准测试再相信输出。希望这篇笔记能帮你把弯路走直一点。本文还有配套的精品资源点击获取
