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

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

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

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

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

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

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

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

<P>这个是利用FDTD(时域有限差分法)得到的式子,0.5是半个步长</P>
页: [1]
查看完整版本: 各位高手们,帮帮小弟,小弟快疯了~!