龙格库塔方法ppt课件

Runge-Kuttua方法和matlab原 理
1
龙格－库塔法（Runge-Kutta） 数值分析中，龙格－库塔法（Runge-Kutta） 是用于模拟常微分方程的解的重要 的一类隐式或显式迭代法。这些技术由数学家卡 尔·龙格和马丁·威尔海姆·库塔于1900年左右发 明。
经典四阶龙格库塔法 龙格库塔法的家族中的一个成员如此常用，以至 于经常被称为“RK4”或者就是 “龙格库塔法”。

2K2

2K3

K4
)
n N 1; N
x1 x0, y1 y0
nN?
Y
输出x1,y1
结束 9
function ff=rk(yy,x0,y0,h,a,b) %yy为y的导函数,x0,y0,为初值,h 为步长,a,b为区间
c=(b-a)/h+1;i1=1; %c为迭代步数；i1为迭代步数累 加值
解 f (x, y) y 2x/ y

yn1 K1 f
yn
xn,
h 6
K
1
yn ,

4K
2

K
3
,


K
2

f

x
n

h, 2
yn

h 2
K
1
,

K
3

f
xn

h,
yn
hK1

2hK 2 .
12
xn Yn
|yn- R-K3 y(xn)|
5
误差分析：
四阶R-K方法的每一步需要计算四 次函数值f，可以证明其局部截断误差 为O(h5).
注意上述公式对于标量或者向 量函数(y可以是向量)都适用。
6
注意的问题
Runge-Kutta法的主要运算在于计算 Ki 的值，即计 f 的值。计算量与可达到的最高精度阶数的关系：
每步须算Ki 的 个数
0. 1.095
1.09544
19
0.0005
误差 0.45e-4
y(xn) 1.0954
0. 1.184
1.18322 0.17e-4
21
0.0009
1.1832
0. 1.266
1.26491 0.15e-4
32
0.0013
1.2649
0. 1.343
1.34165 0.48e-4
44
0.0018
2
3
4
5
6
7 n8
可达到的最高 精度
O(h2) O(h3)
O(h4 )
O(h4 )
O(h5 )
Fra Baidu bibliotek
O(h6 )
O(hn2 )
R-K(高阶)方法不唯一,选择不同的 参数能得到
R不-K同方的法R的-K推公导式是基于Taylor展开 法，因而要求
解具有较好的光滑性，如果光滑性 较差精度可
能不如改进Euler方法,最好采用低 7
z(i1,1)=x;z(i1,2)=k1;z(i1,3)=k2;z(i1,4)=k3
;z(i1,5)=k4;z(i1,6)=y;
i1=i1+1;
end
end
11
例题4
例 求(3阶、4阶R K公式)解初值问题(步长h 0.1)
4

y

y

2x y
(0 x 1)
y(0) 1
1.3416
0. 1.416
1.41422 0.25e-4
54
0.0022
1.4142
0. 1.486
1.48326 0.55e-4
60
0.0028
1.4832
0. 1.552
1.54921 0.14e-4 1.5492
13
75
0.0033
改进Euler法一步需要计算两个函数值 (h=0.1) 四阶Runge-Kutta方法一步需要计算 四个函数值（h=0.2） 总计算量大致相当，但四阶RungeKutta方法精度更高
y2 n1

y(h) n1

17
• 可以通过检查步长折半前后两次计
算结果的偏差
=
(h)
y2 n1

y(h) n1
变步长方
来判断所选取的步长是否合法适
（1）对于给定的精度，如果 ,则反复将步长折半 进行计算直至 为止，这时取折半以后的“ 新值” 作为结果；
y=y0;z=zeros(c,6);
%z生成c行，6列的零矩阵存放结
果；
10
不断迭代运算： for x=a:h:b
if i1<=c k1=feval(yy,x,y); k2=feval(yy,x+h/2,y+(h*k1)/2); k3=feval(yy,x+h/2,y+(h*k2)/2); k4=feval(yy,x+h,y+h*k3); y=y+(h/6)*(k1+2*k2+2*k3+k4);
每跨一步截断误差是C
(
h 2
)5
,
有y(xn+1)-y(nh21)

2C
(
h)5 2
16
步长折半后，误差大约减少为原来的 1 ，即有 16
(h)
y( xn1) y( xn1)
y2 n1
y(h) n1

1 16
事后估计式为
y(
xn1 )

(h)
y2 n1

1 15

(h)
2
四阶RungeKutta方法
3
4
这样，下一个值(yn+1)由现在的值(yn)加上时间间隔 (h)和一个估算的斜率的乘积 决定。该斜率是以下斜率的加权平均： k1是时间段开始时的斜率； k2是时间段中点的斜率，通过欧拉法采用斜率 k1 来 决定 y在点 tn + h/2的值； k3也是中点的斜率，但是这次采用斜率 k2决定 y值； k4是时间段终点的斜率，其 y值用 k3 决定。 当四个斜率取平均时，中点的斜率有更大的权值：
四阶Runge-Kutta方法的MATL
8
四阶R-K方法实现
开始
输入 x0, y0,h, N
x1 x0 h;
K1

f
( x0 ,
y0 ), K2

f
( x0

h 2
,
y0

h 2
K1)
K3
f
( x0

h, 2
y0

h 2
K2 ),
f (x0
h,
y0
hK3)
y1

y0

h 6
(K1
实施方案
• 以经典四阶Runge-Kutta方法为例
设从节点xn出发，先以h为步长求出一个近似值，记为yn(h)1 因为经典格式的局部截断误差为O(h5),因此有
y(xn+1)-y(nh)1 Ch5 其中C与y(5)(x)在[xn, xn1]内的值有关
将步长折半，即取
h 2
为步长从xn跨两步到xn+1，求得一个近似值yn( h21)
14
五、变步长Runge-Kutta 方法 • 从每一步看，步长越小，截断误差
越小；但随着步长的缩小，在一定 求解范围内所要完成的步数就会增 加，步数的增加不但引起计算量的 增大，而且可能导致舍入误差的严 重积累，因此需要选择步长
•如何衡量和检验计算结果的精度 •如何依据所判定的精度来处理步长
15