请高手们帮忙指点迷津,小弟在这里先谢谢~!

young20123 · 发表于 2006-5-18 11:04

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

您需要登录才可以下载或查看，没有账号？我要加入

x

这是我编的一个matlab程序,是用一维FDTD方法对正弦电磁波穿过DNG媒质的平行板时的传播过程模拟,出了问题,但是不知道在哪,希望您能帮指点一下错误. 还有,能不能提供一个PML边界吸收条件的matlab程序,万分谢谢了~! %constants c_0 = 3E8; % Speed of light in free space Nz = 5000; % Number of cells in z-direction Nt = 2000; % Number of time steps N = Nz; % Global number of cells
E = zeros(Nz,1); % E component E2 = zeros(Nz,1); E1 = zeros(Nz,1); Es = zeros(Nz,Nz); % Es is the output for surf function %H = zeros(Nz,1); % H component pulse = 0; % Pulse
%to be set mu_0 = 4.0*pi*1.0e-7; % Permeability of free space eps_0 = 1.0/(c_0*c_0*mu_0); % Permittivty of free space
%c_ref = c_0; % Reference velocity %eps_ref = eps_0; %mu_ref = mu_0;
f_0 = 3*10e10; % Excitation frequency %f_ref = f_0; % Reference frequency 
omega_0 = 2.0*pi*f_0; % Excitation circular frequency %omega_ref = omega_0;
Dt=9.5*10e-14; Dz=3.0*10e-5;
r = 1; 
T_p=1/f_0; omega_p=10e8;
tao=10e8;
m=5; n=10;
% Source position Nz_Source = 0.12*Nz; N_Source = Nz_Source;
for n1 = 1:Nt-1 t=Dt*r*(n1-1); x_on=1-(m*T_p-t)/(m*Tp); x_off=(t-(m+n)*T_p)/(m*T_p); g_on=10*x_on^3-15*x_on^4+6*x_on^5; g_off=1-(10*x_off^3-15*x_off^4+6*x_off^5); %Source function series Source_type = 1; switch Source_type case 1 if t<m*T_p; pulse =g_on*sin(omega_0*t); else pulse = 0; end case 2 if t<=(m+n)*T_p; pulse=sin(omega_0*t); else pulse = 0; end case 3 if t<=(m+n+m)*T_p; pulse=g_off*sin(omega_0*t); else pulse = 0; end case 4 if t>(m+n+m)*T_p; pulse = 0; else pulse = 0; end end E1(N_Source) = E1(N_Source) - r*pulse; E1 = E; H1(N_Source) = H1(N_Source) - r*pulse; H1 = H; J=diff(E,t); J1=diff(E1,t); K=diff(H,t); K1=diff(H1,t); i=2:Nz-1 E(i)=E1(i)-(Dt/(eps_0*Dz))*(H(i+0.5)-H(i-0.5)+0.5*(J1(i+0.5)+J1(i-0.5))*Dz); H(i+0.5)=H1(i+0.5)-(Dt/(mu_0*Dz))*(E1(i+1)-E1(i)+K1(i+0.5)*Dz); K(i+0.5)=((1-0.5*tao*Dt)/(1+0.5*tao*Dt))*K1(i+0.5)+(mu_0*omega_p^2*Dt/(1+0.5*tao*Dt))*H(i+0.5); J(i+0.5)=((1-0.5*tao*Dt)/(1+0.5*tao*Dt))*J1(i+0.5)+(0.5*eps_0*omega_p^2*Dt/(1+0.5*tao*Dt))*(E(i)+E(i+1)); %Boundary E(1) = 0; E(Nz) = 0; % Dirichlet fprintf(fid,'%10.4f %10.4f\n',E(i),t); end plot(E(i),t); hold on

[此贴子已经被作者于2006-5-18 15:31:27编辑过]

suffer · 发表于 2006-5-18 11:12

两个个非常明显的错误
x_on=1-(m*T_p-t)/(m*Tp); 其中Tp应该是T_p吧
H1(N_Source) = H1(N_Source) - r*pulse; H1没有初始值
这种错误自己看一下错误信息应该就能看出来吧
另外给你一个简易，问问题的时候最好不要指名问谁，否则有些人看到了就算知道也未必回答你

suffer · 发表于 2006-5-18 11:12

其他错误自己在检查检查吧

lxq · 发表于 2006-5-18 11:13

<DIV class=quote>以下是引用suffer在2006-5-18 11:12:32的发言： 
两个个非常明显的错误
x_on=1-(m*T_p-t)/(m*Tp); 其中Tp应该是T_p吧
H1(N_Source) = H1(N_Source) - r*pulse; H1没有初始值
这种错误自己看一下错误信息应该就能看出来吧
另外给你一个简易，问问题的时候最好不要指名问谁，否则有些人看到了就算知道也未必回答你</DIV>
 suffer 说的很对啊~~~

young20123 · 发表于 2006-5-18 15:27

谢谢,谢谢,小弟知道了

young20123 · 发表于 2006-5-18 15:30

运行了,就说是 Error: Missing operator, comma, or semicolon. 所以求大家帮忙找错误了

ericlin · 发表于 2006-5-18 15:36

缺了括号或标点符号什么的，用一用字典也行啊。一般有错的地方它都会告诉你再第几行第几列的！！

[此贴子已经被作者于2006-5-18 15:39:39编辑过]

happy · 发表于 2006-5-18 15:38

<DIV class=quote>以下是引用ericlin在2006-5-18 15:36:49的发言： 缺了括号或标点符号什么的，用一用字典也行啊。</DIV>
 呵呵[em17]

young20123 · 发表于 2006-5-18 15:42

弱弱的问一下,字典怎么用啊,小弟刚刚接触matlab

ericlin · 发表于 2006-5-18 15:50

打开金山词霸，开启取词功能，把鼠标放到你不认识的英文字上。 可能你少了逗号或分号吧。 你初学MATLAB就编了这么大一串已经很厉害了，继续努力啊。

[此贴子已经被作者于2006-5-18 15:56:27编辑过]

young20123 · 发表于 2006-5-18 15:52

谢谢大家了

账号		自动登录	找回密码
密码			我要加入

请高手们帮忙指点迷津,小弟在这里先谢谢~!

马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。

回复：（young20123）请Happy老师帮忙指点迷津,谢谢~!...

回复：（young20123）请Happy老师帮忙指点迷津,谢谢~!...

回复：（suffer）回复：（young20123）请Happy老师帮忙...

回复：（ericlin）缺了括号或标点符号什么的，用一用字...