三维pml matlab 程序，有问题，大侠帮忙啊，急！

超过48小时，复制至相关版块，此帖锁定
参与讨论

经过进一步检验，这是修改之后的程序，不过仍存在问题：1、激励源的加法，我的激励源修改后分别加到两个分量上，但距离源点几个网格后仍产生毛刺，说明源的加入还是有问题；2、我的程序是把自由空间和pml层一起做了运算和迭代，其实边角问题已经不存在，不需要考虑了（二维已经验证了），这样编程的好处是缩短程序，全域统一编程。 +~$pkxD"  
现在的程序仍有毛刺，请高手帮我改正，万分感激！！！谢谢！！ CUnBi?Mi  
clear b\S~uFq6  
c=3e8;   d
J{q}U  
muz=4*pi*1e-7;   apgR[=Oy  
epsz=1/(c*c*muz);       weH3\@  
f=1e+9;               HCw,bRxm  
lambda=c/f;             ] @:x<>  
omega=2*pi*f; ,c,@WQ2
:-  
ds=lambda/12;             L;-V Yo#  
dt=ds/(2*c); N6HeZB":  
muz=4*pi*1e-7;   SW}?y%~  
R0=1e-5;   )&j@={0  
8d7 NESYl  
%************************************************ g OK
 
ie=50;                                                                    %x向空间步             "0 $UnR  
je=50;                                                                    %y向空间步       oA?EJ~%  
ke=50;                                                                    %z向空间步 8j)*T9  
ib=ie+1; ly#jl5wmT  
jb=je+1; 8.:WMH`  
kb=ke+1; yoH,4,!G  
is=25;                                                                     %激励源x向位置 Kay\;fXT  
js=25;                                                                     %激励源y向位置 +c$:#9$ |  
ks=25;                                                                     %激励源z向位置 )4TP{tp  
iebc=8;                                                                    %左、右PML层层数  ]0XlI;ah  
jebc=8;                                                                    %上、下PML层层数   b|-S;cw  
kebc=8;                                                                    %前、后PML层层数 m*.+9 6  
ibbc=iebc+1; _:]g:F[ #  
jbbc=jebc+1; tb4^+&.GS  
kbbc=kebc+1; ERy=lP~gV  
nmax=100;                                                                  %时间步 zGNmc7  
dlta=6;                                                                    %观测点距激励源步数 K /$-H#;N  
IGOEqUw*  
YQcaWd(  
%************************************************ _#
qfe  
Exy=zeros(ie,jb,kb); [1nUq!uTm  
Exz=zeros(ie,jb,kb); SBI*[
 
Eyx=zeros(ib,je,kb); Xe&p.v  
Eyz=zeros(ib,je,kb); H?rCIS0  
Ezx=zeros(ib,jb,ke); IN75zn*%  
Ezy=zeros(ib,jb,ke); [RF6mWQ  
Hxy=zeros(ib,je,ke); @H8DGeM
 
Hxz=zeros(ib,je,ke); wjfq"7Q  
Hyx=zeros(ie,jb,ke); V8#NXUg<!  
Hyz=zeros(ie,jb,ke); Iz[ohn!f  
Hzx=zeros(ie,je,kb); Lg~ll$ U  
Hzy=zeros(ie,je,kb); 1obajN  
iK=QP+^VN  
qw 03]a  
fid1=fopen('m1.txt','wt'); 6Yl+IP];i  
fid2=fopen('m2.txt','wt'); /0o#V-E)  
%************************************************ Zo,066'+[.  
sigxmax=-log(R0)/(2*iebc*ds)*3*epsz*c;                                     %x向最大电导率 adPd}rt;  
sigymax=-log(R0)/(2*jebc*ds)*3*epsz*c;                                     %y向最大电导率 < 0YoZSNGj  
sigzmax=-log(R0)/(2*kebc*ds)*3*epsz*c; R.cR:fA  
m(D+!I9  
k{H7+;_  
for i=1:ib                                                                 %x向网格节点处电导率分布 |nfMoUI  
                if      i<=iebc Cu!]-c{  
                    sigx(i)=sigxmax*((iebc+1-i)/iebc)^2;                      %左侧PML层 v[r8-0c  
                elseif i>=ib-iebc+1                                           %右侧PML层 (vp#?-i  
                    sigx(i)=sigxmax*((i-ib+iebc)/iebc)^2; 3m| C8:  
                else a7685Y  
                    sigx(i)=0;                                              %中间自由空间 j^%N:BQ&  
                end     n(`|:h"  
      cax(i)=(1-(dt*sigx(i))/(2*epsz))/(1+(dt*sigx(i))/(2*epsz)); &$ud;r#  
      cbx(i)=(dt/epsz/ds)/(1+(dt*sigx(i))/(2*epsz)); d<6m_!L  
       <_c8F!K)T  
end   %/ctt_p0x  
for i=1:ie                                                                 %x方向网格边中点处电导率分布 fW[ .Q0  
                if      i<=iebc >g m  
                    sigxh(i)=sigxmax*((iebc+1-i-0.5)/iebc)^2;                 %左侧PML层 7Ie=(x8):  
                elseif i>=ie-iebc+1                                           %右侧PML层 W>5[_d  
                    sigxh(i)=sigxmax*((i-ie+iebc-0.5)/iebc)^2; [X91nUz#  
                else ac\([F-  
                    sigxh(i)=0;                                              %中间自由空间 d8iq9AP\o  
                end                  O4Q"2  
      roux(i)=muz*sigxh(i)/epsz;                                              %磁导率       je5[.VTM  
      cpx(i)=(1-(dt*roux(i))/(2*muz))/(1+(dt*roux(i))/(2*muz)); |Sm/s;&c6  
      cqx(i)=(dt/muz/ds)/(1+(dt*roux(i))/(2*muz)); H'JU5nE  
end   o{hX?,4i  
{PR "}x  
A'.=SA2.Y  
.)SR3?  
for j=1:jb                                                                 %y方向网格节点处电导率分布 &B]1 VZUp  
                if      j<=jebc }
V[ORGzox  
                    sigy(j)=sigymax*((jebc+1-j)/jebc)^2;                      %下侧PML层 MT7B'hd  
                elseif j>=jb-jebc+1                                           %上侧PML层 f!{@
{\  
                    sigy(j)=sigymax*((j-jb+jebc)/jebc)^2; 3I(;c ,S  
                else lu8*+.V  
                    sigy(j)=0;                                              %中间自由空间 QYi4A"$`  
                end               RD46@Q`  
      cay(j)=(1-(dt*sigy(j))/(2*epsz))/(1+(dt*sigy(j))/(2*epsz)); w {"1V7|  
      cby(j)=(dt/epsz/ds)/(1+(dt*sigy(j))/(2*epsz)); b;%t*?t  
end   AVm+ 1  
for j=1:je                                                                 %y方向网格边中点处电导率分布                 X`1R&K;z^  
                if      j<=jebc Z/dhp0k  
                    sigyh(j)=sigymax*((jebc+1-j-0.5)/jebc)^2;                 %下侧PML层 )_1 GPS  
                elseif j>=je-jebc+1                                           %上侧PML层 JHXkQz[Jb  
                    sigyh(j)=sigymax*((j-je+jebc-0.5)/jebc)^2; g+5c"Yk+u~  
                else U E$Ix  
                    sigyh(j)=0;                                              %中间自由空间 e,&#,O  
                end tt#dO@G#Fe  
      rouy(j)=muz*sigyh(j)/epsz;                                              %磁导率       a9UXg<4  
      cpy(j)=(1-(dt*rouy(j))/(2*muz))/(1+(dt*rouy(j))/(2*muz)); C$0g2X  
      cqy(j)=(dt/muz/ds)/(1+(dt*rouy(j))/(2*muz));         rOz1tY)l0d  
end   bAbR0)  
S8Y\@C?5  
#wo *2(  
gq"d$Xh$x7  
for k=1:kb                                                                 %z方向网格节点处电导率分布 NSBcYObX  
                if      k<=kebc
:.r_4$F:  
                    sigz(k)=sigzmax*((kebc+1-k)/kebc)^2;                      %前侧PML层 {VKFw=$8  
                elseif k>=kb-kebc+1                                           %后侧PML层 /M+Du,  
                    sigz(k)=sigzmax*((k-kb+kebc)/kebc)^2; * k<@  
                else b\"w
/'XX  
                    sigz(k)=0;                                              %中间自由空间 yYaoA/0  
                end Nke!!A}\|  
      caz(k)=(1-(dt*sigz(k))/(2*epsz))/(1+(dt*sigz(k))/(2*epsz)); U`lK'..  
      cbz(k)=(dt/epsz/ds)/(1+(dt*sigz(k))/(2*epsz));       J/O{x  
end   & +*OV:[;  
for k=1:ke                                                                 %z方向网格边中点处电导率分布                 {}$Zff  
                if      k<=kebc mBE&>}G<  
                    sigzh(k)=sigzmax*((kebc+1-k-0.5)/kebc)^2;                 %前侧PML层 |JP19KFx'B  
                elseif k>=ke-kebc+1                                           %后侧PML层 ,uAp;"YJeV  
                    sigzh(k)=sigzmax*((k-ke+kebc-0.5)/kebc)^2; L
SP p  
                else TL)*onA9  
                    sigzh(k)=0;                                              %中间自由空间 &!OEd]  
                end B]@25  
      rouz(k)=muz*sigzh(k)/epsz;                                              %磁导率     hHGuD2%  
      cpz(k)=(1-(dt*rouz(k))/(2*muz))/(1+(dt*rouz(k))/(2*muz)); </WeB3#6  
      cqz(k)=(dt/muz/ds)/(1+(dt*rouz(k))/(2*muz));         &w#!
 
end   'E+"N'M|  
yu)^s!UY;  
vbVOWX6  
% Ex=Exy+Exz; N0.|Mb"?t  
% Ey=Eyx+Eyz; #c
5jCy}n  
% Ez=Ezy+Ezx; B\v+C!/f|  
.] sJl  
%************************************************ 1
5,JD  
tic jj1\
oyQ8  
}f]Y^>-Ux  
schedule_box=waitbar(0,'Please wait...'); % display the schedule box tq}45{FH3  
                                                                       3+15 yEeA  
for n=1:nmax     m3TR}=n  
pF4Z4?W  
           for i=1:ib BHf$ %?3z,  
               for j=1:je s2#Ia>5!  
                 for k=1:ke ]"lB!O~  
                    if(i==is&&j==js&&k==ks) |1RVm?~i  
                        Hxy(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 IHYLM;@L  
                    else T!8^R|!a6  
                        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)); ATl?./Tu  
                    end *<\K-NSL  
                                         xC,x_:R`  
                 end ;4[[T%&v  
               end ~Ix2O  
           end fEX=csZ86  
             KWZhCS?[(  
           e=WjFnK[x7  
            for i=1:ib %VH,(}i  
               for j=1:je SkE<V0  
                 for k=1:ke aL( hWE  
                   if(i==is&&j==js&&k==ks) 3f] ;y<K
m  
                       Hxz(i,j,k)=0.5*sin(omega*n*dt);                                                                        %源激励 D%abBE1  
                   else o`,~#P|  
                       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)); 8.[F3Tk=  
                   end iN[x *A|h  
%                    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)); B*,)
@h  
               end 68Gywk3]=u  
               end Q-n8~Ey1a  
          end -}9^$}PR  
%           Hx=Hxy+Hxz; 1^4:l!0D  
          for i=1:ie NS~;{d\  
              for j=1:jb D2?H"PH  
                 for k=1:ke ^>?=L\[  
                   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)); /\c'kMAW!  
                 end +yp:douERi  
               end kIVQ2hmv  
           end H*'1bLzq  
          for i=1:ie                                                 pA6KiY&  
               for j=1:jb }=5>h' <  
                 for k=1:ke ,|{`(y/v  
                   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)); RqtBz3v  
               end !Aw^X} C  
               end ]xr0]  
          end WJ25fTsG  
%           Hy=Hyx+Hyz; RZzHl
Z  
          for i=1:ie                                                   y %Q.(  
               for j=1:je %bAQ>E2;m  
                 for k=1:kb [mA-s
l]  
                   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)); A^>@6d $2  
               end qcS.=Cj?)  
               end wQSye*ec  
           end 
kFv*>>X`  
          for i=1:ie                                                   ];OvV ,*  
               for j=1:je <qwf"Ey  
                 for k=1:kb ZpV]X(Px(o  
                   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)); e@Lxduq  
               end -GAF
>  
               end y|2<Vc  
          end
cm(*F0<  
% Hz=Hzx+Hzy; AJbCC  
n^Ca?|} ,  
          for i=1:ie                                                   +e-F`k  
               for j=2:je  % PEC at j=1 & j=jb UX@%1W!8  
                 for k=2:ke  % PEC at k=1 & k=kb  后面相同 9%"7~YCDas  
                   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)); #wI}93E  
                 end 94rSB}b.O  
               end D\AVZ76F1  
          end a%T`c/C  
          for i=1:ie                                               P.'.KZJ:WD  
               for j=2:je u^~7[OkE  
                 for k=2:ke %.Ma_4o Z  
                   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)); V4n~Z+k  
                 end #i[:oC6m:  
               end JaCX}[R  
          end @-'a{hBR  
iT>u&0B-  
            for i=2:ie                                                     )*~A|[  
               for j=1:je V4:/LNq_]  
                 for k=2:ke tWIs |n  
                   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)); Y~WdN<g  
                 end H~a ~'tm  
               end )l7XZ_gw'  
            end JWixY/  
             \8`7E1d  
         .. s-$Wc)l  
7$'AH:K  

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

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

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

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

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

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

如果都没有迭代的话，自己调试一下吧。 ,p(&G_  
Hxy(is,js,ks)是源，那么Eyz和Ezy应该是可以得到值的，你把这几个分量都实时显示出来，看看问题到底在哪里吧。 |S|'o*u  
gVR]z9  
或者你给我个email，我发一个2D+PML的源码给你参考一下。

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

超过48小时，复制至相关版块，此帖锁定 T:Klr=&V  
%KT}Map  
参与讨论