常微分方程初值问题的数值解法

y' (xi ) phy'' (xi ) (h2 )代入（7 -10）得
yi1 y(xi ) h(1 2 ) y' (xi ) 2 ph2 y'' (xi ) (h3 )
又yx在xi 2的二阶泰勒展开式为
y(x)
y(xi )
y' (xi )(x
xi )
y
'' (xi ) 2！
(
y(xi1) h
y(xi )
y' (xi
h)
f
( xi
h, y(xi
h))
所以
y(xi1) y(xi ) hf (xi h, y(xi h)) （7 - 9）
这里y' (xi
h)
f
( xi
h, y(xi
h))
_
k
可以看作
yx在xi , xi1上的平均斜率，有微分中值定理知
（7 - 9）式是精确表达式，但是一般无法求得y' (xi h)
1.183216 1.341641
步长的自动选择
单从第一步看，步长h越小，局部截断误差就越小； 但随着h的缩小，不但引起计算量的增加， 而且也引起舍入误差的严重积累。
所以选取合适的步长h在实际计算中是很重要的。
以四阶为例，从节点xi出发，步长h，求出近似值记为
yh i 1
,
然后折半取步长h
/
2，从xi
h
h
k2 f (xi 2 , yi 2 k1)
同理可得三阶龙格-库塔公式
yi1 yi h 1k1 2k2 3k3
k1
f
xi ,
yi
k2 f xi ph, yi phk1
k3 f xi qh, yi qh rk1 sk2
yi1
k1 f
yi
xi ,
Taylor展开法不仅可得到求数值解的公式，且容易 估计截断误差。
7.2 尤拉方法 1. Euler 公式
初值问题：
y f (x, y) a x b (7 1)
y(x0 ) y0
(7 2)
xn1 y(t)dt xn1 f (t, y(t))dt
xn
xn
y(xn1) y(xn )
,因此，平均斜率k_ 采用不同近似方法，就可以得到
不同的计算公式。
启示
如果在区间内多取几个点的斜率值，然后把它们的线性组合作为平均斜率 的近似值，则有可能构造出更高精度的计算公式，这就是龙格-库塔法的基 本思想
二阶龙格-库塔公式
yi1 yi h(1k1 2k2 )
k1 f (xi , yi ) k2 f (xi ph, yi phk1)
定解问题：微分方程＋定解条件(初值条件、边界条件) 分别称为初值问题和边值问题。
一阶常微分方程的初值问题：
y f (x, y) a x b (7 1)
y(x0 ) y0
(7 2)
函数f (x, y)在带形区域
(x, y) a x b, y
x, y为连续函数，且对任意的y满足李普
逐步计算出y （y x）在xn1点的值： yn1 yn hf (xn, yn ) n 0，1，2，L
用分段的折线逼近函数,此为“折线法”而非“切线 法”，除第一个点是曲线上的切线，其它都不是。
（2）梯形公式
y(xn1) y(xn )
xn1 f (t, y(t))dt
xn
(7 3)
称该方法具有p阶精度或p阶方法。显然，尤拉公式 具有一阶精度，梯形公式具有二阶精度。
3. 改进的尤拉方法
梯形公式虽然提高了精度，但算法复杂。而在实际计 算中只迭代一次，这样建立的预测—校正系统称作改 进的尤拉公式。
预测 校正
yn1 yn hf (xn , yn )
yn1
yn
h[f 2
(xn , yn )
k2
f
x0
h 2
,
y0
h 2
k1
0.918182
k3
f
x0
h 2
,
y0
h 2
k2
0.908637
k4 f x0 h, y0 hk3 0.843239
所以y1
y0
h 6
k1
2k2
2k3
k4
1.183229
在计算y2
k1 f x1, y1
k2
f
x1
h 2
,
y1
h 2
k1
将y(xn1)在xn点 Taylor 展开：
y ( xn1 )
y ( xn
h)
y(xn )
hf
(xn ,
y(xn ))
h2 ( x2n
y( )
xn1)
yn1 yn hf (xn , yn )
Rn1
y(xn1)
yn1
h2 2
y( ) O(h2 )
定义：
数值方法的局部截断误差为：Rn1 O(h p1)
a x0 x1 x2 L xn b 在节点上用离散化方法将连续型微分方程 转化成离散型代数方程即差分方程来求解。
具体作法：利用y(x0 )求出y(x1)的近似值y1，再由 y1求出y2，L ，直到求出yn为止。该算法称为步进 式或递推式算法。
微分方程离散化常用方法
（1）数值积分 用数值积分方法离散化 :
（7-10）
选取适当的 1，,2，p使其截断误差为o h3 .
选取过程
假定对
k1
,
k
做泰勒展开
2
k1 f ( xi , yi ) y '
k2 f (xi p , yi phk1) f (xi ph, yi phk1)
f (xi , yi ) phfx' (xi , yi ) phk1 f y' (xi , yi ) (h2 )
计算两步到xi
1
,
得近似值
yh i 1
,
记
y h / 2 i 1
y(h) i 1
来判断选取步长是否合适，分以下
(7 5)
y(k 1) n1
yn
h(f 2
(xn ,
yn )
f
(
xn
1
,
y(k) n1
)
k
0,1,2,L
2. Euler方法的截断误差
(1) 局部截断误差 在一步中产生的误差而非累积误差：Rn1 y(xn1) yn1
其中，yn1是当yn y(xn )(精确解)时由Euler法求出 的值，即yn无误差。
xn1 f (t, y(t))dt
xn
(7 3)
对右端积分采用取左端点的矩形公式,则有：
xn1 xn
f
(t,
y(t)dt
h
f
(xn ,
y(xn ))
h2 2
y(n )
(n 0,1,L ), xn n xn1
当h充分小，舍去误差项R1
h2 2
y(i )，则：
y(xn1) y(xn ) h f (xn , y(xn ))
xn1 dydx xn1 f (x, y)dx (n 0,1,L )
xn dx
xn
用yn1, yn代替y(xn1), y(xn ), 对右端积分采用 取左端点的矩形公式
则有
xn1 xn
f
(x,
y)dx
hf
(xn ,
yn )
yn1 yn hf (xn , yn ) (n 0,1,L )
f
(xn1, yn1)]
也可表示为
yn1
yn
h[ 2
f
(xn ,
yn )
f
( xn
h,
yn
hf
(xn ,
yn ))]
或者表示为下列平均化形式
yp yc
yn hf (xn , yn hf (xn1
yn ,y
)
p
)
yn1
(
yp
yc
)
/
2
(7 8)
例1：用尤拉公式和改进的尤拉公式解初值问题
（2）数值微分
dy dx
( xn , yn )
y(xn1) y(xn ) xn1 xn
f (xn , y(xn ))
用 h xn1 xn , yn y(xn ), yn1 y(xn1) 代替，
则：
yn1 h
yn
f (xn , yn )
yn1 yn hf (xn , yn ) n 0, 1, 2,L
由于f (x0 , y0 ) 及(x0, y0 ) 已知，必有切线方程。 由点斜式写出切线方程：
dy y y0 (x x0 ) dx (x0 ,y0 ) y0 (x x0 ) f (x0 , y0 )
等步长为h，则x x0 h，可由切线算出y1
y1 y0 hf (x0 , y0 )
h 6
k1
yi
4k2
k3
k2
f
xi
h 2 , yi
h 2
k1
k3 f xi h, yi hk1 2hk2
局部截断误差o(h4)
四阶龙格-库塔公式
yi1
yi
h 6
k1
2k2
2k3
k4
k1
f
xi , yi
k2
f
xi
h 2
,
yi
h 2
k1
k3
f
xi
h 2 , yi
y
y
2x y
(0 x 1)
y(0) 1
解：取步长h 0.1，
尤拉公式为：yn1
yn
h( yn
2 xn yn
)
y
p
yn
h( yn
2 xn yn
)
改进的尤拉公式为： yc
yn
h( yp
2xn1 ) yp
yn1
1 2
(
yp
yc
)
7.3 龙格-库塔（R-K）法
尤拉法
yi1 yi hk1
k3
f
x1
h 2
,
y1
h 2
k2
k4 f x1 h, y1 hk3
所以
y2
y1
h 6
k1
2k2
2k3
k4
1.341667
这样继续下去，计算结果列于表
基点 0 0.2 0.4
改进尤拉法 1
1.184096 1.343360
龙格-库塔法 1
1.183229 1.341667
精确解 1
第7章 常微分方程初值问题的数值解法
7.1 引言 7.2 尤拉方法 7.3 龙格—库塔法 7.4 收敛性和稳定性
7.1 引 言
ODE的求解： 分离变量法、齐次方程的求解、可降阶高阶微分方程 求解——特殊类型的微分方程。
微分方程的近似解法： (1) 近似解析法：逐次逼近法、级数解法 (2) 数值解法：求离散点上的近似值。
Βιβλιοθήκη Baidu
x
xi )2
y''' ()
3！ (x
xi )3
介于x与xi之间
当x
y ( xi1
xi
) y
时，y
1
''' ( )
3！
(
(xi ) hy' (xi )
x
h2 2！y
xi )3
'' (xi )
(h
(h3)
3
)，有
（ ）
比较（）和（ ）的系数即可发现，要使公式具有
二阶精度只要下列条件成立
1 2 1
2
p
1 2
这里有3个未知数，
2 个方程。
7-11
存在无穷多个解
将满足条件（7-11）式到一簇公式（7-10） 统称为二阶龙格-库塔公式
p
1,1
2
1时，二阶龙格- 库塔公式 改进尤拉公式 2
p
1 2
,
1
0, 2
1,二阶龙格- 库塔公式
变形尤拉公式
中点公式
yi1 yi hk2
k1 f (xi , yi )
k1 f xi , yi
局部截断误差O(h2)
考察改进的尤拉法，可以将其改写为：
yi 1
yi
h
1 2
K1
1 2
K
2
K1 f ( xi , yi )
K 2 f ( xi h, yi hK1 )
局部截断误差O(h3)
增加计算f(x,y)在不同点的 值，能否提高局部截断误差的阶？
有微分中值定理，存在 （0,1）使得
用yn1, yn代替y(xn1), y(xn ),得出初值问题（7 -1）（7 - 2） 的离散化形式：
y0 y(x0 )
(7 4)
yn1 yn h( f (xn , yn ) (n 0,1,2,L )
(1) 几何意义 由 (x0, y0 ) 出发取曲线 y y(x)的切线（存在)，则斜率 dy dx (x0 , y0 ) f (x0 , y0 )
积分用梯形公式，且令：yn1 y(xn1), yn y(xn )
则得：yn1
yn
h( 2
f
(xn ,
yn )
f
(xn1, yn1))
Rn1
y(xn1)
yn1
h3 12
y( )
xn xn1
与Euler法结合，形成迭代算法，对n 0，1，2，L
y(0) n1
yn
hf
(xn ,
yn )
（3）在xn附近 y(x)的 Taylor 展开：
y( xn
h)
y(xn ) hy(xn )
h2 2
y(xn ) L
y(xn ) hf
(xn ,
y(xn ))
h2 2
y(xn ) L
取 h 的线性部分，且yn y(xn ) 得 y(xn )的近似值： yn1 yn hf (xn y(xn ) n 0,1, 2,L
h 2
k2
k4 f xi h, yi hk3
局部截断误差o(h5)
例题 用四阶龙格-库塔公式求解初值问题
y'
x
y
2x y
y 0 1
从x 0到x （1 取h=0.2）各节点上的数值解
f x, y y 2x , h 0.2
y
由 x0 0, y0 1
得k1 f x0 , y0 1
希兹(Lipschitz)条件，即存在常数L，使得
f (x, y1) f (x, y2 ) L y1 y2 由常微分方程理论知，初值问题(7 -1)(7 - 2)的 解必存在且唯一。
数值方法的任务：寻求唯一存在的解在离散节点 上的近似值yi(i 0,1, 2,L )。 数值解法的基本思想：在初值问题(7 -1)(7 - 2)解 的区间a x b上插入一系列节点