[讨论]《计算电磁学的数值方法》一书中的程序错误

以下是吕华英老师 编写的《计算电磁学的数值方法》一书中的计算微带天线S11程序 TCS^nBEE  
我输入matlab并运行，结果与给出结果相差甚远， 我找出其中部分错误，但是没有完全找到， X]AbBzy  
请各位大侠帮忙看看 还有什么错误 使其调试成功。 5"8R|NU:\0  
B(U0 ~{7a  
clear @AAkEWo)_  
mu=4*pi*1e-7; e0=8.854*1e-12; ]#Q'~X W  
v=3e8;dt=0.441e-12;dx=0.389e-3;dy=0.400e-3;dz=0.265e-3;r=2.2;m=(1+r)/2; "f3KE=cUm  
A=dt/dz/mu;       nO7#m~  
B=dt/dy/mu; )uIHonXU  
C=dt/dx/mu; et0yS%7+?@  
D=dt/(r*dy)/e0; t4CI+fqy  
E=dt/(r*dz)/e0; i=8){GX4  
F=dt/(r*dx)/e0; vV+>JM6<K  
G=dt/(m*dy)/e0; ky98Bz%  
H=dt/(m*dz)/e0; ]z_C7Y"4BR  
I=dt/(m*dx)/e0; rCFTch"  
J=dt/dy/e0; sKuPV
  
K=dt/dz/e0; L}5IX)#gH
 
L=dt/dx/e0; >5}jM5$  
a=(v*dt/sqrt(r)-dx)/(v*dt/sqrt(r)+dx); QW1d&Gb.(  
b=(v*dt/sqrt(m)-dx)/(v*dt/sqrt(m)+dx); w-(^w9_e  
c=(v*dt-dx)/(v*dt+dx); 
H-Z1i  
d=(v*dt/sqrt(r)-dy)/(v*dt/sqrt(r)+dy); @dp1bkU  
e=(v*dt/sqrt(m)-dy)/(v*dt/sqrt(m)+dy); d*TpHLm  
f=(v*dt-dy)/(v*dt+dy); ,1g*0W^  
g=(v*dt-dz)/(v*dt+dz); =| M[JPr  
%输入初始值 )./.rtP|4  
i=2:62;j=2:102;k=2:18; >!|(n@  
Ex1(i,j,k)=0;Ey1(i,j,k)=0;Ez1(i,j,k)=0; [Q*aJLG  
Ex2(i,j,k)=0;Ey2(i,j,k)=0;Ez2(i,j,k)=0; -(YdK8  
Hx1(i,j,k)=0;Hy1(i,j,k)=0;Hz1(i,j,k)=0; sT"h)I)]*  
Hx2(i,j,k)=0;Hy2(i,j,k)=0;Hz2(i,j,k)=0; 1O]27"9  
%时间迭代(200,400,600,800) Il9pL~u  
for n=0:800 60St99@O  
    %计算磁场 /]*#+;;%  
    i=2:62;j=2:101;k=2:17; ~ E|L4E  
        Hx2(i,j,k)=Hx1(i,j,k)+A*(Ey1(i,j,k+1)-Ey1(i,j,k))-B*(Ez1(i,j+1,k)-Ez1(i,j,k)); t?R=a-ZI  
    i=2:61;j=2:102;k=2:17; $>mTPNF  
        Hy2(i,j,k)=Hy1(i,j,k)+C*(Ez1(i+1,j,k)-Ez1(i,j,k))-A*(Ex1(i,j,k+1)-Ex1(i,j,k)); 8GD!]t#  
    i=2:61;j=2:101;k=2:18; POCFT0R}  
        Hz2(i,j,k)=Hz1(i,j,k)+B*(Ex1(i,j+1,k)-Ex1(i,j,k))-C*(Ey1(i+1,j,k)-Ey1(i,j,k)); zO07X*Bw  
    %计算电场 (6Sf#M  
    i=2:61;j=3:101;k=3:4;%介质层 a,U@ !}K  
        Ex2(i,j,k)=Ex1(i,j,k)+D*(Hz2(i,j,k)-Hz2(i,j-1,k))-E*(Hy2(i,j,k)-Hy2(i,j,k-1)); M=aWL!nJ  
    i=3:61;j=2:101;k=3:4; Sk 10"DB/  
        Ey2(i,j,k)=Ey1(i,j,k)+E*(Hx2(i,j,k)-Hx2(i,j,k-1))-F*(Hz2(i,j,k)-Hz2(i-1,j,k)); Z/@%MEU[zl  
    i=3:61;j=3:101;k=2:4; `nDgwp:b"  
        Ez2(i,j,k)=Ez1(i,j,k)+F*(Hy2(i,j,k)-Hy2(i-1,j,k))-D*(Hx2(i,j,k)-Hx2(i,j-1,k)); 1*Ui=M4  
    i=2:61;j=3:101;k=5;%交界面 >{]mN5  
        Ex2(i,j,k)=Ex1(i,j,k)+G*(Hz2(i,j,k)-Hz2(i,j-1,k))-H*(Hy2(i,j,k)-Hy2(i,j,k-1)); C5jR||  
    i=3:61;j=2:101;k=5; r`T(xJ!)  
        Ey2(i,j,k)=Ey1(i,j,k)+H*(Hx2(i,j,k)-Hx2(i,j,k-1))-I*(Hz2(i,j,k)-Hz2(i-1,j,k)); T c{]w?V  
        Ez2(i,j,k)=Ez1(i,j,k)+F*(Hy2(i,j,k)-Hy2(i-1,j,k))-D*(Hx2(i,j,k)-Hx2(i,j-1,k)); n\Y|0\ B  
         c(<,qWH  
    i=2:61;j=3:101;k=6:17;%空气层 Kzd`|+?'`M  
        Ex2(i,j,k)=Ex1(i,j,k)+J*(Hz2(i,j,k)-Hz2(i,j-1,k))-K*(Hy2(i,j,k)-Hy2(i,j,k-1)); 7VwLyy  
    i=3:61;j=2:101;k=6:17; (iZE}qf7g  
        Ey2(i,j,k)=Ey1(i,j,k)+K*(Hx2(i,j,k)-Hx2(i,j,k-1))-L*(Hz2(i,j,k)-Hz2(i-1,j,k)); X@ Gm:6  
    i=3:61;j=3:101;k=5:17; \PB~6  
        Ez2(i,j,k)=Ez1(i,j,k)+L*(Hy2(i,j,k)-Hy2(i-1,j,k))-J*(Hx2(i,j,k)-Hx2(i,j-1,k)); !B_?_ a  
    %边界1--接地板 [ZP8l'?  
    i=2:61;j=2:102; Ck0R%|  
        Ex2(i,j,2)=0;  Hz2(i,j,2)=0; Rtl;*ZAS  
    i=2:62;j=2:101; =7c1l77z  
        Ey2(i,j,2)=0;  Hz2(i,j,2)=0; Oy!j`  
    %边界2—微带线贴片
~CB6+t>  
    i=21:27;j=2:52; #j{!&4M  
    Ex2(i,j,5)=0;Ey2(i,j,5)=0;  Hz2(i,j,2)=0; 2W=( {e)$  
    %----矩形切片 ZP&"[_  
    i=16:48;j=52:92; = ?hx+-'  
    Ex2(i,j,5)=0;Ey2(i,j,5)=0;  Hz2(i,j,2)=0; $N#f)8v  
    %边界3--左侧吸收边界 G;%Pf9o26  
    j=2:101;k=3:4; K$..#]\TM  
        Ey2(2,j,k)=Ey1(3,j,k)+a*(Ey2(3,j,k)-Ey1(2,j,k)); S6
sw)  
    j=3:101;k=2:4; "A_WU|  
        Ez2(2,j,k)=Ez1(3,j,k)+a*(Ez2(3,j,k)-Ez1(2,j,k)); c=T^)~$$  
    j=2:101; O}!L;?  
        Ey2(2,j,5)=Ey1(3,j,5)+b*(Ey2(3,j,5)-Ey1(2,j,5)); &)p/cOiV  
    j=2:101;k=6:17; ;=.QT  
        Ey2(2,j,k)=Ey1(3,j,k)+c*(Ey2(3,j,k)-Ey1(2,j,k)); $I7/FZP  
    j=3:101;k=5:17; NgnHo\)  
        Ez2(2,j,k)=Ez1(3,j,k)+c*(Ez2(3,j,k)-Ez1(2,j,k)); Svmyg]  
    %边界4--右侧吸收边界 5M{DJ/q  
    j=2:101;k=3:4; t;@VsQ8  
        Ey2(62,j,k)=Ey1(61,j,k)+a*(Ey2(61,j,k)-Ey1(62,j,k)); wxg`[c$:  
    j=3:101;k=2:4;
@:~O  
        Ez2(62,j,k)=Ez1(61,j,k)+a*(Ez2(61,j,k)-Ez1(62,j,k)); Q\~4J1  
    j=2:101; &!{wbm@  
        Ey2(62,j,5)=Ey1(61,j,5)+b*(Ey2(61,j,5)-Ey1(62,j,5)); <M1XG7_I  
    j=2:101;k=6:17; 2>l:: 8Pp  
       Ey2(62,j,k)=Ey1(61,j,k)+c*(Ey2(61,j,k)-Ey1(62,j,k)); wOy1i/oj  
    j=3:101;k=5:17;   1;l&ck-Gg/  
       Ez2(62,j,k)=Ez1(61,j,k)+c*(Ez2(61,j,k)-Ez1(62,j,k)); nJ|8#U7  
    %边界5--后侧吸收边界 y/m^G=Q6g#  
    i=2:61;k=3:4; 61Nj&1Ze  
       Ex2(i,102,k)=Ex1(i,101,k)+d*(Ex2(i,101,k)-Ex1(i,102,k)); q(Y<cJ?X  
    i=3:61;k=2:4; h 
!~u9  
       Ez2(i,102,k)=Ez1(i,101,k)+d*(Ez2(i,101,k)-Ez1(i,102,k)); 7FVu[Qu  
    i=2:61; x hFQjV?V  
       Ex2(i,102,5)=Ex1(i,101,5)+e*(Ex2(i,101,5)-Ex1(i,102,5)); lq$1CI  
    i=2:61;k=6:17; =Po!\[SBU  
       Ex2(i,102,k)=Ex1(i,101,k)+f*(Ex2(i,101,k)-Ex1(i,102,k)); ogX'3L  
    i=3:61;k=5:17; C'hI{4@P  
       Ez2(i,102,k)=Ez1(i,101,k)+f*(Ez2(i,101,k)-Ez1(i,102,k)); r'p;Nj.  
    %边界6--上侧吸收边界 /AR;O4X+  
    i=2:61;j=3:101; jG0{>P#+  
       Ex2(i,j,18)=Ex1(i,j,17)+g*(Ex2(i,j,17)-Ex1(i,j,18)); S< <xlW  
    i=3:61;j=2:101; K'%,dn  
       Ey2(i,j,18)=Ey1(i,j,17)+g*(Ey2(i,j,17)-Ey1(i,j,18)); 9_rNJLj8y  
    %边界7--z方向的棱 W.ud<OKP90  
    k=2:17; d{f3R8~Q.  
       Ez2(2,2,k)=Ez2(2,3,k)+Ez2(3,2,k)-Ez2(3,3,k); HB4Hz0Fa  
       Ez2(62,2,k)=Ez2(62,3,k)+Ez2(61,2,k)-Ez2(61,3,k); g]V}azLr  
       Ez2(2,102,k)=Ez2(2,101,k)+Ez2(3,102,k)-Ez2(3,101,k); h@72eav3+  
       Ez2(62,102,k)=Ez2(62,101,k)+Ez2(61,102,k)-Ez2(61,101,k); EO,;^RtB  
    %边界8--x方向的棱 >,`/ z  
    i=2:61; 0C}7=_?  
       Ex2(i,2,18)=Ex2(i,3,18)+Ex2(i,2,17)-Ex2(i,3,17); W9?Yzl  
       Ex2(i,102,18)=Ex2(i,101,18)+Ex2(i,102,17)-Ex2(i,101,17); ^b`}g  
    %边界9--y方向的棱 ,XW6W&vR;  
    j=2:101; }B_n}<tjD  
       Ey2(2,j,18)=Ey2(3,j,18)+Ey2(2,j,17)-Ey2(3,j,17); 1WPDMLuN  
       Ey2(62,j,18)=Ey2(61,j,18)+Ey2(62,j,17)-Ey2(61,j,17); #(jozl_8  
    %边界9--源平面 '!!w|kd  
     if (n>220), ,sk;|OAI  
         i=2:61;k=3:4; C9
/?B:  
           Ex2(i,2,k)=Ex1(i,3,k)+d*(Ex2(i,3,k)-Ex1(i,2,k)); .BXZ\r`  
         i=3:61;k=2:4; w.cQ|_  
           Ez2(i,2,k)=Ez1(i,3,k)+d*(Ez2(i,3,k)-Ez1(i,2,k)); ?kew[oZ  
         i=2:61; Ho?+?YJ#P  
           Ex2(i,2,5)=Ex1(i,3,5)+e*(Ex2(i,3,5)-Ex1(i,2,5)); 8 F'i5i  
         i=2:61;k=6:17; P sD+?  
           Ex2(i,2,k)=Ex1(i,3,k)+f*(Ex2(i,3,k)-Ex1(i,2,k)); L;xc,"\3  
         i=3:61;k=5:17; ,LXuU8sB  
           Ez2(i,2,k)=Ez1(i,3,k)+f*(Ez2(i,3,k)-Ez1(i,2,k)); XQ<2(}]4  
     else Xu$xO(  
         i=2:61;k=3:4; A^JeB<, 5a  
           Ex2(i,2,k)=Ex1(i,2,k)+D*(2*Hz2(i,2,k))-E*(Hy2(i,2,k)-Hy2(i,2,k-1)); ^6+P&MxM  
         i=2:61; 56*}}B$?  
           Ex2(i,2,5)=Ex1(i,2,5)+G*(2*Hz2(i,2,5))-H*(Hy2(i,2,5)-Hy2(i,2,4)); YizJT0$  
         i=2:61;k=6:17; \qAMs^1-  
           Ex2(i,2,k)=Ex1(i,2,k)+J*(2*Hz2(i,2,k))-K*(Hy2(i,2,k)-Hy2(i,2,k-1)); RC8{QgaI  
       &nb .. vZC2F  
F]W'spF,  

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

建议楼主把程序打包之后发上来，不然拷贝到MATLAB之后会出现乱码。这样别人检查你的错误会很麻烦。

好的  谢谢

以上程序在以下文件中

有没有那位大侠 看出问题了？  讨论讨论一下

我仔细看过这个，书中有数个错误点。 C2e.RTxc  
最重要的是计算S11的方法错了。但是也不知道如何在它上面改。

哎···都是错的··难道就没有一个正确的程序让我这样的新手参考下么····郁闷···