[已解决]三维pml matlab 程序，有问题，大侠帮忙啊，急！

经过进一步检验，这是修改之后的程序，不过仍存在问题：1、激励源的加法，我的激励源修改后分别加到两个分量上，但距离源点几个网格后仍产生毛刺，说明源的加入还是有问题；2、我的程序是把自由空间和pml层一起做了运算和迭代，其实边角问题已经不存在，不需要考虑了（二维已经验证了），这样编程的好处是缩短程序，全域统一编程。 i`!>zl+D  
现在的程序仍有毛刺，请高手帮我改正，万分感激！！！谢谢！！ PD-*rG `  
clear Lv>OBHD  
c=3e8;   Dt'bbX'edw  
muz=4*pi*1e-7;   WA#y&  
epsz=1/(c*c*muz);        $kxu-  
f=1e+9;               j$P`/-N  
lambda=c/f;             $,nidK!"  
omega=2*pi*f; Ru$%gh>v  
ds=lambda/12;             K!v\r
"N  
dt=ds/(2*c); xN!In-v[j;  
muz=4*pi*1e-7;   #d$lN}8  
R0=1e-5;   r>6FJ:Tx  
K#R|GEwr  
%************************************************ <jh=W9.N_  
ie=50;                                                                    %x向空间步             CxrsP.  
je=50;                                                                    %y向空间步      
3s#/d,+  
ke=50;                                                                    %z向空间步  *1*i5c  
ib=ie+1; E=ObfN"ge  
jb=je+1; mJsU7bD`  
kb=ke+1;
X#ud_+6x  
is=25;                                                                     %激励源x向位置 nd{k D>a  
js=25;                                                                     %激励源y向位置 hwSxdT6  
ks=25;                                                                     %激励源z向位置 43*;"w=  
iebc=8;                                                                    %左、右PML层层数 2|${2u`$&y  
jebc=8;                                                                    %上、下PML层层数   Tzfk_h3hE  
kebc=8;                                                                    %前、后PML层层数 ]<u%jTQREd  
ibbc=iebc+1; Tz7|OV_W$  
jbbc=jebc+1; i4)]lWnd  
kbbc=kebc+1; _TbvQY  
nmax=100;                                                                  %时间步 hcR^?  
dlta=6;                                                                    %观测点距激励源步数 LKg9{0Y:  
)g
k tI!  
%Lp#2?*  
%************************************************ J$S*QCo  
Exy=zeros(ie,jb,kb); p\tA&>3-  
Exz=zeros(ie,jb,kb); )4n]n:FjN  
Eyx=zeros(ib,je,kb); Q'aVdJN,  
Eyz=zeros(ib,je,kb); C-H6l6,  
Ezx=zeros(ib,jb,ke); tQ)l4Y 8  
Ezy=zeros(ib,jb,ke); l7(p~+o?h>  
Hxy=zeros(ib,je,ke); [=>[2Ty  
Hxz=zeros(ib,je,ke); 4H`B]Zt7  
Hyx=zeros(ie,jb,ke); XLlJ|xhY-K  
Hyz=zeros(ie,jb,ke); ;33SUgX  
Hzx=zeros(ie,je,kb); ozl>Au  
Hzy=zeros(ie,je,kb); -#r=  
dlMjy$/
T  
H5be5  
fid1=fopen('m1.txt','wt'); sc z8`%  
fid2=fopen('m2.txt','wt'); CUj$ <ay=  
%************************************************ GYV%RD#  
sigxmax=-log(R0)/(2*iebc*ds)*3*epsz*c;                                     %x向最大电导率 ',I$`h  
sigymax=-log(R0)/(2*jebc*ds)*3*epsz*c;                                     %y向最大电导率 xo{z4W  
sigzmax=-log(R0)/(2*kebc*ds)*3*epsz*c; VI)hA ^S  
fBKN?]BdN  
pAyUQe;X#  
for i=1:ib                                                                 %x向网格节点处电导率分布 8L*#zaSAf  
                if      i<=iebc DKG;up0  
                    sigx(i)=sigxmax*((iebc+1-i)/iebc)^2;                      %左侧PML层 /`7G7pQ+  
                elseif i>=ib-iebc+1                                           %右侧PML层 `*cJc6  
                    sigx(i)=sigxmax*((i-ib+iebc)/iebc)^2; )ZpMB  
                else 5[P^O6'  
                    sigx(i)=0;                                              %中间自由空间 BfQ#5  
                end     _OB^ywHn.  
      cax(i)=(1-(dt*sigx(i))/(2*epsz))/(1+(dt*sigx(i))/(2*epsz)); v1 f^gde  
      cbx(i)=(dt/epsz/ds)/(1+(dt*sigx(i))/(2*epsz)); U:8cz=#  
       Iv?1XI=  
end   {T3wOi
  
for i=1:ie                                                                 %x方向网格边中点处电导率分布 afjEN y1  
                if      i<=iebc @Suww@<  
                    sigxh(i)=sigxmax*((iebc+1-i-0.5)/iebc)^2;                 %左侧PML层 6ciA|J'MR  
                elseif i>=ie-iebc+1                                           %右侧PML层 jK2gc^"t  
                    sigxh(i)=sigxmax*((i-ie+iebc-0.5)/iebc)^2; fF.+{-.  
                else Rd|^C$6  
                    sigxh(i)=0;                                              %中间自由空间 y4Nam87;/?  
                end                 ^%qQ)>I=j  
      roux(i)=muz*sigxh(i)/epsz;                                              %磁导率       3Q_)Xs r`  
      cpx(i)=(1-(dt*roux(i))/(2*muz))/(1+(dt*roux(i))/(2*muz)); M+7jJ?n  
      cqx(i)=(dt/muz/ds)/(1+(dt*roux(i))/(2*muz)); u89Q2\z~"M  
end   Hh%|}*f_,  
@`
HW0Y_:  
x/%/MFK)>8  
|C7=$DgwY  
for j=1:jb                                                                 %y方向网格节点处电导率分布 S0;s 7X#c  
                if      j<=jebc ?Xqkf>  
                    sigy(j)=sigymax*((jebc+1-j)/jebc)^2;                      %下侧PML层 Qi|k,1A0  
                elseif j>=jb-jebc+1                                           %上侧PML层 f5'vjWJ30  
                    sigy(j)=sigymax*((j-jb+jebc)/jebc)^2; z/i&Lpr:  
                else 5w</Ga  
                    sigy(j)=0;                                              %中间自由空间 kuv+TN  
                end               cZAf?,>u  
      cay(j)=(1-(dt*sigy(j))/(2*epsz))/(1+(dt*sigy(j))/(2*epsz)); 
BuS[(  
      cby(j)=(dt/epsz/ds)/(1+(dt*sigy(j))/(2*epsz)); +aOX{1w  
end   P} Y .  
for j=1:je                                                                 %y方向网格边中点处电导率分布                 ]F"@+_E  
                if      j<=jebc *++}ll6  
                    sigyh(j)=sigymax*((jebc+1-j-0.5)/jebc)^2;                 %下侧PML层 LKst QP!I  
                elseif j>=je-jebc+1                                           %上侧PML层 J n'SGR  
                    sigyh(j)=sigymax*((j-je+jebc-0.5)/jebc)^2; \Lm`jU(:l  
                else 8aW<lu  
                    sigyh(j)=0;                                              %中间自由空间 PsN_c[+  
                end O*c<m,  
      rouy(j)=muz*sigyh(j)/epsz;                                              %磁导率       )yS8(F0  
      cpy(j)=(1-(dt*rouy(j))/(2*muz))/(1+(dt*rouy(j))/(2*muz)); l$C Y gm  
      cqy(j)=(dt/muz/ds)/(1+(dt*rouy(j))/(2*muz));         !~Kg_*IT  
end   ~P"o
_b6,k  
;V84Dy#b  
9M@,BXOt  
"
nU] 2  
for k=1:kb                                                                 %z方向网格节点处电导率分布 n{"a0O  
                if      k<=kebc w+hpi
5OH  
                    sigz(k)=sigzmax*((kebc+1-k)/kebc)^2;                      %前侧PML层 c=p!2jJ1K~  
                elseif k>=kb-kebc+1                                           %后侧PML层 @Z Dd(xB&  
                    sigz(k)=sigzmax*((k-kb+kebc)/kebc)^2; ]i8t  
                else TmJXkR.5  
                    sigz(k)=0;                                              %中间自由空间 mH )i  
                end Z5[g[Q  
      caz(k)=(1-(dt*sigz(k))/(2*epsz))/(1+(dt*sigz(k))/(2*epsz)); {
}BAQ9|q  
      cbz(k)=(dt/epsz/ds)/(1+(dt*sigz(k))/(2*epsz));       @R;&PR#5  
end  
U=[isi+7  
for k=1:ke                                                                 %z方向网格边中点处电导率分布                 BxB B](  
                if      k<=kebc JG{`tTu  
                    sigzh(k)=sigzmax*((kebc+1-k-0.5)/kebc)^2;                 %前侧PML层 a&B@F]+  
                elseif k>=ke-kebc+1                                           %后侧PML层 t_\&LMD  
                    sigzh(k)=sigzmax*((k-ke+kebc-0.5)/kebc)^2; >*$Xbj*  
                else p,7?rI\N  
                    sigzh(k)=0;                                              %中间自由空间 }w{E<C(M  
                end .L.9e#?3  
      rouz(k)=muz*sigzh(k)/epsz;                                              %磁导率     lF\2a&YRbn  
      cpz(k)=(1-(dt*rouz(k))/(2*muz))/(1+(dt*rouz(k))/(2*muz)); /b410NP5  
      cqz(k)=(dt/muz/ds)/(1+(dt*rouz(k))/(2*muz));         -f"{%<Q  
end   1tlqw  
?3<Y/Vg%c  
7hT@,|(j  
% Ex=Exy+Exz; ;dTxQ_:  
% Ey=Eyx+Eyz; mcxD#+H 3  
% Ez=Ezy+Ezx; gq~"Z[T  
=0SJf 3  
%************************************************ ]^<\a=U  
tic >".@;  
L),bPfz  
schedule_box=waitbar(0,'Please wait...'); % display the schedule box Pxgul7  
                                                                       Op hD_^  
for n=1:nmax     %I@vMs^  
ul!q)cPb{  
           for i=1:ib #{J+BWP\o  
               for j=1:je K+Pa b ?  
                 for k=1:ke )-25?B  
                    if(i==is&&j==js&&k==ks) uYijzHQyD  
                        Hxy(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 ~w}=Oby'y  
                    else mv(/M t  
                        Hxy(i,j,k)=cpy(j).*Hxy(i,j,k)+cqy(j).*(Ezx(i,j,k)+Ezy(i,j,k)-Ezx(i,j+1,k)-Ezy(i,j+1,k)); SZvC4lOn#  
                    end kkXe=f%  
                                         =
]1g*~%  
                 end 6xu%M&ht  
               end 7t+H94KG7  
           end j_S3<wEJ  
             3A\Z]L  
           @
@=,bO  
            for i=1:ib <]*Jhnx/  
               for j=1:je 'h[7AZ&)#  
                 for k=1:ke ZHCrKp  
                   if(i==is&&j==js&&k==ks) CiTWjE?|7  
                       Hxz(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 x!u6LDq0  
                   else )M8@|~~  
                       Hxz(i,j,k)=cpz(k).*Hxz(i,j,k)+cqz(k).*(Eyx(i,j,k+1)+Eyz(i,j,k+1)-Eyx(i,j,k)-Eyz(i,j,k)); {~#d_!(  
                   end ,C|aiSh0-  
%                    Hxz(i,j,k)=cpz(k).*Hxz(i,j,k)+cqz(k).*(Eyx(i,j,k+1)+Eyz(i,j,k+1)-Eyx(i,j,k)-Eyz(i,j,k)); f4vdJ5pV  
               end .`oKd@I*"  
               end TV}=$\D  
          end xM/WS':V  
%           Hx=Hxy+Hxz; S`vw<u4t  
          for i=1:ie qJR8fQ
  
              for j=1:jb U@:l~xJ  
                 for k=1:ke 8yM8O #S  
                   Hyz(i,j,k)=cpz(k).*Hyz(i,j,k)+cqz(k).*(Exy(i,j,k)+Exz(i,j,k)-Exy(i,j,k+1)-Exz(i,j,k+1)); bG+Gg*0p  
                 end r>
G||/Z  
               end BL%3[JQ  
           end zR?1iV.]  
          for i=1:ie                                                 qJw\<7m  
               for j=1:jb %cASk>
^i  
                 for k=1:ke #r]GnC,  
                   Hyx(i,j,k)=cpx(i).*Hyx(i,j,k)+cqx(i).*(Ezx(i+1,j,k)+Ezy(i+1,j,k)-Ezx(i,j,k)-Ezy(i,j,k)); s:I 8~Cc  
               end GE\({V.W  
               end wkT4R\H>  
          end 5'_:>0}  
%           Hy=Hyx+Hyz; <ILi38%Y  
          for i=1:ie                                                   muO;g&  
               for j=1:je K
]&GSro  
                 for k=1:kb <)oxs]<  
                   Hzx(i,j,k)=cpx(i).*Hzx(i,j,k)+cqx(i).*(Eyx(i,j,k)+Eyz(i,j,k)-Eyx(i+1,j,k)-Eyz(i+1,j,k)); o_r{cnu  
               end FV%|*JW[;N  
               end 4+4&}8FH  
           end J$-1odL0Z  
          for i=1:ie                                                   m*v@L4t(1  
               for j=1:je 2SKtdiY  
                 for k=1:kb 3w! NTvp  
                   Hzy(i,j,k)=cpy(j).*Hzy(i,j,k)+cqy(j).*(Exy(i,j+1,k)+Exz(i,j+1,k)-Exy(i,j,k)-Exz(i,j,k)); 2(R{3E4.  
               end >uE<-klv  
               end [ *
!0DW`  
          end B&}lYo  
% Hz=Hzx+Hzy; Zm#,Ike?#  
|^#Z!Hp_Y  
          for i=1:ie                                                   8_3WCbe/  
               for j=2:je  % PEC at j=1 & j=jb NSQ)lSW,;  
                 for k=2:ke  % PEC at k=1 & k=kb  后面相同 z0T6a15f!P  
                   Exy(i,j,k)=cay(j).*Exy(i,j,k)+cby(j).*(Hzx(i,j,k)+Hzy(i,j,k)-Hzx(i,j-1,k)-Hzy(i,j-1,k)); +\J+?jOC4S  
                 end ryTtGx%a  
               end jjg&C9w T  
          end ,Uy~O(Ft  
          for i=1:ie                                               \]Rmq_O  
               for j=2:je .GFKy  
                 for k=2:ke 4hUUQ;xj  
                   Exz(i,j,k)=caz(k).*Exz(i,j,k)+cbz(k).*(Hyx(i,j,k-1)-Hyx(i,j,k)+Hyz(i,j,k-1)-Hyz(i,j,k)); A \Z_br  
                 end fg GTm:   
               end > Qh#pn*  
          end UO-,A j*wW  
iF1zLI<A  
            for i=2:ie                                                     IOl0=+p  
               for j=1:je A$TFa:O|  
                 for k=2:ke w^S]HzMd  
                   Eyz(i,j,k)=caz(k).*Eyz(i,j,k)+cbz(k).*(Hxy(i,j,k)-Hxy(i,j,k-1)+Hxz(i,j,k)-Hxz(i,j,k-1)); 8vUq8[[  
                 end N Q}5'  
               end @>[3[;  
            end s_6Iz^]I  
             ) 3I|6iS  
           .. M,SIs 3  
 FOqD  

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

你说电磁场迭代有大问题，到底是什么大问题，是电磁波传播不稳定？很快就变成无穷大了么？还是计算结果和预期不对？

是除了我的激励那有波形，其他点一律没有产生迭代

第一是PML煤质中损耗电导率的计算在不同的场量(相差半网格)上是不一样的。

分离场的PML，是最繁琐的、代码量最多的PML方法，大体一看代码的长度，就知道算法编写有问题。

分离的场步进计算的系数看着有点别扭啊

如果都没有迭代的话，自己调试一下吧。 =c{/ Z  
Hxy(is,js,ks)是源，那么Eyz和Ezy应该是可以得到值的，你把这几个分量都实时显示出来，看看问题到底在哪里吧。 V~rF`1+5N  
9~8UG (  
或者你给我个email，我发一个2D+PML的源码给你参考一下。

碟代/碟代/碟代/碟代/

版主,最近一直没上网,不好意思,这个程序我还在研究,下个帖子我发最新得程序,但仍未有最好效果,也就是说这个问题还没有解决掉.

你是直接做的3D的吗？