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

经过进一步检验，这是修改之后的程序，不过仍存在问题：1、激励源的加法，我的激励源修改后分别加到两个分量上，但距离源点几个网格后仍产生毛刺，说明源的加入还是有问题；2、我的程序是把自由空间和pml层一起做了运算和迭代，其实边角问题已经不存在，不需要考虑了（二维已经验证了），这样编程的好处是缩短程序，全域统一编程。 {*|y
U"  
现在的程序仍有毛刺，请高手帮我改正，万分感激！！！谢谢！！ p?}Rolk7  
clear /`1zkBj<&  
c=3e8;   GL/\uq  
muz=4*pi*1e-7;   3oSQe"  
epsz=1/(c*c*muz);       PY^Yx$t9  
f=1e+9;               T|E;U  
lambda=c/f;             @1>83-p"X  
omega=2*pi*f;
$Ec;w~e  
ds=lambda/12;             Kg.E~  
dt=ds/(2*c); B82A:t)  
muz=4*pi*1e-7;   :g,
rl\S7  
R0=1e-5;   ,^
+3AT  
D^A_0@  
%************************************************ =Xp3UNXg  
ie=50;                                                                    %x向空间步             KAe) X_R7  
je=50;                                                                    %y向空间步       tHGK<rb  
ke=50;                                                                    %z向空间步 5'o.v^l  
ib=ie+1; 8^^al!0K~  
jb=je+1; iw#luHc
J  
kb=ke+1; mU3UQ j  
is=25;                                                                     %激励源x向位置 S"Efp/-  
js=25;                                                                     %激励源y向位置 ^|8cS0dK]Q  
ks=25;                                                                     %激励源z向位置 0{
j>u`  
iebc=8;                                                                    %左、右PML层层数 $)'{+1  
jebc=8;                                                                    %上、下PML层层数   -du+iOe?  
kebc=8;                                                                    %前、后PML层层数 #>233<  
ibbc=iebc+1; DF|qNX  
jbbc=jebc+1; 46 77uy  
kbbc=kebc+1; [iDa6mcth  
nmax=100;                                                                  %时间步 H fRxgA@  
dlta=6;                                                                    %观测点距激励源步数 cJqPcCq(wn  
'aCnj8B  
kh`X92~  
%************************************************ `xtN+y F  
Exy=zeros(ie,jb,kb); 5X  
Exz=zeros(ie,jb,kb); ALKhZFuz  
Eyx=zeros(ib,je,kb); o$Jk27  
Eyz=zeros(ib,je,kb); M0^r!f>O  
Ezx=zeros(ib,jb,ke); >LW9$[H  
Ezy=zeros(ib,jb,ke); 0xPML}|V  
Hxy=zeros(ib,je,ke); uZqo"  
Hxz=zeros(ib,je,ke); K,So#Ui  
Hyx=zeros(ie,jb,ke); `dj/Uk  
Hyz=zeros(ie,jb,ke); _z}d yp"I  
Hzx=zeros(ie,je,kb); <cl$?].RE!  
Hzy=zeros(ie,je,kb); E&97;VH  
sx/g5?zh  
m,*f6g  
fid1=fopen('m1.txt','wt'); LIR2B"3F  
fid2=fopen('m2.txt','wt'); \O^= Z{3y  
%************************************************ 6!bf,T
]  
sigxmax=-log(R0)/(2*iebc*ds)*3*epsz*c;                                     %x向最大电导率 [*1c.&%(  
sigymax=-log(R0)/(2*jebc*ds)*3*epsz*c;                                     %y向最大电导率 p}j{<y  
sigzmax=-log(R0)/(2*kebc*ds)*3*epsz*c; 9J>DLvl;  
Y&#<{j':  
 1 ft.ZJ  
for i=1:ib                                                                 %x向网格节点处电导率分布 B aXzz  
                if      i<=iebc Y(
&phv&  
                    sigx(i)=sigxmax*((iebc+1-i)/iebc)^2;                      %左侧PML层 "r[Ea|  
                elseif i>=ib-iebc+1                                           %右侧PML层 ln3
.TR*  
                    sigx(i)=sigxmax*((i-ib+iebc)/iebc)^2; UboOIx5:  
                else [%b<%m}L-  
                    sigx(i)=0;                                              %中间自由空间 nrZv>r  
                end     i4->XvC  
      cax(i)=(1-(dt*sigx(i))/(2*epsz))/(1+(dt*sigx(i))/(2*epsz)); Tp9LBF  
      cbx(i)=(dt/epsz/ds)/(1+(dt*sigx(i))/(2*epsz)); n
%ld*EgY  
       pHWol!  
end   @]OI(B  
for i=1:ie                                                                 %x方向网格边中点处电导率分布 T*$uc,  
                if      i<=iebc #Q;#A |EZ  
                    sigxh(i)=sigxmax*((iebc+1-i-0.5)/iebc)^2;                 %左侧PML层 %``FIv15w  
                elseif i>=ie-iebc+1                                           %右侧PML层 :}E*u^v K  
                    sigxh(i)=sigxmax*((i-ie+iebc-0.5)/iebc)^2; l]%|w]i\  
                else jSddjs  
                    sigxh(i)=0;                                              %中间自由空间 `_f3o,5  
                end                 KYlWV<sR  
      roux(i)=muz*sigxh(i)/epsz;                                              %磁导率       OnG!5b  
      cpx(i)=(1-(dt*roux(i))/(2*muz))/(1+(dt*roux(i))/(2*muz)); <1hwXo  
      cqx(i)=(dt/muz/ds)/(1+(dt*roux(i))/(2*muz)); D]4?
UL  
end   wv1?v_4  
~M <4HC  
oiklRf  
4=1lyw  
for j=1:jb                                                                 %y方向网格节点处电导率分布 b?r0n]  
                if      j<=jebc [7$<sN<'  
                    sigy(j)=sigymax*((jebc+1-j)/jebc)^2;                      %下侧PML层 0m?ul%=  
                elseif j>=jb-jebc+1                                           %上侧PML层 uH]^/'8vBd  
                    sigy(j)=sigymax*((j-jb+jebc)/jebc)^2; @m(\f  
                else I{M2nQi  
                    sigy(j)=0;                                              %中间自由空间 xvgIYc{  
                end               3lKIEPf6r  
      cay(j)=(1-(dt*sigy(j))/(2*epsz))/(1+(dt*sigy(j))/(2*epsz)); ?Ww',e  
      cby(j)=(dt/epsz/ds)/(1+(dt*sigy(j))/(2*epsz)); 5xRh'Jkyb  
end   YrB-;R1+  
for j=1:je                                                                 %y方向网格边中点处电导率分布                 !'+t)h9^  
                if      j<=jebc M>0~Ek%3  
                    sigyh(j)=sigymax*((jebc+1-j-0.5)/jebc)^2;                 %下侧PML层 Tvk=NJ  
                elseif j>=je-jebc+1                                           %上侧PML层 TsR20P@  
                    sigyh(j)=sigymax*((j-je+jebc-0.5)/jebc)^2; 0wOgQ
n  
                else Ir]b.6B  
                    sigyh(j)=0;                                              %中间自由空间 Y\j &84  
                end 'dBzv>ngD  
      rouy(j)=muz*sigyh(j)/epsz;                                              %磁导率       u<+;]8[o  
      cpy(j)=(1-(dt*rouy(j))/(2*muz))/(1+(dt*rouy(j))/(2*muz)); | WDX@Q  
      cqy(j)=(dt/muz/ds)/(1+(dt*rouy(j))/(2*muz));         "+|>nA=7  
end   E6n;_{Se/S  
Cu!4ha.e`  
8y+Gvk:  
P~?u2,.E[  
for k=1:kb                                                                 %z方向网格节点处电导率分布 xNjA>S\]W5  
                if      k<=kebc :pNZQX  
                    sigz(k)=sigzmax*((kebc+1-k)/kebc)^2;                      %前侧PML层 #?aR,@n  
                elseif k>=kb-kebc+1                                           %后侧PML层 (L~3nN;rr  
                    sigz(k)=sigzmax*((k-kb+kebc)/kebc)^2; 8o~
\L= l  
                else ]H.+=V;1  
                    sigz(k)=0;                                              %中间自由空间 F.O2;M|x  
                end (spX3n%p  
      caz(k)=(1-(dt*sigz(k))/(2*epsz))/(1+(dt*sigz(k))/(2*epsz)); Jq.26I=  
      cbz(k)=(dt/epsz/ds)/(1+(dt*sigz(k))/(2*epsz));       4$_8#wB1&  
end   Ju:=-5r"'  
for k=1:ke                                                                 %z方向网格边中点处电导率分布                 1-q\C<Q)  
                if      k<=kebc gg6&Fzp  
                    sigzh(k)=sigzmax*((kebc+1-k-0.5)/kebc)^2;                 %前侧PML层 N["(ZSS   
                elseif k>=ke-kebc+1                                           %后侧PML层 {0e5<"i  
                    sigzh(k)=sigzmax*((k-ke+kebc-0.5)/kebc)^2; 2wu 5`Z[E  
                else mV^dIm  
                    sigzh(k)=0;                                              %中间自由空间 1P6~IZVN  
                end &npf %Eub  
      rouz(k)=muz*sigzh(k)/epsz;                                              %磁导率     @ cv`}k  
      cpz(k)=(1-(dt*rouz(k))/(2*muz))/(1+(dt*rouz(k))/(2*muz)); pKp#4Js  
      cqz(k)=(dt/muz/ds)/(1+(dt*rouz(k))/(2*muz));        
SLBKXj|  
end   |rNm_L2  
5ptbz<Xv  
Ef7Kx49I  
% Ex=Exy+Exz; Z5NuLB'  
% Ey=Eyx+Eyz; sX@e1*YE_  
% Ez=Ezy+Ezx; -'ZP_$sA  
D@\97t+  
%************************************************ !WDdq_n*v  
tic c5Offnq'1  
R4{}ZT  
schedule_box=waitbar(0,'Please wait...'); % display the schedule box s2v\R~T  
                                                                       ukri7 n
*  
for n=1:nmax     DNL TJrN  
Q>||HtF$A  
           for i=1:ib )m6=_q5@o  
               for j=1:je M?]ObIM:5  
                 for k=1:ke k"AY7vq@!P  
                    if(i==is&&j==js&&k==ks) +!w?g/dV  
                        Hxy(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 ZtIK"o-|!  
                    else X2o5Hc)l<  
                        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)); 6#dx%TC  
                    end 4CNK ]2  
                                         3fBq~Q  
                 end jQf1h|e  
               end !M}&dW2  
           end Ot v{#bB$  
             0k3
^+#J  
           lJq %me;4m  
            for i=1:ib KX*e2 /0  
               for j=1:je `.><$F  
                 for k=1:ke aIkxN&  
                   if(i==is&&j==js&&k==ks)  eYS  
                       Hxz(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 $|AvT;4  
                   else UY>{e>/H9  
                       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)); P^&+ehp  
                   end bZa?h.IF  
%                    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)); ~PS%^zxyn  
               end E4 JS   
               end pvcf_w`n
 
          end ~~h9yvW7&  
%           Hx=Hxy+Hxz; igz&7U8gg  
          for i=1:ie {'{ssCL  
              for j=1:jb *6k (xL  
                 for k=1:ke "zm.jNn  
                   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)); 6`EyzB%.$  
                 end .llAiv  
               end ujDAs%6MZ  
           end s;$ eq);  
          for i=1:ie                                                 Hjlx,:'M  
               for j=1:jb mB_ba1r  
                 for k=1:ke LG51e7_gFi  
                   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)); `t#C0  
               end %f?#) 01>  
               end zYH6+!VBH#  
          end {K:/(\  
%           Hy=Hyx+Hyz; ;9b?[G  
          for i=1:ie                                                   7rsrC  
               for j=1:je $k}+,tHtJO  
                 for k=1:kb
I /RvU,  
                   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)); i"_JF-IbN  
               end rs\*$20  
               end MJ>(HJY6?%  
           end & yw-y4 =  
          for i=1:ie                                                   4?8GK  
               for j=1:je W*VQ"CW{^]  
                 for k=1:kb XjL( V1  
                   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)); Vd|/]Zj  
               end %#|S  
               end M*@MkN*u&  
          end g[!sGa&  
% Hz=Hzx+Hzy; V GM/ed5-  
Be?mIwc_g  
          for i=1:ie                                                   $^`hu%s,~  
               for j=2:je  % PEC at j=1 & j=jb qOkw6jfluh  
                 for k=2:ke  % PEC at k=1 & k=kb  后面相同 I7]45pF  
                   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)); drF"kTD"7  
                 end r`6XF  
               end D|UDLaz~  
          end <:/V`b3a  
          for i=1:ie                                               P`RM"'Om  
               for j=2:je \#~~,k 6f  
                 for k=2:ke <o p !dS  
                   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)); }8 ,b;Q  
                 end H2|w  
               end v82@']IN  
          end n j1 cqh  
7dxY07yu  
            for i=2:ie                                                     :qw:)i  
               for j=1:je \b~zyt6-  
                 for k=2:ke {T.
$xiR  
                   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)); %K?~$;Z.  
                 end T*LbZ"A  
               end jj.)$|&#`  
            end x4fLe5xv  
             NcqE)"yObo  
           ..  M3u[E  
vO <;Gnh~  

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

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

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

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

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

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

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

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

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

你是直接做的3D的吗？