有 限 差 分 法
流体运动的控制方程多为偏微分方程，在复杂的情况下不存在解析解。但是对于一些简单的情况存在解析解，偏微分方程的解析解可用精确的数学表达式表示，该表达式给出了因变量在整个定义域中的连续变化状况。有限差分法（Finite Difference Method ,FDM ）是数值计算中比较经典的方法，由于其计算格式直观且计算简便，因此被广泛地应用在计算流体力学中。有限差分法首先将求解区域划分为差分网格，变量信息存储在网格节点上，然后将偏微分方程的导数用差商代替，代入微分方程的边界条件，推导出关于网格节点变量的代数方程组，通过求解代数方程组，获得偏微分方程的近似解。偏微分方程被包含离散点未知量的代数方程所替代，这个代数方程能求出离散节点处的变量，这种离散方法叫做有限差分法。
2.1
有 限 差 分 逼 近
2.1.1 有限差分网格
由于有限差分法求解的是网格节点上的未知量值，因此首先介绍有限差分网格。图2.1 – 1是x-y 平面上的矩形差分网格示意图。在x 轴方向的网格间距为△x ，在y 轴方向的网格间距为△y ，网格的交点称为节点，计算变量定义在网格节点上。称△x 和△y 为空间步长，△x 一般不等于△y ，且△x 和△y 也可以不为常数。取各方向等距离的网格，可以大大简化数学模型推导过程，并且经常会取得更加精确的数值解。本章作为计算流体力学入门知识，假设沿坐标轴的各个方向网格间距分别相等，但是并不要求各方向的网格间距一致。例如假设△x 和△y 是定值，但是不要求△x 等于△y 。
在图2.1 - 1中，网格节点在x 方向用i 表示，在y 方向用j 表示。因此，假如（i ，j ）是点P 在图2.1 – 1中的坐标，那么，点P 右边的第一个点的就可以用（i+1，j ）表示；在P 左边的第一个点的就可以用（i —1，j ）表示；点P 上边的第一个点的就可以用（i ，j+1）表示；点P 下边的第一个点的就可以用（i ，j —1）表示。
2.1.2 几种差分近似
将微分方程化为有限差分方程时，最普遍的形式是基于泰勒展开（Taylor expantion ）的。如图2.1 – 2 所示，假如,i j u 表示点（i ，j ）的待求量，那么，点（i+1，j ）的未知量1,i j u +就可以用基于点（i ，j ）的Taylor 展开表示
2
21,,,,2()()()2
i j
i j i j i j x u u u u x x x +???=+?+??
3
3,3()()6
i j
x u x ??++???? （2.1 - 1） 式（2.1 - 1）是1,i j u +的数学精确表达式，假如数值是有限的而且系列是收敛的，并且△x 趋近于零。
通过下面的例子来进一步说明泰勒展开及其计算精度。考虑关于自变量x 的连续方程f(x)，其中所有的微分都针对x 。如在点x 处的函数值f(x)已知，那么，f(x+△x)值可以通过
Taylor 展开从x 处的信息知道，即
222()()()()2!
n
n n x x f f f f x x f x x x x x n ?????+?=+?++???++?????? （2.1 - 2）
式（2.1 - 2）的意义如图2.1 – 2所示。假设知
道f 在x （图2.1 – 2中的点1）处的值，想利用式（2.1 - 2）求解f 在x+△x(图2.1 – 2中的点2)处的值。检查式（2.1 – 2）的右侧，可以看见第一项f(x)不能作为对f(x+△x)的良好近似，除非函数f(x)为常数。一个重要的精度改进是利用点1处的曲线斜率，这就是式（2.1 - 2）的第二项
f
x x
???的作用。为了取得f(x+△x)处的更好的近似，加入了第三项
2
2
2()
2
x f x ???，这是点1和点2之间的曲线曲率。
一般来说，为了取得更高的精度，必须加入更高的项。实际上，当式（2.1 - 2）右端有无穷多项高阶项时，它就变成了f(x+△x)的精确表示。
1. 一阶前差分
下面结合式（2.1 - 1）讨论有限差分的形式，解式（2.1 - 1）的,(
)i j u
x
??，可以得到 2
231,,,,,23()()()()26
i j i j i j i j i j u u x x u u u x x x x +-?????=--+??????? （2.1 - 3） 式中：等式左端的微分项是在点（i ，j ）处取值；右端的第一项，即
1,,i j i j
u u x
+-?是微分项的
有限差分形式；右端剩下的部分是截断误差。所以式（2.1 - 3）可近似写为
1,,,()i j i j i j u u u
x x
+-?≈?? （2.1 - 4） 将式（2.1 - 4）与式（2.1 - 3）对比，可以看出式（2.1 - 4）具有截断误差，并且截断误
差的最低阶是△x 的一次方。因此，称式（2.1 - 4）的有限差分形式为一阶精度。可以将式（2.1 – 4）更加精确的表示为
1,,,()()i j i j i j u u u
O x x x
+-?=+??? （2.1 - 5） 式（2.1 - 5）中，符号()O x ?是表示截断误差，代表“△x 的阶”。在式（2.1 - 5）中，截断误差的阶数被明确的表示出来。从图2.1 – 1可以看出，式（2.1 - 5）中的有限差分格式使用到了点（i ，j ）及其右边的信息，即用到了,i j u 和1,i j u +；没有用到点（i ，j ）左边的信息，所以式（2.1 - 5）中的有限差分格式叫做前差分。因此定义式（2.1 - 5）中,(
)i j u
x
??的一阶有
限差分形式为一阶前差分。图2.1 – 3中将x 方向一阶前差分用到的节点称为差分模块，点周围标明加号或减号是与公式相对应的。
2．一阶后差分
现在写出用,i j u 表示1,i j u - 的Taylor 展开形式：
23
231,,,,,23()()()()()26
i j
i j i j i j i j x x u u u u u x x x x -?????=-?+-+?????? （2.1 - 6）
对于,(
)i j u
x
??，可以得到 ,1,,(
)()i j i j i j u u u
O x x x
--?=+??? （2.1 - 7） 组成式（2.1 - 7）的有限差分的信息来自于点（i ，j ）的左边，也就是说，用到了,i j u 和
1,i j u -，没有用到点（i ，j ）右边的信息，所以等式（2.1 - 7）中的有限差分格式叫做后差分。
而且截断误差中△x 的最低阶数是一阶，所以，称式（2.1- 7）中的有限差分格式为一阶后
差分。相应的差分模块如图2.1 – 4所示。
3．二阶中心差分
在计算流体力学应用中，差分格式仅仅为一阶精度是不够的。为了建立二阶的有限差分格式，可以用式（2.1 - 1）减去式（2.1 - 6），得
3
31,1,,,3()2()2()6
i j i j
i j i j x u u u u x x x +-???-=?++????? （2.1 - 8）
式（2.1 - 8）可以写为
1,1,2,(
)()2i j i j i j u u u
O x x x
+--?=+??? （2.1 - 9） 组成式（2.1 - 9）的有限差分的信息来自于点（i ，j ）的左右两边，也就是说，用到了1,i j u +，
,i j u 和1,i j u -，点（i ，j ）在两个相邻的格点之间。而且，式（2.1 - 8）中的截断误差中△x 的
最低阶数是2()x ?，即二阶。因此，称式（2.1 - 9）的有限差分为二阶中心差分。
同理，y 方向差分的获得使用的是与上述相同的方法，结果也和上面x 方向的等式相似，即
,1,,,1,,1,1
2()(
)()()2i j i j
i j i j i j i j i j u u O y y u u u O y y
y u u O y y +-+--?+??
???-??=+?????
?-+??
???
前差分 （2.1-10）后差分 （2.1-11）
中心差分 （2.1-12）
y 方向相应的差分模块如图2.1 – 6、图2.1 – 7和图2.1 – 8所示。 式（2.1 - 5）、式（2.1 - 7）以及式（2.1 - 9）～式（2.1 - 12）都是对一阶导数项的有限差分形式。
4．二阶导数的中心差分
在流体力学方程中，扩散项或黏性项中存在关于变量的二阶导数。采用下面的Taylor 展开方式可以获得二阶导数的差分形式。
将式（2.1 - 1）和式（2.1 - 6）相加，有
4242
1,1,,,,24()2()()()12
i j i j
i j i j i j x u u u u u x x x +-???+=+?++?????
2,2()i j u
x
??可以表示为 21,,1,2,22
2()()()
i j i j i j i j u u u u
O x x x +--+?=+??? （2.1 - 13） 在式（2.1 - 13）中，右边的第一项是点（i ，j ）在x 方向二阶微分形式的中心差分，从
剩下的项中可以看到中心差分是二阶精度的。x 方向二阶导数的二阶中心差分模块，如图2.1 – 9所示，点周围的加号和减号表示这些点是被加还是被减，例如点（i ，j ）旁边的（—2）表示在形成有限差分形式的时候，这个点被减了两次。
y 方向二阶微分形式可以类似的得出
21,,1,2,22
2()()()
i j i j i j i j u u u u
O y y y +--+?=+??? （2.1 - 14）
式（2.1 - 13）和式（2.1 - 14）是二阶导数的二阶中心差分的例子。
5. 二阶混合导数的差分
对于混合微分形式，例如2/u x y
???，也可以得出适当的有限差分形式。利用式中（2.1 - 1），可以得到
23
234
1,,,,,
23
()() ()()()()()
26
i j i j i j i j i j
x x
u u u u u
x
y y x y x y x y
+
??
?????
=+?+++???????????
（2.1 - 15）利用式（2.1 - 6），可以得到
23
234
1,,,,,
23
()() ()()()()()
26
i j i j i j i j i j
x x
u u u u u
x
y y x y x y x y
-
??
?????
=-?+-+???????????
（2.1 - 16）用式（2.1 - 15）减去式（2.1 - 16），得
3
24
1,1,,,
3
()
()()2()()
6
i j i j i j i j
x
u u u u
x
y y x y x y
+-
?
????
-=?++???
??????
将混合微分
2
,
()
i j
u
x y
?
??
用有限差分表示
1,1,3
24
,,
3
()()
()
()()
212
i j i j
i j i j
u u
x
u u
y y
x y x x y
+-
??
-
?
??
??
=-+???
?????
（2.1 - 17）在式（2.1 - 17）中，右边的第一项含有
u
y
?
?
，计算点为（i+1，j）和（i—1，j）。参考图2.1 – 1，可以看出这两点的
u
y
?
?
都可以用式（2.1 - 12）中的二阶中心差分形式表示，但是首先需要近似格点（i+1，j），然后再近似格点（i—1，j）。具体的说，在式（2.1 - 17）中先将1,
()
i j
u
y+
?
?
表示为
1,11,121,()()2i j i j i j u u u
O y y y
+++-+-?=+??? 然后将1,(
)i j u
y
-??用类似的形式表示为 1,11,121,()()2i j i j i j u u u
O y y y
-+----?=+??? 这样，式（2.1 - 17）变成
21,11,11,11,122
,()(),()4i j i j i j i j i j u u u u u O x y x y x y
+++--+----+???=+???????? （2.1 - 18） 式（2.1 - 18）的截断误差来自于式（2.1 - 17），忽略的最低项是2()O x ?，由于式（2.1
- 12）的中心差分是2()O y ?，因此式（2.1 - 18）的截断误差是22
(),()O x y ??????。式（2.1 -
18）给出了混合微分2,()i j u
x y
???的二阶中心差分，相应的差分模块如图2.1 – 11所示。
6．二阶精度的差分
前面介绍的有限差分形式仅仅是最基本的情况，其他许多差分形式可以用与上面相似的方法表示出来，尤其是可以获得更加精确的三阶精度、四阶精度，甚至更高阶精度的有限差分形式。这些高阶精度的有限差分形式涉及的点比上面
所提到的形式涉及的点多。例如，四阶精度的关于22u
x
??的中
心差分形式是
22,1,,1,2,4
,22
163016()()12()
i j i j i j i j i j i j u u u u u u O x x x ++---+-+-?=+??? （2.1 - 19） 注意到，四阶精度的中心差分形式涉及到五个点，而在式（2.1 - 13）中，2,2()i j u
x
??用到
了3个点的信息，仅具有二阶精度。式（2.1 - 19）可以通过重复使用关于点（i+1，j ），(i ，j)，（i —1，j ）的Taylor 展开得到，强调可以从几乎是无穷多的有限差分形式得到不断增大的精度。在计算流体力学中，根据问题的特点决定差分格式的精度，一般情况下二阶精度就够用了。
7．边界处的差分
上面讨论的是在区域内部差分的离散格式，如果遇到边界，离散点的布置将受到边壁的限制。
在图2.1 – 12中表示了流场的一部分边界，y 轴垂直于边界。点1处于边界上，点2和
点3距边界的距离是△y 和2△y 。希望构造一个关于边界的u
y
??的有限差分近似。构造前差分比较容易得到
211(
)()u u u
O y y y
-?=+??? （2.1 - 20）
式（2.1 - 20）是一阶精度的。然而如何获得一个二阶精度的格式呢？式（2.1 - 12）的中心差分未被使用，因为它要求边界外的一个点，如图2.1 – 12中的点2ˊ。点2ˊ是计算区域外的一点，一般没有这个点的速度信息。在早期计算流体力学中，许多方法试图通过假设2u ˊ=2u 来绕过这个问题，称之为反射边界条
件。在多数情况下，与等式（2.1 - 20）的前差分比较，
是没有物理意义的，而且精度不够。
可以采用多项式展开的方法，设计边界处的二阶精度有限差分是一种不同于前述的Taylor 展开的方法。假设图2.1 – 12 中，边界上的u 值可以用多项式表示为
2u a by cy =++ （2.1 - 21）
将式（2.1 - 21）用于图2.1 – 12 中的点，可以知道
点1处 y=0
1u a = (2.1 - 22)
点2处y=△y
22()u a b y c y =+?+? (2.1 - 23)
点3处y=2△y
23(2)(2)u a b y c y =+?+? (2.1 - 24)
解式（2.1 - 22）～式（2.1 - 24）得到
123
342u u u b y
-+-=
? （2.1 - 25）
利用式（2.1 - 21），对y 求导得到
2u
b cy y
?=+? （2.1 - 26） 考虑到y=0，从式（2.1 - 26）可以得到
1(
)u
b y
?=? （2.1 - 27） 比较式（2.1 - 25）和式（2.1 - 27），得到
123134(
)2u u u u
y y
-+-?=
?? （2.1 - 28） 式（2.1 - 28）是一个一边的边界处的有限差分形式。因为式（2.1 - 28）中所有信息都
是边界一边的节点的，即仅用到了图2.1 – 4中点1上部的信息，所以叫做“一边”。同时，式（2.1 - 28）用到了多项式的形式，即式（2.1 - 21）。这给出了获得有限差分的另一种取得途径。另外，还必须考虑式（2.1 - 28）的精度。这时，再次使用Taylor 展开。考虑点1处的Taylor 展开，得到
23
23111123()()()()26
y y u u u u y u y y y ???=++++?????? （2.1 - 29）
比较等式（2.1 - 29）和等式（2.1 - 21）。假设等式（2.1 - 21）中的多项式表达和泰勒展开的前三项相同，那么等式（2.1 - 21）中含有3()O y ?。检验等式（2.1 - 28）的分子，在这里1u ，2u 和3u 都可以用等式（2.1 - 21）中的多项式表达。因为等式（2.1 - 21）中含有3()O y ?，所以等式（2.1 - 28）的分子也含有3()O y ?。然而，在写等式（2.1 - 29）的微分形式时除了
y ?，所以式（2.1 - 21）应该含有2()O y ?。因此，可以从等式（2.1 - 28）得到
2123134(
)()2u u u u O y y y
-+-?=+??? （2.1 - 30） 式（2.1 - 30）就是在边界上的二阶精度的差分形式。
式（2.1 - 21）和式（2.1 - 30）都被称为单侧差分，因为它们都表示依靠单侧点信息的微分。而且，这些形式是通用的，即它们不仅仅局限于边界处，还可以用于网格的内部。仅仅是在讨论边界差分问题的时候导出了单侧差分，单侧差分对于边界问题有重要用途，但是对于整个计算域内部的问题，它仅仅是多提供了一个解决方法。式（2.1 - 30）表示了二阶精度的单侧有限差分，许多其他的更高阶精度的单侧有限差分可以通过用更多点的信息得到。
2.2 差 分 方 程
在2.1节中讨论了偏微分的有限差分形式。偏微分方程一般含有多个导数项，如对流扩散方程中，非定常项为关于时间的一阶导数项，对流项为空间一阶导数项，扩散项为空间二阶导数项。当偏微分方程中的所有导数项都用有限差分表示之后，得到的数学公式称为差分方程。计算流体力学的有限差分求解方法的要点，在于使用2.1节中推导出的差分形式代替控制方程的偏微分形式，得到关于网格节点上变量的有限差分方程。本节介绍差分方程的求解过程。
2.2.1 差分格式的构造
以一维非恒定热扩散方程为例。
22T T
t x
α??=?? (2.2 - 1)
式（2.2 - 1）中热扩散系数α为常数。将式（2.2
- 1）中的偏微分用有限差分表示，式中有两个相互独立的自变量x 和t ，计算网格如图2.2 – 1 所示。此处，x 方向离散节点编号用i 表示，时间层编号用n 表示。将式（2.2 - 1）中的时间变化项用前差分形式表示：
122()()2
n n n i i i i T T T
T t t t t +-???=-+?????? （2.2 - 2） 时间差分的截断误差为一阶。将式（2.2 - 1）中的二阶空间导数项用中心差分近似：
22411
2242()()()()12
n n n
n n i i i i i T T T x T T x x x +--+???=-+?????? （2.2 – 3）
此时，截断误差和式（2.1 - 13）相同，将式（2.2 - 1）写为
220T T
t x
α??-=?? （2.2 - 4） 将式（2.2 - 2）和式（2.2 - 3）代入式（2.2 - 4），得到
1224112
24(2)()()()0()212n n n n n
n n i i i i i i i T T T T T x T t T t x t x αα++-??--+????-+-++???=????????
（2.2 - 5） 式中方括号内给出了差分方程的截断误差，略去截断误差后得到
1112
(2)()n n n n n
i i i i i T T T T T t x α++---+=?? （2.2 - 6） 式（2.2 - 6）是式（2.2 - 1）的偏微分方程的差分形式，其截断误差是2(,)O t x ??。 有限差分方程与偏微分方程不同，它是代数方程，是在图2.2 – 1 的控制域中的每个节
点上写出来，组成一个代数方程组，通过求解代数方程组得到在所有网格点上的离散解，如
n i T ，1n i T +，1n i T +，11n i T ++，2
n i
T +等。从原则上说，仅仅可以求出T 的偏微分方程的在截断误差之内的近似解。当网格节点数趋于无限时，时间步长和空间步长都趋近于零，差分方程的
截断误差也趋近于零，因此有限差分方程的解就接近于偏微分方程的解。关于差分方程的稳定性和收敛性将在后面讨论。
对于与时间相关的非恒定流动，计算格式主要分为显式格式和隐式格式两大类。下面首先介绍关于时间层离散的显格式。
2.2.2 显式差分格式
结合一维热传导方程式进行讨论。
22T T
t x
α??=?? （2.2 - 7） 将式（2.2 - 7）作为这一节讨论中的模型方程，用这个模型方程就可得到与显式差分方
法和隐式差分方法有关的所有必要点。在2.1节中用前差分来离散
T
t
??，用中心差分来离散22
T
x ??，从而得到微分方程的特殊形式： 1112
(2)()n n n n n
i i i i i T T T T T t x α++---+=?? （2.2 - 8） 整理后可写成
1112
(2)()
n n n n n
i i i i i t T T T T T x α
++-?=+-+? （2.2 - 9） 式（2.2 - 7）是抛物型的偏微分方程，
参考图2.2 – 2 中的有限差分网格，假定n 时层所有网格点上的温度值T 已知，时间推进是指在n+1时层所有网格点上的T 值都可以从n 时层的已知值算出。计算完毕便可得知n+1时层的值。然后可以通过同样的步骤，用n+1时层的已知值来计算n+2时层所有网格点上的T ，进而通过时间步的推进逐步得到更多的解。从式（2.2 - 9）中可以看到，可以通过一个简单机制来完成时间推进。注
意到式（2.2 - 9）中将时层n 的有关项写在
等式右边，而将时层n+1的有关项写在等式左边。
注意到在时间推进的过程中所有n 时层的变量T 都是已知的，而n+1时层变量T 则是需要计算的。式（2.2 - 9）中仅含有一个未知数1n i T +，可以直接由n 时层的已知物理量求解
1n i T +，得到一个只含有一个未知数的方程，余下的计算就比较简单了。例如，考虑图2.2 – 3
所示的网格，在x 轴上取7个点，设点2为中心点，可以得到
1223212
(2)()
n n n n n t
T T T T T x α
+?=+-+? （2.2 - 10） 由于式（2.2 - 10）右边的量均为已知值，可以对12n T +进行直接计算。然后，以格点3作为中心点，式（2.2 - 9）可以写成
1334322
(2)()
n n n n n t
T T T T T x α
+?=+-+? （2.2 - 11） 可以由式（2.2 - 11）右边的已知数直接计算13n T +。
同样，将式（2.2 - 9）应用于格点4、5、6，依次可得到14n T +，15n T +，16n T +。 在显式方法中，每个不同的方程都只含有一个未知数，因此可以很轻易的直接求解。这
种显式方法在图2.2 – 3 中用带有虚线框的有限差分模型可进一步描述，在这里模型在n+1时层上仅含有一个未知量。
如图2.2 – 3 所示，求抛物型偏微分方程的推进解要
求给定边界条件约束，这表示T 在左右边界处的值1T 和7T ，根据预先给定的边界条件在每一时层上均为已知值。
2.2.3 隐式差分格式
式（2.2 - 8）并不是唯一可描述式（2.2 - 7）的差分方
程，它只是最初的偏微分方程差分表达中的一种。作为以上关于显式方法讨论的一个反例，将方程的空间差分按照时层n 和时层n+1之间的平均特性的顺序写在方程右边。
也就是说，将方程表达为如下形式：
1111
11112
111()(22)()222()
n n n n n n n n
i i i i i i i i T T T T T T T T t x α++++++--++--++-=?? （2.2 - 12） 式（2.2 - 12）中使用的空间差分格式称为Crank – Nicolson 格式（C – N 格式）。C – N
格式被广泛使用于求解控制方程为双曲型方程的问题中。在计算流体力学中，C – N 格式以及它的修正形式被频繁地使用。仔细观察式（2.2 - 12），未知量 1n i T +不仅根据n 时层的未
知量1n i T +，n i T ，1n i T -来表达，也根据其他时层n+1的未知量11n i T ++，11n i T +-来表达，式（2.2 - 12）代表了一个有3个未知量的方程，有11n i T ++，1n i T +，11n i T +-，因此，在给定格点i 处的方程并不
是独立的，只通过该方程本身无法求解1n i T +的结果。必须在所有内部格点都列出式（2.2 - 10），最后可以得到一个代数方程组，可联立解出i 取所有值时未知量1n i T +的值。
这是一个隐式方法的实例。根据定义，隐式方
法是指在给定时间层，对于所有格点都列出微分方程联立求出所有未知量的解的方法。因为需要求解大型的代数联立方程组，所以隐式方法通常需要对大型矩阵进行处理。到目前为止，可能会觉得隐式方法涉及比前面讨论的显式方法要复杂的运算。与图2.2 – 3所示的简单显式有限差分模型相对比，图2.2 – 4描述了式（2.2 - 12）的隐式模型的草图，清晰地描绘了时层n+1的3个未知量。
更具体地，以图2.2 – 4 所示的有7个格点的网格为例，将式（2.2 - 12）的项重新排列，将未知项移到等号左边，已知项移到等号右边，结果如下：
11111112222
(2)12()()2()2()
n n n n n
n n i i i i i i i tT tT t T T T t T T x x x x αααα+++-++-?????-+?-++=--????????
(2.2 - 13)
为简单起见，将以下的量用A ，B 和i K 表示
2
2()
t
A x α?=
? ； 2
1()
t
B x α?=+
? ； 112
(2)2()
n n n n i i i i i t
K T T T T x α+-?=--
-+?
可以将式（2.2 - 13）写成如下形式：
111
11n n n i i i i AT BT AT K +++-+-+= (2.2 - 14)
注意到式（2.2 - 14）中的i K 包括了已知的时层n 的特性。因此，i K 在式（2.2 - 14）中是已知量，如图2.2 – 4 所示，对格点2和6依次应用式（2.2 - 14）。
格点2： 1232AT BT AT K -+= （2.2 - 15） 为方便起见，省略上标， 1T ，2T ，3T 代表时层n+1的3个值。如前所述，2K 为已知值。另外，由于格点1和格点7处的边界条件约束，方程中的1T 为已知值，因此，式（2.2 - 15）中与已知量1T 有关的项可以移到等号右边，结果如下：
2321BT AT K AT -+=- （2.2 - 16）
将21K AT -记为2K ˊ，2K ˊ已知，式（2.2 - 16）可以写为
23BT AT -+=2K ˊ （2.2 - 17）
格点3： 2343AT BT AT K -+= （2.2 - 18） 格点4： 3454AT BT AT K -+= （2.2 - 19） 格点5： 4565AT BT AT K -+= （2.2 - 20） 格点6： 5676AT BT AT K -+= （2.2 - 21） 在式（2.2 - 21）中，因为格点7在边界上，由边界条件约束可知7T 为已知量，因此，将式（2.2 - 21）整理后可以写成
56676AT BT K AT K -=-=ˊ （2.2 - 22）
其中6K ˊ为已知量。式（2.2 - 17）～式（2.2 - 22）是以2T ，3T ，4T ，5T ，6T 为未知数的五个方程，这个方程组可以写成如下的矩阵形式：
系数矩阵是三对角阵，三对角阵的定义是仅在式（2.2 - 23）中的三条对角线上的元素为非零元素。式（2.2 - 23）的求解可用追赶法。
这个例子的基础更多的涉及到隐式方法而不是显式方法，式（2.2 - 7）是一个线性偏微分方程，由它得到的是线性差分方程。如果偏微分方程是非线性，例如假定式（2.2 - 7）中的温度扩散系数α是一个温度函数
22()T T
T t x
α??=?? (2.2 - 24) 式（2.2 - 24）是非线性偏微分方程。这里参照式（2.2 - 9）将式（2.2 - 24）写成差分形
式如下：
1112
()
(2)()
n n n n n n
i i i i i i t T T T T T T x α++-?=+-+? （2.2 - 25） 式（2.2 – 25）仍然是单未知量1n i T +的线性方程，因为α在时层n 上赋值，即，α=()n i T α，这里n i T 为已知量。另外，如果在式（2.2 - 24）中使用C – N 方法，右边按照时层n 和n+1
的平均值赋值，即()T α可表示为
11()()2n n
i i T T αα+??+??。得到的差分方程由式（2.2 - 13）给出，在方程中出现α的地方均由11
()()2
n n i i T T αα+??+??代替。显然，这个新的差分方程包含了由独立变量组成的项：11()n n i i T T α++，111()n n i i T T α+++和11
1()n n i i T T α++-。换而言之所得的方
程是一个非线性代数方程。因此隐式求解需要一个大型的非线性方程组的同步解。这是隐式求解一个很大的缺点。为了避免这个问题，通常将差分方程近似地离散化。例如，如果在式（2.2 - 24）中α简单地用n 时层的值来赋值而不是用n 和n+1时层的平均值，那么差分方程中就不会出现非线性的代数项。如此便得到理想的差分方程如方程（2.2 - 13）所示，方程中α的值由()n i T α给定。
对于显式方法而言，x ?一旦选定，t ?也不再是一个可任意选择的独立变量，而限制在一个等于或小于一个由稳定准则限定的范围内取值。如果t ?的取值超过了稳定准则要求的范围，时间推进程序就会迅速发展为不稳定，程序会很快因为数值趋向无限大或者对负数取平方根之类的原因被停止。在很多情况下，t ?为了保持稳定性必须取很小的值，这可能会造成在一个给定的时间段内，完成计算需要很长的程序运行时间。而隐式方法却没有这样的稳定要求。
对于大多数的隐式方法而言，为保持稳定而要求的t ?的取值范围要比相应的显式方法大得多。甚至有的隐式方法是无条件稳定的，也就是说t ?取任意值，不管多大，也能得到一个稳定解。因此，对于隐式方法而言，在给定的时间段内需要的时间步数可以比显式方法少相当一部分，所以，对于有些情况，即使隐式方法由于自身的复杂性，而在每个时间步内需要做更多的计算，实际上由于在给定时间段内减少了相当多的时间步，可能最后在电脑上的程序运行时间比显式方法还要短。
由于可以取较大的t ?，截断误差也较大，使用隐式方法追踪准确的瞬态（变量随时间的变化）可能不如显式方法的精度高。但是，当想要求解稳定状态的时间独立解时，关于时间变量和精度就显得不重要了。
时间步进方法在计算流体力学中经常被使用，是求解流场的最直接的方法。但是，许多更为复杂的流体力学计算的应用中，在流场的某些区域要求密集的网格格点，由于网格要求的推进步长很小，会造成电脑运算需要的时间过长。高Re 数黏性流动的计算中，在流场内剧烈改变发生的地方，要求有密集的空间点，是这个问题的一个实例。这使得上面列出的隐式方法的优点更为突出，在一个适合的网格上使用大的时间步长。因为这个原因，隐式方法称为20世纪80年代计算流体力学应用的主要焦点。但是，由于电脑技术的不断提高，即大型并行处理器的相互连接，显式方法有可能再次成为重点。如此大型的处理器可以完成分为数千个网格的流场的同一时刻的显式计算。对于一个给定的问题是选用显式方法还是隐式方法求解有时并不能分得十分清楚。当需要对此进行选择时，必须要精心选择。
2.3 误 差 和 稳 定 分 析
使用有限差分法求解流体流动时，首先要将计算区域划分为差分网格，选取合适的差分格式。需要说明的是，在上面条件下差分方程的解逼近微分方程的解。如果时间步长t ?的增加超过一定值则有的格式会趋向不稳定。这个允许的时间步长最大值可以通过有限差分格式的稳定性分析得到。关于差分方程的解是否收敛，需要讨论差分格式的相容性、收敛性以及稳定性。
2.3.1 差分方程的相容性
相容性阐述的是差分方程能否代表与之相对应的微分方程。当时间步长t ?和空间步长x ?同时趋近于零时，如果差分方程的截断误差趋近于零，此时差分方程趋近于微分方程，称为差分方程与之对应的微分方程相容。如果不管t ?和x ?以怎样的方式趋近，差分方程的截断误差总是趋向于零，称为差分方程与之对应的微分方程是无条件相容的。
以一维线性热扩散方程为例进行讨论
22T T
t x
α??=?? （2.3 - 1） 取Dufort-Frankel 格式，得到
111111
2
()02n n n n n n
i i i i i i T T T T T T t x
α+-+-+---++-=?? （2.3 - 2） 对该格式各项进行Taylor 展开分析
112()()2n n n i i i T T T
O t t t
+--?=+???
2424112224(22)2()2!4!n n
n n n i i i i
i T T x x u u u x x x x αα+-+????=+++?????????? 2424112224
(22)2()2!4!n n n n n i i i i i T T t t u u u x x t t αα+-+????=+++??????????
将上述三式代入式（2.3 - 2）得到截断误差
1111211
2
2()()2n n n n n n n
n i i i i i i i
i T T T T T T T T R t x t x αα+-+-+---++??=---????
22
22
(/)
O t x t x =?+?+?? （2.3 - 3） 可见，如果/t x ??等于常数，当0t ?→，0x ?→时，截断误差n i R 并不趋近于零，因此差分方程（2.3 - 2）与微分方程（2.3 - 1）是不相容的。
如果选用显式差分格式为
1112
(2)()n n n n n
i i i i i T T T T T t x α++---+=?? （2.3 - 4） 通过截断误差分析，得到2()n i R O t x =?+?，当0t ?→，0x ?→时，截断误差n i R 趋近于零，因此差分方程（2.3 - 4）与微分方程（2.3 - 1）是相容的。
2.3.2 差分方程的收敛性
当时间步长和空间步长趋近于零时，差分方程的解趋近于微分方程的解，这种性质称为
差分方程的收敛性。
差分方程的误差由截断误差和离散误差两部分组成。差分方程的截断误差与相容性有关，差分方程满足相容性，并不一定满足收敛性，相容性只是收敛性的必要条件之一。
讨论差分方程的收敛性，还需要分析离散误差。对于定解域上的任意一点（x ，t ），若用n i A 表示微分方程的精确解，n i D 表示差分方程的精确解（假定计算中不引入初始误差、舍入误差），离散误差的定义为
n n n i i i A D ε=- (2.3 - 5)
当时间步长和空间步长趋近于零时，如果离散误差n i ε也趋近于零，则称差分方程是收敛的。直接证明差分方程的收敛性比较困难，根据Lax 等价定理，可以通过稳定性间接证明收敛性。
Lax 等价定理：对于一个适定的线性初值问题及与它相容的差分格式，其收敛性的充分必要条件是此差分方程的稳定性。也就是说，对于一个适定的线性初值问题，相容性加稳定性等于收敛性。其中微分方程的适定性，是指微分方程的解存在并且唯一。
2.3.3 差分方程的稳定性
差分格式计算是按时层逐渐推进的，n i D 为差分方程的精确解，n i u 为差分方程的数值解，两者往往不相等。这是因为初始时层或边界值一般是由实测或推算而确定的，具有初始误差和边界误差，在计算中又有四舍五入而产生的舍入误差，定义取整误差如下：
n n n i i i e u D =- (2.3 - 6)
若在时层推进的计算中，n i e 不是逐渐增大而导致完全偏离差分方程的精确解，则称差分方程是稳定的，否则是不稳定的。
考虑偏微分方程（2.3 - 1），表示的是一维热扩散方程，将式（2.3 - 6）表示的差分方程数值解n n n i i i u D e =+代入差分方程，得到
111111
2
22()
n n n n n n n n n n i i i i i i i i i i D e D e D e D e D e t x α++++--+--+--++=?? （2.3 - 7） 根据定义，n i D 是差分方程的精确解，精确地满足差分方程。因而得到
111
2
2()n n n n n i i i i i D D D D D t x α++---+=
?? （2.3 - 8） 用式（2.3 - 8）减去式（2.3 - 7），得
111
2
2()
n n n n n i i i i i e e e e e t x α++---+=?? （2.3 - 9） 由式（2.3 - 9），知道误差n i e 通常满足差分方程。在从n 向n+1时层推进时，如果n i e 逐渐增大，则差分格式是不稳定的；如果误差n i e 不增大，则差分格式是稳定的。因此差分格式的稳定条件是：
1/1n n i i e e +≤ （2.3 - 10）
对于以式（2.3 - 1）为例的非恒定一维热扩散问题，取整误差可以假设为随机变量，可
以用下面的傅立叶级数表示
()m jk x m m
e x A e =∑ （2.3 - 11）
式中： m k 为波数；j = 2.3 - 11）代表了一个正弦函数和一个余弦函数，因为cos sin m jk x m m e k x j k x =+。对于长度为L 的一维计算域，划分的网格节点数N+1，可以确定这样的网格划分能够考虑的最小波长为
min 2/L N λ= (2.3 – 12)
那么最大波数为
222/2
m N
k L N L ππ=
= （2.3 - 13）
式（2.3 - 11）中m 的上限应该等于N/2，因此式（2.3 - 11）可以改写为
/2
1
()m N jk x m m e x A e ==∑ （2.3 - 14）
认为变量m A 随时间呈指数发展，误差也随着时间呈现指数增大或减小的趋势，因此有
/2
1(,)m N jk x at m e x t e e ==∑ （2.3 - 15）
这里a 是一个常数，式（2.3 - 15）代表了空间和时间的取整误差的最终合理形式。
由于最初的差分方程（2.3 - 2）是线性的，且取整误差满足这个差分 方程，将式（2.3 - 15）代入式（2.3 - 9），级数的每一项的行为都与该级数一致。因此，只需处理级数的其中一项，有
(,)m jk x at m e x t e e = （2.3 - 16）
现在研究误差随着时间步的变化，并由此得出满足式（2.3 - 9）对时间步长t ?的限制。
将式（2.3 - 16）代入式（2.3 - 9），得
()()
()2
2()
m m m m m jk x jk x jk x x jk x jk x x a t t at at at at e e e e e e e e e e t x α+?-?+?--+=?? （2.3 - 17） 从式（2.3 - 17）得到
2
1(2)()
m m jk x jk x
a t t
e e e x α?-???=+
+-? （2.3 - 18） 由公式cos()2
m m jk x jk x
m e e k x ?-?+?=，式（2.3 - 18）可以写成
[]221cos()1()
a t m t
e k x x α??=+
?-? （2.3 - 19） 由三角公式2
1cos()
sin
22
m m k x k x ?-?=式（2.3 - 19）最终成为 2241sin
()2
a t m k x t e x α???=-
? （2.3 - 20） 由式（2.3 - 16），得
1()m m jk x
n a t t a t i jk x
n at i e e e e e e e
++??== (2.3 - 21) 通过式（2.3 - 20）和式（2.3 - 21），可以得到
12241sin 1()2
n a t i m n i e k x t e e x α+???==-≤? （2.3 - 22） 只有满足式（2.3 - 22），差分方程才有稳定解。在式（2.3 - 22）中，定义因数
22
41sin ()2
m k x t G x α??-
≡? 为放大因数，并用G 表示，给式（2.3 - 22）中的不等式赋值，也就是1G ≤，有两种可能的
情况必须保持同步：
（1）
2241sin 1()2m k x t x α??-
≤? 且 224s i n 0()2
m k x t x α??≥? （2.3 - 23） 由于24/()t x α??总是为正，式（2.3 - 23）的条件始终能够满足。
（2）
2241sin 1()2
m k x t x α??-
≥-? 因此 224sin 11()2m k x t x α??-≤? 要满足上述条件，有
2
1
()2
t
x α?≤
? （2.3 - 24） 式（2.3 - 24）给出了差分方程（2.3 - 2）的解稳定的条件。t ?的允许值必须足够小以满足式（2.3 - 24）。只要2/()1/2t x α??≤，误差就不会随着时间步的推进而增大，其数值解也会表现出稳定的形式。另一方面，如果2/()t x α??>1/2，误差就会逐渐增大，最终使数值推进解在电脑上发生溢出。
2.3.4 波动方程的稳定性分析（CFL 条件）
在流体力学数值模拟中，波动方程对数值格式的要求比较严格，通过分析一次波动方程，获得其稳定条件。一次波动方程如下：
0u u c t x
??+=?? （2.3 - 25） 11
2n n i i u u u x x
+--?=?? （2.3 - 26） 将时间导数用一个简单的前差分代替，那么得到的与式（2.3 - 25）相应的差分方程为
111
2n n n n i i i i u u u u c t x
++---=-?? （2.3 - 27） 称式（2.3 - 27）为Euler 显式格式。但是，冯﹒纽曼稳定性分析在式（2.3 - 27）中的应用说
明不论t ?的值为多少，式（2.3 - 27）将趋向不稳定解，因此式（2.3 - 27）成为无条件不稳定。将时间导数用一次差分代替，n 时层的数值用格点i + 1 和i —1的平均值代替，得到
1111
()2n n n i i i u u u u t t
++--+?=?? (2.3 - 28) 将式（2.3 - 28）代入式（2.3 - 27），有
1111122
n n n n n i i i i i
u u u u t u
c x ++-+-+-?=-? （2.3 - 29）
在上面的方程中使用的差分格式为Lax 方法，在这里用式（2.3 - 28）代表时间导数。
Lax 方法是用最初提出该方法的数学家彼得列克斯（Peter Lax ）的名字命名的。如果假定误差具有如下形式(,)m jk t
at m e x t e e
=，在式（2.3 - 29）中使用这种形式，放大因子变成
cos()sin()at m r m e k x jC k x =?-? （2.3 - 30）
式中：/r C c t x =??。为保证稳定要求1at
e ≤，这里对式（2.3 - 30）应用时得出
1r t
C c
x
?=≤? （2.3 - 31） 在式（2.3 - 31）中，r C 称为库朗数（Courant number ）。这个方程说明/t x c ?≤?时式（2.3 - 29）的数值解稳定。式（2.3 - 31）称为CFL （C o urant – Friedrichs - Lowy ）条件，作为一般性的计算流体力学条件，这是双曲型方程一个很重要的稳定判断准则。
有限差分法及其应用 1有限差分法简介 有限差分法（FDM）是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方程将解域划分为差分网格，用有限个网络节点代替连续的求解域。有限差分法通过泰勒级数展开等方法，把控制方程中的导数用网格节点上的函数值得差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 2有限差分法的数学基础 有限差分法的数学基础是用差分代替微分，用差商代替微商而用差商代替微商的意义是用函数在某区域内的平均变化率来代替函数的真是变化率。而根据泰勒级数展开可以看出，用差商代替微商必然会带来阶段误差，相应的用差分方程代替微分方程也会带来误差，因此，在应用有限差分法进行计算的时候，必须注意差分方程的形式，建立方法及由此产生的误差。 3有限差分解题基本步骤 有限差分法的主要解题步骤如下： 1）建立微分方程 根据问题的性质选择计算区域，建立微分方程式，写出初始条件和边界条件。 2）构建差分格式 首先对求解域进行离散化，确定计算节点，选择网格布局，差分形式和步长；然后以有限差分代替无线微分，以差商代替微商，以差分方程代替微分方程及边界条件。 3）求解差分方程 差分方程通常是一组数量较多的线性代数方程，其求解方法主要包括两种:精确法和近似法。其中精确法又称直接发，主要包括矩阵法，高斯消元法及主元素消元法等；近似法又称间接法，以迭代法为主，主要包括直接迭代法，间接迭代法以及超松弛迭代法。4）精度分析和检验 对所得到的数值进行精度与收敛性分析和检验。 4商用有限差分软件简介 商用有限差分软件主要包括FLAC、UDEC/3DEC和PFC程序，其中，FLAC是一个基于显式有限差分法的连续介质程序，主要用来进行土质、岩石和其他材料的三维结构受力特性模拟和塑性流动分析；UDEC/3DEC是针对岩体不连续问题开发，用于模拟非连续介质在静，动态载荷作用下的反应；PFC是利用显式差分算法和离散元理论开发的微、细观力学程序，它是从介质的基本粒子结构的角度考虑介质的基本力学特性，并认为给定介质在不同应力条件下的基本特征主要取决于粒子之间接粗状态的变化，适用于研究粒状集合体的破裂和破裂发展问题，以及颗粒的流动(大位移）问题。
时域有限差分法的Matlab仿真 关键词: Matlab 矩形波导时域有限差分法 摘要：介绍了时域有限差分法的基本原理，并利用Matlab仿真，对矩形波导谐振腔中的电磁场作了模拟和分析。 关键词：时域有限差分法；Matlab；矩形波导；谐振腔 目前，电磁场的时域计算方法越来越引人注目。时域有限差分（Finite Difference Time Domain，FDTD）法［1］作为一种主要的电磁场时域计算方法，最早是在1966年由K. S. Yee提出的。这种方法通过将Maxwell旋度方程转化为有限差分式而直接在时域求解，通过建立时间离散的递进序列，在相互交织的网格空间中交替计算电场和磁场。经过三十多年的发展，这种方法已经广泛应用到各种电磁问题的分析之中。 Matlab作为一种工程仿真工具得到了广泛应用［2］。用于时域有限差分法，可以简化编程，使研究者的研究重心放在FDTD法本身上，而不必在编程上花费过多的时间。 下面将采用FDTD法，利用Matlab仿真来分析矩形波导谐振腔的电磁场，说明了将二者结合起来的优越性。 1FDTD法基本原理 时域有限差分法的主要思想是把Maxwell方程在空间、时间上离散化，用差分方程代替一阶偏微分方程，求解差分方程组，从而得出各网格单元的场值。FDTD 空间网格单元上电场和磁场各分量的分布如图1所示。 电场和磁场被交叉放置，电场分量位于网格单元每条棱的中心，磁场分量位于网格单元每个面的中心，每个磁场（电场）分量都有4个电场（磁场）分量环绕。这样不仅保证了介质分界面上切向场分量的连续性条件得到自然满足，而且
还允许旋度方程在空间上进行中心差分运算，同时也满足了法拉第电磁感应定律和安培环路积分定律，也可以很恰当地模拟电磁波的实际传播过程。 1.1Maxwell方程的差分形式 旋度方程为： 将其标量化，并将问题空间沿3个轴向分成若干网格单元，用Δx，Δy和Δz 分别表示每个网格单元沿3个轴向的长度，用Δt表示时间步长。网格单元顶点的坐标（x，y，z）可记为： 其中：i，j，k和n为整数。 同时利用二阶精度的中心有限差分式来表示函数对空间和时间的偏导数，即可得到如下FDTD基本差分式： 由于方程式里出现了半个网格和半个时间步，为了便于编程，将上面的差分式改写成如下形式：
课程论文
有限差分法在微分方程中的应用 本学期学习了《微分方程数值解》，本书中有限差分法给我留下的印象比较深刻，下边说说自己在方面的一点理解，请老师指正。 1.有限差分法的基本思想： 当系统的数学模型建立后，我们面对的主要问题就是微分积分方程的求解。基本思想是用离散的只含有限个未知量的差分方程组去近似地代替连续变量的微分方程和定解条件，并把差分方程组的解作为微分方程定解问题的近似解。将原方程及边界条件中的微分用差分来近似，对于方程中的积分用求和或及机械求积公式来近似代替，从而把原微分积分方程和边界条件转化成差分方程组。 2.有限差分法求解偏微分方程的步骤： 区域离散，即把所给偏微分方程的求解区域细分成由有限个格点组成的网格，这些离散点称作网格的节点； 近似替代，即采用有限差分公式替代每一个格点的导数。 逼近求解，换而言之，这一过程可以看作是用一个插值多项式及其微分来代替偏微分方程的解的过程。 从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值进行插值计算来近似得到。理论上，当网格步长趋近于零时，差分方程组的解应该收敛于精确解，但由于机器字节的限制，网格步长不可能也没有必要取得无限小，那么差分法的收敛性或者说算法的稳定性就显得至关重要。因此，在运用有限差分法时，除了要保证精度外，还必须要保证其收敛性。 3.构造差分法的几种形式： 主要草用的是泰勒级数展开的方法。其基本差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等。其中前两种形式为一阶计算精度，后一种为二阶计算精度。
利用有限差分法分析电磁场边界问题 在一个电磁系统中，电场和磁场的计算对于完成该系统的有效设计师极端重要的。例如，在系统中，用一种绝缘材料是导体相互隔离是，就要保证电场强度低于绝缘介质的击穿强度。在磁力开关中，所要求的磁场强弱，应能产生足够大的力来驱动开关。在发射系统中进行天线的有效设计时，关于天线周围介质中电磁场分布的知识显然有实质性的意义。 为了分析电磁场，我们可以从问题所涉及的数学公式入手。依据电磁系统的特性，拉普拉斯方程和泊松方程只能适合于描述静态和准静态（低频）运行条件下的情况。但是，在高频应用中，则必须在时域或频域中求解波动方程，以做到准确地预测电场和磁场，在任何情况下，满足边界条件的一个或多个偏微分方程的解，因此，计算电池系统内部和周围的电场和磁场都是必要的。 对电磁场理论而言，计算电磁场可以为其研究提供进行复杂的数值及解析运算的方法，手段和计算结果；而电磁场理论则为计算电磁场问题提供了电磁规律，数学方程，进而验证计算结果。常用的计算电磁场边值问题的方法主要有两大类，其每一类又包含若干种方法，第一类是解析法；第二类是数值法。对于那些具有最简单的边界条件和几何形状规则的（如矩形、圆形等）问题，可用分离变量法和镜像法求电磁场边值问题的解析解（精确解），但是在许多实际问题中往往由于边界条件过于复杂而无法求得解析解。在这种情况下，一般借助于数值法求解电磁场的数值解。 有限差分法，微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网络来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。 差分运算的基本概念： 有限差分法是指用差分来近似取代微分，从而将微分方程离散成为差分方程组。于是求解边值问题即转换成为求解矩阵方程[5]。 对单元函数 ()x f而言，取变量x的一个增量x?=h，则函数()x f的增量可以表示为 ()x f? = ()h x f+-()x f 称为函数()x f 的差分或一阶差分。函数增量还经常表示为 ()x f? = ? ? ? ? ? + 2 h x f - ? ? ? ? ? - 2 h x f
时域有限差分法发展综述 潘忠 摘要:时域有限差分法(FDTD)是解决复杂电磁问题的有效方法之一，目前FDTD 法的许多重要问题得到了很好的解决，已经发展成为一种成熟的数值计算方法。随着计算机数据处理性能的快速提高和计算机价格的下降，使得FDTD法的应用范围越来越广，而FDTD法本身在应用中又有新的发展.本文介绍并分析了时域有限差分法，对各种条件的应用进行了比较和分析，给出了具有一定参考价值的结论。 关键词:时域有限差分法;研究与发展;比较;分析 A Summary of FDTD and Development at Home and Abroad Zhong Pan Abstract: The finite difference time-domain (FDTD) method is one of the most effective methods to solve electromagnetic problems. Many important questions of FDTD method have been solved well through many scientists’ effort. Now, FDTD method is a mature numerical method. Especially in few years, the range of using FDTD method is becoming wider and wider because of the faster data processing and processing and cheaper price of computer. FDTD method has also been developed during using. FDTD method is introduced and discussed in this paper. The applications of various conditions are compared and analyzed. Finally, some valuable conclusions are drawn. Key words: FDTD; Research and Development; Comparison; Analysis 1966年，K.S.Yee首次提出电磁场数值计算的新方法—时域有限差分法(Finite Difference- Time Domain，简称FDTD)。经历了二十年的发展FDTD法才逐渐走向成熟。上世纪80年代后期以来FDTD法进入了一个新的发展阶段，即由成熟转为被广泛接受和应用的阶段。FDTD法是解决复杂问题的有效方法之一，是一种直接基于时域电磁场微分方程的数值算法，它直接在时域将Maxwell旋度方程用二阶精度的中心差分近似，从而将时域微分方程的求解转换为差分方程的迭代求解。是电磁场和电磁波运动规律和运动过程的计算机模拟。原则上可以求解任意形式的电磁场和电磁波的技术和工程问题，并且对计算机内存容量要求较低、计算速度较快、尤其适用于并行算法。现在FDTD法己被广泛应用于天线的分析与设计、目标电磁散射、电磁兼容、微波电路和光路时域分析、生物电磁剂量学、瞬态电磁场研究等多个领域。
班级：通信13-4 姓名： 学号： 指导教师：徐维 成绩： 电子与信息工程学院 信息与通信工程系
求解金属槽的电位分布 1.实验原理 利用有限差分法和matlab软件解决电位在金属槽中的分布。 有限差分法基本思想是把连续的定解区域用有限个离散点构成的网格来代替,这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似,积分用积分和来近似,于是原微分方程和定解条件就近似地代之以代数方程组,即有限差分方程组,解此方程组就可以得到原问题在离散点上的近似解.然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解.在采用数值计算方法求解偏微分方程时,若将每一处导数由有限差分近似公式替代,从而把求解偏微分方程的问题转换成求解代数方程的问题。 2.有限差分法 方程的定解问题就是在满足某些定解条件下求微分方程的解。在空间区域的边界上要满足的定解条件称为边值条件。如果问题与时间有关，在初始时刻所要满足的定解条件，称为初值条件。不含时间而只带边值条件的定解问题，称为边值问题。与时间有关而只带初值条件的定解问题，称为初值问题。同时带有两种定解条件的问题，称为初值边值混合问题。 定解问题往往不具有解析解，或者其解析解不易计算。所以要采用可行的数值解法。有限差分方法就是一种数值解法，它的基本思想是先把问题的定义域进行网格剖分，然后在网格点上，按适当的数值微分公式把定解问题中的微商换成差商，从而把原问题离散化为差分格式，进而求出数值解。此外，还要研究差分格式的解的存在性和唯一性、解的求法、解法的数值稳定性、差分格式的解与原定解问题的真解的误差估计、差分格式的解当网格大小趋于零时是否趋于真解（即收敛性），等等。 有限差分方法具有简单、灵活以及通用性强等特点，容易在计算机上实现。 2.1有限差分法原理
Efficiency enhancement of homoepitaxial InGaN/GaN light-emitting diodes on free-standing GaN substrate with double embedded SiO2 photonic crystals Tongbo Wei,* Ziqiang Huo, Yonghui Zhang, Haiyang Zheng, Yu Chen, Jiankun Yang, Qiang Hu, Ruifei Duan, Junxi Wang, Yiping Zeng, and Jinmin Li Semiconductor Lighting Technology Research and Development Center, Institute of Semiconductors, Chinese Academy of Sciences, Beijing 100083, China *tbwei@https://www.doczj.com/doc/8b18276572.html, Abstract: Homoepitaxially grown InGaN/GaN light emitting diodes (LEDs) with SiO2 nanodisks embedded in n-GaN and p-GaN as photonic crystal (PhC) structures by nanospherical-lens photolithography are presented and investigated. The introduction of SiO2 nanodisks doesn’t produce the new dislocations and doesn’t also result in the electrical deterioration of PhC LEDs. The light output power of homoepitaxial LEDs with embedded PhC and double PhC at 350 mA current is increased by 29.9% and 47.2%, respectively, compared to that without PhC. The corresponding light radiation patterns in PhC LEDs on GaN substrate show a narrow beam shape due to strong guided light extraction, with a view angle reduction of about 30°. The PhC LEDs are also analyzed in detail by finite-difference time-domain simulation (FDTD) to further reveal the emission characteristics. ?2014 Optical Society of America OCIS codes: (230.0230) Optical devices; (230.3670) Light-emitting diodes; (160.5298) Photonic crystals; (220.4241) Nanostructure fabrication. References and links 1. B. Monemar and B. E. Sernelius, “Defect related issues in the “current roll-off” in InGaN based light emitting diodes,” Appl. Phys. Lett. 91(18), 181103 (2007). 2. G. Verzellesi, D. Saguatti, M. Meneghini, F. Bertazzi, M. Goano, G. Meneghesso, and E. Zanoni, “Efficiency droop in InGaN/GaN blue light-emitting diodes: Physical mechanisms and remedies,” J. Appl. Phys. 114(7), 071101 (2013). 3. K. Akita, T. Kyono, Y. Yoshizumi, H. Kitabayashi, and K. Katayama, “Improvements of external quantum efficiency of InGaN-based blue light-emitting diodes at high current density using GaN substrates,” J. Appl. Phys. 101(3), 033104 (2007). 4. Y. Yang, X. A. Cao, and C. H. Yan, “Rapid efficiency roll-off in high-quality green light-emitting diodes on freestanding GaN substrates,” Appl. Phys. Lett. 94(4), 041117 (2009). 5. C.-L. Chao, R. Xuan, H.-H. Yen, C.-H. Chiu, Y.-H. Fang, Z.-Y. Li, B.-C. Chen, C.-C. Lin, C.-H. Chiu, Y.-D. Guo, J.-F. Chen, and S.-J. Cheng, “Reduction of Efficiency Droop in InGaN Light-Emitting Diode Grown on Self-Separated Freestanding GaN Substrates,” IEEE Photon. Technol. Lett. 23(12), 798–800 (2011). 6. M. J. Cich, R. I. Aldaz, A. Chakraborty, A. David, M. J. Grundmann, A. Tyagi, M. Zhang, F. M. Steranka, and M. R. Krames, “Bulk GaN based violet light-emitting diodes with high efficiency at very high current density,” Appl. Phys. Lett. 101(22), 223509 (2012). 7. X. A. Cao, S. F. LeBoeuf, M. P. D’Evelyn, S. D. Arthur, J. Kretchmer, C. H. Yan, and Z. H. Yang, “Blue and near-ultraviolet light-emitting diodes on free-standing GaN substrates,” Appl. Phys. Lett. 84(21), 4313 (2004). 8. Y. J. Zhao, J. Sonoda, C.-C. Pan, S. Brinkley, I. Koslow, K. Fujito, H. Ohta, S. P. DenBaars, and S. Nakamura, “30-mW-class high-power and high-efficiency blue (1011) semipolar InGaN/GaN light-emitting diodes obtained by backside roughening technique,” Appl. Phys. Express 3, 102101 (2010). 9. Y.-K. Fu, B.-C. Chen, Y.-H. Fang, R.-H. Jiang, Y.-H. Lu, R. Xuan, K.-F. Huang, C.-F. Lin, Y.-K. Su, J.-F. Chen, and C.-Y. Chang, “Study of InGaN-based light-emitting diodes on a roughened backside GaN substrate by a chemical wet-etching process,” IEEE Photon. Technol. Lett. 23(19), 1373–1375 (2011). #209568 - $15.00 USD Received 4 Apr 2014; revised 23 May 2014; accepted 26 May 2014; published 2 Jun 2014 (C) 2014 OSA30 June 2014 | Vol. 22, No. S4 | DOI:10.1364/OE.22.0A1093 | OPTICS EXPRESS A1093
伊犁师范学院硕士研究生 ————期末考核 科目：电磁波有限时域差分方法 姓名：姚伟 学号：1076411203009 学院：电子与信息工程学院 专业：无线电物理
时域有限差分法 1 选题背景 在多种可用的数值方法中，时域有限差分法(FDTD)是一种新近发展起来的可选方法。1966年，K.S.Yee 首次提出电磁场数值计算的新方法—时域有限差分法(Finite Difference- Time Domain ，简称FDTD)。经历了二十年的发展FDTD 法才逐渐走向成熟。上世纪80年代后期以来FDTD 法进入了一个新的发展阶段，即由成熟转为被广泛接受和应用的阶段。FDTD 法是解决复杂问题的有效方法之一，是一种直接基于时域电磁场微分方程的数值算法，它直接在时域将Maxwell 旋度方程用二阶精度的中心差分近似，从而将时域微分方程的求解转换为差分方程的迭代求解。是电磁场和电磁波运动规律和运动过程的计算机模拟。原则上可以求解任意形式的电磁场和电磁波的技术和工程问题，并且对计算机内存容量要求较低、计算速度较快、尤其适用于并行算法。现在FDTD 法己被广泛应用于天线的分析与设计、目标电磁散射、电磁兼容、微波电路和光路时域分析、生物电磁剂量学、瞬态电磁场研究等多个领域[1]。 2 原理分析 2.1 FDTD 的Yee 元胞 E,H 场分量取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理 t t ??=??=??E D H ε t t ??-=??-=??H B E μ 图1 Yee 模型 如图1所示，Yee 单元有以下特点[2]： 1）E 与H 分量在空间交叉放置，相互垂直；每一坐标平面上的E 分量四周由H 分量环绕，H 分量的四周由E 分量环绕；场分量均与坐标轴方向一致。 2）每一个Yee 元胞有8个节点，12条棱边，6个面。棱边上电场分量近似相等，用棱边的中心节点表示，平面上的磁场分量近似相等，用面的中心节点表示。 3）每一场分量自身相距一个空间步长，E 和H 相距半个空间步长 4）每一场分量自身相距一个时间步长，E 和H 相距半个时间步长，电场取n 时刻的值，磁场取n+0.5时刻的值；即：电场n 时刻的值由n-1时刻的值得到，磁场n+0.5时刻的值由n-0.5时刻的值得到；电场n 时刻的旋度对应n+0.5时刻的磁场值，磁场n+0.5时刻的
有限差分法 有限差分法（Finite Differential Method, FDM ） 什么是有限差分法 有限差分法是指用泰勒技术展开式将变量的导数写成变量，在不同时间或空间点值的差分形式的方法。 按时间步长和空间步长将时间和空间区域剖分成若干网格，用未知函数在网格结(节)点上的值所构成的差分近似代替所用偏微分方程中出现的各阶导数，从而把表示变量连续变化关系的偏微分方程离散为有限个代数方程，然后解此线性代数方程组，以求出溶质在各网格结(节)点上不同时刻的浓度。 有限差分法的基本步骤 （1）剖分渗流区，确定离散点。将所研究的水动力弥散区域按某种几何形状(如矩形、任意多边形等)剖分成网络系统。 （2）建立水动力弥散问题的差分方程组。 （3）求解差分方程组。采用各种迭代法，如点逐次超松驰方法(SOR)、线逐次超松驰方法(LSOR)、迭代的交替方向隐式方法(IADI)及强隐式方法(SID)等。 (1) 现在分别对时间（从0时刻到到期日）和股票价格（S max )为可达到的足够高的股票价格） 进行分割，即\triangle S=S_{max}/M,\triangle T/N,这样就分别有N+1个时间段和M+1个股票价格，建立如图(所示的坐标方格，将定解区域网格化，坐标方格上的点（i,j ）对应时刻 和股票价格 ，用变量f i ,j 表示（i,j ）点的期权价格。 2.建立差分格式 （1）内含的有限差分方法 其步骤可分为以下几步：
(1)求前向差分近似:(2) 后向差分格式：(3) 将(2),(3)式平均可更加对称地求出的近似，即 (4) (2)求用前向差分近似： (5) (3)求 (6) (4)将(4),(5),(6)式代入(1)式可得到内含有限差分公式： + b j f i,j?c j f i,j + 1 = f i + 1,j(7) a j f i,j? 1 其中： i=0,1,…，N-1。j=0,1…，M-1 针对看跌期权和看涨期权可分别求出方程的边界条件： 看跌期权： 看涨期权： (5)利用边界条件和(7)式可以给出M-1个联立方程组： + b j f N? 1,j + c j f N? 1,j + 1j=1,2…，M-1 a j f N? 1,j? 1
有限元法与有限差分法的主要区别 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分的空间形式来考虑，可分为中心格式和逆风格式。考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。构造差分的方法有多种形式，目前主要采用的是泰勒级数展开方法。其基本的差分表达式主要有三种形式：一阶向前差分、一阶向后差分、一阶中心差分和二阶中心差分等，其中前两种格式为一阶计算精度，后两种格式为二阶计算精度。通过对时间和空间这几种不同差分格式的组合，可以组合成不同的差分计算格式。有限元方法的基础是变分原理和加权余量法，其基本求解思想是把计算域划分为有限个互不重叠的单元，在每个单元内，选择一些合适的节点作为求解函数的插值点，将微分方程中的变量改写成由各变量或其导数的节点值与所选用的插值函数组成的线性表达式，借助于变分原理或加权余量法，将微分方程离散求解。采用不同的权函数和插值函数形式，便构成不同的有限元方法。有限元方法最早应用于结构力学，后来随着计算机的发展慢慢用于流体力学的数值模拟。在有限元方法中，把计算域离散剖分为有限个互不重叠且相互连接的单元，在每个单元内选择基函数，用单元基函数的线形组合来逼近单元中的真解，整个计算域上总体的基函数可以看为由每个单元基函数组成的，则整个计算域内的解可以看作是由所有单元上的近似解构成。在河道数值模拟中，常见的有限元计算方法是由变分法和加权余量法发展而来的里兹法和伽辽金法、最小二乘法等。根据所采用的权函数和插值函数的不同，有限元方法也分为多种计算格式。从权函数的选择来说，有配置法、矩量法、最小二乘法和伽辽金法，从计算单元网格的形状来划分，有三角形网格、四边形网格和多边形网格，从插值函数的精度来划分，又分为线性插值函数和高次插值函数等。不同的组合同样构成不同的有限元计算格式。对于权函数，伽辽金(Galerkin)法是将权函数取为逼近函数中的基函数；最小二乘法是令权函数等于余量本身，而内积的极小值则为对代求系数的平方误差最小；在配置法中，先在计算域内选取N个配置点。令近似解在选定的N个配置点上严格满足微分方程，即在配置点上令方程余量为0。插值函数一般由不同次幂的多项式组成，但也有采用三角函数或指数函数组成的乘积表示，但最常用的多项式插值函数。有限元插值函数分为两大类，一类只要求插值多项式本身在插值点取已知值，称为拉格朗日(Lagrange)多项式插值；另一种不仅要求插值多项式本身，还要求它的导数值在插值点取已知值，称为哈密特(Hermite)多项式插值。单元坐标有笛卡尔直角坐标系和无因次自然坐标，有对称和不对称等。常采用的无因次坐标是一种局部坐标系，它的定义取决于单元的几何形状，一维看作长度比，二维看作面积比，三维看作体积比。在二维有限元中，三角形单元应用的最早，近来四边形等参元的应用也越来越广。对于二维三角形和四边形电源单元，常采用的插值函数为有La grange插值直角坐标系中的线性插值函数及二阶或更高阶插值函数、面积坐标系中的线性插值函数、二阶或更高阶插值函数等。对于有限元方法，其基本思路和解题步骤可归纳为(1)建立积分方程，根据变分原理或方程余量与权函数正交化原理，建立与微分方程初边值问题等价的积分表达式，这是有限元法的出发点。(2)区域单元剖分，根据求解区域的形状及实际问题的物理特点，将区域剖分为若干相互连接、不重叠的单元。区域单元划分是采用有限元方法的前期准备工作，这部分工作量比较大，除了给计算单元和节点进行编号和确定相互之间的关系之外，还要表示节点的位置坐标，同时还需要列出自然边界和本质边界的节点序号和相应的边界值。(3)确定单元基函数，根据单元中节点数目及对近似解精度的要求，选择满足一定插
% 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 % hagness@https://www.doczj.com/doc/8b18276572.html, % % 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. %
有限差分法 有限差分法有限差分法 finite difference method 微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散 点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函 数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似， 积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差 分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便 可以从离散解得到定解问题在整个区域上的近似解。 有限差分法的主要内容包括：如何根据问题的特点将定解区域作网格剖分；如何把原 微分方程离散化为差分方程组以及如何解此代数方程组。此外为了保证计算过程的可行和 计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分 格式的相容性、收敛性和稳定性。对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。另外，一个差分格 式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。此外，还有一个重要的概念必须考虑，即差分格式的稳定性。因为差分格式的计算过 程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致 差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以 控制的，就认为格式是稳定的。只有在这种情形，差分格式在实际计算中的近似解才可能 任意逼近差分方程的精确解。关于差分格式的构造一般有以下3种方法。最常用的方法是 数值微分法，比如用差商代替微商等。另一方法叫积分插值法，因为在实际问题中得出的 微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。此外还可以用 待定系数法构造一些精度较高的差分格式。 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛运用。该方法 将求解域划分为差分网格，用有限个网格节点代替连续的求解域。有限差分法以Taylor 级数展开等方法，把控制方程中的导数用网格节点上的函数值的差商代替进行离散，从 而建立以网格节点上的值为未知数的代数方程组。该方法是一种直接将微分问题变为代数 问题的近似数值解法，数学概念直观，表达简单，是发展较早且比较成熟的数值方法。 对于有限差分格式，从格式的精度来划分，有一阶格式、二阶格式和高阶格式。从差分 的空间形式来考虑，可分为中心格式和逆风格式。 考虑时间因子的影响，差分格式还可以分为显格式、隐格式、显隐交替格式等。目 前常见的差分格式，主要是上述几种形式的组合，不同的组合构成不同的差分格式。差分 方法主要适用于有结构网格，网格的步长一般根据实际地形的情况和柯朗稳定条件来决定。
时域有限差分法 1 选题背景 在多种可用的数值方法中，时域有限差分法(FDTD)是一种新近发展起来的可选方法。1966年，K.S.Yee 首次提出电磁场数值计算的新方法—时域有限差分法(Finite Difference- Time Domain ，简称FDTD)。经历了二十年的发展FDTD 法才逐渐走向成熟。上世纪80年代后期以来FDTD 法进入了一个新的发展阶段，即由成熟转为被广泛接受和应用的阶段。FDTD 法是解决复杂问题的有效方法之一，是一种直接基于时域电磁场微分方程的数值算法，它直接在时域将Maxwell 旋度方程用二阶精度的中心差分近似，从而将时域微分方程的求解转换为差分方程的迭代求解。是电磁场和电磁波运动规律和运动过程的计算机模拟。原则上可以求解任意形式的电磁场和电磁波的技术和工程问题，并且对计算机内存容量要求较低、计算速度较快、尤其适用于并行算法。现在FDTD 法己被广泛应用于天线的分析与设计、目标电磁散射、电磁兼容、微波电路和光路时域分析、生物电磁剂量学、瞬态电磁场研究等多个领域[1]。 2 原理分析 2.1 FDTD 的Yee 元胞 E,H 场分量取样节点在空间和时间上采取交替排布，利用电生磁，磁生电的原理 t t ??=??=??E D H ε t t ??-=??-=??H B E μ 图1 Yee 模型 如图1所示，Yee 单元有以下特点[2]： 1）E 与H 分量在空间交叉放置，相互垂直；每一坐标平面上的E 分量四周由H 分量环绕，H 分量的四周由E 分量环绕；场分量均与坐标轴方向一致。 2）每一个Yee 元胞有8个节点，12条棱边，6个面。棱边上电场分量近似相等，用棱边的中心节点表示，平面上的磁场分量近似相等，用面的中心节点表示。 3）每一场分量自身相距一个空间步长，E 和H 相距半个空间步长 4）每一场分量自身相距一个时间步长，E 和H 相距半个时间步长，电场取n 时刻的值，磁场取n+0.5时刻的值；即：电场n 时刻的值由n-1时刻的值得到，磁场n+0.5时刻的值由n-0.5时刻的值得到；电场n 时刻的旋度对应n+0.5时刻的磁场值，磁场n+0.5时刻的旋度对应(n+0.5)+0.5时刻的电场值，逐步外推。 5）3个空间方向上的时间步长相等，
有限差分法有限差分法 finite difference method 微分方程和积分微分方程数值解的方法。基本思想是把连续的定解区域用有限个离散点构成的网格来代替，这些离散点称作网格的节点；把连续定解区域上的连续变量的函数用在网格上定义的离散变量函数来近似；把原方程和定解条件中的微商用差商来近似，积分用积分和来近似，于是原微分方程和定解条件就近似地代之以代数方程组，即有限差分方程组，解此方程组就可以得到原问题在离散点上的近似解。然后再利用插值方法便可以从离散解得到定解问题在整个区域上的近似解。 有限差分法的主要内容包括：如何根据问题的特点将定解区域作网格剖分；如何把原微分方程离散化为差分方程组以及如何解此代数方程组。此外为了保证计算过程的可行和计算结果的正确，还需从理论上分析差分方程组的性态，包括解的唯一性、存在性和差分格式的相容性、收敛性和稳定性。对于一个微分方程建立的各种差分格式，为了有实用意义，一个基本要求是它们能够任意逼近微分方程，这就是相容性要求。另外，一个差分格式是否有用，最终要看差分方程的精确解能否任意逼近微分方程的解，这就是收敛性的概念。此外，还有一个重要的概念必须考虑，即差分格式的稳定性。因为差分格式的计算过程是逐层推进的，在计算第n＋1层的近似值时要用到第n层的近似值，直到与初始值有关。前面各层若有舍入误差，必然影响到后面各层的值，如果误差的影响越来越大，以致差分格式的精确解的面貌完全被掩盖，这种格式是不稳定的，相反如果误差的传播是可以控制的，就认为格式是稳定的。只有在这种情形，差分格式在实际计算中的近似解才可能任意逼近差分方程的精确解。关于差分格式的构造一般有以下3种方法。最常用的方法是数值微分法，比如用差商代替微商等。另一方法叫积分插值法，因为在实际问题中得出的微分方程常常反映物理上的某种守恒原理，一般可以通过积分形式来表示。此外还可以用待定系数法构造一些精度较高的差分格式。 有限差分方法(FDM)是计算机数值模拟最早采用的方法，至今仍被广泛