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

#include <stdio.h> VfkQc$/  
#include <math.h> NS){D7T  
#include <fdtd.h> 79Ur1-]/  
#include <string.h> 52.%f+Oa  
#define IE 120 .\mkgAlyaM  
#define JE 120 l`rO)7  
P1dFoQz  
void main(void) 9v?rNJs  
{ m*vz   
 /* z LZHVvL3  
  *葛德彪老师书中的4个参数:ca,cb,cp,cq, 放在结构体param0里面 [i8Ju  
  */ ~WKWx.ul  
 param param0; qflOi8  
 /* E'^ny4gL  
  *葛德彪老师书中的4个参数:ca,cb,cp,cq, 放在结构体param0里面 8f>v[SQ"  
  */ <$Djags,F  
 abcs_2d_mur_te abcs_mur; g5lK&-yu]  
 double hz[IE][JE],ex[IE][JE],ey[IE][JE]; eqk.+~^  
 double hz_last_time[IE][JE],ex_last_time[IE][JE],ey_last_time[IE][JE]; F0ylJ /E  
 /************************************************************************************************ @nOuFX4  
  * 源的设置 
Xa_:B\ic  
  */ {=6CL'_  
 /* : $N43_
Wb  
  *高斯源放置位置 K&BaGrR  
  */ L b-xc]  
 int gauss_pulse_position_x = IE/2; +zZ]Txb(  
 int gauss_pulse_position_y = JE/2; iHeu<3O  
 /* vQ8$C 3  
  *设置高斯脉冲源参数，极大值出现在x=12.0，脉冲宽度为6.0 z-^/<u1p  
  */ ta0;:o?/d  
 double gauss_pulse_t0 = 12.0; }fL ]}&  
 double gauss_pulse_width= 6.0; vDC
bD#.6  
 /* qAi:F=> X  
  *设置正弦源参数，sinusoidal_frequency:正弦源的频率 DWupLJpk;c  
  */ dpcU`$kt  
 double sinusoidal_frequency = 5.0e9; uLr-!T  
 /* \0.!al0  
  *设置源通用参数，pulse_source：源的大小；source_time：源的时间 J~)JsAXAI  
  */ ]T._TZ"  
 double pulse_source = 0; ? kCo/sW  
 int source_time =0; 5h6-aQU[  
 /************************************************************************************************/ !*PX
-  
 /* ce7$# #f  
  * 时间步数 T@mYHKu  
  */ 8@r>`c  
 int time_steps = 200;
%lGT|XrY  
 int n = 0,i = 0,j=0; 9qm'qx  
 /* 85BB{T;  
  * 设置边界 eN0lJ~  
  */ c3(0BSv  
 int left_boundary = 1; zQoJ8i>  
 int right_boundary = IE-1; m O"Rq5  
 int top_boundary = JE-1; jJAr #|  
 int bottom_boundary = 1; {v*X}`.h  
 /* {EJ+   
  * 设置介电常数 "qoJIwl#q  
  */ .p%V]Ka  
 fdtd_init(); aB{OXU}#  
 /* p~OX1RBI  
  * 设置迭代方程参数ca，cb ，cp，cq z`sW5K(A  
  */ Kh{_Bd
N  
 init_param(&param0); )3^#C
D  
 /* :PNhX2F  
  * 设置边界方程的参数 ,qj M1xkL$  
  */ `d4;T|f+=  
 init_abcs_2d_mur_te(&abcs_mur); K0\`0E
^,  
VVqpzDoXG  
 for(j  = 0;j < JE ;j++) &?# YjU"  
 { FpZ5@  
  for(i = 0;i < IE; i++) _ygdv\^Tet  
  { GW2v&Ul7(  
   hz[j] = 0.0; [a7S?%>Bh  
   ex[j] = 0.0; Rp !Rzl<  
   ey[j] = 0.0; ! q+>'Mt  
  } I5E+=.T*ar  
 } Y4N)yMSl"  
* yGlX
[  
 for(n=0;n < time_steps;n++) c=<^pCa9t1  
 { irMBd8WG  
   '![VA8  
  source_time += 1; AmK g;9LS  
  /* \O)u' Bu  
   * 放置高斯脉冲源 J9P\D!  
   */ .}Hs'co  
  //pulse_source = Gauss_Pulse(source_time,gauss_pulse_t0,gauss_pulse_width); vy5I#q(k  
  /* p1s|JI  
   * 放置高斯脉冲源 @B9#Hrc  
   */ \~.elKw<U  
  pulse_source = sinusoidal_source(sinusoidal_frequency,((double)source_time)*dt); V o%GO9b;  
  hz[gauss_pulse_position_x][gauss_pulse_position_y] = pulse_source; 6C!TXV'  
  memcpy(hz_last_time,hz,sizeof(hz)); eBqF@'DQ  
  memcpy(ex_last_time,ex,sizeof(ex)); Mo\nY5  
  memcpy(ey_last_time,ey,sizeof(ey)); Hk?E0.  
P_(<?0l  
  /* 3iw.yR  
   * 计算磁场 R)Dh;XA  
   */ E//*bmww  
  for(j = bottom_boundary+1;j <= top_boundary-1;j++) syk,e4:oA  
  { lHO.pN`2  
   for(i = left_boundary+1;i <= right_boundary-1;i++) :Yn{:%p  
   { EUS^Gtc  
    hz[j] = param0.cp*hz[j] - param0.cq*((ey[i+1][j]- ey[j])/dx - (ex[j+1] - ex[j])/dy); 0oo_m6ie&  
   } GK
&R.R]  
  } 1-PFM-  
  /* !J(6E:,b#  
   * 计算电场 Kir|in)r0  
   */ [Vj|fy4  
  for(j = bottom_boundary+1;j <= top_boundary;j++) Vk%W4P"l  
  { !X 0 (4^  
   for(i = left_boundary;i <= right_boundary;i++) Qm4cuV-0{  
   { mV}eMw  
    ex[j] = param0.ca*ex[j] + param0.cb*(hz[j] - hz[j-1])/dy;  /UtSZ(  
   } 'grb@+w(
 
  }    :-~x~ah-  
  for(j = bottom_boundary;j <= top_boundary;j++) 4Yd$RP  
  { p7[&H/  
   for(i = left_boundary+1;i <= right_boundary;i++) 0gr#<(  
   { }3R
:7N`,|  
    ey[j] = param0.ca*ey[j] - param0.cb*(hz[j] - hz[i-1][j])/dx; |`Be(  
   } *2Q x69`  
  } Pz-=Eq  
  /* 8xHjdQr  
   * 加入Mur二阶边界条件 yn04[PN2  
   */ x5,|kJ9S  
  for(j = bottom_boundary+1;j <= top_boundary-1;j++) ]#dZLm_  
  { _#w5hXcu  
   //x=1，左边界 =<U'Jtu6'  
   hz[left_boundary ][j] = hz_last_time[left_boundary +1][j] 9HAK  
         + abcs_mur.left_h1*(hz[left_boundary +1][j] - hz_last_time[left_boundary ][j]) 1wW4bg 5  
         + 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]); Jo9c|\4  
   //x=h，右边界 mWuhXY^Q  
   hz[right_boundary][j] = hz_last_time[right_boundary-1][j] !%QbE[Kl>  
         + abcs_mur.left_h1*(hz[right_boundary-1][j] - hz_last_time[right_boundary][j]) BXr._y, cr  
         + 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]); '>>@I~<\  
  } n;k B_i*l  
     for(i = left_boundary+1;i <= right_boundary-1;i++) ]CgZt'h{  
  { :U-yO 9!j  
   //y=1，下边界 ;h#CT#R2  
   hz[bottom_boundary] = hz_last_time[bottom_boundary+1] ~AR0,lak  
         + abcs_mur.left_h1*(hz[bottom_boundary+1] - hz_last_time[bottom_boundary]) \p\rPfY{>  
         + 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]); [Rzn>  
   //y=h，上边界 uU1q?|4  
   hz[ top_boundary ] = hz_last_time[  top_boundary-1 ] Dm`gzGl  
         + abcs_mur.left_h1*(hz[  top_boundary-1 ] - hz_last_time[ top_boundary  ]) "~T06!F45  
         + 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 ]); >k(AQW5?  
  } Job&qW9W`  
  //角点的处理 N~B'gJJDx  
  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]); rm
,h\  
  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]); s_76)7  
  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]); hYh~[Kr^@^  
  hz[right_boundary][top_boundary   ] = hz_last_time[right_boundary-1][ top_boundary-1  ] + abcs_mur.C2*(hz[right_boun .. cX5tx]  
@4jPaqa(  

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

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

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

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

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