怎么在matlab中解圆形波导的本征值的理论值

Question

怎么在matlab中解半圆形波导的本征值的理论值。目前有一个这样的波导

以下是现有的matlab代码，框架来源于@bit_ljf 的文章

clc
clear
%矩阵的读入
p=importdata('p.mat');
trig=importdata('trig.mat');
edges=find_edges(trig);
%矩阵创建
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);
%获取顶点编号：
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

%循环组装刚度矩阵和循环矩阵
%矩阵组装第二步，组装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
%处理边界条件并进行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(edges)
    A_tm(edges(1,i), :) = zeros(1, pn);% 将对应行全部设为0
    A_tm(:, edges(1,i)) = zeros(pn, 1);% 将对应列全部设为0
    B_tm(edges(1,i), :) = zeros(1, pn);% 将对应行全部设为0
    B_tm(:, edges(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;
n=299;
kc_2_tm=zeros(n,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
kc_tm=sqrt(kc_2_tm);
%根据理论公式对矩形波导进行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 );
%输出图像
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波理论值' );


function edges = find_edges(trig)
    % 初始化空的集合来存储所有的边
    edgeSet = containers.Map('KeyType', 'char', 'ValueType', 'double');

    % 遍历每个三角形
    for i = 1:size(trig, 2)
        % 获取三角形的三个顶点
        v1 = trig(1, i);
        v2 = trig(2, i);
        v3 = trig(3, i);

        % 提取三角形的三条边
        edges = [
            sort([v1, v2]);
            sort([v2, v3]);
            sort([v3, v1])
        ];

        % 将每条边加入集合
        for j = 1:3
            edgeKey = mat2str(edges(j, :));
            if ~isKey(edgeSet, edgeKey)
                edgeSet(edgeKey) = 1;
            end
        end
    end

    % 将集合中的边转化为矩阵
    edgeKeys = keys(edgeSet);
    numEdges = length(edgeKeys);
    edges = zeros(2, numEdges);

    for i = 1:numEdges
        edge = str2num(edgeKeys{i}); 
        edges(:, i) = edge(:);
    end
end

Answer 1

引自免费微信小程序：皆我百晓生

要解决半圆形波导的本征值问题，我们首先需要理解波导的数学模型和边界条件。对于圆形或半圆形波导，问题可以归结为求解亥姆霍兹方程在特定边界条件下的本征值。对于半圆形波导，我们可以将其视为圆柱坐标系下的问题，并利用分离变量法来简化方程。

下面是一个简化的步骤和示例代码，展示如何在MATLAB中求解半圆形波导的本征值。我们将采用有限元方法(FEM)来离散亥姆霍兹方程，并通过eigs函数来求解本征值问题。

步骤概述：

1. 定义几何：生成半圆形波导的网格。
2. 建立方程：设置亥姆霍兹方程。
3. 边界条件：应用适当的边界条件。
4. 离散化：使用有限元方法离散化方程。
5. 求解本征值：使用eigs函数求解本征值问题。

示例代码：

% 定义波导参数
radius = 1; % 半圆半径
length = 2; % 波导长度
lambda_min = 0.05; % 最小波长，用于确定网格密度

% 创建半圆几何
g = multisphere(radius);
g = rotate(g,[1 0 0],pi/2);
g = translate(g,[0 -radius 0]);
g = geomesh(g,'Hmax',lambda_min/2);
[p,e,t] = initmesh(g,'Hmax',lambda_min/2);

% 创建刚度矩阵和质量矩阵
A = assempde(b,p,e,t,1,0,1);
M = assema(p,e,t,1,0);

% 应用边界条件
c = applyBoundaryCondition('dirichlet','Edge',[1:4],'u',0);
[A,M] = assempde(c,[],[],1,0,1,p,e,t);

% 求解本征值问题
[V,D] = eigs(A,M,20,'sm');
D = diag(D);

% 输出本征值
disp('本征值:');
disp(D);

% 绘制前几个模式
figure;
for i = 1:5
    pdeplot(p,e,t,'XYData',V(:,i),'Contour','on');
    title(['Mode ',num2str(i)]);
    drawnow;
end

注意事项：

• 这段代码假定你有MATLAB的PDE Toolbox，它提供了构建和求解偏微分方程的工具。
• multisphere, initmesh, assempde, assema, 和 applyBoundaryCondition 是PDE Toolbox中的函数，用于定义几何、网格生成、组装刚度和质量矩阵以及应用边界条件。
• 使用eigs函数求解本征值时，可以选择求解最小的实数本征值，这对于波导模式分析非常有用。

这段代码提供了一个基础框架，你可能需要根据具体的波导尺寸和材料属性进行调整。此外，为了得到更准确的结果，你可能需要细化网格并检查收敛性。