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

以下是吕华英老师 编写的《计算电磁学的数值方法》一书中的计算微带天线S11程序 2L'vB1`  
我输入matlab并运行，结果与给出结果相差甚远， 我找出其中部分错误，但是没有完全找到， ~;pv&s5}  
请各位大侠帮忙看看 还有什么错误 使其调试成功。 .7!n%Ks  
Xg](V.B6  
clear RnA>oKc  
mu=4*pi*1e-7; e0=8.854*1e-12; :tv:46+s=  
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; 3e<FlH{  
A=dt/dz/mu;       |#r[{2sS  
B=dt/dy/mu; ?&H1C4   
C=dt/dx/mu; -RSPYQjz  
D=dt/(r*dy)/e0; {u1Rc/Lw  
E=dt/(r*dz)/e0; m _0D^e7#  
F=dt/(r*dx)/e0; +:s]>R eDa  
G=dt/(m*dy)/e0; jf_0I
E  
H=dt/(m*dz)/e0; 7H1 ii   
I=dt/(m*dx)/e0; N<xf=a+j  
J=dt/dy/e0; 5K~kzRL$r  
K=dt/dz/e0; UP@a ?w  
L=dt/dx/e0; jmcb-=ts  
a=(v*dt/sqrt(r)-dx)/(v*dt/sqrt(r)+dx); X CjYm  
b=(v*dt/sqrt(m)-dx)/(v*dt/sqrt(m)+dx); HhmC
+3w.7  
c=(v*dt-dx)/(v*dt+dx);  QTN _Z#'  
d=(v*dt/sqrt(r)-dy)/(v*dt/sqrt(r)+dy); b :Knc$  
e=(v*dt/sqrt(m)-dy)/(v*dt/sqrt(m)+dy); V ifQ@  
f=(v*dt-dy)/(v*dt+dy); k#4%d1O}  
g=(v*dt-dz)/(v*dt+dz); 2hjR'6h"Y  
%输入初始值 :oH~{EQ  
i=2:62;j=2:102;k=2:18; >UP{=`  
Ex1(i,j,k)=0;Ey1(i,j,k)=0;Ez1(i,j,k)=0; ?H c~ 3  
Ex2(i,j,k)=0;Ey2(i,j,k)=0;Ez2(i,j,k)=0; *mc]Oa  
Hx1(i,j,k)=0;Hy1(i,j,k)=0;Hz1(i,j,k)=0; S<
hj6A  
Hx2(i,j,k)=0;Hy2(i,j,k)=0;Hz2(i,j,k)=0; T@n-^B!Xq  
%时间迭代(200,400,600,800) ;FO1b*  
for n=0:800 C3Hq&TVf/  
    %计算磁场 
S_z}h  
    i=2:62;j=2:101;k=2:17; ?D].Za^km  
        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)); <eObQ[mQ  
    i=2:61;j=2:102;k=2:17; ?:Y0#Btj  
        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)); \~(ww3e  
    i=2:61;j=2:101;k=2:18; m"2KAq61  
        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)); g[rxKn\Z  
    %计算电场 iaqhP7!  
    i=2:61;j=3:101;k=3:4;%介质层 h48JpZ"  
        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)); Wp0e?bK_  
    i=3:61;j=2:101;k=3:4; {\OIowa  
        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)); X[frL)k]  
    i=3:61;j=3:101;k=2:4; n\nC.|_G@  
        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)); 4{:W5eT!/  
    i=2:61;j=3:101;k=5;%交界面 `Z3Qx~fx  
        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)); e7{n=M  
    i=3:61;j=2:101;k=5; I[~EQ{Iz  
        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)); Q]';1#J\  
        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)); AnbY<&OC1  
         ZWC-<QO"<  
    i=2:61;j=3:101;k=6:17;%空气层 6\MJv
g\;  
        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)); X(-e-:B4;  
    i=3:61;j=2:101;k=6:17; |P_\l,f8`  
        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)); xZ51iD$  
    i=3:61;j=3:101;k=5:17; z2m%L0  
        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)); cT8`l!RD<  
    %边界1--接地板 t|gEMDGa3  
    i=2:61;j=2:102; @quNVx(y  
        Ex2(i,j,2)=0;  Hz2(i,j,2)=0; +F &,,s"&  
    i=2:62;j=2:101; \haJe~  
        Ey2(i,j,2)=0;  Hz2(i,j,2)=0; q;T3bxp+  
    %边界2—微带线贴片 vrtK~5K  
    i=21:27;j=2:52; &&0,;r,-)  
    Ex2(i,j,5)=0;Ey2(i,j,5)=0;  Hz2(i,j,2)=0; k)W8%=R  
    %----矩形切片 [<)/ c>Y  
    i=16:48;j=52:92; VrF(0,-Z`3  
    Ex2(i,j,5)=0;Ey2(i,j,5)=0;  Hz2(i,j,2)=0; 05(lh<C  
    %边界3--左侧吸收边界 ?|8QL9Q"|  
    j=2:101;k=3:4; C+r<DC3  
        Ey2(2,j,k)=Ey1(3,j,k)+a*(Ey2(3,j,k)-Ey1(2,j,k)); x?J- {6k  
    j=3:101;k=2:4; 5Rv6+d  
        Ez2(2,j,k)=Ez1(3,j,k)+a*(Ez2(3,j,k)-Ez1(2,j,k)); >IW0YIQy,  
    j=2:101; {iP^51fy  
        Ey2(2,j,5)=Ey1(3,j,5)+b*(Ey2(3,j,5)-Ey1(2,j,5)); apxZ}  
    j=2:101;k=6:17; U9D4bn D  
        Ey2(2,j,k)=Ey1(3,j,k)+c*(Ey2(3,j,k)-Ey1(2,j,k)); `<Z5/;a5W  
    j=3:101;k=5:17; UFxQ-GV4  
        Ez2(2,j,k)=Ez1(3,j,k)+c*(Ez2(3,j,k)-Ez1(2,j,k)); <a-I-~  
    %边界4--右侧吸收边界 tylMJ$ 9*.  
    j=2:101;k=3:4; x%ZgLvdp,  
        Ey2(62,j,k)=Ey1(61,j,k)+a*(Ey2(61,j,k)-Ey1(62,j,k)); U!:Q|':=h  
    j=3:101;k=2:4; ckqU2ETpD}  
        Ez2(62,j,k)=Ez1(61,j,k)+a*(Ez2(61,j,k)-Ez1(62,j,k)); G?LPj*=$?  
    j=2:101; Hno:"k?  
        Ey2(62,j,5)=Ey1(61,j,5)+b*(Ey2(61,j,5)-Ey1(62,j,5)); u,:GJU  
    j=2:101;k=6:17; h[D"O6 y  
       Ey2(62,j,k)=Ey1(61,j,k)+c*(Ey2(61,j,k)-Ey1(62,j,k)); v535LwFW  
    j=3:101;k=5:17;   ^h+<Q%'a'  
       Ez2(62,j,k)=Ez1(61,j,k)+c*(Ez2(61,j,k)-Ez1(62,j,k)); JM5w`=  
    %边界5--后侧吸收边界 _y}]j;e8>{  
    i=2:61;k=3:4; 6#xP[hlR[  
       Ex2(i,102,k)=Ex1(i,101,k)+d*(Ex2(i,101,k)-Ex1(i,102,k)); Q 'R@'W9  
    i=3:61;k=2:4; u>Z0u
g6x  
       Ez2(i,102,k)=Ez1(i,101,k)+d*(Ez2(i,101,k)-Ez1(i,102,k)); IqK??KSC  
    i=2:61; 3K2`1+kBVG  
       Ex2(i,102,5)=Ex1(i,101,5)+e*(Ex2(i,101,5)-Ex1(i,102,5)); *P_ 3A:_  
    i=2:61;k=6:17; eRC /Pr  
       Ex2(i,102,k)=Ex1(i,101,k)+f*(Ex2(i,101,k)-Ex1(i,102,k)); $_-f}E  
    i=3:61;k=5:17; ~X^L3=!vf  
       Ez2(i,102,k)=Ez1(i,101,k)+f*(Ez2(i,101,k)-Ez1(i,102,k)); kji*7a?y  
    %边界6--上侧吸收边界 .Od.lxz"mp  
    i=2:61;j=3:101; AL/q6PWi  
       Ex2(i,j,18)=Ex1(i,j,17)+g*(Ex2(i,j,17)-Ex1(i,j,18)); ^CUeq"GYoZ  
    i=3:61;j=2:101; .6%-Il  
       Ey2(i,j,18)=Ey1(i,j,17)+g*(Ey2(i,j,17)-Ey1(i,j,18)); j2T Z`Z?a^  
    %边界7--z方向的棱 @[0zZX2EE  
    k=2:17; >9{Gdq[gyr  
       Ez2(2,2,k)=Ez2(2,3,k)+Ez2(3,2,k)-Ez2(3,3,k); !>B|z=  
       Ez2(62,2,k)=Ez2(62,3,k)+Ez2(61,2,k)-Ez2(61,3,k); 1F*gPhm  
       Ez2(2,102,k)=Ez2(2,101,k)+Ez2(3,102,k)-Ez2(3,101,k); }&d@6m]  
       Ez2(62,102,k)=Ez2(62,101,k)+Ez2(61,102,k)-Ez2(61,101,k); hKw4[wB]  
    %边界8--x方向的棱 8(UUc>g  
    i=2:61; 4P@Ak7iL(V  
       Ex2(i,2,18)=Ex2(i,3,18)+Ex2(i,2,17)-Ex2(i,3,17); ':8yp|A|  
       Ex2(i,102,18)=Ex2(i,101,18)+Ex2(i,102,17)-Ex2(i,101,17); 8\m_.e  
    %边界9--y方向的棱 B$aA=+<S  
    j=2:101; #C1u~db  
       Ey2(2,j,18)=Ey2(3,j,18)+Ey2(2,j,17)-Ey2(3,j,17); eK\1cs  
       Ey2(62,j,18)=Ey2(61,j,18)+Ey2(62,j,17)-Ey2(61,j,17); 2f1WT g)  
    %边界9--源平面 Vx@JP93|  
     if (n>220), +K4d(!Sb  
         i=2:61;k=3:4; #%U5,[<a8  
           Ex2(i,2,k)=Ex1(i,3,k)+d*(Ex2(i,3,k)-Ex1(i,2,k)); -W(O~AK  
         i=3:61;k=2:4; WL4{_X  
           Ez2(i,2,k)=Ez1(i,3,k)+d*(Ez2(i,3,k)-Ez1(i,2,k)); \2#>@6Sqrl  
         i=2:61; `;-K
/)/x  
           Ex2(i,2,5)=Ex1(i,3,5)+e*(Ex2(i,3,5)-Ex1(i,2,5)); MXY[t  
         i=2:61;k=6:17;
ANEW^\  
           Ex2(i,2,k)=Ex1(i,3,k)+f*(Ex2(i,3,k)-Ex1(i,2,k)); 'qS&7 W(  
         i=3:61;k=5:17; 3]BK*OqJ  
           Ez2(i,2,k)=Ez1(i,3,k)+f*(Ez2(i,3,k)-Ez1(i,2,k)); V`Z-m-V~1  
     else w-?_U7'  
         i=2:61;k=3:4; @b\/\\{  
           Ex2(i,2,k)=Ex1(i,2,k)+D*(2*Hz2(i,2,k))-E*(Hy2(i,2,k)-Hy2(i,2,k-1)); n7`R+4/s  
         i=2:61; umrfA  
           Ex2(i,2,5)=Ex1(i,2,5)+G*(2*Hz2(i,2,5))-H*(Hy2(i,2,5)-Hy2(i,2,4)); <rc?EV  
         i=2:61;k=6:17; oLEqy  
           Ex2(i,2,k)=Ex1(i,2,k)+J*(2*Hz2(i,2,k))-K*(Hy2(i,2,k)-Hy2(i,2,k-1)); OD!b*Iy|  
       &nb .. 2xvTijO0  
Jrd:6
Z  

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

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

好的  谢谢

以上程序在以下文件中

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

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

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