二维TE波CNDG-FDTD方法引入高斯源,求数值解误差,显示收敛阶并画图matlab实现
时间: 2023-06-19 17:09:51 浏览: 52
二维TE波CNDG-FDTD方法是一种求解电磁波传播问题的数值方法,而高斯源是一种常用的激励源,可以用来模拟电磁波在空间中的传播。在这里,我们介绍如何引入高斯源,并求解数值解误差、显示收敛阶以及画图。
首先,我们需要了解二维TE波CNDG-FDTD方法的基本原理。在这种方法中,我们将空间划分成一个个小网格,然后通过有限差分的方式来近似求解电磁场的变化。具体来说,我们用中心差分法来近似求解电场和磁场在时间和空间上的变化,并在每个网格点上引入一个时间步长dt和空间步长dx。这样,我们就可以通过迭代计算,得到电磁场在空间中的分布情况。
在引入高斯源的时候,我们需要在计算过程中添加一个激励项,用来描述高斯源的影响。具体来说,我们可以在计算区域的中心位置引入一个高斯源,其具体形式为:
$E_x(x,y,t) = E_{x0}e^{-\frac{(x-x_0)^2+(y-y_0)^2}{2\sigma^2}}sin(\omega t)$
其中,$E_{x0}$表示高斯源的幅值,$(x_0,y_0)$表示高斯源的位置,$\sigma$表示高斯源的宽度,$\omega$表示高斯源的频率。在实现中,我们可以将高斯源看作一个时间上的脉冲信号,通过每个时间步长来逐渐加入高斯源的影响。
为了求解数值解误差和显示收敛阶,我们需要将计算区域划分成不同的网格大小,并分别计算出每个网格大小下的数值解和精确解之间的误差。然后,我们可以通过比较不同网格大小下的误差来计算出收敛阶。具体来说,如果我们将网格大小按照$d$的倍数逐渐减小,那么误差应该会按照$d^p$的方式逐渐减小,其中$p$就是收敛阶。
最后,我们可以使用Matlab来实现上述过程,并通过绘制图表来显示结果。具体实现过程如下:
1. 首先,我们需要定义计算区域的大小和网格大小,以及高斯源的参数。
```matlab
% 定义计算区域大小和网格大小
nx = 100; % x方向网格数
ny = 100; % y方向网格数
dx = 0.01; % x方向网格大小
dy = 0.01; % y方向网格大小
dt = dx / (2 * 10^8); % 时间步长
% 定义高斯源参数
Ex0 = 1; % 高斯源幅值
sigma = 0.02; % 高斯源宽度
x0 = 0.5; % 高斯源x方向位置
y0 = 0.5; % 高斯源y方向位置
omega = 2 * pi * 10^9; % 高斯源频率
```
2. 接着,我们需要初始化电场和磁场的数值,并按照时间步长迭代计算电磁场的变化。
```matlab
% 初始化电场和磁场的数值
Ex = zeros(nx, ny);
Ey = zeros(nx, ny);
Hx = zeros(nx, ny);
Hy = zeros(nx, ny);
% 迭代计算电磁场的变化
for n = 1:1000
% 计算电场在时间和空间上的变化
for i = 2:nx-1
for j = 2:ny-1
Ex(i,j) = Ex(i,j) + (dt/(dx*dy)) * (Hy(i,j)-Hy(i-1,j)-Hx(i,j)+Hx(i,j-1));
Ey(i,j) = Ey(i,j) - (dt/(dx*dy)) * (Hy(i,j)-Hy(i-1,j)-Hx(i,j)+Hx(i,j-1));
end
end
% 在中心位置引入高斯源
Ex(round(nx/2), round(ny/2)) = Ex(round(nx/2), round(ny/2)) + Ex0 * exp(-((dx*(round(nx/2)-x0))^2+(dy*(round(ny/2)-y0))^2)/(2*sigma^2)) * sin(omega*n*dt);
% 计算磁场在时间和空间上的变化
for i = 2:nx-1
for j = 2:ny-1
Hx(i,j) = Hx(i,j) - (dt/(dx*dy)) * (Ey(i,j+1)-Ey(i,j)+Ex(i,j)-Ex(i,j+1));
Hy(i,j) = Hy(i,j) + (dt/(dx*dy)) * (Ex(i+1,j)-Ex(i,j)-Ey(i,j)+Ey(i,j-1));
end
end
end
```
3. 接着,我们可以计算数值解和精确解之间的误差,并显示收敛阶。在这里,我们将网格大小按照2的倍数逐渐减小。
```matlab
% 定义不同网格大小
dxs = [0.02, 0.01, 0.005, 0.0025];
% 计算不同网格大小下的误差
for i = 1:length(dxs)
% 定义当前网格大小
dx = dxs(i);
dy = dx;
dt = dx / (2 * 10^8);
% 初始化电场和磁场的数值
Ex = zeros(nx, ny);
Ey = zeros(nx, ny);
Hx = zeros(nx, ny);
Hy = zeros(nx, ny);
% 迭代计算电磁场的变化
for n = 1:1000
for i = 2:nx-1
for j = 2:ny-1
Ex(i,j) = Ex(i,j) + (dt/(dx*dy)) * (Hy(i,j)-Hy(i-1,j)-Hx(i,j)+Hx(i,j-1));
Ey(i,j) = Ey(i,j) - (dt/(dx*dy)) * (Hy(i,j)-Hy(i-1,j)-Hx(i,j)+Hx(i,j-1));
end
end
% 在中心位置引入高斯源
Ex(round(nx/2), round(ny/2)) = Ex(round(nx/2), round(ny/2)) + Ex0 * exp(-((dx*(round(nx/2)-x0))^2+(dy*(round(ny/2)-y0))^2)/(2*sigma^2)) * sin(omega*n*dt);
for i = 2:nx-1
for j = 2:ny-1
Hx(i,j) = Hx(i,j) - (dt/(dx*dy)) * (Ey(i,j+1)-Ey(i,j)+Ex(i,j)-Ex(i,j+1));
Hy(i,j) = Hy(i,j) + (dt/(dx*dy)) * (Ex(i+1,j)-Ex(i,j)-Ey(i,j)+Ey(i,j-1));
end
end
end
% 计算数值解和精确解之间的误差
error(i) = norm(Ex - exact_solution) / norm(exact_solution);
% 显示收敛阶
if i > 1
order(i-1) = log(error(i-1)/error(i)) / log(2);
end
end
```
4. 最后,我们可以使用Matlab中的plot函数来绘制误差随网格大小变化的图表。
```matlab
% 绘制误差随网格大小变化的图表
figure(1);
loglog(dxs, error, '-o');
xlabel('dx');
ylabel('error');
title('Error vs. Grid Size');
grid on;
% 绘制收敛阶随网格大小变化的图表
figure(2);
plot(dxs(2:end), order, '-o');
xlabel('dx');
ylabel('order');
title('Order vs. Grid Size');
grid on;
```
这样,我们就可以得到误差随网格大小变化和收敛阶随网格大小变化的图表了。需要注意的是,为了计算精确解,我们需要知道电磁场在计算区域边界处的边界条件。在这里,我们假设边界条件为电场和磁场都为零。
相关推荐
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![](https://csdnimg.cn/download_wenku/file_type_ask_c1.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)
![zip](https://img-home.csdnimg.cn/images/20210720083736.png)