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

小弟刚学fdtd不久，照葛德彪那本书天线那张做了个二维柱坐标系下的TE波线电流源，激励源采用双指数函数，吸收边界采用一阶Mur #UP,;W  
感觉激励源的加入还是有很大问题 `UzCq06rJ1  
另外步长改为1的话出来的波形明显震荡 F~11 _  
请各位高手指点一下 'Gm!Jblo@  
RMs1{64:  
tPMgZ  
j[Hg]  
1{Sx V  
clear all jHUz`.8B  
close all btkMY<o7  
clc .1l[
l5$  
'RN"yMv7l  
d-_93
 
%-------------------------------------------------- GtGyY0  
%FDTD法模拟柱坐标二维TE波的传播(Mur吸收边界条件) t's5~  
%-------------------------------------------------- bN<
O<x1j  
\dQ2[Ek  
Jf8'N ot  
%初始化 ^O,6(@>  
c=3e8;                                                       %光速 ]2u7?l  
delta=10; 8tB{rK,  
% delta=1e-3;                                                  %空间步长 _YD<Q@  
dt=delta/(2*c);                                              %时间步长 .*clY  
u0=(4*pi)*1e-7;                                              %磁导率 .WPuQZ!  
e0=(1/(36*pi))*1e-9;                                         %介电常数 QvyUd%e'5A  
I0=15e3;                                                     %基电流 [jG uO%  
n0=200;                                                      %空间网格数 a'L7y%  
% n1=33;                                                       %源的位置 ' ]H#0.  
abi=4*dt/(pi*e0*delta^2);                                    %激励源系数 f{oxF?|89  
CA=1; j+ T\c2d  
CB=dt/(e0*delta);                                            %Er与Ez中的系数 _^,[
wD  
CP=1; cmC&s'/8`D  
CQ=dt/(u0*delta);                                            %Hphi中的系数 d,E/9y\e  
C1=(1/(2*delta)-1/(2*c*dt)-1/(8*delta*(n0)))/(1/(2*delta)+1/(2*c*dt)+1/(8*delta*(n0)));                              %一阶边界条件中的系数 kB!M[[t  
C2=(-1/(2*delta)+1/(2*c*dt)-1/(8*delta*(n0)))/(1/(2*delta)+1/(2*c*dt)+1/(8*delta*(n0))); s^x ,S  
C3=(1/(2*delta)+1/(2*c*dt)-1/(8*delta*(n0)))/(1/(2*delta)+1/(2*c*dt)+1/(8*delta*(n0))); *jqPKK/  
% C1=-(4*n0+1)/(12*n0+1);                              %一阶边界条件中的系数 '!2  
% C2=-(-4*n0+1)/(12*n0+1); q0xE&[C[M  
% C3=-(-12*n0+1)/(12*n0+1); e>s.mH6A  
% C1=-(4*(n0+1)+1)/(12*(n0+1)+1);                              %一阶边界条件中的系数 u6CMRZ$  
% C2=-(-4*(n0+1)+1)/(12*(n0+1)+1); zv3<i (  
% C3=-(-4*(n0+1)+1)/(-12*(n0+1)+1); ABS BtH ?  
C4=(1.3e8*dt-delta)/(1.3e8*dt+delta); 1tK6lrhj  
Er=zeros(n0+1,n0+1); =V4_DJ(&  
Ez=zeros(n0+1,n0+1); 8~XI7g'5x  
Hphi=zeros(n0,n0); TQyFF/K  
Hphi_last=zeros(n0,n0);                                 %大型矩阵维数的预先确定 B3i=pcef  
RdpQJ)3F  
62KW HB9S  
okW)s*7  
IcNIuv  
for k=1:2000 ,_7tRkn  
hJL0M!  
3hpz.ISk  
    %迭代求电场 Fb`7aFIf  
)SO1P6  
.LVOaxT  
> & lg  
    %求Er #'{PYr  
    for j=2:n0 l2_E6U"  
        for i=1:n0 EEnTq  
            Er(i,j)=CA*Er(i,j)-CB*(Hphi(i,j)-Hphi(i,j-1)); (Su2\x  
        end 0yEyt7 ~@  
    end                                                %Er循环结束 )a2m<"  
    %求Ez "}Sid+)<  
    for j=1:n0 */@bNT9BgO  
        for i=2:n0 K,T]Fuy  
            Ez(i,j)=CA*Ez(i,j)+CB*(delta*i*Hphi(i,j)-delta*(i-1)*Hphi(i-1,j))/(delta*(i-0.5)); '}[L sU  
        end /$N#_Xblr  
    end                                                %Ez循环结束 k?*DBXJv  
    %激励源 rD)v%vvr&`  
    t=k*dt; bJ5z??  
    I=1.1*I0*(exp(-9.2e4*t)-exp(-0.5e7*t)); ur_"m+  
    for j=1:n0 b|C,b"$N0  
        Ez(1,j)=CA*Ez(1,j)+4*CB*Hphi(1,j)/delta-abi*I; x'PjP1  
    end N4JL.(m){I  
za 4B+&JJ  
7|?@\ZE  
    %迭代求解磁场 x#xO {  
    Hphi_last=Hphi; GA@Q:n8UuR  
    for j=1:n0 iPi'5g(a   
        for i=1:n0 "r(pK@h  
            Hphi(i,j)=CP*Hphi(i,j)+CQ*(Ez(i+1,j)-Ez(i,j))-CQ*(Er(i,j+1)-Er(i,j)); p!5JO4F$  
        end D +%k1  
    end                                                   %磁场循环结束 [(ib9_`A'1  
Q,+*u%/u  
T<~[vjA  
IF21T  
    %Mur边界条件 G"R>aw  
tfU3 6PR  
    % r=h T;e(Q,!H  
    for j=1:n0 ${H&Q*  
        Hphi(n0,j)=C1*Hphi(n0-1,j)+C2*Hphi_last(n0,j)+C3*Hphi_last(n0-1,j); pz
p"NKxi  
    end BN>$LL  
    %z=0,h $@!&ML  
    for i=1:n0 /bj`%Q.n  
.. NNrZb?  
AZ@
Zo'  

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