简介：面向激光器原理学习与数值仿真的入门者，这份主动调Q固体激光器MATLAB仿真资料围绕四能级系统、粒子数反转、声光调Q及速率方程建模展开，可用于理解脉冲形成、参数敏感性分析与结果解读。压缩包为rar格式，共2个文件，均以.m脚本形式提供，包体约919B，便于直接打开运行或按需修改参数。已有2102人学习，说明其在相关课程实验与自学场景中有一定参考价值。读者可借助脚本复现主动调Q过程，调整泵浦速率、增益与损耗等变量观察脉冲宽度、峰值功率和能量效率变化，并参考附带的被动调Q例程对比两种调Q方式的差异，从而加深对激光器设计与优化的理解。

1. 主动调Q固体激光器仿真：先看它能不能替我算出那 20 ns

调Q脉冲到底能压到多窄、峰值功率够不够烧穿一片靶材、重复频率提上去之后单脉冲能量掉多少——这些问题在画结构图之前就该有个数。这套主动调Q固体激光器matlab仿真文件做的事很直接：把腔内光子数密度和反转粒子数密度的耦合方程丢给求解器，跑出纳秒级的巨脉冲波形，再给出脉宽、峰值功率、单脉冲能量三个数。它适合三类人：做激光器方案前期评估、需要在几种腔型和泵浦功率之间快速取舍的；课程设计或论文里要交一张像样的脉冲波形的；以及实验台上调不出脉冲、想先在线把参数走一遍再回实验室的。代码本身不长，值钱的是参数表里每个数字的物理口径——口径错了，曲线画得再漂亮也是假的。

2. 速率方程怎么落：两个物理量和一张参数表撑起整个模型

2.1 主动调Q的时域节奏：几百微秒储能，几十纳秒放完

主动调Q的过程拆成三段看最清楚。第一段是泵浦储能期，Q开关保持在关断状态，腔内插入损耗被人为抬高，阈值反转粒子数密度被顶到一个很高的位置，谐振腔“想起来也起不来”，上能级只能老老实实攒粒子数。这一段的时间尺度由上能级寿命和泵浦速率决定，常见在几十到几百微秒量级。

第二段是开关打开后的瞬态。腔内损耗在几个纳秒内塌下去，阈值跟着塌，而此时的反转粒子数密度是阈值的几倍，净增益为正且远大于损耗，腔内自发辐射在往返过程中被反复放大，光子数密度呈指数爬升。第三段是增益饱和与能量泄放，光子数密度涨到一定程度后，受激辐射消耗反转粒子的速度超过泵浦补充速度，增益被拉回阈值以下，光子数密度冲高回落，形成一个宽度几十纳秒的巨脉冲。

这三段的时间尺度差了三到四个数量级。储能段几百微秒，脉冲段几十纳秒，同一个常微分方程组里既有慢变过程又有快变过程，这就是刚性方程最典型的特征，也直接决定了后面求解器不能随便选。很多人第一次跑这个仿真，卡就卡在这里：拿默认的 ode45 一路跑下去，要么跑到天荒地老，要么步长被压到极小之后报步长不足。

2.2 状态变量选两个就够：光子数密度与反转粒子数密度

速率方程的写法有很多版本，有的把上下能级粒子数分开写，有的加上下能级吸收项做准三能级处理。对于主动调Q这个场景，只要不是研究极低掺杂或者高温工况下的准三能级效应，四能级近似的两变量版本足够覆盖工程判断，也是这类仿真文件里最常见的结构。

两个状态量分别是腔内平均光子数密度 φ（单位 1/m³）和增益介质内的反转粒子数密度 N（单位 1/m³）。前者描述腔里有多少光子，后者描述还能放出多少能量。方程写作：

dφ/dt = c·σ_eff·N·φ - φ/τc + β·N/τf dN/dt = Rp - N/τf - c·σ_eff·N·φ

第一项是受激放大，第二项是腔损耗，第三项是自发辐射往激光模式里的种子贡献。反转粒子的方程里，泵浦项往里灌，自发辐射衰减往外漏，受激辐射消耗往下掉。三条腿，少一条模型就失真。

这里有个必须提前讲清楚的口径问题。增益只在长度为 l 的介质里产生，而光子在长度为 L 的腔里往返，严格推导要把 l/L 这个填充因子带上。教学代码里常见的处理是把截面折算成有效截面 σ_eff = σ·l/L，相当于把填充因子吸收进参数里。好处是两个方程的受激项系数完全一致，能量账天然守恒——每产生一个光子就精确消耗一个反转粒子。代价是 σ_eff 不再是一个纯材料参数，改腔长时它要跟着重算。这一点如果不写清楚，后面做腔长扫描就会得到一组没法解释的曲线。

提示：如果你的模型需要单独研究填充因子或者非均匀增益分布，就不要用折算截面，老老实实把 l/L 显式写进受激项，并接受非守恒带来的误差；工程判断场景下折算写法够用且更稳。

2.3 一张参数表：把器件手册上的数字翻译成模型能吃的口径

参数表是这个仿真文件里最该被认真对待的部分。代码可以照抄，参数抄错了整条曲线就是废的。下面这张表是一套能跑出毫焦级脉冲的典型口径，读者可以按自己手上的器件替换，但单位必须统一到 SI。

两个阈值可以当场算出来做校核。基础阈值反转粒子数密度 N_th = 1/(c·σ_eff·τc)，代进去约 9.4e23 /m³；关断状态下的阈值约 7.4e24 /m³。泵浦 200 μs 后实际攒到的反转粒子数密度 N = Rp·τf·(1-exp(-t/τf))，约 3.3e24 /m³，低于关断阈值、高于开门阈值，这正是想要的工作点：关着的时候不出脉冲，一开就出。

如果算出来 N 比关断阈值还高，说明泵浦太猛或者开关损耗不够，仿真会告诉你“关断期就漏光了”，这时候要么降泵浦、要么加大 δQ。这一步手算花两分钟，能省掉后面半小时的排错。

2.4 Q开关的数学表达：阶跃、双曲正切还是实测驱动波形

主动调Q的“主动”体现在开关是外部信号驱动的，模型里必须有一个随时间变化的损耗项。最简单的写法是理想阶跃：t < t_open 时损耗为 δ0+δQ，t ≥ t_open 时降到 δ0。写起来一行，跑出来脉冲前沿会偏陡，脉宽偏窄，和实测对不上。

第二种是双曲正切过渡，用 gQ(t) = 0.5·(1 + tanh((t - t_open)/t_rise)) 表示开关的开启程度，损耗写成 δ(t) = δ0 + δQ·(1 - gQ(t))。t_rise 对应驱动电路的上升沿，声光Q开关常见在几十到几百纳秒，电光Q开关可以做到几纳秒。这个写法平滑、可导，求解器不会在阶跃点附近反复缩步长，是这类仿真文件里最常用的形式。

第三种是直接把示波器采到的驱动波形插值进来，用 interp1 或 timeseries 喂给求解器。适合已经做完实验、想反过来标定模型参数的场景。缺点是对波形质量敏感，噪声会让求解器频繁触发事件检测。

我一般先用 tanh 版本跑通，把 t_rise 设成器件手册上的典型值，等曲线和预期对上了，再换实测波形做精细比对。上来就插值实测数据的，十个里有八个会因为波形毛刺把 ode15s 逼到 MaxStep 下限，然后误以为是模型本身有问题。

3. 让 ode15s 跑出波形：分段求解、开关函数与后处理脚本

3.1 参数脚本：把所有物理量收进一个结构体

仿真文件最容易腐化的地方是参数散落在十几个脚本里，改一个泵浦功率要翻三个文件。我的习惯是单独一个参数文件，所有量挂在同一个结构体上，主脚本只认这个结构体，不认工作区里的裸变量。

%% params_laser.m —— 全部物理量集中在这里，改参数只动这一个文件 p.c = 3e8; % 光速 m/s p.sigma = 2.8e-23; % 受激发射截面 m^2 p.l = 5e-3; % 增益介质长度 m p.L = 0.30; % 腔光学长度 m p.sigma_eff = p.sigma * p.l / p.L; % 折算到腔长的有效截面 p.tau_f = 230e-6; % 上能级寿命 s p.R = 0.80; % 输出镜反射率 p.delta0 = -0.5*log(p.R) + 0.02; % 基础往返损耗（含散射） p.deltaQ = 0.90; % Q开关关断时的附加损耗 p.tau_c = p.L / (p.c * p.delta0); % 基础光子寿命 s p.Rp = 2.5e28; % 泵浦速率密度 1/(m^3*s) p.beta = 1e-6; % 自发辐射耦合因子 p.phi0 = 1e10; % 初始光子数密度 1/m^3 p.t_open = 200e-6; % 开关打开时刻 s p.t_rise = 10e-9; % 开关过渡时间 s p.t_span = 500e-9; % 脉冲段观察窗口 s p.A_mode = 1e-6; % 模场面积 m^2

逻辑上就是一次“把物理量翻译成代码量”的动作。σ 和 l、L 一起折算成 σ_eff，后面方程里不再出现长度比；τc 由损损耗反算，不手工填数字，避免改反射率忘了同步改寿命。δQ 取 0.9 是一个能让关断阈值抬到储能之上的常用量级，实际值取决于开关的衍射效率和腔内光束尺寸，声光器件通常在 0.5 到 1.5 之间。

p.phi0给 1e10 是数值上的种子，避免零解陷阱——方程在 φ=0 处是平衡点，初值给零的话求解器会告诉你“一切正常”，然后给你一条死线。

3.2 速率方程函数：开关函数写在一个地方

function dydt = rate_eq(t, y, p) % y(1): 腔内光子数密度 phi [1/m^3] % y(2): 反转粒子数密度 N [1/m^3] phi = max(y(1), 0); N = max(y(2), 0); % Q开关开启程度：0 为关断（高损耗），1 为全开（低损耗） gQ = 0.5 * (1 + tanh((t - p.t_open) / p.t_rise)); % 时变损耗与时变光子寿命 delta = p.delta0 + p.deltaQ * (1 - gQ); tau_c = p.L / (p.c * delta); gain = p.c * p.sigma_eff * N; % 增益系数 1/s dphi = gain * phi - phi / tau_c + p.beta * N / p.tau_f; dN = p.Rp - N / p.tau_f - gain * phi; dydt = [dphi; dN]; end

几处值得说。max(...,0)是数值保护，求解器在自适应步长时会试探性地把状态量推到负值，截断能避免 NaN 沿链条扩散，代价是极小概率掩盖真实的不物理振荡——如果发现截断被频繁触发，说明步长设太大了。

开关函数只在这一处出现，损耗和光子寿命都由它派生。这样改开关模型时只改一行，不会出现“损耗改了但寿命没改”的错配。gain在两个方程里用同一个符号，保证光子产生和反转消耗严格对应，这是前面折算截面写法换来的好处。

p.beta * N / p.tau_f是自发辐射种子项，量级很小但在脉冲起振的最初几十皮秒里决定谁先起来。把它删掉也不会报错，只是脉冲建立时间会往后挪，脉宽略微变窄。

3.3 分段求解：泵浦段走解析，脉冲段交给 ode15s

刚性方程的经典翻车点是把两百微秒的储能过程和二十纳秒的脉冲放在同一个求解区间里跑。储能段解析解是现成的，没必要让求解器去啃，直接算完当初始条件喂给第二段，又快又稳。

%% solve_pulse.m —— 泵浦段解析、脉冲段数值 run('params_laser.m'); % ---- 第一段：0 ~ t_open 泵浦储能，无激光输出 ---- t_pump = linspace(0, p.t_open, 400); N_pump = p.Rp * p.tau_f * (1 - exp(-t_pump / p.tau_f)); phi_pump = zeros(size(t_pump)); N_end = N_pump(end); % 开门瞬间的反转粒子数密度 % ---- 第二段：开关打开后的巨脉冲 ---- y0 = [p.phi0; N_end]; opts = odeset('RelTol', 1e-8, ... 'AbsTol', [1e6 1e16], ... 'MaxStep', 1e-10); [t_p, y_p] = ode15s(@(t,y) rate_eq(t + p.t_open, y, p), ... [0 p.t_span], y0, opts); t_all = [t_pump(:); t_p + p.t_open]; phi_all = [phi_pump(:); y_p(:,1)]; N_all = [N_pump(:); y_p(:,2)];

第一段用解析解，前提是关断期确实没有激光输出，也就是自发辐射放大不到阈值。前面算过 N_end ≈ 3.3e24 低于关断阈值 7.4e24，这个前提成立。如果你的参数下 N_end 超过关断阈值，解析解就不适用了，那时候必须整段数值求解，或者至少把第一段也交给 ode15s。

第二段用 ode15s，MaxStep卡在 1e-10（0.1 ns），对应 τc 的百分之一量级，够解析脉冲前沿。AbsTol给两个分量不同的量级，光子数密度 1e6、反转粒子数密度 1e16，这是按各自稳态值的相对精度折算的，两个分量量级差十几个数量级时统一 AbsTol 会让小的那个被过度求解。

匿名函数里写的是t + p.t_open，把相对时间平移回绝对时间。这一行是血泪经验：开关函数里比较的是绝对时刻，如果忘了加偏移，开关会在 t=0 就打开，你会得到一条在储能段就出光的荒谬曲线，而且怎么看都看不出哪里错了。

3.4 后处理：脉宽、峰值功率、单脉冲能量三个数怎么取

%% post_process.m h_nu = 1.87e-19; % 1064 nm 光子能量 J V_cav = p.A_mode * p.L; % 腔内模式体积 m^3 delta_o = -0.5*log(p.R); % 输出耦合损耗 tau_co = p.L / (p.c * delta_o); % 输出耦合对应的时间常数 s t_ns = t_p * 1e9; P = y_p(:,1) * V_cav * h_nu / tau_co; % 输出功率 W [Pmax, i_pk] = max(P); idx = find(P > 0.5*Pmax); fwhm = t_ns(idx(end)) - t_ns(idx(1)); % 脉宽 ns E = trapz(t_p, P); % 单脉冲能量 J fprintf('峰值功率 %.1f kW\n', Pmax/1e3); fprintf('脉宽 %.2f ns\n', fwhm); fprintf('单脉冲能量 %.3f mJ\n', E*1e3);

输出功率的换算走的是“腔内光子数 × 光子能量 / 输出耦合时间常数”。τ_co 用输出耦合损耗单独算，而不是用总损耗 τc，因为只有从输出镜出去的那部分才算有效输出，腔内散射和开关吸收不算。这个细节在做效率评估时会差出百分之十几。

脉宽用半高全宽，find取第一个和最后一个超过半高的点，对单峰脉冲没问题。如果扫描参数时出现了多峰或者脉冲分裂，这个写法会把几个峰一起框进去，得先做峰值检测再分段算。

单脉冲能量用 trapz 对功率积分，和后面第 6 章的储能法互相校验，两条路算出来的数应该只差几个百分点。差得离谱，一定是单位或者体积口径出了问题。

4. 参数扫描怎么扫：重复频率、泵浦功率、开关时间的三条曲线

4.1 重复频率扫描：把单脉冲脚本包成函数

单次求解只告诉你一个工作点，实际选型要看趋势。第一步是把整个求解流程封成一个函数，输入参数结构体，输出四个标量：单脉冲能量、峰值功率、脉宽、开门时的初始反转粒子数密度。

function [E, Ppk, tau, N0] = one_pulse(p) % 输入 p 为参数结构体，返回单脉冲能量(J)、峰值功率(W)、脉宽(ns)、初始反转密度 T = 1 / p.f_rep; % 泵浦周期 s N0 = p.Rp * p.tau_f * (1 - exp(-T / p.tau_f));% 周期末储能 y0 = [p.phi0; N0]; opts = odeset('RelTol',1e-8,'AbsTol',[1e6 1e16],'MaxStep',1e-10); [t, y] = ode15s(@(t,y) rate_eq(t + p.t_open, y, p), [0 p.t_span], y0, opts); h_nu = 1.87e-19; V_cav = p.A_mode * p.L; tau_o = p.L / (p.c * (-0.5*log(p.R))); P = y(:,1) * V_cav * h_nu / tau_o; [Ppk, ~] = max(P); Ppk = Ppk; idx = find(P > 0.5*max(P)); tau = (t(idx(end)) - t(idx(1))) * 1e9; E = trapz(t, P); end

关键改动是把储能公式里的时间从固定的 t_open 换成泵浦周期 T。重复频率越高，周期越短，两次脉冲之间攒回来的反转粒子数越少，单脉冲能量随之下降。这个下降不是线性的，当 T 接近 τf 时开始明显，T 远小于 τf 时能量基本正比于 T。

%% scan_reprate.m f_list = [1 2 5 10 20 50] * 1e3; % 1 kHz 到 50 kHz res = zeros(numel(f_list), 4); for k = 1:numel(f_list) p.f_rep = f_list(k); p.t_open = 1 / f_list(k) * 0.98; % 留 2% 余量给脉冲建立 [res(k,1), res(k,2), res(k,3), res(k,4)] = one_pulse(p); end

t_open取周期的 98% 是个工程习惯，留一点余量让脉冲在下一个泵浦周期开始前完全泄放完，否则残余光子会带进下一轮，扫描曲线尾部会出现不自然的抖动。

跑完这条曲线通常能看到三段：低频段能量饱和在一个平台上，中频段线性下降，高频段掉得比线性还快，因为开门时的反转粒子数密度已经接近阈值，增益开关效应被削弱。这条曲线的拐点位置直接决定你这台器件能跑多高的重频还保持可用的单脉冲能量。

4.2 泵浦功率扫描：阈值附近别用线性外推

泵浦速率 Rp 和泵浦功率近似成正比，扫 Rp 相当于扫泵浦电流。这条曲线上最值得看的是阈值附近的行为：Rp 低于某个值时开门后根本不出脉冲，因为 N0 没到阈值；刚过阈值时脉宽很宽、峰值功率很低；继续加大泵浦，能量快速上升然后进入饱和。

%% scan_pump.m Rp_list = linspace(0.6, 2.0, 15) * 2.5e28; for k = 1:numel(Rp_list) p.Rp = Rp_list(k); [E(k), Ppk(k), tau(k), N0(k)] = one_pulse(p); end

扫描时把 N0 一起存下来很有用。把 N0 和阈值 N_th = 1/(c·σ_eff·τc) 画在同一张图上，能直观看出哪些泵浦点工作在“刚过阈值”的脆弱区。这个区域里脉宽对参数极其敏感，仿真结果的可信度也最低，因为实际器件的自发辐射种子、模式竞争、热效应都会把结果推到另一条路上。工程上一般不把工作点选在这里。

另外提醒一句，Rp 和泵浦电功率之间还隔着一个量子效率和一个模式重叠因子，仿真里给的 Rp 是已经折算过的“有效泵浦速率”，别拿电功率直接除以光子能量去对，那样算出来的 Rp 会偏大好几倍，曲线整体左移，阈值对不上。

4.3 开关过渡时间：驱动电路的上升沿值多少钱

t_rise 这个参数在模型里只占一行，在实物上对应的是驱动电路的上升沿和开关的响应速度。扫描它的代价极低，收益却很直接。

%% scan_rise.m rise_list = [2 5 10 20 50 100] * 1e-9; for k = 1:numel(rise_list) p.t_rise = rise_list(k); [E(k), Ppk(k), tau(k), ~] = one_pulse(p); end

典型结果是这样：t_rise 从 2 ns 加到 50 ns，脉宽会明显变宽，峰值功率下降，但单脉冲能量基本不变——能量由储能决定，开关只决定这些能量在多长时间里放出来。这个结论对选型很有指导意义：如果指标卡的是峰值功率，那驱动电路的上升沿值得花钱；如果卡的是单脉冲能量，上升沿够用就行，钱花在泵浦上更划算。

再往下压到几纳秒，会看到脉宽不再明显变窄，因为此时脉宽已经被腔光子寿命和增益饱和的动力学限制住了，开关再快也没用。这个“饱和点”是这个模型给出的最有价值的信息之一，它告诉你这套腔型和增益介质的脉宽下限在哪。

4.4 扫描结果的整理：别把结构体数组堆在工作区

十几组参数扫下来，工作区里会堆满匿名变量。建议每条扫描曲线直接落成一个表格文件，字段名写清楚，后面做图或者写报告直接读表。

T = table(f_list(:)/1e3, res(:,1)*1e3, res(:,2)/1e3, res(:,3), ... 'VariableNames', {'RepRate_kHz','Energy_mJ','PeakPower_kW','PulseWidth_ns'}); writetable(T, 'scan_reprate.csv');

字段名带上单位，是给自己留的后悔药。三个月后回头翻数据，谁也不记得那一列到底是 0.003 秒还是 3 毫秒。

5. 排查清单：脉冲消失、能量翻倍、波形拖尾的五种典型翻车

5.1 曲线是一条死线，从头到尾没有脉冲

现象：主脚本跑完不报错，画出来的光子数密度是一条贴着零的直线，或者只有极缓慢的爬升，看不到任何尖峰。

原因通常有三个。一是开门时的反转粒子数密度 N0 低于开门后的阈值 N_th，增益压根没起来，这种情况在泵浦速率填小了或者泵浦时间填短了的时候最常见。二是开关函数写反了，tanh 的符号一颠倒，开关就在 t_open 之前打开、之后反而关断，腔损耗始终没降下来。三是初始光子数密度给了 0，方程在零解处稳定，求解器认为“一直不动就是正确解”。

解决：先手算一次 N0 和 N_th 做对比，这两个数对不上就不要往下跑。再把 gQ 单独画出来看时间轴对不对，一条从 0 平滑升到 1 的曲线是正常的。最后确认 phi0 给了非零种子，量级 1e8 到 1e12 都可以，差别只在脉冲建立时间上。

5.2 ode45 跑到一半卡住，或者报步长不足

现象：换成 ode45 之后，前几秒进度条几乎不动，或者跑出 “Unable to meet integration tolerances without reducing the step size” 的警告，解在中途断掉。

原因就是刚性。储能段的时间常数是 τf 的 230 μs 量级，脉冲段的特征时间是 τc 的 7.6 ns 量级，两者差了三万倍。显式求解器为了保证稳定，步长会被脉冲段逼到纳秒以下，于是一共要走几十万步才能爬完两百微秒的储能段。

解决：换 ode15s 或 ode23s 这类刚性求解器，同时把泵浦段用解析解切出去，只让求解器处理脉冲段。再加一条 MaxStep 上限，防止求解器在脉冲峰值附近一步跨过整个上升沿。三步做完，同样的算例从几分钟降到几秒。

5.3 积分出来的脉冲能量比储能还大

现象：用 trapz 积出来的单脉冲能量是十几个毫焦，手算储能只有一两个毫焦，能量凭空多出来一个量级。

原因基本锁定在单位或者体积口径上。最常见的三种：截面填的是 cm² 没换成 m²，差 1e4

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