二维TE波CNDG-FDTD方法引入高斯源,求数值解与解析解误差,显示收敛阶并画图matlab实现
时间: 2023-06-20 10:05:35 浏览: 154
CNDG-FDTD方法是一种结合了传统的有限差分时域(FDTD)方法和Curl-Nodal Discontinuous Galerkin(CNDG)方法的求解电磁波问题的方法。在二维情况下,TE波表示电场只在z方向有分量,磁场只在x、y方向有分量。在本文中,我们将引入高斯源,并利用Matlab实现求解数值解与解析解误差、显示收敛阶并画图。
假设我们要求解的是一个长宽分别为Lx和Ly的矩形区域内的TE波。首先,我们需要将区域网格化,将整个区域划分为N_x \times N_y个小矩形,每个小矩形的长宽分别为\Delta x和\Delta y。我们假设在该区域内引入一个高斯源,其表达式为:
$$
S(x,y,t)=S_0 exp\left(-\frac{(x-x_0)^2+(y-y_0)^2}{2\sigma^2}\right)cos(\omega t)
$$
其中,S_0是源的振幅,(x_0,y_0)是源的位置,\sigma是源的宽度,\omega是源的频率。
接下来,我们可以利用CNDG-FDTD方法求解TE波方程:
$$
\nabla \times \nabla \times \mathbf{E} - \mu \epsilon \frac{\partial^2 \mathbf{E}}{\partial t^2} = \mu \mathbf{J}_s
$$
其中,\mathbf{E}表示电场,\mu和\epsilon分别表示介质的磁导率和电介质常数,\mathbf{J}_s表示高斯源产生的电流密度。我们可以将该方程离散化为:
$$
\mathbf{C} \frac{\partial \mathbf{E}}{\partial t} + \mathbf{G} \mathbf{E} = \mathbf{S}
$$
其中,\mathbf{C}和\mathbf{G}是系数矩阵,\mathbf{S}是源项矩阵。我们可以利用FDTD方法求解\mathbf{C}和\mathbf{G},利用CNDG方法求解\mathbf{S}。
最后,我们可以利用解析解求出TE波的精确解,与数值解做比较,计算误差,并绘制误差随网格大小变化的图像,以显示收敛阶。具体实现可参考以下Matlab代码:
```
% 定义常数
Lx = 1; % 区域长度
Ly = 1; % 区域宽度
Nx = 100; % x方向网格数
Ny = 100; % y方向网格数
dx = Lx / Nx; % x方向网格间距
dy = Ly / Ny; % y方向网格间距
T = 1e-8; % 模拟总时间
dt = dx / sqrt(2); % 时间间隔
Nt = T / dt; % 时间步数
x = linspace(0, Lx, Nx+1); % x方向网格
y = linspace(0, Ly, Ny+1); % y方向网格
omega = 2 * pi * 1e9; % 频率
mu = pi * 4e-7; % 磁导率
epsilon = 8.85e-12; % 电介质常数
sigma = dx / 4; % 高斯源宽度
S0 = 1; % 高斯源振幅
x0 = Lx / 2; % 高斯源x方向位置
y0 = Ly / 2; % 高斯源y方向位置
% 初始化系数矩阵和源项矩阵
C = zeros(Nx*Ny, Nx*Ny);
G = zeros(Nx*Ny, Nx*Ny);
S = zeros(Nx*Ny, 1);
% 构造系数矩阵和源项矩阵
for i = 1:Nx
for j = 1:Ny
% 计算节点编号
n = (j-1)*Nx + i;
% 构造C矩阵
if i == 1 || i == Nx || j == 1 || j == Ny
C(n, n) = 1;
else
C(n, (j-2)*Nx+i) = 1/(dy*dy);
C(n, j*Nx+i) = 1/(dy*dy);
C(n, (j-1)*Nx+i-1) = 1/(dx*dx);
C(n, (j-1)*Nx+i+1) = 1/(dx*dx);
C(n, n) = -2/(dx*dx) - 2/(dy*dy);
end
% 构造G矩阵
G(n, n) = mu * epsilon / (dt*dt);
% 构造S矩阵
xj = (i-0.5) * dx;
yj = (j-0.5) * dy;
r = sqrt((xj-x0)^2 + (yj-y0)^2);
S(n) = S0 * exp(-r^2/(2*sigma*sigma)) * cos(omega*dt*Nt);
end
end
% 求解电场
E = zeros(Nx*Ny, Nt);
E(:, 1) = 0;
E(:, 2) = dt * inv(C) * (S - G * E(:, 1));
for n = 3:Nt
E(:, n) = inv(C) * (2*S - G * E(:, n-1) - E(:, n-2));
end
% 计算解析解
E_analytical = zeros(Nx, Ny, Nt);
for n = 1:Nt
for i = 1:Nx
for j = 1:Ny
xj = (i-0.5) * dx;
yj = (j-0.5) * dy;
r = sqrt((xj-x0)^2 + (yj-y0)^2);
E_analytical(i, j, n) = S0 * exp(-r^2/(2*sigma*sigma)) * cos(omega*dt*n);
end
end
end
% 计算误差
error = zeros(1, Nt);
for n = 1:Nt
error(n) = norm(E(:, n) - reshape(E_analytical, Nx*Ny, Nt)(:, n));
end
% 计算收敛阶
order = log(error(2:end)./error(1:end-1))/log(sqrt(dx*dx+dy*dy));
order = [order(1) order];
% 绘制误差随时间变化图像
figure;
plot(dt*(0:Nt-1), error);
xlabel('时间(秒)');
ylabel('误差');
% 绘制误差随网格大小变化图像
figure;
loglog(sqrt(dx*dx+dy*dy).^(1:Nt), error);
xlabel('网格大小');
ylabel('误差');
% 绘制收敛阶随时间变化图像
figure;
plot(dt*(1:Nt), order);
xlabel('时间(秒)');
ylabel('收敛阶');
```
需要注意的是,本文中的TE波方程是在时域中求解的,因此需要将频率离散化为角频率,并在源项中加入cos(\omega t)项。此外,我们还需要注意在计算解析解时的节点位置和电场方向。
阅读全文