[求助]请帮忙看下我做的fdtd程序

小弟刚学fdtd不久，照葛德彪那本书天线那张做了个二维柱坐标系下的TE波线电流源，激励源采用双指数函数，吸收边界采用一阶Mur x)eF{%QB  
感觉激励源的加入还是有很大问题 2/A*\  
另外步长改为1的话出来的波形明显震荡 H{i|?
a)  
请各位高手指点一下 y<nPZ<h  
RS&BS;  
n,s7!z/  
-@]b7J?`k  
SvkCx>6/G  
clear all ^C~R)M:C  
close all <2<2[F5Q%  
clc ^jRX6  
&M3KJ I0L  
"Vl4=W)u  
%-------------------------------------------------- a @3s71  
%FDTD法模拟柱坐标二维TE波的传播(Mur吸收边界条件) aY.cx1"  
%-------------------------------------------------- Yd~J(  
pJ1\@G  
! N!pvK;  
%初始化 \I}EWI  
c=3e8;                                                       %光速 ':tdb$h  
delta=10; }:$ot18  
% delta=1e-3;                                                  %空间步长 hP.Km%C)0n  
dt=delta/(2*c);                                              %时间步长 ='+I dn#5  
u0=(4*pi)*1e-7;                                              %磁导率 \JR^uJ{Y  
e0=(1/(36*pi))*1e-9;                                         %介电常数 KTot40osj  
I0=15e3;                                                     %基电流 ;-"q;&1e  
n0=200;                                                      %空间网格数 >)diXe}j  
% n1=33;                                                       %源的位置 d)Z&_v<|  
abi=4*dt/(pi*e0*delta^2);                                    %激励源系数 >/ A'G  
CA=1; |w
#~v%w  
CB=dt/(e0*delta);                                            %Er与Ez中的系数 'xnI Nu  
CP=1; tEo-Mj5:  
CQ=dt/(u0*delta);                                            %Hphi中的系数 ue+{djz[4  
C1=(1/(2*delta)-1/(2*c*dt)-1/(8*delta*(n0)))/(1/(2*delta)+1/(2*c*dt)+1/(8*delta*(n0)));                              %一阶边界条件中的系数 :HrFbq  
C2=(-1/(2*delta)+1/(2*c*dt)-1/(8*delta*(n0)))/(1/(2*delta)+1/(2*c*dt)+1/(8*delta*(n0))); 4e|N^h*!  
C3=(1/(2*delta)+1/(2*c*dt)-1/(8*delta*(n0)))/(1/(2*delta)+1/(2*c*dt)+1/(8*delta*(n0))); ?>V>6cDQ  
% C1=-(4*n0+1)/(12*n0+1);                              %一阶边界条件中的系数 A*/8j\{n  
% C2=-(-4*n0+1)/(12*n0+1); nnT#S  
% C3=-(-12*n0+1)/(12*n0+1); 2,g4yXws5  
% C1=-(4*(n0+1)+1)/(12*(n0+1)+1);                              %一阶边界条件中的系数 @'M"c q  
% C2=-(-4*(n0+1)+1)/(12*(n0+1)+1); a-FI`Dv  
% C3=-(-4*(n0+1)+1)/(-12*(n0+1)+1); R)SY#*Y  
C4=(1.3e8*dt-delta)/(1.3e8*dt+delta); E>5p7=Or;"  
Er=zeros(n0+1,n0+1); |dqESl,2  
Ez=zeros(n0+1,n0+1); T2rBH]5  
Hphi=zeros(n0,n0); iV#A-9  
Hphi_last=zeros(n0,n0);                                 %大型矩阵维数的预先确定 kQd|qZ=:w  
i0+e3!QU  
I#;dS!W"'  
^
LB]  
L{r4hL [  
for k=1:2000  CDuA2e  
L$);50E  
N`?/kubD  
    %迭代求电场 eYlI};  
+zLw%WD[l  
~a_X 7
 
#6D>e~>n  
    %求Er K*,,j\Q.  
    for j=2:n0 4]p#9`j  
        for i=1:n0 f9kd&#O&  
            Er(i,j)=CA*Er(i,j)-CB*(Hphi(i,j)-Hphi(i,j-1)); +|X`cmnuU  
        end b$*2bSdv0<  
    end                                                %Er循环结束 q`L)^In"  
    %求Ez [dFcxzM-N  
    for j=1:n0 o^"OKHU,S0  
        for i=2:n0 r1vS~ 4Z  
            Ez(i,j)=CA*Ez(i,j)+CB*(delta*i*Hphi(i,j)-delta*(i-1)*Hphi(i-1,j))/(delta*(i-0.5)); rMjb,2*rC7  
        end ns\I Y<Yo  
    end                                                %Ez循环结束 {dRZ2U3  
    %激励源 $- %um  
    t=k*dt; j%bC9UkE3  
    I=1.1*I0*(exp(-9.2e4*t)-exp(-0.5e7*t)); `T1bY9O.  
    for j=1:n0  |xg#Q`O  
        Ez(1,j)=CA*Ez(1,j)+4*CB*Hphi(1,j)/delta-abi*I; j[w=pF,o  
    end 4<E <sD  
Ck%if  
-
Hw3rv3o  
    %迭代求解磁场 sSGXd=":  
    Hphi_last=Hphi; L%\Wt1\[  
    for j=1:n0  #$2
/<  
        for i=1:n0 BaIuOZ@,  
            Hphi(i,j)=CP*Hphi(i,j)+CQ*(Ez(i+1,j)-Ez(i,j))-CQ*(Er(i,j+1)-Er(i,j)); 8qw{e`c  
        end QZ;DZMP  
    end                                                   %磁场循环结束 _ gYj@ %  
> cWE@P  
Q/4ICgo4  
0} &/n>F  
    %Mur边界条件 , .;0xyc  
ohe[rV>EX  
    % r=h %t!S 7UD  
    for j=1:n0 #Y   
        Hphi(n0,j)=C1*Hphi(n0-1,j)+C2*Hphi_last(n0,j)+C3*Hphi_last(n0-1,j); ZADMtsk  
    end k%O3\q  
    %z=0,h >35w"a7S  
    for i=1:n0 V'Kied+  
.. mDbTOtD  
3)?WSOsL:  

未注册仅能浏览部分内容，查看全部内容及附件请先 登录 或 注册