S
sai1ormoon
Guest
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
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