2983阅读
6回复

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

UID ：17376

注册：2008-08-31
登录：2015-03-22
发帖：45
等级：仿真新人

0楼发表于: 2010-12-06 21:28:27

以下是吕华英老师编写的《计算电磁学的数值方法》一书中的计算微带天线S11程序
我输入matlab并运行，结果与给出结果相差甚远，我找出其中部分错误，但是没有完全找到，
请各位大侠帮忙看看还有什么错误使其调试成功。

clear
mu=4*pi*1e-7; e0=8.854*1e-12;
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;
A=dt/dz/mu;
B=dt/dy/mu;
C=dt/dx/mu;
D=dt/(r*dy)/e0;
E=dt/(r*dz)/e0;
F=dt/(r*dx)/e0;
G=dt/(m*dy)/e0;
H=dt/(m*dz)/e0;
I=dt/(m*dx)/e0;
J=dt/dy/e0;
K=dt/dz/e0;
L=dt/dx/e0;
a=(v*dt/sqrt(r)-dx)/(v*dt/sqrt(r)+dx);
b=(v*dt/sqrt(m)-dx)/(v*dt/sqrt(m)+dx);
c=(v*dt-dx)/(v*dt+dx);
d=(v*dt/sqrt(r)-dy)/(v*dt/sqrt(r)+dy);
e=(v*dt/sqrt(m)-dy)/(v*dt/sqrt(m)+dy);
f=(v*dt-dy)/(v*dt+dy);
g=(v*dt-dz)/(v*dt+dz);
%输入初始值
i=2:62;j=2:102;k=2:18;
Ex1(i,j,k)=0;Ey1(i,j,k)=0;Ez1(i,j,k)=0;
Ex2(i,j,k)=0;Ey2(i,j,k)=0;Ez2(i,j,k)=0;
Hx1(i,j,k)=0;Hy1(i,j,k)=0;Hz1(i,j,k)=0;
Hx2(i,j,k)=0;Hy2(i,j,k)=0;Hz2(i,j,k)=0;
%时间迭代(200,400,600,800)
for n=0:800
    %计算磁场
    i=2:62;j=2:101;k=2:17;
        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));
    i=2:61;j=2:102;k=2:17;
        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));
    i=2:61;j=2:101;k=2:18;
        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));
    %计算电场
    i=2:61;j=3:101;k=3:4;%介质层
        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));
    i=3:61;j=2:101;k=3:4;
        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));
    i=3:61;j=3:101;k=2:4;
        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));
    i=2:61;j=3:101;k=5;%交界面
        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));
    i=3:61;j=2:101;k=5;
        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));
        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));

    i=2:61;j=3:101;k=6:17;%空气层
        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));
    i=3:61;j=2:101;k=6:17;
        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));
    i=3:61;j=3:101;k=5:17;
        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));
    %边界1--接地板
    i=2:61;j=2:102;
        Ex2(i,j,2)=0;  Hz2(i,j,2)=0;
    i=2:62;j=2:101;
        Ey2(i,j,2)=0;  Hz2(i,j,2)=0;
    %边界2—微带线贴片
    i=21:27;j=2:52;
    Ex2(i,j,5)=0;Ey2(i,j,5)=0;  Hz2(i,j,2)=0;
    %----矩形切片
    i=16:48;j=52:92;
    Ex2(i,j,5)=0;Ey2(i,j,5)=0;  Hz2(i,j,2)=0;
    %边界3--左侧吸收边界
    j=2:101;k=3:4;
        Ey2(2,j,k)=Ey1(3,j,k)+a*(Ey2(3,j,k)-Ey1(2,j,k));
    j=3:101;k=2:4;
        Ez2(2,j,k)=Ez1(3,j,k)+a*(Ez2(3,j,k)-Ez1(2,j,k));
    j=2:101;
        Ey2(2,j,5)=Ey1(3,j,5)+b*(Ey2(3,j,5)-Ey1(2,j,5));
    j=2:101;k=6:17;
        Ey2(2,j,k)=Ey1(3,j,k)+c*(Ey2(3,j,k)-Ey1(2,j,k));
    j=3:101;k=5:17;
        Ez2(2,j,k)=Ez1(3,j,k)+c*(Ez2(3,j,k)-Ez1(2,j,k));
    %边界4--右侧吸收边界
    j=2:101;k=3:4;
        Ey2(62,j,k)=Ey1(61,j,k)+a*(Ey2(61,j,k)-Ey1(62,j,k));
    j=3:101;k=2:4;
        Ez2(62,j,k)=Ez1(61,j,k)+a*(Ez2(61,j,k)-Ez1(62,j,k));
    j=2:101;
        Ey2(62,j,5)=Ey1(61,j,5)+b*(Ey2(61,j,5)-Ey1(62,j,5));
    j=2:101;k=6:17;
       Ey2(62,j,k)=Ey1(61,j,k)+c*(Ey2(61,j,k)-Ey1(62,j,k));
    j=3:101;k=5:17;
       Ez2(62,j,k)=Ez1(61,j,k)+c*(Ez2(61,j,k)-Ez1(62,j,k));
    %边界5--后侧吸收边界
    i=2:61;k=3:4;
       Ex2(i,102,k)=Ex1(i,101,k)+d*(Ex2(i,101,k)-Ex1(i,102,k));
    i=3:61;k=2:4;
       Ez2(i,102,k)=Ez1(i,101,k)+d*(Ez2(i,101,k)-Ez1(i,102,k));
    i=2:61;
       Ex2(i,102,5)=Ex1(i,101,5)+e*(Ex2(i,101,5)-Ex1(i,102,5));
    i=2:61;k=6:17;
       Ex2(i,102,k)=Ex1(i,101,k)+f*(Ex2(i,101,k)-Ex1(i,102,k));
    i=3:61;k=5:17;
       Ez2(i,102,k)=Ez1(i,101,k)+f*(Ez2(i,101,k)-Ez1(i,102,k));
    %边界6--上侧吸收边界
    i=2:61;j=3:101;
       Ex2(i,j,18)=Ex1(i,j,17)+g*(Ex2(i,j,17)-Ex1(i,j,18));
    i=3:61;j=2:101;
       Ey2(i,j,18)=Ey1(i,j,17)+g*(Ey2(i,j,17)-Ey1(i,j,18));
    %边界7--z方向的棱
    k=2:17;
       Ez2(2,2,k)=Ez2(2,3,k)+Ez2(3,2,k)-Ez2(3,3,k);
       Ez2(62,2,k)=Ez2(62,3,k)+Ez2(61,2,k)-Ez2(61,3,k);
       Ez2(2,102,k)=Ez2(2,101,k)+Ez2(3,102,k)-Ez2(3,101,k);
       Ez2(62,102,k)=Ez2(62,101,k)+Ez2(61,102,k)-Ez2(61,101,k);
    %边界8--x方向的棱
    i=2:61;
       Ex2(i,2,18)=Ex2(i,3,18)+Ex2(i,2,17)-Ex2(i,3,17);
       Ex2(i,102,18)=Ex2(i,101,18)+Ex2(i,102,17)-Ex2(i,101,17);
    %边界9--y方向的棱
    j=2:101;
       Ey2(2,j,18)=Ey2(3,j,18)+Ey2(2,j,17)-Ey2(3,j,17);
       Ey2(62,j,18)=Ey2(61,j,18)+Ey2(62,j,17)-Ey2(61,j,17);
    %边界9--源平面
     if (n>220),
         i=2:61;k=3:4;
           Ex2(i,2,k)=Ex1(i,3,k)+d*(Ex2(i,3,k)-Ex1(i,2,k));
         i=3:61;k=2:4;
           Ez2(i,2,k)=Ez1(i,3,k)+d*(Ez2(i,3,k)-Ez1(i,2,k));
         i=2:61;
           Ex2(i,2,5)=Ex1(i,3,5)+e*(Ex2(i,3,5)-Ex1(i,2,5));
         i=2:61;k=6:17;
           Ex2(i,2,k)=Ex1(i,3,k)+f*(Ex2(i,3,k)-Ex1(i,2,k));
         i=3:61;k=5:17;
           Ez2(i,2,k)=Ez1(i,3,k)+f*(Ez2(i,3,k)-Ez1(i,2,k));
     else
         i=2:61;k=3:4;
           Ex2(i,2,k)=Ex1(i,2,k)+D*(2*Hz2(i,2,k))-E*(Hy2(i,2,k)-Hy2(i,2,k-1));
         i=2:61;
           Ex2(i,2,5)=Ex1(i,2,5)+G*(2*Hz2(i,2,5))-H*(Hy2(i,2,5)-Hy2(i,2,4));
         i=2:61;k=6:17;
           Ex2(i,2,k)=Ex1(i,2,k)+J*(2*Hz2(i,2,k))-K*(Hy2(i,2,k)-Hy2(i,2,k-1));
       &nb ..

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

共1条评分

cc81372365

rf币 +1

积极参与论坛交流，欢迎继续参与本贴交流！

2010-12-06

离线cc81372365

UID ：27752