二维TE波FDTD方法引入高斯源作图,显示收敛阶matlab实现
时间: 2023-06-25 21:06:53 浏览: 87
二维TE波的FDTD方法可以通过引入高斯源来模拟电磁波在介质中的传播。以下是基于MATLAB的二维TE波FDTD方法引入高斯源的实现步骤:
1. 定义模拟区域和网格
```matlab
% 定义模拟区域大小和网格大小
L = 1e-2; % 模拟区域边长,单位:m
dx = 1e-3; % 网格大小,单位:m
dy = dx;
Nx = round(L/dx); % x方向网格数
Ny = round(L/dy); % y方向网格数
```
2. 定义模拟参数
```matlab
% 定义模拟的频率和时间步长
f = 1e9; % 频率,单位:Hz
lambda = 3e8/f; % 波长,单位:m
dt = dx/(2*3e8); % 时间步长,单位:s
nsteps = 1000; % 模拟步数
```
3. 初始化场值和介质参数
```matlab
% 初始化场值和介质参数
ex = zeros(Nx,Ny);
ey = zeros(Nx,Ny);
ez = zeros(Nx,Ny);
hx = zeros(Nx,Ny);
hy = zeros(Nx,Ny);
hz = zeros(Nx,Ny);
% 定义介质的电导率、介电常数和磁导率
sigma = ones(Nx,Ny)*0.01; % 电导率,单位:S/m
epsilon = ones(Nx,Ny)*8.85e-12; % 介电常数,单位:F/m
mu = ones(Nx,Ny)*4*pi*1e-7; % 磁导率,单位:H/m
```
4. 定义高斯源
```matlab
% 定义高斯源
cx = Nx/2; % x方向中心位置
cy = Ny/2; % y方向中心位置
s = 5; % 高斯函数的标准差
t0 = 50; % 高斯函数的峰值时间
gauss = zeros(Nx,Ny);
for i = 1:Nx
for j = 1:Ny
gauss(i,j) = exp(-((i-cx)^2 + (j-cy)^2)/(2*s^2))*sin(2*pi*f*(n*dt-t0));
end
end
```
5. 更新场值
```matlab
% 更新场值
for n = 1:nsteps
% 更新ex、ey、ez
for i = 2:Nx-1
for j = 2:Ny-1
ex(i,j) = ex(i,j) + dt/(epsilon(i,j)*dx)*(hz(i,j)-hz(i,j-1));
ey(i,j) = ey(i,j) - dt/(epsilon(i,j)*dy)*(hz(i,j)-hz(i-1,j));
ez(i,j) = ez(i,j) + dt/epsilon(i,j)*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1));
end
end
% 更新hx、hy、hz
for i = 1:Nx-1
for j = 1:Ny-1
hx(i,j) = hx(i,j) - dt/(mu(i,j)*dy)*(ez(i+1,j)-ez(i,j));
hy(i,j) = hy(i,j) + dt/(mu(i,j)*dx)*(ez(i,j+1)-ez(i,j));
hz(i,j) = hz(i,j) + dt/mu(i,j)*(ex(i,j+1)-ex(i,j)-ey(i+1,j)+ey(i,j));
end
end
% 引入高斯源
ez(Nx/2,Ny/2) = ez(Nx/2,Ny/2) + gauss(Nx/2,Ny/2);
% 显示场值
imagesc(ez');
xlabel('x');
ylabel('y');
title(['Step ' num2str(n) ' of ' num2str(nsteps)]);
colorbar;
drawnow;
end
```
6. 计算收敛阶
```matlab
% 计算收敛阶
N = 2:2:Nx;
error = zeros(size(N));
for k = 1:length(N)
Nx = N(k);
Ny = Nx;
dx = L/Nx;
dy = dx;
dt = dx/(2*3e8);
ex = zeros(Nx,Ny);
ey = zeros(Nx,Ny);
ez = zeros(Nx,Ny);
hx = zeros(Nx,Ny);
hy = zeros(Nx,Ny);
hz = zeros(Nx,Ny);
sigma = ones(Nx,Ny)*0.01;
epsilon = ones(Nx,Ny)*8.85e-12;
mu = ones(Nx,Ny)*4*pi*1e-7;
gauss = zeros(Nx,Ny);
cx = Nx/2;
cy = Ny/2;
s = 5;
t0 = 50;
for i = 1:Nx
for j = 1:Ny
gauss(i,j) = exp(-((i-cx)^2 + (j-cy)^2)/(2*s^2))*sin(2*pi*f*(n*dt-t0));
end
end
for n = 1:nsteps
for i = 2:Nx-1
for j = 2:Ny-1
ex(i,j) = ex(i,j) + dt/(epsilon(i,j)*dx)*(hz(i,j)-hz(i,j-1));
ey(i,j) = ey(i,j) - dt/(epsilon(i,j)*dy)*(hz(i,j)-hz(i-1,j));
ez(i,j) = ez(i,j) + dt/epsilon(i,j)*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1));
end
end
for i = 1:Nx-1
for j = 1:Ny-1
hx(i,j) = hx(i,j) - dt/(mu(i,j)*dy)*(ez(i+1,j)-ez(i,j));
hy(i,j) = hy(i,j) + dt/(mu(i,j)*dx)*(ez(i,j+1)-ez(i,j));
hz(i,j) = hz(i,j) + dt/mu(i,j)*(ex(i,j+1)-ex(i,j)-ey(i+1,j)+ey(i,j));
end
end
ez(Nx/2,Ny/2) = ez(Nx/2,Ny/2) + gauss(Nx/2,Ny/2);
end
error(k) = norm(ez(:)-ez_old(:));
ez_old = ez;
end
order = diff(log(error))./diff(log(N));
```
以上是二维TE波FDTD方法引入高斯源作图,显示收敛阶的MATLAB实现。