简介：面向非线性光学仿真与数值计算的完整工程资料包，适合光学工程、物理及相关专业的学生与研究人员，侧重二次谐波产生、参量下转换、四波混频等典型非线性效应，覆盖光场传播、相位匹配、材料极化率与转换效率计算等核心问题。压缩包内共包含五百七十七个文件，核心文件类型包括MATLAB脚本、C++源文件与头文件、预编译的MEX内核，以及Visual Studio工程配置、CHM帮助文档、数据文件与图示文件；全部内容压缩后约五点二三MB，目录结构清晰，便于按功能模块检索和学习。当前已有八百五十三人学习或下载。压缩包内提供了信号计算、最小误差内核修正、二次谐波求解等仿真代码与MEX工具，并配有工程配置和说明文档，可帮助学习者将非线性光学原理转化为可运行的程序实践，也适合在此基础上扩展新的数值实验，深入理解高功率光场与材料相互作用时的各种非线性效应。

1. 拿到 Nonlinear-Optics 代码包，先搞清楚它在算什么：相位匹配比非线性方程更要命

做光学仿真的人第一次打开 Nonlinear-Optics 这类代码包，第一反应多半是去找非线性薛定谔方程的求解器。我的经验是反过来的：这类项目里最容易翻车的，从来不是方程本身解不动，而是你根本没把相位匹配算对。相位匹配没对准，再高阶的积分器、再细的网格都白搭。这篇文章把非线性光学仿真最常用的一条技术路线讲透：从 I 类倍频切入，用耦合波方程加 RK4 把转换效率算出来，再逐一扫角度、温度这些调谐参数，最后落在代码和事务怎么组织上。适合正在做倍频、OPO、OPA 或者激光频率变换的工程师和研究生，照着改参数就能复现。

2. 先把物理模型立住：从 I 类倍频到三波互作用的耦合波方程

2.1 为什么先啃 I 类倍频：三波互作用里最省事的试金石

倍频（SHG）是非线性光学里最经典、也最容易做数值验证的过程：一束 1064 nm 的基频光打进晶体，出来一束 532 nm 的倍频光。它背后是三波互作用最简形式——两个基频光子湮灭，产生一个倍频光子，信道只有一个，极化率只看 χ⁽²⁾。

我一般建议新手或者新项目先从 I 类倍频入手，理由有三个。第一，实验数据遍地都是，BBO、LBO、KTP 的倍频效率曲线在文献里能查到一大堆，仿真对得上对不上，一眼就能看出来。第二，它的耦合波方程只有两个，数值上是最低门槛，却已经包含了后面 OPO、OPA 的几乎所有核心要素：三个波、一个相位失配量、一组耦合系数。第三，从 I 类倍频出发，把方程从两个扩到三个，就是完整的光参量放大（OPA），再加一个谐振腔边界条件就是 OPO。所以这一章把倍频的物理模型啃下来，后面全是机械扩展。

很多人会问，为什么不用 FDTD 或者非线性薛定谔方程来做。FDTD 确实能算非线性过程，但代价是三维网格加时间步进，一个毫米级晶体的计算量非常可观，而且它给的是场分布，不是工程师最关心的转换效率和匹配带宽。连续波或者纳秒脉冲场景下，慢变包络近似下的耦合波方程是行业里最常用的做法，几行 RK4 就能解决，没必要上重型求解器。至于非线性薛定谔方程，那是处理自聚焦、自相位调制这类效应时才需要的，倍频场景里如果泵浦功率还没到损伤阈值，这些效应通常可以放一边。

2.2 耦合波方程怎么落到可数值求解的形式

从麦克斯韦方程组出发，做三个假设：平面波、无损耗晶体、慢变包络近似。慢变包络的意思是，场的振幅包络在波长尺度上变化很慢，所以可以把快振荡的相位项单独拎出来处理。做完这些假设，二阶非线性极化项会以驱动项的形式出现在波动方程里，最后整理成一组一阶常微分方程，也就是耦合波方程。

对 I 类倍频（两个基频光子合成一个倍频光子），场振幅的演化可以写成：

dE1/dz = i (ω1·deff)/(n1·c) · E2·E1*·exp(−iΔkz)

dE2/dz = i (ω2·deff)/(2·n2·c) · E1²·exp(iΔkz)

E1 是基频场，E2 是倍频场，deff 是有效非线性系数，Δk = k2 − 2k1 是相位失配量。第一个方程描述泵浦光在转换过程中的损耗，第二个方程描述倍频光的增长。注意第二个方程前面那个 1/2 因子，它来源于 E² 展开后 2ω 分量前面的系数，丢了它，能量守恒直接破功，效率算出超过 100% 就是从这里来的。

这套方程本身没有解析解（除了小信号近似下），所以它的落地形式必然是数值积分。把实部和虚部分开，或者直接用复数场推进，沿着晶体长度 z 一步一步算，每一小步里场和场之间的耦合按上面的右端项更新。参数上的关键点在于，deff 的量纲是 m/V，n1、n2 是折射率，ω1、ω2 是圆频率，这几个量一旦单位混了，结果会差出几个数量级，后面第 5 章专门讲这个坑。

2.3 相位匹配条件与失配量的物理直觉

为什么一定要相位匹配？因为介质有正常色散，基频光和倍频光在晶体里的折射率不一样，走的速度不一样。两列波在前进过程中不断累积相位差，转换出来的倍频光和后面新产生的倍频光相位对不上，就会发生能量倒流——先转过去的又转回来。这就是所谓 back-conversion，效率曲线长到一定长度后不升反降的原因。

所以 Δk 是整件事的灵魂。Δk 为零时，转换效率随长度近似二次方增长；Δk 不为零时，效率变成振荡的，峰值被压低。真实晶体里 Δk 几乎不可能天然为零，因为任何材料的色散都让 k2 大于 2k1。工程师做的事，就是通过角度、温度、周期极化这三种手段把 Δk 压回零附近。

相位匹配不是一个点，而是一个窗口。窗口的宽容度由晶体长度决定，长度越长，允许的 Δk 偏差越小，角度或温度的容忍度就越窄。这个直觉在调参时非常重要——你扫角度扫出来的峰，宽度往往只有零点几度，不是程序写错了，是物理本来就这个脾气。

3. 用 Python 把倍频效率算出来：一套可直接改的 RK4 脚本

3.1 最小脚本：RK4 积分两个耦合方程

先定义参数区，把所有物理量统一到 SI，避免中途换算：

import numpy as np # 物理常数 c = 2.99792458e8 # 光速, m/s # 晶体与光束参数（BBO I类，基频 1064 nm -> 倍频 532 nm） wl1 = 1.064e-6 # 基频波长, m wl2 = wl1 / 2.0 # 倍频波长, m n1 = 1.655 # 基频 o 光折射率（近似，精确值查 Sellmeier 方程） n2 = 1.556 # 倍频 e 光在匹配角附近的折射率（近似） deff = 2.0e-12 # 有效非线性系数, m/V（数量级近似） I0 = 500e4 # 峰值功率密度, W/m^2，即 500 MW/cm^2 Lc = 3.0e-3 # 晶体长度, m w1 = 2 * np.pi * c / wl1 w2 = 2 * np.pi * c / wl2 k1 = 2 * np.pi * n1 / wl1 k2 = 2 * np.pi * n2 / wl2 dk = k2 - 2 * k1 # 相位失配量, 1/m # 振幅归一化：让 |A|^2 直接对应功率密度（W/m^2） A10 = np.sqrt(I0) A20 = 0.0 # 两个方程的耦合系数 k1_cpl = w1 * deff / (n1 * c) k2_cpl = w2 * deff / (2 * n2 * c)

这套初始化的思路是，不把场振幅定义成电场强度 V/m，而是把 |A|² 直接映射到功率密度。好处是末端效率可以直接用 |A2|²/(|A1|²+|A2|²) 算，不用再转换单位。坏处是耦合系数的表达式跟着变，但这个变是透明的，所有量都在 SI 下，不容易错。

然后是 RK4 积分器：

def shg_rk4(z_grid, A1_0, A2_0, dk, k1_cpl, k2_cpl): # 恒定步长 RK4，每一步推进两个复振幅 A1, A2 = A1_0, A2_0 for i in range(len(z_grid) - 1): h = z_grid[i+1] - z_grid[i] def f(z, y): a1, a2 = y[0], y[1] return np.array([ -1j * k1_cpl * a2 * np.conj(a1) * np.exp(-1j * dk * z), -1j * k2_cpl * a1 * a1 * np.exp(1j * dk * z), ]) y = np.array([A1, A2], dtype=complex) k1v = f(z_grid[i], y) k2v = f(z_grid[i] + h/2, y + h*k1v/2) k3v = f(z_grid[i] + h/2, y + h*k2v/2) k4v = f(z_grid[i] + h, y + h*k3v) y = y + h/6 * (k1v + 2*k2v + 2*k3v + k4v) A1, A2 = y[0], y[1] return A1, A2

主流程：

z_grid = np.linspace(0, Lc, 4000) A1_end, A2_end = shg_rk4(z_grid, A10, A20, dk, k1_cpl, k2_cpl) eta = abs(A2_end)**2 / (abs(A1_end)**2 + abs(A2_end)**2) print(f"晶体末端倍频效率: {eta*100:.2f}%")

逻辑说明：右端项里 exp(−iΔkz) 这个快振荡项不能忽略，它代表两波之间的相位滑移，RK4 的优势就在处理这种带振荡相位的耦合方程。为什么不用欧拉法？因为欧拉法的局部截断误差是 O(h²)，对 3 mm 晶体、4000 步这种尺度，误差累积会让效率曲线出现锯齿。RK4 的局部误差是 O(h⁵)，同样步数下精度高得多。

参数方面，步数 4000 不是拍脑袋。基频在晶体里的波长约 0.64 μm，3 mm 晶体约 4700 个波长，步长取 0.75 μm 左右，能保证每个波长周期内有接近一个采样点，包络近似才站得住。如果你要扫大 Δk，步长上限要用 π/|Δk| 来约束，后面第 5 章细说。

3.2 从小信号到泵浦耗尽：效率曲线怎么看

把这段代码跑起来，你会看到效率随晶体长度并不是一路涨到 100% 就停下来。小信号区（转换效率低于 10%）效率近似正比于 L²，这是解析解能覆盖的范围。继续涨到 50% 以上，泵浦开始被明显消耗，曲线斜率变缓。再往后可能出现效率回落，这就是 back-conversion——倍频光又差频回去了，能量倒流。

很多人第一次看到效率回落以为自己 bug 了。不是，是 Δk 不为零或者长度超过相干长度时，三波之间的相位关系反转了。这个现象在真实实验里也看得到，不是数值伪影。判断标准很简单：把 Δk 设成零，看效率是否单调饱和并趋近 100%，如果这个行为对了，说明程序逻辑没问题，回落现象是物理，不用修。

另一个验证手段是能量守恒检验。每走一步，|A1|² + |A2|² 应该保持常数，误差在 1e-6 量级以内。如果这条不成立，优先怀疑耦合系数写错了。这是我每次跑完新参数必做的一件事，比盯着效率曲线本身更有用。

3.3 参数表：哪些量要调，调了会怎样

注意 dk 很少是自由参数。它由角度、温度、周期极化决定，所以实际项目里你调的是晶体角度或者烘烤温度，而不是直接改 dk。下一章就讲这两个最常用的调谐手段，以及它们各自带来的副作用。

4. 相位匹配调谐的落地参数：角度、温度与走离效应

4.1 角度调谐：把 θ 扫一遍，找到 I 类匹配峰与带宽

BBO 是负单轴晶体，I 类匹配配置是 oo→e，也就是两个基频光子都是寻常光，倍频光子是非寻常光。非寻常光折射率随传播方向与光轴夹角 θ 变化，关系是：

n_e(θ) = n_o·n_eo / sqrt(n_o²·sin²θ + n_eo²·cos²θ)

n_eo 是主轴方向上的非寻常光折射率。θ 一变，k2 跟着变，Δk 就能扫过零。所以角度调谐的本质是，用几何关系补偿材料色散。下面这段代码把角度扫描接进上一章的求解器：

# 角度扫描：在 20°~30° 范围内扫 BBO 的相位匹配角 n_o = 1.655 n_eo = 1.556 theta_deg = np.linspace(20, 30, 2001) theta = np.deg2rad(theta_deg) n2_theta = (n_o * n_eo) / np.sqrt( n_o**2 * np.sin(theta)**2 + n_eo**2 * np.cos(theta)**2 ) # 每个角度对应一个新的 dk，重跑 RK4 eta_scan = [] for n2_i in n2_theta: dk_i = 2*np.pi*n2_i/wl2 - 2*(2*np.pi*n1/wl1) A1_end, A2_end = shg_rk4( z_grid, A10, A20, dk_i, k1_cpl, k2_cpl ) eta_scan.append(abs(A2_end)**2 / (abs(A1_end)**2 + abs(A2_end)**2)) eta_scan = np.array(eta_scan) best_idx = np.argmax(eta_scan) print(f"最佳匹配角: {theta_deg[best_idx]:.2f}°") print(f"峰值效率: {eta_scan[best_idx]*100:.2f}%")

跑出来的峰值角度，文献上 BBO 对 1064 nm 的 I 类匹配角一般在 22.8° 上下，你的 Sellmeier 常数版本不同会差个零点几度，这很正常。真正要注意的是峰的宽度：半高全宽通常只有零点几度，稍微偏一点效率就掉一半。这是因为 3 mm 晶体的相干长度有限，Δk 稍大，整个晶体长度内累积的相位滑移就超过 π。

角度容忍度这个参数直接决定实验里要用什么精度的旋转台。如果仿真显示容忍度只有 0.1°，那手动调节基本没戏，得上带刻度的步进电机。这一条信息比效率绝对值更有工程价值。

4.2 温度调谐与 90° 非临界相位匹配

角度调谐有个天然副作用，就是倍频光相对基频光会走离，下一节细说。如果想躲开走离，就得让传播方向垂直于光轴，也就是 θ = 90°，这时非寻常光和寻常光的传播方向一致，但 Δk 通常不为零。怎么办？换温度。以 LiNbO₃ 为例，它的折射率随温度变化很明显，升温可以把 Δk 压回零，这就是所谓 90° 非临界相位匹配。

在脚本里做温度扫描，比角度扫描更简单：找一个温度依赖的 Sellmeier 方程，算出每个温度下的 n1(T)、n2(T)，然后套同一套 RK4 循环。折射率对温度的变化率大约在 10⁻⁵/K 到 10⁻⁴/K 量级，换算成温度带宽，通常也是几度的量级，比角度容忍度宽松得多。

另一个完全不同的思路是周期极化，比如 PPLN。它不靠温度也不靠角度，而是让晶体的自发极化方向周期翻转，每过一个相干长度就把相位误差掰回来。极化周期 Λ = 2π/Δk，设计好周期之后，可以在很宽的温度范围里做到非临界匹配，同时用上晶体最大的有效非线性系数。仿真上，周期极化只是在右端项里把 dk 换成 dk − 2π/Λ，代码改动极小，但物理效果完全不同。

4.3 走离效应：为什么临界匹配的晶体不能开大口径

走离是角度调谐绕不开的代价。非寻常光在晶体里传播时，能量传播方向（坡印廷矢量）和波矢方向不一致，夹角叫走离角。BBO 的走离角大约 3°，看起来不大，但基频光和倍频光在横向上会慢慢分开，重叠区域变小，有效作用长度变短。

数值上怎么处理？一种常见做法是把晶体长度 L 换成有效长度 L_eff = 光束直径 / tan(ρ)，然后再跑同一套方程。比如光束直径 1 mm，tan(3°) 约 0.052，有效长度约 19 mm，比晶体长，走离可以忽略；但光束直径 100 μm 时，有效长度只有约 1.9 mm，这时候再不建模走离，效率会被高估不少。

所以选型的时候有一个基本权衡：

如果你的应用要求大光束、高效率、长晶体，优先考虑非临界或者周期极化方案；如果只是小口径、短晶体、追求简单，BBO 角度调谐完全够用。这个取舍在做系统设计时比盲目上高非线性系数晶体重要得多。

5. 非线性光学仿真常踩的 5 个坑与排查路径

5.1 坑一：效率曲线在中段出现锯齿或振荡

现象：eta 随晶体长度不是光滑上升，而是出现周期性的锯齿，甚至效率在某个位置跳变。

原因：步长不够。当 Δk 很大时，右端项里的 exp(−iΔkz) 振荡非常快，步长必须小于 π/|Δk|，否则采样点落在振荡周期的错误位置上，RK4 会把振荡混叠成低频误差。

解决：先算出 π/|Δk|，然后把步长设成它的 1/5 到 1/10。比如 |Δk| = 10⁶ /m 时，π/|Δk| 约 3 μm，步长取 0.3 μm，对 3 mm 晶体就是 10000 步。不要嫌多，这个场景下步数不是性能瓶颈。

5.2 坑二：泵浦耗尽后效率超过 100%

现象：|A2|²/(|A1|²+|A2|²) 算出来 1.2，能量无中生有。

原因：第二个方程的 1/2 因子丢了。倍频方程里 E1² 本身带有 2ω 分量系数，很多人在推导时只记了形式，把 ω2·deff/(n2·c) 直接写上，漏掉分母里的 2。这个错误在小信号区（效率 < 10%）几乎看不出来，一旦进入泵浦耗尽区就爆发。

解决：在 RK4 循环里加一行断言：abs(A1)**2 + abs(A2)**2 的极差小于 1e-6，不满足就直接打印错误退出。这个校验应该成为所有耦合波方程的标配，换晶体、换波长之后都先跑一遍。

5.3 坑三：角度扫描整个区间都找不到峰值

现象：把角度从 10° 扫到 40°，效率全程接近 0，或者峰值出现在明显不合理的位置。

原因：折射率数据用错了。最常见的是把文献里某一温度下的 Sellmeier 常数直接拿来算 1064 nm，没注意那个常数只覆盖某个波长范围；另一个常见问题是偏振配置搞反，oo→e 写成了 oe→o，计算出来的 Δk 永远不为零。

解决：先在脚本里打印 Δk 的扫描范围。如果 Δk 全程没有跨过零，那不管扫多少个角度都不可能找到峰值。然后回去核对 Sellmeier 常数的波长适用范围和温度条件。我在换晶体参数时，总是先用小信号解析解算一个粗略匹配角，再上 RK4，匹配角差太远就先修数据，不急着看图。

5.4 坑四：单位制混乱，结果差几个数量级

现象：同样的物理参数，不同人跑出来的效率差一万倍，而且谁都觉得自己没错。

原因：cm 和 m 混用是最常见的。功率密度有人用 W/cm²，有人用 W/m²，折射率倒是无量纲，但波长一会儿 μm 一会儿 m，最后 k 的量级差出 4 个次方。pm/V 和 m/V 的有效非线性系数也经常被直接套用。

解决：脚本头部单独做一个换算区，所有物理量一次性转到 SI，后续变量只引用不换算。比如功率密度统一写成 W/m²，在参数区用 500e4 表示 500 MW/cm²，但注释里写明原始工程单位。这个习惯能省掉大量排错时间，属于血泪经验。

5.5 坑五：仿真结果和 SNLO 对不上

现象：同样的晶体、波长、功率密度，自己的脚本算出来效率 60%，SNLO 给 45%，带宽也有偏差。

原因：SNLO 这类工程软件内部做了平面波近似之外的修正，包括走离平均、孔径效应、聚焦参数、高斯光束的重叠积分。你的耦合波方程是理想平面波，两者算的本来就不是同一个物理对象。这不叫 bug，叫模型边界差异。

解决：先把场景设成小信号、薄晶体、低效率，这个区间里解析解、你的数值解、SNLO 三者应该高度一致。对不上就先查自己的数据；对得上，再逐步提高功率、加长晶体，观察两者分叉。分叉本身就是信息：它告诉你什么时候平面波模型失效，什么时候必须上更完整的模型。盲目追求和 SNLO 完全一致，很多时候是在浪费时间调一个本来就不该一致的参数。

6. 从倍频走向 OPO、OPA：脚本扩展与 master 分支代码组织

倍频的方程组只有两个，扩展到 OPA 就是再加一个闲频光方程，结构完全一样：

# OPA 三波耦合方程骨架（信号光、闲频光、泵浦光） def rhs_opa(z, y, kappa_s, kappa_i, kappa_p, dk): As, Ai, Ap = y[0], y[1], y[2] return np.array([ -1j * kappa_s * Ap * np.conj(Ai) * np.exp(-1j * dk * z), -1j * kappa_i * Ap * np.conj(As) * np.exp(-1j * dk * z), -1j * kappa_p * As * Ai * np.exp(1j * dk * z), ])

这里的 dk = kp − ks − ki，符号定义和倍频场景一致。OPA 的物理行为和倍频完全不同：泵浦功率超过阈值后，信号光和闲频光会指数增长，增益曲线对波长非常敏感，所以它的仿真重点从效率计算变成了增益谱和阈值扫描。脚本骨架不变，RK4 循环不用动，只换右端项和初值，这就是我推荐先做倍频的原因——后面全是复用。

代码组织上，我现在的习惯是把这类 Nonlinear-Optics 项目按 src、sims、results 三块分，master 分支只保留可复现的版本。有人问 master 上的代码怎么切到 dev，git checkout -b dev 一条命令的事，但在科学计算项目里，真正该做的不是频繁开分支，而是每次跑批前在 config 文件里固定参数、跑完打 tag。tag 里记的是“这个结果对应哪一组晶体参数、哪个温度、哪个功率密度”，比分支名可靠得多。我自己曾经吃过没打 tag 的亏：同一份代码改了 Sellmeier 常数忘了记录，两周后对着旧结果怎么都对不上。现在每个仿真任务新建一个 config 文件，task 名和 tag 同步，后悔药算是吃上了。

现在拿到一个非线性频率变换题目，我第一件事不是找求解器，而是先查这批晶体的折射率数据对不对、匹配角落在哪个区间、走离影响大不大，这些确认了再写 RK4。这个习惯救了我好几次，希望帮到你。

本文还有配套的精品资源，点击获取