第8章 常微分方程边值问题的数值解法
8.1 引 言
第7章介绍了求解常微分方程初值问题的常用的数值方法；本章将介绍常微分方程的边值问题的数值方法。
只含边界条件(boundary-value condition)作为定解条件的常微分方程求解问题称为常微分方程的边值问题(boundary-value problem). 为简明起见，我们以二阶边值问题为
则边值问题(8.1.1)有唯一解。
推论 若线性边值问题
()()()()()(),,
(),()y x p x y x q x y x f x a x b y a y b αβ'''=++≤≤??
==?
(8.1.2) 满足
(1) (),()p x q x 和()f x 在[,]a b 上连续； (2) 在[,]a b 上, ()0q x >, 则边值问题(8.1.1)有唯一解。
求边值问题的近似解，有三类基本方法：
(1) 差分法(difference method)，也就是用差商代替微分方程及边界条件中的导数，最终化为代数方程求解；
(2) 有限元法(finite element method)；
(3) 把边值问题转化为初值问题，然后用求初值问题的方法求解。
8.2 差分法
8.2.1 一类特殊类型二阶线性常微分方程的边值问题的差分法
设二阶线性常微分方程的边值问题为
(8.2.1)(8.2.2)
()()()(),,(),(),
y x q x y x f x a x b y a y b αβ''-=<
==?
其中(),()q x f x 在[,]a b 上连续，且()0q x ≥.
用差分法解微分方程边值问题的过程是：
(i) 把求解区间[,]a b 分成若干个等距或不等距的小区间，称之为单元；
(ii) 构造逼近微分方程边值问题的差分格式. 构造差分格式的方法有差分法, 积分插值法及变分插值法；本节采用差分法构造差分格式；
(iii) 讨论差分解存在的唯一性、收敛性及稳定性；最后求解差分方程. 现在来建立相应于二阶线性常微分方程的边值问题(8.2.1), (8.2.2)的差分方程. ( i ) 把区间[,]I a b =N 等分，即得到区间[,]I a b =的一个网格剖分：
011N N a x x x x b -=<<
<<=,
其中分点(0,1,
,)i x a ih i N =+=，并称之为网格节点(grid nodes)；步长b a N
h -=. ( ii ) 将二阶常微分方程(8.2.2)在节点i x 处离散化：在部节点
(1,2,
,1)i x i N =-处用数值微分公式
2(4)
1112
()2()()()(),12
i i i i i i i i y x y x y x h y x y x x h ξξ+---+''=-<< (8.2.3) 代替方程(8.2.2)中()i y x ''，得
112
()2()()
()()()()i i i i i i i y x y x y x q x y x f x R x h +--+-=+,
(8.2.4) 其中2(4)
()()12
i i h R x y ξ=. 当h 充分小时，略去式(8.2.4)中的()i R x ，便得到方程(8.2.1)的近似方程
11
2
2i i i i i i y y y q y f h +--+-=,
(8.2.5) 其中(),()i i i i q q x f f x ==, 11,,i i i y y y +-分别是11(),(),()i i i y x y x y x +-的近似值, 称式(8.2.5)为差分方程(difference equation)，而()i R x 称为差分方程(8.2.5)逼近方程
(8.2.2)的截断误差(truncation error). 边界条件(8.7.2)写成
0,.
N y y αβ==
(8.2.6)
于是方程(8.2.5), (8.2.6)合在一起就是关于1N +个未知量01,,,N y y y ，以及1N +个
方程式的线性方程组：
22112122
11222111(2),(2),1,2,,1,(2).i i i i i N N N N q h y y h f y q h y y h f i N y q h y h f αβ-+----?-++=-?-++==-??-+=-?
(8.2.7)
这个方程组就称为逼近边值问题(8.2.1), (8.2.2)的差分方程组(system of difference equations)或差分格式(difference scheme)，写成矩阵形式
2211122
2222233322222
2
2
111(2)1
1(2)11(2)1
1(2)11(2)N N N N N N y q h h f y q h h f y q h h f y q h h f y q h h f αβ------????-+-????????-+????
????????-+=?????????
?????????-+??????-+-????????????
. (8.2.8)
用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.7)或(8.2.8), 其解
01,,,N y y y 称为边值问题(8.2.1), (8.2.2)的差分解(difference solution). 由于
(8.2.5)是用二阶中心差商代替方程(8.2.1)中的二阶微商得到的，所以也称式(8.2.7)为中心差分格式(centered-difference scheme).
( iii ) 讨论差分方程组(8.2.7)或(8.2.8)的解是否收敛到边值问题(8.2.1), (8.2.2)的解，估计误差.
对于差分方程组(8.2.7)，我们自然关心它是否有唯一解；此外，当网格无限加密，或当0h →时，差分解i y 是否收敛到微分方程的解()i y x . 为此介绍下列极值原理：
定理8.2.1 (极值原理) 设01,,
,N y y y 是给定的一组不全相等的数，设
11
2
2(),0,1,2,
,i i i i i i i y y y l y q y q i N h +--+=
-≥=.
(8.2.9)
(1) 若()0,1,2,,i l y i N ≥=, 则{}0N
i i y =中非负的最大值只能是0y 或N y ； (2) 若()0,1,2,
,i l y i N ≤=, 则{}0N
i i y =中非正的最小值只能是0y 或N y .
证 只证(1) ()0i l y ≥的情形，而(2) ()0i l y ≤的情形可类似证明. 用反证法. 记{}0max i i N
M y ≤≤=，假设0M ≥, 且在121,,
,N y y y -中达到. 因为i y 不
全相等，所以总可以找到某个00(11)i i N ≤≤-，使0i y M =，而01i y -和01i y +中至少有一个是小于M 的. 此时
000000
0011
2
2
2()2.i i i i i i i i y y y l y q y h M M M q M q M h
+--+=
--+<-=-
因为00,0i q M ≥≥，所以0()0i l y <, 这与假设矛盾，故M 只能是0y 或N y . 证毕！
推论 差分方程组(8.2.7)或(8.2.8)的解存在且唯一. 证明 只要证明齐次方程组
112
02()0,0,1,2,,1,
0,0
i i i i i i i N y y y l y q y q i N h y y +--+?
=-=≥=-??
?==? (8.2.10)
只有零解就可以了. 由定理8.7.1知，上述齐次方程组的解01,,
,N y y y 的非负的最大值
和非正的最小值只能是0y 或N y . 而00,0N y y ==，于是0,1,2,
,.i y i N == 证毕！
利用定理8.2.1还可以证明差分解的收敛性及误差估计. 这里只给出结果： 定理8.2.2 设i y 是差分方程组(8.2.7)的解，而()i y x 是边值问题(8.2.1), (8.2.2)的解()y x 在i x 上的值，其中0,1,
,i N =. 则有
2
24()(),96i i i M h y x y b a ε=-≤-
(8.2.11)
其中(4)
4max ()a x b
M y
x ≤≤=.
显然当0h →时，()0i i i y x y ε=-→. 这表明当0h →时，差分方程组(8.2.7)或(8.2.8)的解收敛到原边值问题(8.7.1), (8.7.2)的解.
例8.2.1 取步长0.1h =，用差分法解边值问题
43,
01,
(0)(1)0,
y y x x y y ''-=≤≤??
==?
并将结果与精确解(
)()2222()3434x
x y x e
e e e x --=---进行比较.
解 因为110N h ==，()4,()3q x f x x ==, 由式(8.2.7)得差分格式
221222
1122
89(240.1)30.10.1,(240.1)30.1,2,3,,8,(240.1)30.10.9,
i i i i y y y y y x i y y -+?-+?+=???-+?+=?=??-+?=???
0100y y ==, 00.1,1,2,
,9i x ih i i =+==, 其结果列于表8.2.1.
从表8.2.1可以看出, 差分方法的计算结果的精度还是比较高的. 若要得到更精确的数值解，可用缩小步长h 的方法来实现.
8.2.2 一般二阶线性常微分方程边值问题的差分法
对一般的二阶微分方程边值问题
1212
()()()()()(),,()(),
()(),
y x p x y x q x y x f x a x b y a y a y b y b αααβββ'''++=<
'+=??'+=? (8.2.12) 假定其解存在唯一.
为求解的近似值，类似于前面的做法，
( i ) 把区间[,]I a b =N 等分，即得到区间[,]I a b =的一个网格剖分：
011N N a x x x x b -=<<
<<=,
其中分点(0,1,
,)i x a ih i N =+=，步长b a N
h -=. ( ii ) 对式(8.2.12)中的二阶导数仍用数值微分公式
2(4)
1112()2()()()(),12i i i i i i i i
y x y x y x h y x y x x h ξξ+---+''=-<<
代替，而对一阶导数，为了保证略去的逼近误差为2
()O h ，则用3点数值微分公式；另外为了保证插，在2个端点所用的3点数值微分公式与网格点所用的公式不同，即
2
1112012000022
212()()()(),,1,2,,1,263()4()()()(),,23()4()3()()(),.23i i i i i i i N N N N N N N N y x y x h y x y x x i N h y x y x y x h y x y x x h y x y x y x h y x y x x h ξξξξξξ+-----?-''''=-<<=-??
-+-?
''''=
+<
?-+''''=
+<
(8.2.13)
略去误差，并用()i y x 的近似值i y 代替()i y x ，(),(),()i i i i i i p p x q q x f f x ===，便得到差分方程组
1
111221*********
(2)(),1,2,,1,
2(34),
2(43),2i i i i i i i i i N N N N p y y y y y q y f i N h h
y y y y h y y y y h αααβββ-++---?-++-+==-??
?
+-+-=??
?+-+=??
(8.2.14)
其中(),(),(),1,2,,1i i i i i i q q x p p x f f x i N ====-, i y 是()i y x 的近似值. 整理
得
120212222
1122
2121(23)42,(2)2(2)(2)2,1,2,,1,4(32)2.
i i i i i i i N N N h y y y h hp y h q y hp y h f i N y y h y h αααααβββββ-+---+-=??---++==-??-++=? (8.2.15)
解差分方程组(8.2.15)，便得边值问题(8.2.12)的差分解01,,,N y y y .
特别地, 若12121,0,1,0ααββ====，则式(8.2.12)中的边界条件是第一类边值
条件：(),();y a y b αβ==此时方程组(7.7.16)为
2211121122
11221211112(2)(2)2(2),
(2)2(2)(2)2,2,3,,2,(2)2(2)2(2).i i i i i i i N N N N N N h q y hp y h f hp hp y h q y hp y h f i N hp y h q y h f hp αβ-+------?--++=--?---++==-??---=-+?
(8.2.16)
方程组(8.2.16)是三对角方程组，用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.16)，便得边值问题(8.2.12)的差分解01,,
,N y y y .
( iii ) 讨论差分方程组(8.2.16)的解是否收敛到微分方程的解，估计误差. 这里就不再详细介绍.
例8.2.2 取步长/16h π=，用差分法求下列边值问题的近似解，并将结果与精确解进行比较
.
精确解是1
()(sin 3cos )10
y x x x =-
+. 解 因为(20)8N h π=-=，()1,()2,()cos p x q x f x x =-=-=, 由式(8.2.17)得差分格式
()()()()()()()()()()()()()21222
112
2
2122216(2)216(1)216cos 16216(1)(0.3),216(1)2216(2)216(1)216cos 16,2,3,,6,216(1)2216(2)216cos 7i i i N N y y
y y y i i y y πππππππππππππ-+--??--?-++?-?????
?=--?-?-??????-?---?-++?-??????????==??-?---?-??????=()()16216(1)(0.1),ππ???
?
???????
?-+?-?-?????
080.3,0.1y y =-=-, 00.1,1,2,
,9i x ih i i =+==, 其结果列于表8.2.2.
8.3 有限元法
有限元法(finite element method)是求解微分方程定解问题的有效方法之一，它特别适用在几何、物理上比较复杂的问题. 有限元法首先成功地应用于结构力学和固体力学，以后又应用于流体力学、物理学和其他工程科学. 为简明起见，本节以线性两点边值问题为例介绍有限元法.
考虑线性两点边值问题
()(8.3.1)(8.3.2)
()()()()(),,(),(),
Ly p x y x q x y x f x a x b y a y b αβ?''?=-+=≤≤?==??
其中1
()0,()0,C [,]p x q x p a b >≥∈, ,C[,]q f a b ∈.
此微分方程描述了长度为b a -的可变交叉截面(表示为()q x )的横梁在应力()p x 和
()f x 下的偏差()y x .
8.3.1 等价性定理
记{
}
2
2
1C [,]()C [,],(),()a b y y y x a b y a y b αβ==∈==, 引进积分
()2
2()()[()]
()()2()()d b
a
I y p x y x q x y x f x y x x '=+-?
. (8.3.3)
任取2
1()C [,]y y x a b =∈，就有一个积分值()I y 与之对应，因此()I y 是一个泛函
(functional)，即函数的函数. 因为这里是,y y '的二次函数，因此称()I y 为二次泛函.
对泛函(8.3.3)有如下变分问题(variation problem)：求函数2
1C [,]y a b ∈，使得对任意2
1C [,]y a b ∈, 均有
()()I y I y ≥, (8.3.4)
即()I y 在y 处达到极小, 并称y 为变分问题(8.3.4)的解.
可以证明：
定理8.3.1(等价性定理) y 是边值问题(8.3.1), (8.3.2)的解的充分必要条件是y 使泛函()I y 在21C [,]a b 上达到极小，即y 是变分问题(8.3.4)在2
1C [,]a b 上的解.
证 (充分性) 设21C [,]y a b ∈是变分问题()I y 的解；即y 使泛函()I y 在2
1C [,]a b 上达到极小，证明y 必是边值问题(8.3.1), (8.3.2)的解.
设()x η是2
C [,]a b 任意一个满足()()0a b ηη==的函数，则函数
2
1()()()C [,]y x y x x a b αη=+∈，
其中α为参数. 因为y 使得()I y 达到极小，所以()()I y I y αη+≥，即积分
()2
2()()[()()]
()[()()]2()[()()]b
a
I y p x y x x q x y x x f x y x x dx
αηαηαηαη''+=+++-+?
作为α的函数，在0α=处取极小值()I y ，故
d
()0d I y ααηα=+=. (8.3.5) 计算上式，得
()()()()
()0
2
2
(8.d
()d d ()[()()]()[()()]2()[()()]d d 2()[()()]()2()[()()]()2()()d 2()()()()()()()()d .
b
a
b a
b
a
I y p x y x x q x y x x f x y x x x p x y x x x q x y x x x f x x x p x y x x q x y x x f x x x ααααηααηαηαηα
αηηαηηηηηη===+''=+++-+'''=+++-''=+-???
3.6)
利用分部积分法计算积分
[][]()()()d ()()d ()
()()()()()()d ()()()d ,
b
b
a
a
b b
a a
b
a
p x y x x x p x y x x p x y x x x p x y x x x p x y x x ηηηηη'''='''=-''=-?
???
代入式(8.3.6)，得
()0
(8.3.7)d
()2()()()()()()d 0.
d b a I y p x y x q x y x f x x x ααηηα'=????'+=-+-=?
因为()x η是任意函数，所以必有
()()()()()()0p x y x q x y x f x ''-+-≡. (8.3.8)
否则，若在[,]a b 上某点0x 处有
()00000()()()()()0p x y x q x y x f x ''-+-≠，
不妨设
()00000()()()()()0p x y x q x y x f x ''-+->，
则由函数的连续性知，在包含0x 的某一区间00[,]a b 上有
()()()()()()0p x y x q x y x f x ''-+->.
作
0022
00000,
[,]\[,],()()(),.
x a b a b x x a x b a x b η∈??=?--≤≤?? 显然2
()C [,]x a b η∈，且()()0a b ηη==，但
()()()()()()()d b
a p x y x q x y x f x x x η??''-+-????
?
()0
0()()()()()()d 0b a p x y x q x y x f x x x η??''=-+->????
?，
这与式(8.3.7)矛盾. 于是式(8.3.8)成立，即变分问题(8.3.4)的解y 满足微分方程(8.3.1), 且(),()y a y b αβ==故它是边值问题(8.3.1), (8.3.2)的解.
(必要性) 设()y y x =是边值问题(8.3.1), (8.3.2)的解，证明y 是变分问题(8.3.4)的解；即：y 使泛函()I y 在2
1C [,]a b 上达到极小.
因为()y y x =满足方程(8.3.1)，所以()()()()()()0p x y x q x y x f x ''-+≡.
设任意2
1()C [,]y y x a b =∈，则函数()()()x y x y x η=-满足条件()()0a b ηη==，
且2
()C [,]x a b η∈. 于是
()()[]()2
2
2
2
2
2
()()()()
()[()()]()[()()]2()[()()]d ()[()]()[()]2()()d 2()()()()()()()()d ()[()]()[()]d b
a
b
a b
a
a
I y I y I y I y p x y x x q x y x x f x y x x x p x y x q x y x f x y x x
p x y x x q x y x x f x x x p x x q x x x
ηηηηηηηηη-=+-''=+++-+'-+-''=+-++?
??()()()222
22()()()()()()d ()[()]()[()]d ()[()]
()[()]d .
b
b b
a a b
a
p x y x q x y x f x x x p x x q x x x
p x x q x x x ηηηηη??'''=--+++????
'=+????
因为()0,()0p x q x >≥，所以当()0x η≠时，
()2
2()[()]
()[()]d 0b
a
p x x q x x x ηη'+>?, 即
()()0I y I y ->.
只有当()0x η≡时，()()0I y I y -=. 这说明y 使泛函()I y 在2
1C [,]a b 上达到极小. 证毕！
定理8.3.2 边值问题(8.3.1), (8.3.2)存在唯一解.
证明 用反证法. 若12(),()y x y x 都是边值问题(8.3.1), (8.3.2)的解，且不相等，则由定理8.3.1知，它们都使泛函()I y 在2
1C [,]a b 上达到极小，因而
12()()I y I y > 且 21()()I y I y >，
矛盾！因此边值问题(8.3.1), (8.3.2)的解是唯一的.
由边值问题解的唯一性，不难推出边值问题(8.3.1), (8.3.2)解的存在性(留给读者自行推导).
8.3.2 有限元法
等价性定理说明，边值问题(8.3.1), (8.3.2)的解可化为变分问题(8.3.4)的求解问题. 有限元法就是求变分问题近似解的一种有效方法. 下面给出其解题过程：
第1步 对求解区间进行网格剖分
01,i n a x x x x b =<<<<<=
区间1[,]i i i I x x -=称为单元，长度1(1,2,
,)i i i h x x i n -=-=称为步长，1max i i n
h h ≤≤=. 若
(1,2,
,)i h h i n ==，则称上述网格剖分为均匀剖分.
给定剖分后，泛函(8.3.3)可以写成
()2
2()()[()]
()()2()()d b
a
I y p x y x q x y x f x y x x '=+-?
()1
2
2
1
1
()[()]
()()2()()d i
i n
n
x i x i i p x y x q x y x f x y x x
S -=='=+-∑∑?
记
. (8.3.9)
第2步 构造试探函数空间。为了计算积分(8.3.9)，最简单的近似方法是将分段线性函数的集合作为试探函数空间。设01,,,n y y y 分别是边值问题(8.3.1), (8.3.2)的解()
y x 在节点01,,
,n x x x 处的值，用分段线性插值
1
11(),[,].i i i i i i i i i
x x x x y y x y y x I x x h h -----==
+∈= (8.3.10) 近似1()[,]i i i y x x I x x -∈=，
，并称式(8.3.10)为单元形状函数(element shape function). 为了将线性插值(8.3.10)标准化，令
1
i i
x x t h --=
， 则将1[,]i i i I x x -=变到t 轴上的单元[0,1]. 记01()1,()N t t N t t =-=，则式(8.3.10)可写成
011()(),[0,1].i i y N t y N t y t -=+∈ (8.3.11)
第3步 建立有限元方程组. 将式(8.3.10)代入泛函(8.3.9)，有
()1
1
2
2
1
1
()()()[()]()[()]d 2()()d i
i
i i n
n
x x x x i i I y I y p x y x q x y x x f x y x x --=='≈=+-∑∑?
?
.
由式(8.3.11)知
2
1
2111011011
10110
11
122111
1
1
1()
()()()[()()]d 2()[()()]d ()()(1)2()({n
i i i i i i i i i i i n
i i i i i i n
i i i i i i i i i i i i y y I y p x th h q x th N t y N t y t
h h f x th N t y N t y t
h p x th h q x th t y h p x th h q x ----=--=----=--??-=++++ ???
-++??=+++-??+-++∑?∑?∑?11122
11)(1)()()d }i i i i i i i i i i i th t t y y h p x th h q x th t y t
-----??+-????++++??
1
110
1
2()[(1)]d .n
i i i i i i h f x th t y ty t --=-+-+∑? (8.3.12)
式(8.3.12)右端第1个求和号的第i 项(对应第i 个单元)是关于1,i i y y -的二次型，可以写成
()T ()()()i i i Y K Y , (8.3.13)
其中
()T 1(,)i i i y y -=Y ，()()
1,1
1,()
()(),1,i i i i i i
i i i i i i i K K K K -
---??= ? ???
K , (8.3.14) ()i K 称为单元刚度矩阵(element stiffness matrix)，
1()121,11101()12
,1101
()()11,,1110()()(1)d ,
()()d ,
()()(1)d .i i i i
i i i i i i i i i i i i i i i i i i i i i i i i i i K h p x th h q x th t t K h p x th h q x th t t K K h p x th h q x th t t t -------------???=+++-??????=+++???????==-+++-???
??? (8.3.15) 而式(8.3.12)的第2个求和号的第i 项可以写成
()T ()()i i b Y (8.3.16)
其中
()()()T
12(,)i i i b b =b ，()T 1(,)i i i y y -=Y ，
1
1
()()1
12
10
()(1)d ,()d .i i i i i i i i b h f x th t t b h f x th t t --=+-=+??
于是式(8.3.12)中求和号的项可以写成
()T ()()()T ()()2()i i i i i i S =-Y K Y b Y . (8.3.17)
再令T 01(,,
,)n y y y =Y 及2(1)n ?+矩阵
()1001000,0
0010
0i i i -??=
???
列列
C 则()
()i i =Y
C Y .
于是式(8.3.17)又可写成
()T ()()()T ()()()2()()i i i i i i S =-C Y K C Y b C Y
(
)
T
T
()T
()
()
()T
()()2()i i i i i =-Y C K C Y C b
Y . (8.3.18)
把式(8.3.18)代入式(8.3.12)右端求和号，得
()T T
()T ()()()T ()111()()2()n
n n i i i i i i i i i I y S ===????==- ? ?????
∑∑∑Y C K C Y C b Y . (8.3.19)
记
()
()
1,1
1,()T
()
()
()()1
1,1,0
000000000
000000()00000000000
0i i n
n i i i i
i i i i i i i i i i i K K K K ---==-??
?????????
?==??????????????
∑∑K C K C
(1)
(1)0001(1)(1)(2)(2)101111
12(2)(2)(3)(3)21
2222
23(3)(3)(4)(4)
32
3333
34
(1)(1)()
()
1,2
1,11,1
1,()(),1
,n n n n n n n n n n n n n n n n n n K K K K K K K K K K K K K K K K K K K K ----------????+????+??=+??????+?????
?
， (8.3.20)
()T ()()()
T 121
1
1,()(0,
,0,,,0,
,0)n n
i i i i i i i i
b b ==-==∑∑b C b
(1)(1)(2)(2)(3)(3)
(4)(1)()()T
1212121212(,,,,,,)n n n b b b b b b b b b b -=++++, (8.3.21)
则式(8.3.19)化为
T T ()2.I y =-Y KY b Y (8.3.22)
并称K 为总刚度矩阵(total stiffness matrix)，称b 为右端向量.
因为()y x 是使()I y 取极小值的函数，所求的01,,,n y y y 自然使式(8.3.22)的右边取
极小值，所以应有
T T d
(2)0,0,1,2,,d i
i n y -==Y KY b Y . (8.3.23)
记T 12(),(,,
,)ij n n n K b b b ?==K b , 则式(8.8.23)为
0000
d
2d 2n n n r j j r r r r j r i
n n
r i r i j j i
r j K y y b y y K y K y b =====????
-??
???????=+-∑∑∑∑∑
020,1,2,,1,n ir r i r K y b i n =??
=-==- ???
∑
或
,0,1,
,,n
ir
r i r K
y b i n ===∑ (8.3.24)
得方程组
=KY b . (8.3.25)
因为0,n y y αβ==是已知的，不能任意选取，所以不能要求式(8.3.23)对0,n y y 也成立. 因此方程组(8.3.24)或(8.3.25)中应当去掉首末2个方程，并把其他方程中含有的0,n y y 改为已知量，所得方程组就是未知量121,,
,n y y y -满足的代数方程组
(1)(2)
(2)1111112(2)(2)(3)(3)
2212222
23
(2)(2)(1)
(1)
22,32,22,2
2,1(1)
(1)()11,2
1,11,1n n n n n n n n n n n n n n n n n n n n n n n y K K K y K K K K y K K K K y K K K ----------------------??+??
????+????????????+????????+??
?
? (1)(2)(1)
2110(2)(3)
21(2)(1)
21
(1)()()2
11,n n n n n n n b b K b b b b b b K αβ----??+-??+??
?
?=??+????
+-??. (8.3.26) 方程组(8.3.25)或(8.3.26)称为有限元方程组(finite element system of equations).
用第2章介绍的解三对角方程组的追赶法求解有限元方程组(8.3.26)，可解出
121,,,n y y y -，即得变分问题(8.3.4)的近似解，也就是边值问题(8.3.1), (8.3.2)解的近
似值. 这种方法称为有限单元法(finite element method)或有限元法.
例8.3.1 用有限元法求下列边值问题的近似解，并将结果与精确解进行比较.
2(())4()481,01,
(0)(1)0,
xy x y x x x x y y ''?-+=-+≤≤?
==? 取0.2h =，精确解是2
()y x x x =-.
解 因为15N h ==，2
(),()4,()481p x x q x f x x x ===-+, 由式(8.3.26)得有限元方程组
12
123
234342.53333 1.366670.082667,1.36667 4.53333 2.366670.306667,2.36667 6.53333 3.366670.466667,3.366678.533330.562667.y y y y y y y y y y -=-??-+-=-??
-+-=-?
?-+=-?
其结果列于表8.3.1.
上面虽然是对边值问题(8.3.1), (8.3.2)介绍的有限元解法，但其解题步骤对一般的微分方程定解问题也是适用的.
对所给微分方程定解问题，首先找出相应微分方程的变分形式，然后进行下列步骤： 第1步 对定义区域(或定义区间)进行网格剖分, 将定义区域(或定义区间)剖分为若干个小单元(一维是区间；多维是区域，如矩形、三角形等)；
第2步 构造试探空间; 即选择适当的插值函数类.
第3步 建立有限元方程组. 计算单元刚度矩阵及右端向量，再形成总刚度矩阵及总右端项，写出有限元方程组；结合定解条件，求解有限元方程组.
注 从形式上看，用有限元法解微分方程定解问题很繁，但是有限元法的求解步骤规，便于在计算机上实现，并且总刚度矩阵是三对角对称正定矩阵，因此有限元方程组可用第2章介绍的解三对角方程组的追赶法求解. 有限元法最主要的优点是可以处理相当复杂的区域上的边值问题。
8.4 打靶法
解常微分方程边值问题的打靶法(shooting method)，也称为尝试法，其基本思想是把边值问题转化为初值问题来解：首先作出一些只满足一端边值条件的解，然后再从这些解中找出适合另一端边值条件的解. 下面以二阶常微分方程带第一类边界条件的边值问题
(8.4.1)(8.4.2)
()(,(),()),
,(),().
y x f x y x y x a x b y a y b αβ'''=<
==?
为例介绍常微分方程边值问题的打靶法.
7.6节曾介绍过二阶常微分方程初值问题(7.6.10)
()(,(),()),,
(),(),
y t f t y t y t a t b y a y a αα'''=<≤??
''==?
的求解方法. 将上式中的t 改为x ，α'改为s ，得
()(,(),()),,
(),(),
y x f x y x y x a x b y a y a s α'''=<≤??
'==? (8.4.3)
设初值问题(8.4.3)的解为1()y x ，显然1()y x 依赖于s ，即11()(,)y x y x s =. 为了求解边值问题(8.4.1), (8.4.2)，必须求s s =，使之满足11()(,)y b y b s β==.
下面介绍2种方法来求s .
方法1 根据实际问题情况选一个1s ，解初值问题
1
(,,),
,(),(),y f x y y a x b y a y a s α'''=<≤??
'==? (8.4.4) 得到的解仍记为1()y x .
若1()y b β=或1()y b βε-<(ε为事先给定的精度)，则1()y x 就是边值问题(8.4.1), (8.4.2)的解.
否则，根据11()y b β=与β的误差，将1s 修改为2s ；例如取211
s s β
β=, 再解初值问题
2()(,(),()),,
(),(),
y x f x y x y x a x b y a y a s α'''=<≤??
'==? (8.4.5)
得到解2()y x .
若22()y b β=满足2ββ=或2ββε-<，则2()y x 就是边值问题(8.4.1), (8.4.2)
的解.
否则，再将2s 适当修改为3s ; 例如可用1122(,),(,)s s ββ作线性插值求3s ：
21
31121
()s s s s ββββ-=+
--,
然后解初值问题
3()(,(),()),,
(),(),
y x f x y x y x a x b y a y a s α'''=<≤??
'==? (8.4.6)
的解. 如此继续下去，直到达到精度要求为止.
方法2 求s 另一种方法是求函数()(,)F s y b s β=-的一个零点s . 对于每一个自变量s ，通过解初值问题(8.4.4)，可解出(,)y x s . 计算
()(,)F s y b s β=-, (8.4.7)
然后采用第3章介绍的求方程根的方法求()F s 的零点.
首先寻找(1)
(2)
,s s
，使(1)(2)
()0,()0F s F s
<>，则在区间(1)(2)(,)s s 或(2)(1)(,)s s 至
少有()F s 的一个零点. 然后可用二分法求s . 也可用Newton 迭代公式
1()
()
k k k k F s s s F s +=-
' 求s 的近似值.
几何解释：微分方程边值问题(8.4.1), (8.4.2)的解()y x 是一条通过
(,()),(,())A a y a B b y b 两点的曲线(如图8.4.1). 初值问题(8.4.4)的解1(,)y x s 是一条
通过点(,())A a y a 、斜率为1()y a s '=的曲线(如图8.4.1). 初值问题(8.4.5)的解2(,)y x s 是一条通过点(,())A a y a 、斜率为2()y a s '=的曲线(如图8.4.1). 这有点象射击者从定点
A 向目标
B 射击一样. 根据经验以某一角度(斜率为1s )试射一次. 如果与目标相差太大，
调整射击角度(斜率为2s )，再进行射击. 如此继续进行下去，直到击中目标，或击中的点与B 的误差在允许的围之。
图8.4.1
参考文献[1]提供的资料还讨论了选择初始值()y a s =的重要性，并介绍了多重打靶法，这里就不作详细介绍，有兴趣的读者可参看相关资料.
8.5 算法程序
8.5.1 用差分法解二阶线性常微分方程的边值问题
%用差分法解二阶线性常微分方程的边值问题
%[a,b]是区间, Step 是步长, Alpha, Beta 是初值 % f, q 分别是式(8.2.1)中的f (x ), q (x )
function Y = DiffMethod(f, q, a, b, Step, Alpha, Beta) N = (b-a) / Step; X = a : Step : b; A = zeros( N-1 ); for i = 1 : N-1
A(i,i) = -1 * ( 2+feval(q, X(i+1)) * Step^2 ); if i ~= N-1
A(i,i+1) = 1; A(i+1,i) = 1; end end
B(1) = Step^2 * feval(f, X(2)) - Alpha; B(2:N-2) = Step^2 * feval(f, X(3:N-1)); B(N-1) = Step^2 * feval(f, X(N)) - Beta; B = B';
[L,U] = lu( A ); Y = U \ ( L \ B ) ; for i = 1 : length(Y)
sprintf('%10.7f',Y(i)) end end
例8.5.1 取0.1h =，用差分法求下列边值问题的近似解
()4(()),01,
(0)0,(1) 2.y x y x x x y y ''=-≤≤??
==?
解 在MATLAB 命令窗口输入：
f = inline('-4.*x'); q= inline('4'); DiffMethod(f, q, 0, 1, 0.1, 0, 2)
回车，可的结果。
8.5.2 用差分法解一般二阶常微分方程的边值问题
%用差分法解一般二阶常微分方程的边值问题
% [a,b]是区间, Step 是步长, Alpha, Beta 是初值 % f, p, q 分别是式(8.2.12)中的f (x ), p (x ), q (x )
function Y = DiffMethod_2(f, p, q, a, b, Step, Alpha, Beta) N = (b-a) / Step; X = a : Step : b; A = zeros( N-1 );
A(1, 1) = -2 * ( 2 - Step^2 * feval(p, X(2)) ); A(1,2) = 2 + Step * feval(f, X(3)); for i = 2 : N-2
A(i,i) = -2 * ( 2 - feval(p, X(i+1)) * Step^2 ); A(i,i-1)= 2 - Step * feval(f, X(i+1)); A(i-1,i) = 2 + Step * feval(f, X(i+1)); end
A(N-1,N-2) = 2 - Step * feval(f, X(N));
A(N-1,N-1) = -2 * ( 2 - feval(p, X(N)) * Step^2 );
B(1) = 2*Step^2 * feval(q, X(2)) - ( 2 - Step * feval(f, X(2)) ) * Alpha; B(2:N-2) = 2 * Step^2 * feval(q, X(3:N-1));
B(N-1) = 2*Step^2 * feval(q, X(N)) - ( 2 + Step * feval(f, X(N)) ) * Beta; B = B';
[L,U] = lu( A ); Y = U \ ( L \ B ) ; for i = 1 : length(Y)
sprintf('%10.7f',Y(i)) end end
例8.5.2 取/16h π=，用差分法求下列边值问题的近似解
解 在MATLAB 命令窗口输入：
f = inline('cos(x)'); p = inline('-1'); q = inline('-2'); DiffMethod_2(f, p, q, 0, pi/2, pi/16, -0.3, -0.1) 回车，可得结果。
8.5.3 用有限元法解二阶常微分方程的边值问题
%用有限元法解二阶常微分方程的边值问题 %[a,b]是区间, h 是步长, Alpha, Beta 是初值
% f, p, q 分别是式(8.3.1)中的f (x ), p (x ), q (x )
function Y = FiniElem(f, p, q, a , b, Step, Alpha, Beta) N = length(Step) - 1; X(1) = a;
for i = 1 : N+1
X(i+1) = X(i) + Step(i); end syms t;
ff = -1/Step(N+1) * feval(f, X(N)+t*Step(N+1)) + ...
Step(N+1)* feval(p, X(N)+t*Step(N+1)) * (1-t)*t;
Knnn = int(ff, t, 0, 1);
syms t;
ff = -1/Step(2) * feval(f, X(1)+t*Step(2)) + ...
Step(2)*feval(p, X(1)+t*Step(2))*(1-t)*t;
K101 = int(ff, t, 0, 1);
for i = 2 : N
syms t;
ff = Step(i+1) * feval(q, X(i)+t*Step(i+1)) * (1-t);
B1(i) = int(ff, t, 0, 1);
end
for i = 1 : N-1
syms t;
ff = Step(i+1) * feval(q, X(i)+t*Step(i+1)) * t;
B2(i) = int(ff, t, 0, 1);
end
B(1) = B2(1) + B1(2) - Alpha*K101;
B(N-1) = B2(N-1) + B1(N) - Beta*Knnn;
for i = 2: N-2
B(i) = B2(i) + B1(i+1);
end
B=B';
A=zeros(N-1);
for i = 1 : N-1
syms t;
ff = 1/Step(i+1) * feval(f, X(i)+t*Step(i+1)) + ...
Step(i+1)* feval(p, X(i)+t*Step(i+1)) * t^2;
K1 = int(ff, t, 0, 1);
syms t;
ff = 1/Step(i+2) * feval(f, X(i+1)+t*Step(i+2)) + ...
Step(i+2) *feval(p, X(i+1)+t*Step(i+2)) * (1-t)^2;
K2 = int(ff, t, 0, 1);
A(i,i) = eval(K1+K2);
end
for i = 1 : N-2
syms t;
ff = -1/Step(i+2) *feval(f, X(i+1)+t*Step(i+2)) + ...
Step(i+2) * feval(p, X(i+1)+t*Step(i+2)) * (1-t) * t; K1 = int(ff, t, 0, 1);
A(i,i+1) = eval(K1);
A(i+1,i) = eval(K1);
end
[L,U]=lu(A);
Y = eval( U \ (L\B) );
小 结
本章介绍了解微分方程边值问题的差分方法、有限元法和打靶法。前两种方法是解微分方程边值问题最常用的方法，它们最终得到的方程组往往是三对角方程组，可通过追赶法来解。有限元法最主要的优点是可以处理相当复杂的区域上的边值问题。而且差分方法和有限元法还是求解偏微分方程的两类重要的方法.
习 题 八
1. 用差分法求下列边值问题的近似解，并将结果与精确解进行比较。 (1) ()4(),01,
(0)0,(1)2,
y x y x x y y ''=-≤≤??
==?
取0.1h =，精确解是222
2
()x x
e e y x x e e
---=+-.
(2) ()2()(),02,
(0)0,(2)4,
x y x y x y x xe x x y y '''?=-+-≤≤?==-?
取0.2h =，精确解是315()22
63
x x
x y x x e xe e x =
-+--. 2. 用有限元法求下列边值问题的近似解，并将结果与精确解进行比较。
取0.2h =，精确解是()()()
1132234()cos cos y x x x x πππ=-+.