687阅读
2回复

[求助]求救高手指导

UID ：41338

注册：2009-09-10
登录：2010-11-04
发帖：125
等级：仿真二级

0楼发表于: 2009-09-21 17:59:50

这是我写的一个FDTD的程序（MATLAB，参考葛得彪、闫玉波的书籍）中关于UPML的部分的，我实在是不知道该如何去修改才能免去发散的情况了，
只好来这里求救，希望能有高手来救命了！！！
感激不尽啊
另外，附件附上我的整个的程序，有感兴趣的朋友请多指教

sigmax=5/(150*pi*dx*sqrt(1));%暂时这样设计
epsrmax=1.75;
d=10;                                             %吸收UPML层厚度
sig=zeros(ie,je,ke,3);                    %最后的3分别用来表示x、y、z分量
epsr=ones(ie,je,ke,3);
%%%%%%%先是构造PML层的本构参数
epsr(itmin-1-d:itmin-1,jtmin-1-d:jtmax+1+d,ktmin-1-d:ktmax+1+d,2:3)=1;
epsr(itmax+1:itmax+1+d,jtmin-1-d:jtmax+1+d,ktmin-1-d:ktmax+1+d,2:3)=1;
for p=itmin-1-d:itmin-1                                                  %垂直于x轴的sigma,epsr迭代
    sig(p,jtmin-1-d:jtmax+1+d,ktmin-1-d:ktmax+1+d,1)=sigmax*(p-itmin)^4/d^4;
    epsr(p,jtmin-1-d:jtmax+1+d,ktmin-1-d:ktmax+1+d,1)=1+(epsrmax-1)*(p-itmin)^4/d^4;
end
for p=itmax+1:itmax+1+d
sig(p,jtmin-1-d:jtmax+1+d,ktmin-1-d:ktmax+1+d,1)=sigmax*(p-itmax)^4/d^4;
    epsr(p,jtmin-1-d:jtmax+1+d,ktmin-1-d:ktmax+1+d,1)=1+(epsrmax-1)*(p-itmax)^4/d^4;
end

epsr(itmin-1-d:itmax+1+d,jtmin-1-d:jtmax+1+d,ktmin-1-d:ktmin-1,1:2)=1;    %垂直于z轴
epsr(itmin-1-d:itmax+1+d,jtmin-1-d:jtmax+1+d,ktmax+1:ktmax+1+d,1:2)=1;
for k=ktmin-1-d:ktmin-1
    sig(itmin-1-d:itmax+1+d,jtmin-1-d:jtmax+1+d,k,3)=sigmax*(k-ktmin)^4/d^4;
    epsr(itmin-1-d:itmax+1+d,jtmin-1-d:jtmax+1+d,k,3)=1+(epsrmax-1)*(k-ktmin)^4/d^4;
end
for k=ktmax+1:ktmax+1+d
    sig(itmin-1-d:itmax+1+d,jtmin-1-d:jtmax+1+d,k,3)=sigmax*(k-ktmax)^4/d^4;
    epsr(itmin-1-d:itmax+1+d,jtmin-1-d:jtmax+1+d,k,3)=1+(epsrmax-1)*(k-ktmax)^4/d^4;
end

epsr(itmin-1-d:itmax+1+d,jtmin-1-d:jtmin-1,ktmin-1-d:ktmax+1+d,1)=1;      %垂直于y轴
epsr(itmin-1-d:itmax+1+d,jtmin-1-d:jtmin-1,ktmin-1-d:ktmax+1+d,3)=1;
epsr(itmin-1-d:itmax+1+d,jtmax+1:jtmax+1+d,ktmin-1-d:ktmax+1+d,3)=1;
epsr(itmin-1-d:itmax+1+d,jtmax+1:jtmax+1+d,ktmin-1-d:ktmax+1+d,1)=1;
for j=jtmin-1-d:jtmin-1
    sig(itmin-1-d:itmax+1+d,j,ktmin-1-d:ktmax+1+d,2)=sigmax*(j-jtmin)^4/d^4;
    epsr(itmin-1-d:itmax+1+d,j,ktmin-1-d:ktmax+1+d,2)=1+(epsrmax-1)*(j-jtmin)^4/d^4;
end
for j=jtmax+1:jtmax+1+d
    sig(itmin-1-d:itmax+1+d,j,ktmin-1-d:ktmax+1+d,2)=sigmax*(j-jtmax)^4/d^4;
    epsr(itmin-1-d:itmax+1+d,j,ktmin-1-d:ktmax+1+d,2)=1+(epsrmax-1)*(j-jtmax)^4/d^4;
end
%%%%%%%%%%棱边处理
epsr(itmin-1-d,jtmin-1-d,ktmin-1-d:ktmax+1+d,3)=1;      %；平行于z轴
sig(itmin-1-d,jtmin-1-d,ktmin-1-d:ktmax+1+d,3)=0;
epsr(itmin-1-d,jtmax+1+d,ktmin-1-d:ktmax+1+d,3)=1;
sig(itmin-1-d,jtmax+1+d,ktmin-1-d:ktmax+1+d,3)=0;
epsr(itmax+1+d,jtmax+1+d,ktmin-1-d:ktmax+1+d,3)=1;
sig(itmax+1+d,jtmax+1+d,ktmin-1-d:ktmax+1+d,3)=0;
epsr(itmax+1+d,jtmin-1-d,ktmin-1-d:ktmax+1+d,3)=1;
sig(itmax+1+d,jtmin-1-d,ktmin-1-d:ktmax+1+d,3)=0;
                                                      %平行于x轴
epsr(itmin-1-d:itmax+1+d,jtmin-1-d,ktmax+1+d,1)=1;
sig(itmin-1-d:itmax+1+d,jtmin-1-d,ktmax+1+d,1)=0;
epsr(itmin-1-d:itmax+1+d,jtmax+1+d,ktmax+1+d,1)=1;
sig(itmin-1-d:itmax+1+d,jtmax+1+d,ktmax+1+d,1)=0;
epsr(itmin-1-d:itmax+1+d,jtmax+1+d,ktmin-1-d,1)=1;
sig(itmin-1-d:itmax+1+d,jtmax+1+d,ktmin-1-d,1)=0;
epsr(itmin-1-d:itmax+1+d,jtmin-1-d,ktmin-1-d,1)=1;
sig(itmin-1-d:itmax+1+d,jtmin-1-d,ktmin-1-d,1)=0;

epsr(itmin-1-d,jtmin-1-d:jtmax+1+d,ktmin-1-d,2)=1;         %平行于y轴
sig(itmin-1-d,jtmin-1-d:jtmax+1+d,ktmin-1-d,2)=0;
epsr(itmin-1-d,jtmin-1-d:jtmax+1+d,ktmax+1+d,2)=1;
sig(itmin-1-d,jtmin-1-d:jtmax+1+d,ktmax+1+d,2)=0;
epsr(itmax+1+d,jtmin-1-d:jtmax+1+d,ktmax+1+d,2)=1;
sig(itmax+1+d,jtmin-1-d:jtmax+1+d,ktmax+1+d,2)=0;
epsr(itmax+1+d,jtmin-1-d:jtmax+1+d,ktmin-1-d,2)=1;
sig(itmax+1+d,jtmin-1-d:jtmax+1+d,ktmin-1-d,2)=0;

D=zeros(ie,je,ke,3);                                                             %这是书中引入的几个量
B=zeros(ie,je,ke,3);
CA=zeros(ie,je,ke,3);                %CA CB 就是进行 D B的迭代式所用的系数
CB=zeros(ie,je,ke,3);                 %因为  公式的形式一致所以可以共用

DB=zeros(ie,je,ke,3);           %DB DC 还有下面的DA都是进行 e   和 h 分量时所使用
DC=zeros(ie,je,ke,3);           %的系数，但与书上不同  的是，分母中都没有相对介电、Q=itmax-itmin+2+d;              % 磁电常数。这两个可以在迭代时分别引入
for p=itmin-1-d:itmin-1          %CA CB DB DC 的1 2 3 分量对应的是公式上分子的x y z
    for t=1:3;                             %的分量。
        if t+2>3                      %在使用这些系数是注意要对应于迭代公式的对应x y z 的系数
            o=t-1;
        else                               %从这里开始进行系数的计算，往下的了两个循环也是。公
            o=t+2;                         %式都一致的，所以循环的内容也一样
        end
        CA(p,:,:,t)=(epsr(p,:,:,t)/dt-sig(p,:,:,t)/2/epsz)./(epsr(p,:,:,t)/dt+sig(p,:,:,t)/2/epsz);
        CB(p,:,:,t)=1/(epsr(p,:,:,t)/dt+sig(p,:,:,t)/2/epsz);
        DB(p,:,:,t)=(epsr(p,:,:,t)/dt+sig(p,:,:,t)/2/epsz)./(epsr(p,:,:,o)/dt+...
sig(p,:,:,o)/2/epsz);
        DC(p,:,:,t)=(epsr(p,:,:,t)/dt-sig(p,:,:,t)/2/epsz)./(epsr(p,:,:,o)/dt+sig(p,:,:,o)/2/epsz);

        CA(p+Q,:,:,t)=(epsr(p+Q,:,:,t)/dt-sig(p+Q,:,:,t)/2/epsz)./(epsr(p+Q,:,:,t)/dt+...
sig(p+Q,:,:,t)/2/epsz);
        CB(p+Q,:,:,t)=1/(epsr(p+Q,:,:,t)/dt+sig(p+Q,:,:,t)/2/epsz);
        DB(p+Q,:,:,t)=(epsr(p+Q,:,:,t)/dt+sig(p+Q,:,:,t)/2/epsz)./(epsr(p+Q,:,:,o)/dt+...
sig(p+Q,:,:,o)/2/epsz);
        DC(p+Q,:,:,t)=(epsr(p+Q,:,:,t)/dt-sig(p+Q,:,:,t)/2/epsz)./(epsr(p+Q,:,:,o)/dt+...
sig(p+Q,:,:,o)/2/epsz);
    end
end
Q=jtmax-jtmin+2+d;
for j=jtmin-1-d:jtmin-1
    for t=1:3
        if t+2>3
            o=t-1;
        else
            o=t+2;
        end
        CA(:,j,:,t)=(epsr(:,j,:,t)/dt-sig(:,j,:,t)/2/epsz)./(epsr(:,j,:,t)/dt+sig(:,j,:,t)/2/epsz);
        CB(:,j,:,t)=1/(epsr(:,j,:,t)/dt+sig(:,j,:,t)/2/epsz);
        DB(:,j,:,t)=(epsr(:,j,:,t)/dt+sig(:,j,:,t)/2/epsz)./(epsr(:,j,:,o)/dt+...
sig(:,j,:,o)/2/epsz);
        DC(:,j,:,t)=(epsr(:,j,:,t)/dt-sig(:,j,:,t)/2/epsz)./(epsr(:,j,:,o)/dt+sig(:,j,:,o)/2/epsz);
        CA(:,j+Q,:,t)=(epsr(:,j+Q,:,t)/dt-sig(:,j+Q,:,t)/2/epsz)./(epsr(:,j+Q,:,t)/dt+...
sig(:,j+Q,:,t)/2/epsz);
        CB(:,j+Q,:,t)=1/(epsr(:,j+Q,:,t)/dt+sig(:,j+Q,:,t)/2/epsz);
        DB(:,j+Q,:,t)=(epsr(:,j+Q,:,t)/dt+sig(:,j+Q,:,t)/2/epsz)./(epsr(:,j+Q,:,o)/dt+...
sig(:,j+Q,:,o)/2/epsz);
        DC(:,j+Q,:,t)=(epsr(:,j+Q,:,t)/dt-sig(:,j+Q,:,t)/2/epsz)./(epsr(:,j+Q,:,o)/dt+...
sig(:,j+Q,:,o)/2/epsz);
    end
end
Q=ktmax-ktmin+2+d;
for k=ktmin-1-d:ktmin-1
    for t=1:3
        if t+2>3
            o=t-1;
        else
            o=t+2;
        end
        CA(:,:,k,t)=(epsr(:,:,k,t)/dt-sig(:,:,k,t)/2/epsz)./(epsr(:,:,k,t)/dt+sig(:,:,k,t)/2/epsz);
        CB(:,:,k,t)=1./(epsr(:,:,k,t)/dt+sig(:,:,k,t)/2/epsz);
        DB(:,:,k,t)=(epsr(:,:,k,t)/dt+sig(:,:,k,t)/2/epsz)./(epsr(:,:,k,o)/dt+sig(:,:,k,o)/2/epsz);
        DC(:,:,k,t)=(epsr(:,:,k,t)/dt-sig(:,:,k,t)/2/epsz)./(epsr(:,:,k,o)/dt+sig(:,:,k,o)/2/epsz);
        CA(:,:,k+Q,t)=(epsr(:,:,k+Q,t)/dt-sig(:,:,k+Q,t)/2/epsz)./(epsr(:,:,k+Q,t)/dt+...
sig(:,:,k+Q,t)/2/epsz);
        CB(:,:,k+Q,t)=1./(epsr(:,:,k+Q,t)/dt+sig(:,:,k+Q,t)/2/epsz);
        DB(:,:,k+Q,t)=(epsr(:,:,k+Q,t)/dt+sig(:,:,k+Q,t)/2/epsz)./(epsr(:,:,k+Q,o)/dt+...
sig(:,:,k+Q,o)/2/epsz);
        DC(:,:,k+Q,t)=(epsr(:,:,k+Q,t)/dt-sig(:,:,k+Q,t)/2/epsz)./(epsr(:,:,k+Q,o)/dt+...
sig(:,:,k+Q,o)/2/epsz);
    end
end
DA=CA;
A=zeros(ie,je,ke,3);
M=zeros(ie,je,ke,3);
A=D;                        %用于存储D B的n时刻的值，以便在进行e 和 h分量迭代时使用
M=B;
for p=itmin-1-d:itmin-1                           %从这里开始UPML的循环计算、
    for j=jtmin-1-d:jtmax+1+d                   %往下的循环也是，只是区域不同而已
        for k=ktmin-1-d:ktmax+1+d
    D(p,j,k,1)=CA(p,j,k,2)*D(p,j,k,1)+CB(p,j,k,2)*(hz(p,j,k)-hz(p,j-1,k)-hy(p,j,k)+...
hy(p,j,k-1))/dx;
    D(p,j,k,2)=CA(p,j,k,3)*D(p,j,k,2)+CB(p,j,k,3)*(hx(p,j,k)-hx(p,j,k-1)-hz(p,j,k)+...
hy(p-1,j,k))/dx;
    D(p,j,k,3)=CA(p,j,k,1)*D(p,j,k,3)+CB(p,j,k,1)*(hy(p,j,k)-hy(p-1,j,k)-hx(p,j,k)+...
hx(p,j-1,k-1))/dx;

ex(p,j,k)=DA(p,j,k,3)*ex(p,j,k)+DB(p,j,k,1)*D(p,j,k,1)/epsz-DC(p,j,k,1)*A(p,j,k,1)/epsz;
    ey(p,j,k)=DA(p,j,k,1)*ey(p,j,k)+DB(p,j,k,2)*D(p,j,k,2)/epsz-DC(p,j,k,1)*A(p,j,k,2)/epsz;
    ez(p,j,k)=DA(p,j,k,2)*ez(p,j,k)+DB(p,j,k,3)*D(p,j,k,3)/epsz-DC(p,j,k,1)*A(p,j,k,3)/epsz;

    B(p,j,k,1)=CA(p,j,k,2)*B(p,j,k,1)-CB(p,j,k,2)*((ez(p,j+1,k))-ez(p,j,k)-ey(p,j,k+1)+...
ey(p,j,k))/dx;
    B(p,j,k,2)=CA(p,j,k,3)*B(p,j,k,3)-CB(p,j,k,3)*((ex(p,j,k+1))-ex(p,j,k)-ez(p+1,j,k)+...
ez(p,j,k))/dx;
    B(p,j,k,3)=CA(p,j,k,1)*B(p,j,k,2)-CB(p,j,k,1)*((ey(p+1,j,k))-ey(p,j,k)-ex(p,j+1,k)+...
ex(p,j,k))/dx;

    hx(p,j,k)=DA(p,j,k,3)*hx(p,j,k)+DB(p,j,k,1)*B(p,j,k,1)/muz-DC(p,j,k,1)*M(p,j,k,1)/muz;
    hy(p,j,k)=DA(p,j,k,1)*hy(p,j,k)+DB(p,j,k,2)*B(p,j,k,2)/muz-DC(p,j,k,2)*M(p,j,k,2)/muz;
    hz(p,j,k)=DA(p,j,k,2)*hz(p,j,k)+DB(p,j,k,3)*B(p,j,k,3)/muz-DC(p,j,k,3)*M(p,j,k,3)/muz;
        end
    end
end
for p=itmax+1:itmax+1+d
    for j=jtmin-1-d:jtmax+1+d
        for k=ktmin-1-d:ktmax+1+d
    D(p,j,k,1)=CA(p,j,k,2)*D(p,j,k,1)+CB(p,j,k,2)*(hz(p,j,k)-hz(p,j-1,k)-hy(p,j,k)+...
hy(p,j,k-1))/dx;
    D(p,j,k,2)=CA(p,j,k,3)*D(p,j,k,2)+CB(p,j,k,3)*(hx(p,j,k)-hx(p,j,k-1)-hz(p,j,k)+...
hy(p-1,j,k))/dx;
    D(p,j,k,3)=CA(p,j,k,1)*D(p,j,k,3)+CB(p,j,k,1)*(hy(p,j,k)-hy(p-1,j,k)-hx(p,j,k)+...
hx(p,j-1,k-1))/dx;

ex(p,j,k)=DA(p,j,k,3)*ex(p,j,k)+DB(p,j,k,1)*D(p,j,k,1)/epsz-DC(p,j,k,1)*A(p,j,k,1)/epsz;
    ey(p,j,k)=DA(p,j,k,1)*ey(p,j,k)+DB(p,j,k,2)*D(p,j,k,2)/epsz-DC(p,j,k,1)*A(p,j,k,2)/epsz;
    ez(p,j,k)=DA(p,j,k,2)*ez(p,j,k)+DB(p,j,k,3)*D(p,j,k,3)/epsz-DC(p,j,k,1)*A(p,j,k,3)/epsz;

    B(p,j,k,1)=CA(p,j,k,2)*B(p,j,k,1)-CB(p,j,k,2)*((ez(p,j+1,k))-ez(p,j,k)-ey(p,j,k+1)+...
ey(p,j,k))/dx;
    B(p,j,k,2)=CA(p,j,k,3)*B(p,j,k,3)-CB(p,j,k,3)*((ex(p,j,k+1))-ex(p,j,k)-ez(p+1,j,k)+...
ez(p,j,k))/dx;
    B(p,j,k,3)=CA(p,j,k,1)*B(p,j,k,2)-CB(p,j,k,1)*((ey(p+1,j,k))-ey(p,j,k)-ex(p,j+1,k)+...
ex(p,j,k))/dx;

    hx(p,j,k)=DA(p,j,k,3)*hx(p,j,k)+DB(p,j,k,1)*B(p,j,k,1)/muz-DC(p,j,k,1)*M(p,j,k,1)/muz;
    hy(p,j,k)=DA(p,j,k,1)*hy(p,j,k)+DB(p,j,k,2)*B(p,j,k,2)/muz-DC(p,j,k,2)*M(p,j,k,2)/muz;
    hz(p,j,k)=DA(p,j,k,2)*hz(p,j,k)+DB(p,j,k,3)*B(p,j,k,3)/muz-DC(p,j,k,3)*M(p,j,k,3)/muz;
        end
    end
end
for p=itmin-1-d:itmax+1+d
    for j=jtmin-1-d:jtmin-1
        for k=ktmin-1-d:ktmax+1+d
    D(p,j,k,1)=CA(p,j,k,2)*D(p,j,k,1)+CB(p,j,k,2)*(hz(p,j,k)-hz(p,j-1,k)-hy(p,j,k)+...
hy(p,j,k-1))/dx;
    D(p,j,k,2)=CA(p,j,k,3)*D(p,j,k,2)+CB(p,j,k,3)*(hx(p,j,k)-hx(p,j,k-1)-hz(p,j,k)+...
hy(p-1,j,k))/dx;
    D(p,j,k,3)=CA(p,j,k,1)*D(p,j,k,3)+CB(p,j,k,1)*(hy(p,j,k)-hy(p-1,j,k)-hx(p,j,k)+...
hx(p,j-1,k-1))/dx;

ex(p,j,k)=DA(p,j,k,3)*ex(p,j,k)+DB(p,j,k,1)*D(p,j,k,1)/epsz-DC(p,j,k,1)*A(p,j,k,1)/epsz;
    ey(p,j,k)=DA(p,j,k,1)*ey(p,j,k)+DB(p,j,k,2)*D(p,j,k,2)/epsz-DC(p,j,k,1)*A(p,j,k,2)/epsz;
    ez(p,j,k)=DA(p,j,k,2)*ez(p,j,k)+DB(p,j,k,3)*D(p,j,k,3)/epsz-DC(p,j,k,1)*A(p,j,k,3)/epsz;

    B(p,j,k,1)=CA(p,j,k,2)*B(p,j,k,1)-CB(p,j,k,2)*((ez(p,j+1,k))-ez(p,j,k)-ey(p,j,k+1)+ey(p,j,k))/dx;
    B(p,j,k,2)=CA(p,j,k,3)*B(p,j,k,3)-CB(p,j,k,3)*((ex(p,j,k+1))-ex(p,j,k)-ez(p+1,j,k)+ez(p,j,k))/dx;
    B(p,j,k,3)=CA(p,j,k,1)*B(p,j,k,2)-CB(p,j,k,1)*((ey(p+1,j,k))-ey(p,j,k)-ex(p,j+1,k)+ex(p,j,k))/dx;

    hx(p,j,k)=DA(p,j,k,3)*hx(p,j,k)+DB(p,j,k,1)*B(p,j,k,1)/muz-DC(p,j,k,1)*M(p,j,k,1)/muz;
    hy(p,j,k)=DA(p,j,k,1)*hy(p,j,k)+DB(p,j,k,2)*B(p,j,k,2)/muz-DC(p,j,k,2)*M(p,j,k,2)/muz;
    hz(p,j,k)=DA(p,j,k,2)*hz(p,j,k)+DB(p,j,k,3)*B(p,j,k,3)/muz-DC(p,j,k,3)*M(p,j,k,3)/muz;
        end
    end
end
for p=itmin-1-d:itmax+1+d
    for j=jtmax+1:jtmax+1+d
        for k=ktmin-1-d:ktmax+1+d
    D(p,j,k,1)=CA(p,j,k,2)*D(p,j,k,1)+CB(p,j,k,2)*(hz(p,j,k)-hz(p,j-1,k)-hy(p,j,k)+...
hy(p,j,k-1))/dx;
    D(p,j,k,2)=CA(p,j,k,3)*D(p,j,k,2)+CB(p,j,k,3)*(hx(p,j,k)-hx(p,j,k-1)-hz(p,j,k)+...
hy(p-1,j,k))/dx;
    D(p,j,k,3)=CA(p,j,k,1)*D(p,j,k,3)+CB(p,j,k,1)*(hy(p,j,k)-hy(p-1,j,k)-hx(p,j,k)+...
hx(p,j-1,k-1))/dx;

ex(p,j,k)=DA(p,j,k,3)*ex(p,j,k)+DB(p,j,k,1)*D(p,j,k,1)/epsz-DC(p,j,k,1)*A(p,j,k,1)/epsz;
    ey(p,j,k)=DA(p,j,k,1)*ey(p,j,k)+DB(p,j,k,2)*D(p,j,k,2)/epsz-DC(p,j,k,1)*A(p,j,k,2)/epsz;
    ez(p,j,k)=DA(p,j,k,2)*ez(p,j,k)+DB(p,j,k,3)*D(p,j,k,3)/epsz-DC(p,j,k,1)*A(p,j,k,3)/epsz;

    B(p,j,k,1)=CA(p,j,k,2)*B(p,j,k,1)-CB(p,j,k,2)*((ez(p,j+1,k))-ez(p,j,k)-ey(p,j,k+1)+...
ey(p,j,k))/dx;
    B(p,j,k,2)=CA(p,j,k,3)*B(p,j,k,3)-CB(p,j,k,3)*((ex(p,j,k+1))-ex(p,j,k)-ez(p+1,j,k)+...
ez(p,j,k))/dx;
    B(p,j,k,3)=CA(p,j,k,1)*B(p,j,k,2)-CB(p,j,k,1)*((ey(p+1,j,k))-ey(p,j,k)-ex(p,j+1,k)+...
ex(p,j,k))/dx;

    hx(p,j,k)=DA(p,j,k,3)*hx(p,j,k)+DB(p,j,k,1)*B(p,j,k,1)/muz-DC(p,j,k,1)*M(p,j,k,1)/muz;
    hy(p,j,k)=DA(p,j,k,1)*hy(p,j,k)+DB(p,j,k,2)*B(p,j,k,2)/muz-DC(p,j,k,2)*M(p,j,k,2)/muz;
    hz(p,j,k)=DA(p,j,k,2)*hz(p,j,k)+DB(p,j,k,3)*B(p,j,k,3)/muz-DC(p,j,k,3)*M(p,j,k,3)/muz;
        end
    end
end
for p=itmin-1-d:itmax+1+d
    for j=jtmin-1-d:jtmax+1+d
        for k=ktmin-1-d:ktmin-1
    D(p,j,k,1)=CA(p,j,k,2)*D(p,j,k,1)+CB(p,j,k,2)*(hz(p,j,k)-hz(p,j-1,k)-hy(p,j,k)+...
hy(p,j,k-1))/dx;
    D(p,j,k,2)=CA(p,j,k,3)*D(p,j,k,2)+CB(p,j,k,3)*(hx(p,j,k)-hx(p,j,k-1)-hz(p,j,k)+...
hy(p-1,j,k))/dx;
    D(p,j,k,3)=CA(p,j,k,1)*D(p,j,k,3)+CB(p,j,k,1)*(hy(p,j,k)-hy(p-1,j,k)-hx(p,j,k)+...
hx(p,j-1,k-1))/dx;

ex(p,j,k)=DA(p,j,k,3)*ex(p,j,k)+DB(p,j,k,1)*D(p,j,k,1)/epsz-DC(p,j,k,1)*A(p,j,k,1)/epsz;
    ey(p,j,k)=DA(p,j,k,1)*ey(p,j,k)+DB(p,j,k,2)*D(p,j,k,2)/epsz-DC(p,j,k,1)*A(p,j,k,2)/epsz;
    ez(p,j,k)=DA(p,j,k,2)*ez(p,j,k)+DB(p,j,k,3)*D(p,j,k,3)/epsz-DC(p,j,k,1)*A(p,j,k,3)/epsz;

    B(p,j,k,1)=CA(p,j,k,2)*B(p,j,k,1)-CB(p,j,k,2)*((ez(p,j+1,k))-ez(p,j,k)-ey(p,j,k+1)+...
ey(p,j,k))/dx;
    B(p,j,k,2)=CA( ..

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

描述:我写的程序
附件: