最近做耦合问题又需要用到有限差分，就把这个放在草稿箱里的有限差分整理一遍发出来就当复习了！
 
1.有限差分法的引入与介绍
 
如果用一句话来概括有限差分法，那就是利用用网格节点逼近导数并建立代数方程组来求解各个网格上节点的值。在系列一里面我们知道要把一个无限、连续的问题通过区域的划分变成有限、离散的问题。在有限差分里的区域划分是等距的网格。
 
 
 
有限差分：利用用网格节点逼近导数，并建立有限个未知数的代数方程组来求解各个网格上节点的值。
 
 
上面这句话很显然就产生了两个问题，如何用网格节点逼近导数？（可以用高数里面学的泰勒展示）逼近哪里的导数？（这里面的区别就产生了向前向后和CN）。解决了上述两个问题，也就产生了不同的差分格式。
 
 
 
1.离散区域并给出差分格式
 
 
获得差分格式其实并不是一个很复杂的过程，随便几种差商格式组合、泰勒展开都能得到很多的差分格式，但有限差分最核心的地方并不在于差分格式构造，而是在于分析得到的差分格式靠谱不靠谱。
 
而这个分析格式是否靠谱，这就像兄弟们找对象，对方能不能逼近自己的理想型，这个叫收敛性，现实的遇到的那个人（数值解）和自己的理想型（方程的真解）越接近当然越好；在一起后遇到小矛盾（对应过来是小扰动或者截断、舍入误差）两个人的处理方式和后续的发展，这个叫稳定性，不管多小的矛盾都要闹成大矛盾，这种就是不稳定的，当然得淘汰掉ಠ_ಠ。
 
收敛性这个没有什么好说的，在本系列第一篇文章中有说到。至于稳定性，电脑计算小数时四舍五入造成的误差（舍入误差）一定存在，就要考虑他在逐层求解的时候是否会变得越来越大。常用的分析方法有矩阵法和fourier方法。
 
相容性是在得到一个点的值后，他和初始定义的值能不能相等，这个是相容性，一般是在第四步代数方程组求解完在验证。
 
以上一起就是有限差分法最核心的理论分析：
 
 
 
2.差分解的唯一性、收敛性计稳定性的讨论
 
 
接下来，这一步就是在前面的离散里面得到了一堆离散的关系表达式，如果这个关系表达式是显式的，而第一层（初边值条件）上的每一个离散点是已知的，这样可以用第一层离散点值来直接表示出第二层上所有的离散点的值。依次这样算下去，所有层上的离散点都能被求出来。当网格足够密的时候，也就是离散点足够多的时候，可以将其近似看作求出了$U$在区域内的解。
 
对隐式差分格式的计算，一般来说是一个三对角矩阵线性代数方程组的求解。这个问题就要从差分理论的战场转移到线性代数求解算法和计算复杂度优化上了：在数学上可以用追赶法和SOR方法来求解优化，在计算机上可以用稀疏矩阵的格式来存储来优化（因为三对角阵的0元素比较多）
 
这就是理论分析过后进行实际求解并且和计算机结合起来的一步：
 
 
 
3.求解代数方程组
 
 
最后一步当然就是通过程序的结果看看他是否和我们当时理论所计算的误差精度、稳定性是否一样啦
 
 
 
4.验证程序结果和理论结果的一致性
 
 
2.几个常见的差分格式
 
2.1 一维问题
 
1.对求解域做网格剖分
 
我们将区间[a,b]划分为等距的n个部分，对应的有n+1个节点，取这些节点为：
 
$x_0,x_1...x_i...x_n$
 
相邻节点间距为$h=\frac{b-a}{n}$，记为步长。
 
2.几种常见的差分格式
 
 
 
差分格式：用离散网格点的组合（这里也就是差分或者差商）代替方程中的偏导数的格式
 
 
在一维情况，取等距模板点$x_i$、$x_{i-1}$、$x_{i+1}$，两点间距都为 $h$，常见的差分格式有：
 
向前差分格式：即两个离散点的差商是用来代替前面那个点的导数
 $u_x=\frac{x_{i+1}-x_i}{h}$
向后差分格式：即两个离散点的差商是用来代替后面那个点的导数
 $u_x=\frac{x_i-x_{i-1}}{h}$
中心差分格式 ：这时候要三个离散点，最边上两个点的差商来代替中间那个点的导数。
 $u_x=\frac{x_{i+1}-x_{i-1}}{2h}$
 这时候也能刻画二阶导数，这就是二阶中心差分格式
 $u_{xx}=\frac{x_{i+1}+x_{i-1}-2x_i}{h^2}$
 
热传导方程中C-N格式（克兰克－尼科尔森方法）必须得拿来介绍一下。因为是在热传导方程中，所研究的对象 $u$ 就是描述一个区域内的温度如何随时间变化，所以他含有对时间的偏导数。
 
CN格式的核心思想简单理解就是在时间层 $t=n+1/2$处的导数用中心差分格式，在$xj$用二阶中心差分格式。他的思路来源是：既然(j,n+1)(j,n)两点的差商可以代替t = n和t = n+1处的导数，那么能不能用这两个点的差商来代替中间点$t=n+1/2$的导数呢？当然是可以的，并且代替后，因为$t=n+1/2$这个时刻虽然在现实意义本身是存在的，但是他并没有出现在我们离散的网格点里面，我们就不能把它当作存在来处理，所以就用中心差分格式把这个时刻的导数用(j,n+1)(j,n)两点的差商来表示，这就是CN格式的由来。
 
同理还可以代替其他点比如$t=n+1/12$处的导数，代替不同点的导数截断误差也会有差别。
 
 
注1：在二维情况下，根据不同的目标问题，还可以取更多的模板点来构造差分格式以获得更好的数值解，例如椭圆型问题中取六个点的六点对称格式、激波管问题中引进人工粘性项的Lax-Wendroff格式等。
 （待更新）
 
综上，所以可以看出来有限差分的步骤并不难，理论的基础相比有限元法要简单很多，以后的文章再利用具体的代表算例来求解。

 
原标题：本科 研究生《时域有限差分法(基于MATLAB)》
 
"
 

 

 
时域有限差分法(FDTD)是一种时域电磁算法，参数设置灵活，对复杂介质的模拟具有先天优势。该算法自1966年由Yee提出以来发展迅速，获得了广泛应用。FDTD方法将电场和磁场分别在空间和时间上交错采样，将麦克斯韦方程组转化为差分方程，表述十分简单，容易理解。但是在具体编程实现时，涉及多个维度、多种场量，处理起来非常繁琐。
 
本书采用MATLAB语言编程实现FDTD，充分利用MATLAB向量化编程的特点，将复杂的运算在尽量短的代码内完成，大大简化了编程。对于初学者，这是一本很好的入门教材; 对于已经具有一定基础的学者，本书也能够给予一定的参考。
 
本书共10章。
 
第1章主要介绍MATLAB的一些编程技巧。市面上已经有大量关于MATLAB的教材，本书不再详述，而是只挑选一些与本书的代码密切相关的内容进行讲解， 如向量化运算、维度拓展等。
 
第2章介绍电磁波基础理论。该章中的一些内容可作为理论基础，应用到后续章节的算法中; 另一些内容求出了典型问题的解析解，其结果可以作为验证算法的依据。
 
第3章介绍了FDTD的网格划分方法以及时间推进方法，讨论了空间和时间步长对仿真的影响。
 
第4章和第5章分别介绍了Mur吸收边界条件和完全匹配层(PML)吸收边界条件。研究开域问题时，由于计算机内存有限，只能计算有限区域的场，因此必须在截断边界处加以处理，吸收外向行波，以模拟无限大空间。
 
第6章介绍各种激励源的特点。
 
第7章介绍连接边界条件使用，通过连接边界入射波引入到总场区。
 
第8章介绍远场外推方法。在很多问题中，人们更关心的是远场情况，如计算RCS。FDTD的优势之一就是通过脉冲响应的傅立叶变换得到整个频域上的解，因此这一章以瞬态场的外推为主。
 
第9章介绍了色散介质的处理方法，包括递归卷积法、Z变换方法以及辅助微分方程法。
 
第10章介绍了周期边界的处理方法，以垂直入射情形为主，通过光子晶体、频率选择表面等算例对算法进行了验证。
 
责任编辑：

 
霍文晓
 
(青岛农业大学理学与信息科学学院，山东 青岛 266109)
 
【摘要】为了提高教学效果，在教学过程中引入仿真教学方式。结果表明，利用仿真软件进行演示，能够形象的反映有限差分法的解题过程，并得到电位分布图，加强了学生对抽象理论的理解。
 
教育期刊网 http://www.jyqkw.com
 
关键词有限差分法；MATLAB；仿真分析
 
在电磁场理论中，已知场量在场域边界上的值，求场域中的场分布称为边值问题。通常将静态场边值问题的求解简化成：在一定边界条件下对位函数的泊松方程或拉普拉斯方程的求解[1]。在电磁场与电磁波课程教学中，边值问题的求解既是重点又是难点。教材中主要讲了三种方法：镜像法、分离变量法和有限差分法。随着计算机技术的发展和模拟软件的进步，有限差分法得到了迅速发展和广泛应用。因此，为了与实际接轨，在课堂讲授中，我们将有限差分法作为边值问题这部分的重点内容。并采用理论讲解与模拟演示的教学方法，同时提高学生对理论知识的理解和应用能力。
 
1　有限差分法的原理
 
有限差分法的基本思想是将场域划分成网格，把求解场域内连续的场分布用求解网格节点上的离散的数值解来代替，即用网格节点的差分方程近似代替场域内的偏微分方程来求解。
 
1.1　位函数的差分方程
 
在一个边界为L的二维无源区域S内，电位函数φ(x,y)满足拉普拉斯方程和边界条件为：
 

 
通常将场域分成足够小的正方形网格,网格线之间的距离为h,节点(xi,yi)处的电位φi,j可由其周围直接相邻的四个节点的电位表示，即二维拉普拉斯方程的差分形式。
 

 
同时将边界条件进行离散化，成为边界节点上的已知数值。在这些已知节点条件下，求解各节点的差分方程，整个区域中的节点上电位值即可求出。
 
1.2　差分方程的求解方法
 
在求解实际问题时，为了达到足够的精度，需将网格划分的充分细，节点的个数很多，建立的差分方程数量大，一一求解工作量大。因此如果节点数量较多，通常使用迭代法。
 
1.2.1　简单迭代法
 

 
进行反复迭代(k=0,1,2,…)。若当第N次迭代以后，所有内节点的相邻两次迭代值之间的最大误差不超过允许范围，则终止迭代，并将最后一次迭代的结果作为内节点上电位的最终数值解。
 
1.2.2　超松弛迭代法
 
简单迭代法的收敛速度较慢，为了加快收敛速度，实际中常采用超松弛迭代法[3-4]。迭代公式为
 

 
式中：α称为加速收敛因子，其取值范围是1≤α＜2，当α≥2时，迭代过程将不收敛。
 
加速收敛因子α有一个最佳取值问题，但随具体问题而异。对于第一类边值问题，若求解区域为矩形场域，且由正方形网格分割(每边结点数分别为m和n)，则最佳收敛因子α可按下式计算。
 

 
2　模拟演示
 
以一个正方形截面的无限长金属盒为例，演示用MATLAB对有限差分法的仿真。盒子的两侧及底的电位为零，顶部电位为100V，求盒内的电位分布。
 
由于场域内不存在电荷，其电位分布必满足拉普拉斯方程。将正方形区域划分成10×10的网格。
 
2.1　简单迭代法仿真结果
 
为简单起见，将场域内部节点上的电位初始值全部取为零，利用式(3)求出各内部节点电位值的一次解φij(1)。原来零次解中的各节点电位值将被一次解中的相应电位值所取代。重复上述步骤，令每一个内部节点上的第k+1次解电位值等于该节点周围四个相邻节点(或边界点)第k次解电位值的算术平均值。直到相邻两次的迭代值相差不超过设定的误差范围(1e-6)后，退出迭代。
 
仿真结果：迭代次数为150次，迭代后各节点的电位如图1所示。
 

 
根据节点电位画出电位分布曲线如图2所示。
 

 
由得出的数值解可以看出，金属盒内点电位分布是越靠中间电位越高，越靠近金属盒顶部电位越高，这是由于金属盒底部和两边的电位都为零，而顶部最高。由此表明此方法计算出的电位值，符合金属盒内的电位分布情况。
 
2.2　超松弛迭代法仿真结果
 
超松弛迭代法的关键在于收敛因此的取值，合适的收敛因子可以减少迭代次数。
 
由式(5)可知，当网格节数为10×10时，收敛因子取值为1.56迭代次数最少。通过MATLAB仿真，可得收敛因子与迭代次数的关系，如表1所示。
 

 
从结果可知，收敛因子选择1.56迭代次数最少。与理论计算所得结果相等。
 
通过简单迭代法和超松弛迭代法对比发现，两种方法求出的数值解相同，得到的等电位线分布一样，但超松弛迭代法比简单迭代法收敛速度更快，迭代次数更少，计算时间更短。
 
3　结论
 
利用MATLAB软件对有限差分法分析边值问题进行了仿真分析，不但让学生对有限差分法这种分析方法有了直接的接触和了解，同时对边值问题的处理方法和结果有了更深的认识。通过这种教学方法，使抽象的电磁场问题变得直观、形象，既可以活跃课堂气氛，同时可以加深学生对理论知识的理解，并进一步为将来接触时域的有限差分法打下了良好的基础。对培养学生学习的主动性和积极性有着重要的作用。
 
教育期刊网 http://www.jyqkw.com
 
参考文献
 
［1］谢处方,饶克谨.电磁场与电磁波[M]．4版.北京:高等教育出版社,2006:128-165.
 
［2］王洁,陈超波.基于MATLAB的静态场边值问题有限差分法的研究[J].微计算机应用,2010(03).
 
［3］何红雨.电磁场数值计算法与MATLAB实现[M].武汉:华中科技大学出版社,2004:4210.
 
［4］赵德奎,刘勇.MATLAB在有限差分法数值计算中的应用[J].四川理工学院学报:自然科学版,2005,18(4):61-64.
 
［责任编辑：汤静］

 
章 MATLAB编程技巧
 
1.1 MATLAB基本操作
 
1.2 向量化编程方法
 
1.3 数组的自动扩展
 
1.4 计算结果可视化
 
1.5 MATLAB编程原则
 
复习思考题
 
第2章 电磁学基础
 
2.1 场论
 
2.1.1 场的基本概念
 
2.1.2 矢量微分算子
 
2.1.3 其他坐标系下矢量微分算子表达形式
 
2.1.4 时谐场
 
2.2 麦克斯韦方程组
 
2.2.1 麦克斯韦方程组的形式
 
2.2.2 本构关系
 
2.2.3 无源区域电磁场波动方程
 
2.2.4 平面电磁波的传播特性
 
2.2.5 电磁波的极化
 
2.3 电磁辐射
 
2.3.1 格林函数
 
2.3.2 位函数
 
2.4 媒质中电磁波的传播
 
2.4.1 导体中电磁波的传播
 
2.4.2 一般媒质中电磁波的传播
 
2.5 电磁波的反射和透射
 
2.5.1 电磁场边界条件
 
2.5.2 相位匹配条件
 
2.5.3 反射系数和透射系数
 
2.5.4 垂直入射时的匹配条件
 
2.6 电磁散射
 
2.6.1 无限长导体圆柱
 
2.6.2 无限长介质圆柱
 
2.6.3 理想导体球
 
复习思考题
 
第3章 FDTD差分格式
 
3.1 麦克斯韦方程组的改写
 
3.2 Yee元胞形式和时间推进
 
3.2.1 一维情形
 
3.2.2 损耗的处理
 
3.2.3 二维情形
 
3.2.4 三维情形
 
3.3 空间步长和时间步长稳定性
 
3.4 Courant-Friedrichs-Lewy稳定性条件
 
3.5 FDTD差分的数值色散和各向异性特性
 
3.6 指数差分格式
 
复习思考题
 
第4章 Mur吸收边界条件
 
4.1 Engqist-Majda吸收边界条件
 
4.1.1 一维情形
 
4.1.2 二维情形
 
4.1.3 三维情形
 
4.2 一阶和二阶近似
 
4.2.1 一阶近似
 
4.2.2 二阶近似
 
4.3 FDTD差分形式
 
4.3.1 一维情形
 
4.3.2 二维情形
 
4.3.3 三维情形
 
4.4 吸收效果评估
 
复习思考题
 
第5章 完全匹配层
 
5.1 完全匹配层电磁参数
 
5.1.1 Berenger完全匹配层
 
5.1.2 各向异性完全匹配层
 
5.1.3 坐标伸缩完全匹配层
 
5.1.4 三种完全匹配层的等价证明
 
5.2 Berenger完全匹配层FDTD实现
 
5.2.1 一维情形
 
5.2.2 二维情形
 
5.2.3 三维情形
 
5.3 各向异性完全匹配层FDTD实现
 
5.3.1 一维情形
 
5.3.2 二维情形
 
5.3.3 三维情形
 
5.4 坐标伸缩完全匹配层FDTD实现
 
5.4.1 一维情形
 
5.4.2 二维情形
 
5.4.3 三维情形
 
5.5 几种PML吸收边界的比较
 
复习思考题
 
第6章 FDTD常用激励源
 
6.1 几种随时间变化的源
 
6.1.1 时谐场源
 
6.1.2 脉冲源
 
6.2 时谐场振幅和相位的提取
 
复习思考题
 
第7章 连接边界条件
 
7.1 总场和散射场边界入射波的引入
 
7.1.1 一维总场边界条件
 
7.1.2 二维总场边界条件
 
7.1.3 三维总场边界条件
 
7.2 入射波加入的FDTD实现方法
 
7.2.1 入射波加入一维总场边界
 
7.2.2 入射波加入二维总场边界
 
7.2.3 入射波加入三维总场边界
 
复习思考题
 
第8章 远场外推
 
8.1 远场外推等效原理
 
8.2 二维时谐场的外推
 
8.3 三维时谐场的外推
 
8.4 二维瞬态场的外推
 
8.4.1 基本公式
 
8.4.2 外推的MATLAB实现
 
8.5 三维瞬态场的外推
 
8.5.1 基本公式
 
8.5.2 外推的MATLAB实现
 
复习思考题
 
第9章 色散介质
 
9.1 色散介质基本模型
 
9.1.1 频域模型
 
9.1.2 时域模型
 
9.2 递归卷积法(RC)
 
9.2.1 基本公式
 
9.2.2 Debye介质
 
9.2.3 Drude介质
 
9.2.4 Lorentz介质
 
9.3 Z变换方法
 
9.3.1 Z变换形式
 
9.3.2 Debye介质
 
9.3.3 Drude介质
 
9.3.4 Lorentz介质
 
9.4 辅助微分方程法(ADE)
 
9.4.1 Debye介质
 
9.4.2 Drude介质
 
9.4.3 Lorentz介质
 
9.5 介质板反射系数计算
 
9.6 等离子体光子晶体的仿真
 
9.7 双负介质的仿真
 
复习思考题
 
0章 周期边界条件
 
10.1 二维周期结构
 
10.1.1 二维周期边界条件
 
10.1.2 二维周期结构TM波模拟
 
10.1.3 二维光子晶体仿真
 
10.1.4 二维等离子体光子晶体仿真
 
10.2 频率选择表面
 
复习思考题
 
附录 变量名命名规范
 
参考文献显示全部>>隐藏全部>>