上一页 下一页湘潭大学数学与计算科学学院 1
第八章 常微分方程数值解
§ 1 引 论上一页 下一页湘潭大学数学与计算科学学院 2
解析法求解常微分方程的初值问题
( 0 ) 1
yy
y
如又由 得 0 1,1C e C(0 ) 1y?
xye?初值问题解为很多时候解析解求不出来,如 22(,)f x y x y
xy ce?由 得yy
00
(,),,| |,
( ),,
y f x y a x b y
y a y y R
(1.1)
上一页 下一页湘潭大学数学与计算科学学院 3
常微分方程的初值问题
00
(,),,| |,
( ),,
y f x y a x b y
y a y y R
(1.1)
其中 f 为,xy 的已知函数,0y 为给定的初值,
[,]a b R?,G
[,],G a b R
为简便起见，我们将区域,记为 即设,f G R?为连续映射,若存在常数 L>0使得不等式
1 2 1 2| (,) (,) | | |,f x y f x y L y y
上一页 下一页湘潭大学数学与计算科学学院 4
设,f G R?为连续映射,若存在常数 L>0使得不等式
1 2 1 2| (,) (,) | | |,f x y f x y L y y
对一切 12(,),(,)x y x y G?都成立，
则称 f (x,y)在 G上 关于 y满足 Lipschitz条件,
而式中的常数 L称为 Lipschitz常数,
一切在 G上关于 y满足 Lipschitz条件的连续映射 f
所构成的集合记为?，
而相应的初值问题（ 1.1）构成的问题类记为?.
上一页 下一页湘潭大学数学与计算科学学院 5
定理 1?中的任何初值问题在 [a,b]上有连续可微的解存在并且惟一,
定义 1 初值问题 (1.1)称为在 [a,b]上是 适定的,
0,0,k如果存在常数 0,使得对于任何的正数
(,)f x y 0,y及任给的函数 和常数 当
00| |,yy| (,) (,) |,(,)f x y f x y x y G
时初值问题
0
(,),[,]
( ),
z f x z x a b
z a y
有解 ()zx 存在，且不等式 | ( ) ( ) |y x z x k 对任给
[,]x a b? 都成立,
上一页 下一页湘潭大学数学与计算科学学院 6
定理 2?中的任何初值问题在 [a,b]上是适定的,
上一页 下一页湘潭大学数学与计算科学学院 7
§ 2 Euler方法一,Euler方法二、误差分析三,Euler方法的 收敛性和稳定性上一页 下一页湘潭大学数学与计算科学学院 8
1,Euler方法
0 1 2,Na x x x x b
记：
1ii
bah x x
N?
因为：
(等距剖分 )
( ) ( ) (,( ) )xhxy x h y x f y d(积分方程 )
令：,mxx? 有：
上一页 下一页湘潭大学数学与计算科学学院 9
( ) ( ) (,( ) ),m m m m my x h y x h f x y x R
1 (,( ) ) (,( ) )m
m
x
m m mxR f x y x dx h f x y x
1( ) ( ) (,( ) ),m m m my x y x h f x y x
1 (,),0,1,,1m m m my y h f x y m N
—— Euler方法截去 mR 有：
由于,00()y x y?(已知 ),
又称 Euler折线法,
可得递推关系：
上一页 下一页湘潭大学数学与计算科学学院 10
欧拉方法的几何意义：
0 1 2
n
x x x X
x
0y
Xy
ny
1xy
2xy
1y
2y
h步长
Euler方法的几何意义上一页 下一页湘潭大学数学与计算科学学院 11
2,误差分析称为局部截断误差，mR
计算时 的误差,()my x h?
记：
( ),m m my x y
1 2 1 2| (,) (,) | | |,f x y f x y K x x
有：
.m?估计关于假设 (,)f x y x 满足 Lipschitz条件,
精确值时，
()mmy y x? 为它表示当上一页 下一页湘潭大学数学与计算科学学院 12
1 | (,( ) ) (,( ) ) |m
m
x
mx f x y x f x y x dx
1 | (,( ) ) (,( ) ) |m
m
x
m m mx f x y x f x y x d x
11| | | ( ) ( ) |mmxx
mmK x x d x L y x y x d x
12 / 2 | ( ( ) ) | ( )m
m
x
m m mxKh L y x x x x x d x?
2( ) / 2K L M h
其中：
[,] [,]0 1,m a x | ( ) | m a x | (,( ) ) |,x a b x a bM y x f x y x
1| | | [ (,( ) ) (,( ) ) ] |m
m
x
m m mxR f x y x f x y x d x
上一页 下一页湘潭大学数学与计算科学学院 13
几何分析：
Euler公式的误差
ix 1ix?
iy
1iyx? 1iy?
y y x?
( ),i i i iY y x x f x y
||mRR?
记 2( ) / 2R K LM h，则有上一页 下一页湘潭大学数学与计算科学学院 14
整体截断误差,
( ),m m my x y
由：
( ) ( ) (,( ) ),m m m m my x h y x h f x y x R
1 (,),m m m my y hf x y
1 [ (,( ) ) (,) ],m m m m m m mh f x y x f x y R
上一页 下一页湘潭大学数学与计算科学学院 15
从而有：
1| | | | | |,m m mh L R
对任一 0,1,,1,mN
1| | ( 1 ) | |mm h L R
2 2( 1 ) | | ( 1 )mh L h L R R
1
0
0
( 1 ) | | ( 1 )
m
mj
j
h L R h L?
有：
上一页 下一页湘潭大学数学与计算科学学院 16
0( 1 ) | | [ ( 1 ) 1 ]
mm Rh L h L
hL
( ) ( )
0| | ( 1 )
L b a L b aRee
hL?
于是便得 Euler方法的 整体截断误差界
( ) ( )
0| | | | ( ) ( 1 ),2
L b a L b a
m
hKe M e
L
(*)
()1 12 L b am hM eL
0 0 0 0y y x若,
1 ( ),,mM K L M O h这 里 =可 记 为说明 Euler方法的整体截断误差与 h同阶。
上一页 下一页湘潭大学数学与计算科学学院 17
因 此 难 以 使 用 上 述 结 果 。
3 y x M M） 实 际 上,的 估 计 式 中 的 是 不 知 道 的 。
2,） 方 法 的 精 度 由 整 体 误 差 决 定 但 局 部 误 差 容 易 估 计 。
注意,对于 Euler方法
1 ） 一 般 来 说,整 体 误 差 比 局 部 误 差 低 一 阶 。
上一页 下一页湘潭大学数学与计算科学学院 18
定理 5 设 f (x,y)属于?且关于 x满足 Lipschitz条件，
其 Lipschitz常数为 K,且当 时，0h?
0 ( ),y y a?
则 Euler方法一致收敛于真解 ( ),myx
成立,
并且 有估计式
( ) ( )
0| | | | ( ) ( 1 )2
L b a L b a
m
hKe M e
L
(*)
上一页 下一页湘潭大学数学与计算科学学院 19
2.2 改进的 Euler方法
00
(,),,| |,
( ),,
y f x y a x b y
y a y y R
等价于 积分方程：
( ) ( ) (,( ) ),xhxy x h y x f y d
微分方程：
上一页 下一页湘潭大学数学与计算科学学院 20
( ) ( ) (,( ) )xhxy x h y x f y d
有,(令 )mxx?
( ) ( )mmy x h y x (,( ) )m
m
xh
x f y d
11[ (,) (,) ] / 2m m m m m my h f x y f x y R
1 (,( ) ) [ (,( ) )
2
m
m
x
m m mx
hR f x y x dx f x y x
11(,( ) ) ],mmf x y x
去掉,mR
1( ) ( )mmy x y x 11[ (,) (,) ] / 2.m m m mh f x y f x y
有：
上一页 下一页湘潭大学数学与计算科学学院 21
设 为 的近似值,my ()myx
1mmyy 11[ (,) (,) ] / 2,m m m mh f x y f x y
称为 改进的 Euler方法,
称为改进的 Euler方法的 局部截断误差,mR
误差分析,仍记 ( ),m m my x y
注意：
1
3
1( ) [ ( ) ( ) ] ( ),2 1 2
m
m
x
m m mx
hhy x d x y x y x y x h
0 1,
则：
上一页 下一页湘潭大学数学与计算科学学院 22
于是：
若记
[,]m a x | ( ) |,x a bM y x
3 ( ) / 1 2,mmR h y x h
3| | / 1 2,mR h M?
整体截断误差的阶由局部截断误差的阶来决定,
可见改进的 Euler方法误差比 Euler方法要 高一阶,
则有上一页 下一页湘潭大学数学与计算科学学院 23
三,Euler方法的 收敛性和稳定性
立是差分方程的解，若成初值问题的解是设
my
O D Exy,
0 1li m m a x 0mmh mn y y x
的。也称此数值方法是收敛则称差分方程是收敛的不能无限缩小。以实际上对扩大了对解的影响。所误差的积累可能时离散点增多。舍入当
h
h 0?
hOE u l er 方法是收敛的，且阶为结论,
注意,
收敛性上一页 下一页湘潭大学数学与计算科学学院 24
的变化情况。
导致公式的差分方程解稳定性描述初值的改变稳定性
,2,1,
,0
,
,,,,
00
0
000
mvucvu
abmhhh
vu
vuhc
mm
mm
时当满足：的解它们使对任给的初值定义方法是稳定的。则称 E u l er
上一页 下一页湘潭大学数学与计算科学学院 25
结论,
方法是稳定的。则条件，满足若
E u l er
L i p s ch i zyxf,
1,m m m mu u h f x u
mmm vu证明：记
mhL?1
1,m m m mv v h f x v
mmmmmm vxfuxfh,, 1
() 0,L b ae
011 mhL
上一页 下一页湘潭大学数学与计算科学学院 26
计算问题：
1mmyy 11[ (,) (,) ] / 2,m m m mh f x y f x y
隐式计算格式由迭代法去完成,
1 1 1(,) (,) 022m m m m m m
hhy f x y y f x y
将上式变形为记
1 1 1 1( ) (,) (,)22m m m m m m m
hhy y f x y y f x y?
求 1my? 即求隐式方程
1( ) 0my的根。
上一页 下一页湘潭大学数学与计算科学学院 27
总结通过对 Euler方法的讨论可以看到,微分方程数值方法的研究应包括以下方面
1.数值计算公式的构造 ;
2.方法稳定性,收敛性的研究 ;
3.方法的误差估计 ;
4.方法的实现等,
上一页 下一页湘潭大学数学与计算科学学院 28
欧拉公式的改进：
隐式欧拉法向后差商近似导数
h
xyxyxy )()()( 01
1
))(,()( 1101 xyxfhyxy
)1,...,0(),( 111 niyxfhyy iiii
由于未知数 yi+1同时出现在等式的两边，不能直接得到，故称为 隐式 欧拉公式，而前者称为 显式 欧拉公式。
一般先用显式计算一个初值，再 迭代 求解。
隐式 欧拉法的局部截断误差：
11 )( iii yxyR )()( 322 hOxy ih
即隐式欧拉公式具有 1 阶精度。
上一页 下一页湘潭大学数学与计算科学学院 29
梯形公式 — 显、隐式两种算法的 平均
)1,.,,,0()],(),([2 111 niyxfyxfhyy iiiiii
注,的确有局部截断误差,
即梯形公式具有 2 阶精度，比欧拉方法有了进步。
但注意到该公式是 隐式 公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。
)()( 311 hOyxyR iii
中点欧拉公式中心差商近似导数
h
xyxyxy
2
)()()( 02
1
x0 x2x1
))(,(2)()( 1102 xyxfhxyxy
1,...,1),(211 niyxfhyy iiii
假设,则可以导出即中点公式具有 2 阶精度。
)(),( 11 iiii xyyxyy )()( 311 hOyxyR iii
需要 2个初值 y0和 y1来启动递推过程，这样的算法称为 双步法,
而前面的三种算法都是 单步法 。
上一页 下一页湘潭大学数学与计算科学学院 30
方 法
显式欧拉隐式欧拉梯形公式中点公式简单 精度低稳定性最好 精度低,计算量大精度提高 计算量大精度提高,显式 多一个初值,可能影响精度上一页 下一页湘潭大学数学与计算科学学院 31
改进欧拉法
Step 1,先用 显式 欧拉公式作 预测,算出 ),(1 iiii yxfhyy
Step 2,再将 代入 隐式 梯形公式的右边作 校正,得到1?iy
)],(),([2 111 iiiiii yxfyxfhyy
注,此法亦称为 预测 -校正法 。可以证明该算法具有 2 阶精度，同时可以看到它是个 单步 递推格式，比隐式公式的迭代求解过程 简单 。后面将看到，它的 稳定性高 于显式欧拉法。
)1,.,,,0(),(,),(2 11 niyxfhyxfyxfhyy iiiiiiii
上一页 下一页湘潭大学数学与计算科学学院 32
§ 3 龙格 - 库塔法建立高精度的单步递推格式。
单步递推法的 基本思想 是从 ( xi,yi ) 点出发，以 某一斜率 沿直线达到 ( xi+1,yi+1 ) 点。欧拉法及其各种变形所能达到的最高精度为 2阶 。
考察改进的欧拉法，可以将其改写为：
),(
),(
2
1
2
1
12
1
211
hKyhxfK
yxfK
KKhyy
ii
ii
ii
斜率一定取 K1 K2
的 平均值 吗？
步长一定是一个 h 吗？
上一页 下一页湘潭大学数学与计算科学学院 33
首先希望能确定系数?1,?2,p，使得到的算法格式有 2阶精度，即在 的前提假设下，使得)(
ii xyy?
)()( 311 hOyxyR iii
Step 1,将 K2 在 ( xi,yi ) 点作 Taylor 展开
)(),(),(),(
),(
2
1
12
hOyxfp h Kyxp h fyxf
p h KyphxfK
iiyiixii
ii
)()()( 2hOxyphxy ii
将改进欧拉法推广为：
),(
),(
][
12
1
22111
phKyphxfK
yxfK
KKhyy
ii
ii
ii
),(),(),(
),(),(
),()(
yxfyxfyxf
dx
dy
yxfyxf
yxf
dx
d
xy
yx
yx
Step 2,将 K2 代入第 1式，得到
)()()()(
)]()()([)(
32
221
2
211
hOxyphxyhy
hOxyphxyxyhyy
iii
iiiii
上一页 下一页湘潭大学数学与计算科学学院 34
Step 3,将 yi+1 与 y( xi+1 ) 在 xi 点的 泰勒 展开作比较
)()()()( 322211 hOxyphxyhyy iiii
)()(2)()()( 3
2
1 hOxy
hxyhxyxy
iiii
要求,则必须有,)()( 3
11 hOyxyR iii
2
1,1
221 p
这里有 个未知数,个方程。
3
2
存在 无穷多个解 。所有满足上式的格式统称为 2阶龙格 - 库塔格式 。
2
1,1
21p
注意到,就是改进的欧拉法。
Q,为获得更高的精度，应该如何进一步推广？
上一页 下一页湘潭大学数学与计算科学学院 35
其中?i ( i = 1,…,m ),?i ( i
= 2,…,m ) 和?ij ( i = 2,…,
m; j = 1,…,i?1 ) 均为待定系数，确定这些系数的步骤与前面相似。
)...,(
......
),(
),(
),(
]...[
112211
23213133
12122
1
22111
mm mmmmim
ii
ii
ii
mmii
hKhKhKyhxfK
hKhKyhxfK
hKyhxfK
yxfK
KKKhyy
最常用为四级 4阶 经典龙格 -库塔法,
),(
),(
),(
),(
)22(
34
2223
1222
1
432161
hKyhxfK
KyxfK
KyxfK
yxfK
KKKKyy
ii
h
i
h
i
h
i
h
i
ii
h
ii
上一页 下一页湘潭大学数学与计算科学学院 36
注：
龙格 -库塔法 的主要运算在于计算 Ki 的值，即计算 f 的值。 Butcher 于 1965年给出了计算量与可达到的最高精度阶数的关系：
753
可达到的最高精度
642每步须算 Ki 的个数
)( 2hO )( 3hO )( 4hO )( 5hO )( 6hO)( 4hO )( 2?nhO
8?n
由于龙格 -库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用 低阶算法 而将步长 h取小 。
上一页 下一页湘潭大学数学与计算科学学院 37
§ 3 收敛性与稳定性
收敛性定义 若某算法对于任意固定的 x = xi = x0 + i h,当 h?0
( 同时 i) 时有 yi? y( xi ),则称该算法是 收敛 的 。
例,就初值问题 考察欧拉显式格式的收敛性。
0)0( yy
yy?
解,该问题的精确解为 xeyxy?
0)(?
欧拉公式为
iiii yhyhyy )1(1 0)1( yhy ii
对任意固定的 x = xi = i h,有
ii xhhxi hyhyy ])1[()1( /10/0
eh hh /10 )1(l i m
)(0 ix xyey i
上一页 下一页湘潭大学数学与计算科学学院 38
定义 若某算法在计算过程中任一步产生的误差在以后的计算中都 逐步衰减,则称该算法是 绝对稳定的 。
一般分析时为简单起见，只考虑 试验方程
yy
常数，可以是复数当步长取为 h 时，将某算法应用于上式，并假设只在初值产生误差,则若此误差以后逐步衰减，就称该算法相对于 绝对稳定,的全体构成 绝对稳定区域 。
我们称 算法 A 比算法 B 稳定,就是指 A 的绝对稳定区域比 B
的 大 。
000 yy
h? h? h
上一页 下一页湘潭大学数学与计算科学学院 39
例,考察显式欧拉法
011 )1( yhyhyy iiii
000 yy 011 )1( yhy ii
01111 )1( iiii hyy
由此可见，要保证初始误差?0 以后逐步衰减，
必须满足：
hh
1|1| h
0-1-2 Re
Img
例,考察隐式欧拉法
11 iii yhyy?
ii yhy
1
1
1 0
1
1 1
1
i
i h
可见绝对稳定区域为,1|1| h
210 Re
Img
注,一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。
上一页 下一页湘潭大学数学与计算科学学院 40
例,隐式龙格 -库塔法
),.,,,1(
).,,,(
].,,[
11
111
mj
hKhKyhxfK
KKhyy
mmjjijij
mmii
而 显式 1~ 4 阶方法的绝对稳定区域为
)2,2( 11
11
KhyhxfK
hKyy
ii
ii其中 2阶方法 的绝对稳定区域为
0 Re
Img
k=1
k=2
k=3
k=4
-1-2-3
-
-
-
1
2
3
Re
Img
无条件稳定
第八章 常微分方程数值解
§ 1 引 论上一页 下一页湘潭大学数学与计算科学学院 2
解析法求解常微分方程的初值问题
( 0 ) 1
yy
y
如又由 得 0 1,1C e C(0 ) 1y?
xye?初值问题解为很多时候解析解求不出来,如 22(,)f x y x y
xy ce?由 得yy
00
(,),,| |,
( ),,
y f x y a x b y
y a y y R
(1.1)
上一页 下一页湘潭大学数学与计算科学学院 3
常微分方程的初值问题
00
(,),,| |,
( ),,
y f x y a x b y
y a y y R
(1.1)
其中 f 为,xy 的已知函数,0y 为给定的初值,
[,]a b R?,G
[,],G a b R
为简便起见，我们将区域,记为 即设,f G R?为连续映射,若存在常数 L>0使得不等式
1 2 1 2| (,) (,) | | |,f x y f x y L y y
上一页 下一页湘潭大学数学与计算科学学院 4
设,f G R?为连续映射,若存在常数 L>0使得不等式
1 2 1 2| (,) (,) | | |,f x y f x y L y y
对一切 12(,),(,)x y x y G?都成立，
则称 f (x,y)在 G上 关于 y满足 Lipschitz条件,
而式中的常数 L称为 Lipschitz常数,
一切在 G上关于 y满足 Lipschitz条件的连续映射 f
所构成的集合记为?，
而相应的初值问题（ 1.1）构成的问题类记为?.
上一页 下一页湘潭大学数学与计算科学学院 5
定理 1?中的任何初值问题在 [a,b]上有连续可微的解存在并且惟一,
定义 1 初值问题 (1.1)称为在 [a,b]上是 适定的,
0,0,k如果存在常数 0,使得对于任何的正数
(,)f x y 0,y及任给的函数 和常数 当
00| |,yy| (,) (,) |,(,)f x y f x y x y G
时初值问题
0
(,),[,]
( ),
z f x z x a b
z a y
有解 ()zx 存在，且不等式 | ( ) ( ) |y x z x k 对任给
[,]x a b? 都成立,
上一页 下一页湘潭大学数学与计算科学学院 6
定理 2?中的任何初值问题在 [a,b]上是适定的,
上一页 下一页湘潭大学数学与计算科学学院 7
§ 2 Euler方法一,Euler方法二、误差分析三,Euler方法的 收敛性和稳定性上一页 下一页湘潭大学数学与计算科学学院 8
1,Euler方法
0 1 2,Na x x x x b
记：
1ii
bah x x
N?
因为：
(等距剖分 )
( ) ( ) (,( ) )xhxy x h y x f y d(积分方程 )
令：,mxx? 有：
上一页 下一页湘潭大学数学与计算科学学院 9
( ) ( ) (,( ) ),m m m m my x h y x h f x y x R
1 (,( ) ) (,( ) )m
m
x
m m mxR f x y x dx h f x y x
1( ) ( ) (,( ) ),m m m my x y x h f x y x
1 (,),0,1,,1m m m my y h f x y m N
—— Euler方法截去 mR 有：
由于,00()y x y?(已知 ),
又称 Euler折线法,
可得递推关系：
上一页 下一页湘潭大学数学与计算科学学院 10
欧拉方法的几何意义：
0 1 2
n
x x x X
x
0y
Xy
ny
1xy
2xy
1y
2y
h步长
Euler方法的几何意义上一页 下一页湘潭大学数学与计算科学学院 11
2,误差分析称为局部截断误差，mR
计算时 的误差,()my x h?
记：
( ),m m my x y
1 2 1 2| (,) (,) | | |,f x y f x y K x x
有：
.m?估计关于假设 (,)f x y x 满足 Lipschitz条件,
精确值时，
()mmy y x? 为它表示当上一页 下一页湘潭大学数学与计算科学学院 12
1 | (,( ) ) (,( ) ) |m
m
x
mx f x y x f x y x dx
1 | (,( ) ) (,( ) ) |m
m
x
m m mx f x y x f x y x d x
11| | | ( ) ( ) |mmxx
mmK x x d x L y x y x d x
12 / 2 | ( ( ) ) | ( )m
m
x
m m mxKh L y x x x x x d x?
2( ) / 2K L M h
其中：
[,] [,]0 1,m a x | ( ) | m a x | (,( ) ) |,x a b x a bM y x f x y x
1| | | [ (,( ) ) (,( ) ) ] |m
m
x
m m mxR f x y x f x y x d x
上一页 下一页湘潭大学数学与计算科学学院 13
几何分析：
Euler公式的误差
ix 1ix?
iy
1iyx? 1iy?
y y x?
( ),i i i iY y x x f x y
||mRR?
记 2( ) / 2R K LM h，则有上一页 下一页湘潭大学数学与计算科学学院 14
整体截断误差,
( ),m m my x y
由：
( ) ( ) (,( ) ),m m m m my x h y x h f x y x R
1 (,),m m m my y hf x y
1 [ (,( ) ) (,) ],m m m m m m mh f x y x f x y R
上一页 下一页湘潭大学数学与计算科学学院 15
从而有：
1| | | | | |,m m mh L R
对任一 0,1,,1,mN
1| | ( 1 ) | |mm h L R
2 2( 1 ) | | ( 1 )mh L h L R R
1
0
0
( 1 ) | | ( 1 )
m
mj
j
h L R h L?
有：
上一页 下一页湘潭大学数学与计算科学学院 16
0( 1 ) | | [ ( 1 ) 1 ]
mm Rh L h L
hL
( ) ( )
0| | ( 1 )
L b a L b aRee
hL?
于是便得 Euler方法的 整体截断误差界
( ) ( )
0| | | | ( ) ( 1 ),2
L b a L b a
m
hKe M e
L
(*)
()1 12 L b am hM eL
0 0 0 0y y x若,
1 ( ),,mM K L M O h这 里 =可 记 为说明 Euler方法的整体截断误差与 h同阶。
上一页 下一页湘潭大学数学与计算科学学院 17
因 此 难 以 使 用 上 述 结 果 。
3 y x M M） 实 际 上,的 估 计 式 中 的 是 不 知 道 的 。
2,） 方 法 的 精 度 由 整 体 误 差 决 定 但 局 部 误 差 容 易 估 计 。
注意,对于 Euler方法
1 ） 一 般 来 说,整 体 误 差 比 局 部 误 差 低 一 阶 。
上一页 下一页湘潭大学数学与计算科学学院 18
定理 5 设 f (x,y)属于?且关于 x满足 Lipschitz条件，
其 Lipschitz常数为 K,且当 时，0h?
0 ( ),y y a?
则 Euler方法一致收敛于真解 ( ),myx
成立,
并且 有估计式
( ) ( )
0| | | | ( ) ( 1 )2
L b a L b a
m
hKe M e
L
(*)
上一页 下一页湘潭大学数学与计算科学学院 19
2.2 改进的 Euler方法
00
(,),,| |,
( ),,
y f x y a x b y
y a y y R
等价于 积分方程：
( ) ( ) (,( ) ),xhxy x h y x f y d
微分方程：
上一页 下一页湘潭大学数学与计算科学学院 20
( ) ( ) (,( ) )xhxy x h y x f y d
有,(令 )mxx?
( ) ( )mmy x h y x (,( ) )m
m
xh
x f y d
11[ (,) (,) ] / 2m m m m m my h f x y f x y R
1 (,( ) ) [ (,( ) )
2
m
m
x
m m mx
hR f x y x dx f x y x
11(,( ) ) ],mmf x y x
去掉,mR
1( ) ( )mmy x y x 11[ (,) (,) ] / 2.m m m mh f x y f x y
有：
上一页 下一页湘潭大学数学与计算科学学院 21
设 为 的近似值,my ()myx
1mmyy 11[ (,) (,) ] / 2,m m m mh f x y f x y
称为 改进的 Euler方法,
称为改进的 Euler方法的 局部截断误差,mR
误差分析,仍记 ( ),m m my x y
注意：
1
3
1( ) [ ( ) ( ) ] ( ),2 1 2
m
m
x
m m mx
hhy x d x y x y x y x h
0 1,
则：
上一页 下一页湘潭大学数学与计算科学学院 22
于是：
若记
[,]m a x | ( ) |,x a bM y x
3 ( ) / 1 2,mmR h y x h
3| | / 1 2,mR h M?
整体截断误差的阶由局部截断误差的阶来决定,
可见改进的 Euler方法误差比 Euler方法要 高一阶,
则有上一页 下一页湘潭大学数学与计算科学学院 23
三,Euler方法的 收敛性和稳定性
立是差分方程的解，若成初值问题的解是设
my
O D Exy,
0 1li m m a x 0mmh mn y y x
的。也称此数值方法是收敛则称差分方程是收敛的不能无限缩小。以实际上对扩大了对解的影响。所误差的积累可能时离散点增多。舍入当
h
h 0?
hOE u l er 方法是收敛的，且阶为结论,
注意,
收敛性上一页 下一页湘潭大学数学与计算科学学院 24
的变化情况。
导致公式的差分方程解稳定性描述初值的改变稳定性
,2,1,
,0
,
,,,,
00
0
000
mvucvu
abmhhh
vu
vuhc
mm
mm
时当满足：的解它们使对任给的初值定义方法是稳定的。则称 E u l er
上一页 下一页湘潭大学数学与计算科学学院 25
结论,
方法是稳定的。则条件，满足若
E u l er
L i p s ch i zyxf,
1,m m m mu u h f x u
mmm vu证明：记
mhL?1
1,m m m mv v h f x v
mmmmmm vxfuxfh,, 1
() 0,L b ae
011 mhL
上一页 下一页湘潭大学数学与计算科学学院 26
计算问题：
1mmyy 11[ (,) (,) ] / 2,m m m mh f x y f x y
隐式计算格式由迭代法去完成,
1 1 1(,) (,) 022m m m m m m
hhy f x y y f x y
将上式变形为记
1 1 1 1( ) (,) (,)22m m m m m m m
hhy y f x y y f x y?
求 1my? 即求隐式方程
1( ) 0my的根。
上一页 下一页湘潭大学数学与计算科学学院 27
总结通过对 Euler方法的讨论可以看到,微分方程数值方法的研究应包括以下方面
1.数值计算公式的构造 ;
2.方法稳定性,收敛性的研究 ;
3.方法的误差估计 ;
4.方法的实现等,
上一页 下一页湘潭大学数学与计算科学学院 28
欧拉公式的改进：
隐式欧拉法向后差商近似导数
h
xyxyxy )()()( 01
1
))(,()( 1101 xyxfhyxy
)1,...,0(),( 111 niyxfhyy iiii
由于未知数 yi+1同时出现在等式的两边，不能直接得到，故称为 隐式 欧拉公式，而前者称为 显式 欧拉公式。
一般先用显式计算一个初值，再 迭代 求解。
隐式 欧拉法的局部截断误差：
11 )( iii yxyR )()( 322 hOxy ih
即隐式欧拉公式具有 1 阶精度。
上一页 下一页湘潭大学数学与计算科学学院 29
梯形公式 — 显、隐式两种算法的 平均
)1,.,,,0()],(),([2 111 niyxfyxfhyy iiiiii
注,的确有局部截断误差,
即梯形公式具有 2 阶精度，比欧拉方法有了进步。
但注意到该公式是 隐式 公式，计算时不得不用到迭代法，其迭代收敛性与欧拉公式相似。
)()( 311 hOyxyR iii
中点欧拉公式中心差商近似导数
h
xyxyxy
2
)()()( 02
1
x0 x2x1
))(,(2)()( 1102 xyxfhxyxy
1,...,1),(211 niyxfhyy iiii
假设,则可以导出即中点公式具有 2 阶精度。
)(),( 11 iiii xyyxyy )()( 311 hOyxyR iii
需要 2个初值 y0和 y1来启动递推过程，这样的算法称为 双步法,
而前面的三种算法都是 单步法 。
上一页 下一页湘潭大学数学与计算科学学院 30
方 法
显式欧拉隐式欧拉梯形公式中点公式简单 精度低稳定性最好 精度低,计算量大精度提高 计算量大精度提高,显式 多一个初值,可能影响精度上一页 下一页湘潭大学数学与计算科学学院 31
改进欧拉法
Step 1,先用 显式 欧拉公式作 预测,算出 ),(1 iiii yxfhyy
Step 2,再将 代入 隐式 梯形公式的右边作 校正,得到1?iy
)],(),([2 111 iiiiii yxfyxfhyy
注,此法亦称为 预测 -校正法 。可以证明该算法具有 2 阶精度，同时可以看到它是个 单步 递推格式，比隐式公式的迭代求解过程 简单 。后面将看到，它的 稳定性高 于显式欧拉法。
)1,.,,,0(),(,),(2 11 niyxfhyxfyxfhyy iiiiiiii
上一页 下一页湘潭大学数学与计算科学学院 32
§ 3 龙格 - 库塔法建立高精度的单步递推格式。
单步递推法的 基本思想 是从 ( xi,yi ) 点出发，以 某一斜率 沿直线达到 ( xi+1,yi+1 ) 点。欧拉法及其各种变形所能达到的最高精度为 2阶 。
考察改进的欧拉法，可以将其改写为：
),(
),(
2
1
2
1
12
1
211
hKyhxfK
yxfK
KKhyy
ii
ii
ii
斜率一定取 K1 K2
的 平均值 吗？
步长一定是一个 h 吗？
上一页 下一页湘潭大学数学与计算科学学院 33
首先希望能确定系数?1,?2,p，使得到的算法格式有 2阶精度，即在 的前提假设下，使得)(
ii xyy?
)()( 311 hOyxyR iii
Step 1,将 K2 在 ( xi,yi ) 点作 Taylor 展开
)(),(),(),(
),(
2
1
12
hOyxfp h Kyxp h fyxf
p h KyphxfK
iiyiixii
ii
)()()( 2hOxyphxy ii
将改进欧拉法推广为：
),(
),(
][
12
1
22111
phKyphxfK
yxfK
KKhyy
ii
ii
ii
),(),(),(
),(),(
),()(
yxfyxfyxf
dx
dy
yxfyxf
yxf
dx
d
xy
yx
yx
Step 2,将 K2 代入第 1式，得到
)()()()(
)]()()([)(
32
221
2
211
hOxyphxyhy
hOxyphxyxyhyy
iii
iiiii
上一页 下一页湘潭大学数学与计算科学学院 34
Step 3,将 yi+1 与 y( xi+1 ) 在 xi 点的 泰勒 展开作比较
)()()()( 322211 hOxyphxyhyy iiii
)()(2)()()( 3
2
1 hOxy
hxyhxyxy
iiii
要求,则必须有,)()( 3
11 hOyxyR iii
2
1,1
221 p
这里有 个未知数,个方程。
3
2
存在 无穷多个解 。所有满足上式的格式统称为 2阶龙格 - 库塔格式 。
2
1,1
21p
注意到,就是改进的欧拉法。
Q,为获得更高的精度，应该如何进一步推广？
上一页 下一页湘潭大学数学与计算科学学院 35
其中?i ( i = 1,…,m ),?i ( i
= 2,…,m ) 和?ij ( i = 2,…,
m; j = 1,…,i?1 ) 均为待定系数，确定这些系数的步骤与前面相似。
)...,(
......
),(
),(
),(
]...[
112211
23213133
12122
1
22111
mm mmmmim
ii
ii
ii
mmii
hKhKhKyhxfK
hKhKyhxfK
hKyhxfK
yxfK
KKKhyy
最常用为四级 4阶 经典龙格 -库塔法,
),(
),(
),(
),(
)22(
34
2223
1222
1
432161
hKyhxfK
KyxfK
KyxfK
yxfK
KKKKyy
ii
h
i
h
i
h
i
h
i
ii
h
ii
上一页 下一页湘潭大学数学与计算科学学院 36
注：
龙格 -库塔法 的主要运算在于计算 Ki 的值，即计算 f 的值。 Butcher 于 1965年给出了计算量与可达到的最高精度阶数的关系：
753
可达到的最高精度
642每步须算 Ki 的个数
)( 2hO )( 3hO )( 4hO )( 5hO )( 6hO)( 4hO )( 2?nhO
8?n
由于龙格 -库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响。对于光滑性不太好的解，最好采用 低阶算法 而将步长 h取小 。
上一页 下一页湘潭大学数学与计算科学学院 37
§ 3 收敛性与稳定性
收敛性定义 若某算法对于任意固定的 x = xi = x0 + i h,当 h?0
( 同时 i) 时有 yi? y( xi ),则称该算法是 收敛 的 。
例,就初值问题 考察欧拉显式格式的收敛性。
0)0( yy
yy?
解,该问题的精确解为 xeyxy?
0)(?
欧拉公式为
iiii yhyhyy )1(1 0)1( yhy ii
对任意固定的 x = xi = i h,有
ii xhhxi hyhyy ])1[()1( /10/0
eh hh /10 )1(l i m
)(0 ix xyey i
上一页 下一页湘潭大学数学与计算科学学院 38
定义 若某算法在计算过程中任一步产生的误差在以后的计算中都 逐步衰减,则称该算法是 绝对稳定的 。
一般分析时为简单起见，只考虑 试验方程
yy
常数，可以是复数当步长取为 h 时，将某算法应用于上式，并假设只在初值产生误差,则若此误差以后逐步衰减，就称该算法相对于 绝对稳定,的全体构成 绝对稳定区域 。
我们称 算法 A 比算法 B 稳定,就是指 A 的绝对稳定区域比 B
的 大 。
000 yy
h? h? h
上一页 下一页湘潭大学数学与计算科学学院 39
例,考察显式欧拉法
011 )1( yhyhyy iiii
000 yy 011 )1( yhy ii
01111 )1( iiii hyy
由此可见，要保证初始误差?0 以后逐步衰减，
必须满足：
hh
1|1| h
0-1-2 Re
Img
例,考察隐式欧拉法
11 iii yhyy?
ii yhy
1
1
1 0
1
1 1
1
i
i h
可见绝对稳定区域为,1|1| h
210 Re
Img
注,一般来说，隐式欧拉法的绝对稳定性比同阶的显式法的好。
上一页 下一页湘潭大学数学与计算科学学院 40
例,隐式龙格 -库塔法
),.,,,1(
).,,,(
].,,[
11
111
mj
hKhKyhxfK
KKhyy
mmjjijij
mmii
而 显式 1~ 4 阶方法的绝对稳定区域为
)2,2( 11
11
KhyhxfK
hKyy
ii
ii其中 2阶方法 的绝对稳定区域为
0 Re
Img
k=1
k=2
k=3
k=4
-1-2-3
-
-
-
1
2
3
Re
Img
无条件稳定
