规则与不规则波导本征值问题有限元求解（利用Matlab pdetool）

规则波导的有限元分析

求解原理

上一部分我们使用一维当中的平行电容求解，现在我们开始解决规则波导的求解

1. 实验原理与一维有限元分析原理相同，但是本征值问题并不是形如的求解，而是，除此之外，在区域离散化、基函数选取中也不相同。

2. 区域离散化
   将二维平面进行离散化，将其划分为若干个三角形，将三角形的顶点坐标、个数、边进行存储；
3. 选取插值基函数
   插值基函数选取为，通过插值条件判断：
   

4. 可得出具体插值基函数为：

   

   

   

5. 建立并组装刚度矩阵
   由下列方程建立矩阵
   
   
   再通过循环进行组装矩阵
   
   
6. 处理边界条件并求解
         由于TM波边界条件为0，需要在刚度矩阵设立边界条件，即在边界处设为0，TE波的边界条件已经被包含在泛函当中，无需再次设立。
   得到设立边界条件后的刚度矩阵与结果矩阵后可通过求解出，即该方程的特征值；
7. 求解精确解
          因为TE波与TM波的矩形空波导计算公式为可画出TE波与TM波前十个本征模的传播常数不进行对比

8. 可以通过数值解与精确解的一一比对，来确保数值解的精确性。

求解代码

在使用matlab求解时，使用pdetool进行网格划分

输入之后会弹出以下页面

我们可以按照波导的形状进行绘图

在绘图之后点击三角进行网格划分

最后点击mesh，底部有Export mesh

导出三角形顶点坐标p、每个单元节点的节点编号e、与三角形单元编号的trig三个矩阵。

点击OK，可以关闭pdetool页面，将文件保存为.m格式

1. 通过计算p与trig矩阵大小来确定刚度矩阵与结果矩阵的大小
         其中在矩阵创建时可以使用稀疏矩阵，但是稀疏矩阵带来的缺点是搜索时效率下降，而稠密矩阵中存储太多0使空间利用率较差​​​​​​​

   Np=size(p);%p is mesh of points
   Nt=size(trig);%trig is number of triangles
   %e is the edge of triangle
   pn=Np(2);
   tn=Nt(2);
   %矩阵组装中的第一步，将Ae与Be求出后组装
   b=zeros(3,1);
   c=zeros(3,1);
   Ae=zeros(tn,3,3);
   Be=zeros(tn,3,3);

2.  获取顶点编号，得到矩阵
   ​​​​​​​​​​​​​​

   for k=1:tn
       i1=trig(1,k); i2=trig(2,k); i3=trig(3,k);
       x1=p(1,i1);   x2=p(1,i2);   x3=p(1,i3);
       y1=p(2,i1);   y2=p(2,i2);   y3=p(2,i3);
       %求得三角形面积
       S = x3 * (y1 - y2) + x1 * (y2 - y3) + x2 * (y3 - y1);
       %从数学角度来讲是要求梯度
       b(1)=y2-y3;b(2)=y3-y1;b(3)=y1-y2;
       c(1)=x3-x2;c(2)=x1-x3;c(3)=x2-x1;
       %求梯度的内积
       for i=1:3
           for j=1:3
               Ae(k,i,j)=(b(i)*b(j)+c(i)*c(j))/(S);
               Be(k,i,j)=(1+~(i-j))*S/12;
           end
       end
   end

3. 循环组装刚度矩阵与结果矩阵

   %矩阵组装第二步，组装A与B
   A = zeros( pn , pn );
   B = zeros( pn , pn );
   for k=1:tn
       for i=1:3
           for j=1:3
               %求刚度矩阵A
               A(trig(i,k),trig(j,k)) = A(trig(i,k),trig(j,k)) + Ae(k,i,j);
               B(trig(i,k),trig(j,k)) = B(trig(i,k),trig(j,k)) + Be(k,i,j);
           end
       end
   end

    
4.  处理边界条件并进行TE波与TM波求解

   %求解TE波，TE波不需要进行边界条件处理
   [Hz,kc_TE]=eigs(A,B,10,'smallestreal');
   kc_te=sqrt(kc_TE);
   kc_2_te =  sort( eig( B\A ) );
   kc_Te = sqrt( kc_2_te );
   %求解TM波，TM波需要将边界置为0
   A_tm=A;
   B_tm=B;
   for i=1:length(e)
       A_tm(e(1,i), :) = zeros(1, pn);% 将对应行全部设为0
       A_tm(:, e(1,i)) = zeros(pn, 1);% 将对应列全部设为0
       B_tm(e(1,i), :) = zeros(1, pn);% 将对应行全部设为0
       B_tm(:, e(1,i)) = zeros(pn, 1);% 将对应列全部设为0
   end
   [Ez , D] = eig( pinv(B_tm) * A_tm );%求特征值和特征向量
   [kc_2_tm0 , I] = sort( diag(D) );%kc_2_tm是kc方的升序，其中可能含有0；I为索引
   Ez_tm = Ez(:,I); %对特征向量排序，使之对应,即为Ez
   j = 1;
   for i = 1:pn
       if(kc_2_tm0(i) == 0)
           continue;
       else
           kc_2_tm(j,1) = kc_2_tm0(i);    %kc_2_tm中存的是去掉0的kc方的值
           j = j+1;
       end
   end

5.  根据理论公式对举行波导进行TE波与TM波计算

   %计算TE波的理论值
   k=1;
   for m = 0:5
       for n = 0:5
           if( m==0 && n==0 )
               continue;
           else
               kc_th_2_te( k ) = ( m*pi / 0.4572).^2 + ( n*pi / 0.2286).^2;
               k = k+1;
           end
       end     
   end
   kc_th_2_te = sort( kc_th_2_te );
   
   %计算TM波的理论值
   k=1;
   for m = 1:5
       for n = 1:5
           kc_th_2_tm( k )=( m*pi / 0.4572).^2 + ( n*pi / 0.2286).^2;
           k = k+1;
       end     
   end
   kc_th_2_tm = sort( kc_th_2_tm );
   

6. 最后输出图像

   x = 1:10;
   plot( x, kc_2_te(2:11) );
   title('前10个本征模Kc平方的对比结果');
   hold on;
   plot( x, kc_2_tm(1:10) );
   hold on;
   plot( x, kc_th_2_te(1:10) );
   hold on;
   plot( x, kc_th_2_tm(1:10) );
   legend('TE波程序计算值', 'TM波程序计算值', 'TE波理论值', 'TM波理论值' );

    

不规则波导的求解

求解过程

求解代码与规则波导相同，需要注意的是，我们只需要将第一步画边界时画出自己需要的边界即可，利用区域的加减等方法即可。