平行平面腔自再现模的matlab数值计算

《平行平面腔自再现模的matlab数值计算》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、平行平面腔自再现模的数值解法石鹏（2010302749）目前，平行平面腔在中等以上激光器中仍普遍应用，但是平行平面腔的自再现模积分方程至今尚未得到精确的解析解。本文就条形腔和矩形对积分方程做数值计算。计算的思想基于Fox-Li数值迭代法，计算用matlab软件实现，结果可以很直观地用图形表示出来。1.矩形腔的分离变量法先来看谐振腔的自再现模。所谓自再现，就是指不管初始分布的具体特征如何，只要经过足够多次的往返渡越后，所生成的场分布都将明显带有衍射的痕迹，而且，具有这种特点的场分布将不再受到衍射作用的影响，形成一种稳定的场分布。设为经过次

2、渡越后在某一镜面上的光场分布，根据菲涅耳—基尔霍夫衍射积分公式，可以得到光波再渡越一次腔长距离后的光场分布与之间的关系为：（1）按照自再现模的概念，和应满足关系：。光场分布的模长为振幅分布，相角为相位分布。对于矩形腔，的计算公式为：按幂级数展开近似为：代入（1），可以得到积分方程的核为：可以看出，方形镜平行平面腔的自再现模的积分方程对、两个坐标是对称的。令，（2）这样就可以对积分方程进行分离变量，和的方程具有相同的形式：（3）这就是条形腔的积分方程，求解矩形腔的自再现模的关键就是解这个方程。2.积分方程的数值解法积分方程（1），（3）在

3、数学上称为第二类弗里德霍姆方程。这儿使用Fox-Li数值迭代的方法来求解。方程（3）可以表示为以下形式：（4）（5）其中为积分的核函数。方程（4）的含义是场分布经过一次传输以后得到新的场分布。根据自再现模的性质，不管初始分布具有什么样的形式，用方程（4）经过足够多次迭代得到的光场分布都具有相同的形式。现在任意给一个初始分布，来解方程（4）。先将方程（4）离散化：（6）再次迭代就要用到，把代入上式，令分别等于1,2,3，计算的前几项，；；，可以发现（6）式具有矩阵的形式，把它写出来：（7）把上式简化表示为，可以发现矩阵表示有特殊的含义：在

4、（4）中，是一个关于和的函数，而在（7）中是一个数值矩阵。传输一次就是给左乘一个，于是，初始分布经过q次传输以后就是给初始分布左乘一个：（8）这就是积分方程（4）的数值方法的迭代公式。可以发现，用积分方程迭代的重点就是计算核矩阵。我把它做成了一个函数文件kernel：functionK=kernel(Lambda,L,x)%KERNEL计算积分方程的核Kk=2*pi/Lambda;[X,Y]=meshgrid(x);K=sqrt(1i/(Lambda*L)*exp(-1i*k*L))*exp(-1i*k*(Y-X).^2/(2*L));

5、endKernel函数的输入参数为波长、腔长和一个把坐标轴离散化后的向量x。其中用到了一个生成网格数据的函数meshgrid，它可以得到一个横向向右递增的矩阵X和一个纵向向下递增的矩阵Y，用来计算矩阵。这种算法的结果得到了一个数值矩阵，并且把函数积分的N次迭代转化为对数值矩阵求N次方。matlab计算矩阵的速度比计算函数的速度要快，这就增加了程序的效率，优化了算法。这儿得到的中并没有步长2a/n，因为n是在函数外部定义的，所以在函数外部调用比较方便。3.条形腔的自再现模方程（3）就是条形腔的自再现模的积分方程，经过前面的讨论，它的解法就

6、是利用公式（8）来进行迭代。计算条形腔的振幅和相位分布的函数为bar1，下面给出bar1函数的一段代码：M=500;T=300;x=linspace(-a,a,M);x1=ones(M,1);K=kernel(Lambda,L,x);y=(2*a/(M-1)*K)^T*x1;y0=max(abs(y));y1=abs(y)/y0;y2a=angle(y);y2=xiugaiAngle(y2a,M);其中，M为离散程度，T为传播次数，即迭代次数。代码中有一个全一列向量x1，这就是我给的初始分布，全一表示它是一个平面波。第三行生成了一个需要

7、的K矩阵，真正进行迭代计算的是第四行。步长2*a/(M-1)在这儿调用，迭代了T次。这儿输出的y就是最终的光场分布。最后得到的y1为振幅分布，y2为相位分布。在计算相位分布时没有直接输出y2a，而是用到了一个函数xiugaiAngle，可以称它为修改相位函数。因为我们发现在计算相位时，angle函数的值域为，这就意味着当相位超过时angle函数会自动跳到，这样的话画出来的相位分布就会有一个个“坑”，如右图：要做的修改就是找出这些“坑”并且给它加上2。xiugaiAngle函数如下：functiony=xiugaiAngle(x,M)y=

8、x;x1=[x(1);x];x2=[x;x(M)];xa=x2-x1;y1=find(xa<-2*pi+0.2);y2=find(xa>2*pi-0.2);n=length(y1);forii=1:ny=[