一维电磁波的时域有限差分法(FTDT)

Question

我尝试编写一个代码来求解以下耦合偏微分电磁波方程：

​

​

代码使用使用Yee算法的时域有限差分，该算法可在以下两个在线文档中阅读：

http://www.eecs.wsu.edu/~schneidj/ufdtd/ufdtd.pdf

http://www.eecs.wsu.edu/~schneidj/ufdtd/chap3.pdf

我希望我的左手边界处的信号源是一个正弦波，参数如下：

Ex(0,t) = E0 sin(2πft)      for 0 ≤ t ≤ 1/f

代码使用每个循环更新波的电场和磁场属性。

我的初始代码如下：

%FDTD Yee algorithm to solve coupled EM wave equations
clear
clc


G=50;           %Specify size of the grid
f=10^3;         %choose initial frequency of wave in hertz
e=1;            %specify permitivity and permeability (normalised condition)
u=1;            
Nt=150;         %time steps
E0=1;           %Electric Field initial amplitude

%Specify step sizes using corruant condition
c=3*10^8;
dx=0.01;
dt=dx/2*c;

%make constant terms
c1=-dt/(dx*e);
c2=-dt/(dx*u);

%create vgector place holders
Ex=zeros(1,G);
Hy=zeros(1,G);

%create updating loop
M=moviein(Nt);
for t=1:Nt

    % Spatial Ex

    for k=2:G-1
        Ex(k)=Ex(k)+c1*(Hy(k)-Hy(k-1));
    end
    Ex(G)=0; %PEC boundary condition
    %E Source at LHS boundary 
    Ex(1)=E0*sin(2*pi*f*t);
    %Spatial Hy
    for n=1:G-1
        Hy(n)=Hy(n)+c2*(Ex(n)-Ex(n+1));
    end
    Hy(G)=0; %PMC boundary condition

plot(Ex);
M(:,t) = getframe;
end
movie(M,1);

基本上，我想创建一个更新的电影，显示正弦波传播到右侧边界，编码为完美的电导体；因此反射波，然后传播回左侧边界，编码为完美绝缘体；吸收波。

我所遇到的问题如下：

1)我不确定如何正确地实现我想要的源码。它看起来不是纯正弦的。

2)我编码的波开始传播，但在模拟的大部分时间里它很快就消失了。我不知道为什么会发生这种情况

3)我不太了解运行电影模拟，并且随着解决方案的运行，情节振荡。我怎么才能阻止它呢？

Answer 1

您的波衰减是因为差分方程没有正确实现；相反：

Ex(k)=Ex(k)+c1*(Hy(k)-Hy(k-1));

你应该使用

Ex1(k)=Ex(k)+c1*(Hy(k)-Hy(k-1));

而不是：

Hy(n)=Hy(n)+c2*(Ex(n)-Ex(n+1));

您应该使用：

Hy1(n)=Hy(n)+c2*(Ex(n)-Ex(n+1));

并且，在循环结束时，更新当前的"dataframe"：

Hy = Hy1;
Ey = Ey1;

(您还应注意边界条件)。

如果您希望固定的绘图框在数据更改时不会更改，请首先创建一个可在其中绘制的axis，其中定义了xmin/max和ymin/max，请参见http://www.mathworks.com/help/matlab/ref/axis.html

Answer 2

您应该将Courant数设置为更接近于1，例如0.995。因此，delta_t =0.995*delta/c。另外，假设delta_x使用公制单位，则e和u应该使用公制单位。我不知道具体使用的编码语言，但在c或c++中，不需要中间变量Ey1等。此外，为了准确起见，每个波长至少应该有10个样本(最好是60个)。因此波长= 60*delta_x，因此频率大致等于10的9次方。另外，我认为正弦源应该是E0 * sin(2* pi *f*t* delta_t)。您需要调整您的常量，然后重试