计算电磁学之FDTD算法的MATLAB语言实现
矩形谐振腔电磁场的FDTD分析和Matlab仿真
矩形谐振腔电磁场的FDTD分析和Matlab仿真摘要：目前，电磁场的时域计算方法越来越引人注目。
这种方法已经广泛应用到各种电磁问题的分析之中。
而将Matlab作为一种仿真工具，用于时域有限差分法，可以简化编程，使研究者重心放在FDTD本身上，而不必在编程上花费过多的时间。
本课题通过用FDTD方法计算矩形谐振腔电磁场分布，并用Matlab 进行仿真。
关键词：时域有限差分法，Matlab仿真，矩形谐振腔1.引言时域有限差分法（Finite-Dfference Time-Domain Method）是求解电磁问题的一种数值技术，是在1966年由K.S.Yee第一次提出的。
FDTD法直接将有限差分式代替麦克斯韦时域场旋度方程中的微分式，得到关于场分量的有限差分式，用具有相同电参量的空间网格去模拟被研究体，选取合适的场始值和计算空间的边界条件，可以得到包括时间变量的麦克斯韦方程的四维数值解，通过傅里叶变换可求得三维空间的频域解。
时域有限差分法突出的优点是所需的存储量及计算时间与N成正比，使得很多复杂的电磁场计算问题成为可能，用时域有限差分法容易模拟各种复杂的结构，使得用其他方法不能解决的问题有了新的处理方法。
本文主要讨论如何用Matlab语言来编写FDTD的吸收边界条件以及编程时应注意的问题。
2时域有限差分法的基本理论2.1 时域有限差分法的简介1966年K.S.Yee首次提出了一种电磁场数值计算的新方法——时域有限差分（Finite-Dfference Time-Domain Method）方法。
对电磁场E、H分量在时间和空间上采取交替抽样的离散方式，每一个E（或H）场分量四周有四个H（或E）场分量环绕，应用这种离散方式将含时间变量的麦克斯韦旋度方程转化为一组差分方程，并在时间轴上逐步推进地求解空间电磁场。
Yee提出的这种抽样方式后来被称为Yee元胞。
FDTD方法是求解麦克斯韦方程的直接时域方法。
在计算中将空间某一样本点的电场（或磁场）与周围格点的磁场（或电场）直接相关联，且介质参数已赋值给空间每一个元胞，因此这一方法可以处理复杂形状目标和非均匀介质物体的电磁散射、辐射等问题。
基于MATLAB的FDTD算法编程
FD TD ar ithm etic programm ing ba sed on M ATLAB language
ZHAO J ia (Com pu ter Engineering D epartm en t, Guangx i U n iversity
of T echno logy, L iuzhou 545006, Ch ina) Abstract: T he basic p rincip le of FD TD w as in troduced and the 22D TM m ode FD TD arithm etic exp ression w as derived. T he p rogramm ing m ethod based on M A TLAB w as illu strated w ith exam p le. Key words: M A TLAB; FD TD; Yee A rithm etic; P rogramm ing
自由空间的光速。∃x、∃y 和 ∃ t 的选取必须使它们的关系满足 Cournan t 稳定性条件。尽管时间和空间步长越
小, FD TD 运算的精度会越高, 但太细的网格会使计算机运算时间太长, 同时, 在网格分辨率 (即波长与空间
步长之比) 小到一定程度时, 数值误差基本就不变了。对于采用M A TLAB 编程, 需要特别注意时间步长和空
211 算例
考察区域为正方形的二维自由空间, 区域的边
界为理想电导体, 为了激励外向辐射波, 设在区域的
中心有电场分量 E z 随时间按高斯变化, 求解该区
域内电场随时间的变化, 场域如图 2 所示。
212 M ATLAB 编程方法及结果
基于M A TLAB 对本算例 FD TD 算法进行编程 求 解可归纳为以下几个步骤: (1) 正确建立二维
FDTD算法与Matlab仿真在电磁场与电磁波教学中的应用
FDTD算法与Matlab仿真在电磁场与电磁波教学中的应用黎杨;程莉
【期刊名称】《高师理科学刊》
【年(卷),期】2022(42)11
【摘要】麦克斯韦方程组是电磁场与电磁波课程中最核心的一个知识点,物理量的抽象性、数学运算的复杂程度和预言产生电磁波等这些方程组的性质学生学习起来特别困难.从麦克斯韦方程组的微分形式出发,使用时域有限差分法(FDTD)推演时变电场和时变磁场在Yee胞中的空间和时间维度更新过程,借助于Matlab编程实现,可视化演示了电磁波的产生和传播过程,使学生更深入、更直观地理解“由麦克斯韦方程组预言电磁波的产生”这一重要性质.
【总页数】5页(P86-90)
【作者】黎杨;程莉
【作者单位】武汉工程大学电气信息学院
【正文语种】中文
【中图分类】O441;G642.0
【相关文献】
1.FDTD算法在"电磁场与电磁波"课程教学中的应用
2.基于Matlab语言实现电磁场中的FDTD算法编程
3.应用MATLAB设计电磁场与电磁波模拟仿真实验
4.Matlab软件在电磁场与电磁波可视化教学中的应用
5.Matlab在“电磁场与电磁波”课程可视化教学中的探索——以平面电磁波在介质交界面处的反射与透射分析为例
因版权原因，仅展示原文概要，查看原文内容请购买。
MATLAB编程在FDTD算法中的应用
MATLAB编程在FDTD算法中的应用摘要：文章介绍了时域有限差分（FDTD）法的基本原理，推导了二维TM 波的FDTD表达式。
给出了MATLAB语言编程的步骤和应注意的问题，并结合算例阐述了基于MATLAB仿真的基本方法,最后得出用MATLAB语言对FDTD 算法仿真的几点结论。
关键词：MATLAB；FDTD；编程时域有限差分（FDTD）法是在21世纪60年代由K.S.Yee首先提出并用于求解电磁场散射问题，其主要思路是在空间轴和时间轴上对场量进行离散，并用中心差分代替偏微分，这就将麦克斯韦方程组转化为了差分方程，通过在时间轴和空间轴上采取蛙跳法（leapfrog）逐步推进地求解，最终求得在一定边值与初值条件下的空间场解。
随着计算机技术的发展，FDTD的应用越来越多，对于FDTD算法的编程求解，最常用的程序语言有VC和FORTRUN，而MATLAB 作为一种可视化效果好的软件，在FDTD计算中可视化程度较高，并具有能显示动态场效果的特点。
文章采用FDTD法对高压开关柜内超高频电磁波的辐射和传播特性进行仿真，仿真中将对二维TM电磁波进行FDTD表达式的推导，并结合FDTD算法边界条件的特点，用MATLAB语言进行编程。
1二维TM电磁波FDTD算法1.1算式推导在自由空间中，对于二维TM电磁波，，MAXWELL的两个旋度方程可以分解为：＝-0 （1）＝0（2）０=-（3）构造二维TM波FDTD cell如图1所示：按照FDTD元胞对以上三式中的偏导用中心差商代替，分别可得：=-0 （4）=0 （5）０=-（6）由于方程中出现了半网格和半时间步，为了便于编程，可以将上面差分式改为如下FDTD算式：Hx（i,j,k+1）=Hx（i,j,k）-[EZ（i,j+1,n）-EZ（i,j,n）]（7）Hy（i,j,n+1）=Hy（i,j,n）+[EZ（i+1,j,n）-EZ（i,j,n）]（8）Ez（i,j,n+1）=Ez（i,j,n）+[-] （9）1.2数值稳定性的条件FDTD方法是以一组有限差分方程来替代MAXWELL旋度方程来进行数值计算的方法，在执行形如FDTD算法时，随着时间步长的增长，如何保证该算法的稳定性是一个重要问题。
matlab实现一维fdtd仿真
麦克斯韦方程组为
(2—1—1)
（2—1—2）
各项同性线性介质中的本构关系为
（2—1—3）
在直角坐标系中，（2—1—1）、（2—1—2）式写为
（2—1—4）
以及
（2—1—5）下面我们考虑（2—1—4）、（2—1—5）式的FDTD差分离散。令 代表 或 直角坐标系中某一分量，在时间和空间域中的离散取以下符号表示：
FDTD在电磁研究的多个领域获得广泛应用，其中有：辐射天线的分析，微波器件和导行波结构的研究，散射和雷达截面计算，周期机构分析，电子封装和电磁兼容分析，核电磁脉冲的传播和散射，微光学元器件中光的传播和颜射特性。
二、FDTD模拟电磁波传播的原理
麦克斯韦方程组是支配宏观电磁现象的一组基本方程。这组方程既可以写成微分形式，又可以写成积分形式。FDTD方法是由微分形式的麦克斯韦旋度方程出发进行差分离散从而得到一组时域推进公式。首先给出麦克斯韦旋度方程及其在直角坐标系中的FDTD离散形式，由三维情况到一维情况，即一维电磁波传播的理论基础。
c = 1/sqrt(mu_0*eps_0); %光速
%定义问题空间和参数
domain_size = 1; %空间问题的长度
dx = 1e-3; %单元网格大小
dt = 3e-12; %步长时间的持续
number_of_time_steps = 2000; %迭代次数
nx = round(domain_size/dx); %在问题空间的单元数目
3.2参数设置
程序设计了电流屏在z轴上传播，产生的电场和磁场，此电流屏放置在问题空间的中心，与两个相互平行的完善导体平面平行，整个区间充满空气。激励源在两板之间传播和反射。
无论是简单目标还是复杂目标，在进行FDTD离散时网格尺寸的确定，除了受计算资源的限制不可能取得很小外，还需要考虑以下几个因素：
矩形谐振腔电磁场的FDTD分析和Matlab仿真
矩形谐振腔电磁场的FDTD分析和Matlab仿真摘要：目前，电磁场的时域计算方法越来越引人注目。
这种方法已经广泛应用到各种电磁问题的分析之中。
而将Matlab作为一种仿真工具，用于时域有限差分法，可以简化编程，使研究者重心放在FDTD本身上，而不必在编程上花费过多的时间。
本课题通过用FDTD方法计算矩形谐振腔电磁场分布，并用Matlab 进行仿真。
关键词：时域有限差分法，Matlab仿真，矩形谐振腔1.引言时域有限差分法（Finite-Dfference Time-Domain Method）是求解电磁问题的一种数值技术，是在1966年由K.S.Yee第一次提出的。
FDTD法直接将有限差分式代替麦克斯韦时域场旋度方程中的微分式，得到关于场分量的有限差分式，用具有相同电参量的空间网格去模拟被研究体，选取合适的场始值和计算空间的边界条件，可以得到包括时间变量的麦克斯韦方程的四维数值解，通过傅里叶变换可求得三维空间的频域解。
时域有限差分法突出的优点是所需的存储量及计算时间与N成正比，使得很多复杂的电磁场计算问题成为可能，用时域有限差分法容易模拟各种复杂的结构，使得用其他方法不能解决的问题有了新的处理方法。
本文主要讨论如何用Matlab语言来编写FDTD的吸收边界条件以及编程时应注意的问题。
2时域有限差分法的基本理论2.1 时域有限差分法的简介1966年K.S.Yee首次提出了一种电磁场数值计算的新方法——时域有限差分（Finite-Dfference Time-Domain Method）方法。
对电磁场E、H分量在时间和空间上采取交替抽样的离散方式，每一个E（或H）场分量四周有四个H（或E）场分量环绕，应用这种离散方式将含时间变量的麦克斯韦旋度方程转化为一组差分方程，并在时间轴上逐步推进地求解空间电磁场。
Yee提出的这种抽样方式后来被称为Yee元胞。
FDTD方法是求解麦克斯韦方程的直接时域方法。
在计算中将空间某一样本点的电场（或磁场）与周围格点的磁场（或电场）直接相关联，且介质参数已赋值给空间每一个元胞，因此这一方法可以处理复杂形状目标和非均匀介质物体的电磁散射、辐射等问题。
matlab模拟的电磁学时域有限差分法 pdf
matlab模拟的电磁学时域有限差分法 pdf电磁学时域有限差分法（FDTD）是一种基于数值模拟的电磁场计算方法，它使用有限差分来近似微分方程。
该方法广泛用于电磁学、电波传播、微波技术、光学等领域，以求解电磁场分布和场的辐射、散射等问题。
而在这个领域中，MATLAB是非常流行的工具之一。
本文将围绕“MATLAB模拟的电磁学时域有限差分法”这一主题，从以下几个方面进行阐述：1.时域有限差分法的基础概念在FDTD方法中，将时域中的Maxwell方程组转化为差分形式，使得可以在计算机上进行数值解法。
通过在空间和时间上的离散，可以得到电磁场在时域内的各种分布，进而求得特定情况下的电磁场变化。
2.MATLAB中的FDTD仿真在MATLAB中，我们可以使用PDE工具箱中的电磁学模块来实现FDTD仿真。
通过选择适当的几何形状和边界条件，可以利用该工具箱演示电磁场的传输、反射、折射、透射等现象。
同时，MATLAB中还提供了不同的场分量计算和可视化工具，以便用户可以更好地理解电磁场分布。
3.MATLAB代码实现以下是一些MATLAB代码示例，展示了FDTD模拟的基础实现方法。
代码中的示例模拟了平面波在一个矩形和圆形障碍物上的传播情况。
% 1. Square obstaclegridSize = 200; % Grid sizemaxTime = 600; % Maximum time (in steps)imp0 = 377.0; % Impedance of free spacecourantNumber = 0.5; % Courant numbercdtds = ones(gridSize,gridSize); % Courant number in space% (not variable in this example)Ez = zeros(gridSize, gridSize); % Define EzHy = zeros(gridSize, gridSize); % Define Hy% Simulation loopfor n = 1:maxTime% Update magnetic fieldHy(:,1:end-1) = Hy(:,1:end-1) + ...(Ez(:,2:end) - Ez(:,1:end-1)) .*cdtds(:,1:end-1) / imp0;% Update electric fieldEz(2:end-1,2:end-1) = Ez(2:end-1,2:end-1) + ...(Hy(2:end-1,2:end-1) - Hy(1:end-2,2:end-1)) .* cdtds(2:end-1,2:end-1) .* imp0;end% 2. Circular obstacleradius = 50;xAxis = [-100:99];[X,Y] = meshgrid(xAxis);obstacle = sqrt((X-50).^2 + (Y).^2) < radius;gridSize = length(xAxis); % Grid sizemaxTime = 500; % Maximum time (in steps)imp0 = 377.0; % Impedance of free space courantNumber = 0.5; % Courant numbercdtds = ones(gridSize,gridSize); % Courant number in space% (not variable in this example)Ez = zeros(gridSize, gridSize); % Define EzHy = zeros(gridSize, gridSize); % Define Hy% Simulation loopfor n = 1:maxTime% Update magnetic fieldHy(:,1:end-1) = Hy(:,1:end-1) + ...(Ez(:,2:end) - Ez(:,1:end-1)) .*cdtds(:,1:end-1) / imp0;% Update electric field, with obstacleEz(2:end-1,2:end-1) = Ez(2:end-1,2:end-1) + ...(Hy(2:end-1,2:end-1) - Hy(1:end-2,2:end-1)) .* cdtds(2:end-1,2:end-1) .* imp0;Ez(obstacle) = 0;end以上代码仅供参考，不同条件下的模拟需要适当修改，以便获得特定的模拟结果。
matlab模拟的电磁学时域有限差分法
matlab模拟的电磁学时域有限差分法时域有限差分法（FDTD）是一种计算电磁波传播及散射的数值模拟方法。
它是基于麦克斯韦方程组进行仿真的一种方法，而且从计算电磁波传播的实质上来看，FDTD方法是一种求解时域麦克斯韦方程的有限差分方法。
在FDTD方法中，我们将区域空间离散化，并定义电场、磁场等量的格点值。
然后，根据麦克斯韦方程组的时域形式，在各个时刻进行场量的更新。
FDTD方法在实践应用中具有计算时间和空间复杂度低，且适用于复杂的结构和非线性介质等特点，所以在电磁学数值仿真中应用广泛。
我们可以用MATLAB来进行FDTD的电磁学仿真，下面详细介绍MATLAB的使用步骤：1. 建立空间离散化格点在仿真开始前，需要先根据空间大小和仿真目的来建立离散化格点。
对于一个一维的结构，我们可以用以下代码来建立：x = linspace(0,1,N); %建立离散化空间格点Ex = zeros(1,N); %电场，长度为N的全0数组Hy = zeros(1,N); %磁场，长度为N的全0数组其中N为获取离散化格点数量的参数，x为离散化空间格点，Ex和Hy为电场和磁场。
2. 定义电场和磁场边界条件在进行仿真时，需要了解仿真的边界情况并将其定义成特殊的边界条件。
例如，仿真空间内可能存在各种元件、环境等，这些都会对电场和磁场的性质产生影响。
所以，我们需要用特殊边界条件来约束仿真空间内电场和磁场的行为。
在FDTD中，通常采用数值反射边界条件（DNG Boundary）来进行仿真。
例如，在这个边界条件下，在仿真空间内部设置经典的电场边界条件：场强等于零；并在仿真空间外部添加一层基质，该基质的介电常数和磁导率均为负值，并且在该基质中场的强度和方向均反向。
相当于在仿真空间外设置一个虚拟折射界面，能够将场边界反射。
我们设定如下代码：M = 20; % 反射界面层数Ex_low_M1 = 0; %反射界面边界条件Ex_high_M1 = 0; %反射界面边界条件for i = 1:MEx_low_M2(i) = Ex_high_M2(i-1); %反转反射界面内的电场贡献Ex_high_M2(i) = Ex_low_M2(i-1); %反转反射界面内的电场贡献end3. 计算电场的场值FDTD仿真中最核心的内容就是判断时刻要计算的电场场值。
计算电磁学之FDTD算法的MATLAB语言实现要点
South China Normal University课程设计实验报告课程名称：计算电磁学指导老师：专业班级： 2014级电路与系统姓名：学号：FDTD算法的MATLAB语言实现摘要：时域有限差分(FDTD)算法是K.S.Yee于1966年提出的直接对麦克斯韦方程作差分处理,用来解决电磁脉冲在电磁介质中传播和反射问题的算法。
其基本思想是:FDTD计算域空间节点采用Yee元胞的方法,同时电场和磁场节点空间与时间上都采用交错抽样；把整个计算域划分成包括散射体的总场区以及只有反射波的散射场区,这两个区域是以连接边界相连接,最外边是采用特殊的吸收边界,同时在这两个边界之间有个输出边界,用于近、远场转换;在连接边界上采用连接边界条件加入入射波,从而使得入射波限制在总场区域；在吸收边界上采用吸收边界条件,尽量消除反射波在吸收边界上的非物理性反射波。
本文主要结合FDTD算法边界条件特点,在特定的参数设置下，用MATLAB语言进行编程，在二维自由空间TEz网格中，实现脉冲平面波。
关键词：FDTD；MATLAB；算法1 绪论1.1 课程设计背景与意义20世纪60年代以来，随着计算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（FEM）、矩量法(MM)和单矩法等；属于时域技术方面的时域有限差分法(FDTD)、传输线矩阵法(TLM)和时域积分方程法等。
其中FDTD是一种已经获得广泛应用并且有很大发展前景的时域数值计算方法。
时域有限差分（FDTD）方法于1966年由K.S.Yee提出并迅速发展，且获得广泛应用。
K.S.Yee用后来被称作Yee氏网格的空间离散方式，把含时间变量的Maxwell旋度方程转化为差分方程，并成功地模拟了电磁脉冲与理想导体作用的时域响应。
但是由于当时理论的不成熟和计算机软硬件条件的限制，该方法并未得到相应的发展。
20世纪80年代中期以后，随着上述两个条件限制的逐步解除，FDTD便凭借其特有的优势得以迅速发展。
fdtd solutions matlab代码
以下是关于FDTD（Finite Difference Time Domain）方法的Matlab代码。
1. FDTD方法简介FDTD方法是一种数值分析电磁场问题的方法，最早应用于求解Maxwell方程组。
该方法的基本思想是将时间和空间分割为离散网格，并利用差分法求解Maxwell方程组。
FDTD方法广泛应用于天线设计、电磁兼容性分析、光学器件仿真等领域。
2. FDTD方法的Matlab代码以下是一个简单的一维FDTD方法的Matlab代码示例：```matlab定义常数c = 3e8; 光速dx = 0.01; 空间步长dt = dx/(2*c); 时间步长初始化场量Ez = zeros(1,1000); 电场Hy = zeros(1,1000); 磁场模拟时间步进for n = 1:1000更新磁场for i = 1:999Hy(i) = Hy(i) + (Ez(i+1) - Ez(i))/(c*dx);end更新电场for i = 2:1000Ez(i) = Ez(i) + (Hy(i) - Hy(i-1))*(c*dt/dx);endend绘制结果figure;plot(Ez);xlabel('空间步长');ylabel('电场强度');title('FDTD方法仿真结果');```3. 代码解释- 代码首先定义了常数c（光速）、空间步长dx和时间步长dt。
- 然后初始化了电场Ez和磁场Hy的空间网格。
- 在时间步进的循环中，首先更新磁场Hy，然后更新电场Ez。
- 最后绘制了Ez随空间步长的变化图。
这是一个非常简单的一维FDTD方法的Matlab代码示例，实际应用中可能会有更复杂的三维情况，需要考虑更多的边界条件和介质性质。
4. FDTD方法的应用FDTD方法在天线设计、电磁兼容性分析、光学器件仿真等领域有着广泛的应用。
通过编写相应的Matlab代码，可以对复杂的电磁场问题进行数值分析，得到定量的仿真结果，为工程设计和科学研究提供重要的参考。
计算电磁学之FDTD算法的MATLAB语言实现
South China Normal University课程设计实验报告课程名称：计算电磁学指导老师：专业班级： 2014级电路与系统姓名：学号：FDTD算法的MATLAB语言实现摘要：时域有限差分(FDTD)算法是K.S.Yee于1966年提出的直接对麦克斯韦方程作差分处理,用来解决电磁脉冲在电磁介质中传播和反射问题的算法。
其基本思想是:FDTD计算域空间节点采用Yee元胞的方法,同时电场和磁场节点空间与时间上都采用交错抽样；把整个计算域划分成包括散射体的总场区以及只有反射波的散射场区,这两个区域是以连接边界相连接,最外边是采用特殊的吸收边界,同时在这两个边界之间有个输出边界,用于近、远场转换;在连接边界上采用连接边界条件加入入射波,从而使得入射波限制在总场区域；在吸收边界上采用吸收边界条件,尽量消除反射波在吸收边界上的非物理性反射波。
本文主要结合FDTD算法边界条件特点,在特定的参数设置下，用MATLAB语言进行编程，在二维自由空间TEz网格中，实现脉冲平面波。
关键词：FDTD；MATLAB；算法1 绪论1.1 课程设计背景与意义20世纪60年代以来，随着计算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（FEM）、矩量法(MM)和单矩法等；属于时域技术方面的时域有限差分法(FDTD)、传输线矩阵法(TLM)和时域积分方程法等。
其中FDTD是一种已经获得广泛应用并且有很大发展前景的时域数值计算方法。
时域有限差分（FDTD）方法于1966年由K.S.Y ee提出并迅速发展，且获得广泛应用。
K.S.Y ee用后来被称作Y ee氏网格的空间离散方式，把含时间变量的Maxwell旋度方程转化为差分方程，并成功地模拟了电磁脉冲与理想导体作用的时域响应。
但是由于当时理论的不成熟和计算机软硬件条件的限制，该方法并未得到相应的发展。
20世纪80年代中期以后，随着上述两个条件限制的逐步解除，FDTD便凭借其特有的优势得以迅速发展。
基于MATLAB的FDTD算法编程
（） ３
方程中出现了半个网格和半个时间步 ， 为便于编程 ， 可将上面差分式改写为如下 Ｆ Ｔ Ｄ Ｄ算式 ：
Ｈ ， ＋ １ － Ｈ ， ） （， ｋ ） （， 愚 一
ｔ ，
￣Ｅ ｌ 一 ＿］ ｙ ｊ ， Ｅ，） Ｅ， （， （＋ ） 『
［ ， 一Ｅ （，，） Ｅ （＋１ ，） ］
，
８７ １ｚ
ａ ￡ Ｈ
’
Ｂ Ｅ． ａ Ｈ
。 。 一
８ Ｈ
ａ ｖ
３ 一 一 ｙ
构 造 二维 Ｔ 波 ＹＥ ｅ 如 图 ｌ 示 ： Ｍ Ｅｃｌ ｌ 所
按 Ｙ Ｅ元胞对上式偏导用中心差商代替， Ｅ 可得 ：
生
Ａｙ
一
生１ — １ １ — １ ． ＋
△￡
— — — — — — — ～ ： ：： — — — — — — — — — — — — — 一
（） ２ Leabharlann ｏ —Ｅ ，＋ ） ， 丢 （， 吉一 （ ） Ｅ ， 一
— — — — — — — — 一
日（ ｌ， 丢一 十 ）Ｈ（ ＋， 专 日（『＋ ） ， ＋ ）Ｈ（ ， 号 ， ｌ＋ 一 ， ＋ ， ｊ ， ＿
时间轴和空间轴上采取蛙跳法 （ ａｆ ｇ ： 步推进地求解 ， １ ｐｒ ） ｅ ｏ ￣ 最终求出一定边值与初值条件下的空间场解 。随
着计 算机 技术 的发 展 ， 近年来 Ｆ Ｄ计 算技 术 也得到 了越 来越 多 的应用 。 ＤＴ 对于 Ｆ Ｔ 算法 的编程 求解 ， Ｄ Ｄ 最常 用 的 有 ＶＣ和 Ｆ ＴＲ ＯＲ ＵＮ， ＭＡＴ ＡＢ作 为一 种可 视化 效果 很好 的科学 计算 软件 ， Ｆ Ｔ 而 Ｌ 在 Ｄ Ｄ计算 中能充分 发 挥 编程 简单 、 可视 化 程度 高 、 显示 动态 场 能
一维等离子体FDTD的Matlab源代码(两种方法)
维等离子体FDTD的Matlab 源代码（两种方法）%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% %%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1D %%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% %%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%% %%%%%%%%%%%%%%%%%%%%%%%%% 初始化clear;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%% 系统参数TimeT=3000;% 迭代次数KE=2000;% 网格树木kc=450;% 源的位置kpstart=500;% 等离子体开始位置kpstop=1000;% 等离子体终止位置%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%% 物理参数c0=3e8;% 真空中波速zdelta=1e-9;% 网格大小dt=zdelta/(2*c0);% 时间间隔f=900e12;%Gause 脉冲的载频d=3e-15% 脉冲底座宽度t0=2.25/f;% 脉冲中心时间u0=57e12% 碰撞频率fpe=2000e12;% 等离子体频率wpe=2*pi*fpe;% 等离子体圆频率epsz=1/(4*pi*9*10A9); % 真空介电常数mu=1/(c0A2*epsz);% 磁常数ex_low_m1=0;ex_low_m2=0;ex_high_m1=0;ex_high_m2=0;aO=2*uO/dt+(2/dtF2;a1=-8/(dtf2;a2=-2*u0/dt+(2/dt)A2;b0=wpe A2+2*u0/dt+(2/dt)A2;b1=2*wpeA2-8/(dt)A2;b2=wpeA2-2*uO/dt+(2/dt)A2;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 初始化电磁场Ex=zeros(1,KE);Ex_Pre=zeros(1,KE);Hy=zeros(1,KE);Hy_Pre=zeros(1,KE);Dx=zeros(1,KE);Dx_Pre=zeros(1,KE);Sx1=zeros(1,KE);Sx2=zeros(1,KE);Sx3=zeros(1,KE);Sx=zeros(1,KE);Dx=Ex;Dx1=Ex;Dx2=Ex;Dx3=Ex;Ex1=Ex;Ex2=Ex; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%% 开始计算for T=1:TimeT%%% 保存前一时间的电磁场Ex_Pre=Ex;Hy_Pre=Hy;%%%% 中间差分计算Dxfor i=2:KEDx(i)=Dx(i)-(dt/zdelta)*(Hy(i)-Hy(i-1));end%%%%%%%% 力叭源Dx(kc)=cos(2*pi*f*T*dt)*exp(-4*pi*((T*dt-tO)/d)A2);% FDTD Main Fu nction Jobs toWorkers % %**************************************************************** *******% 3-D FDTD code with PEC boundaries%*********************************************************************** %% This MATLAB M-file implements the finite-difference time-domain %solution of Maxwell's curl equations over a three-dimensional% Cartesian space lattice comprised of uniform cubic grid cells.%% To illustrate the algorithm, an air-filled rectangular cavity% resonator is modeled. The length, width, and height ofthe% cavity are X cm (x-direction), Y cm (y-direction), and% Z cm (z-direction), respectively.%% The computational domain is truncated using PEC boundary% conditions:% ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes% ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes% ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes% These PEC boundaries form the outer lossless walls of the cavity.%% The cavity is excited by an additive current source oriented% along the z-direction. The source waveform is a differentiated% Gaussian pulse given by% J(t)=-J0*(t-t0)*exp(-(t-t0)A2 /tau A2),% where tau=50 ps. The FWHM spectral bandwidth of this zero-dc- % content pulse is approximately 7 GHz. The grid resolution% (dx = 2 mm) was chosen to provide at least 10 samples per% wavelength up through 15 GHz.%% To execute this M-file, type "fdtd3D" at the MATLAB prompt.% This M-file displays the FDTD-computed Ez fields at every other% time step, and records those frames in a movie matrix, M, which% is played at the end of the simulation using the "movie" command.%%*********************************************************************** function [Ex,Ey,Ez]=FDTD3D_Main(handles)global SimRunStop% if ~isdir('C:\MATLAB7\work\cavity\figures')% mkdir 'C:\MATLAB7\work\cavity\figures' % end%*********************************************************************** %Grid Partition%***********************************************************************p.ip = get(handles.XdirPar,'Value');p.jp = get(handles.YdirPar,'Value');p.PN = get(handles.PartNo,'Value');%***********************************************************************% Grid Dimensons%***********************************************************************ib=ie+1;jb=je+1; kb=ke+1;%*********************************************************************** %All Domains Fields Ini.%***********************************************************************Ex=zeros(ie,jb,kb);Ey=zeros(ib,je,kb);Ez=zeros(ib,jb,ke);Hx=zeros(ib,je,ke);Hy=zeros(ie,jb,ke);Hz=zeros(ie,je,kb);%*********************************************************************** %Fundamental constants%***********************************************************************=2.99792458e8; %speed of light in free spaceparam.muz=4.0*pi*1.0e-7; %permeability of free space param.epsz =1.0/(**param.muz); %permittivity of free space%***********************************************************************% Gridparameters %***********************************************************************ie =get(handles.xslider,'Value');je = %number of grid cells in x-direction %number of grid cells in y-param.is =get(handles.xsource,'Value'); param.js =get(handles.ysource,'Value'); %location of z-directed current source%location of z-directed current sourceparam.kobs = floor(ke/ 2); %Surface of observationparam.dx = get(handles.CellSize,'Value');; %space increment of cubic lattice param.dt = param.dx/(2.0*); %time stepparam.nmax = get(handles.TimeStep,'Value');; %total number of time steps%***********************************************************************% Differentiated Gaussian pulseexcitation %***********************************************************************param.rtau=get(handles.tau,'Value')*100e-12; param.tau=param.rtau/param.dt;param.ndelay=3*param.tau;param.Amp = get(handles.sourceamp,'Value')*10e11; param.srcconst=-param.dt*param.Amp;%***********************************************************************% Materialparameters %***********************************************************************param.eps= get(handles.epsilon,'Value');param.sig= get(handles.sigma,'Value');%***********************************************************************% Updating coefficients%***********************************************************************param.ca=(1.0-(param.dt*param.sig)/(2.0*param.epsz*param.eps))/(1.0+(param.dt*param.sig)/( 2.0*param.epsz* param.eps));param.cb=(param.dt/param.epsz/param.eps/param.dx)/(1.0+(param.dt*param.sig)/(2.0*param.epsz*param.eps));param.da=1.0; param.db=param.dt/param.muz/param.dx;%***********************************************************************% Calling FDTD Algorithm%*********************************************************************** ex=zeros(ib,jb,kb); ey=zeros(ib,jb,kb);ez=zeros(ib,jb,kb);hx=zeros(ib,jb,kb);hy=zeros(ib,jb,kb);hz=zeros(ib,jb,kb);[X,Y,Z] = meshgrid(1:ib,1:jb,1:kb); % Grid coordinatesPsim = zeros(param.nmax,1);Panl = zeros(param.nmax,1);if ((p.ip == 1)&(p.jp == 0))x = ceil(ie/p.PN)param.a = [1:x-1:ie-x];param.b = [x+1:x-1:ie];param.c = [1,1];param.d = [je,je];m2 = 1;for n=1:1:param.nmaxfor m1=1:1:p.PN[ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p);[hx,hy,hz] =Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p); end[Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n);field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p); Dyn_FFT pause(0.01);endelseif ((p.ip == 0)&(p.jp == 1))y = ceil(je/p.PN);param.c = [1:y-1:je-y];param.d = [y+1:y-1:je];param.a = [1,1];param.b = [ie,ie];m1 = 1;for n=1:1:param.nmaxfor m2=1:1:p.PN[ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p);[hx,hy,hz] =Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p); end[Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n);field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p); pause(0.01);endelseif ((p.ip == 1)&(p.jp == 1))x = ceil(ie/p.PN); param.a = [1:x-1:ie-x];param.b = [x+1:x-1:ie];y = ceil(je/p.PN);param.c = [1:y-1:je-y];param.d = [y+1:y-1:je];for n=1:1:param.nmaxfor m2=1:1:p.PN for m1=1:1:p.PN[ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p); [hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);endend[Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n);field_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);pause(0.01);endelseparam.a = 1;param.b = ie;param.c = 1; param.d = je;m1 = 1;m2=1;for n=1:1:param.nmax[ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,m1,m2,p); [hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,m1,m2,p);SimRunStop = get(handles.Stop_sim,'value');if SimRunStop == 1 h = warndlg('Simulation Run is Stopped byUser !','!! Warning !!'); waitfor(h);break;end [Psim(n),Panl(n)] = Cavity_Power(param,handles,ex,ey,ez,n); if n>=2Panl(n)=Panl(n) + Panl(n-1); endfield_viz(param,handles,ex,ey,ez,X,Y,Z,n,Psim,Panl,p);pause(0.01);end end%文件 2% Cavity Field Viz. %function [] = field_viz(param,handles,ex,ey,ez,X,Y ,Z,n,Psim,Panl,p)%***********************************************************************% Cross-Section initialization%***********************************************************************tview = squeeze(ez(:,:,param.kobs));sview = squeeze(ez(:,param.js,:));ax1 = handles.axes1;ax2 = handles.axes2;ax3 = handles.axes3;%*********************************************************************** % Visualize fields %***********************************************************************timestep=int2str(n);ezview=squeeze(ez(:,:,param.kobs)); exview=squeeze(ex(:,:,param.kobs));eyview=squeeze(ey(:,:,param.kobs));xmin = min(X(:)); xmax = max(X(:)); ymin = min(Y(:)); ymax = max(Y(:)); zmin = min(Z(:));daspect([2,2,1])xrange = linspace(xmin,xmax,8);yrange = linspace(ymin,ymax,8);zrange = 3:4:15;[cx cy cz] = meshgrid(xrange,yrange,zrange);% sview=squeeze(ez(:,param.js,:));rect = [-50 -35 360 210];rectp = [-50 -40 350 260];axes(ax1);axis tight set(gca,'nextplot','replacechildren');E_total = sqrt(ex.A2 + ey.A2 + ez. A2);etview=squeeze(E_total(:,:,param.kobs));sview = squeeze(E_total(:,param.js,:));surf(etview');shading interp;caxis([-1.0 1.0]);colorbar;axis image;title(['Total E-Field, time step = ',timestep],'fontname','Times Roman','fontsize',12,'FontWeight','BOLD');xlabel('i coordinate');ylabel('j coordinate');set(gca,'fontname','Times New Roman','fontsize',10);% F1 = getframe(gca,rect);% M1 = frame2im(F1);% filename =fullfile('C:\MATLAB7\work\cavity\figures',strcat('XY_ETotal',num2str(n),'.jpeg')); % imwrite(M1,filename,'jpeg')axes(ax2);axis tight set(gca,'nextplot','replacechildren');surf(sview');shading interp;caxis([-1.0 1.0]);colorbar;axis image;title(['Ez(i,j=13,k), time step = ',timestep],'fontname','TimesRoman','fontsize',12,'FontWeight','BOLD');xlabel('i coordinate');ylabel('k coordinate');set(gca,'fontname','Times New Roman','fontsize',10);% F2 = getframe(gca,rect);% M2 = frame2im(F2);% filename = New Newfullfile('C:\MATLAB7\work\cavity\figures',strcat('XZ_ETotal',num2str(n),'.jpeg'));% imwrite(M2,filename,'jpeg')%***********************************************************************% Cavity Power - Analiticexpression %***********************************************************************axes(ax3);% axis tight% set(gca,'nextplot','replacechildren');t = [1:1:param.nmax];Psim = 1e3*Psim./max(Psim);Panl = 1e3*Panl./max(Panl);semilogy(t,Psim,'b.-',t,Panl,'rx-');str(1) = {'Normalized Analytic Vs. Simulated Power'};str(2) = {'as function of time-steps'}; title(str,'fontname','Times NewRoman','fontsize',12,'FontWeight','BOLD' ); xlabel('Time Step');ylabel('Log(Power)'); legend('Simulation','Analytic','location','SouthEast');set(gca,'fontname','Times New Roman','fontsize',10);% F3 = getframe(gca,rectp);% M3 = frame2im(F3);% filename = fullfile('C:\MATLAB7\work\cavity\figures',strcat('Power',num2str(n),'.jpeg')); % imwrite(M3,filename,'jpeg')%文件 3% Source waveform functionfunction [source]=waveform(param,handles,n)ax1 = handles.axes1;type = get(handles.sourcetype,'value');amp = get(handles.sourceamp,'value')*1e4;f = get(handles.Frequency,'value')*1e9;switch typecase 2source = param.srcconst*(n-param.ndelay)*exp(-((n-param.ndelay)A2 /param.tau^));case 1source = param.srcconst*sin(param.dt*n*2*pi*f);case 3sourceparam.srcconst*exp(-((n-param.ndelay)A2 /param.tau A2))*sin(2*pi*f*(n-param.ndelay)*param.dt) Jotherwisesource = param.srcconst*(n-param.ndelay)*exp(-((n-param.ndelay)A2 /param.tauA2)); end%文件 4function [hx,hy,hz] = Hfields(param,hx,hy,hz,ex,ey,ez,ie,je,ke,ib,jb,kb,n,x,y,p)a = param.a(x);b = param.b(x);c = param.c(y);d = param.d(y);hx(a+1:b,c:d,1:ke)=...hx(a+1:b,c:d,1:ke)+...param.db*(ey(a+1:b,c:d,2:kb)-...ey(a+1:b,c:d,1:ke)+...ez(a+1:b,c:d,1:ke)-...ez(a+1:b,c+1:d+1,1:ke));hy(a:b,c+1:d,1:ke)=...hy(a:b,c+1:d,1:ke)+...param.db*(ex(a:b,c+1:d,1:ke)-...ex(a:b,c+1:d,2:kb)+...ez(a+1:b+1,c+1:d,1:ke)-...ez(a:b,c+1:d,1:ke));hz(a:b,c:d,2:ke)=...hz(a:b,c:d,2:ke)+...param.db*(ex(a:b,c+1:d+1,2:ke)-...ex(a:b,c:d,2:ke)+...ey(a:b,c:d,2:ke)-...ey(a+1:b+1,c:d,2:ke));%文件 5 function[ex,ey,ez]=Efields(param,handles,ex,ey,ez,hx,hy,hz,ie,je,ke,ib,jb,kb,n,x,y,p) a =param.a(x);b = param.b(x);c = param.c(y);d = param.d(y);if (param.is>=a)&(param.is<=b)|(param.js>=c)&(param.js<=d) source =waveform(param,handles,n);elsesource = 0;endex(a:b,c+1:d,2:ke)=...param.ca*ex(a:b,c+1:d,2:ke)+... param.cb*(hz(a:b,c+1:d,2:ke)-... hz(a:b,c:d-1,2:ke)+...hy(a:b,c+1:d,1:ke-1)-...hy(a:b,c+1:d,2:ke));ey(a+1:b,c:d,2:ke)=...param.ca*ey(a+1:b,c:d,2:ke)+... param.cb*(hx(a+1:b,c:d,2:ke)-...hx(a+1:b,c:d,1:ke-1)+...hz(a:b-1,c:d,2:ke)-...hz(a+1:b,c:d,2:ke));ez(a+1:b,c+1:d,1:ke)=...param.ca*ez(a+1:b,c+1:d,1:ke)+... param.cb*(hx(a+1:b,c:d-1,1:ke)-...hx(a+1:b,c+1:d,1:ke)+... hy(a+1:b,c+1:d,1:ke)-...hy(a:b-1,c+1:d,1:ke));ez(param.is,param.js,1:ke)=ez(param.is,param.js,1:ke)+source;function FDTDonedimensionpipei(L,d,T) %version1.0 终端匹配%FDTDonedimensionpipei(6,0.18,0.5e-9) t0=3*T;c=3e8;u=4*pi*1e-7; e=8.8541878e-12;dz=T*c/10;Nz=L/dz;Nz=fix(Nz);dt=dz/2/c;Ex=zeros(1,Nz+1);B=zeros(1,Nz+1);Hy=zeros(1,Nz);Nt=2*Nz;for n=0:Ntt=n*dt;F=exp(-(t-tO)A2./P2);Ex(1)=F;for k=1:NzHy(k)=Hy(k)+dt./u.*(Ex(k)-Ex(k+1))./dz;endfor k=1:Nz-1Ex(k+1)=Ex(k+1)+dt./e.*(Hy(k)-Hy(k+1))./dz;endEx(1)=B(2)+(c*dt-dz)./(c*dt+dz).*(Ex(2)-B(1));Ex(Nz+1)=B(Nz)+(c*dt-dz)./(c*dt+dz).*(Ex(Nz)-B(Nz+1));Vref1=d.*Ex(Nz-300);Vref2=d.*Ex(Nz-100);plot(t,Vref1,'s');hold on;plot(t,Vref2,'rx');hold on;B=Ex;endfunction FDTD_debug% Speed of light in free space% Number of cells in z-direction % Number of time steps % Global number of cells% E component% Es is the output for surf function% H component mu_0 = 4.0*pi*1.0e-7; % Permeability of free spaceeps_0 = 1.0/(c_0*c_0*mu_0); % Permittivty of free space c_ref =c_0; % Reference velocity eps_ref = eps_0;mu_ref = mu_0;f_0 = 10e9;% Excitation frequency f_ref = f_0; % Reference frequency omega_0 = 2.0*pi*f_0; omega_ref = omega_0; lambda_ref = c_ref/f_ref; Dx_ref = lambda_ref/ 20; Dz = Dx_ref;nDz = Dz/Dx_ref;Z = Nz*Dz;r = 1; % Normalized timestep Dt_ref =r*Dx_ref/c_ref; % Referencetime step Dt = Dt_ref;% Source positionNz_Source = 0.5*Nz;N_Source = Nz_Source;for n = 1:Nt-1t = Dt_ref*r*(n-1); % Actual time %Source function seriesSource_type = 1;switch Source_typecase 1 % modified source functionncycles = 2;if t < ncycles*2.0*pi/(omega_0)pulse = -0.5*( 1.0 - cos(omega_0*t/ncycles) ) * sin(omega_0*t); elsepulse = 0;endcase 2 % sigle cos source function if t < 2.0*pi/(omega_0)pulse = 8*c_0A2*Dt_refA2*mu_ref*omega_0*cos(omega_0*t);elsepulse = 0;endcase 3 % Gaussian pulse if t < Dt_ref*r*50pulse = -40*c_0*(t-t*25 /(n-1))*exp(-(t-t*25 /(n-1))人2 /2/(50/2.3548)人%constants c_0 = 3E8; Nz = 100; Nt = 200; N = Nz; E = zeros(Nz,1);E2 = zeros(Nz,1);E1 = zeros(Nz,1); Es = zeros(Nz,Nz); %H =zeros(Nz,1); pulse = 0;% Pulse %to be set % Excitation circular frequency% Reference wavelength % Reference cells width2/("(门-1))人2); elsepulse = 0; end end E2(N_Source) = E2(N_Source) - r*pulse; E1 = E2;E2 = E;m = 2:Nz-1;E(m) = r A2*(E2(m+1) - 2*E2(m) + E2(m-1)) + 2*E2(m) - E1(m);%BoundaryE(1) = 0; E(Nz) = 0; % Dirichlet%E(1) = E(2);E(Nz) = E(Nz-1); % Neumann%E(Nz) = E2(Nz-1) + (r-1)/(r+1)*(E(Nz-1)-E2(Nz)); % Mur ABC right side z = Nz%E(1) = E2(2) + (r-1)/(r+1)*(E(2)-E2(1)); % Mur ABC left side z = 0%display%choice***********************************************display = 3; % choice: '1' = line plot; '2' = imagesc; '3' = surf switch display case 1i = 1:Nz;plot(i, E(i));axis([0 Nz -4 4]);case 2A = rot90(E); imagesc(A);case 3i = 1:Nz;for j = 1:NzEs(i,j) = E(i);endEs = rot90(Es);j = 1:Nz;surf(i,j,Es);axis([0 Nz 0 Nz -4 4]) otherwiseA = rot90(E); imagesc(A);end pause(0.005);end% 3-D FDTD code with PECboundaries %*********************************************************************** %% Program author: Susan C. Hagness% Department of Electrical and Computer Engineering% University of Wisconsin-Madison% 1415 Engineering Drive% Madison, WI 53706-1691% 608-265-5739% %% Date of this version: February 2000%% This MATLAB M-file implements the finite-difference time-domain % solution of Maxwell's curl equations over a three-dimensional% Cartesian space lattice comprised of uniform cubic grid cells. %% To illustrate the algorithm, an air-filled rectangular cavity % resonator is modeled. The length, width, and height of the% cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and % 2.0 cm (z-direction), respectively.%% The computational domain is truncated using PEC boundary% conditions:% ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes% ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes% ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes% These PEC boundaries form the outer lossless walls of the cavity. %% The cavity is excited by an additive current source oriented% along the z-direction. The source waveform is a differentiated% Gaussian pulse given by% J(t)=-J0*(t-t0)*exp(-(t-t0F2 /tau A2),% where tau=50 ps. The FWHM spectral bandwidth of this zero-dc- % content pulse is approximately 7 GHz. The grid resolution% (dx = 2 mm) was chosen to provide at least 10 samples per % wavelength up through 15 GHz. %% To execute this M-file, type "fdtd3D" at the MATLAB prompt.% This M-file displays the FDTD-computed Ez fields at every other% time step, and records those frames in a movie matrix, M, which % is played at the end of the simulation using the "movie" command.%clear %*********************************************************************** % Fundamental constants%***********************************************************************cc=2.99792458e8; %speed of light in free space muz=4.0*pi*1.0e-7; %permeability of freespace epsz=1.0/(cc*cc*muz); %permittivity of free space%*********************************************************************** % Grid parameters %***********************************************************************ie=50; %number of grid cells in x-directionje=24; %number of grid cells in y-directionke=10; %number of grid cells in z-directionib=ie+1;jb=je+1;kb=ke+1;is=26; %location of z-directed current sourcejs=13; %location of z-directed current sourcekobs=5;dx=0.002; %space increment of cubic latticedt=dx/(2.0*cc); %time stepnmax=500; %total number of time steps%*********************************************************************** % Differentiated Gaussian pulse excitation%***********************************************************************rtau=50.0e-12;tau=rtau/dt;ndelay=3*tau;srcconst=-dt*3.0e+11;%*********************************************************************** % Material parameters%*********************************************************************** eps=1.0;sig=0.0;%***********************************************************************% Updating coefficients%***********************************************************************ca=(1.0-(dt*sig)/(2.0*epsz*eps))/(1.0+(dt*sig)/(2.0*epsz*eps));cb=(dt/epsz/eps/dx)/(1.0+(dt*sig)/(2.0*epsz*eps));da=1.0;db=dt/muz/dx;%***********************************************************************% Field arrays%***********************************************************************ex=zeros(ie,jb,kb);ey=zeros(ib,je,kb);ez=zeros(ib,jb,ke);hx=zeros(ib,je,ke);hy=zeros(ie,jb,ke);hz=zeros(ie,je,kb);%***********************************************************************% Movie initialization%***********************************************************************tview(:,:)=ez(:,:,kobs);sview(:,:)=ez(:,js,:);subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview');shading flat;caxis([-1.0 1.0]);colorbar;axis image;title(['Ez(i,j,k=5), time step = 0']);xlabel('i coordinate');ylabel('j coordinate');subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview');shading flat;caxis([-1.0 1.0]);colorbar;axis image; title(['Ez(i,j=13,k), time step = 0']); xlabel('i coordinate'); ylabel('k coordinate');rect=get(gcf,'Position');rect(1:2)=[0 0];M=moviein(nmax/ 2,gcf,rect);%*********************************************************************** % BEGIN TIME-STEPPING LOOP%***********************************************************************for n=1:nmax%*********************************************************************** % Update electric fields%***********************************************************************ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+...cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+... hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke));ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+... cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+... hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke));ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+... cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+... hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke));ez(is,js,1:ke)=ez(is,js,1:ke)+...srcconst*(n-ndelay)*exp(-((n-ndelay)A2 /tau A2));%*********************************************************************** % Update magnetic fields%***********************************************************************hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+... db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+...ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke));hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+... db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+...ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke));hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+... db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+...ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke));%*********************************************************************** % Visualize fields %***********************************************************************if mod(n,2)==0;timestep=int2str(n); tview(:,:)=ez(:,:,kobs); sview(:,:)=ez(:,js,:);subplot('position',[0.15 0.45 0.7 0.45]),pcolor(tview'); shading flat;caxis([-1.0 1.0]);colorbar;axis image;title(['Ez(i,j,k=5), time step = ',timestep]); xlabel('i coordinate'); ylabel('jcoordinate');subplot('position',[0.15 0.10 0.7 0.25]),pcolor(sview'); shading flat;caxis([-1.0 1.0]);colorbar;axis image;title(['Ez(i,j=13,k), time step = ',timestep]); xlabel('i coordinate'); ylabel('kcoordinate');nn=n/2;M(:,nn)=getframe(gcf,rect);end;%*********************************************************************** % END TIME-STEPPING LOOP%***********************************************************************endmovie(gcf,M,0,10,rect);。
FDTD中的MATlAB编程基础
• 2． pause函数：暂停程序的执行。
• 调用格式： pause(延迟秒数) • 注：如果省略延迟时间，直接使用pause，则将暂停程序，直到 用户按任一键后程序继续执行。
3.Drawnow函数： 将还未处理完的图像实时的显示出来。当代码执行时间
长，需要反复执行plot时，Matlab程序不会马上把图像画到figure上，这时，要 想实时看到图像的每一步变化情况，需要使用这个语句。
x、y、z 是向量时，plot3 命令的使用
t=0:0.1:8*pi; plot3(sin(t),cos(t),t) title(’绘制螺旋线’) %用命令 title 对图形主题进行标注 xlabel(’sin(t)’,’FontWeight’,’bold’,’FontAngle’,’italic’) ylabel(’cos(t)’,’FontWeight’,’bold’,’FontAngle’,’italic’) zlabel(’t’,’FontWeight’,’bold’,’FontAngle’,’italic’) %命令 zlabel 用来指定 z 轴的数据名称 grid on
常见矩阵生成函数
zeros(m,n) ones(m,n) eye(m,n) diag(X) tril(A) triu(A) 生成一个 m 行 n 列的零矩阵，m=n 时可简写为 zeros(n) 生成一个 m 行 n 列的元素全为 1 的矩阵, m=n 时可写为 ones(n) 生成一个主对角线全为 1 的 m 行 n 列矩阵, m=n 时可简写为 eye(n)，即为 n 维单位矩阵 若 X 是矩阵，则 diag(X) 为 X 的主对角线向量 若 X 是向量，diag(X) 产生以 X 为主对角线的对角矩阵 提取一个矩阵的下三角部分 提取一个矩阵的上三角部分
二维fdtd球坐标matlab代码
二维FDTD球坐标Matlab代码1. 简介在电磁学领域，有关于二维FDTD（时域有限差分）球坐标的Matlab 代码是一个重要的课题。
FDTD方法是一种数值求解Maxwell方程组的方法，通过网格化的形式将Maxwell方程组离散化为差分方程，进而求解电磁场在空间和时间上的变化。
而二维FDTD球坐标的应用则常见于球形结构的电磁场仿真，比如天线、雷达系统等。
本文将详细介绍如何使用Matlab编写二维FDTD球坐标的代码，并结合实例进行讲解。
2. 基本原理FDTD方法是一种求解Maxwell方程组的数值求解方法，它通过将Maxwell方程组离散化为差分方程，并采用逐步推进的方式求解电磁场在空间和时间上的变化。
在二维空间中，我们可以将电磁场的分布用网格进行离散化，通过更新电场和磁场的值来模拟电磁波在空间中的传播。
球坐标系是一种直角坐标系的补充，它在处理具有球对称性的问题时具有优势。
在二维FDTD球坐标中，我们通常会将空间离散化为一系列的环形网格，再结合球坐标系的变换关系来更新电场和磁场的值。
通过这种方式，我们可以较为准确地模拟球形结构上的电磁场分布。
3. Matlab代码实现下面我们将介绍如何使用Matlab编写二维FDTD球坐标的代码。
我们需要定义球坐标系下的电磁场更新方程，然后通过迭代来更新电场和磁场的数值，最终得到电磁场在空间和时间上的分布。
3.1 球坐标系下的电磁场更新方程在球坐标系下，电磁场的更新方程如下：（此处省略具体公式，具体可根据需要补充）3.2 迭代更新电场和磁场在Matlab中，我们可以通过嵌套循环来进行电场和磁场的迭代更新。
我们需要初始化电场和磁场的数值，然后通过迭代更新它们的值，最终得到电磁场在空间和时间上的分布。
```matlab此处为Matlab代码示例，具体实现可根据需要编写```4. 实例分析接下来，我们将通过一个实例来展示二维FDTD球坐标的Matlab代码。
假设我们希望模拟一个球形天线的电磁场分布，我们可以利用上述的Matlab代码来实现。
电磁场FDTD算法以及仿真图
function [ output_args ] = Untitled2( input_args )%UNTITLED2 Summary of this function goes here% Detailed explanation goes here%******************************************************************** ***% 3-D FDTD code with PEC boundaries%******************************************************************** ***%% Program author: Susan C. Hagness% Department of Electrical and Computer Engineering% University of Wisconsin-Madison% 1415 Engineering Drive% Madison, WI 53706-1691% 608-265-5739%*****************.edu%% Date of this version: February 2000%% This MATLAB M-file implements the finite-difference time-domain% solution of Maxwell's curl equations over a three-dimensional% Cartesian space lattice comprised of uniform cubic grid cells.%% To illustrate the algorithm, an air-filled rectangular cavity% resonator is modeled. The length, width, and height of the% cavity are 10.0 cm (x-direction), 4.8 cm (y-direction), and% 2.0 cm (z-direction), respectively.%% The computational domain is truncated using PEC boundary% conditions:% ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes% ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes% ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes% These PEC boundaries form the outer lossless walls of the cavity. %% The cavity is excited by an additive current source oriented% along the z-direction. The source waveform is a differentiated% Gaussian pulse given by% J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),% where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-% content pulse is approximately 7 GHz. The grid resolution% (dx = 2 mm) was chosen to provide at least 10 samples per% wavelength up through 15 GHz.%% To execute this M-file, type "fdtd3D" at the MATLAB prompt.% This M-file displays the FDTD-computed Ez fields at every other% time step, and records those frames in a movie matrix, M, which% is played at the end of the simulation using the "movie" command.%这个MATLAB的m文件实现了麦克斯韦旋度方程的有限差分方法在三维笛卡尔空间点阵组成的统一的立方网格细胞。
FDTD算法的Matlab程序
F D T D算法的M a t l a b程序(总7页)本页仅作为文档封面，使用时可以删除This document is for reference only-rar21year.March%*********************************************************************** % 3-D FDTD code with PEC boundaries%*********************************************************************** %% Program author: Susan C. Hagness% Department of Electrical and Computer Engineering% University of Wisconsin-Madison% 1415 Engineering Drive% Madison, WI 53706-1691% 608-265-5739%%% Date of this version: February 2000%% This MATLAB M-file implements the finite-difference time-domain% solution of Maxwell's curl equations over a three-dimensional% Cartesian space lattice comprised of uniform cubic grid cells.%% To illustrate the algorithm, an air-filled rectangular cavity% resonator is modeled. The length, width, and height of the% cavity are cm (x-direction), cm (y-direction), and% cm (z-direction), respectively.%% The computational domain is truncated using PEC boundary% conditions:% ex(i,j,k)=0 on the j=1, j=jb, k=1, and k=kb planes% ey(i,j,k)=0 on the i=1, i=ib, k=1, and k=kb planes% ez(i,j,k)=0 on the i=1, i=ib, j=1, and j=jb planes% These PEC boundaries form the outer lossless walls of the cavity.%% The cavity is excited by an additive current source oriented% along the z-direction. The source waveform is a differentiated% Gaussian pulse given by% J(t)=-J0*(t-t0)*exp(-(t-t0)^2/tau^2),% where tau=50 ps. The FWHM spectral bandwidth of this zero-dc-% content pulse is approximately 7 GHz. The grid resolution% (dx = 2 mm) was chosento provide at least 10 samples per% wavelength up through 15 GHz.%% To execute this M-file, type "fdtd3D" at the MATLAB prompt.% This M-file displays the FDTD-computed Ez fields at every other% time step, and records those frames in a movie matrix, M, which% is played at the end of the simulation using the "movie" command.%%***********************************************************************clear%*********************************************************************** % Fundamental constants%***********************************************************************cc=; %speed of light in free spacemuz=*pi*; %permeability of free spaceepsz=(cc*cc*muz); %permittivity of free space%*********************************************************************** % Grid parameters%***********************************************************************ie=50; %number of grid cells in x-directionje=24; %number ofgrid cells in y-directionke=10; %number of grid cells in z-directionib=ie+1;jb=je+1;kb=ke+1;is=26; %location of z-directed current sourcejs=13; %location of z-directed current sourcekobs=5;dx=; %space increment of cubic latticedt=dx/*cc); %time stepnmax=500; %total number of time steps%*********************************************************************** % Differentiated Gaussian pulse excitation%***********************************************************************rtau=;tau=rtau/dt;ndelay=3*tau;srcconst=-dt*+11;%*********************************************************************** % Material parameters%***********************************************************************eps=;sig=;%*********************************************************************** %Updating coefficients%***********************************************************************ca=(dt*sig)/*epsz*eps))/+(dt*sig)/*epsz*eps));cb=(dt/epsz/eps/dx)/+(dt*sig)/*epsz*eps));da=;db=dt/muz/dx;%*********************************************************************** % Field arrays%***********************************************************************ex=zeros(ie,jb,kb);ey=zeros(ib,je,kb);ez=zeros(ib,jb,ke);hx=zeros(ib,je,ke);hy=zeros(ie,jb,ke);hz=zeros(ie,je,kb);%*********************************************************************** % Movie initialization%***********************************************************************tview(:,:)=ez(:,:,kobs);sview(:,:)=ez(:,js,:);subplot('position',[ ]),pcolor(tview');shading flat;*****is([ ]);colorbar;axis image;title(['Ez(i,j,k=5), time step = 0']);xlabel('i coordinate');ylabel('j coordinate');subplot('position',[ ]),pcolor(sview');shading flat;*****is([ ]);colorbar;axis image;title(['Ez(i,j=13,k), time step = 0']);xlabel('icoordinate');ylabel('k coordinate');rect=get(gcf,'Position');rect(1:2)=[0 0];M=moviein(nmax/2,gcf,rect);%*********************************************************************** % BEGIN TIME-STEPPING LOOP%*********************************************************************** for n=1:nmax%*********************************************************************** % Update electric fields%***********************************************************************ex(1:ie,2:je,2:ke)=ca*ex(1:ie,2:je,2:ke)+...cb*(hz(1:ie,2:je,2:ke)-hz(1:ie,1:je-1,2:ke)+...hy(1:ie,2:je,1:ke-1)-hy(1:ie,2:je,2:ke));ey(2:ie,1:je,2:ke)=ca*ey(2:ie,1:je,2:ke)+...cb*(hx(2:ie,1:je,2:ke)-hx(2:ie,1:je,1:ke-1)+...hz(1:ie-1,1:je,2:ke)-hz(2:ie,1:je,2:ke));ez(2:ie,2:je,1:ke)=ca*ez(2:ie,2:je,1:ke)+...cb*(hx(2:ie,1:je-1,1:ke)-hx(2:ie,2:je,1:ke)+...hy(2:ie,2:je,1:ke)-hy(1:ie-1,2:je,1:ke));ez(is,js,1:ke)=ez(is,js,1:ke)+...srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));%*********************************************************************** % Update magnetic fields%***********************************************************************hx(2:ie,1:je,1:ke)=hx(2:ie,1:je,1:ke)+...db*(ey(2:ie,1:je,2:kb)-ey(2:ie,1:je,1:ke)+...ez(2:ie,1:je,1:ke)-ez(2:ie,2:jb,1:ke));hy(1:ie,2:je,1:ke)=hy(1:ie,2:je,1:ke)+...db*(ex(1:ie,2:je,1:ke)-ex(1:ie,2:je,2:kb)+...ez(2:ib,2:je,1:ke)-ez(1:ie,2:je,1:ke));hz(1:ie,1:je,2:ke)=hz(1:ie,1:je,2:ke)+...db*(ex(1:ie,2:jb,2:ke)-ex(1:ie,1:je,2:ke)+...ey(1:ie,1:je,2:ke)-ey(2:ib,1:je,2:ke));%*********************************************************************** % Visualize fields%*********************************************************************** if mod(n,2)==0;timestep=int2str(n);tview(:,:)=ez(:,:,kobs);sview(:,:)=ez(:,js,:);subplot('position',[ ]),pcolor(tview');shading flat;*****is([ ]);colorbar;axis image;title(['Ez(i,j,k=5), time step = ',timestep]);xlabel('i coordinate');ylabel('j coordinate');subplot('position',[ ]),pcolor(sview');shading flat;*****is([ ]);colorbar;axis image;title(['Ez(i,j=13,k), time step = ',timestep]);xlabel('i coordinate');ylabel('k coordinate');nn=n/2;M(:,nn)=getframe(gcf,rect);end;%*********************************************************************** % END TIME-STEPPING LOOP%*********************************************************************** endmovie(gcf,M,0,10,rect);。
用Matlab语言实现电磁场中FDTD法编程
n+1
( Q0 ) + f
n- 1
( P0 ) ] -
2
[f f
Δt - δ c [f δ
n
( P0 ) ] +
( P1 ) + f
n
(c Δt ) 2 ・ Δt + δ ) 2δ( c n ( P3 ) + f ( P4 ) +
n
n
( Q1 ) + f
( Q2 ) + f
( Q3 ) + f
一般来说 , 对于一个非扩频系统 , 时延越小越 好 ,时延越小 、 衰落周期越大 , 频率选择性够衰落越 轻。 但对于一个扩频系统 , 当时延小于并且越接近 一个码片周期时 , 首先影响解扩 , 不能进行多径分 离 ,在一个码片内的时延与码片宽度之比越大 ,影响 越大 。
完全匹配层为有耗介质 ,进入 PML 层的透射波 将迅速衰减 ,它对于入射波有很好的吸收效果 。在 不同条件下的具体推倒公式可参考文献 [ 2 ] 。
2 用 Matlab 方法实现编程
由于计算机容量的限制 ,FDTD 计算只能在有限 区域进行 。首先 , 根据计算目标的不同可以选择不 同的边界条件 。一般来说 , 如果对计算要求的精度 不高 ,且计算机性能较差的情况 , 选用 Mur 边界条 件 ; 如果对计算要求的精度比较高 ,且计算机硬件条 件较好 ,则可选用 PML 边界条件 。然后 ,根据 FDTD 的计算步骤进行以下准备工作 : ① FDTD 区的划分 。对于散射问题 ,FDTD 计算 区域划分为总场区和散射场区 。对于辐射问题 , 激 励源直接加到辐射天线上 , 整个 FDTD 计算区域为 辐射场 。可以借助等效原理应用计算区域内的近场 数据实现计算区域以外远场的外推 ; ② 求入射场在总场边界上的等效电磁场切向 分量 ,并写成差分格式 ,如上面列出的 2 种边界条件 的差分格式 ; ③将计算区域 ( 总场区和散射场区 ) 用网格离 散化 ,并用 FDTD 方程的差分离散形式按时间步写 出程序 ; ④根据选定的边界 条 件 写 出 边 界 和 角 点 的 程序 。 以真空中三维空间点电源的辐射为例 ; 如图 1 所示 。设 FDTD 计算区域为 100 3 100 3 100 个元胞 , 激励源为电偶极子 。 在三维 FDTD 计算中 , 电磁场节点可以区分为 3 种情况 : ①计算区域内部 ( 不包括截断边界) 节点 ;
FDTD算法的Matlab语言实现
FDTD算法的Matlab语言实现
张亮;寇晓艳
【期刊名称】《延安大学学报（自然科学版）》
【年(卷),期】2009(28)2
【摘要】分析了FDTD算法的基本原理及两种典型边界条件的算法特点,给出了Matlab语言编程的步骤和应注意的问题,并给出了实际的仿真结果,最后得出用Matlab语言对FDTD算法编程的几点结论.
【总页数】3页(P57-59)
【作者】张亮;寇晓艳
【作者单位】延安大学,信息学院,陕西,延安,716000;延安市医疗办,陕西,延
安,716000
【正文语种】中文
【中图分类】TP312
【相关文献】
1.模糊聚类分析算法的MATLAB语言实现 [J], 郭珉
2.BP算法用MATLAB语言的实现方法 [J], 黄君冉;钱东平;张嘏伟;程浩
3.用Matlab语言实现电磁场中FDTD法编程 [J], 田甜
4.基于Matlab语言实现电磁场中的FDTD算法编程 [J], 郑木生
5.费希尔算法的Matlab语言实现 [J], 陈晶晶;王艳
因版权原因，仅展示原文概要，查看原文内容请购买。
基于MATLAB的FDTD算法编程
基于MATLAB的FDTD算法编程
赵嘉
【期刊名称】《广西工学院学报》
【年(卷),期】2006(017)004
【摘要】介绍了时域有限差分(FDTD)法的基本原理,推导了二维TM模Yee算法的FDFD表达式,并结合算例阐述了基于MATLAB编程的基本方法.
【总页数】4页(P43-46)
【作者】赵嘉
【作者单位】广西工学院,计算机工程系,广西,柳州,545006
【正文语种】中文
【中图分类】TP321
【相关文献】
1.基于MATLAB交互式软件包实现对FDTD算法的改进 [J], 赵辉
2.基于FDTD技术计算目标远场RCS方法的Matlab实现 [J], 崔方;沈卫东
3.基于Matlab语言实现电磁场中的FDTD算法编程 [J], 郑木生
4.基于MATLAB的时域有限差分法(FDTD)的研究 [J], 李伟;赵春晖;徐娜
5.基于MATLAB的FDTD两种常见吸收边界条件的编程 [J], 赵金昌;童玲
因版权原因，仅展示原文概要，查看原文内容请购买。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
South China Normal University课程设计实验报告课程名称：计算电磁学指导老师：专业班级： 2014级电路与系统姓名：学号：FDTD算法的MATLAB语言实现摘要：时域有限差分(FDTD)算法是K.S.Yee于1966年提出的直接对麦克斯韦方程作差分处理,用来解决电磁脉冲在电磁介质中传播和反射问题的算法。
其基本思想是:FDTD计算域空间节点采用Yee元胞的方法,同时电场和磁场节点空间与时间上都采用交错抽样；把整个计算域划分成包括散射体的总场区以及只有反射波的散射场区,这两个区域是以连接边界相连接,最外边是采用特殊的吸收边界,同时在这两个边界之间有个输出边界,用于近、远场转换;在连接边界上采用连接边界条件加入入射波,从而使得入射波限制在总场区域；在吸收边界上采用吸收边界条件,尽量消除反射波在吸收边界上的非物理性反射波。
本文主要结合FDTD算法边界条件特点,在特定的参数设置下，用MATLAB语言进行编程，在二维自由空间TEz网格中，实现脉冲平面波。
关键词：FDTD；MATLAB；算法1 绪论1.1 课程设计背景与意义20世纪60年代以来，随着计算机技术的发展，一些电磁场的数值计算方法逐步发展起来，并得到广泛应用，其中主要有：属于频域技术的有限元法（FEM）、矩量法(MM)和单矩法等；属于时域技术方面的时域有限差分法(FDTD)、传输线矩阵法(TLM)和时域积分方程法等。
其中FDTD是一种已经获得广泛应用并且有很大发展前景的时域数值计算方法。
时域有限差分（FDTD）方法于1966年由K.S.Yee提出并迅速发展，且获得广泛应用。
K.S.Yee用后来被称作Yee氏网格的空间离散方式，把含时间变量的Maxwell旋度方程转化为差分方程，并成功地模拟了电磁脉冲与理想导体作用的时域响应。
但是由于当时理论的不成熟和计算机软硬件条件的限制，该方法并未得到相应的发展。
20世纪80年代中期以后，随着上述两个条件限制的逐步解除，FDTD便凭借其特有的优势得以迅速发展。
它能方便、精确地预测实际工程中的大量复杂电磁问题，应用范围几乎涉及所有电磁领域，成为电磁工程界和理论界研究的一个热点。
目前，FDTD日趋成熟，并成为分析大部分实际电磁问题的首选方法。
1.2 时域有限差分法的发展与应用经过四十多年的发展，FDTD已发展成为一种成熟的数值计算方法。
在发展过程中，几乎都是围绕几个重要问题展开的，即数值稳定性、计算精度、数值色散、激励源技术以及开域电磁问题的吸收边界条件等。
数值稳定和计算精度对任何一种数值计算方法都是至关重要的。
其中激励源的设计和引入也是FDTD的一个重要任务。
目前，应用最广泛的激励源引入技术是总场/散射场体系。
对于散射问题，通常在FDTD计算空间中引入连接边界，它将整个计算空间划分为内部的总场区和外部的散射场区，如图1.1。
利用Huygens原理，可以在连接边界处引入入射场，使入射场的加入变得简单易行。
图1.1总场/散射场区2 FDTD法基本原理2.1 Maxwell方程和Yee氏算法根据电磁场基本方程组的微分形式，若在无源空间，其空间中的煤质是各向同性、线性和均匀的，即煤质的参数不随时间变化且各向同性，则Maxwell旋度方程可写成：式中，E是电场强度，单位为伏/米（V/m）；H是磁场强度，单位为安/米（A/m）；ε表示介质介电系数，单位为法拉/米（F/m）；μ表示磁导系数，单位为亨利/米（H/m）；σ表示介质电导率，单位为西门子/米（S/m）；σ表示导磁率，单位m为欧姆/米（Ω/m）。
K.S.Yee在1966年建立了如图2.1所示的空间网格，这就是著名的Yee氏元胞网格。
图2.1 Yee氏网格及其电磁场分量分布电场和磁场被交叉放置,电场分量位于网格单元每条棱的中心,磁场分量位于网格单元每个面的中心,每个磁场(电场)分量都有4个电场(磁场)分量环绕。
这样不仅保证了介质分界面上切向场分量的连续性条件得到自然满足,而且还允许旋度方程在空间上进行中心差分运算,同时也满足了法拉第电磁感应定律和安培环路积分定律,也可以很恰当地模拟电磁波的实际传播过程。
2.2 时域有限差分法相关技术2.2.1数值稳定性问题FDTD方程是一种显式差分方程，在执行时，存在一个重要的问题：即算法的稳定性问题。
这种不稳定性表现为在解显式方程时，随着时间步数的继续增加，计算结果也将无限制地增加。
Taflove于1975年对Yee氏差分格式的稳定性进行了讨论，并导出了对时间步长的限制条件。
数值解是否稳定主要取决于时间步长Δt与空间步长Δx、Δy、Δz的关系。
对于非均匀媒质构成的计算空间选用如下的稳定性条件：上式是空间和时间离散之间应当满足的关系，又称为Courant稳定性条件。
若采用均匀立方体网格：Δx=Δy=Δz=Δs其中，v为计算空间中的电磁波的最大速度。
2.2.2数值色散FDTD方程组是对Maxwell旋度方程进行差分近似，在进行数值计算时，将会在计算网格中引起数字波模的色散，即在FDTD网格中，电磁波的相速与频率有关，电磁波的相速度随波长、传播方向及变量离散化的情况不同而改变。
这种关系由非物理因素引起，且色散将导致非物理因素引起的脉冲波形畸变、人为的各向异性和虚假折射等现象。
为获得理想的色散关系，问题空间分割应按照小于正常网格的原则进行。
一般选取的最大空间步长为Δmax=λmin/20，λmin为所研究范围内电磁波的最小波长。
由上分析说明，数值色散在用FDTD法分析电磁场传播中的影响是不可能避免的，但我们可以尽可能的减小数值色散的影响。
2.2.3离散网格的确定无论是简单目标还是复杂目标，在进行FDTD离散时网格尺寸的确定，除了受计算资源的限制不可能取得很小外，还需要考虑以下几个因素：1．目标离散精确度的要求。
网格应当足够小以便能精确模拟目标几何形状和电磁参数。
2．FDTD方法本身的要求。
主要是考虑色散误差的影响。
设网格为立方体Δx=Δy=Δz=δ，所关心频段的频率上限为fmax ，对应波长为fmin，则考虑FDTD的数值色散要求δ≤λmin/N通常N≥10。
上式是根据已知所关心频率上限情况下来确定FDTD网格尺寸的；反之，若给定，则FDTD计算结果可用的上限频率也随之确定。
2.3 吸收边界条件由时域有限差分法的基本原理可知，在利用时域有限差分法研究电磁场时，需在全部问题空间建立Yee氏网格空间，并存储每个单元网格上任一时间步的六个场分量用于下一时间步的计算。
而在对于辐射、散射这类开放系统的实际研究中，不可能有无限大的存储空间。
因此，必须在某处将网格空间截断，且在截断边界网格点处运用特殊的场分量计算方法，使得向边界面行进的波在边界处保持“外向行进”特性、无明显的反射现象，并且不会使内部空间的场产生畸变，从而用有限网格空间模拟电磁波在无界空间中传播的情况。
具有这种功能的边界条件称之为吸收边界条件，或辐射边界条件，或网格截断条件，如图2.2所示图 2.2 附加截断边界使计算区域变为有限域从FDTD的基本差分方程组可以看出，在截断边界面上切向场分量的计算需要利用计算空间以外的电磁场分量，因此FDTD基本差分方程对这些截断边界面上的场分量失效。
如何处理截断边界上的场分量，使之与需要考虑的无限空间有尽量小的差异，是FDTD中必须很好解决的一个重要问题。
实际上，这是要求在误差可容忍的范围内，计算空间中的外向波能够顺利通过截断边界面而不引起波的明显反射，使有限计算空间的数值模拟与实际情况趋于一致，对外向波而言，就像在无限大空间中传播一样。
所以，需要一种截断边界网格处的特殊计算方法，它不仅要保证边界场计算的必要精度，而且还要大大消除非物理因素引起的波反射，使得用有限的网格空间就能模拟电磁波在无限空间中的传播。
但是如果处理不当，截断边界面可能造成较大反射，构成数值模拟误差的一部分，甚至可能造成算法不稳定。
加于截断边界场分量符合上述要求的算法就称为吸收边界条件（Absorbing Boundary Conditions）。
2.4 FDTD计算所需时间步的估计为了使计算达到稳定,通常计算所需要时间步按照电磁波往返穿越FDTD计算区对角线3～5次来估计。
若FDTD计算区总元胞数为N，则对角线上元胞约为（三维)。
按照Courant稳定条件，设计算中心区Δt=δ/（2c），即穿越对角线一次需要时间步为。
总计算时间步约为需。
对于二维情况则约为。
或者说，计算时间步大约等于FDTD计算区对角线上元胞数目的12～20倍。
实际上，计算所需时间步还与散射体具体形状、结构有关。
图2.3给出了应用FDTD分析电磁场问题时的程序流程图图2.3 FDTD 程序流程图3 课程设计要求及MATLAB语言实现3.1课程设计要求➢采用总场/散射场技术，在二维自由空间TEz网格中，实现脉冲平面波。
➢总场区域100*100个网格。
➢周边散射场区域20个网格。
➢散射场区域外采用PML吸收边界条件。
➢平面波0度角入射。
使用高斯脉冲，中心频率为1GHz,带宽为1/e。
➢采用正方形Yee网格，离散步长Δ=1.5cm，时间离散步长为Δ=Δ/2c.3.2 MATLAB程序%%%%%%%% 此程序实现自由空间的平面波 %%%%%%%%%%%%%%%% 此程序采用总场/散射场技术 %%%%%%%%%%%%%%%% 此程序使用PML设置吸收边界条件 %%%%%%%%KE=128;%真空区域网格数IE=KE;JE=KE;T=0;NSTEP=300;%迭代次数e=2.71828;%%%%%%%% 初始化 %%%%%%%%dz=zeros(KE,KE);ez=zeros(KE,KE);hx=zeros(KE,KE);hy=zeros(KE,KE);ihx=zeros(KE,KE);ihy=zeros(KE,KE);ga=ones(KE,KE);ez_inc=zeros(1,JE);hx_inc=zeros(1,JE);%%%%%%%% PML 参数 %%%%%%%%gi2=ones(1,IE);gi3=gi2;fi1=zeros(1,IE);fi2=ones(1,IE);fi3=fi2;gj2=ones(1,JE);gj3=gj2;fj1=zeros(1,JE);fj2=ones(1,JE);fj3=fj2;for i=1:8 %设置PML层中的参数xnum=npml-i;xd=npml;xxn=xnum/xd;xn=0.33*xxn^3;gi2(i)=1.0/(1.0+xn);gi2(IE-1-i)=1.0/(1.0+xn);gi3(i)=(1.0-xn)/(1.0+xn);gi3(IE-1-i)=(1.0-xn)/(1.0+xn);xxn=(xnum-0.5)/xd;xn=0.25*xxn^3;fi1(i)=xn;fi1(IE-2-i)=xn;fi2(i)=1.0/(1.0+xn);fi2(IE-2-i)=1.0/(1.0+xn);fi3(i)=(1.0-xn)/(1.0+xn);fi3(IE-2-i)=(1.0-xn)/(1.0+xn); endfor j=1:8xnum=npml-j;xd=npml;xxn=xnum/xd;xn=0.33*xxn^3;gj2(j)=1.0/(1.0+xn);gj2(JE-1-j)=1.0/(1.0+xn);gj3(j)=(1.0-xn)/(1.0+xn);gj3(JE-1-j)=(1.0-xn)/(1.0+xn);xxn=(xnum-0.5)/xd;xn=0.25*xxn^3;fj1(j)=xn;fj1(JE-2-j)=xn;fj2(j)=1.0/(1.0+xn);fj2(JE-2-j)=1.0/(1.0+xn);fj3(j)=(1-xn)/(1.0+xn);fj3(JE-2-j)=(1.0-xn)/(1.0+xn); end%%%%%%%% 总场/散射场边界 %%%%%%%%ib=IE-ia-1;ja=28;jb=JE-ja-1;ez_inc_low_m1=0;ez_inc_low_m2=0;ez_inc_high_m1=0;ez_inc_high_m2=0;%%%%%%%% 信号源 %%%%%%%%t0=80;%脉冲高度spread=12;%脉冲宽度%%%%%%%% 网格 %%%%%%%%ddx=0.015; %%%%%%%% 离散步长 %%%%%%%%dt=ddx/(2*3e8); %%%%%%%% 计算时间离散步长 %%%%%%%%epsz=8.851*1e-12; %%%%%%%% 自由空间的介电常数 %%%%%%%%%%%%%% 迭代求解电场与磁场 %%%%%%%%for n=1:NSTEPa=1;T=T+1;for j=2:JEez_inc(j)= ez_inc(j)+0.5*( hx_inc(j-1)- hx_inc(j));end%%%%%%%% 入射波缓冲区 %%%%%%%%ez_inc(1) =ez_inc_low_m2;ez_inc_low_m2=ez_inc_low_m1;ez_inc_low_m1=ez_inc(2);ez_inc(JE-1) =ez_inc_high_m2;ez_inc_high_m2=ez_inc_high_m1;ez_inc_high_m1=ez_inc(JE-2);%%%%%%%% 计算 dx 区域 %%%%%%%%for i=2:KEfor j=2:KEdz(i,j)=gi3(i)*gj3(j)*dz(i,j)+gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx (i,j)+hx(i,j-1));endend%%%%%%%% 输入高斯脉冲信号源 %%%%%%%%pulse=exp(-0.5*(t0-T)^2/spread^2); %%%%%%%%ez_inc(3)=pulse;%%%%%%%% 入射波 DZ 参数 %%%%%%%%for i=ia:ibdz(i,ja)= dz(i,ja)+0.5*hx_inc(ja-1);dz(i,jb)= dz(i,jb)-0.5*hx_inc(jb);end%%%%%%%% 计算电场dx区域 %%%%%%%%for i=1:KEfor j=1:KEez(i,j)=ga(i,j)*dz(i,j);endend%%%%%%%% 设置电场 EZ 边缘为0, 作为PML的参数%%%%%%%%for j=1:JE-1ez(1,j)=0;ez(IE,j)=0;endfor i=1:IE-1ez(i,1)=0;ez(i,JE)=0;endfor j=1:JE-1hx_inc(j)=hx_inc(j)+0.5*( ez_inc(j)-ez_inc(j+1));end%%%%%%%% 计算磁场hx区域 %%%%%%%%for i=1:KEfor j=1:KE-1curl_e=ez(i,j)-ez(i,j+1);ihx(i,j)=ihx(i,j)+fi1(i)*curl_e;hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5* (curl_e+ihx(i,j)); endend%%%%%%%% 入射磁场Hx区域参数 %%%%%%%%for i=ia:ibhx(i,ja-1)= hx(i,ja-1)+0.5*ez_inc(ja);hx(i,jb)= hx(i,jb)-0.5*ez_inc(jb);end%%%%%%%% 计算磁场hy区域%%%%%%%%for i=1:KE-1for j=1:KEcurl_e=ez(i+1,j)-ez(i,j);ihy(i,j)=ihy(i,j)+fj1(j)*curl_e;hy(i,j)=fi3(i)*hy(i,j)+fi2(i)*0.5* (curl_e+ihy(i,j)); endend%%%%%%%% 入射磁场hy区域参数 %%%%%%%%for j=ja:jbhy(ia-1,j)= hy(ia-1,j)-0.5*ez_inc(j);hy(ib,j)= hy(ib,j)+0.5*ez_inc(j);end%%%%%%%% 图形展示 %%%%%%%%m=moviein(500);x=1:KE;y=1:KE;[X,Y]=meshgrid(x,y);surf(X,Y,ez);axis([0 KE 0 KE -1 1 ]);set(gcf,'Color', 'white', 'Number', 'off', 'Name', sprintf('Simulation FDTD 2D, Iteration = %i', n));title( sprintf('t = %.3f nsec',n*dt*1e9))xlabel('x');ylabel('y');zlabel('Ez');m(a)=getframe(gcf);a=a+1;end下图为该程序实现的效果图：4.结语以上结合FDTD和MATLAB在自由空间中实现了平面波，所编MATLAB程序简洁明了，运行效率也较高。