第三章 有限差分法函数()f x ，x 为定义在区间[]a b ,上的连续变 量。
将区间[]a b ,等分成n 份，令()h b a n =-称为 步长，x 在这些离散点处的取值为x a ih i =+ ()i n =01,,,称为节点。
函数()f x 在这些节点处的差值()()()()()()f x h f x f x f x h f x h f x h i i i i i i +---+--⎧⎨⎪⎩⎪(5-1)分别称为一阶向前、向后和中心差分，可以用它 们作为函数()f x 在x i 处的微分近似值。
这些差分 与相应x 区间的比值()()[]()()[]()()[]1112h f x h f x h f x f x h h f x h f x h i i i i i i +---+--⎧⎨⎪⎪⎪⎩⎪⎪⎪ (5-2) 分别称为一阶向前、向后和中心差商，可以用它 们作为函数()f x 在x i 处的导数近似值。
完全类似地可以定义高阶差商，例如常用的二阶中心差商()()()[]122hf x h f x f x h i i i +-+- (5-3)可以作为函数()f x 在x i 处的二阶导数近似值。
§3.1 常微分方程初值问题的差分解法考虑电学中的一个问题：如图5-1。
研究 电容器上的电荷随时间的变化规律。
图5-1 RC 放电回路这个问题对应的微分方程及其定解条件为：d d Q tQ RC QQ t =-=⎧⎨⎪⎩⎪=00(5-4) 这是一阶微分方程的初值问题，它的解析解为 Q Q e t RC =-0 (5-5)一、欧拉(Euler )折线法求解下列普遍形式的一阶微分方程的初值 问题：()[]()'=∈=⎧⎨⎪⎩⎪y f x y x a b y a y ,,0(5-6) 首先，将区间[]a b ,等分n 份，取值a x x xb n =<<<=01 ，步长h x x i i =-+1。
然后，用一阶向前差商近似一阶导数，即()()()()[]y x y x hy x f x y x i i i i i +-≈'=1, (5-7) 简记()y x y i i ≈，则式(5-7)可以写成差分格式：()y y h f x y i i i i +=+⋅1, ()i n =-011,,,(5-8)此即向前欧拉差分格式。
这是一个递推计算格式， 从区间左端点即式(5-6)中的初始条件出发，按式 (5-8)依次可以算到区间右端点，得到的 y y y n 12,,, 就是原方程解()y x 的近似值。
应用式(5-8)计算RC 放电方程(5-4)，按SI 单 位制，取Q 010=，RC =8，时间步长h =1，计 算结果如下：从表5-1可见，时间较小时，计算值与解析 值比较一致，而随着时间的增加，两者之差有增 加的趋势，其原因可以从向前欧拉差分格式(5-8) 的几何意义得到解释。
图5-2 向前欧拉差分格式的几何意义如图5-2，一开始从y 0求y 1时，公式为 ()y y h f x y 1000=+⋅,也可以将其写为()()y y f x y x x 100010-=-,这是过点()P x y 000,、斜率为()f x y 00,的直线方 程，由y 0求y 1的过程相当于从点()P x y 000,沿该 直线求点()P x y 111,的过程，直线上的y 1与曲线上 的()y x 1一般是不相等的，这就是由于采用直线近 似曲线产生的计算误差，并且由此往后的每一步 计算都是如此，欧拉法实际上是用一条折线近似 原来的曲线，所以又称为欧拉折线法。
关于欧拉法的局部截断误差，可以利用级数 展开的方法进行讨论，在x i 点处作泰勒级数展开：()()()()y x y x h y x h y i i i +=+'+''122ξ []ξ∈a b , (5-9)与欧拉差分格式(5-8)比较可知，局部截断误差为()()()22112h O y h y x y R i i =''=-=++ξ (5-10) 记号()O h 2表示局部截断误差是h 的二次方数量级，此时称该计算具有一阶精度。
一般情况下， 如果局部截断误差为()O h m +1，则称有m 阶精度。
如果采用一阶向后差商近似一阶导数，即()()()()[]y x y x hy x f x y x i i i i i -≈'=-1, (5-11) 并将所有角标加1，可以得到向后欧拉差分格式：()y y h f x y i i i i +++=+⋅111, ()i n =-011,,, (5-12)计算时，要从式(5-12)中解出y i +1。
如果称向前 欧拉差分格式(5-8)为显式格式的话，那么向后 欧拉差分格式(5-12)就是隐式格式。
向后欧拉差分格式(5-12)的局部截断误差也 是()O h 2。
我们还可以采用一阶中心差商近似一阶导 数，即()()()()[]y x y x hy x f x y x i i i i i +--≈'=112, (5-13) 可以得到中心欧拉差分格式：()y y h f x y i i i i +-=+⋅112, ()i n =-011,,,(5-14)如果我们将式(5-8)和式(5-12)称为单步格式的话， 那么式(5-14)就是两步格式。
中心欧拉差分格式的局部截断误差为()O h 3，比前两者精度高一阶。
二、Lax 等价定理关于微分方程的适定性，差分方程的相容性、 收敛性和稳定性。
如图5-3所示。
图5-3 微分方程的适定性与差分方程的相容性、收敛性和稳定性之间的关联微分方程的适定性：定解问题解的存在性、唯一性和稳定性统称 为定解问题的适定性。
满足存在性、唯一性和稳 定性的定解问题称为适定的，否则就是不适定的。
相容性、收敛性和稳定性是差分方程的三个 最基本的性质。
相容性是指差分方程是否逼近原微分方程的问题，它不涉及差分方程的解。
如果将微分方程 的解代入到差分方程之中所得到的局部截断误差 随计算步长趋于零时而趋于零，则称差分方程和 原微分方程是相容的。
收敛性是指差分方程的解收敛于原微分方程 的解。
设原微分方程的准确解为U ，差分方程的 准确解为u ，如果当计算步长趋于零时，U u -→0，则称差分方程收敛于原微分方程。
其中⋅是任意一种范数。
稳定性是指差分方程对于初始误差的抗扰性， 它和微分方程无关，是差分方程本身的一个特性。
稳定性概念较为严格的叙述为：如果对于初始误 差ε0，存在与计算步长无关的常数L ，使得任意 一步上的误差εi ，在一定范数下满足不等式： εεi L ≤0，则称差分格式是稳定的。
Lax 等价定理：定理5-1：如果线性微分方程的初值问题是 适定的，差分格式对它的逼近是相容的，则差分 格式的收敛性和稳定性是等价的。
在具体计算工作中，相容性是容易满足的， 只要保证了差分格式的稳定性，利用Lax 等价 定理，我们也就保证了差分格式的收敛性，这 是Lax 等价定理的重要意义所在。
三、龙格--库塔(Runge--Kutta )法欧拉法是在小区间的左端点求斜率建立直线 方程以近似曲线；为了减少由此造成的误差，可 以用小区间左、右端点斜率的平均值建立直线方 程以近似曲线，这就是梯形法。
该方法的斜率为：()()[]()()[]K y x y x f x y f x y i i i i i i ='+'=++++1212111,,(5-15)计算公式为：()()[]y y hf x y f x y i i i i i i +++=++1112,, (5-16)这是一个单步、隐式格式，截断误差为()O h 3。
当y i +1不方便以显式求出时，通常可以采用迭代 法求解，即由第i 步计算第 i +1步，采用如下 迭代格式：()()()()[]y y h f x y f x y i k i i i i i k ++++=++11112,, (5-17)为了得到显式格式，也可以对式(5-16)右边 的y i +1进行一级近似，这就是二阶龙格--库塔法：()y y hK K i i +=++1122(5-18)其中：()()K f x y K f x y hK i i i i 1211==+⎧⎨⎪⎩⎪+,, (5-19) 为了进一步提高计算精度，还可以对小区间 左、右端点以及中点的斜率近似值进行加权平均， 得到应用广泛的四阶龙格--库塔法：()y y hK K K K i i +=++++11234622 (5-20)其中：()()⎪⎪⎪⎩⎪⎪⎪⎨⎧++=⎪⎭⎫ ⎝⎛++=⎪⎭⎫ ⎝⎛++==3423121,2,22,2,K h y h x f K K h y h x f K K h y h x f K y x f K i i ii ii i i (5-21)四阶龙格--库塔法的截断误差为()O h 5。
对于一般计算要求，四阶龙格--库塔格式精度 已足够高，而且计算稳定、易于编程。
四、一阶微分方程组初值问题的数值解法为简单，讨论只有两个方程组成的一阶微分 方程组的初值问题：()()()()'='===⎧⎨⎪⎪⎩⎪⎪y f x y z z g x y z y x y z x z,,,,0000 (5-22) 为了求解式(5-22)，只需将y 和z 看作一个 矢量的两个分量，则前面所述的各种差分解法 就可直接推广到方程组中加以应用。
以四阶龙 格--库塔法为例，方程组(5-22)的计算格式为：()()y y h K K K K z z h L L L L i i i i ++=++++=++++⎧⎨⎪⎪⎩⎪⎪1123411234622622 (5-23)其中：()()K f x y z L g x y z K f x h y h K z h L L g x h y h K z h L K f x h y h K z h L L g x h y h K z h L K f x i i i i i i i i ii i ii i i i i i 112112113223224222222222222===+++⎛⎝ ⎫⎭⎪=+++⎛⎝ ⎫⎭⎪=+++⎛⎝⎫⎭⎪=+++⎛⎝⎫⎭⎪=,,,,,,,,,,,,()()i i i i i i h y hK z hL L g x h y hK z hL+++=+++⎧⎨⎪⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎪,,,,33433 (5-24)下面是几个可以用此种方法求解的物理问题。
[例题5-1] 洛伦兹(Lorenz )吸引子()d d d d d d xt x y yt xz r x y ztx y bz =--=-+-=-⎧⎨⎪⎪⎪⎩⎪⎪⎪σ (5-25)其中，σ、r 、b 是来自流体力学中的参量。