проблеми> про Аді FDTD 1D програми

[sai1ormoon]

привіт:
Я намагаюся код ADI (змінних напрямків неявних) FDTD,
так що я код 1D ADIFDTD у вільному просторі в першу чергу.

Моя програма моделювання гауссовского імпульсу поширюються у вільному просторі і до н.е. використовували періодичні BC (хвилі, propage направо виявиться з лівого боку).

У своїй програмі, пульсової хвилі можуть поширюватися назовні від точкового джерела, і крок за часом (DT) може більше, ніж оригінальні (близько 1,25 рази).

Проте своєрідний помилка відбувається, хвиля "розпаду", коли вона поширюється.

Я дійсно не знаю, де не так.Хто-небудь знає які помилки можуть привести його?

Велике спасибі.

це мій код

Код:%

% 1D FDTD зв'язок випробування

% Одному вимірі хвилі хвилі поширюються в напрямку х і поляризації в напрямку осі г (

Ez% Hy)

% До sai1ormoon

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%

% Аді алгоритм

%

% Т-> т 1 / 2

% Оновлення Ez (неявна)

Hy% оновлення (явним)

Безліч періодичних% до н.е.

%

% T 1 / 2 -> т 1

% Оновлення Ez (прив'язки)

% Додати джерело

Hy% оновлення (явним)

Безліч періодичних% до н.е.

%

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%очистити всі;

% Сітках

Nx = 100;

% Ez поля і поля Hy

Ег = нулі (1, Nx 1);

Hy = нулі (1, Nx);% Фізичний параметр

CC = 2.99792458e8;

muz0 = 4 * пі * 1e-7;

epsz0 = 1/cc/cc/muz0;% Джерел лямбда-параметра і частоти

лямбда = 500E-9;

Частотний = сс / лямбда;

хвильове число = частота * 2 * пі / см;розмір сітки%

ДХ = lambda/20.0;

% Куранта умова

Куранта = 0,8;

DT = ах / CC / Куранта;%% Для ЕГ неявної у вільному просторі

AEZ нулів = (1, Nx 1);

Без нулів = (1, Nx 1);

CEZ = нулі (1, Nx 1);

G = нулі (1, Nx 1);

P = нулі (1, Nx 1);(Р = 1: Nx 1)

AEZ (р) =- DT ^ ^ 2/4/epsz0/muz0/dx 2;

Без (р) = 1 2 * DT ^ ^ 2/4/epsz0/muz0/dx 2;

CEZ (р) =- DT ^ ^ 2/4/epsz0/muz0/dx 2;

кінець% FDTD основний цикл

для (T = 1:600)Оновлення Ez% т-> т 1 / 2 (передбачувані)

для (E = 2: Nx)

P (е) = Е, (е) * dt/2/epsz0/dx (Hy (е)-Hy (е-1));

кінець

% * Вирішити Ez = P є трьох-діагональні матриці, ЕГ полі, вирішувач

Алгоритм% від Чисельні рецепти, щоб вирішити трехдіагональной матриці '

Без ставка = (2);

Ez (2) = P (2) / ставки;

для (I = 3: Nx)

G (я) = CEZ (я-1) / ставки;

Без ставка = (я)-AEZ (я) * G (я);

Ez (я) = (P (я)-AEZ (я) * Ez (я-1)) / ставки;

кінець

для (я = Nx-1: -1:2)

Ez (я) = Е, (я)-G (я 1) * Ez (я 1);

кінець% Оновлення Hy т-> т 1 / 2 (явним)

для (е = 1: (Nx-1))

Hy (е) = Hy (е) * dt/2/muz0/dx (Ez (е 1)-Ez (е));

кінець% Спрощення кордону за допомогою періодичної граничне умова

Ег (1) = Ez (Nx);

Hy (Nx) = Hy (1);% Оновлення Ez т 1 / 2 -> т 1для (E = 2: Nx)

Ez (е) = Е, (е) * dt/2/epsz0/dx (Hy (е)-Hy (е-1));

кінецьточки% джерело

Ez (50) = Е, (50) ехр (- ((T-50) / 13) ^ 2);

% Оновлення Hy т 1 / 2 -> т 1

для (е = 1: (Nx-1))

Hy (е) = Hy (е) * dt/2/muz0/dx (Ez (е 1)-Ez (е));

кінець% Періодичної кордоні

Ег (1) = Ez (Nx);

Hy (Nx) = Hy (1);% Дільниці

якщо мода (т, 5) == 0;rtime = num2str (кругла (т * dt/1.0e-15));сюжетні (2,1,1), земельну ділянку (Ez), осі ([0 (Nx) -1 1]);

Назва (час ['=', rtime, 'фс']);

ylabel ('EZ');сюжетні (2,1,2);

ділянка (Hy);

осі ([0 NX-3.0E-3 3.0E-3]);

Назва (час ['=', rtime, 'нс']);

xlabel ('х (у метрах)');

ylabel ('HY');

пауза (0,1);

кінецькінець

 

[albede]

Я думаю, періодичні граничне умова не є правильним.
Ег (1) = Ez (Nx) -> Е, (1) = Ez (Nx 1);