[求助]帮小弟看个问题，是fdtd 二维TE波，mur吸收边界

#include <stdio.h> y+:<  
#include <math.h> 0/KNXz  
#include <fdtd.h> n? s4"N6  
#include <string.h> dbZPt~S'$  
#define IE 120 K0I-7/L  
#define JE 120 )kUq2-r  
>Ft jrEB  
void main(void) 3!$rp- !<)  
{ E/GI:}YUy_  
 /* b'-g
y0  
  *葛德彪老师书中的4个参数:ca,cb,cp,cq, 放在结构体param0里面 ~`N|sI,  
  */ hV8A<VT  
 param param0; [=E<iPl  
 /* cRg$~rYd  
  *葛德彪老师书中的4个参数:ca,cb,cp,cq, 放在结构体param0里面 % 9WWBxS  
  */ 8'cDK[L  
 abcs_2d_mur_te abcs_mur; \W@?revK  
 double hz[IE][JE],ex[IE][JE],ey[IE][JE]; mI`dZ3h  
 double hz_last_time[IE][JE],ex_last_time[IE][JE],ey_last_time[IE][JE]; ]c5GG!E-g  
 /************************************************************************************************ %)aDh }  
  * 源的设置 \aW5V:?  
  */ "J{,P9P6  
 /* x@k9]6/zs  
  *高斯源放置位置 4t8 Hy  
  */ o!H"~5Trv!  
 int gauss_pulse_position_x = IE/2; fCV
SVn"o  
 int gauss_pulse_position_y = JE/2; x`eYCi  
 /* *smo{!0Gg  
  *设置高斯脉冲源参数，极大值出现在x=12.0，脉冲宽度为6.0 (~#PzE:  
  */ P7
Y[?='v  
 double gauss_pulse_t0 = 12.0; cL WM]\Y  
 double gauss_pulse_width= 6.0; 3S5Qq
Am  
 /* Q'=!1^&  
  *设置正弦源参数，sinusoidal_frequency:正弦源的频率 i([A8C_A  
  */ *@Qt*
f  
 double sinusoidal_frequency = 5.0e9; ZK'I$p]b  
 /* >$"bwr}'4B  
  *设置源通用参数，pulse_source：源的大小；source_time：源的时间 ,vhR99g{  
  */ pI^n("|  
 double pulse_source = 0; X>wQYIi  
 int source_time =0; uw<Ruy  
 /************************************************************************************************/ \dc`}}Lc  
 /* ]CL70+[^9  
  * 时间步数 gh 0\9;h  
  */ LQQhn{[D  
 int time_steps = 200; T@S+5(  
 int n = 0,i = 0,j=0; z`_N|iEd  
 /* 7PvuKAv?k  
  * 设置边界 '",5Bu#C  
  */ ,,iQG' *  
 int left_boundary = 1; ?
QMs<  
 int right_boundary = IE-1; !{3pp  
 int top_boundary = JE-1; 7xd}J(l  
 int bottom_boundary = 1; L'6zs:i  
 /* E|fPI u  
  * 设置介电常数  :D/R  
  */ Z2@&4_P  
 fdtd_init(); WMC6dD_6e  
 /* [,e_2<  
  * 设置迭代方程参数ca，cb ，cp，cq }gn0bCJy  
  */ eX$Biv1N  
 init_param(&param0); IWQ8e$N  
 /* ``\H'^{B  
  * 设置边界方程的参数 kR_[p._  
  */ ~gAx  
 init_abcs_2d_mur_te(&abcs_mur); O,;SA  
P~*fZ)\}F@  
 for(j  = 0;j < JE ;j++) lubS{3<  
 { z DK+8  
  for(i = 0;i < IE; i++) e'?(`yW>  
  { 9(I4x]`  
   hz[j] = 0.0; lk'RWy"pw  
   ex[j] = 0.0; l-npz)EM  
   ey[j] = 0.0; Ar$LA"vu4  
  } ~lL($rE  
 } p*'?(o:=  
s-DtkO  
 for(n=0;n < time_steps;n++) OL=bhZ  
 { #nd,cn  
   s6 ^JgdW  
  source_time += 1; WgL!@
g  
  /* &KB{,:)?  
   * 放置高斯脉冲源 D *tBbV  
   */ 5u!
cA4e"  
  //pulse_source = Gauss_Pulse(source_time,gauss_pulse_t0,gauss_pulse_width); .KD0
7  
  /* }qg!U
m0  
   * 放置高斯脉冲源 bd9c/>&  
   */ MWuVV=rd8a  
  pulse_source = sinusoidal_source(sinusoidal_frequency,((double)source_time)*dt); Tv KX8m"  
  hz[gauss_pulse_position_x][gauss_pulse_position_y] = pulse_source; +l#2u#e  
  memcpy(hz_last_time,hz,sizeof(hz)); - t+Mh.  
  memcpy(ex_last_time,ex,sizeof(ex));
,8 .`;  
  memcpy(ey_last_time,ey,sizeof(ey)); =
osj}(  
~+g5?y  
  /* 4H:WpW*r  
   * 计算磁场 TvP# /qGgG  
   */ #-`lLI:w0  
  for(j = bottom_boundary+1;j <= top_boundary-1;j++) BOG )JaDW  
  { prWid3}  
   for(i = left_boundary+1;i <= right_boundary-1;i++) _Dv^~e1c  
   { nk{1z\D{  
    hz[j] = param0.cp*hz[j] - param0.cq*((ey[i+1][j]- ey[j])/dx - (ex[j+1] - ex[j])/dy);
v7R&9kU{  
   } @PI\.y_w  
  } E5dXu5+ye  
  /* 0Nfj}sXCWE  
   * 计算电场 Ob6vg
^#  
   */ Q9'p2@Z  
  for(j = bottom_boundary+1;j <= top_boundary;j++)  96BMJE'  
  { _`\INZe-G  
   for(i = left_boundary;i <= right_boundary;i++) 4j*}|@x  
   { Zry>s0  
    ex[j] = param0.ca*ex[j] + param0.cb*(hz[j] - hz[j-1])/dy; )<|TEp4r-  
   } o?3R HP47  
  }    )eFK@goGeb  
  for(j = bottom_boundary;j <= top_boundary;j++) y/FisX  
  { g[$B90  
   for(i = left_boundary+1;i <= right_boundary;i++) o6r4tpiR5  
   { Y:a(y*y<  
    ey[j] = param0.ca*ey[j] - param0.cb*(hz[j] - hz[i-1][j])/dx; j0GI[#  
   } Vh<`MS0X  
  } 1m0':n Vdu  
  /* JjmL6(*ui  
   * 加入Mur二阶边界条件 TQou.'+v  
   */ 4|(?Wt)5  
  for(j = bottom_boundary+1;j <= top_boundary-1;j++) :U9R 1^}A  
  { ![ QQF|  
   //x=1，左边界 3%} Ma,  
   hz[left_boundary ][j] = hz_last_time[left_boundary +1][j] xEBjf
n  
         + abcs_mur.left_h1*(hz[left_boundary +1][j] - hz_last_time[left_boundary ][j]) \x!>5Z Y  
         + abcs_mur.left_h2*(ex_last_time[left_boundary ][j] - ex_last_time[left_boundary ][j-1] + ex_last_time[left_boundary+1][j] - ex_last_time[left_boundary+1][j-1]); %L/=heBBd  
   //x=h，右边界 NR*SEbUU*  
   hz[right_boundary][j] = hz_last_time[right_boundary-1][j] bm4W,  
         + abcs_mur.left_h1*(hz[right_boundary-1][j] - hz_last_time[right_boundary][j]) cNVdGY%&  
         + abcs_mur.left_h2*(ex_last_time[right_boundary][j] - ex_last_time[right_boundary][j-1] + ex_last_time[right_boundary-1][j] - ex_last_time[right_boundary-1][j-1]); X_qXH5^%  
  } ^WNJQg'  
     for(i = left_boundary+1;i <= right_boundary-1;i++) =6sXZ"_Tw  
  { S|[UEU3FpB  
   //y=1，下边界 ST3qg6Cq2J  
   hz[bottom_boundary] = hz_last_time[bottom_boundary+1] (TE2t7ab|M  
         + abcs_mur.left_h1*(hz[bottom_boundary+1] - hz_last_time[bottom_boundary]) ?Z=v&d[o)  
         + abcs_mur.left_h2*(ey_last_time[bottom_boundary] - ey_last_time[i-1][bottom_boundary] + ey_last_time[bottom_boundary+1] - ey_last_time[i-1][bottom_boundary+1]); H@zk8]_P  
   //y=h，上边界 o;[oy#aWl_  
   hz[ top_boundary ] = hz_last_time[  top_boundary-1 ] X=3@M_Jzo  
         + abcs_mur.left_h1*(hz[  top_boundary-1 ] - hz_last_time[ top_boundary  ]) WdT|xf.Q&  
         + abcs_mur.left_h2*(ey_last_time[  top_boundary ] - ey_last_time[i-1][  top_boundary ] + ey_last_time[  top_boundary-1 ] - ey_last_time[i-1][  top_boundary-1 ]); /\#5\dHj  
  } (Ka#6   
  //角点的处理 c;^J!e  
  hz[left_boundary ][bottom_boundary] = hz_last_time[left_boundary+1 ][bottom_boundary+1] + abcs_mur.C2*(hz[left_boundary+1 ][bottom_boundary+1] - hz_last_time[left_boundary ][bottom_boundary+1]); |(SW  
  hz[right_boundary][bottom_boundary] = hz_last_time[right_boundary-1][bottom_boundary+1] + abcs_mur.C2*(hz[right_boundary-1][bottom_boundary+1] - hz_last_time[right_boundary][bottom_boundary]); B4}XK=)  
  hz[left_boundary ][top_boundary   ] = hz_last_time[left_boundary+1 ][ top_boundary-1  ] + abcs_mur.C2*(hz[left_boundary+1 ][ top_boundary-1  ] - hz_last_time[left_boundary ][top_boundary]); dc05,Bz  
  hz[right_boundary][top_boundary   ] = hz_last_time[right_boundary-1][ top_boundary-1  ] + abcs_mur.C2*(hz[right_boun .. Bag#An1  
^6+x0[13  

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

谁来帮小弟看看，肯定是边界条件问题！主要看一下边界！

个人认为是你对二阶摩尔吸收边界的公式理解有问题，要分清楚是n时刻的E值还是n+1时刻的，可以引入中间变量来保存上一时刻的电场值

借贴问个问题，为什么边界条件没有E的分量呢，是不是不用考虑呢

你在程序中定义个中间变量试试看？