数值分析第九讲(2)

0 yi 1 yi 1 k yi k h0 fi 1 fi 1 k fi k
或写为：

j 0
k
j
yi j h j f i j ,
j 0
k
i 0,1,2,, n k.
9
(9.15)
结束
其中αj,βj(j=0,1,…,k)都是实常数,且αk≠0, ｜α0｜+｜β0｜≠0, fi+j= f (xi+j,yi+j), j=0,1,…,k, 由(9.15)可看出要计算yi+k,要利用它前面的k个值yi,yi+1,…,yi+k-1,又 因为(9.15)关于yi+j和 f i+j都是线性组合,所以这一类方法都称为线 性k步法.欧拉法,隐式欧拉法和梯形法都是线性一步法,欧拉中点 公式是线性二步法.
的形式,通常称为改进的欧拉法. (2)取c1=0,c2=1,a2=1/2 得:
这公式又称中点公式.我们还可以构造其他的二阶R-K法. 6 结束
yi 1 yi hK 2 K1 f ( x i , yi ) h h K 2 f x i , yi K1 2 2
p ( xi , yi , h) c j K j j 1 K1 f ( xi , yi ) j 1 K j f ( xi a j h, yi h b jl K l ) l 1
j 2,3,, p
其中cj,aj,bjl是待定参数, aj和bjl满足
yi 1 y( xi ) c1 c2 hy( xi ) a2c2h2 y( xi ) O(h3 )
而Taylor展开式 代入(9.10) 代入(9.11)
h2 y ( xi 1 ) y ( xi ) hy( xi ) y( xi ) O(h3 ) 2
(9.13)
9.3.2 四阶 R-K 法 用类似的方法可以确定三级和四级R-K法的参数,构造出三阶和 四阶的R-K法.但最常用的是四阶R-K法,四阶R-K法也不只一个, 下面给出的是最常用的四阶经典的R-K公式：
h yi 1 yi 6 K1 2 K 2 2 K 3 K 4 K f ( x , y ) i i 1 K f x h , y h K i i 1 2 2 2 h h K f x i , yi K 2 3 2 2 K f x h, y hK i i 3 4
令L( y( x), h) j y( x jh) h j y( x jh) ,
j 0
k


把y(x+jh)和y’(x+jh)作Taylor展开：
jh y( x jh) y( x) jhy( x)
2!
2百度文库
jh y( x)
p 1
y( x) y( x) f x f y
所以
f ( xi a2h, y( xi ) a2hy( xi )) y( xi ) a2hy( xi ) O(h2 )
( xi , yi , h) c1 c2 y( xi ) a2hy( xi ) O(h2 )
公式(9.19)称为米尔尼(Milne)公式,它的局部截断误差为：
14 5 ( 5) R h y ( xi ) O ( h 6 ) 45
它也可以写成
4h yi 1 yi 3 (2 f i 2 f i 1 2 f i ) 3
14 5 ( 5 ) R h y ( ) 45


i 0,1,2,, n 1 (9.14)
7
结束
例3 用经典的四阶R-K法计算例2题目,取步长为0.2,且与准确值比 较. 计算结果列入表9-3：
可见即使用h=0.2计算,也比一阶和二阶方法精度好得多
8
结束
§9.4 线性多步法
单步法只利用前一步的结果,只要给出初值,就能开始计算,但也因 为它只利用前一步的值,为了提高精度就要计算一些非结点处的函 数值,增加了计算量.R-K法就是通过这一途径提高精度的.下面介 绍的线性多步法,在求yi+1时,不仅用到yi的值, 还用到前若干步的 yi-1,…,yi-k的值,这些值都是已知的,因此可在计算量增加不多的情况 下提高精度. 9.4.1 用待定系数法构造线性多步法 线性多步法的一般形式是：
15
i 2,3, , n 1 (9.20)
(9.21)
结束
例5 试确定下列公式的系数和局部截断误差
yi+1=a yi+b yi-2+h (c fi+1+d fi+e fi-1) 解 按(9.15)和(9.16)可知k=3,α0= - b,α1=0,α2= - a,α3=1,
4 结束
二式相减,得局部截断误差
Ri 1 y ( xi 1 ) yi 1 1 2 1 c1 c2 hy ( xi ) a2 c2 h y ( xi ) O(h3 ) 2 1 令1 c1 c2 0, a 2 c2 0 得: 2 c1 c 2 1 a c 1 2 2 2
就是一种p阶的线性k步方法,Cp+1hp+1y(p+1)(xi)称为局部截断误差 的主项.当βk=0时,是显式方法,当βk≠0时是隐式方法. 下面构造几个实用的线性多步法公式 例4 形如：
yi+4= -α0yi+h(β1fi+1+β2fi+2+β3fi+3)的线性4步法公式,试确定 α0,β1,β2,β3并求其局部截断误差主项.
只要c1,c2,a2满足以上方程,就得到一个二阶的R-K法. 这是一个不定方程,有无穷多解.比如： 5 结束
(1)取c1=c2=1/2,a2=1得
h yi 1 yi 2 K1 K 2 (9.12) K1 f ( x i , yi ) K f ( x h, y hK ) i i 1 2 这实际上是 (9.9)公式,即梯形公式的预估-校正公式只迭代一次
k
C2 1 22 k 1 22 kk
1 2 k2 2 3
C
1 6 1
C1 1 22 kk 0 1 k
2 3
2

k2 6
2
由二元函数的泰勒展开：
f ( xi a2 h, y ( xi ) a2 hy( xi )) O( h 2 ) f ( xi , y ( xi )) a2 hf x a2 hy( xi ) f y
其中所有的偏导数都是它们在点(xi,y(xi))的值,下同 3 结束
又由于： y( x) f ( x, y ( x))
此时
( xi , yi , h) c1 f ( xi , y( xi )) c2 f ( xi a2 h, y( xi ) a2 hf ( xi , y( xi ))) c1 y( xi ) c2 f ( xi a2 h, y( xi ) a2 hy( xi )) (9.11)
§9.3 龙格-库塔(Runge-Kutta)法
欧拉方法是显式的一步法,使用方便,但精度较低.本节将构造出 高精度的显式一步法：龙格-库塔法,简称R-K法. 9.3.1 二阶R-K法 欧拉法的公式为： yi+1=yi+h f (xi ,yi) i=0,1,2,…,n-1 决定其精度的是函数f (xi ,yi).如能改进这个函数,就可能提高公 式的精度.为此把公式改写成： yi+1=yi+h (xi ,yi,h) i=0,1,2,…,n-1 (9.10) 选择函数 (xi ,yi,h),一种方法是用若干个点的函数值的线性组 合代替 (xi ,yi,h),如： 1 结束
p
jh y( x jh) y( x) jhy( x)
10
p!
y ( p ) ( x) ,
( p 1)!
y ( p ) ( x) ,
结束
于是, L( y ( x), h) j p j j p 1 j p ( p ) h y ( x) j y ( x) j j j hy( x) p! ( p 1)! j 0 C0 y ( x) C1hy( x) C p h p y ( p ) ( x) , 其中 C0 0 1 k
j 0 j 0
k 1
k
C p 1h p 1 y ( p 1) ( x) O ( h p 2 )
即
y ( xi k ) j y ( xi j ) h j f ( xi j , y ( xi j ))
j 0 j 0
k 1
k
C p 1h p 1 y ( p 1) ( x) O ( h p 2 )
13
结束
解: 由(9.15)知α1=α2=α3=0,α4=1,β0=β4=0,因为有4个待定系数,由 (9.16)写出前4个方程：
C0 0 1 0 C 4 ( ) 0 1 2 3 1 42 C2 ( 1 2 2 3 3 ) 0 2! 43 1 2 2 2 32 3 C3 0 3! 2!
j 0
k


C p 1h p 1 y ( p 1) ( x ) O ( h p 2 )
将x=xi代入上式,设αk=1,并注意到y(xi+jh)=f(xi+jh,y(xi+jh)),可推 出
y ( xi k ) j y ( xi jh) h j f ( xi jh, y ( xi jh))
12
(9.17)
结束
设yi=y(xi),yi+1=y(xi+1),…,yi+k-1=y(xi+k-1),记左端为
yi+k,并舍去最后两项,(9.17)变为：
yi k j yi j h j f ( xi j , yi j )
j 0 j 0
k 1
k
(9.18)
k 1 2 1 2
k2 2
k


(9.16)
结束
Cp 1 2 2 k k / p! 1 2 2 k k /( p 1)!
p p p1 p1


11
若选择αj,βj,使C0=C1=…=Cp=0,Cp+1≠0,则
L( y ( x ), h) j y ( x jh) h j y ( x jh)
a j b jl
l 1
j 1
j 2,3,, p
以上方法称为p级R-K法,选择cj,aj和bjl,可能使以上方法为p阶 方法.显然欧拉法就是一阶R-K法.
2
结束
二级R-K法的形式是:
y i 1 y i hc1K 1 c 2 K 2 K1 f ( xi , yi ) K f ( x a h, y a hK ) i 2 i 2 1 2
解之,得α0=-1,β1=β3=8/3,β2=-4/3 故所求的公式为
yi 4
4h yi (2 f i 1 f i 2 2 f i 3 ) 3
14
(9.19)
结束
将α0,β1,β2,β3代入C4
1 23 2 33 3 44 C4 0 4! 3! 再求C5 1 2 4 2 34 3 45 14 C5 5! 4! 45