Python系统仿真其二：半导体激光器速率方程

速率方程由复光场方程和激光腔内的载流子密度方程组成，是描述系统的动态特性的时变方程。设激光器输出的场振幅和相位分别为 A A A和 ∅ \emptyset ∅，载流子密度为 N N N，参考Ohtsubo的半导体激光器专著，设计相应程序。由于方程比较复杂，因此这里不对方程作过多介绍，感兴趣的读者可以自己去阅读原文。1985年哈肯从两能级原子结构出发，写出带有环形谐振腔的FP激光器的时变电场，并将时变电场、偏振方程与反转粒子数的方程通过线性代换和化简，将其转化成Lorenz方程形式(Lorenz-Haken方程)。这是非常有意思的一种想法，毕竟半导体激光器能够沿着hopf分岔的路径进入混沌。后来经过几年的不断研究，重写了载流子密度和场方程，舍弃了偏振方程才最终得到了如今的半导体激光器速率方程。 d A d t = 1 2 [ g ( N − N 0 ) 1 + ϵ A 2 − τ P − 1 ] ∗ A + k f τ L ∗ A ( t − τ f ) cos ⁡ θ \frac{dA}{dt}=\frac{1}{2}\left[\frac{g\left(N-N_0\right)}{1+\epsilon A^2}-\tau_P^{-1}\right]\ast A+\frac{k_f}{\tau_L}\ast A\left(t-\tau_f\right)\cos{\theta} dtdA​=21​[1+ϵA2g(N−N0​)​−τP−1​]∗A+τL​kf​​∗A(t−τf​)cosθ d ∅ d t = α 2 [ g ( N − N 0 ) 1 + ϵ A 2 − τ P − 1 ] − k f τ L ∗ A ( t − τ f ) A sin ⁡ θ      \frac{d\emptyset}{dt}=\frac{\alpha}{2}\left[\frac{g\left(N-N_0\right)}{1+\epsilon A^2}-\tau_P^{-1}\right]-\frac{k_f}{\tau_L}\ast\frac{A\left(t-\tau_f\right)}{A}\sin{\theta}\ \ \ \ dtd∅​=2α​[1+ϵA2g(N−N0​)​−τP−1​]−τL​kf​​∗AA(t−τf​)​sinθ     d N d t = I q V − N τ N − g ( N − N 0 ) A 2 \frac{dN}{dt}=\frac{I}{qV}-\frac{N}{\tau_N}-{g\left(N-N_0\right)A}^2 dtdN​=qVI​−τN​N​−g(N−N0​)A2其中公式的第一项描述内腔场变化，第二项是反馈光，表示反馈光对输出的振幅和相位的影响，根据半导体激光器不同的内腔，方程有不同的形式。方程是时滞的，因此无法普通求解，这里用到前面写的MATLAB混沌系统仿真其一：Lorenz系统和Rossler系统的龙格库塔四阶算法求解。仅是传递的参数变得更复杂，但是在Python上实现的思路也是一样的。因此也非常简单。

程序设计主要按照以下思路，龙格库塔四阶算法求解实际上就是不断地迭代，与Python系统仿真其一：logistics递归方程求解logistics的方法相同。 求解初值和调用函数的核心代码如下

由于不像MATLAB需要传递一个函数句柄，因此直接设计一个描述L&K方程的类内函数在RK4函数内不断调用迭代即可，参考MATLAB混沌系统仿真其一：Lorenz系统和Rossler系统，这里不做过多介绍。在不同反馈强度下，半导体激光器呈现不同的输出： 有了输出，不断改变反馈强度的值，暴力迭代即可得到分岔图。