有限差分法

有限差分法以变量离散取值后对应的函数值来近似微分方程中独立变量的连续取值。
在有限差分方法中，我们放弃了微分方程中独立变量可以取连续值的特征，而关注独立变量离 散取值后对应的函数值。但是从原则上说，这种方法仍然可以达到任意满意的计算精度。因为方程的 连续数值解可以通过减小独立变量离散取值的间格，或者通过离散点上的函数值插值计算来近似得 到。这种方法是随着计算机的诞生和应用而发展起来的。其计算格式和程序的设计都比较直观和简单， 因而，它的实际应用已经构成了计算数学和计算物理的重要组成部分。
f ( xi + h) − f ( xi − h).
节点 xi 的一阶向前差分 节点 xi 的一阶向后差分 节点 xi 的一阶中心差分
2
与 xi 点相邻两点的泰勒展开式可以写为
f (xi
− h) =
f
( xi
)
−
hf
'
( xi
)
+
h2 2
f
''
( xi
)
−
h3 3!
f
'''
( xi
)
+
h4 4!
f
−
φ
0
)
⎤ ⎥
⎦
+
f 0φ0
=
q0
（4.2.18）
如果在 x 和 y 方向的步长分别相等， 即 h1 = h3 = hx 和 h2 = h4 = hy 时，则上式化为
φ1
− 2φ 0 hx 2
+φ3
+ φ2
− 2φ 0 hy 2
+φ4
+
f 0φ 0
=
q0
，
（4.2.19）
一般可以用角标来表示节点的标记，将上式写为
⎞ ⎟
的差分表达式：
⎝ ∂y 2 ⎠ 0
⎛ ⎜ ⎝
∂
2
φ
⎞ ⎟
∂y 2 ⎠ 0
≈
2
h4 (φ 2
− φ 0 ) + h2 (φ 4 h2h4 (h2 + h4 )
−φ0)
.
当采用等步长 h2 = h4 = hy 时， 有
⎛ ⎜
∂
2φ
⎞ ⎟
≈ φ 2 − 2φ 0 + φ 4 .
⎝ ∂y 2 ⎠ 0
hy 2
⎛ ∂φ ⎞ ⎝⎜ ∂x ⎠⎟ 0
≈
φ1 −φ0 h1
.
或
（4.2.7）
⎛ ∂φ ⎞ ⎝⎜ ∂x ⎠⎟ 0
≈
φ0 −φ3 h3
.
显然这种单侧差商的误差较大。
（4.2.8）
如果要寻求更精确的差分格式，我们可以引入待定常数 α , β ，由 φ1 和 φ3 的泰勒展开，构造出如下的 关系式：
α (φ1
− φ0
'''' (xi )−...... .
f (xi + h) =
f
( xi
)
+
hf
' ( xi
)
+
h2 2
f '' ( xi ) +
h3 3!
f
''' ( xi ) +
h4 4!
f
'''' (xi )+...... .
(4.1.1) (4.1.2)
（4.1.1）-(4.1.2)，并忽略 h 的平方和更高阶的项得到一阶微分的中心差商表示：
f1 + f2
+
f3 + h2
f4
− 4 f0
−
2h 2 4!
∂4f ( ∂x4
+
∂4f ∂y 4
).
(4.1.7)
j+1 (i-1, j+1)
(i, j+1) 2
(i+1, j+1)
h2
j
(i-1, j)
3
0
(i, j)
1 (i+1, j)
h4
h3
4
(i+1, j-1)
j-1
(i-1, j-1) (i, j-1)
第四章 有限差分方法
4.1 引言
有限差分法：数值求解常微分方程或偏微分方程的方法。
物理学和其他学科领域的许多问题在被分析研究之后, 往往可以归结为常微分方程或偏微分方 程的求解问题。一般说来，处理一个特定的物理问题，除了需要知道它满足的数学方程外，还应当同 时知道这个问题的定解条件，然后才能设计出行之有效的计算方法来求解。
..., x = x + h = b , 然
n+1
n
后求出 f(x)在这些点上的近似值。显然步长 h 越小，近似解的精度就越好。与节点 xi 相邻的节点有
xi − h 和 xi + h ，因此在 xi 点可以构造如下形式的差值：
f (xi + h) − f (xi ), f (xi ) − f (xi − h),
⎞ ⎟ ⎠
0
≈
2
h3 (φ 1
ຫໍສະໝຸດ Baidu
−φ0) + h1h3 (h1
h1 (φ 3 + h3 )
−
φ0)
.
当用等步长 h1 = h3 = hx 时，上式成为
（4.2.14）
10
⎛ ⎜
∂
2φ
⎞ ⎟
≈ φ1 − 2φ 0 + φ 3 .
⎝ ∂x2 ⎠ 0
hx 2
用完全相同的计算方法，我们可以推导出
⎛ ⎜
∂
2φ
我们现在考虑方程(4.2.1)中 p 为常数的二维的情况，我们可以得到下面的方程：
∂ 2φ + ∂ 2φ + f (x, y)φ = q(x, y). ∂x 2 ∂y 2
设函数φ 在区域 D 内满足方程（4.2.6）式(区域 D 的边界为 G)。
(4.2.6)
7
区域 D 的离散化：
即通过任意的网络划分方法把区域 D 离散为许许多多的小单元。原则上讲这种网格分割是可以任 意的，但是在实际应用中，常常是根据边界 G 的形状，采用最简单，最有规律，和边界的拟合程度最 佳的方法来分割。常用的有正方形分割法和矩形分割法（如图 4.2.1）。有时也用三角形分割法（见 图 4.2.2）。对圆形区域，应用图（4.2.3）所示的极网络格式也许更方便些。这些网络单元通常称 为元素，网络点称为节点。
f ' (xi ) ≈
f (xi + h) − f (xi 2h
− h)
.
(4.1.3)
利用(4.1.1)和(4.1.2)式我们还可以得到一阶微分的向前，向后一阶差商表示
f ' (xi ) ≈
f (xi + h) − h
f (xi ) ,
f ' (xi ) ≈
f (xi ) − f (xi − h) h
假定某方程形式上可以写为：
Lφ = q.
其中 L 为含偏微商的算符．
(4.2.1)
它的边界条件一般可写为：
φ
|G
+ g1(s)
∂φ ∂n
|G =
g2 (s).
G 表示场域 D 的边界， g1(s), g2 (s) 为边界上 s 点的逐点函数。
(4.2.2)
6
三类边界条件：
（1）第一类边界条件，或称为狄利克莱(Dirichlet)问题 (g1 = 0, g2 ≠ 0) 。
h1
j-1
j
i+1
图 4.1.1 节点 0 及邻近节点
4
对微分方程数值求解的误差的来源：
（1）方法误差（或截断误差）。这是由于采用的计算方法所引起的误差。例如上面我们介绍的差商
表示中，采用的泰勒展开式展开到第 n + 1项时的截断误差阶数为O(hn+1)。具体方法的误差阶数取决于
在离散化时的近似阶数。因此若改进算法就可以减小截断误差。 （2）舍入误差（或计算误差）。这是由于计算机的有限字长而造成数据在计算机中的表示出现误差。 在计算机运算的过程中，随着运算次数的增加舍入误差会积累得很大。如果在多次运算后，舍入误差 的精度影响是有限的，那么这个算法是稳定的，否则是不稳定的。不稳定的算法是不能用的。
)
+
β (φ3
− φ0 )
=
⎜⎛ ⎝
∂φ ∂x
⎟⎞ ⎠0
(αh1
−
βh3 )
+
1 2
⎜⎜⎝⎛
∂ 2φ ∂x 2
⎟⎟⎠⎞0
(αh12
+
βh32 )
+
....
令⎛ ⎜
∂
2φ
⎞ ⎟
项的系数为零，则得到α 和 β
之间应当满足
⎝ ∂x 2 ⎠ 0
α
=
−
h32 h12
β.
（4.2.9） （4.2.10）
9
将公式（4.2.10）带入（4.2.9），并舍去高阶项，得到
本书中我们将略去对差分法稳定性和收敛性理论的讨论，尽管这方面的内容是相当重要的。以下 的讨论中所讲到的各种差分格式，我们均假定求解方法满足稳定性和收敛性的要求。
5
4.2 有限差分法和偏微分方程
利用上节所介绍的微分的差分表示，我们就很容易地将微分方程离散化为差分方程组的形式。但 是由差分方程所得的解完全取决于待求微分方程的特性。正如我们在物理上所知道的，边界条件的情 况变化将会引起差分方程组的不同。在求解微分方程中，我们会遇到两类问题：一类是初始值问题； 另一类是边值条件的问题。在初始值问题中，部分边界上的函数值和部分的函数偏导值是给定的。通 常在这类问题中的独立变量之一是时间 t。在边界值问题中，边界上的信息是给定的。本书中我们仅 讨论后一类问题。
y
y
D
D
0
x0
x
4.2.1 求解区域的矩形分割。 4.2.2 求解区域的三角形分割。 4.2.3 求解区域的极网络分割。
8
用节点上的函数值来表示节点上偏导的数值。
设区域内部某节点 0 附近的各节点如图 4.1.1 所示。这里我们取步长 h 不相等的最一般情况。以
φ0,φ1,φ2,φ3,φ4 分别代表在节点 0，1，2，3，4 处φ 的函数值。如前所述，0 点的一阶偏导数可以通过先 前或向后的差商，由 1 和 3 节点近似写出
（4.2.15） （4.2.16） （4.2.17）
将公式（4.2.14）和（4.2.16）两式代入方程（4.2.13），我们就得到该方程的差分表达式为
(∇ 2φ )0
=
⎡ 2⎢
h3
(φ1
⎣
−φ0) + h1h3 (h1
h1 (φ3 + h3 )
− φ0 )
+
h4 (φ2 − φ0 ) + h2 (φ4 h2 h4 (h2 + h4 )
1
有限差分法的具体操作分为两个部分：
(1)用差分代替微分方程中的微分，将连续变化的变量离散化，从而得到差分方程组的数学形式； (2)求解差分方程组。 在第一步中，我们通过所谓的网络分割法，将函数定义域分成大量相邻而不重合的子区域。通 常采用的是规则的分割方式。这样可以便于计算机自动实现和减少计算的复杂性。网络线划分的交点 称为节点。若与某个节点 P 相邻的节点都是定义在场域内的节点，则 P 点称为正则节点；反之，若节 点 P 有处在定义域外的相邻节点，则 P 点称为非正则节点。 在第二步中，数值求解的关键就是要应 用适当的计算方法，求得特定问题在所有这些节点上的离散近似值。
上述差分步骤应用于偏微分：
例如，对于 f
=
f
(x,
y)
的情况，拉普拉斯算符在
0
点作用在此函数上的值
⎛ ⎜
∇
2
⎝
f
=
⎛ ⎜ ⎝
∂ ∂
2f x2
+∂2f ∂ y2
⎞ ⎟
⎞ ⎟
，也可
⎠⎠
以用临近的点上的函数值来表示出来。(见图 4.1.1, 且 h1 = h2 = h3 = h4 = h 时)
∇2 f ≈
⎛⎝⎜
∂φ ∂x
⎞⎠⎟
0
的另一个差分表达式
⎛⎝⎜
∂φ ∂x
⎞⎠⎟
0
≈
h32 (φ1
−φ0) − h1h3 (h1
h12 (φ 3 + h3 )
−φ0)
.
当选用等步距 h1 = h3 = hx 时，上式成为
（4.2.11）
⎛ ∂φ ⎞ ⎝⎜ ∂x ⎠⎟ 0
≈
φ1 −φ3 2hx
.
（4.2.12）
（此为中心差商表达式。）
(4.2.4)
L ≡ −∇( p∇) + f .
(4.2.5)
公式中的 p 和 f 是给定的函数。我们将会得到一类很重要的微分方程。它是在流体力学、等离子物 理、天体物理…等学科中，势函数起关键作用的许多问题当中的基本方程。当 p=1, f=0 时，我们得 到（4.2.1）式的特殊情况—泊松(Poisson)方程。
二阶偏导数的差分表达式（“五点格式”或“菱形格式”）
在（4.2.9）式中，如果令
⎛ ⎝⎜
∂φ ∂x
⎞ ⎠⎟
0
的系数为零，则有
α
和
β
间存在关系式：
α = h3 β. h1
（4.2.13）
将上式（4.2.13）代入（4.2.9）式中，并忽略 h 三阶以上的高次项，则得到表达式：
⎛ ⎜ ⎝
∂ 2φ ∂x 2
.
(4.1.4) (4.1.5)
将(4.1.1)和(4.1.2)式相加，忽略 h 的立方及更高阶的项得到二阶微分的中心差商表示
f '' ( xi ) ≈
f (xi + h) − 2 f (xi ) + h2
f (xi − h)
.
(4.1.6)
3
利用(4.1.3) ~ (4.1.6)式，我们就可以构造出微分方程的差分格式。这里要指出的是：在构造 差分格式时，究竟应该选择向前，向后还是中间差分或差商来代替微分方程中的微分或微商，应当根 据由此得到的差分方程解的稳定性和收敛性来考虑。同时兼顾到差分格式的简单和求解的方便。
有限差分法的差分格式：
一个函数在 x 点上的一阶和二阶微商，可以近似地用它所临近的两点上的函数值的差分来表示。
如对一个单变量函数 f(x)，x 为定义在区间[a,b]的连续变量。以步长 h= Δ x 将[a,b]区间离散化，我
们得到一系列节点 x = a , 1
x=x +h,
2
1
x =x
3
2
+ h = a + 2Δx ,
φ |G = g(s).
(4.2.3)
（2）第二类边界条件，或称诺伊曼(Neumann)问题 (g1 ≠ 0, g2 = 0) 。
∂φ ∂n
|G
=
g ( s ).
（3）第三类边界条件，或称混合问题 (g1 ≠ 0, g2 ≠ 0) 。
对于算符 L 为斯杜-刘维尔(Sturm-Liouville)算符的特定情况，即