程序结果：
1、
时域有限差分法
问题：在 TM 模式下，编写 matlab 程序，使用 FDTD 方法分析点脉冲源激励情
况下，正方形 MUR 边界的性能。
分析：对 TM 波，FDTD 公式为：
Hx
n 1 2 1 n Ezn (i, j 1) Ezn (i, j ) 1 1 2 (i , j ) cp (m) H x (i , j ) cq (m) 2 2 y
N N ' k L( J sz ) i L[ i (r r i )] i H 02(k | r r i |)dl 4 li i 1 i 1
若以狄拉克函数作为权函数切令它定义在各分段的中点 ( x j , y j ) 上，即令
W j ( x x j , y y j ) …………………………………………………… （1-2）
Matlab 程序如下：
yps=8.85*10^-12;%真空中介电常数 u=4*pi*10^-7; c=3*10^8;%真空中光速 L=0.1; %TM 波约定： %Hx(i,j)=Hx(i*deltaX,(j+1/2)*deltaY); %Hy(i,j)=Hy((i+1/2)*deltaX,j*deltaY); %Ez(i,j)=Ez(i*deltaX,j*deltaY); %区域设定 N=10 %将正方形区域均分成为 N*N 个小正方区域 %初始电磁场值 %旧时刻场量值 Ez1=zeros(N+1,N+1);
《电磁场问题数值算法》学习经典习题
吕波
矩量法
问题：求解金属圆柱在水平入射均匀平面波照射下，面电流分布及由此产生的远 区散射场分布。 分析：可将面电流密度用脉冲函数展开，即：
N J sz i i (r r i ) ……………………………………………………（1-1） i 1
它表示把导体横截面的周界分成 N 个小段， 且设每一个小段中的面电流密度 i 作 均匀分布。这样算子方程为：
Ezn 1 (i , j ) Ezn (i 1, j 1) Ezn 1 (i , j ) Ezn (i 1, j 1) Ezn 1 (i , j ) Ezn (i 1, j 1) Ezn 1 (i , j ) Ezn (i 1, j 1) Ezn 1 (i 1, j 1) Ezn (i , j ) （左下角） ct 2 Ezn 1 (i 1, j 1) Ezn (i , j ) （右上角） ct 2 Ezn 1 (i 1, j 1) Ezn (i , j ) （左上角） ct 2 Ezn 1 (i 1, j 1) Ezn (i , j ) （右下角） ct 2 ct 2 ct 2 ct 2 ct 2
[g]中个元素可按下面的分析定之。设自由空间中有一均匀平面波沿 r ' 方向传播，
r ' 和 x 轴之间的夹角为θi，则在空间任意点（x,y）处的电场强度可表示为： Ez( i ) e jk ( x cosi y sin i )
导体表面上匹配点处的电场强度应为
Ez( i ) ( x 以转化为代数方程组
[l ][ ] [ g ] ……………………………………………………………（1-3）
式中 （1-4） g j ( x x j , y y j ), Ezi Ezi ( x j , y j ) ………………………………
k (2) lij l ji ( x x j , y y j ), L i (r r i ) H0 (k ( xi x j ) 2 ( yi y j ) 2dli （1-5） 4 li
%计算 l 矩阵各个元素 D=1.781; for wx=1:N % this is for fufill the matrix for wy=1:N x(wx)=r*cos(pi/N+((wx-1)*2*pi)/N); x(wy)=r*cos(pi/N+((wy-1)*2*pi)/N); y(wx)=r*sin(pi/N+((wx-1)*2*pi)/N); y(wy)=r*sin(pi/N+((wy-1)*2*pi)/N); if wx==wy l(wx,wy)=(k*yita*detaL/4)*(1-j*2*(log(D*k*detaL/4)-1)); else l(wx,wy)=(k*yita/4)*besselh(0,2,k*sqrt((x(wx)-x(wy))^2+(y(wx)-y(wy))^2))*detaL; l(wy,wx)=(k*yita/4)*besselh(0,2,k*sqrt((x(wx)-x(wy))^2+(y(wx)-y(wy))^2))*detaL; end end end;

Hx
n
其中，真空中 ca(m)=cp(m)=1; cb(m)
0
t
; cq (m)
0
t
二维 Mur 吸收边界条件的 FDTD 形式如下（一阶） ： ct n 1 Ez (i 1, j ) Ezn (i , j ) Ezn 1 (i , j ) Ezn (i 1, j ) （左边界） ct ct n 1 Ezn 1 (i , j ) Ezn (i 1, j ) Ez (i 1, j ) Ezn (i , j ) （右边界） ct ct n 1 Ezn 1 (i , j ) Ezn (i , j 1) Ez (i , j 1) Ezn (i , j ) （上边界） ct ct n 1 Ezn 1 (i , j ) Ezn (i , j 1) Ez (i , j 1) Ezn (i , j ) （下边界） ct 二维 Mur 吸收边界角点的处理如下：
Hx1=zeros(N+1,N); Hy1=zeros(N,N+1); %迭代后新时刻场量值 Ez2=zeros(N+1,N+1); Hx2=zeros(N+1,N); Hy2=zeros(N,N+1); %设置脉冲源 tao=2*10^-10;%高斯脉冲的宽度 t0=tao/2%高斯脉冲峰值出现的时刻 deltaL=L/N;%空间间隔要求 deltaT=tao/(4*N);%时间间隔要求 t=0;%可以修改 t 的值一观察不同时刻的电磁场 %迭代过程 for dc=1:1000%迭代的次数 ca=1;%迭代系数 cb=deltaT/yps; cp=1; cq=deltaT/u; %激励源加入 It=exp(-4*pi*(t-t0).^2/(tao^2)); %非边界场值迭代 clear i j; for i=1:(N+1) for j=1:N Hx2(i,j)=cp*Hx1(i,j)-cq*(Ez1(i,j+1)-Ez1(i,j))/deltaL; end end clear i j; for i=1:N for j=1:(N+1) Hy2(i,j)=cp*Hy1(i,j)+cq*(Ez1(i+1,j)-Ez1(i,j))/deltaL; end end clear i j; for i=2:N for j=2:N
jk ( x j cos i y j sin i )
（j=1,2,…,N）……………………（1-8）
Matlab 源程序如下：
%矩量法计算圆柱散射 %copyright by lvbo % clear; yita=120*pi%自由空间波阻抗 freq=5*10^8%频率为 500MHz k=(2*pi*freq)/(3*10^8)%该频率的传播常数 N=100%将圆柱导体横截面的周界 l 分成 N 个小段 r=0.1%圆柱导体的半径为 0.1 detaL=2*pi*r/N%每个小段的长度 2003,6,24
当 i j 时，采用近似式
lij l ji k (2) H 0 (k (xi x j )2 ( yi y j )2 )li ………………………（1-6） 4
当 i j 时，
lii k li 4 Dk li 2 1 j (ln 4 1) ……………………………………（1-7）
%给 g 矩阵赋值 %入射波以 0 度(水平)入射圆柱导体
for px=1:N g(px)=exp(j*k*(x(px)*cos(0)+y(px)*sin(0))); end
%求解电流展开系数 afa=g*inv(l);
%求解面电流密度 for py=1:N J(py)=afa(py); magJ(py)=abs(J(py)); phasJ(py)=angle(J(py)); end subplot(2,2,1); polar(pi/N:2*pi/N:2*pi,magJ); title('面电流密度幅度分布'); subplot(2,2,2); plot(phasJ); title('面电流密度相位分布'); %求解远区散射场 fai=pi/N:2*pi/N:2*pi; p=1000%所求远区场与柱体的距离 Es=0; for pj=1:N Es=Es+(k*yita*exp(-j*k*p+0.75*pi)/sqrt(j*8*k*p))*afa(pj)*detaL*exp(-j*k*(x(pj)*cos(fai)+y(pj) *sin(fai))); end subplot(2,2,3); polar(fai,abs(Es)/max(abs(Es))); title('电场幅度随离开柱体距离的变化趋势'); subplot(2,2,4); plot(1:N,angle(Es)); title('电场相位随角度变化趋势');