problema sobre el programa IDA FDTD 1D

[sai1ormoon]

hola:
Trato de código IDA (dirección alterna implícito) FDTD,
así que el código 1D ADIFDTD en el espacio libre en primer lugar.

Mi programa de simulación de un pulso gaussiano se propagan en el espacio libre y el periódico se utiliza BC BC (la ola que propage aparecerá a la derecha de la izquierda).

En mi programa, la onda de pulso puede propagar hacia el exterior desde el punto de origen, y el paso de tiempo (dt) más grande que puede original (alrededor de 1,25 veces).

Sin embargo, una peculiar error ocurre, la ola de "decadencia" cuando se propaga.

Realmente no sé que
está mal.¿Alguien sabe qué tipo de error que puede resultar?

muchas gracias.

este es mi código

Código:%

% IDA FDTD 1D prueba

% Una dimensión ola ola propagar hacia x y polarizar en dirección z (

Ez Hy%)

% Por sai1ormoon

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

% El algoritmo es adi

%

% T> t 1 / 2

% Actualización Ez (implícito)

% Actualización Hy (explícito)

% Fijado periódico BC

%

% T 1 / 2 -> t 1

% Actualización Ez (explícito)

% Añadir fuente

% Actualización Hy (explícito)

% Fijado periódico BC

%

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%claro todos;

% redes

Nx = 100;

Ez% Hy campo y campo

Ez ceros = (1, Nx 1);

Hy ceros = (1, Nx);% parámetro físico

cc = 2.99792458e8;

muz0 = 4 * pi * 1e-7;

epsz0 = 1/cc/cc/muz0;% lambda fuente parámetro y la frecuencia

lambda = 500e-9;

freq = cc / lambda;

wavenum = frecuencia * 2 * pi / cc;% tamaño de la red

dx = lambda/20.0;

courant condición%

courant = 0,8;

dx = dt / cc / courant;%% para Ez implícita en el espacio libre

AEZ ceros = (1, Nx 1);

Bez ceros = (1, Nx 1);

Cez ceros = (1, Nx 1);

G = ceros (1, Nx 1);

P = ceros (1, Nx 1);para (p = 1: Nx 1)

AEZ (p) =-dt 2/4/epsz0/muz0/dx ^ ^ 2;

Bez (p) = 1 2 * dt 2/4/epsz0/muz0/dx ^ ^ 2;

Cez (p) =-dt 2/4/epsz0/muz0/dx ^ ^ 2;

finalFDTD% bucle principal

for (t = 1:600)% actualización Ez-T> t 1 / 2 (implícito)

para (e = 2: Nx)

P (E) = Ez (e) dt/2/epsz0/dx * (Hy (e)-Hy (e-1));

final

% A * resolver Ez = P es un tri-diagonal de la matriz, es Ez campo, el solucionador

% algoritmo de 'Recetas numéricos para resolver un tridiagonal matriz'

Bez apuesta = (2);

Ez (2) = P (2) / apuesta;

for (i = 3: Nx)

G (i) = Cez (i-1) / apuesta;

apuesta = Bez (i)-AEZ (i) * G (i);

Ez (i) = (P (i)-AEZ (i) * Ez (i-1)) / apuesta;

final

for (i = Nx-1: -1:2)

Ez (i) = Ez (i)-G (i 1) * Ez (i 1);

final% actualización de Hy-T> t 1 / 2 (explícita)

para (e = 1: (NX-1))

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

final% simplificar frontera utilizando periódico condición límite

Ez (1) = Ez (Nx);

Hy (Nx) = Hy (1);Ez actualización% t 1 / 2 -> t 1para (e = 2: Nx)

Ez (e) = Ez (e) dt/2/epsz0/dx * (Hy (e)-Hy (e-1));

finalpunto fuente%

Ez (50) Ez = (50) exp (- ((t-50) / 13) ^ 2);

Hy% actualizar t 1 / 2 -> t 1

para (e = 1: (NX-1))

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

final% periódico frontera

Ez (1) = Ez (Nx);

Hy (Nx) = Hy (1);parcela%

si mod (t, 5) == 0;rtime = num2str (ronda (t * dt/1.0e-15));subparcela (2,1,1), parcela (Ez), eje ([0 (Nx) -1 1]);

título ([ 'tiempo =', rtime, "fs"]);

ylabel ( 'EZ');subparcela (2,1,2);

parcela (Hy);

axis ([0-Nx 3.0e-3 3.0e-3]);

título ([ 'tiempo =', rtime, 'ns']);

xlabel ( 'x (m)');

ylabel ( 'IA');

pausa (0,1);

finalfinal

 

[albede]

Supongo que la condición de frontera periódico no es correcta.
Ez (1) = Ez (Nx) -> Ez (1) = Ez (Nx 1);