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