FDTD中的MATlAB编程基础
Matlab编程的基础知识详解
Matlab编程的基础知识详解一、引言Matlab是一种高效且强大的数值计算软件，被广泛应用于科学、工程和金融等领域。
本文将详细介绍Matlab编程的基础知识，包括变量、数据类型、数组和矩阵操作、控制流程和函数等方面的内容。
二、变量和数据类型在Matlab中，变量是用来存储数据的容器。
在定义变量时，需要为其指定一个名称，并给其赋予一个值。
Matlab中常用的数据类型包括数值类型、字符型和逻辑型。
数值类型包括整型（int）、浮点型（double）、复数型（complex）等。
字符型用于存储文本信息，逻辑型用于存储逻辑值（true/false）。
变量可以通过赋值运算符“=”进行赋值操作。
例如，可以使用语句“x = 10”将变量x的值设置为10。
三、数组和矩阵操作在Matlab中，数组是一个包含相同类型数据的集合。
矩阵是特殊的数组，是一个二维表格，其中的元素可以通过行和列的索引进行访问。
Matlab提供了丰富的数组和矩阵操作函数，用于对数据进行变换、运算和统计分析。
例如，可以使用“size”函数获取数组的大小，使用“transpose”函数进行矩阵转置，使用“reshape”函数改变矩阵的形状等。
四、控制流程控制流程用于控制程序的执行流程，包括条件判断和循环结构。
条件判断使用“if-else”语句，用于在特定条件下执行不同的代码块。
例如，可以使用“if x>0”判断变量x是否大于0，如果成立则执行相应的代码块，否则执行其他代码块。
循环结构用于重复执行特定的代码块。
常用的循环结构有“for”循环和“while”循环。
例如，可以使用“for i=1:10”循环语句执行一个代码块10次。
五、函数函数是一段具有特定功能的代码块，可以反复利用。
在Matlab中，可以使用内置函数或自定义函数。
使用内置函数可以实现诸如数学运算、数据分析和图形绘制等功能。
例如，可以使用“sin”函数计算正弦值，使用“mean”函数计算平均值。
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编程基础入门
MATLAB编程基础入门MATLAB是一种常用于科学计算和数据分析的高级编程语言和环境。
它提供了丰富的工具集，使得处理数值计算、绘制图形以及实现算法变得更加便捷。
本文将为初学者介绍MATLAB的基础知识和编程技巧，以帮助读者快速入门。
1. MATLAB的安装和启动首先，我们需要到MathWorks官网上下载并安装MATLAB。
安装完成后，双击MATLAB图标即可启动软件。
MATLAB的主界面分为命令窗口、编辑器和工作空间等几个主要部分，用户可以通过这些界面进行编程和运行程序。
2. MATLAB的基本语法MATLAB的基本语法与其他编程语言有所不同。
在MATLAB中，不需要声明变量的类型，只需要直接给变量赋值即可。
例如：```a = 10;b = 3.14;c = 'Hello, MATLAB!';```MATLAB中还有一些特殊变量和函数，比如`pi`表示圆周率，`sin`表示正弦函数。
使用这些特殊变量和函数可以实现更加高效的数值计算和数据处理。
3. MATLAB的基本操作MATLAB提供了丰富的操作符和函数，可以用于数值计算、矩阵运算、图形绘制等。
下面是一些常用操作的示例：3.1 数值计算```a = 5;b = 3;c = a + b; % 加法运算d = a * b; % 乘法运算e = sqrt(a); % 开方运算```3.2 矩阵运算```A = [1 2 3; 4 5 6; 7 8 9]; % 创建一个3x3的矩阵B = [10 11 12; 13 14 15; 16 17 18];C = A + B; % 矩阵相加D = A * B; % 矩阵相乘```3.3 图形绘制```x = linspace(0, 2*pi, 100); % 在0到2π之间生成100个等间隔的点y = sin(x);plot(x, y); % 绘制正弦函数图像xlabel('x'); % 设置x轴标签ylabel('y'); % 设置y轴标签title('Sin Function'); % 设置图像标题```4. MATLAB的程序控制MATLAB提供了丰富的控制结构，可以用于实现条件判断和循环等功能。
MATLAB编程的基础与实践
MATLAB编程的基础与实践MATLAB是一种常用的科学与工程计算软件，它提供了强大的数学计算、数据可视化和算法开发等功能。
在科学研究和工程设计中，MATLAB有着广泛的应用。
本文将介绍MATLAB编程的基础知识和实践技巧，并给出一些实用例子。
一、MATLAB的基础语法MATLAB的基础语法与众多编程语言类似，主要包括数据类型、变量、运算符和控制结构等。
下面我们逐一介绍。
1.数据类型MATLAB支持常见的数据类型，包括数值型、字符型、逻辑型和结构型。
其中，数值型包括整型和浮点型，用于存储数值数据；字符型用于存储文本数据；逻辑型用于存储True或False两种取值；结构型用于存储复杂数据结构。
2.变量在MATLAB中，变量指的是存储数据的内存单元。
变量必须以字母或下划线开头，可以包含数字和字母，不能使用 MATLAB 关键字作为变量名。
变量名的长度可以任意，但是只有前面几个字符才会被识别。
3.运算符MATLAB支持常见的算术运算符、关系运算符和逻辑运算符。
其中，算术运算符包括加减乘除和幂运算；关系运算符包括等于、大于、小于等；逻辑运算符包括与、或和非等。
4.控制结构在MATLAB中，控制结构包括条件语句和循环语句。
条件语句用于实现根据条件决定执行不同代码块，常见的有if语句和switch语句。
循环语句用于实现重复执行某一段代码，常见的有for循环和while循环。
二、MATLAB的实用技巧除了基础语法之外，MATLAB编程需要掌握一些实用技巧，以提高编程效率和实现更复杂的功能。
1.函数的定义和调用在MATLAB中，函数是实现特定功能的代码块。
通过定义函数，我们可以将复杂的计算拆分成多个可重复使用的小模块，从而提高编程效率。
函数的定义包括函数名、参数和主体代码。
函数的调用可以通过在命令行窗口中输入函数名和参数值来实现。
2.数组和矩阵操作MATLAB是一个优秀的数值计算软件，因此数组和矩阵的操作是编程中常见的任务。
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编程基础教程在计算机科学和工程领域，MATLAB（Matrix Laboratory）是一个广泛用于数值分析和科学计算的高级编程语言和环境。
它的强大功能和简洁的语法使得它成为许多科学家和工程师的首选工具。
本文将介绍MATLAB的基础知识，帮助读者快速入门并进行简单的编程。
1. MATLAB的安装与环境配置首先，我们需要下载并安装MATLAB软件。
MATLAB可以在官方网站上免费获取到，并提供不同的版本供选择。
安装程序非常简单，只需按照向导的指示进行操作即可完成安装。
安装完成后，我们需要进行一些环境配置。
首先启动MATLAB软件，然后选择合适的工作目录。
工作目录是我们存储和管理MATLAB文件的地方。
选择一个方便和易于查找的目录，并将其设置为工作目录。
接下来，我们还可以对编辑器的外观和功能进行自定义设置，以适应个人的需要。
2. MATLAB的基本语法和语句在MATLAB中，所有的操作都是通过输入命令来完成的。
MATLAB的命令由一个或多个关键字组成，可以用于执行各种操作，包括数值计算、数据可视化和文件处理等。
下面是一些常用的MATLAB命令示例：- disp('Hello, world!')：显示一个文本消息- a = 1 + 2：将1与2相加，并将结果保存到变量a中- b = sqrt(9)：计算9的平方根，并将结果保存到变量b中- c = linspace(1, 10, 10)：生成一个由1到10的10个等间距数字组成的向量，并将结果保存到变量c中MATLAB还提供了丰富的数学函数和运算符，可以进行各种数值计算操作。
例如，可以使用'+'运算符进行加法运算，使用'-'运算符进行减法运算，使用'*'运算符进行乘法运算，使用'/'运算符进行除法运算等。
此外，MATLAB还提供了一些特殊的函数，如sin、cos、exp、log等，用于实现各种数学运算。
MATLAB编程基础
MATLAB编程基础MATLAB是一种高级的科学计算软件，广泛应用于工程、科学和数学领域。
作为一种强大的计算工具，掌握MATLAB编程基础对于提高科学研究和工程设计的效率至关重要。
一、MATLAB介绍与安装MATLAB是Matrix Laboratory的缩写，其最大的特点是擅长矩阵运算和数据可视化。
作为一种商业软件，需要购买正版授权才能使用。
安装MATLAB非常简单，只需按照官方指引进行安装，然后输入授权码即可启动。
同时，MATLAB还提供了独立开发环境MATLAB Integrated Development Environment（简称MATLAB IDE），方便用户进行编写、调试和运行代码。
二、MATLAB基本语法MATLAB使用类似英语的简明语法，不需要像其他编程语言那样严格遵循特定的规范。
基本的MATLAB语法包括变量声明、赋值、运算符、条件语句和循环语句。
例如，使用"="符号可以将数值赋给变量，并使用"=="符号判断两个数值是否相等。
此外，还可以使用"+、-、*、/"等运算符进行数值计算。
三、MATLAB函数与脚本在MATLAB中，函数是一种可重复使用的代码块，用于封装特定的计算任务。
通过定义函数，可以编写更加模块化的程序。
在函数中，可以设置输入参数和输出参数，并且可以在函数中调用其他函数进行嵌套计算。
与函数不同，脚本是一种顺序执行的代码片段，它可以直接在MATLAB命令窗口中运行。
脚本通常用于简单的计算和快速验证算法。
四、MATLAB图形界面与数据可视化MATLAB具有强大的图形界面工具，可用于数据可视化和图形绘制。
通过简单的MATLAB命令，可以绘制二维和三维的曲线、散点图、柱状图等。
同时，还可以自定义绘图风格、添加标签和标题，使得图形更加美观而富有表达力。
图形界面工具还支持交互式操作，用户可以通过鼠标和键盘进行数据点的选取、拖动和缩放。
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代码，可以对复杂的电磁场问题进行数值分析，得到定量的仿真结果，为工程设计和科学研究提供重要的参考。
基于MATLAB的FDTD算法编程
（） ３
方程中出现了半个网格和半个时间步 ， 为便于编程 ， 可将上面差分式改写为如下 Ｆ Ｔ Ｄ Ｄ算式 ：
Ｈ ， ＋ １ － Ｈ ， ） （， ｋ ） （， 愚 一
ｔ ，
￣Ｅ ｌ 一 ＿］ ｙ ｊ ， Ｅ，） Ｅ， （， （＋ ） 『
［ ， 一Ｅ （，，） Ｅ （＋１ ，） ］
，
８７ １ｚ
ａ ￡ Ｈ
’
Ｂ Ｅ． ａ Ｈ
。 。 一
８ Ｈ
ａ ｖ
３ 一 一 ｙ
构 造 二维 Ｔ 波 ＹＥ ｅ 如 图 ｌ 示 ： Ｍ Ｅｃｌ ｌ 所
按 Ｙ Ｅ元胞对上式偏导用中心差商代替， Ｅ 可得 ：
生
Ａｙ
一
生１ — １ １ — １ ． ＋
△￡
— — — — — — — ～ ： ：： — — — — — — — — — — — — — 一
（） ２ Leabharlann ｏ —Ｅ ，＋ ） ， 丢 （， 吉一 （ ） Ｅ ， 一
— — — — — — — — 一
日（ ｌ， 丢一 十 ）Ｈ（ ＋， 专 日（『＋ ） ， ＋ ）Ｈ（ ， 号 ， ｌ＋ 一 ， ＋ ， ｊ ， ＿
时间轴和空间轴上采取蛙跳法 （ ａｆ ｇ ： 步推进地求解 ， １ ｐｒ ） ｅ ｏ ￣ 最终求出一定边值与初值条件下的空间场解 。随
着计 算机 技术 的发 展 ， 近年来 Ｆ Ｄ计 算技 术 也得到 了越 来越 多 的应用 。 ＤＴ 对于 Ｆ Ｔ 算法 的编程 求解 ， Ｄ Ｄ 最常 用 的 有 ＶＣ和 Ｆ ＴＲ ＯＲ ＵＮ， ＭＡＴ ＡＢ作 为一 种可 视化 效果 很好 的科学 计算 软件 ， Ｆ Ｔ 而 Ｌ 在 Ｄ Ｄ计算 中能充分 发 挥 编程 简单 、 可视 化 程度 高 、 显示 动态 场 能
matlab编程入门基础
例： A=input('Please input A: ')
例： name=input('What''s your name? ')
输入字符串时必须带单引号
单引号的输出
最新课件
20
disp
数据的输出：disp
disp(X)
00
0
0
01
0
1
10
0
1
11
1
1
非 异或
~A Xor(A,B)
1
0
1
1
0
1
0
0
在 Matlab 中，0 表示 “假”，非零表示 “真”
最新课件
14
逻辑运算
相关函数
any(x)
如果向量 X 中存在非零元素，则返回 1， 否则返回 0
all(x)
如果向量 X 中所有元素都非零，则返回 1， 否则返回 0
"|"与“||”同理。
最新课件
11
逻辑运算
逻辑运算符
A&&B 首先判断A的逻辑值，如果A的值为假，就可以判断整 个表达式的值为假，就不需要再判断B的值。这种用法非常 有用，如果A是一个计算量较小的函数，B是一个计算量较大 的函数，那么首先判断A对减少计算量是有好处的。
另外这也可以防止类似被0除的错误：
在 Matlab 程序设计中，要充分利用 Matlab 数据结构的 特点，提高编程效率
最新课件
2
主要内容
M 文件介绍
Matlab 编程基础
算术运算、关系运算、逻辑运算 控制结构：顺序结构、选择结构、循环结构
二维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的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*10^9); % 真空介电常数mu=1/(c0^2*epsz);%磁常数ex_low_m1=0;ex_low_m2=0;ex_high_m1=0;ex_high_m2=0;a0=2*u0/dt+(2/dt)^2;a1=-8/(dt)^2;a2=-2*u0/dt+(2/dt)^2;b0=wpe^2+2*u0/dt+(2/dt)^2;b1=2*wpe^2-8/(dt)^2;b2=wpe^2-2*u0/dt+(2/dt)^2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%初始化电磁场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-t0)/d)^2);%%%计算电场Exfor i=1:kpstart-1Ex(i)=Dx(i)/epsz;endfor i=kpstop+1:KEEx(i)=Dx(i)/epsz;endDx3=Dx2;Dx2=Dx1;Dx1=Dx;for i=kpstart:kpstopEx(i)=(1/b0)*(a0*Dx1(i)+a1*Dx2(i)+a2*Dx3(i)-b1*Ex1(i)-b2*Ex 2(i));endEx2=Ex1;Ex1=Ex;Sx3=Sx2;Sx2=Sx1;Ex(1)=ex_low_m2;ex_low_m2=ex_low_m1;ex_low_m1=Ex(2);Ex(KE)=ex_high_m2;ex_high_m2=ex_high_m1;ex_high_m1=Ex(KE-1);%%%%%%%%%%%%%%%%%%计算磁场for i=1:KE-1Hy(i)=Hy(i)-(dt/(mu*zdelta))*(Ex(i+1)-Ex(i));endplot(Ex);grid on;pause(0.01);end% FDTD Main Function Jobs to Workers %%*********************************************************************** % 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 of the% 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)^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.%%*********************************************************************** 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%***********************************************************************ie = get(handles.xslider,'Value'); %number of grid cells in x-directionje = get(handles.yslider,'Value'); %number of grid cells in y-directionke = get(handles.zslider,'Value'); %number of grid cells in z-directionib=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 spaceparam.epsz = 1.0/(**param.muz); %permittivity of free space%*********************************************************************** % Grid parameters%***********************************************************************param.is = get(handles.xsource,'Value'); %location of z-directed current source param.js = get(handles.ysource,'Value'); %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 pulse excitation%***********************************************************************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;%***********************************************************************% Material parameters%***********************************************************************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.e psz*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_FFTpause(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.PNfor 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 == 1h = warndlg('Simulation Run is Stopped by User !','!! 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);endend%文件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 tightset(gca,'nextplot','replacechildren');E_total = sqrt(ex.^2 + ey.^2 + ez.^2);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 New 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 tightset(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','Times New Roman','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 = fullfile('C:\MATLAB7\work\cavity\figures',strcat('XZ_ETotal',num2str(n),'.jpeg')); % imwrite(M2,filename,'jpeg')%***********************************************************************% Cavity Power - Analitic expression%***********************************************************************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 New Roman','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)^2/param.tau^2));case 1source = param.srcconst*sin(param.dt*n*2*pi*f);case 3source = param.srcconst*exp(-((n-param.ndelay)^2/param.tau^2))*sin(2*pi*f*(n-param.ndelay)*param.dt) ;otherwisesource = param.srcconst*(n-param.ndelay)*exp(-((n-param.ndelay)^2/param.tau^2)); 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));%文件5function [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-t0).^2./T^2);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_debugc_0 = 3E8; % Speed of light in free spaceNz = 100; % Number of cells in z-directionNt = 200; % Number of time stepsN = Nz; % Global number of cellsE = zeros(Nz,1); % E componentE2 = zeros(Nz,1);E1 = zeros(Nz,1);Es = zeros(Nz,Nz); % Es is the output for surf function%H = zeros(Nz,1); % H componentpulse = 0; % Pulse%to be setmu_0 = 4.0*pi*1.0e-7; % Permeability of free spaceeps_0 = 1.0/(c_0*c_0*mu_0); % Permittivty of free spacec_ref = c_0; % Reference velocityeps_ref = eps_0;mu_ref = mu_0;f_0 = 10e9; % Excitation frequencyf_ref = f_0; % Reference frequencyomega_0 = 2.0*pi*f_0; % Excitation circular frequencyomega_ref = omega_0;lambda_ref = c_ref/f_ref; % Reference wavelengthDx_ref = lambda_ref/20; % Reference cells widthDz = Dx_ref;nDz = Dz/Dx_ref;Z = Nz*Dz;r = 1; % Normalized time stepDt_ref = r*Dx_ref/c_ref; % Reference time stepDt = 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;case 2 % sigle cos source functionif t < 2.0*pi/(omega_0)pulse = 8*c_0^2*Dt_ref^2*mu_ref*omega_0*cos(omega_0*t);elsepulse = 0;endcase 3 % Gaussian pulseif t < Dt_ref*r*50pulse = -40*c_0*(t-t*25/(n-1))*exp(-(t-t*25/(n-1))^2/2/(50/2.3548)^2/(t/(n-1))^2);elsepulse = 0;endendE2(N_Source) = E2(N_Source) - r*pulse;E1 = E2;E2 = E;m = 2:Nz-1;E(m) = r^2*(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' = surfswitch displaycase 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);endpause(0.005);end%***********************************************************************% 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 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.%%***********************************************************************clear%*********************************************************************** % Fundamental constants%***********************************************************************cc=2.99792458e8; %speed of light in free spacemuz=4.0*pi*1.0e-7; %permeability of free spaceepsz=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)^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',[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('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 = ',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法编程
使 用 Ｍ ｆｂ语 言 编 写 的 计 算 在 真 空 中 的 三 维 空 间 点 电源 辐 射 电 场 的 算 例 ，记 述 了 区 域 划 分 、 写 有 限 差 分 式 、 编 写 程 ａａ ｌ 序 ３个 主 要 计 算 步 骤 。 列 出 了 在 Ｆ Ｔ 网 格 划 分 中 一 些 参 数 的 选 取 以 及 Ｍ ｔ ｂＰ Ｅ 工 具 箱 的 网 格 生 成 器 。 得 出 用 ＤＤ ａａ Ｄ ｌ Ｍ ｔｂ语 言对 Ｆ Ｔ 算 法 编 程 的 几 点 结 论 。 ａａ ｌ ＤＤ 关键词 Ｆ Ｔ 边 界 条 件 ； ａａ 言 ； 法 编 程 Ｄ Ｄ； Ｍ ｄｂ语 算
Ｋｅ ｒ ｓ Ｆ ｙ ｗｏ ｄ ＤＴＤ ； ｕ ｄ ｒ ｏ ｄ ｔ ｎ； ｄａ ａ ｇ ａ ｅ；ｒｔ ｍｅｉ ｒｇａ ｏ ｂ ｎ ａ ｃ ｎ ｉｉ ｙ ｏ Ｍａ ｂ ｌｎ ｕ ｇ ａｉ ｈ ｔｃ ｐｏ ｍｍｉｇ ｒ ｎ
０ 引言
时域 有 限 差 分 法 （ ｉ ｔ Ｄｆ ｒ ｃ ｉ ｅＤ ｍ ｉ Ｆｎ ｅ ｉｅｅ ｅＴｍ ． ｏ ａ ｉ． ｆ ｎ ｎ
Ｍｅ ｏ ） ｔ ｄ 是求解 电磁 问题 的一种 数值 技术 ， 在 １６ ｈ 是 ９６ 年 由 Ｋ． ． ｅ 第一 次 提 出 的 。Ｆ Ｔ ＳＹ ｅ Ｄ Ｄ法 直接 将 有 限 差分式 代替麦 克 斯 韦时 域 场旋 度 方 程 中 的微 分式 ， 得到关 于场分 量 的有 限 差 分式 ， 具 有 相 同 电参 量 用 的空 间网格去模 拟 被 研 究 体 ， 取 合 适 的场 始 值 和 选 计算 空 间的边界 条 件 ， 以得 到包 括 时 间变 量 的麦 可 克斯 韦方 程 的四维 数值 解 ， 通过 傅 里 叶变 换 可求 得 三 维空 间的频 域解 。时域 有限差 分法突 出 的优 点是 所 需 的存 储量 及计 算时 间 与 Ⅳ 成正 比 ， 得很 多复 使 杂 的电磁场 计算 问题成 为可 能 ， 时域 有 限差 分 法 用 容易模 拟各 种复 杂 的结 构 ， 得 用 其他 方法 不 能 解 使 决 的问题有 了新 的处 理方法 。
matlab的dstep函数用法
一、dstep函数概述Matlab中的dstep函数是一个用于求离散系统单位阶跃响应的函数。
它可以帮助用户分析和验证离散系统在单位阶跃输入下的响应情况。
在控制系统工程中，dstep函数是一个非常重要的工具，能够帮助工程师进行模拟、分析和设计。
二、dstep函数的基本语法dstep函数的基本语法如下：[R, L, X] = dstep(SYS)其中，SYS是一个离散系统的状态空间模型，R是单位阶跃响应的输出向量，L是单位阶跃响应的时间向量，X是单位阶跃响应的状态向量。
三、dstep函数的参数说明1. SYS：离散系统的状态空间模型，可以通过ss函数创建。
离散系统的状态空间模型包括状态矩阵A、输入矩阵B、输出矩阵C和传递矩阵D。
2. R：单位阶跃响应的输出向量，包括系统对单位阶跃输入的响应值。
3. L：单位阶跃响应的时间向量，表示单位时间内系统的响应情况。
4. X：单位阶跃响应的状态向量，包括系统的状态响应值。
五、dstep函数的使用示例下面是一个使用dstep函数的简单示例，假设有一个离散系统的状态空间模型如下：A = [1, 0.1; 0, 1]B = [0.01; 0.1]C = [1, 0]D = 0我们可以通过以下代码求取该系统的单位阶跃响应：SYS = ss(A, B, C, D)[R, L, X] = dstep(SYS)六、dstep函数的注意事项1. 在使用dstep函数时，需要输入一个有效的离散系统状态空间模型，确保该模型符合实际工程系统的特性。
2. dstep函数计算的是离散系统对单位阶跃输入的响应，因此在实际工程应用中，需要结合其他方法和工具对系统进行综合分析和设计。
七、总结通过上述介绍，我们了解了Matlab中dstep函数的基本用法和注意事项。
在工程领域中，dstep函数是一个非常有用的工具，能够帮助工程师分析离散系统的响应特性，从而进行模拟、分析和设计工作。
当然，在实际应用中，也需要结合其他工具和方法，进行全面的系统分析和设计。
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算法编程
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编程基础
• 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 为主对角线的对角矩阵 提取一个矩阵的下三角部分 提取一个矩阵的上三角部分
Matlab编程入门（一）：编程基础
Matlab编程⼊门（⼀）：编程基础上学期学了⼀些matlab的知识，这学期再⽤时竟然发现已经忘得差不多了(┬＿┬)于是决定重新开始并将它们记录下来，也⽅便⾃⼰以后查漏补缺！M⽂件编程脚本⽂件 matlab有⾃⼰的命令⾏窗⼝，对于简单的命令，可以直接在命令⾏窗⼝输⼊，但随着命令⾏的增加或者命令本⾝复杂度的增加，再使⽤命令⾏就显得有些不便了，这时就需要脚本⽂件了。
可以说，脚本⽂件是matlab指令集合的封装。
函数⽂件 函数⽂件以function开始，end结束，这也是区别于脚本⽂件的地⽅。
在function后⾯接着定义输出参数，函数名和输⼊参数，⽐如： function [x,y,z] = math_count(a,b,c) x，y，z是输出参数，以⽅框括起来，math_count是函数名，a，b，c是输⼊参数，以圆括号括起来。
也可以没有参数，⽐如： function printresults(x,y) printresults是函数名，x和y是输⼊参数，没有输出参数。
数据类型 matlab共有6中基本数据类型，分别是数值类型、逻辑类型、字符串、函数句柄、结构体和单元数组。
这⾥我们简单地介绍前四种。
1. 数值类型 基本的数值类型包括整数类型和浮点数类型，额外的数值类型还有复数类型、⽆穷量（Inf）和⾮数值量（NaN)，后⾯的两种算是matlab的特⾊类型，当然要记录⼀下啦！复数类型 复数包括实部和虚部两部分，matlab中默认使⽤i和j作为复数的虚部标志。
创建复数时，可以直接利⽤复数形式进⾏输⼊或使⽤函数complex。
z=1+2i %利⽤复数形式进⾏输⼊ 输出结果： z=complex(2,3) %利⽤complex（x,y）函数进⾏输⼊ 输出结果： z=complex（2）%利⽤complex（x），如果x是实数，z=x+0i，如果x是复数，z=x 输出结果：⽆穷量（Inf）和⾮数值量（NaN) Matlab中使⽤Inf和-Inf分别代表正⽆穷量和负⽆穷量，NaN代表⾮数值量。
计算电磁学之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便凭借其特有的优势得以迅速发展。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
3. 坐标轴的调整 实现坐标系的调整的命令是 axis 函数。 调用格式为： axis([xmin,xmax,ymin,ymax,zmin,zmax])
坐标的最小值（ xmin,ymin,zmin）必须小于相应的最大值 （ xmax,ymax,zmax），否则会出错。
自动坐标系与用 axis 函数调整后的坐标系的比较。
用plot绘图在确定自变量的取值间隔时，一般采用平均间隔，有时会因 某处 间距太大，而不能反映出函数的变化情况。fplot是绘制函数 y=f(x)图形的专用命令,它的数据点是自适应产生的，对那些导数变化较 大的函数，用 fplot 函数绘出的曲线比等分取点所画出的曲线更加接近 真实。 fplot 函数命令的调用格式为： [X,Y]=fplot(‘fun’,lims） fun：函数名字符串； lims：定义 x 的取值区间，lims=[xmin,xmax]；
例
3.1416
3.1416 3.14159265358979 3.1416e+000
format long e
format short g format long g format compact format loose
长格式e方式
短格式g方式 长格式g方式 压缩格式 自由格式
3.141592653589793e+000
用 mesh 命令绘制上例中的网格曲面。
[X,Y] = meshgrid(-2:.2:2, -2:.2:2); Z = X .* exp(-X.^2 - Y.^2); mesh(Z)
与 mesh 相关的 另外两个函数是 meshc 和 meshz ，它们的调用形 式与 mesh 相同。
c. 三维表面命令 surf 函数 surf 可实现对网格曲面片进行着色，将网格曲面转化 为实曲面。surf 命令的调用格式与 mesh 相同。
Matlab 中的数默认是双精度实数，表示方法同 C 语言 3, -9, 0.4, 1.603e-12, 3.23e+20
浮点运算的相对误差为 eps 浮点数表示范围为：10-308 ～ 10308
输出格式
Matlab 以双精度执行所有的运算，运算结果可以在 屏幕上输出，同时赋给指定变量；若无指定变量，则系 统会自动将结果赋给变量 “ans”
四、常用命令
• 1. input函数：用于向计算机输入一个参数。
• 调用格式： A=input(提示信息，选项)； • 注：‘s’选项，则允许用户输入一个字符串。 • 例如想输入一个人的姓名，可采用命令 • xm=input('What''s your name:','s') •
• • • • • •
常见矩阵生成函数
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 为主对角线的对角矩阵 提取一个矩阵的下三角部分 提取一个矩阵的上三角部分
清除当前工作空间中的变量
clear
清除当前工作空间中的所有变量 清除指定的变量
clear A x
三、建立矩阵的函数
• 常用函数有：
• • • • • • eye(size(A)) 产生与A矩阵同阶的单位矩阵 zeros(m,n) 产生0矩阵 ones(m,n) 产生幺矩阵 rand (m,n) 产生随机元素的矩阵 Size(a) 返回包含两个元素的向量。 Length(a) 返回向量的长度。
Matlab图形可视化的几个命令
A、二维平面图形与坐标系
1. 几个基本的绘图命令
a. 线性坐标曲线 plot 函数命令 plot 是 MATLAB 二维曲线绘图中最简单、最重 要、使用最广泛的一个线性绘图函数。它可以生成线段、 曲线和参数方程曲线的函数图形。 命令格式: plot(X,Y) plot(x1,y1,x2,y2,…)：综合调用方式
• 二维函数曲线专用命令
ezplot
2. 图形窗口的分割
有时需要在一个图形窗口中显示几幅图，以便对几个函数进行直观、 便捷的比较。由于每个绘图命令在绘制数据图像时都会将已有图形 覆盖掉，而用 hold 命令不能实现同时显示几个不同坐标尺寸下的图 形，用 figure 命令再创窗口又很难同时比较由不同的数据绘得的图 像。
FDTD中的MATlAB编程基础
一、系统预定义变量
pi ： 圆周率 ，其值为 imag(log(-1)) inf，Inf ：无穷大
nan，NaN ：Not-a-Number，一个不定值，如 0/0 eps ：浮点运算相对精度
i，j ：虚部单位，即
1
应尽量避免给系统预定义变量重新赋值！
s 字符串可以是三种类型的符号之一，也可以是线型与颜色和定点标记 与颜色的组合； 如果没有 s 参数，plot 将使用缺省设置（实线，前七种颜色顺序着色） 绘制曲线； 在当前坐标系中绘图时，每调入一次绘图函数，MATLAB将擦掉坐标 系中已有的图形对象。可以用 hold on 命令在一个坐标系中增加新的图 形对象。注意MATLAB会根据新图形的大小，重新改变坐标系的比例。
用不同的线型和标注来绘制两条曲线。
t1=0:0.1:2*pi; t2=0:0.1:6; y1=sin(t1); y2=sqrt(t2); plot(t1,y1,':hb',t2,y2,'--g')
线型和颜色 plot 函数可以设置曲线的线段类型、定点标记和线段颜色。
常用的线段、颜色与定点标记参数
• 二维函数曲线专用命令 fplot
3.1416 3.14159265358979
format + / format bank / format rat / format hex (详情查看联机帮助)
二、变量的读取
将数据文件中的变量载入当前工作空间
load mydata 载入数据文件中的所有变量
load mydata A x 从数据文件中提取指定变量
subplot(2,1,1) t=0:0.1:4*pi; y=sin(t); plot(t,y) subplot(2,1,2) t=0:0.1:4*pi; y=sin(t); plot(t,y) axis([0,max(t),min(y),max(y)])
B、三维绘图 1.
三维曲线绘图命令 三维函数 plot3主要用来表现单参数的三维曲线，与二维绘图函数 plot 相比，只 多了第三维数据。 其调用格式为： plot3(X1,Y1,Z1,s1,X2,Y2,Z2,s2,…) 参数的含义如下： Xn、Yn、Zn：第一到三维数据，是尺寸相等的向量/矩阵； s、s1、s2：是字符串，用来设置线型、颜色、数据点标记。
实现在同一个窗口中同时显示多个图像的命令subplot。 使用格式为： subplot(m，n，i) 其含义为 ：把图形窗口分割为 m 行 n 列子窗口，然后选 定第 i 个窗口为当前窗口。
subplot 命令不仅用于二维图形，对三维图形一样适用。其本质是将 figure 窗口分为几个区域，再在每个区域内分别绘图。
Matlab 中数的输出格式可以通过 format 命令指定
format 只改变变量的输出格式， 但不会影响变量的值！
各种 format short format long format short e
解释
短格式（缺省显示格式），同short
短格式（缺省显示格式），只显示5位 长格式，双精度数15位，单精度数7位 短格式e方式（科学计数格式）
• 2． pause函数：暂停程序的执行。
• 调用格式： pause(延迟秒数) • 注：如果省略延迟时间，直接使用pause，则将暂停程序，直到 用户按任一键后程序继续执行。
3.Drawnow函数： 将还未处理完的图像实时的显示出来。当代码执行时间
长，需要反复执行plot时，Matlab程序不会马上把图像画到figure上，这时，要 想实时看到图像的每一步变化情况，需要使用这个语句。
rand(m,n)
randn(m,n)
产生 0～1 间均匀分布的随机矩阵 m=n 时简写为 rand(n)
产生均值为0，方差为1的标准正态分布随机矩阵 m=n 时简写为 randn(n)
其它特殊矩阵生成函数：magic、hilb、pascal
可利用冒号提取矩阵 的整行或整列。 例：>> A(1, :) >> A(:, 1:3) >> A(:, :)
【例2】 求一元二次方程a2 +bx+c=0的根。
a=input('a=?'); b=input('b=?'); c=input('c=?'); d=b*b-4*a*c; x=[(-b+sqrt(d))/(2*a),(-b-sqrt(d))/(2*a)] 将该程序以aa.m文件存盘，然后运行aa.m文件。
用命令 plot(x,y)绘制函数 y=cos(x)在两个周期内的图形。
x=0:0.01:2*pi; y=cos(x); plot(x,y)
在同一图形窗口中用命令 plot(x,y)绘出正弦余弦函数的图形。
x=0:0.01:2*pi; y=[sin(x);cos(x)]; plot(x,y)
调用格式：plot(x,y,s) ，s 为类型说明参数，是字符串。