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

经过进一步检验，这是修改之后的程序，不过仍存在问题：1、激励源的加法，我的激励源修改后分别加到两个分量上，但距离源点几个网格后仍产生毛刺，说明源的加入还是有问题；2、我的程序是把自由空间和pml层一起做了运算和迭代，其实边角问题已经不存在，不需要考虑了（二维已经验证了），这样编程的好处是缩短程序，全域统一编程。 4(aesZ8h  
现在的程序仍有毛刺，请高手帮我改正，万分感激！！！谢谢！！ \aZ(@eF@@Q  
clear uAj

GR  
c=3e8;   ^Nsl5  
muz=4*pi*1e-7;   u=E?N:I~F  
epsz=1/(c*c*muz);       (v4  
f=1e+9;               TLSy+x_gX  
lambda=c/f;             5&X  
omega=2*pi*f; uCu,'F,6Y  
ds=lambda/12;             =(n'#mV  
dt=ds/(2*c); -/gS s<"  
muz=4*pi*1e-7;   )%y~{j+M  
R0=1e-5;   y9*H  
9uS7G*  
%************************************************ MW`a>'0t?  
ie=50;                                                                    %x向空间步             6ZG)`u".("  
je=50;                                                                    %y向空间步       u4<r$[]V  
ke=50;                                                                    %z向空间步 qz0v1057#  
ib=ie+1; _).'SU)>  
jb=je+1; R
NA03  
kb=ke+1; L|q<Bpz  
is=25;                                                                     %激励源x向位置 #Y5I_:k  
js=25;                                                                     %激励源y向位置 :):Y6)giBD  
ks=25;                                                                     %激励源z向位置 gw*d"~A  
iebc=8;                                                                    %左、右PML层层数 +1)C&:  
jebc=8;                                                                    %上、下PML层层数   <;6])  
kebc=8;                                                                    %前、后PML层层数 l#~FeD  
ibbc=iebc+1; GVl u4  
jbbc=jebc+1; Vw{Ys6q  
kbbc=kebc+1; -8tA~;p
 
nmax=100;                                                                  %时间步 QhGg^h%6  
dlta=6;                                                                    %观测点距激励源步数 *}_/:\v  
j&Hn`G  
Z@[,"{Sn  
%************************************************ H~ZSw7!M8  
Exy=zeros(ie,jb,kb); `>sqP aD  
Exz=zeros(ie,jb,kb); h+74W0 $  
Eyx=zeros(ib,je,kb); YjX=@  
Eyz=zeros(ib,je,kb); 5dgBSL$A}]  
Ezx=zeros(ib,jb,ke); w@: ]]R  
Ezy=zeros(ib,jb,ke); hL`zV  
Hxy=zeros(ib,je,ke); W1@;94Sb~  
Hxz=zeros(ib,je,ke); X#3<hN*v  
Hyx=zeros(ie,jb,ke); H*\[:tPa  
Hyz=zeros(ie,jb,ke); Pj.~|5gnf  
Hzx=zeros(ie,je,kb); y:A0!75  
Hzy=zeros(ie,je,kb); yX-xVvlv@  
R{
[Q+y'E  
OpL6Y+<  
fid1=fopen('m1.txt','wt'); =wj~6:Bf  
fid2=fopen('m2.txt','wt'); @qC:% |>  
%************************************************ P+b^;+\1s  
sigxmax=-log(R0)/(2*iebc*ds)*3*epsz*c;                                     %x向最大电导率 b}4/4Z.  
sigymax=-log(R0)/(2*jebc*ds)*3*epsz*c;                                     %y向最大电导率 {cv;S2  
sigzmax=-log(R0)/(2*kebc*ds)*3*epsz*c; k_a'a)`$6  
xF/DYXC{8  
aqM_t  
for i=1:ib                                                                 %x向网格节点处电导率分布 
*n*y!z  
                if      i<=iebc ya_'
Oz!C  
                    sigx(i)=sigxmax*((iebc+1-i)/iebc)^2;                      %左侧PML层 gPwp [  
                elseif i>=ib-iebc+1                                           %右侧PML层 x>J3tp$2  
                    sigx(i)=sigxmax*((i-ib+iebc)/iebc)^2; t3GK{X  
                else Hxl,U>za#  
                    sigx(i)=0;                                              %中间自由空间 +E [bLz^  
                end     DrEtnt  
      cax(i)=(1-(dt*sigx(i))/(2*epsz))/(1+(dt*sigx(i))/(2*epsz)); @}?D<O8#"#  
      cbx(i)=(dt/epsz/ds)/(1+(dt*sigx(i))/(2*epsz)); iCK$ o_`?  
       zOE6;c81  
end   xI<dBg|]+  
for i=1:ie                                                                 %x方向网格边中点处电导率分布 s =5H.q%PV  
                if      i<=iebc `e9uSF:9C  
                    sigxh(i)=sigxmax*((iebc+1-i-0.5)/iebc)^2;                 %左侧PML层 ^97ZH)Ww  
                elseif i>=ie-iebc+1                                           %右侧PML层 \Zv =?\  
                    sigxh(i)=sigxmax*((i-ie+iebc-0.5)/iebc)^2; Az-!LAu9 R  
                else - X_w&  
                    sigxh(i)=0;                                              %中间自由空间 _Y|kX2l S@  
                end                 [;=ky<K0E  
      roux(i)=muz*sigxh(i)/epsz;                                              %磁导率       2wBU@T1  
      cpx(i)=(1-(dt*roux(i))/(2*muz))/(1+(dt*roux(i))/(2*muz)); 'sF563kE  
      cqx(i)=(dt/muz/ds)/(1+(dt*roux(i))/(2*muz)); 0l6iv[qu5w  
end   YxtkI:C?
  
,>UmKrYo  
{T.Vu]L80  
h(8;7}K  
for j=1:jb                                                                 %y方向网格节点处电导率分布 `az`?`i7  
                if      j<=jebc `_'Dj>  
                    sigy(j)=sigymax*((jebc+1-j)/jebc)^2;                      %下侧PML层 Uu X"AFy~\  
                elseif j>=jb-jebc+1                                           %上侧PML层 ^d9
raYE`'  
                    sigy(j)=sigymax*((j-jb+jebc)/jebc)^2; 2SJh
6U  
                else :'dc=C  
                    sigy(j)=0;                                              %中间自由空间 .V3Dql@z"  
                end               8e@JvAaa$  
      cay(j)=(1-(dt*sigy(j))/(2*epsz))/(1+(dt*sigy(j))/(2*epsz)); B<XPu=|  
      cby(j)=(dt/epsz/ds)/(1+(dt*sigy(j))/(2*epsz)); / 3k\kkv!  
end   YgiGI <U  
for j=1:je                                                                 %y方向网格边中点处电导率分布                 Mak9qaWqF>  
                if      j<=jebc 4LjSDgA  
                    sigyh(j)=sigymax*((jebc+1-j-0.5)/jebc)^2;                 %下侧PML层 E`4=C@NN+,  
                elseif j>=je-jebc+1                                           %上侧PML层 9AYe,R  
                    sigyh(j)=sigymax*((j-je+jebc-0.5)/jebc)^2; "J5Pwvs-  
                else Mdrv/x{  
                    sigyh(j)=0;                                              %中间自由空间 fE*I
+pe  
                end ?03Zy3/  
      rouy(j)=muz*sigyh(j)/epsz;                                              %磁导率       lDU:EJ&DHE  
      cpy(j)=(1-(dt*rouy(j))/(2*muz))/(1+(dt*rouy(j))/(2*muz)); V 3]p3  
      cqy(j)=(dt/muz/ds)/(1+(dt*rouy(j))/(2*muz));         1Kh?JH  
end   mG~y8nUtp  
B%@!\D#  
a1`cI5
n  
erP>P  
for k=1:kb                                                                 %z方向网格节点处电导率分布 t4k'9Y:\Q  
                if      k<=kebc arLl8G[  
                    sigz(k)=sigzmax*((kebc+1-k)/kebc)^2;                      %前侧PML层 Hm*vKFhz  
                elseif k>=kb-kebc+1                                           %后侧PML层 ptv4v[gQ  
                    sigz(k)=sigzmax*((k-kb+kebc)/kebc)^2; vEe  
                else <e&QTyb  
                    sigz(k)=0;                                              %中间自由空间 IJc#)J.2A  
                end 1:!rw,Jzl`  
      caz(k)=(1-(dt*sigz(k))/(2*epsz))/(1+(dt*sigz(k))/(2*epsz)); _<a)\UR  
      cbz(k)=(dt/epsz/ds)/(1+(dt*sigz(k))/(2*epsz));       (YJAT  
end   Gr@{p"./z  
for k=1:ke                                                                 %z方向网格边中点处电导率分布                 >1U@NK)HfY  
                if      k<=kebc 'h^DI`  
                    sigzh(k)=sigzmax*((kebc+1-k-0.5)/kebc)^2;                 %前侧PML层 =+h!JgY/L  
                elseif k>=ke-kebc+1                                           %后侧PML层 B^(rUR  
                    sigzh(k)=sigzmax*((k-ke+kebc-0.5)/kebc)^2; G`#
gV"PlC  
                else N&ql(#r  
                    sigzh(k)=0;                                              %中间自由空间 7=.VqC^  
                end e-cb?.WU?  
      rouz(k)=muz*sigzh(k)/epsz;                                              %磁导率     tN3 {7'\7  
      cpz(k)=(1-(dt*rouz(k))/(2*muz))/(1+(dt*rouz(k))/(2*muz)); /mwr1GU  
      cqz(k)=(dt/muz/ds)/(1+(dt*rouz(k))/(2*muz));         'B5J.Xe:  
end   ba9<(0`  
-fx88  
')<FLCFwT  
% Ex=Exy+Exz; \ui^ d  
% Ey=Eyx+Eyz; 9BD|uU;0  
% Ez=Ezy+Ezx; YaZt+WA  
DsW`V~T  
%************************************************ 'HWgvmw(  
tic A]?O&m|  
g**%J Xo  
schedule_box=waitbar(0,'Please wait...'); % display the schedule box 7 3z Y^x  
                                                                       0bxvM  
for n=1:nmax     j'*p  
M y"!j,Up  
           for i=1:ib ZU.)K
>'  
               for j=1:je -mHhB(Td'  
                 for k=1:ke WQLL[{mhS  
                    if(i==is&&j==js&&k==ks) '}`hY1v  
                        Hxy(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 @ KPv&UB  
                    else Mto~ /  
                        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)); :5K~/=6x  
                    end n{Qh8"  
                                         K}=8:BaUL  
                 end sHTePEJ_h  
               end w52HN;Jm  
           end q3-;}+  
             D^s0EW-E  
           }Q=se[((  
            for i=1:ib <oaBh)=7  
               for j=1:je sR>;h /  
                 for k=1:ke Px7g\[]  
                   if(i==is&&j==js&&k==ks) Z^`=!n-V  
                       Hxz(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 ms\/=96F  
                   else +Q!xEfpO;  
                       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)); ;k&k#>L!K  
                   end lcT+$4zk.  
%                    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 ,ZNh1xt  
               end i)=89?8  
               end FV6he[,  
          end y]pN=<*h5  
%           Hx=Hxy+Hxz; KvJP(!{  
          for i=1:ie :IT U0%;!+  
              for j=1:jb `{|}LFS>  
                 for k=1:ke `?o1cf A  
                   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)); e=]oh$]  
                 end /-K dCp~  
               end A[ECa{v  
           end x\bRj>%(  
          for i=1:ie                                                 ck
jVa\  
               for j=1:jb 3Ra\2(bR  
                 for k=1:ke Ddl% V7  
                   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)); W3V{Xk|  
               end D@(M+u9/%  
               end ({_:^$E\  
          end )Kk(P/s  
%           Hy=Hyx+Hyz; d!QD vO  
          for i=1:ie                                                   Z:Y.":[ Qi  
               for j=1:je V0'p1J tD  
                 for k=1:kb -KwL9J4u  
                   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)); H=o-ScA  
               end [OoH5dD  
               end KYRm Ui#  
           end .xz,pn}  
          for i=1:ie                                                   E2h;hr;W  
               for j=1:je p-yOiG8b}  
                 for k=1:kb u};]LX\E  
                   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)); P$6f+{  
               end &^3~=$   
               end [%iUg\'7d  
          end &X]=Qpl  
% Hz=Hzx+Hzy; `/wq3+?  
@Uj_+c q  
          for i=1:ie                                                   X|1_0  
               for j=2:je  % PEC at j=1 & j=jb qtwT#z;Y  
                 for k=2:ke  % PEC at k=1 & k=kb  后面相同 Vi?~0.Z
%  
                   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)); &v auLp  
                 end rGn5Q
V  
               end f0mH|tI`  
          end %I|+_ z&x  
          for i=1:ie                                               O^R
:_vb3I  
               for j=2:je lGnql1(  
                 for k=2:ke ]~ #+b>  
                   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)); 3dxnh,]&@  
                 end IGK_1@tq  
               end K*_{Rs0P  
          end  ny  
_5U%'\5s  
            for i=2:ie                                                     ,+NE:_  
               for j=1:je -0KbdHIKb'  
                 for k=2:ke ycl>git]  
                   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)); (@WDvgi(  
                 end X6Nm!od'  
               end 0(hv#C4  
            end S QM(8*:X  
             H81.p  
           .. R%
)2(\  
iA%' ;V  

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

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

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

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

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

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

如果都没有迭代的话，自己调试一下吧。 7m1KR#j  
Hxy(is,js,ks)是源，那么Eyz和Ezy应该是可以得到值的，你把这几个分量都实时显示出来，看看问题到底在哪里吧。 !/I0i8T  
zAScRg$:?  
或者你给我个email，我发一个2D+PML的源码给你参考一下。

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

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

你是直接做的3D的吗？