使用差分法计算矩形波导前20截止频率（matlab实现）
一、任务要求
(a) Writea finite difference (FD) code to compute the cutoff frequencies of the firsttwenty TM modes of a rectangular waveguide with dimensions 1 cm by 2 cm.Compare your results to the analytical cutoff frequencies. How rapidly do theFD results approach the exact results as the number of grid.points increases?
(b) Same as (a), but forTE modes.
(c) Which case (TE or TM) has moreaccurate results
comparing with the analytical results? Why?
二、波导问题分析
求解波导问题相当于求解亥姆霍兹方程：
（1）
上式中，对于TM模式， ，对于TE模式， ，其中 是波数。
二者边界条件为：
1）对于TM波，
2）对于TE波，
在上述方程的本征值问题中， 和 都是待定，截止波长 ，对截止波长的每一个值都有一个本征函数对应，代表一个传播模式。
使用有限差分法，用正方形来剖分波导横截面得到：
(2)
式中: 为网格尺寸。
可以得到N个方程
求解A矩阵的特征值可得到截止波数即： 。
三、matlab代码实现
1) FD for TM
对于一个20mm*10mm的波导，
网格划分图：
边缘节点为固定矩阵，TM边界为0 ，网格大小为a/nx
通过计算得到前20个截止频率为：
由上表可以看出：点数越多，计算越精确。当网格长度是波长的0.1倍时，可以计算出较为准确值。
2） FD for TE
边缘节点为虚拟节点，网格大小为a/nx-1
TE前20截止频率
由上表可以看出：点数越多，计算越精确。但TE计算没TM较为准确，是由于在边界条件上多了一阶近似。
Matlab：代码
FD：TM
clc
clear all
a=20*1e-3; %长边20mmb=10*1e-3; %短边10mm
c=2.99792458*1e8; %光速
nx=10; %x方向剖分网格
ny=nx*0.5; %y方向剖分网格保证为正方形
h=a/nx; %网格大小
N=(nx-1)*(ny-1); %未知数个数
A=zeros(N,N);
k=0 ; %用于填充m
y1=zeros(2,20); %前20个截止频率比较
for j=1:ny-1
for i=1:nx-1
k=k+1
NL(j,i) = k;
end %NL用于i、j和m的对应关系
end
%%%%%填充A矩阵%%%%%%%
for m=1:N
for n=1:N
[jm,im]=find(NL==m)
[jn,in]=find(NL==n);
if m==n
A(m,n)=4;
elseif im==in&jm==jn+1
A(m,n)=-1;
elseif im==in&jm==jn-1
A(m,n)=-1;
elseif im==in+1&jm==jn
A(m,n)=-1;
elseif im==in-1&jm==jn
A(m,n)=-1;
else
A(m,n)=0;
end
end
end
%%%%%%%填充A矩阵%%%%%%%
la=eig(A); %求特征值la
la=round(la,4); %保留小数点后4位
la=unique(la); %去除相同元素
Kc_fd=sqrt(la)./h;
lamda_fd=2*pi./Kc_fd;
f_fd=c./lamda_fd;
B1=f_fd;
for i=1:20
[p,q]=min(B1); %取截止频率最小前20
B1(q)=inf;
y1(1,i)=p; %fd
end
%%%%%%%理论计算截止频率 TM11-88 %%%%%%%%%%%
for m=1:8
for n=1:8
Kc_cal(m,n)=sqrt( (m*pi/a)^2 +(n*pi/b)^2 )
lambda_cal(m,n)=2*pi/Kc_cal(m,n)
f_cal(m,n)=c/lambda_cal(m,n);
end
end
B2=sort(f_cal(:),1,'ascend')
B2=round(B2,4) %保留小数点后4位
B2=unique(B2) %去除相同元素
for i=1:20
[p,q]=min(B2);
B2(q)=inf;
y1(2,i)=p; %cal
end
y1=y1*1e-9; %换算为GHz
y1=round(y1,4); %保留小数点后4位
disp(y1) % y1第一行为FD计算值 第二行为理论计算值
FD:TE
clc
clear all
a=20*1e-3; %长边20mm
b=10*1e-3; %短边10mm
c=2.99792458*1e8; %光速
nx=10; %x方向剖分网格
ny=nx*0.5; %y方向剖分网格
h=a/nx; %网格大小
N=(nx-1)*(ny-1); %未知数个数
A=zeros(N,N);
k=0 ; %用于填充m
y1=zeros(2,20); %前20个截止频率比较
for j=1:ny-1
for i=1:nx-1
k=k+1
NL(j,i) = k;
end %NL用于i、j和m的对应关系
end
%%%%%填充A矩阵%%%%%%%
for m=1:N
for n=1:N
[jm,im]=find(NL==m)
[jn,in]=find(NL==n);
if m==n
A(m,n)=4;
elseif im==in&jm==jn+1
A(m,n)=-1;
elseif im==in&jm==jn-1
A(m,n)=-1;
elseif im==in+1&jm==jn
A(m,n)=-1;
elseif im==in-1&jm==jn
A(m,n)=-1;
else
A(m,n)=0;
end
%%%%%%%%%%% TE边界条件 %%%%%%%%%
if im==in&jm==jn-1&jm==1 %下边界
A(m,n)=-2;
elseif im==in&jm==jn+1&jm==(ny-1) %上边界
A(m,n)=-2;
elseif im==in-1&jm==jn& im==1 %左边界
A(m,n)=-2;
elseif im==in+1&jm==jn&im==(nx-1) %右边界
A(m,n)=-2;
end
end
end
%%%%%%%填充A矩阵%%%%%%%
la=eig(A); %求特征值la
la=round(la,4); %保留小数点后4位
la=unique(la); %去除相同元素
Kc_fd=sqrt(la)./h;
lamda_fd=2*pi./Kc_fd;
f_fd=c./lamda_fd;
B1=f_fd;
B1(find(B1==0))=[]; %去除0
for i=1:20
[p,q]=min(B1); %取截止频率最小前20
B1(q)=inf;
y1(1,i)=p; %fd
end
%%%%%%%理论计算截止频率 TE10- TE87 %%%%%%%%%%%
for m=1:8
for n=1:8
Kc_cal(m,n)=sqrt( (m*pi/a)^2 +((n-1)*pi/b)^2 )
lambda_cal(m,n)=2*pi/Kc_cal(m,n)
f_cal(m,n)=c/lambda_cal(m,n);
end
end
B2=sort(f_cal(:),1,'ascend')
B2=round(B2,4) %保留小数点后4位
B2=unique(B2) %去除相同元素
for i=1:20
[p,q]=min(B2);
B2(q)=inf;
y1(2,i)=p; %cal
end
y1=y1*1e-9; %换算为GHz
y1=round(y1,4); %保留小数点后4位
disp(y1) % y1第一行为FD计算值 第二行为理论计算值
新人贴，好紧张呐