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

#include <stdio.h> `b# w3 2  
#include <math.h> ?N&"WL^|  
#include <fdtd.h> !@p@u;djJ  
#include <string.h> [ wr0TbtV  
#define IE 120 Xp4pN{
he  
#define JE 120 u`R  
xa5I{<<U  
void main(void) AIN_.=]"?  
{ P_lk40X  
 /* /J wQ5  
  *葛德彪老师书中的4个参数:ca,cb,cp,cq, 放在结构体param0里面 ! FhN(L[=j  
  */ gV$Lfkz  
 param param0; w3fi2B&q  
 /* $ [M8G  
  *葛德彪老师书中的4个参数:ca,cb,cp,cq, 放在结构体param0里面 %(,JBa:G  
  */ m:ITyQ+  
 abcs_2d_mur_te abcs_mur; E.}T.
St  
 double hz[IE][JE],ex[IE][JE],ey[IE][JE]; q#c\  
 double hz_last_time[IE][JE],ex_last_time[IE][JE],ey_last_time[IE][JE]; | t3_E  
 /************************************************************************************************ U3ED3) D  
  * 源的设置 rF>:pS,`&  
  */ 7f~.Qus  
 /* p`T7Y\\#!  
  *高斯源放置位置 <+ >y GPp  
  */ CdC&y}u  
 int gauss_pulse_position_x = IE/2; ERCW5b[RT  
 int gauss_pulse_position_y = JE/2; rP^2MH"  
 /* `?2S4lN/  
  *设置高斯脉冲源参数，极大值出现在x=12.0，脉冲宽度为6.0 MJ4+|riB  
  */ Mpb|qGi!  
 double gauss_pulse_t0 = 12.0; KokmylHu  
 double gauss_pulse_width= 6.0; :&9#p%/  
 /* .y#@~H($  
  *设置正弦源参数，sinusoidal_frequency:正弦源的频率 f.,S-1D]h  
  */ <eQS16  
 double sinusoidal_frequency = 5.0e9; Nq9@^ E-{M  
 /* (VU: &.  
  *设置源通用参数，pulse_source：源的大小；source_time：源的时间 `hY%<L sI  
  */
HwiG~'Ah9  
 double pulse_source = 0; 1
g^N7YF  
 int source_time =0; *^}(LoPZ  
 /************************************************************************************************/ )LKutN?tBy  
 /* {LrezE4  
  * 时间步数 {[NBTT9&  
  */ X3<<f`X  
 int time_steps = 200; +p>tO\mo  
 int n = 0,i = 0,j=0; !1-:1Whz8  
 /* xpRQ"6  
  * 设置边界 )Uo)3FAn
 
  */ {6 C!^ 5  
 int left_boundary = 1; Jd7+~isu~  
 int right_boundary = IE-1; /2~qm/%Q  
 int top_boundary = JE-1; e#<%`\qH  
 int bottom_boundary = 1; = 0
Z}s  
 /* doW_vu  
  * 设置介电常数 j2k,)MHu!x  
  */ Rm&

i"  
 fdtd_init(); 118A6qyi  
 /* @'7'3+ c  
  * 设置迭代方程参数ca，cb ，cp，cq R*'rg-
d  
  */ M(jSv  
 init_param(&param0); |9@?8\  
 /* _@ev(B  
  * 设置边界方程的参数 6zaO$  
  */ W(9-XlYKE  
 init_abcs_2d_mur_te(&abcs_mur); }[l
d=9p(  
;
vy"i
 
 for(j  = 0;j < JE ;j++) L(fOe3 v  
 { 3ktjMVy\  
  for(i = 0;i < IE; i++) |'aGj  
  { L 1H!o!*  
   hz[j] = 0.0; SRRqIQz  
   ex[j] = 0.0; |~Z.l  
   ey[j] = 0.0; @aAB#,  
  } }i^]uW*h  
 } zn@<>o8hU  
}~DlOvsq  
 for(n=0;n < time_steps;n++) Pv|g.hH9m  
 { +C' u!^)  
   o:W>7~$jr=  
  source_time += 1;
@@-n/9>vs  
  /* "jb`KBH%"  
   * 放置高斯脉冲源 TWZ**S-  
   */ ;:4&nJ*qG  
  //pulse_source = Gauss_Pulse(source_time,gauss_pulse_t0,gauss_pulse_width); EZ:pcnL{  
  /* m(i84~  
   * 放置高斯脉冲源 1hp`.!3]H  
   */ 2!Yq9,`  
  pulse_source = sinusoidal_source(sinusoidal_frequency,((double)source_time)*dt); ]w`)"{j5m  
  hz[gauss_pulse_position_x][gauss_pulse_position_y] = pulse_source; ikX"f?Q;S2  
  memcpy(hz_last_time,hz,sizeof(hz)); o$;t  
  memcpy(ex_last_time,ex,sizeof(ex)); ^~9fQJNs  
  memcpy(ey_last_time,ey,sizeof(ey)); q^; SZ^yW5  
1LYz X;H1  
  /* ZP~Mgz{f  
   * 计算磁场 [ R  
   */ e=Ko4Ao2y  
  for(j = bottom_boundary+1;j <= top_boundary-1;j++) 3qujz)o  
  { cvxY
uP~  
   for(i = left_boundary+1;i <= right_boundary-1;i++) T$rhz)_q  
   { VwvL  
    hz[j] = param0.cp*hz[j] - param0.cq*((ey[i+1][j]- ey[j])/dx - (ex[j+1] - ex[j])/dy); BbsgZ4  
   } FjR/_GPo6  
  } @6h,#8#  
  /* >9dzl#  
   * 计算电场 DcYL8u  
   */  >M~1{  
  for(j = bottom_boundary+1;j <= top_boundary;j++) J=7.-R|t  
  { iZY4+
X  
   for(i = left_boundary;i <= right_boundary;i++) !ij R  
   { \]<eLw-v  
    ex[j] = param0.ca*ex[j] + param0.cb*(hz[j] - hz[j-1])/dy; 5|O~  
   } J>Pc@,y  
  }    niC ;WK  
  for(j = bottom_boundary;j <= top_boundary;j++) f#mNx  
  { 7>.d*?eao\  
   for(i = left_boundary+1;i <= right_boundary;i++) vz6No%8X  
   { J;HYGu:  
    ey[j] = param0.ca*ey[j] - param0.cb*(hz[j] - hz[i-1][j])/dx; ]bxBo  
   }
YYNh| 2  
  } !ZNirvk  
  /* 2_T2?weD5  
   * 加入Mur二阶边界条件 8u8-:c%{  
   */ <<'%2q5  
  for(j = bottom_boundary+1;j <= top_boundary-1;j++) lCafsIB  
  { CWT#1L=  
   //x=1，左边界 `]~1pc  
   hz[left_boundary ][j] = hz_last_time[left_boundary +1][j] H=lzW_(  
         + abcs_mur.left_h1*(hz[left_boundary +1][j] - hz_last_time[left_boundary ][j]) d
e=){.7Y  
         + 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]); ^^,
cnDlm  
   //x=h，右边界 d?:KEi-<7  
   hz[right_boundary][j] = hz_last_time[right_boundary-1][j] T<? (KW  
         + abcs_mur.left_h1*(hz[right_boundary-1][j] - hz_last_time[right_boundary][j]) D/[;Y<X#V  
         + 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]); =hRo#]{(K  
  } S)LvYOOB@  
     for(i = left_boundary+1;i <= right_boundary-1;i++) 5 w-Pq&q  
  { ^?M# |>  
   //y=1，下边界 (aC~0 #4  
   hz[bottom_boundary] = hz_last_time[bottom_boundary+1] =/5^/vwgY  
         + abcs_mur.left_h1*(hz[bottom_boundary+1] - hz_last_time[bottom_boundary]) <+y%k~("  
         + 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]); lXpbAW  
   //y=h，上边界 cN% r\  
   hz[ top_boundary ] = hz_last_time[  top_boundary-1 ] i?wEd!=w  
         + abcs_mur.left_h1*(hz[  top_boundary-1 ] - hz_last_time[ top_boundary  ]) b:WA}x V  
         + 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 ]); 8:t!m>(*  
  } 2# 72B  
  //角点的处理 .;\uh$c  
  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]); pO?v$Rjl  
  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]); `T\_Wje(  
  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]); knfEbH  
  hz[right_boundary][top_boundary   ] = hz_last_time[right_boundary-1][ top_boundary-1  ] + abcs_mur.C2*(hz[right_boun .. ?e{hidg  
CdZ. T/x  

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

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

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

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

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