龙格库塔公式的推导

(3-5-22)
8
r=4时

同样可以得到常用的四阶RK公式：
h y n1 y n 6 (k1 2k 2 2k3 k 4 ) k f (t , y ) n n 1 h h k 2 f (t n , y n k1 ) 2 2 k f (t h , y h k ) n n 2 3 2 2 k 4 f (t n h, y n hk3 )
s 1 y () limsG( s)U ( s) lim s 1 2 2 s 0 s 0 ( s 1) s
18

由 G( z ) 可得：
z 1 y () lim G ( z )U ( z ) z 1 z z 1 K z ( z 1) Tz l im T 2 2 z 1 z ( z e ) ( z 1) K zT (1 e T ) 2
(3-5-23)
9
最优步长控制策略

由此作出判断： 1、若 en 0，则本步积分成功，现确 定下一步的最大步长 hn1 。 假定 hn 1 足够小，则 (t n hn1 ) (t n )， 下一步误差为： (tn＋1 ) hnk1 (tn ) hnk1 en 1 yn 1 1 yn 1

s G( s) 2 s 2s 1
( s 1)
（2）
p1 1, p2 1, q1 0, n 2, m 1
e T , p2 e T , q 1 p1
17
（3）
K z ( z 1) G( z ) T 2 (z e )
（4）对斜坡函数 u(t ) t , G(s) 具有非零且 有限的稳态值：
龙格库塔公式的推导

为了避免计算各阶导数和偏导数，将式 (3.5.12)写成 r y (t h) y (t ) h b k (3.5.13)

i 1
i i

其中r称为阶数，bi 待定系数，k i 由下式 决定 i 1 ki f (t ci h , y (t ) h a j k j ) (3.5.14)
1/ k
0 e n
1/ k
hn
(3-5-38)
12

2、若 en 0 ， 则本步失败，按式(3.5.38) 求出一个积分步长，它表示重新积分的 本步步长，再算一遍，即：
0 hn hn en
1/ k
(3-5-39)
13

（3－5－17）
3

将(3-5-17)式代入(3-5-16)式，再将式 (3.5.16)代入式(3.5.13)，得
y (t h) y (t ) h b1k1 b2 k 2 f (t , y ) f (t , y ) y (t ) h b1 f (t , y ) b2 f (t , y ) c2b2 h a1b2 k1h t y 2 f 2 f y (t ) (b1 b2 )hf (t , y ) b2c2 h a1b2 h f (t , y ) t y
(3-5-18)
4

将式 (3.5.18) 与式 (3.5.12) 逐项比较，按 照对应项系数相等比较可得：
b1 b2 1 b2 c2 1 / 2 a b 1 / 2 1 2
(3-5-19)
5


式(3.5.19)是一个不定方程，它有无穷多 个解。 取 a1 1/ 2, b1 0, b2 1, c2 1/ 2 ，可得
由于假定了 hn1 足够小，因此 (t n ) 基本不 变，故必须限制步长的缩小与放大，一 般限制的最大放缩系数为10，即要求：
0.1hn hn1 10hn
14



有关最优步长的控制，除此方法之外，还有吉 尔(Gear)法等。 采用最优步长控制后、 f 计算量有明显减少， 但上述两种控制方法对于函数中含有间断特性 的情况不适合。 因为在间断点附近会出现步长频繁放大、缩小 的振荡现象，由于最优步长控制法是以本步误 差外推下一步步长，因此振荡现象更为严重。
i 1, 2, 3,
j 1
r
c1 0
1


下面针对r的取值进行讨论： c1 0, k1 f (t , y) r=1，此时，有 因此有：
y(t h) y(t ) b1hf (t , y) （3－5－15）

取 b1 1 即得一阶RK公式，它就是欧拉 公式，因此可以说欧拉公式是RK公式的 特例。
试用欧拉法作为预估公式，梯形法作为校正公式， 写出5步仿真过程。
21
dy dt y t y 0 0 t0 0 h 0.5

s 2、已知 G s 2 s 3s 2
，
试用根匹配法建立该系统的频域模 型G z 。

Leabharlann Baidu
22
19
(1 e ) （5） 确定 K z T
T
2
0 （6）附加一个零点，为简单起见，令 q2 即零点匹配在z平面的原点。故有：
(1 e T ) 2 z ( z 1) G( z ) T ( z e T ) 2
附加这个零点不影响 k z
20
课堂测验

1、已知一微分方程模型模型为：
2
r＝2时


此时有： k1 f (t , y ) （3－5－16） k 2 f (t c2 h, y (t ) a1k1h) 将 f (t c2 h, y(t ) a1k1h) 在点 (t , y ) 展成泰勒 级数
f (t c2 h, y (t ) a1k1h) f (t , y ) c2 h f (t , y ) a1k1h f (t , y ) t y
15
根匹配法



关于最后一步，附加零点的说明： mn 因为一般的传递函数模型有： 因此在s平面上有n-m个零点在负无穷远 处。不妨假设均在 处，由此可见， T 0 在z平面上尚有n-m个零点在 e 处，即尚有n-m个零点在z平面的原点。
16
一个例子
设 ，试利用根匹配法， 求与之相匹配的离散化模型，步骤如下： （1） G(s) s 2 ，求得
y n1 y n hk 2 k1 f (t n , y n ) h h k 2 f (t n , y n k1 ) 2 2
(3-5-20)
6

取 a1 1, b1 b2 1/ 2, c2 1
，可得：
h y n1 y n 2 (k1 k 2 ) k1 f (t n , y n ) k f (t h, y hk ) n n 1 2
10

为使 en1 0 ，即：
(t n ) h
k n 1
yn 1

0
则有
0 ( y n 1) hn1 (t ) n
1/ k
11

将式(3.5.37)代入上式得:
hn1 0 h e n
k n

(3-5-21)
显然式(3.5.21)正好是前面介绍的改进 的欧拉公式。
7
r＝3时

按照前面的推导方法可以得到常用的三 阶RK公式 ：
h y n1 y n 4 (k1 3k 3 ) k f (t , y ) n n 1 k f (t h , y h k ) n n 1 2 3 3 2h 2h k 3 f (t n , y n k 2 ) 3 3