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

平行平面腔自再现模的数值解法

石鹏（2010302749）

目前，平行平面腔在中等以上激光器中仍普遍应用，但是平行平面腔的自再现模积分方程至今尚未得到精确的解析解。本文就条形腔和矩形对积分方程做数值计算。计算的思想基于Fox-Li数值迭代法，计算用matlab软件实现，结果可以很直观地用图形表示出来。

1.矩形腔的分离变量法

先来看谐振腔的自再现模。所谓自再现，就是指不管初始分布u1的具体特征如何，只要经过足够多次的往返渡越后，所生成的场分布都将明显带有衍射的痕迹，而且，具有这种特点的场分布将不再受到衍射作用的影响，形成一种稳定的场分布。

设uq?x?,y??为经过q次渡越后在某一镜面上的光场分布，根据菲涅耳—基尔霍夫衍射积分公式，可以得到光波再渡越一次腔长距离后的光场分布uq?1与uq之间的关系为：

ikuq?1?x,y??4???u?x?,y??qe?ik???1?cos??ds? （1）

按照自再现模的概念，uq?1和uq应满足关系：uq?1??uq。光场分布uq的模长为振幅分布，相角为相位分布。

对于矩形腔，?的计算公式为：???x?x????y?y??222?L2 x?x???按幂级数展开近似为：??L?2Ly?y????2L2

??x?x??2?y?y??2??ik??2L2L???ikL??代入（1），可以得到积分方程的核为：K?x,y,x?,y???ie?Le

可以看出，方形镜平行平面腔的自再现模的积分方程对x、y两个坐标是对称的。令

umn?x,y??um?x?un?y?，?mn??m?n （2）

这样就可以对积分方程进行分离变量，x和y的方程具有相同的形式：

?x?x??22Li?ikLa?ik?mum?x??e?e?a?Lum?x??dx? （3）

这就是条形腔的积分方程，求解矩形腔的自再现模的关键就是解这个方程。

2.积分方程的数值解法

积分方程（1），（3）在数学上称为第二类弗里德霍姆方程。这儿使用Fox-Li数值迭

代的方法来求解。方程（3）可以表示为以下形式：

uq?1?x???K?x?x??uq?x??dx? （4）

?aaK?x?x???i?ikLee?L2x?x????ik2L （5）

其中K?x?x??为积分的核函数。方程（4）的含义是场分布uq?x?经过一次传输以后得到新的场分布uq?1?x?。根据自再现模的性质，不管初始分布u1?x?具有什么样的形式，用方程（4）经过足够多次迭代得到的光场分布都具有相同的形式。

现在任意给一个初始分布，来解方程（4）。先将方程（4）离散化：

2anuq?1?x??Kx?xi?uqxi? （6） ?ni?1再次迭代就要用到uq?1xi?，把xi?代入上式，令i分别等于1,2,3，计算uq?1xi?的前

n2an2a几项，uq?1?x1??Kx1?xi?uqxi?；uq?1?x2??Kx2?xi?uqxi?；??ni?1ni?1????????????????2anuq?1?x3??Kx3?xi?uqxi?，可以发现（6）式具有矩阵的形式，把它写出来： ?ni?1?????uq?1?x1???K?x1?x1?K?x1?x2????ux??2a?q?12???K?x2?x1?K?x2?x2??…?n?……???K?x?x?K?x?x??u?x???n1n2??q?1n?…K?x1?xn???uq?x1?????…K?x2?xn???uq?x2??（7）

??…?……???u?x???…K?xn?xn????qn?把上式简化表示为Uq?1?KUq，可以发现矩阵表示有特殊的含义：在（4）中，K?x?x??是一个关于x和x?的函数，而在（7）中K是一个数值矩阵。传输一次就是给Uq左乘一个K，于是，初始分布经过q次传输以后就是给初始分布左乘一个K：

qUq?KqU1 （8）

这就是积分方程（4）的数值方法的迭代公式。可以发现，用积分方程迭代的重点就是

计算核矩阵K。我把它做成了一个函数文件kernel ：

function K = kernel( Lambda , L , x ) %KERNEL 计算积分方程的核K k = 2*pi/Lambda ;

[ X , Y ] = meshgrid(x) ;

K = sqrt( 1i/(Lambda*L) * exp(-1i*k*L) ) * exp(-1i*k * (Y-X).^2 / (2*L) ) ; end

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

3.条形腔的自再现模

方程（3）就是条形腔的自再现模的积分方程，经过前面的讨论，它的解法就是利用公式（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，这就是我给的初始分布，全一表示它是一个平面波。第三行生成了一个需要的K矩阵，真正进行迭代计算的是第四行。步长2*a/(M-1)在这儿调用，迭代了T次。这儿输出的y就是最终的光场分布。最后得到的y1为振幅分布，y2为相位分布。

在计算相位分布时没有直接输出y2a，而是用到了一个函数xiugaiAngle，可以称它为修改相位函数。因为我们发现在计算相位时，angle函数的值域为?-?，??，这就意味着当相位超过?时angle函数会自动跳到-?，这样的话画出来的相位分布就会有一个个“坑”，如右图：

要做的修改就是找出这些“坑”并且给它加上2? 。xiugaiAngle函数如下：

function y = xiugaiAngle( x , M ) y = 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 ) ; for ii = 1 : n

y = [ y(1:y1(ii)-1) ; y(y1(ii):y2(ii)-1)+2*pi ; y(y2(ii):M) ] ; end end

用find函数找到“坑”的位置，用for循环给它加上2? 。

4.矩形腔的自再现模

从方程（1）由（2）到方程（3），说明了矩形腔的分离变量法。现在解出了方程（3），那么就可以由方程（2）到方程（1），计算矩形腔的自再现模。计算用的函数是rectangleq2 。

下面一段代码是计算矩形腔的振幅分布的：

Z = xa*conj(ya)' ;

Z1 = abs(Z)/max( max(abs(Z)) ) ;

其中xa和ya分别代表了x方向的光场分布和y方向上的光场分布。它们都是列向量，对ya转置（matlab中的算符 ‘ 表示取一个矩阵的共轭转置，还要再求一次共轭才能得到转置），代码中的乘法是矩阵乘法，也就是说乘得的结果是一个矩阵，也就是由公式（2）得到的光场分布。

由于xiugaiAngle函数是一维的，所以不能用上面得到的Z来计算相角分布。相角分布的计算由下面一段代码实现：

xa2 = angle(xa) ; ya2 = angle(ya) ;

xb2 = xiugaiAngle( xa2 , M ) ; yb2 = xiugaiAngle( ya2 , M ) ; [ X2 , Y2 ] = meshgrid( xb2 , yb2 ) ; Z2 = X2 + Y2 ;

复数求积后相角是相加的，上面第四行代码就是实现相角相加。

参考文献：

[1]李东 激光原理与技术课件

[2]俞宽新 激光原理与机关技术 北京工业大学出版社 149-164 [3]孙蓬 matlab基础教程

[4]吴鹏 MATLAB高效编程技巧与应用：25个案例分析 北京航空航天大学出版社 98-99

说明：一共有7个m文件：

bar000.m 脚本文件 是这道题的第一个成功的m文件，画的是条形腔的振幅分布

d = 1 ; k = 2*pi/d ; a = 25*d ; L = 100*d ; M = 300 ; T = 300 ;

x = linspace(-a,a,M) ; x1 = ones( M , 1 ) ; [ X , Y ] = meshgrid(x) ;

K1 = sqrt( 1i/(d*L) * exp(-1i*k*L) ) ; K = K1 * exp(-1i*k * (Y-X).^2 / (2*L) ) ;

y1 = (2*a/(M-1) * K)^T *x1 ; y = abs(y1)/max(abs(y1)) ; plot(x,y)

bar1.m 函数文件 用来计算并绘制条形腔的振幅和相位分布

function bar1( Lambda , a , L )

or1 绘制条形腔的振幅分布和相位分布 % 输入参数为波长，镜面宽度，腔长

M = 500 ; T = 300 ; %M为离散程度，T为传播次数 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 ) ;

na_L = int2str(L) ; na_a = int2str(a) ;

na_M = int2str(M) ; na_T = int2str(T) ; namez = [ '条形腔振幅分布' ] ; namex = [ '条形腔相位分布' ] ;

name_xq = [ 'a：' na_a '倍波长' 'L：' na_L '倍波长' ] ; name_xh = [ '离散程度：' na_M '迭代次数：' na_T ] ;

figure(1) ; plot(x,y1)

title( [ namez char(10) name_xq char(10) name_xh] ) ; saveas(gcf,[ namez na_a '_' na_L '.jpg' ]) figure(2) ; plot(x,y2)

title( [ namex char(10) name_xq char(10) name_xh] ) ; saveas(gcf,[ namex na_a '_' na_L '.jpg' ]) end

barq.m 脚本文件 调用bar1.m来实现条形腔光场的计算

disp('绘制条形腔的振幅分布和相位分布') ; Lambda = input('请输入波长 λ = ') ; a = input('请输入镜面宽度 2a = ') /2 ; L = input('请输入腔长 L = ') ; bar1( Lambda , a , L ) ;

kernel.m 函数文件 计算积分核 K

function K = kernel( Lambda , L , x ) %KERNEL 计算积分方程的核 K % 可以输入参数为波长和腔长

k = 2*pi/Lambda ;

[ X , Y ] = meshgrid(x) ;

K = sqrt( 1i/(Lambda*L) * exp(-1i*k*L) ) * exp(-1i*k * (Y-X).^2 / (2*L) ) ; End

rectangleq.m 脚本文件 调用rectangle2.m计算矩形腔的光场分布

disp('绘制矩形腔的振幅分布和相位分布') ; Lambda = input('请输入波长 λ = ') ; a = input('请输入镜面长度 2a = ') /2 ; b = input('请输入镜面宽度 2b = ') /2 ; L = input('请输入腔长 L = ') ; rectangleq2( Lambda , a , b , L )

rectangleq2.m 函数文件 计算矩形腔的光场分布

function rectangleq2( Lambda , a1 , b1 , L1 )

%RECTANGLE2 绘制矩形腔的振幅分布和相位分布 % 修补相位分布的缺陷

a = a1/Lambda ; b = b1/Lambda ; L = L1/Lambda ; M = 500 ; T = 300 ;

x = linspace(-a,a,M) ; y = linspace(-b,b,M) ; x0 = ones( M , 1 ) ;

[ X , Y ] = meshgrid( x , y ) ;

Kx = kernel( Lambda , L , x ) ; xa = (2*a/(M-1) * Kx)^T *x0 ;

Ky = kernel( Lambda , L , y ) ; ya = (2*b/(M-1) * Ky)^T *x0 ;

Z = xa*conj(ya)' ;

Z1 = abs(Z)/max( max(abs(Z)) ) ;

%计算相位分布

xa2 = angle(xa) ; ya2 = angle(ya) ;

xb2 = xiugaiAngle( xa2 , M ) ; yb2 = xiugaiAngle( ya2 , M ) ; [ X2 , Y2 ] = meshgrid( xb2 , yb2 ) ; Z2 = X2 + Y2 ;

na_L = int2str(L) ;

na_a = int2str(a) ; na_b = int2str(b) ; na_M = int2str(M) ; na_T = int2str(T) ; namez = '矩形腔振幅分布' ; namex = '矩形腔相位分布' ;

name_xq = [ 'a：' na_a '倍波长' ' b：' na_b '倍波长' 'L：' na_L '倍波长' ] name_xh = [ '离散程度：' na_M '迭代次数：' na_T ] ;

figure(1) ; mesh( X , Y , Z1 )

title( [ namez char(10) name_xq char(10) name_xh] ) ; saveas(gcf,[ namez na_a '×' na_b '.jpg' ]) figure(2) ; mesh( X , Y , Z2 )

title( [ namex char(10) name_xq char(10) name_xh] ) ; saveas(gcf,[ namex na_a '×' na_b '.jpg' ]) end

xiugaiAngle.m 函数文件 修改相位函数

function y = xiugaiAngle( x , M )

%UNTITLED2 Summary of this function goes here % Detailed explanation goes here

y = 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 ) ; for ii = 1 : n

y = [ y(1:y1(ii)-1) ; y(y1(ii):y2(ii)-1)+2*pi ; y(y2(ii):M) ] ;

; end end

有11张图片，图片名就是图片的说明

百度搜索“77cn”或“免费范文网”即可找到本站免费阅读全部范文。收藏本站方便下次阅读，免费范文网，提供经典小说综合文库平行平面腔自再现模的matlab数值计算在线全文阅读。