FDTD продолжает переполняться квази границ и некоторые другие вопросы

Question

У меня есть назначение FDTD, оно не очень хорошо. Всякий раз, когда я добавляю источник и позволяю моделированию немного поработать (Примечание: если вы тестируете, вы должны использовать отладчик и пошагово выполнять приращения цикла, в противном случае произойдет сильное переполнение).

В большинстве симуляций, которые я видел, люди разбивают Ex и Ey на отдельные части, я складываю их в одну большую матрицу E и использую смещение индекса (видно в Ep (u, v + 1) и E (u). + 1, v) рассчитать Ey и Ex независимо).

Источники, которые я использовал для справки: http://fdtd.wikispaces.com/ а также http://www.mathworks.com/matlabcentral/fileexchange/21000-tiny-fdtd-v1-0 Это акустика, но она хорошо работает.

    close all;


    %% Some user modifiable parameters
    mu0 = pi*4E-7; % pH/µm
    e0 =   8.854187E-12; % Picofarads/micron
    c = 1/sqrt(mu0*e0);

    cellsizeX = 100;      % Size of Yee-Cell in microns
    cellsizeY = 100;      % Size of Yee-Cell in microns
    numX = 100; % Number of cells in X direction
    numY = 100; % Number of cells in Y direction
    lambda = 700*10^-9;
    dx = lambda/20;
    dy = lambda/20;
    dt =  (c*sqrt(dx^-2+dy^-2))^-1;

    t0 = 100;       %index time of gaussian pulse peak
    width = 10;     %peakyness of gaussian pulse

    %% Initialise the H and E array
    H = zeros(2*numX, 2*numY);       
    Hp = zeros(2*numX, 2*numY);
    E = zeros(2*numX+1,2*numY+1);    
    Ep = zeros(2*numX+1,2*numY+1);    
    Etemp = zeros(2*numX+1,2*numY+1);
    Htemp = zeros(2*numX, 2*numY);
    P = zeros(2*numX+1,2*numY+1);
    Pp = zeros(2*numX+1,2*numY+1);

    % Scaling factors for H and E fields

    CEx = dt/(dx*mu0);
    CEy = dt/(dy*mu0);
    CHx = dt/(dy*e0);
    CHy = dt/(dx*e0);
    x = 2:2:2*numX; %2*numX-2;

    %Initialize Permibilities
    Perm = ones(2*numX+1,2*numY+1);


    %% FDTD loop
for n = 1:1500;

    if(n < 200)
        E(numX-10:2:numX+10,numY+1) =  1E-19*sin(2*pi*428E12*n*dt);    %exp(-.5 * ((n-t0)/width).^2); %insert hard source
    end;

    for u = 2:2:2*numX-1
        for v = 2:2:2*numY-1
            Hp(u,v) = H(u,v) + (dt/mu0)*( -(E(u+1,v) - E(u-1,v))/dx) + ( E(u,v+1) - E(u,v-1)/dy ); % Solving for Hplus
            Ep(u,v+1) = E(u,v+1) + (dt/(dy*e0))*(Hp(u,v+2) - Hp(u, v)); % Solving for Ex plus
            Ep(u+1, v) = E(u+1, v) - (dt/(dx*e0))*(Hp(u+2, v) - Hp(u, v)); % Solving for Ey plus

        end
    end;


    % Dirichlet Boundary Conditions
    Ep(1,:) = 0;
    Ep(:,1) = 0;
    Ep(2*numX+1,:) = 0;
    Ep(:,2*numX+1) = 0;

    % Plotting
    surf(Ep); shading interp; lighting phong; colormap hot; axis off; zlim([0 1]);
    set(gcf,'Color', [0 0 0], 'Number', 'on', 'Name', sprintf('FDTD Project, step = %i', n));
    %title(sprintf('Time = %.2f microsec',n*dt*1e12),'Color',[1 0 0],'FontSize', 22); 
    drawnow;

    E = Ep;
    H = Hp;



end;

end;

Answer 1

2 вещи:

1 - это Ep(:,2*numX+1,:) = 0; правильно? (т.е. это должно быть Ep(:,2*numX+1) = 0;?)

2 - в самом конце, я думаю, вы можете изменить код на

if (max(max(Ep)) > 1e-3)
      E = Ep/max(max(Ep));
else E = Ep;
end
if (max(max(Hp)) > 1e-3)
    H = Hp/max(max(Hp));
else H = Hp;
end

как max (Ep) даст вам матрицу 1x100. 1e-3 для защиты, вы можете понизить его до любого значения по вашему выбору.