常微分方程边值问题的数值解法
第8章 常微分方程边值问题的数值解法
引 言
第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.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), 的差分方程. ( i ) 把区间[,]I a b =N 等分，即得到区间[,]I a b =的一个网格剖分：
011N N a x x x x b -=<<<<=L ,
其中分点(0,1,,)i x a ih i N =+=L ，并称之为网格节点(grid nodes)；步长b a N
h -=. ( ii ) 将二阶常微分方程(8.2.2)在节点i x 处离散化：在内部节点
(1,2,,1)i x i N =-L 处用数值微分公式
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 ，便得到方程的近似方程
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 称为差分方程逼近方程的截断误差
(truncation error). 边界条件写成
0,.
N y y αβ==
(8.2.6)
于是方程(8.2.5), 合在一起就是关于1N +个未知量01,,,N y y y L ，以及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 αβ-+----?-++=-?-++==-??-+=-?
L (8.2.7) 这个方程组就称为逼近边值问题(8.2.1), 的差分方程组(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 αβ------????-+-????????-+????
????????-+=???
??????
?????????-+??????-+-????????????
M O O O M . (8.2.8)
用第2章介绍的解三对角方程组的追赶法求解差分方程组(8.2.7)或 其解01,,,N y y y L 称为边值问题 的差分解(difference solution). 由于是用二阶中心差商代替方程中的二阶微商得到的，所以也称式为中心差分格式(centered-difference scheme).
( iii ) 讨论差分方程组(8.2.7)或的解是否收敛到边值问题 的解，估计误差. 对于差分方程组(8.2.7)，我们自然关心它是否有唯一解；此外，当网格无限加密，或当0h →时，差分解i y 是否收敛到微分方程的解()i y x . 为此介绍下列极值原理：
定理8.2.1 (极值原理) 设01,,,N y y y L 是给定的一组不全相等的数，设
11
22(),0,1,2,,i i i i i i i y y y l y q y q i N h +--+=
-≥=L .
(8.2.9) (1) 若()0,1,2,,i l y i N ≥=L , 则{}0N
i i y =中非负的最大值只能是0y 或N y ； (2) 若()0,1,2,,i l y i N ≤=L , 则{}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 -L 中达到. 因为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)或的解存在且唯一. 证明 只要证明齐次方程组
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 +--+?
=-=≥=-??
?==?L (8.2.10) 只有零解就可以了. 由定理8.7.1知，上述齐次方程组的解01,,,N y y y L 的非负的最大值和非正的最小值只能是0y 或N y . 而00,0N y y ==，于是0,1,2,,.i y i N ==L 证毕！
利用定理8.2.1还可以证明差分解的收敛性及误差估计. 这里只给出结果： 定理8.2.2 设i y 是差分方程组的解，而()i y x 是边值问题 的解()y x 在i x 上的值，其中0,1,,i N =L . 则有
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.1 取步长0.1h =，用差分法解边值问题
43,
01,
(0)(1)0,
y y x x y y ''-=≤≤??
==?
并将结果与精确解()()222
2()3434x x
y x e e e
e x --=---进行比较.
解 因为110N h ==，()4,()3q x f x x ==, 由式(8.2.7)得差分格式
221222
112289(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 -+?-+?+=???-+?+=?=??-+?=???
L 0100y y ==, 00.1,1,2,,9i x ih i i =+==L , 其结果列于表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 -=<<<<=L ,
其中分点(0,1,,)i x a ih i N =+=L ，步长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 ξξξξξξ+-----?-''''=-<<=-??
-+-?
''''=
+<
?-+''''=
+<
L (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
1112210
01221211
(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 αααβββ-++---?-++-+==-??
?
+-+-=??
?+-+=??
L (8.2.14) 其中(),(),(),1,2,,1i i i i i i q q x p p x f f x i N ====-L , 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 αααααβββββ-+---+-=??---++==-??-++=?L (8.2.15) 解差分方程组(8.2.15)，便得边值问题的差分解01,,,N y y y L .
特别地, 若12121,0,1,0ααββ====，则式(8.2.12)中的边界条件是第一类边值条件：(),();y a y b αβ==此时方程组为
2211121122
1122
1211112(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 αβ-+------?--++=--?---++==-??---=-+?L
(8.2.16)
方程组(8.2.16)是三对角方程组，用第2章介绍的解三对角方程组的追赶法求解差分方程组，便得边值问题的差分解01,,,N y y y L .
( 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 πππππππππππππ-+--??--?-++?-?????
?=--?-?-??????-?---?-++?-??????????==??-?---?-??????=L ()()16216(1)(0.1),ππ???
?
???????
?-+?-?-?????
080.3,0.1y y =-=-, 00.1,1,2,,9i x ih i i =+==L , 其结果列于表8.2.2.
有限元法
有限元法(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 %是边值问题 的解的充分必要条件是y
%使泛函()I y 在21C [,]a b 上达到极小，即y %是变分问题在2
1C [,]a b 上的解.
证 (充分性) 设21C [,]y a b ∈%是变分问题()I y 的解；即y %使泛函()I y 在2
1C [,]a b 上
达到极小，证明y %必是边值问题(8.3.1), 的解.
设()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)矛盾. 于是式成立，即变分问题的解y %满足微分方程 且(),()y a y b αβ==%%故它是边值问题 的解.
(必要性) 设()y y x =%%是边值问题(8.3.1), 的解，证明y %是变分问题的解；即：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 边值问题 存在唯一解.
证明 用反证法. 若12(),()y x y x 都是边值问题(8.3.1), 的解，且不相等，则由定理知，它们都使泛函()I y 在2
1C [,]a b 上达到极小，因而
12()()I y I y > 且 21()()I y I y >，
矛盾！因此边值问题(8.3.1), 的解是唯一的.
由边值问题解的唯一性，不难推出边值问题(8.3.1), 解的存在性(留给读者自行推导).
8.3.2 有限元法
等价性定理说明，边值问题(8.3.1), 的解可化为变分问题的求解问题. 有限元法就是
求变分问题近似解的一种有效方法. 下面给出其解题过程：
第1步 对求解区间进行网格剖分
01,i n a x x x x b =<<<<<=L L
区间1[,]i i i I x x -=称为单元，长度1(1,2,,)i i i h x x i n -=-=L 称为步长，1max i i n
h h ≤≤=. 若
(1,2,,)i h h i n ==L ，则称上述网格剖分为均匀剖分.
给定剖分后，泛函(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 L 分别是边值问题 的解()y x 在节点
01,,,n x x x L 处的值，用分段线性插值
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)代入泛函，有
()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 ，
11
()()11210
()(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 =L Y 及2(1)n ?+矩阵
()1001000,000100i i i -??=
???
L L L
L
列列
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)代入式右端求和号内，得
()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
00i i n
n i i i i
i i i i i i i i i i i K K K K ---==-??
?????????
?==??????????????
∑∑L L M M M M M M L L L L L L L L M M M M M L
L
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 ----------??
??+????+??=+??????+?????
?
O
O O
， (8.3.20)
()T ()()()
T 121
1
1,()(0,,0,,,0,,0)n n
i i i i i i i i
b b ==-==∑∑L L 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 -=++++L , (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 L 自然使式(8.3.22)的右边取极小值，所以应有
T T d
(2)0,0,1,2,,d i
i n y -==L Y KY b Y . (8.3.23) 记T 12(),(,,,)ij n n n K b b b ?==L 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 =??
=-==- ???
∑L
或
,0,1,,,n
ir
r i r K
y b i n ===∑L (8.3.24)
得方程组
=KY b . (8.3.25)
因为0,n y y αβ==是已知的，不能任意选取，所以不能要求式(8.3.23)对0,n y y 也成立. 因此方程组或中应当去掉首末2个方程，并把其他方程中含有的0,n y y 改为已知量，所得方程组就是未知量121,,,n y y y -L 满足的代数方程组
(1)(2)
(2)1111112(2)(2)(3)(3)2212222
23(2)(2)(1)
(1)
22,3
2,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 ----------------------??+??
????+????????????+????????+??
?
?M O
O O
(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 αβ----??+-??+??
??=??+????
+-??M
. (8.3.26)
方程组(8.3.25)或称为有限元方程组(finite element system of equations).
用第2章介绍的解三对角方程组的追赶法求解有限元方程组(8.3.26)，可解出
121,,,n y y y -L ，即得变分问题的近似解，也就是边值问题 解的近似值. 这种方法称为有限
单元法(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), 介绍的有限元解法，但其解题步骤对一般的微分方程定解问题也是适用的.
对所给微分方程定解问题，首先找出相应微分方程的变分形式，然后进行下列步骤： 第1步 对定义区域(或定义区间)进行网格剖分, 将定义区域(或定义区间)剖分为若干个小单元(一维是区间；多维是区域，如矩形、三角形等)；
第2步 构造试探空间; 即选择适当的插值函数类.
第3步 建立有限元方程组. 计算单元刚度矩阵及右端向量，再形成总刚度矩阵及总右端项，写出有限元方程组；结合定解条件，求解有限元方程组.
注 从形式上看，用有限元法解微分方程定解问题很繁，但是有限元法的求解步骤规范，便于在计算机上实现，并且总刚度矩阵是三对角对称正定矩阵，因此有限元方程组可用第2章介绍的解三对角方程组的追赶法求解. 有限元法最主要的优点是可以处理相当复杂的区域上的边值问题。
打靶法
解常微分方程边值问题的打靶法(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.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 =. 为了求解边
值问题 ，必须求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), 的解.
否则，根据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), 的解.
否则，再将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), 的解()y x 是一条通过(,()),(,())A a y a B b y b 两点的曲线(如图 初值问题的解1(,)y x s 是一条通过点(,())A a y a 、斜率为1()y a s '=的曲线(如图 初值问题的解2(,)y x s 是一条通过点(,())A a y a 、斜率为2()y a s '=的曲线(如图 这有点象射击者从定点A 向目标B 射击一样. 根据经验以某一角度(斜率为1s )试射一次. 如果与目标相差太大，调整射击角度(斜率为2s )，再进行射击. 如此继续进行下去，直到击中目标，或击中的点与B 的误差在允许的范围之内。
图8.4.1
参考文献[1]提供的资料还讨论了选择初始值()y a s =的重要性，并介绍了多重打靶法，这里就不作详细介绍，有兴趣的读者可参看相关资料.
算法程序
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('%',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, 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('%',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, , 回车，可得结果。
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
小 结
本章介绍了解微分方程边值问题的差分方法、有限元法和打靶法。前两种方法是解微分方程边值问题最常用的方法，它们最终得到的方程组往往是三对角方程组，可通过追赶法来解。有限元法最主要的优点是可以处理相当复杂的区域上的边值问题。而且差分方法和有限元法还是求解偏微分方程的两类重要的方法.
习 题 八
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()2263
x x
x y x x e xe e x =
-+--. 2. 用有限元法求下列边值问题的近似解，并将结果与精确解进行比较。
取0.2h =，精确解是(
)()()
1132234
()cos cos y x x x x πππ=-+. (2) ()()()(2),01,
(0)(1)0,
x x x e y x e y x x x e x y y ?
''?-+=+-≤≤??==?
取0.1h =，精确解是()(1)(1)x
y x x e
-=--.
3. 证明：通过某种变量代换，可将边值问题
()()()()()(),01,
(0),(1)p x y x q x y x f x x y y αβ
?''?-+=≤≤?
==?? 化为下列形式的边值问题
常微分方程边值问题与不动点定论文
目录 引言 (1) 1预备知识 (2) 定义1.1（奇异Sturm-Liouville边值问题的正解） (2) 引理1.1.1 (2) 定义1.2（凸集的概念) (3) 定义1.3锥的定义 (3) 定义1.4（全连续算子的概念） (3) 1.5 （常微分边值问题的定义） (4) 定义1.6混合单调算子得定义） (4) 2 常微分方程边值问题正解得存在性 (5) 2.1 奇异Sturm-Liouville常微分边值问题的正解存在在 (5) 子 (8) 2.2 一类二阶边值问题的存在性 (9) 3一类混合单调算子应用 (11) 3.1一类混合单调算子的存在唯一性?........................ 错误!未定义书签。 3.2 求常微分边值问题的例题 (13) 结束语 (15) 参考文献 (15) 致 (16)
常微分方程边值问题与不动点定 （数学与统计学院 11级数学与应用数学2班）指导教师：攀峰 引言 从历史上看在有了微积分这个概念以后，紧接着出现了常微分方程。发展初期是属于“求通解”得时代，当人们从初期的热潮中结束要从维尔证明了卡帝方程中是一定不会存在一般性的初等解的时候开始的，并且柯西紧接着又提出了初值问题，常微分方程开始从重视“求通解”转向重视“求定解”的历史时代。 大学我们都学习了常微分方程这门学科，如果要研究它的定解问题，我们首先就会知道是常微分方程的初值问题。然而，在科学技术、生产实际问题中，我们还是提出了另一类定解问题-边值问题。对于常微分方程边值问题，伟大的科学家最早在解决二阶线性微分方程时，提出了分离变量法。[]1.在牛顿时期，科学家们已经提出过常微分的边值问题，牛顿也对常微分边值问题进行过研究，并且在1666年10月牛顿已经在这个领域取得了很大的成就，但是由于种种原因当时并没有整理成论文，所以没有及时出版。但在1687年他终于把在常微分方程上研究的成果发表了，虽然不是在数学著作中，却是他的一本力学著作中（《自然哲学的数学原理》）。 在微积分刚创立时期，雅克.伯努利来自瑞士的科学家提出了远著文明的问题-悬链线问题，紧着的地二年著名数学家莱布尼兹就给出了正确的解答，通过对绳子上个点受力分析，建立了以下方程 这个方程满足的定解条件是y（a）=α；y(b)=β.这是一个典型的常微分方程的边值问题。从这开始，常微分边值问题已经是科学家研究微分方程是不可或缺的工具，我就简单列举几个例子：（比如种族的生态系统；梁的非线性震动）等。对于怎么研究它，
常微分方程边值问题的数值解法
第8章 常微分方程边值问题的数值解法 引 言 第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.1 一类特殊类型二阶线性常微分方程的边值问题的差分法 设二阶线性常微分方程的边值问题为 (8.2.1)(8.2.2) ()()()(),,(),(), y x q x y x f x a x b y a y b αβ''-=<
常微分方程组(边值)
常微分方程组边值问题解法 打靶法Shooting Method （shooting.m ） %打靶法求常微分方程的边值问题 function [x,a,b,n]=shooting(fun,x0,xn,eps) if nargin<3 eps=1e-3; end x1=x0+rand; [a,b]=ode45(fun,[0,10],[0,x0]'); c0=b(length(b),1); [a,b]=ode45(fun,[0,10],[0,x1]'); c1=b(length(b),1); x2=x1-(c1-xn)*(x1-x0)/(c1-c0); n=1; while (norm(c1-xn)>=eps & norm(x2-x1)>=eps) x0=x1;x1=x2; [a,b]=ode45(fun,[0,10],[0,x0]'); c0=b(length(b),1); [a,b]=ode45(fun,[0,10],[0,x1]'); c1=b(length(b),1) x2=x1-(c1-xn)*(x1-x0)/(c1-c0); n=n+1; end x=x2; 应用打靶法求解下列边值问题： ()()??? ????==- =010004822y y y dx y d 解：将其转化为常微分方程组的初值问题
()????? ? ?????==-==t dx dy y y y dx dy y dx dy x 0011221 048 命令： x0=[0:0.1:10]; y0=32*((cos(5)-1)/sin(5)*sin(x0/2)-cos(x0/2)+1); 真实解 plot(x0,y0,'r') hold on [x,y]=ode45('odebvp',[0,10],[0,2]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,5]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,8]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,10]'); plot(x,y(:,1))
偏微分方程边值问题的数值解法论文
求解偏微分方程的边值问题 本实验学习使用MATLAB 的图形用户命令pdetool 来求解偏微分方程的边值问题。这个工具是用有限元方法来求解的，而且采用三角元。我们用个例题来说明它的用法。 一、MATLAB 支持的偏微分方程类型 考虑平面有界区域D 上的二阶椭圆型PDE 边值问题： ()c u u f α-??+=g (1.1) 其中 (1) , (2) a,f D c x y ?????=? ????? 是上的已知函数(3)是标量或22的函数方阵 未知函数为(,) (,)u x y x y D ∈。它的边界条件分为三类： （1）Direchlet 条件： hu f = (1.2) （2）Neumann 条件： ()n c u qu g ?+=g (1.3) （3）混合边界条件：在边界D ?上部分为Direchlet 条件，另外部分为Neumann 条件。 其中,,,,h r q g c 是定义在边界D ?的已知函数，另外c 也可以是一个2*2的函数矩阵，n 是沿边界的外法线的单位向量。 在使用pdetool 时要向它提供这些已知参数。 二、例题 例题1 用pdetool 求解 22D 1 D: 10u x y u ??-?=+≤??=?? (1.4)
解：首先在MATLAB 的工作命令行中键入pdetool ，按回牟键确定，于是出现PDE Toolbox 窗口，选Genenic Scalar模式． ( l ）画区域圆 单击椭圆工具按钮，大致在（0，0）位置单击鼠标右键，拖拉鼠标到适当位置松开。为了保证所绘制的圆是标准的单位园，在所绘园上双击，打开 Object Dialog 对话框，精确地输入
第十一章 常微分方程边值问题的数值解法汇总
第十一章 常微分方程边值问题的数值解法 工程技术与科学实验中提出的大量问题是常微分方程边值问题．本章将研究常微分方程边值问题的数值求解方法.主要介绍三种边界条件下的定解问题和两大类求解边值问题的数值方法,打靶法算法和有限差分方法. 11.1 引言 在很多实际问题中都会遇到求解常微分方程边值问题. 考虑如下形式的二阶常微分方程 ),,(y y x f y '='', b x a <<, (11.1.1) 在如下三种边界条件下的定解问题: 第一种边界条件: α=)(a y , β=)(b y (11.1.2) 第二种边界条件: α=')(a y , β=')(b y (11.1.2) 第三种边界条件: ? ? ?=-'=-'101 0)()()()(b b y b y a a y a y βα, (11.1.13) 其中0 0, ,00000>+≥≥b a b a . 常微分方程边值问题有很多不同解法, 本书仅介绍打靶方法和有限差分方法. 11.2 打靶法 对于二阶非线性边值问题 ()()().,,βα==≤≤'=''b y a y b x a y y x f y ，，， (11.2.1) 打靶法近似于使用初值求解的情况． 我们需要利用一个如下形式问题初值解的序列： ()()v a w a w b x a w w x f w ='=≤≤'='')(,,，，，α， （11.2.2） 引进参数v 以近似原边界值问题的解．选择参数k v v =，以使： ()()β==∞ →b y v b w k k ,lim ， （11.2.3）
其中),(k v x w 定义为初值问题（11.2.2）在k v v =时的解，同时()x y 定义为边值问题（11.2.1）的解． 首先定义参数0v ，沿着如下初值问题解的曲线，可以求出点),(αa 对应的初始正视图 ()()v a w a w b x a w w x f w ='=≤≤'='')(,,，，，α． （11.2.4） 如果),(0v b w 不严格收敛于β，那么我们选择1v 等值以修正近似值，直到),(0v b w 严格逼近β． 为了取得合适的参数k v ，现在假定边值问题（11.2.1）有唯一解，如果),(v x w 定义为初始问题（11.2.2）的解，那么v 可由下式确定： 0),(=-βv b w ． （11.2.5） 由于这是一个非线性方程，我们可以利用Newton 法求解．首先选择初始值0v ，然后由下式生成序列 ),)(()),((111----- =k k k k v b dv dw v b w v v β，此处),(),)(( 11--=k k v b dv dw v b dv dw ， （11.2.6） 同时要求求得),)(( 1-k v b dv dw ，因为),(v b w 的表达式未知，所以求解这个有一点难度；我们只能得到这么一系列的值。 ，，，),(),(),(),(1210-??k v b w v b w v b w v b w 假如我们如下改写初值问题（11.2.2），使其强调解对x 和v 的依赖性 ()()v v a w v a w b x a v x w v x w x f w ='=≤≤'=''),(,),(),,(,，，，α，（11.2.7） 保留初始记号以显式与x 的微分相关．既然要求当k v v =时),)((v b dv dw 的值，那么我们需要求出表达式（11.2.7）关于v 的偏导数．过程如下： )),(),,(,(),(v x w v x w x v f v x v w '??=?''? ),()),(),,(,()),(),,(,(v x v w v x w v x w x w f v x v x w v x w x x f ??'??+??'??= ) ,()),(),,(,(v x v w v x w v x w x w f ?'?''??+ 又因为x 跟v 相互独立，所以当b x a ≤≤上式如下；
常微分方程数值解法的误差分析教材
淮北师范大学 2013届学士学位论文 常微分方程数值解法的误差分析 学院、专业数学科学学院数学与应用数学 研究方向计算数学 学生姓名李娜 学号 20091101070 指导教师姓名陈昊 指导教师职称讲师 年月日
常微分方程数值解法的误差分析 李娜 （淮北师范大学数学科学学院，淮北，235000） 摘要 自然界与工程技术中的很多现象，往往归结为常微分方程定解问题。许多偏微分方程问题也可以化为常微分方程问题来近似求解。因此，研究常微分方程的数值解法是有实际应用意义的。数值解法是一种离散化的数学方法，可以求出函数的精确解在自变量一系列离散点处的近似值。随着计算机计算能力的增强以及数值计算方法的发展，常微分方程的数值求解方法越来越多，比较成熟的有Euler 法、后退Euler法、梯形方法、Runge—Kutta方法、投影法和多步法，等等．本文将对这些解的误差进行分析,以求能够得到求解常微分数值解的精度更好的方法。 关键词：常微分方程, 数值解法, 单步法, 线性多步法, 局部截断误差
Error Analysis of Numerical Method for Solving the Ordinary Differential Equation Li Na (School of Mathematical Science, Huaibei Normal University, Huaibei, 235000) Abstract In nature and engineering have many phenomena , definite solution of the problem often boils down to ordinary differential equations. So study the numerical solution of ordinary differential equations is practical significance. The numerical method is a discrete mathematical methods, and exact solution of the function can be obtained in the approximation of a series of discrete points of the argument.With the enhanced computing power and the development of numerical methods,ordinary differential equations have more and more numerical solution,there are some mature methods. Such as Euler method, backward Euler method, trapezoidal method, Runge-Kutta method, projection method and multi-step method and so on.Therefore, numerical solution of differential equation is of great practical significance. Through this paper, error of these solutions will be analyzed in order to get a the accuracy better way to solve the numerical solution of ordinary differential. Keywords:Ordinary differential equations, numerical solution methods, s ingle ste p methods, l inear multi-step methods, local truncation error
常微分方程数值解法
i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。
常微分方程组(边值)
常微分方程组边值问题解法 打靶法Shooti ng Method (shoot in g.m ) % 丁靶法求常微分方程的边值问题 function [x,a,b ,n]=shooti ng(fu n, xO,x n, eps) if nargin<3 eps=1e-3; end x1=x0+ra nd; [a,b]=ode45(fu n, [0,10],[0,x0]'); c0=b(le ngth(b),1); [a,b]=ode45(fu n, [0,10],[0,x1]'); c1=b(le ngth(b),1); x2=x1-(c1-x n)*(x1-x0)/(c1-c0); n=1; while (no rm(c1-x n)>=eps & no rm(x2-x1)>=eps) x0=x1;x 仁x2; [a,b]=ode45(fu n,[ 0,10],[0,x0]'); cO=b(le ngth(b),1); [a,b]=ode45(fu n,[ 0,10],[0,x1]'); c1= b(le ngth(b),1) x2=x1-(c1-x n)*(x1-x0)/(c1-c0); n=n+1; end x=x2; 应用打靶法求解下列边值问题: y 10 0 解：将其转化为常微分方程组的初值问题
命令: xO=[O:O.1:1O]; y0=32*((cos(5)-1)/si n( 5)*si n(x0/2)-cos(x0/2)+1); plot(xO,yO,'r') hold on [x,y]=ode45('odebvp',[0,10],[0,2]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,5]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,8]'); plot(x,y(:,1)) [x,y]=ode45('odebvp',[0,10],[0,10]'); plot(x,y(:,1)) dy i dx y 2 dy 2 dx y i 0 y 4 y o dy dx X0 真实解 30 ' 12^4567^9 10
常微分方程初值问题数值解法
常微分方程初值问题数值解法 朱欲辉 （浙江海洋学院数理信息学院, 浙江舟山316004） [摘要]：在常微分方程的课程中讨论的都是对一些典型方程求解析解的方法.然而在生产实 际和科学研究中所遇到的问题往往很复杂, 在很多情况下都不可能给出解的解析表达式. 本篇文章详细介绍了常微分方程初值问题的一些数值方法, 导出了若干种数值方法, 如Euler法、改进的Euler法、Runge－Kutta法以及线性多步法中的Adams显隐式公式和预测校正 公式, 并且对其稳定性及收敛性作了理论分析. 最后给出了数值例子, 分别用不同的方法计算出近似解, 从得出的结果对比各种方法的优缺点. [关键词]：常微分方程;初值问题; 数值方法; 收敛性; 稳定性; 误差估计 Numerical Method for Initial-Value Problems Zhu Yuhui (School of Mathematics, Physics, and Information Science, Zhejiang Ocean University, Zhoushan, Zhejiang 316004) [Abstract]:In the course about ordinary differential equations, the methods for analytic solutions of some typical equations are often discussed. However, in scientific research, the problems are very complex and the analytic solutions about these problems can’t be e xpressed explicitly. In this paper, some numerical methods for the initial-value problems are introduced. these methods include Euler method, improved Euler method, Runge-Kutta method and some linear multistep method (e.g. Adams formula and predicted-corrected formula). The stability and convergence about the methods are presented. Some numerical examples are give to demonstrate the effectiveness and accuracy of theoretical analysis. [Keywords]:Ordinary differential equation; Initial-value problem; Numerical method; Convergence; Stability;Error estimate
两点边值问题的两种数值解法
常微分方程组两点边值问题的数值解法 ----张亚苗2011年9月 3)1(1)0(04===-''y y y y 可化为微分方程组3 )1(1)0(41221==='='y y y y y y 方法一：配置法 Matlab 程序： function bvcollation clc solinit = bvpinit(linspace(0,1,20),[100 600]);% sol = bvp4c(@twoode,@twobc,solinit); x = linspace(0,1,20); y = deval(sol,x); y' plot(x,y(1,:),x,y(2,:)); end %微分方程组 function dydx = twoode(x,y) dydx = [ y(2) 4*y(1)]; end %边值条件 function res = twobc(ya,yb) res = [ ya(1)-1 yb(1)-3]; end 运行结果： 1.0000 -0.4203 0.9834 -0.2117 0.9777 -0.0055 0.9828 0.2007 0.9988 0.4091 1.0259 0.6220 1.0644 0.8419 1.1147 1.0710 1.1774 1.3121 1.2531 1.5677 1.3427 1.8407 1.4472 2.1341 1.5678 2.4512 1.7057 2.7954 1.8626 3.1707 2.0401 3.5811 2.2402 4.0313 2.4652 4.5261 2.7175 5.0712 3.0000 5.6724
常微分方程数值解法
第八章 常微分方程数值解法 考核知识点： 欧拉法，改进欧拉法，龙格-库塔法，单步法的收敛性与稳定性。 考核要求： 1． 解欧拉法，改进欧拉法的基本思想；熟练掌握用欧拉法，改进欧拉法、求微 分方程近似解的方法。 2． 了解龙格-库塔法的基本思想；掌握用龙格-库塔法求微分方程近似解的方 法。 3． 了解单步法的收敛性、稳定性与绝对稳定性。 例1 用欧拉法，预估——校正法求一阶微分方程初值问题 ? ??=-='1)0(y y x y ，在0=x （0，1）0.2近似解 解 （1）用1.0=h 欧拉法计算公式 n n n n n n x y y x y y 1.09.0)(1.01+=-+=+，1.0=n 计算得 9.01=y 82.01.01.09.09.02=?+?=y （2）用预估——校正法计算公式 1,0)(05.01.09.0)0(111)0(1=???-+-+=+=++++n y x y x y y x y y n n n n n n n n n 计算得 91.01=y ，83805.02=y 例2 已知一阶初值问题 ???=-='1 )0(5y y y 求使欧拉法绝对稳定的步长n 值。 解 由欧拉法公式 n n n n y h y h y y )51(51-=-=+ n n y h y ~)51(~1-=+
相减得01)51()51(e h e h e n n n -==-=-Λ 当 151≤-h 时，4.00≤ 二阶常微分方程边值问题的数值解法 摘要 求解微分方程数值解的方法是多种多样的，它本身已形成一个独立的研究方向，其要点是对微分方程定解问题进行离散化．本文以研究二阶常微分方程边值问题的数值解法为目标，综合所学相关知识和二阶常微分方程的相关理论，通过对此类方程的数值解法的研究，系统的复习并进一步加深对二阶常微分方成的数值解法的理解，为下一步更加深入的学习和研究奠定基础． 对于二阶常微分方程的边值问题，我们总结了两种常用的数值方法：打靶法和有限差分法．在本文中我们主要探讨关于有限差分法的数值解法．构造差分格式主要有两种途径：基于数值积分的构造方法和基于Taylor展开的构造方法．后一种更为灵活，它在构造差分格式的同时还可以得到关于截断误差的估计．在本文中对差分方法列出了详细的计算步骤和Matlab 程序代码，通过具体的算例对这种方法的优缺点进行了细致的比较．在第一章中，本文将系统地介绍二阶常微分方程和差分法的一些背景材料．在第二章中，本文将通过Taylor展开分别求得二阶常微分方程边值问题数值解的差分格式．在第三章中，在第二章的基础上利用Matlab求解具体算例，并进行误差分析． 关键词：常微分方程,边值问题，差分法，Taylor展开，数值解 The Numerical Solutions of Second-Order Ordinary Differential Equations with the Boundary Value Problems ABSTRACT The numerical solutions for solving differential equations are various. It formed an independent research branch. The key point is the discretization of the definite solution problems of differential equations. The goal of this paper is the numerical methods for solving second-order ordinary differential equations with the boundary value problems. This paper introduces the mathematics knowledge with the theory of finite difference. Through solving the problems, reviewing what have been learned systematically and understanding the ideas and methods of the finite difference method in a deeper layer, we can establish a foundation for the future learning. 实验4 常微分方程的数值解 【实验目的】 1．掌握用MATLAB软件求微分方程初值问题数值解的方法； 2．通过实例用微分方程模型解决简化的实际问题； 3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。 【实验内容】 题3 小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。 模型及其求解 火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。 在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式： a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8 dh/dt=v 又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下： function [ dx ] = rocket( t,x ) a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8; dx=[x(2);a]; end ts=0:1:60; x0=[0,0]; [t,x]=ode45(@rocket,ts,x0); h=x(:,1); v=x(:,2); a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8; [t,h,v,a]; 数据如下： t h v a 0 0 0 13.06 1.00 6.57 13.19 13.30 2.00 26.44 26.58 1 3.45 3.00 59.76 40.06 13.50 4.00 106.57 53.54 13.43 5.00 16 6.79 66.89 13.26 6.00 240.27 80.02 12.99 7.00 326.72 92.83 12.61 8.00 425.79 105.22 12.15 9.00 536.99 117.11 11.62 10.00 659.80 128.43 11.02 11.00 793.63 139.14 10.38 12.00 937.85 149.18 9.71 13.00 1091.79 158.55 9.02 14.00 1254.71 167.23 8.33 15.00 1425.93 175.22 7.65 16.00 1604.83 182.55 6.99 17.00 1790.78 189.22 6.36 18.00 1983.13 195.27 5.76 19.00 2181.24 200.75 5.21 20.00 2384.47 205.70 4.69 21.00 2592.36 210.18 4.22 22.00 2804.52 214.19 3.79 23.00 3020.56 217.79 3.41 24.00 3240.08 221.01 3.07 25.00 3462.65 223.92 2.77 26.00 3687.88 226.56 2.50 27.00 3915.58 228.97 2.27 Matlab 求解常微分方程边值问题的方法：bvp4c 函数 常微分方程的边值问题，即boundary value problems ，简称BVP 问题，是指表达形式为 (,)((),())0'=??=?y f x y g y a y b 或(,,)((),(),)0'=??=? y f x y p g y a y b p 的方程组（p 是未知参数），在MA TLAB 中使用积分器bvp4c 来求解。 [命令函数] bvp4c [调用格式] sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) sol 为一结构体，sol.x 、sol.y 、sol.yp 分别是所选择的网格点及其对应的y(x)与y'(x)数值; bvp4c 为带边值条件常微分方程积分器的函数命令；odefun 为描述微分方程组的函数文件；bcfun 为计算边界条件g(f(a),f(b),p)=0的函数文件；solinit 为一结构体，solinit.x 与solinit.y 分别是初始网格的有序节点与初始估计值，边界值条件分别对应a=solinit.x(l)和b=solinit.x(end)； options 为bvpset 命令设定的可选函数，可采用系统默认值；p1, p2…为未知参数。 例 求常微分方程0''+=y y 在(0)2=y 与(4)2=-y 时的数值解。 [解题过程] 仍使用常用方法改变方程的形式： 令1=y y ，21'=y y ，则原方程等价于标准形式的方程组1221 ?'=??'=-??y y y y ； 将其写为函数文件twoode.m ； 同时写出边界条件函数对应文件twobc.m ； 分别使用结构solinit 和命令bvp4c 确定y-x 的关系； 作出y-x 的关系曲线图。 [算例代码] solinit =bvpinit(linspace(0,4,5),[1 0]); % linspace(0,4,5)为初始网格，[1,0]为初始估计值 sol=bvp4c(@twoode,@twobc,solinit); % twoode 与twobc 分别为微分方程与边界条件的函数，solinit 为结构 x=linspace(0,4); %确定x 范围 y=deval(sol,x); %确定y 范围 plot(x,y(1,:)); %画出y-x 的图形 %定义twoode 函数（下述代码另存为工作目录下的twoode.m 文件） function dydx= twoode(x,y) %微分方程函数的定义 dydx =[y(2) -abs(y(1))]; %定义twobc 函数（下述代码另存为工作目录下的twobc.m 文件） function res= twobc(ya,yb); %边界条件函数的定义 res=[ya(1);yb(1)+2]; -293- 第十五章 常微分方程的解法 建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如 22x y dx dy +=。于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。 §1 常微分方程的离散化 下面主要讨论一阶常微分方程的初值问题，其一般形式是 ?????=≤≤=0 )() ,(y a y b x a y x f dx dy （1） 在下面的讨论中，我们总假定函数),(y x f 连续，且关于y 满足李普希兹(Lipschitz)条 件，即存在常数L ，使得 |||),(),(|y y L y x f y x f ?≤? 这样，由常微分方程理论知，初值问题(1)的解必定存在唯一。 所谓数值解法，就是求问题（1）的解)(x y 在若干点 b x x x x a N =<<<<=L 210 处的近似值),,2,1(N n y n L =的方法，),,2,1(N n y n L =称为问题（1）的数值解， n n n x x h ?=+1称为由n x 到1+n x 的步长。今后如无特别说明，我们总取步长为常量h 。 建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： （i ）用差商近似导数 若用向前差商h x y x y n n ) ()(1?+代替)('n x y 代入（1）中的微分方程，则得 )1,,1,0())(,() ()(1?=≈?+N n x y x f h x y x y n n n n L 化简得 ))(,()()(1n n n n x y x hf x y x y +≈+ 如果用)(n x y 的近似值n y 代入上式右端，所得结果作为)(1+n x y 的近似值，记为1+n y ， 则有 )1,,1,0() ,(1?=+=+N n y x hf y y n n n n L （2） 这样，问题（1）的近似解可通过求解下述问题 ?? ?=?=+=+) () 1,,1,0(),(01a y y N n y x hf y y n n n n L （3） 得到，按式（3）由初值0y 可逐次算出N y y y ,,,21L 。式（3）是个离散化的问题，称为差分方程初值问题。 i.常微分方程初值问题数值解法 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法--差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<<<<=L （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1,,1n n n n u u hf t u n N +=+=-L 方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t L 上的差分解1,,N u u L 。(完整版)二阶常微分方程边值问题的数值解法毕业论文
常微分方程的数值解
Matlab求解常微分方程边值问题的方法
15第十五章 常微分方程的解法
(整理)常微分方程数值解法