Transcript 浏览观看

第十一章 常微分方程初值问题的
数值解法
11.1 欧拉方法
11.2 计算公式的误差分析
11.3 龙格－库塔方法
11.4 向一阶方程组与高阶方程的推广
问题：
 y  f ( x, y )
求解 
x  [ x0 , b]
 y( x0 )  y0
(11 1)
数值求解方法：
对区间[ x0 , b] 作等距分割, 生成节点xi  x0  ih (i  0,1,, n),
b  x0
, 然后逐个求解出节点上的函数值 y ( x1 ), y ( x2 ),,
N
y ( xN ) 的近似值
h
11.1 欧拉方法
11.1.1 欧拉公式与改进欧拉公式
y ( xi 1 )  y ( xi )
 y( xi )  f ( xi , y ( xi ))
h
得 : y ( xi 1 )  y ( xi )  h f ( xi , y ( xi ))
由
由此可建立计算公式
 yi 1  yi  h f ( xi , yi )

 y0  y ( x0 )
(i  0,1,  , N  1)
(11 2)
这称为欧拉公式
例11.1 以 h=0.1为步长，用欧拉法求常微分方程初值问题
 y   xe x  y

 y (0)  1
x  [0,1]
的数值解 , 并与精确解 y ( x) 
1 2
( x  2)e  x 比较
2
y ( xi 1 )  y ( xi )
 y( xi 1 )  f ( xi 1 , y ( xi 1 ))
h
得 : y ( xi 1 )  y ( xi )  h f ( xi 1 , y ( xi 1 ))
类似地,由
由此可建立另一公式
 yi 1  yi  h f ( xi 1 , yi 1 )

 y0  y ( x0 )
(i  0,1,  , N  1)
(11 3)
这称为后退欧拉公式
后退欧拉公式是一个隐式公式，通常采用迭代法求解。
11.1.2
梯形公式与改进欧拉公式
欧拉公式与后退欧拉公式也可采用积分近似的方法推出
对微分方程积分y  f ( x, y ) 得
y ( xi 1 )  y ( xi )  
xi 1
xi
y( x)dx  
xi 1
xi
f ( x, y ( x))dx
从而把微分方程问题转化为积分方程问题
y ( xi 1 )  y ( xi )  
xi 1
xi
由
xi 1
xi
f ( x, y ( x))dx
f ( x, y ( x))dx  
xi 1
xi
f ( xi , y ( xi ))dx  h f ( xi , y ( xi ))
得y ( xi 1 )  y ( xi )  h f ( xi , y ( xi ))
从而可导出欧拉公式
(11 4)
由
xi 1
xi
f ( x, y ( x))dx  
xi 1
xi
f ( xi 1 , y ( xi 1 ))dx  h f ( xi 1 , y ( xi 1 ))
得y ( xi 1 )  y ( xi )  h f ( xi 1 , y ( xi 1 ))
从而可导出后退欧拉公式
xi 1
h
由  f ( x, y ( x))dx  [ f ( xi , y ( xi ))  f ( xi 1 , y ( xi 1 ))]
xi
2
h
得 y ( xi 1 )  y ( xi )  [ f ( xi , y ( xi ))  f ( xi 1 , y ( xi 1 ))]
2
从而导出梯形公式
h

y

y

[ f ( xi , yi )  f ( xi 1 , yi 1 )]
 i 1
i
2

 y0  y ( x0 )
梯形公式也是隐式单步法公式
(6  5)
用梯形公式计算时，通常取欧拉公式的解作为迭代初值进行迭
代计算，即采用下式
 yi(01)  yi  h f ( xi , yi )

(k  0,1, )
 ( k 1)
h
(k )
 yi 1  yi  [ f ( xi , yi )  f ( xi 1 , yi 1 )]
2

迭代计算到某个k , 使得 | yi(k1)  yi(k11) |  时, 取 yi(k1) 作为函数
值 y( xi 1 )的近似值 yi 1
简单地, 直接 取 yi(11) 代替 yi 1则得到显式计算公式
yi 1  yi  h f ( xi , yi )
~

(i  0,1,, N  1)

h
~
 yi 1  yi  2 [ f ( xi , yi )  f ( xi 1 , yi 1 )]
这称为改进欧拉公式
例11.2 仍取步长h = 0.1，采用改进欧拉法重新计算例 6.1 的
常微分方程初值问题。
解
这时改进欧拉公式为
~
yi 1  yi  h ( xi e  xi  yi )


h
 xi
 xi 1
~
y

y

[(
x
e

y
)

(
x
e

yi 1 )]
 i 1
i
i
i
i 1
2

计算结果见表6-2（书125页）
(i  0,1,,9)
11.2 计算公式的误差分析
定义11.1 若 yi+1 是 yi=y(xi) 从计算得到近似解，则称
y(xi+1) － yi+1为所用公式的局部截断误差，简称为截断误差。
定理11.1 若单步法 yi+1 = yi+h (xi , yi , h) 的局部截断
误差为 O (h p+1) ，且增量函数 (x , y , h) 关于 y 满足李普希
兹条件，即存在常数 L>0，使对  y, ~
y 成立不等式
|  ( x, y, h)   ( x, ~
y , h) | L | y  ~
y|
则其整体截断误差 y(xi)－ yi=O(hp)
截断误差的估计(基本假设： yi = y( xi ) )
设 y(x)C 3 [x0 , b] , 则
h2
y( xi 1 )  y( xi  h)  y( xi )  hy( xi ) 
y( xi )  O(h3 ) (11 7)
2
（1）对欧拉公式，有
yi 1  y ( xi )  h f ( xi , y ( xi ))  y ( xi )  hy( xi )
故
h2
y ( xi 1 )  yi 1 
y( xi )  O(h 3 )  O(h 2 )
2
(11 8)
因此，欧拉公式的局部截断误差为 O (h2)
（2）对后退欧拉公式，有
yi 1  y ( xi )  h f ( xi 1 , y ( xi 1 ))
故
h2
y ( xi 1 )  yi 1   y( xi )  O(h3 )  O(h 2 )
2
因此，后退欧拉公式的局部截断误差为 O (h2)
(11 9)
（3）对梯形公式，注意到其公式可改写为
1
yi 1  [ yi  hf ( xi , yi )]  [ yi  h f ( xi 1 , yi 1 )]
2
故由式（6-9）和（6-9）得
y ( xi 1 )  y i 1
1
  y ( xi 1 )  [ y i  h f ( xi , y i )]
2
 y ( xi 1 )  [ y i  h f ( xi 1 , y i 1 )]

2
1  h2
h
3
3 
 [ y ( xi )  O (h )]  [
y ( xi )  O(h )]
2 2
2

 O(h 3 )
因此，梯形公式的局部截断误差为 O ( h3 )
(6  10)
（4）对改进欧拉公式，有
yi 1  y( xi )  h y( xi )
~


h
~
y

[
f
(
x
,
y
)

f
(
x
,
i 1
i
i
i 1 yi 1 )]

2

而由 y  f ( x, y) 得 y  f x( x, y)  yf y( x, y) ，故有
f ( xi 1 , ~
yi 1 )  f ( xi  h , y ( xi )  h y( xi ))
 f ( xi , y ( xi ))  [h f x  h y f y ]( xi , y ( xi ))  O(h 2 )
 y( xi )  h y( xi )  O(h 2 )
h2
所以 yi 1  y ( xi )  h y( xi ) 
y( xi )  O(h 3 )
2
与式（6-7）比较得 y(xi+1) －yi+1 = O ( h3 )
因此，改进欧拉公式的局部截断误差为 O ( h3 )
定义11.2 若一种求解常微分方程初值问题的数值计算
方法的局部截断误差为 O ( hp+1 ) ，则称该方法为 p阶精度，
或称该方法为 p阶方法。
由此定义知，欧拉方法与后退欧拉方法为一阶精度，梯
形法与改进欧拉方法为二阶精度。
11.3 龙格-库塔方法
由中值定理，有
y( xi 1 )  y( xi )  ( xi 1  xi ) y( )
 h f ( , y( )),   ( xi , xi 1 )
因此，以上介绍的各种单步法本质上都是对平均斜
率 f( , y( )) 进行近似，龙格-库塔据之提出了适当选取若
干点上的斜率值作近似以构造高精度计算公式的方法，其
基本思想是基于泰勒展式的待定系数法。
11.3.1
二阶R-K公式
问题：建立二阶精度的计算格式形为
 yi 1  yi  h(1 K1  2 K 2 )

 K1  f ( xi , yi )
 K  f ( x  ah, y  bhK )
i
i
1
 2
(11  12)
这里1 , 2 , a , b 为待定系数
解
在 y(xi) = yi 的假设下，有
K1  f ( xi , y ( xi ))  y( xi )
K 2  f ( xi  ah, y ( xi )  bhK1 )
 f ( xi , y ( xi ))  [ahf x  bhK1 f y ]( xi , y ( xi ))  O(h 2 )
 y( xi )  h [a f x  byf y ]( xi , y ( xi ))  O(h 2 )
故
yi 1  y( xi )  (1  2 )hy( xi )  h2[a2 f x  b2 y f y]( xi , y ( xi ))  O(h3 )
yi 1  y( xi )  (1  2 )hy( xi )  h2[a2 f x  b2 y f y]( xi , y ( xi ))  O(h3 )
2
3
而 y ( xi 1 )  y ( xi )  hy( xi )  (h / 2) y( xi )  O(h )
1
2 1


 y ( xi )  hy ( xi )  h [ f x  y f y ]( xi , y ( xi ))  O(h 3 )
2
2
根据格式为二阶精度，即 y(xi+1) －yi+1 = O(h3) 比较两式系数得
1  2  1

(11 13)
a2  1 / 2
b  1 / 2
 2
系数满足(6-13)的形为(6-12)计算格式统称为二阶R-K公式。
当令1=1/2时，解得 2=1/2 ，a=b=1，即为改进欧拉公式。若
令 1=0，解得 2=1，a=b=1/2，则得另一计算公式
 yi 1  yi  hK2

 K1  f ( xi , yi )
 K  f ( x  h / 2, y  hK / 2)
i
i
1
 2
(11 14)
变形欧拉公式
四阶 R-K 公式
1965年，Butcher研究发现显式R-K公式的精度与需要组
合的斜率值的个数具有如下关系
11.3.2
7 n
8
精度阶
1 2 3 4 4 5 6 n-2
可见，超过四阶精度的R-K公式效率并不高，实际计算通
每一步需计算的 f 值的个数
1
2
3
4
5
6
常选用如下四阶格式
h

 yi 1  yi  6 ( K1  2 K 2  2 K 3  K 4 )

 K1  f ( xi , yi )

 K 2  f ( xi  h / 2, yi  hK1 / 2)
 K  f ( x  h / 2, y  hK / 2)
i
i
2
 3
 K 4  f ( xi  h, yi  hK3 )

(11 15)
经典R-K公式
例11.3 取步长h = 0.2，采用经典R-K法计算例 6.1 的常微
分方程初值问题。
解 这时经典R-K公式为
h

y

y

( K1  2 K 2  2 K 3  K 4 )
i
 i 1
6

x
 K1  xi e i  yi

h ( xi  h / 2 )
h

 ( y i  K1 )
 K 2  ( xi  )e
2
2

h  ( xi  h / 2 )
h

 ( yi  K 2 )
 K 3  ( xi  2 )e
2

 ( xi  h )
K

(
x

h
)
e
 ( yi  hK3 )
 4
i
取 h=0.2 计算得到表6-4(书133页）。
与例6.1和例6.2比较可见，用经典R-K法计算得到的解比用
欧拉法和改进欧拉法所得到的解精确得多。
11.3.3 步长的自动选择
为便于进行事后误差估计，实际计算时通常采用步长减半
算法。
对于 p 阶精度的计算格式，当取步长为 h 时，记 y [ h ] 为
从 y(xi) 计算得到的 y (xi+1) (xi+1= xi+h) 的近似解，则有
y ( xi 1 )  yi[h1]  Ch p 1
y ( xi 1 )  y
~ h 
 2C  
2
[ h / 2]
i 1
(C为常数)
p 1
~
(C为常数)
~
h 较小时, C  C 故有
y ( xi 1 )  yi[h1/ 2]
1
 p
[h]
y ( xi 1 )  yi 1
2
所以 y ( xi 1 )  y
[ h / 2]
i 1
i 1
1
 p
( yi[h1/ 2]  yi[h1] )
2 1
记  | yi[h1/ 2]  yi[h1] ，则对给定的精度要求
 ，可根据  按如
|
下方式调整步长：
（1）若 > ，则把步长逐次减半计算，直至  < 为止，
这时最终得到的解即为满足精度要求的近似解。
（2）若 < ，则把步长逐次加倍计算，直至   为止，
这时取前一次步长计算所得到的解作为满足精度要求的近似解。