5.2龙格库塔法

第五章 常微分方程的数值解法

5.2 龙格-库塔法

一、教学目标及基本要求

通过对本节课的学习，使学生掌握常微分方程、常微分方程方程组的数值解法。

二、教学内容及学时分配

本节课主要介绍常微分方程的数值解法。具体内容如下： 讲授内容：龙格库塔方法、收敛性与稳定性

三、教学重点难点

1．教学重点：龙格库塔方法、收敛性与稳定性。 2. 教学难点：收敛性与稳定性。

四、教学中应注意的问题

多媒体课堂教学为主。适当提问，加深学生对概念的理解。

五、正文

龙格－库塔方法

引言 龙格－库塔方法的基本思想

'11()()

()()()(,())n n n n y x y x y y x y x hf y h

ξξξ++-=⇒=+

令*(,())K f y ξξ=称为区间1[,]n n x x +上的平均斜率，只要对*K 提供一种算法，即可推导出一种计算公式。

欧拉公式只是取点n x 的斜率作为区间1[,]n n x x +的平均斜率*K ，精度自然很低。

考察改进的欧拉公式

11211

()22

n n y y h k k +=++1(,)n n k f x y =21(,)n n k f x h y hk =++

它利用n x ，1n x +两点的斜率取算术平均，1n x +处斜率通过已知信息n y 用欧拉

公式预报得到。

可以考虑设法在1[,]n n x x +多取几个点的斜率值，将它们加权平均作为区间

1[,]n n x x +的平均斜率*K 。这就是龙格－库塔方法的基本思想。

1、二阶龙格－库塔方法

考察1[,]n n x x +内一点,01n p n x x hp p +=+<≤，希望通过,n n p x x +两个点的斜率

12,K K 加权平均得到*K ，即令112((1))n n y y h K K λλ+=+-+

取1(,)n n K f x y =，如何预报n p x +处斜率2K ？

仿照改进的欧拉公式，先用欧拉公式预测()n p y x +的值n p y +：

1n p n y y phK +=+

然后用n p y +计算2(,)n p n p K f x y ++=，从而得

112121[(1)](,)(,)

n n n n n p n y y h K K K f x y K f x y phK λλ++=+-+==+

适当选取,p λ，使上述公式具有较高得精度。假定()n n y y x =，分别将12

,K K 泰勒展开：

'1(,)()n n n K f x y y x ==

212

'

''

2

(,)

(,)[(,)(,)(,)]()()()()

n p n n n x n n n n y n n n n K f x y pkK f x y ph f x y f x y f x y O h y x phy x O h +=+=+++=++

代入得

'2''31()()()()n n n n y y x hy x ph y x O h λ+=+++

按泰勒展开2'

''

31()()()()2

n n n n h y y x hy x y x O h +=+++ 比较得，只要1

2

p λ=

，公式截断误差为3()O h 特别，当1,1/2p λ==，就是改进的欧拉公式， 改取1,1/2p λ==，

12n n y y hk +=+，1(,)n n k f x y =，21(,)22

n n h h

k f x y k =++

称为变形的欧拉公式，亦称中点公式。 2、三阶龙格－库塔方法

为提高精度，设除n p x +外再考察一点,01n q n x x hq p q +=+<≤≤，用三个点

,,n n p n q x x x ++的斜率123,,K K K 加权平均得出*123(1)K K K K λμλμ=--++

1123((1))n n y y h K K K λμλμ+=+--++

为预测n q x +的斜率值3K ，将12,K K 加权平均得出[,]n n q x x +上的平均斜率，从而得

12((1))n q n y y qh K K αα+=+-+

然后用n q y +计算3(,)n q n q K f x y ++= 设计的计算公式具有如下形式

1123121312[(1)](,)(,)

(,((1)))

n n n n n p n n q n y y h K K K K f x y K f x y phK K f x y qh K K λμλμαα+++=+--++==+=+-+

运用泰勒展开选择参数,,,,p q λμα 下列库塔公式是其中一种

112312112

312[4]

6

(,)(,)

2(,(2))

n n n n n n n q n h

y y K K K K f x y h

K f x y K K f x y h K K +++=+++==+=+-+

3、四阶龙格－库塔方法

继续上述过程，得出四阶方法，下列是经典格式

11234121123122

413[22]

6

(,)(,)2(,)2

(,)

n n n n n n n n n n h

y y K K K K K f x y h K f x y K h K f x

y K K f x y hK ++

+

+=++++==+=+=+

流程图P105 例P105

1211121312222413332(,)2()

2(,)()22()

2

2()

2(,)()22()

22()

(,)()()

n n n n n

n n n n n n n n n n n n n n n x K f x y y y h x h h K f x y K y K h y K h x h h K f x y K y K h y K x h K f x y hK y hK y hK +++==-

+=+=+-

++=+=+-

++=+=+-

+ 4、变步长龙格－库塔方法

步长越小，截断误差越小，但步长太小，计算量大，舍入误差也大，需要数值解法选择步长。步长选择考虑的问题： 1）如何衡量和检验计算结果的精度； 2）如何依据判定的精度来处理步长。

以经典龙格－库塔方法为例，先以步长h 求得一个近似值记为1h

n y +，经典龙格－库塔方法的截断误差为5()O h ，故有5

11

()()h n n y x y CO h ++-≈ 步长折半/2

511

()2(())2

h n n h y x y CO ++-≈ /2

/2/2/211

11111111()11()()()()1615

h h h h h h n n n n n n n n h

n n y x y y x y y y y y y x y δ++++++++++-≈⇒-≈-⇒=-- 区分两种情况处理步长：

1）δε>，继续将步长折半，直至δε<为止，取步长折半后的新值作为结果；