各位高手们,帮帮小弟,小弟快疯了~!

young20123 · 发表于 2006-5-26 10:02

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

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

x

下面是小弟编的程,出问题,解决不了了,哪位好心的高手帮忙看一下. function FDTD_debug
%constants c_0 = 3E8; % Speed of light in free space Nz = 100; % Number of cells in z-direction Nt = 200; % Number of time steps N = Nz; % Global number of cells
E = zeros(Nz,1); % E component E1 = zeros(Nz,1); Es = zeros(Nz,Nz); % Es is the output for surf function H = zeros(Nz,1); % H component H1 = zeros(Nz,1); J = zeros(Nz,1); J1 = zeros(Nz,1); K = zeros(Nz,1); K1 = zeros(Nz,1); 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 = 10e9; % Excitation frequency f_ref = f_0; % Reference frequency 
omega_0 = 2.0*pi*f_0; % Excitation circular frequency omega_ref = omega_0; omega_p=10e8;
tao=10e8;
lambda_ref = c_ref/f_ref; % Reference wavelength
Dx_ref = lambda_ref/20; % Reference cells width Dz = Dx_ref; nDz = Dz/Dx_ref;
Z = Nz*Dz; r = 1; % Normalized time step Dt_ref = r*Dx_ref/c_ref; % Reference time step Dt = Dt_ref;
% Source position Nz_Source = 0.5*Nz; N_Source = Nz_Source;
for n = 1:Nt-1 t = Dt_ref*r*(n-1); % Actual time %Source function series Source_type = 1; switch Source_type case 1 % modified source function ncycles = 2; if t < ncycles*2.0*pi/(omega_0) pulse = -0.5*( 1.0 - cos(omega_0*t/ncycles) ) * sin(omega_0*t); else pulse = 0; end case 2 % sigle cos source function if t < 2.0*pi/(omega_0) pulse = 8*c_0^2*Dt_ref^2*mu_ref*omega_0*cos(omega_0*t); else pulse = 0; end case 3 % Gaussian pulse if t < Dt_ref*r*50 pulse = -40*c_0*(t-t*25/(n-1))*exp(-(t-t*25/(n-1))^2/2/(50/2.3548)^2/(t/(n-1))^2); 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; J1(N_Source) = J1(N_Source) - r*pulse; J1 = J; K1(N_Source) = K1(N_Source) - r*pulse; K1 = K; b=2:Nz-1 E(b)=E1(b)-(Dt/(eps_0*Dz))*(H(b+0.5)-H(b-0.5)+0.5*(J1(b+0.5)+J1(b-0.5))*Dz); H(b+0.5)=H1(b+0.5)-(Dt/(mu_0*Dz))*(E1(b+1)-E1(b)+K1(b+0.5)*Dz); K(b+0.5)=((1-0.5*tao*Dt)/(1+0.5*tao*Dt))*K1(b+0.5)+(mu_0*omega_p^2*Dt/(1+0.5*tao*Dt))*H(b+0.5); J(b+0.5)=((1-0.5*tao*Dt)/(1+0.5*tao*Dt))*J1(b+0.5)+(0.5*eps_0*omega_p^2*Dt/(1+0.5*tao*Dt))*(E(b)+E(b+1)); %Boundary E(1) = 0; E(Nz) = 0; % Dirichlet %E(1) = E(2);E(Nz) = E(Nz-1); % Neumann %E(Nz) = E2(Nz-1) + (r-1)/(r+1)*(E(Nz-1)-E2(Nz)); % Mur ABC @ right side z = Nz %E(1) = E2(2) + (r-1)/(r+1)*(E(2)-E2(1)); % Mur ABC @ left side z = 0 %display %choice*********************************************** display = 3; % choice: '1' = line plot; '2' = imagesc; '3' = surf switch display case 1 i = 1:Nz; plot(i, E(i)); axis([0 Nz -4 4]); case 2 A = rot90(E); imagesc(A); case 3 i = 1:Nz; for j = 1:Nz Es(i,j) = E(i); end Es = rot90(Es); j = 1:Nz; surf(i,j,Es); axis([0 Nz 0 Nz -4 4]) otherwise A = rot90(E); imagesc(A); end pause(0.005); end ************************************************************************************************************* ??? Subscript indices must either be real positive integers or logicals.
Error in ==> E:\matlab6.5\work\FDTD1Dme.m On line 91 ==> E(b)=E1(b)-(Dt/(eps_0*Dz))*(E(b+0.5)-E(b-0.5)+0.5*(E1(b+0.5)+E1(b-0.5))*Dz);[em06]

young20123 · 发表于 2006-5-26 10:04

我就不明白那4个式子为什么出问题,哪位高手帮忙解答一下

yangzj · 发表于 2006-5-26 10:55

H(b+0.5)这里面怎么能出现小数呢,只能是正整数

young20123 · 发表于 2006-5-26 16:58

这个是利用FDTD(时域有限差分法)得到的式子,0.5是半个步长

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

各位高手们,帮帮小弟,小弟快疯了~!

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

浏览过的版块