简介：等离子体高次谐波产生（HHG）的Matlab仿真源代码，专为光学、原子物理与超快激光领域的研究人员和学生设计，可用于从量子力学层面完整模拟强激光与气体相互作用的非线性过程。代码覆盖喷嘴出口超音速气流、等离子体动态离子、自由电子量子响应、相位匹配与吸收等核心物理环节，并配有用户友好的GUI，允许读者自由调节输入参数即时观察对谐波幅度、电离率等结果的影响。资源共451个文件，以Matlab脚本（.m）和计算数据文件（.txt）为主，辅以PDF说明文档、Markdown笔记及图型文件，压缩包整体29.75MB。所有计算过程均输出至.txt保存，相同参数二次运行时能显著缩短仿真周期，项目中已内置常用运行结果，便于对照学习。尤其适合科研人员用于机理验证、参数扫描或教学演示，已有716人浏览学习，可作为高次谐波数值模拟的重要参考代码库。

1. 高次谐波仿真：为什么一个MATLAB仓库值得你花一个下午跑通

高次谐波产生（High-Harmonic Generation, HHG）是强激光与原子、等离子体相互作用时出现的极端非线性光学现象，谐波能延伸至极紫外甚至软X射线波段，也是阿秒脉冲产生的核心物理机制。这份MATLAB仿真源代码把含时薛定谔方程（TDSE）的数值解算完整封装，跑通一次就能输出拿到谐波谱，看到典型的“平台区+截止频率”结构，整个系统从激光脉冲构造、波函数演化到谐波谱提取全部覆盖。它不只是一段物理演示代码，更是一套可运行、可调参、可换原子的数值仿真系统。适合两类人：一是在学强场物理、超快光学的研究生，想复现教科书上的经典结果；二是做数值仿真开发的人，想参考split-operator算法在MATLAB里的工程化写法。仓库不依赖Simulink，纯脚本就能跑，MATLAB版本要求不高，本地环境干净，跑通一次后面改参就非常顺手。

2. 物理模型与代码结构：高次谐波仿真的黑匣子拆开看

2.1 三步模型：隧穿电离、经典加速、复合辐射

高次谐波产生的物理图像，在强场物理里最常用的解释是“三步模型”（Three-Step Model），也叫再碰撞模型。第一步，强激光场把原子势垒压窄，电子通过隧穿电离逃出束缚；第二步，自由电子在激光电场中做经典加速运动，电场反向时电子被拉回母核附近；第三步，电子与母核复合，把从激光场中获得的动能以光子形式辐射出去。这个辐射光子的能量可以远高于激光光子本身，于是产生了高次谐波。

三步模型是半个经典的图像，而代码里实际求解的是含时薛定谔方程，是真正的量子力学处理。三步模型的价值在于它给出了谐波谱截止频率的半经典公式：谐波阶次约等于（离子化能 + 3.17倍有质动力能量）除以激光光子能量。跑完仿真后，你可以用这个公式验证数值结果——这个交叉验证本身就是理解这份代码的关键路径。如果数值结果偏离半经典公式超过20%，大概率是数值参数没配对，而不是物理模型错了，这是排查时首先要想清楚的。

这里有一个容易困惑的点：为什么代码里不直接做经典轨迹采样，而要费劲解TDSE？因为三步模型只描述了一个电子的经典轨迹，但真实物理是大量电子波包相干叠加，量子干涉效应（比如谐波谱平台区的干涉结构）只有波函数演化才能复现。TDSE解算能同时给出电离、加速、复合三个过程的完整量子信息，代价就是计算量上了一个台阶。

2.2 数值方法选择：为什么用split-operator而非直接差分

如果你对TDSE直接做有限差分，时间步长要压到很小才能保证稳定，因为动能算符是二阶导，显式格式的稳定性条件相当苛刻，步长不合适当场发散给你看。隐式格式（比如Crank-Nicolson）稳定性好一些，但每个时间步都要解一个大型线性方程组，在MATLAB里就意味着频繁的大矩阵计算，效率感人。Split-operator（分裂算符法）的思路完全不同：把演化算符exp(-iHΔt)近似拆成exp(-iVΔt/2)·exp(-iTΔt)·exp(-iVΔt/2)，其中T是动能项，V是势能项。势能项在坐标表象下是对角矩阵，动能项在动量表象下是对角矩阵，两次傅里叶变换就能完成一步演化，全程不需要解任何线性方程组。

这个方案在MATLAB里特别顺手，因为fft是内建函数，矩阵化操作一次搞定一个时间步。代码里每个时间步的核心操作就是一次fft和一次ifft，向量化程度很高，不需要写任何显式循环。代价是边界条件处理得更小心——周期性边界意味着电子波包跑出计算区域会从另一边穿进来，必须在边界上施加吸收掩蔽函数，这部分后面会详细讲。

踩坑提示：split-operator不是唯一的解题思路，也有人用Chebyshev多项式展开演化算符，精度更高但代码复杂度上去了。这个仓库采用分裂算符是合理的工程权衡——在“够用”的精度下拿到最快的速度和最简单的代码。对于学生做课题或者工程师做预研，这个选择完全够用。

2.3 目录结构与核心文件职责

从资源打包的目录结构看，这是一个典型的单脚本驱动型MATLAB工程，没有复杂的类层次，核心逻辑集中在几个函数文件里。

. ├── main.m # 主脚本，定义物理参数并驱动整个仿真流程 ├── init_wavefunction.m # 初始化基态波函数，用虚时间演化求目标原子基态 ├── laser_pulse.m # 构造激光场，sin²包络 + 载波相位 ├── propagate.m # 分裂算符演化，核心时间步进函数 ├── extract_hhg.m # 计算偶极矩加速度并提取谐波谱 ├── parameters.m # 集中管理所有物理参数 └── utils/ ├── absorbing_boundary.m # 吸收边界掩蔽函数，防波函数边界反射 └── plot_spectrum.m # 绘图工具，半对数坐标画谐波谱

每个文件的职责边界很清晰：parameters.m只负责定义参数结构体params，不执行计算；main.m只负责调用流程，不写具体算法；propagate.m只做一步波函数推进，不关心物理参数含义。这种分层方式对二次开发很友好——你想改激光参数就动parameters.m，想换数值方法就只动propagate.m，互不牵连。

其中init_wavefunction.m是最容易被低估的一个文件。它不是随便给一个高斯波包就开始演化，而是先做虚时间演化（imaginary time propagation）求出目标原子势下的基态波函数。这一步意义很大：用真实基态作为初始态，电离动力学才是物理上正确的，否则会混入大量人为的瞬态响应噪声。我第一次跑类似代码时偷懒直接给高斯波包，结果谐波谱上全是非物理的毛刺峰，后来加上虚时间演化才拿到干净的平台区，这个血泪经验值得记一下。

3. 核心模块拆解：一个波函数从初始化到谐波谱输出的完整链路

3.1 激光脉冲构造：sin²包络与载波相位的实现

function E = laser_pulse(t, params) % 构造沿x方向偏振的激光电场 % t : 时间向量（原子单位制） % params: 结构体，包含omega0（载波频率）、t0（峰值时刻）、tau（脉宽） % 返回 : 电场强度数组，单位与原子单位制下的场强一致 % sin²包络：在[t0-2*tau, t0+2*tau]区间内非零，平滑起止避免频谱泄漏 envelope = sin(pi * (t - t0) / (4 * tau)) .^ 2; envelope(t < t0 - 2*tau | t > t0 + 2*tau) = 0; % 载波：cos(omega0*t)表示偏振轴上的电场振荡 % 如果改成sin，则载波包络相位（CEP）不同，谐波奇偶次比例会变化 carrier = cos(params.omega0 * t); % 叠加：峰值场强E0决定有质动力能Up，从而决定截止频率 E = params.E0 * envelope .* carrier; end

这里有几个细节值得注意。sin²包络的两端被强制置零，这是为了保证激光脉冲在时域上是紧支撑的，避免因为直接截断造成频谱上出现非物理的高频分量。载波用cos而不是sin是有讲究的——两者只差一个π/2的相位，也就是CEP不同。在少周期激光脉冲条件下，CEP会显著影响谐波谱的奇偶次分量比例和截止附近的精细结构。如果你发现跑出来的谱和文献对不上，第一件事就是检查载波用的是cos还是sin。

参数E0的单位是原子单位制下的电场强度。1 a.u.的电场强度约等于5.14e9 V/cm，实验上常用的1e14 W/cm²激光强度换算过来大约对应0.053 a.u.。这样换算，恰好落在氢原子强场电离的典型区间，这也是为什么代码的默认参数靠谱的原因——在物理上是对的。你可以用这个换算关系来调参，比如想模拟1.5e14 W/cm²就把E0设成0.053×sqrt(1.5)≈0.065。

3.2 波函数演化：split-operator时间步进的MATLAB实现

function psi = propagate(psi0, V, dt, nsteps, dx, mask) % 用split-operator方法推进波函数 % psi0 : 初始波函数（复数组，长度N） % V : 势能数组（对角势，含原子势+激光-偶极耦合项） % dt : 时间步长（原子单位制，通常取0.05~0.2 a.u.） % nsteps: 总步数 % dx : 空间网格间距 % mask : 吸收边界掩蔽数组，每步演化后乘上以吸收边界泄露 N = length(psi0); psi = psi0; % 动量空间波数向量，配合fft使用：fft的输出是0~2pi对应的k % 通过fftshift把它变成对称的负频到正频排列 k = (2*pi/(N*dx)) * (-N/2 : N/2-1); Texp = exp(-1i * 0.5 * k.^2 * dt); % 动能传播子（动量空间） for j = 1:nsteps % 第一半势能步：坐标空间逐点相乘 psi = exp(-1i * 0.5 * V * dt) .* psi; % 动能步：变换到动量空间，乘动能传播子 psi = fftshift(fft(psi)); psi = Texp .* psi; % 变换回坐标空间 psi = ifft(ifftshift(psi)); % 第二半势能步 psi = exp(-1i * 0.5 * V * dt) .* psi; % 吸收边界：每步都乘掩蔽函数，防止边界反射 psi = mask .* psi; end end

这段代码的核心思想是：把演化拆成三步，每一步都在对应的对角表象里做逐点相乘。split-operator的精度是二阶的，误差正比于dt的三次方，所以把dt从0.1减到0.05，精度会提升约8倍，但计算量翻倍。实际跑的时候，dt取0.05~0.1 a.u.之间最合适，再小是纯浪费算力，再大就有发散风险。我一般先按dt=0.1跑通，确认结果合理后再加密时间步。

这里有两个实现细节极容易翻车。第一，fftshift和ifftshift必须配对使用。MATLAB的fft默认期望零频率在第一个元素，而k向量是从负到正的对称排列，必须先shift对齐才能逐点相乘，乘完再shift回来。搞反了或者漏了shift，动能传播子会作用在错误的频率成分上，直接导致波函数畸变。第二，势能项和动能项的相位都是作用在复数波函数上的，所以psi全程都是复数组，不能用real()取实部。

还有一个性能优化点：V在激光场中其实是随时间变化的，因为激光-偶极耦合项-x·E(t)的E(t)随时间震荡。常见做法是提前把V(t)堆叠成一个nsteps×N的大矩阵，在propagate外面一次性构建好，函数内部每步取一行来用。这样能省掉大量重复构建势能的时间。如果内存不够，退而求其次每步调用一次laser_pulse计算当前时刻的电场，但要接受大约30%的性能损失。

3.3 偶极矩提取与谐波谱计算

function [spectrum, freq] = extract_hhg(dipole, dt) % 从偶极矩时间序列提取谐波谱 % dipole: 偶极矩期望值序列，长度nsteps % dt : 时间步长 % 返回 : 频谱强度与频率向量 % 对偶极矩做二阶时间导数（中心差分），得到加速度 % 谐波辐射功率正比于偶极矩加速度的模平方，而不是偶极矩本身 acc = zeros(size(dipole)); acc(2:end-1) = (dipole(3:end) - 2*dipole(2:end-1) + dipole(1:end-2)) / dt^2; % 加汉宁窗抑制时间窗截断带来的频谱泄漏 N = length(acc); w = 0.5 - 0.5*cos(2*pi*(0:N-1)/(N-1)); signal = acc .* w; % 傅里叶变换并取模平方 f = fft(signal); freq = (0:floor(N/2)-1) / (N*dt); % 只取正频率部分 spectrum = abs(f(1:floor(N/2))).^2; end

谐波谱的正确算法是：先在每步演化后计算偶极矩期望值d(t) = ∫ψ*(x,t)·x·ψ(x,t)dx，然后对时间序列求二阶导得到偶极加速度，再对加速度做傅里叶变换。为什么不用偶极矩直接做FFT？因为偶极矩本身包含一个随激光场线性振荡的缓变成分，直接变换得到的是巨幅低频峰，谐波区的弱信号全被淹没了，根本看不清平台区。加速度相当于对低频成分做了二阶高通滤波，谐波信号才能浮现出来。

汉宁窗的作用是抑制时间窗截断导致的频谱泄漏。仿真时间窗有限，等于是把无限长的信号乘了一个矩形窗，频谱上会叠sinc函数的边瓣，把相近的谐波峰糊在一起。加汉宁窗能把边瓣压低到约-31dB，谐波峰干净很多。如果你的谐波谱上每个主峰两侧都有可疑的小鼓包，先怀疑是不是忘了加窗。

3.4 基态初始化：虚时间演化的妙用

function psi0 = init_wavefunction(V_atom, dx, N) % 用虚时间演化求原子势下的基态波函数 % 原理：把演化算符中的时间t换成虚数τ = -i*t， % 较高的能量本征态会以更快的速率指数衰减，最终只剩基态 psi = exp(-( (1:N) - N/2 ).^2 / (2*10^2) )'; % 初始高斯猜测 psi = psi / norm(psi); % 归一化 tau = 0.1; % 虚时间步长 n_iter = 500; % 迭代步数 for k = 1:n_iter % 一步虚时间演化：exp(-H*tau)作用于波函数 psi = ifft( fft(psi) .* exp(-0.5*(2*pi*((0:N-1)-N/2)/(N*dx)).^2 * tau) ); psi = exp(-V_atom * tau) .* psi; psi = psi / norm(psi); % 每步归一化，防止整体衰减为零 end psi0 = psi; end

虚时间演化的原理很巧妙：把薛定谔方程中的t换成纯虚数τ = -i·t，原来的振荡相位演化exp(-i·E·t)变成了指数衰减exp(-E·τ)。能量越高的本征态衰减越快，经过足够多步虚时间演化后，所有激发态成分都被压制，只剩下能量最低的基态。每步的归一化是必须的，否则波函数模方会指数衰减到零，整个算法就失效了。

初始化基态的精度直接影响后续物理过程的正确性。如果基态能量和理论值误差超过1%，谐波谱平台区的干涉结构会出现不规则的扭曲。这个函数的收敛判据很简单：算完后再用哈密顿量作用一次，看期望能量和前一步差多少，如果相对变化小于1e-6就可以停。实际跑的时候500步虚时间通常足够，但如果你换了更重的原子（势阱更深），可能需要增加到1000步。

4. 运行与调参：从零开始跑通一个真实案例

4.1 环境确认与三步启动法

这份代码是纯MATLAB脚本，不涉及Simulink、App Designer等组件。理论上R2016b之后的版本都能跑，因为用到的核心函数只有fft、ifft、trapz、abs等基础函数，连Signal Processing Toolbox都不需要。我第一次跑的时候用的是R2020a，整个过程没有任何工具箱依赖报错。

开始之前把环境确认和目录安排好：

1. 解压到不含空格的路径：D:\HHG\simulation（中文路径和空格偶尔会在fft线程池调度上出问题，避开最省心） 2. 在MATLAB里cd到该目录，先运行parameters.m，工作区出现params结构体 3. 运行main.m，观察命令行输出的进度信息和弹出的两个图形窗口

如果一切正常，你会看到两个图：一个是波函数模方分布的时间演变（能看到波包在激光电场作用下振荡、部分逃离原子势阱），另一个是谐波谱曲线（横轴为谐波阶次，纵轴为强度，用log坐标表示）。第一次跑通后，重点不是看图有多漂亮，而是验证一个关键特征——谐波谱在低阶区快速下降后进入一个平台区，然后在某个阶次附近急剧截止。这个截止位置就是你要和半经典公式对照的数字。

4.2 核心参数对照表

第一次跑通之后，最重要的验证动作是核对谐波谱的截止位置。如果你看到的谱是全频谱强度单调下降、没有平台区，先检查E0是否太低——低于0.02 a.u.时隧穿电离概率很低，谐波产生效率大幅下降，平台区就不明显。反过来E0高于0.1 a.u.时，电子可能在电场峰值前就已经完全电离，后续的复合过程几乎没有，谐波也会被抑制。找到合适的E0区间，是调参的第一步。

4.3 单参数扫描：改变激光强度看截止频率移动

% 参数扫描脚本：scan_E0.m % 目标：对比不同峰值场强下的谐波截止位置 E0_list = [0.04, 0.06, 0.08, 0.10]; % 四组场强 cutoffs = zeros(size(E0_list)); figure; for idx = 1:length(E0_list) params.E0 = E0_list(idx); % 基态不依赖激光场，只依赖原子势，所以只需要初始化一次 psi0 = init_wavefunction(params.V_atom, params.dx, params.N); % 构建含激光场的完整势能矩阵并推进 V_total = build_effective_potential(params); % 原子势+偶极耦合项 psi = propagate(psi0, V_total, params.dt, params.nsteps, params.dx, params.mask); % 提取偶极矩序列（实际代码中在propagate内同步记录） dipole = record_dipole_moment(psi, params); [spectrum, freq] = extract_hhg(dipole, params.dt); % 半经典截止公式：N_cutoff = (Ip + 3.17*Up) / omega0 Up = params.E0^2 / (4*params.omega0^2); % 有质动力能 Nc = (params.Ip + 3.17*Up) / params.omega0; cutoffs(idx) = Nc; % 半对数坐标绘制谐波谱，横轴用谐波阶次 loglog(freq/params.omega0, spectrum); hold on; end legend(arrayfun(@(x) sprintf('E0=%.2f', x), E0_list, 'UniformOutput', false)); xlabel('Harmonic Order'); ylabel('Intensity (a.u.)');

跑参数扫描时有个关键的性能提示：如果不做任何优化，一组参数（nsteps=4000，网格N=1024）在普通桌面上约需1~2分钟，四组就是接近10分钟甚至更久。实际的快慢强烈依赖CPU单核性能和MATLAB的fftw线程配置。建议先把一组参数跑完并记录tic/toc耗时，再决定是否缩小网格或减少光学周期数量。如果每组超过5分钟，先把N减到512，dx增加一倍，快速看趋势，锁定最优区域后再用高分辨率精跑。

另外注意扫描脚本中基态只初始化了一次——这是合理的，因为基态波函数只依赖原子势，和激光场强无关。如果你不小心把init_wavefunction放进了循环里，计算量会凭空多出4倍，而且没有任何收益。这也是调试脚本时最常看到的浪费。

5. 避坑指南：数值发散、边界反射与性能瓶颈的排查记录

5.1 波函数NaN：时间步长与空间网格不匹配

现象：波函数模方在某个时间步突然出现尖峰，数值变成NaN或Inf，谐波谱直接崩掉，命令行抛错。

原因：dt和dx不匹配。split-operator方法的稳定性依赖动量空间最高频率与相位增量的乘积不过大。动量空间最高波数是π/dx，动能传播子在这个波数上的相位是π²·dt/(2·dx²)，这个值一旦超过视觉安全阈值，高频谱成分会被过度相移，叠加出非物理振荡，最终数值爆炸。

解决：先把dx算清楚。如果你把计算区域L设成500 a.u.，网格N还是1024，那么dx约等于0.98 a.u.，这时候dt必须压到0.02 a.u.才有救。我一般先固定L和N，从dt=0.1开始每跑100步检查一次max(abs(psi))，出现异常就减半，用这个方式快速找到临界dt。经验上dx在0.3~0.5之间时dt取0.05~0.1最稳定，两种情况必须同时收严，只收一个没用。

5.2 边界反射伪峰：吸收边界宽度不够

现象：波函数在边界附近出现驻波条纹，谐波谱背景上出现一系列等间距的伪峰，波长还和激光波长有关。

原因：周期性边界条件让泄露的波包从另一侧穿回来，与真实物理波包发生干涉。吸收边界mask如果宽度不够，电子在每一步时间内能穿透mask抵达真实边界，下一步就从另一侧冒出来了。这不是物理效应，纯粹是数值边界处理不当。

解决：mask用cos²窗口，宽度取整个计算区域的10%~20%，而不是随便给一个窄过渡带。同时增加一个预处理步骤——在激光到来前先让波函数虚时间演化几百步，把被mask扭曲的基态残余消耗掉。我踩过一次这个坑：mask宽度只设了5%，跑出来的谐波谱低频区有一堆可疑的等间距振荡，排查了大半天才发现是边界问题而非物理信号。从那以后mask宽度低于10%我都不开跑。

5.3 性能瓶颈：多线程反而不快

现象：运行时间比预期多出好几倍，CPU占用率却上不去，运行过程中偶发小停顿。

原因：MATLAB的fft在多线程下会自适应选择变换算法，但网格N=1024这种小数组的开销里，线程调度和内存同步反而占了相当比例，多线程收益趋近于零甚至为负。另一个常见瓶颈是每步都重新构建包含cos计算的V数组——4000步乘以数组长度N的cos计算量不可小觑。

解决：小网格仿真用maxNumCompThreads(1)强制单线程，这是最快见效的手段。大网格（N=4096以上）则相反，默认多线程就有收益，不要手动限制。关于势能矩阵，提前分配一个nsteps×N的大矩阵一次性构建，propagate内部逐行取出，能省约30%的总时间。内存不够时退而求其次，在propagate里按步调一个轻量的激光电场函数，但要把重复计算的开销还回去。

5.4 换原子后平台区消失：电离能变了，参数没跟上

现象：把氢原子换成氩原子类似的参数之后，谐波谱的平台区明显变窄，甚至完全消失，只剩下快速衰减的少量低阶谐波。

原因：不同原子的电离能Ip不同，同样的场强下隧穿电离率不同，波函数空间延展也不同。氩的电离能约0.58 a.u.比氢的0.5高，但关键在于直接用氢的E0时，Keldysh参数γ=omega0·sqrt(2·Ip)/E0大于1，进入了多光子电离区，而非隧穿区。多光子区的谐波产率低得可怜，平台区自然看不到。

解决：换原子前先用Keldysh参数算一下当前场强是否在隧穿区域。γ<1才是隧穿电离主导，谐波截止才会偏离半经典公式。理想区间是0.3~0.6。如果γ偏大，就把E0往上调或者把omega0调小（用更长的波长）。这个检查在换原子时是强制动作，不查就直接跑，十有八九得到一条没法解释的谱线。

5.5 谐波峰形不对称：时间窗截断了物理过程

现象：谐波谱的峰形左右不对称，平台区出现起伏的调制条纹，和文献里的光滑曲线差异明显。

原因：激光脉冲只有8个光学周期，整个谐波产生过程被限制在极短的时间窗内。如果偶极矩时间序列在窗的两端没有被衰减到接近零，FFT的矩形窗等效作用会产生严重的频谱泄漏，sinc函数边瓣叠加在真实信号上。

解决：确保偶极矩序列在脉冲开始前和结束后都有至少1~2个光学周期长的“静默区”，让物理信号在窗两端自然衰减到近零。另一个有效手段是加窗函数宽度更缓的平坦窗，比如Hann窗的变体或Nuttall窗，边瓣抑制更好。检查手段很简单：取偶极矩序列看两端值的大小，如果和峰值只差一个数量级，那频谱泄漏一定在污染你的结果。

6. 二次开发：从复现到改造的实用思路与验证方法

6.1 换原子种类：调整软核势与电离能

在init_wavefunction.m里把原子势从纯库仑势换成软核势V(r) = -Z/sqrt(r² + a²)，其中a是软核参数，用于避免r=0处的奇点导致数值不稳定。氢原子取Z=1、a≈0.5，基态能量恰好落在-0.5 a.u.附近。换氦原子时改成Z=1.34、a≈0.9，基态能量约-0.9 a.u.。改完这两个参数后，main.m里parameters.Ip也要同步更新，谐波谱的截止公式里要用新Ip值。

换原子后的验证方法是评估基态能量误差。在虚时间演化收敛后，计算哈密顿量期望值E_exp = sum(conj(psi).(Hpsi))*dx，和理论值对比，误差小于1%就算合格。如果误差偏大，先缩短虚时间步长τ，再看收敛曲线是否单调下降。这个步骤能避免把错误的基态带入后续动力学计算，是换原子的第一道质量关卡。

6.2 加椭偏光场：从一维到二维的改动思路

如果你想把线偏振光改成椭偏光，需要在laser_pulse.m里返回两个电场分量Ex和Ey，分别对应两个偏振方向的电场。对应的偶极耦合项变成-x·Ex(t) - y·Ey(t)，波函数必须从一维升级到二维，动力学方程里加一个y方向的空间导数。此时动能传播子要做二维FFT，也就是fft2/ifft2，波数向量也要相应变成两个方向的网格。

二维计算的代价增长很猛。256×256的网格在普通机器上勉强能跑，但时间步数需要按比例减少，否则计算量爆炸。我试过一次，算到一半发现风扇狂转，果断把光学周期数减半，才在可接受的时间内拿到结果。建议第一次改造成椭偏光时，先用64×64的粗网格跑通流程验证物理行为，再逐步加密，不要一上来就追求分辨率。

6.3 普适验证手段：概率守恒检验

这是我认为整个二次开发周期里性价比最高的检查动作。在无吸收理想情形下，波函数模方的总积分应该是1，边界吸收会带来概率流失，但流失应该在合理范围内。实现方式是：在propagate的循环里每100步计算一次sum(abs(psi).^2)*dx，记录整个演化过程中的变化。

% 在main.m中添加概率守恒检查 prob = zeros(1, params.nsteps/100 + 1); for j = 1:params.nsteps % 正常propagate演化 psi = one_step(psi, V(j,:), dt, dx, mask); if mod(j, 100) == 0 prob(j/100 + 1) = sum(abs(psi).^2) * params.dx; end end % 如果总概率衰减超过5%，说明边界吸收过强或物理区域不够大 if prob(end) / prob(1) < 0.95 warning('Probability loss exceeds 5%%: check mask width and L'); end

如果总概率衰减超过5%，说明吸收边界过强把物理信号也吃掉了，或者计算区域不够大导致真实波包抵达边界。这个检验在改任何参数之后都强制跑一遍，比看谐波谱暴露问题快得多——概率异常比谱形异常早发现至少几百个时间步。另外，概率守恒检验也是检查数值发散预警信号的好方法，发散总是先让概率变为异常值，再在谐波谱上体现出来。

从那以后我每次跑这类仿真，都会把“概率守恒检查”写进脚本的第一步。改完任何参数，先跑50步看概率有没有掉，确认稳定再跑完整流程，这个习惯帮我挡掉过至少三次因参数改动引入的隐性错误。希望帮到你——如果你第一次跑通这份代码，能在截止频率那条线上复现出半经典公式预测的位置，说明你的每一步踩对了。

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