计算方法龙格库塔方法.

其中 yi f , yi f (xi , yi )x f x ff y
yi ( f x ff y )x f xx 2 f xy f f yy f 2 f x f y ( f y )2 f
2020/3/18
2
显然p=1时, y i+1=y i+hf(xi,y i)
它即为我们熟悉的Euler方法。
14
注意到，a2 b21 1, c1就是c2二阶12龙格 - 库塔公式，也就
是改进的欧拉法。
y
K2
yi 1
yi
Biblioteka Baidu
1 2
K1
K2
K1 hf (xi , yi )
y y(x)
K1
K
K2 hf (xi h, yi K1)
xi
xi1
x
因此，凡满足条件式有一簇形如上式的计算格
式，这些格式统称为二阶龙格—库塔格式。因此改
y(xi )
hf
(xi ,
yi )
h2 2!
fx ff y O(h3)
将 K2 hy(xi a2h) hf (xi a2h, yi b21K1) 在x=xi处进行Taylor展开：
2020/3/18
12
K2 h f (xi , yi) a2hfx b21K1 f y O(h2)
h2 2
y(xi ) O(h3)
yi hf (xi , yi )
h2 2
f x(xi , yi ) f (xi , yi ) f y (xi , yi ) O(h3 )
类似地，若取前P+1项作为y(xi+1)的近似值，便得到
yi1
yi
hyi
h2 2!
yi
hP P!
yi(P)
P阶泰勒方法
例1. 使用高阶R-K方法计算初值问题
例如：取前两项可得到
y(xi1) y(xi ) hy(xi ) O(h2 )
y(xi ) hf (xi , y(xi )) O(h2 )
yi hf (xi , yi ) O(h2 )
2020/3/18
1
若取前三项，可得到截断误差为O(h3)的公式
y(xi1)
y(xi ) hy(xi )
2020/3/18
22
类似于前面的讨论，把K2、K3、K4分别在xi点 展成h的幂级数，代入线性组合式中，将得到的 公式与y(xi+1)在xi点上的泰勒展开式比较，使其 两式右端直到h4的系数相等，经过较复杂的解方 程过程便可得到关于ci，ai，bij的一组特解
a2=a3=b21=b32=1/2 b31=b41=b42=0 a4=b43=1 c1=c4=1/6 c2=c3=1/3
确定系数 c1、c2、a2、b21 ，可得到有2阶精度的算法格式
2020/3/18
11
因此
y(xi1) y(xi ) K
y(xi ) (c1K1 c2K2)
将y(xi+1)在x=xi处进行Taylor展开：
y(xi1)
y(xi ) hy(xi )
h2 2!
y(xi ) O(h3)
3
Runge-Kutta 方法是一种高精度的单步法,简称R-K法 9.4.1 龙格-库塔(R-K)法的基本思想
Euler公式可改写成
yi1 yi K K hf (xi , yi )
则 yi+1 的 表 达 式 与 y(xi+1) 的 Taylor 展 开 式 的 前 两 项 完全相同,即局部截断误差为O(h2)。
a2c2h2 fx b21c2h2 ff y O(h3)
2020/3/18
13
令 y(xi1) yi1 对应项的系数相等，得到
c1 c2 1 ,
a2c2
1 2
,
b21c2
1 2
这里有 4 个未知 数，3 个方程。
存在无穷多个解。所有满足上式的格式统称为2阶 龙格 - 库塔格式。
2020/3/18
yi
1
yi
1 6
(K1
4K2
K3)
K1
K
2
hf (xi , yi ) h
hf (xi 2
,
yi
1 2
K1 )
K3 hf (xi h, yi K1 2K 2 )
2020/3/18
20
9.4.4 四阶(经典)龙格—库塔法
如果需要再提高精度，用类似上述的处理 方法，只需在区间[xi,xi+1]上用四个点处的斜 率加权平均作为平均斜率K*的近似值，构成一 系列四阶龙格—库塔公式。具有四阶精度，即 局部截断误差是O(h5)。
a2c2
a3c3
1 2
,
a22c2
a32c3
1 3
,
b221c2
(b31
b32 )2 c3
1 3
a2b32c2
1, 6
b21c2
(b31
b32 )c3
1 2
a2b21c2
a3c3 (b31
b32 )
1 3
b21b32c3
1 6
2020/3/18
19
参数的选择不唯一，从而构成一类不同的三阶R-K 公式，下面给出一种常用的三阶R-K公式，形似 simpson公式：
h f (xi , yi ) a2hfx b21hff y O(h3)
K1 =hf(xi, yi)
yi1 yi (c1K1 c2K2)
y(xi ) c1hf (xi, yi )
c2h f (xi , yi ) a2hfx b21hff y O(h3)
y(xi ) (c1 c2)hf (xi, yi )
答案是否定的！无论四个参数怎样选择,都不能使公式 的局部截断误差提高到三阶。
这说明每一步计算两个函数值的二阶R-K方法最高阶为 二阶。若要获得更高阶得数值方法,就必须增加计算函数 值的次数。
2020/3/18
17
9.4.3 三阶龙格—库塔法
为进一步提高精度，在区间[xi, xi+1]上除两点xi和
2020/3/18
23
➢四阶(经典)Runge-Kutta方法
yi
1
yi
1 6
(K1
2K2
2K3
K4)
K1 hf (xi , yi )
K 2
hf
( xi
h 2
,
yi
1 2
K1 )
K 3
hf
( xi
h 2
,
yi
1 2
K2)
K4 hf (xi h, yi K3 )
2020/3/18
24
hf (xi, yi )
则上式化为
如下图
y
y y(x)
yi1 yi hf (xi , yi )
K
K
即Euler方法
xi
Euler方法也称为一阶Runge-Kutta方法
xi1
x
2020/3/18
10
9.4.2 二阶龙格—库塔法 在[xi, xi+1]上取两点xi和xi+a2= xi +a2h,以该两点处
y(xi1) y(xi ) y(i )(xi1 xi )
其中i (xi , xi1)
即
y(xi1) y(xi ) hy(i )
2020/3/18
8
引入记号 y(xi1) y(xi ) K
K hy(i) hf i, y(i)
yi1 yi K
K可以认为是y y(x)在区间[xi , xi1]上的平均斜率
推导过程与前面类似，由于过程复杂，这 里从略，只介绍最常用的一种四阶经典龙格— 库塔公式。
2020/3/18
21
设 yi+1=yi+c1K1+c2K2+c3K3+c4K4
这里K1、K2、K3、K4为四个不同点上的函数值， 分别设其为
K1=hf (xi, yi) K2=hf (xi+a2h, yi+b21K1) K3=hf (xi+a3h, yi+b31K1+b32K2) K4=hf (xi+a4h, yi+b41K1+b42K2+b43K3) 其中c1、c2、c3、c4、a2、a3、a4、b21、b31、 b32、b41、b42、b43均为待定系数。
y
只要使用适当的方法求 出y(x)在区 间[xi , xi1]上平均斜率的近似值 K
就可得到相应的Runge-Kutta方法
y y(x)
K
2020/3/18
xi
xi1
x
9
如果以y(x)在xi处的斜率作为y(x)在[xi , xi1]上的平均斜率K
即 K hy(xi ) hf [xi, y(xi )]
的斜率值K1和K2的加权平均(或称为线性组合)来求取 平均斜率k*的近似值K，即
K c1K1 c2K2
式中:K1为xi点处的切线斜率值 K1 =hf(xi, yi)=hy'(xi) K2为xi +a2h点处的切线斜率值,比照改进的欧拉
法,将xi+a2视为xi+1，即可得 K2 hf (xi a2h, yi b21K1)
2020/3/18
5
于是可考虑用函数f(x,y)在若干点上的函数值的 线性组合来构造近似公式，构造时要求近似公式在 (xi,yi)处的Taylor展开式与解y(x)在xi处的Taylor展开式 的前面几项重合，从而使近似公式达到所需要的阶数。 既避免求高阶导数,又提高了计算方法精度的阶数。或 者说,在[xi,xi+1]这一步内多计算几个点的斜率值，然 后将其进行加权平均作为平均斜率，则可构造出更高 精度的计算格式，这就是龙格—库塔（Runge-Kutta） 法的基本思想。
xi
1为h区间
2
x的i , x中i1点。也称中点公式。
Q:为获得更高的精度，应该如何进一步推广？
2020/3/18
16
二级R-K方法是显式单步式,每前进一步需要计算两个
函数值。由上面的讨论可知,适当选择四个参数c1,c2,a2,
b21，可使每步计算两次函数值的二阶R-K方法达到二阶 精度。能否在计算函数值次数不变的情况下,通过选择不 同的参数值,使得二阶R-K方法的精度再提高呢?
称为P阶龙格－库塔方法。
其中ai,bij,ci为待定参数，要求上式yi+1在点(xi,yi)处作 Tailor展开，通过相同项的系数确定参数。
2020/3/18
7
Runge-Kutta方法的推导思想
对于常微分方程的初值问题
y f (x, y) a x b
y(a)
y0
的解y=y(x)，在区间[xi, xi+1]上使用微分中值定理，有
§ 9.4 龙格－库塔方法
得到高精度方法的一个直接想法是利用Taylor展开
假设式
y' =f(x,y) (a≤x≤b)
中的 f(x,y) 充分光滑,将y(xi+1)在x i点作Taylor展开，若 取右端不同的有限项作为y(xi+1)的近似值，就可得到 计算y(xi+1)的各种不同截断误差的数值公式。
当p≥2时,要利用泰勒方法就需要计算f(x,y)的高阶 微商。这个计算量是很大的，尤其当f(x,y)较复杂 时，其高阶导数会很复杂。因此，利用泰勒公式 构造高阶公式是不实用的。但是泰勒级数展开法 的基本思想是许多数值方法的基础。
R-K方法不是直接使用Taylor级数,而是利用它的思想
2020/3/18
2020/3/18
6
一般龙格－库塔方法的形式为
yi1 yi c1K1 c2K2 cpK p K1 hf (xi, yi ) K2 hf (xi a2h, yi b21K1) • •• •• •• •• •• K p hf (xi aph, yi bp1K1 bp, p1K p1)
K2
y y(x)
K3
K1
K
K3 hf (xi a3h, yi b31K1 b32K2 )
2020/3/18
xi
xia2
xia3
x
18
同理推导二阶公式，将y(xi+1)和yi+1在x=xi处进 行Taylor展开，使局部截断误差达到O(h4)，使对
应项的系数相等，得到系数方程组：
c1 c2 c3 1
2020/3/18
4
同理，改进Euler公式可改写成
yi
1
yi
1 2
K1
1 2
K2
K1 hf ( xi , yi )
K2 hf ( xi h, yi K1)
局部截断误差为O(h3)
上述两组公式在形式上共同点:都是用f(x,y)在某 些点上值的线性组合得出y(xi+1)的近似值yi+1, 且增 加计算的次数f(x,y)的次数,可提高截断误差的阶。如 欧拉法:每步计算一次f(x,y)的值,为一阶方法。改进欧 拉法需计算两次f(x,y)的值，为二阶方法。
xi+a2= xi +a2h,以外，再增加一点xi+a3= xi +a3h ，用这 三点处的斜率值K1、K2和K3的加权平均得出平均斜率
K*的近似值K，这时计算格式具有形式:
y
yi1 yi c1K1 c2K2 c3K3
KK12
hf hf
(xi , ( xi
yi ) a2h,
yi
b21K1)
进的欧拉格式是众多的二阶龙格—库塔法中的一种
特殊格式。
2020/3/18
15
若取 a2 b21 塔公式。
1, 2
c1， 就0, 是c2另 一1 种形式的二阶龙格 - 库
yi1 yi K2
K1 hf (xi , yi )
K2
hf
( xi
1 2
h,
yi
1 2
K1)
此计算公式称为变形的二阶龙格—库塔法。式中