变步长的龙格库塔法
合集下载
82第二节 龙格—库塔法
h k y x n O hk 1 k! 2
k
(1)
h h k 若令 yn1 y xn hy xn y xn y xn (2) 2! k! 则 y xn1 yn1 O hk 1
y0 k1 2 k2 hf x0 h 2, y0 k1 2
y0 k3
k4 hf x0 h, y0 k3
y0 k2 2 k3 hf x0 h 2, y0 k2 2
k
x1 x0 h y1 y0 k
数学学院 信息与计算科学系
0.1832292
0.1584376
数学学院 信息与计算科学系
接上图
0.4 0.5 0.5 0.6 0.6 1.341667 1.416026 1.412676 1.482627 1.483281 0.0745394 0.0710094 0.0708400 0.0673253
0.1416245
数学学院 信息与计算科学系
数学学院 信息与计算科学系
由表8-4可见，虽然四阶龙格-库塔方法每步要 计算四次 f 的值，但以h=0.2为步长ቤተ መጻሕፍቲ ባይዱ计算结果就
有5 位有效数字，而欧拉法与预估计-校正方法以
h=0.1为步长的计算结果才具有2 位与3 位有效数字.
如果步长 h 也取0.2，则结果的精度会更低.
即公式(2)为k 阶方法.
数学学院 信息与计算科学系
二、龙格-库塔方法(R-K方法)
R-K方法不是通过求导数的方法构造近似公式, 而是通过计算不同点上的函数值, 并对这些函数值作 线性组合, 构造近似公式, 再把近似公式与解的泰勒 展开式进行比较, 使前面的若干项相同 , 从而使近似 公式达到一定的阶数.
k
(1)
h h k 若令 yn1 y xn hy xn y xn y xn (2) 2! k! 则 y xn1 yn1 O hk 1
y0 k1 2 k2 hf x0 h 2, y0 k1 2
y0 k3
k4 hf x0 h, y0 k3
y0 k2 2 k3 hf x0 h 2, y0 k2 2
k
x1 x0 h y1 y0 k
数学学院 信息与计算科学系
0.1832292
0.1584376
数学学院 信息与计算科学系
接上图
0.4 0.5 0.5 0.6 0.6 1.341667 1.416026 1.412676 1.482627 1.483281 0.0745394 0.0710094 0.0708400 0.0673253
0.1416245
数学学院 信息与计算科学系
数学学院 信息与计算科学系
由表8-4可见，虽然四阶龙格-库塔方法每步要 计算四次 f 的值，但以h=0.2为步长ቤተ መጻሕፍቲ ባይዱ计算结果就
有5 位有效数字，而欧拉法与预估计-校正方法以
h=0.1为步长的计算结果才具有2 位与3 位有效数字.
如果步长 h 也取0.2，则结果的精度会更低.
即公式(2)为k 阶方法.
数学学院 信息与计算科学系
二、龙格-库塔方法(R-K方法)
R-K方法不是通过求导数的方法构造近似公式, 而是通过计算不同点上的函数值, 并对这些函数值作 线性组合, 构造近似公式, 再把近似公式与解的泰勒 展开式进行比较, 使前面的若干项相同 , 从而使近似 公式达到一定的阶数.
第三部分龙格-库塔方法
内江师范学院数学与信息科学学院 吴开腾 制作
于是有
其中
y ( xn +1 ) − y ( xn ) = y '(ξ ), ξ ∈ ( xn , xn +1 ) h y ( xn +1 ) = y ( xn ) + hf (ξ , y (ξ ))
k * = f (ξ , y (ξ )) 称作区间 [ xn , xn +1 ] 上的平均斜率。 上的平均斜率 平均斜率。 问题：计算近似值y ( xn +1 ) 的关键是如何选择算法确定平均斜率 k *
（15）
f ( xn +1 , yn + h ( − k1 + 2 k 2 ))
内江师范学院数学与信息科学学院 吴开腾 制作
注释1 可以用Taylor展示证明格式（14） 注释1：可以用Taylor展示证明格式（14）具有三阶精 展示证明格式
度，并且还可以用类似的方法得到四阶及其以上的更高 阶精度的Runge-Kutta格式 阶精度的Runge-Kutta格式。 Runge 格式。
内江师范学院数学与信息科学学院 吴开腾 制作
h yn + （ k1 + 2 k 2 + 2 k3 +k 4） 6 f ( xn , y n ) h f ( x 1 , yn + k1 ) n+ 2 2 h f ( x 1 , yn + k 2 ) n+ 2 2 f ( xn +1 , yn + hk3 ) （16）
四阶龙格－ 四阶龙格－库塔格式计算结果
xn
0.1 0.2 0.3 0.4 0.5
yn
欧拉格式计算结果 xn yn y ( xn )
变步长的龙格库塔法
yn1
yn
h(1 2
K1
1 2
K2 )
K1
f (xn, yn )
K2 f (xn h, yn hK1)
一般地，RK方法设近似公式为
p
yn1 yn h ci Ki
K1
i 1
f ( xn , yn )
i 1
K
i
f ( xn ai h, yn h bij K j )
j 1
(i 2, 3
选择步长时要考虑两个问题： 1.如何衡量和检验计算结果的精度？ 2.如何根据所获得的精度处理步长？
HW: p.201 #6-8
结束语
谢谢大家聆听！！！
17
yi
2 2
hK2
1
2 2
hK3)
§2 Runge-Kutta Method
注：
龙格-库塔法的主要运算在于计算 Ki 的值，即计算 f 的
值。Butcher 于1965年给出了计算量与可达到的最高精 度阶数的关系：
每步须算Ki 的个数 2 3
4
5
6
可达到的最高精度 O (h2 ) O (h3 ) O(h4 ) O(h4 ) O(h5 )
7
O(h6 )
n8
O(hn2 )
由于龙格-库塔法的导出基于泰勒展开，故精度主要受
解函数的光滑性影响。对于光滑性不太好的解，最好 采用低阶算法而将步长h 取小。
➢ 变步长的Runge—Kutta Method
§2 Runge-Kutta Method
R n 1chp 1y(p 1) xn
Q: 由局部截断误差可以看出，步长 h 越小，局部截断 误差越小；但步长减小，在一定求解范围（区间）内 要完成的步数就增加了，步数增加会引起计算量增大， 导致舍入误差积累。因此要选取适当的步长。
《龙格-库塔方法》课件
优点
适用于求解一般和刚性常微分方程，并且具有较好的稳定性和收敛性。
龙格-库塔方法的原理
1
基本思路
通过将微分方程转化为差分方程，使用逼近公式计算数值解。
2
推导过程
使用泰勒级数展开和积分逼近，推导出不同阶次的龙格-库塔方法公式。
龙格-库塔方法的算法流程
1
步骤分解
包括选择步长、计算斜率、更新数值、迭代等一系列操作。
《龙格-库塔方法》PPT课 件
欢迎大家来到本次的PPT课件《龙格-库塔方法》。在本课程中，我们将一起 探索这一数值计算方法的原理、算法流程、应用领域以及注意事项。
什么是龙格-库塔方法？法。
特点
通过多步骤迭代计算，可以得到精确性较高的数值解。
2
具体实现过程
将算法流程转化为编程语言中的相应代码实现。
图解龙格-库塔方法
龙格-库塔方法示意图
图中展示了使用龙格-库塔方法进行迭代计算的示意 图。
对比其他数值方法的表现
通过与其他数值方法的比较，可以更好地理解龙格库塔方法的优势。
龙格-库塔方法的应用领域
常见应用范围
广泛应用于科学研究、工程计算、物理模拟等领域。
实际计算案例
通过实际案例，展示龙格-库塔方法在解决实际问题中的应用。
注意事项
1 龙格-库塔方法的收敛性和稳定性
在使用龙格-库塔方法时，需要注意数值解的收敛性和算法的稳定性。
2 常见误差和偏差的处理方法
介绍常见的误差和偏差来源，并提供相应的处理方法。
总结
龙格-库塔方法的优势和不足
总结龙格-库塔方法的优点和局限性，指出其适 用性和改进空间。
对未来数值计算的展望
展望未来龙格-库塔方法在数值计算领域的发展 和应用前景。
适用于求解一般和刚性常微分方程，并且具有较好的稳定性和收敛性。
龙格-库塔方法的原理
1
基本思路
通过将微分方程转化为差分方程，使用逼近公式计算数值解。
2
推导过程
使用泰勒级数展开和积分逼近，推导出不同阶次的龙格-库塔方法公式。
龙格-库塔方法的算法流程
1
步骤分解
包括选择步长、计算斜率、更新数值、迭代等一系列操作。
《龙格-库塔方法》PPT课 件
欢迎大家来到本次的PPT课件《龙格-库塔方法》。在本课程中，我们将一起 探索这一数值计算方法的原理、算法流程、应用领域以及注意事项。
什么是龙格-库塔方法？法。
特点
通过多步骤迭代计算，可以得到精确性较高的数值解。
2
具体实现过程
将算法流程转化为编程语言中的相应代码实现。
图解龙格-库塔方法
龙格-库塔方法示意图
图中展示了使用龙格-库塔方法进行迭代计算的示意 图。
对比其他数值方法的表现
通过与其他数值方法的比较，可以更好地理解龙格库塔方法的优势。
龙格-库塔方法的应用领域
常见应用范围
广泛应用于科学研究、工程计算、物理模拟等领域。
实际计算案例
通过实际案例，展示龙格-库塔方法在解决实际问题中的应用。
注意事项
1 龙格-库塔方法的收敛性和稳定性
在使用龙格-库塔方法时，需要注意数值解的收敛性和算法的稳定性。
2 常见误差和偏差的处理方法
介绍常见的误差和偏差来源，并提供相应的处理方法。
总结
龙格-库塔方法的优势和不足
总结龙格-库塔方法的优点和局限性，指出其适 用性和改进空间。
对未来数值计算的展望
展望未来龙格-库塔方法在数值计算领域的发展 和应用前景。
Runge-Kutta算法
Yn 1 Yn h(c1 K1 c2 K 2 c3 K 3 ) K F (t , Y ) 1 n n K 2 F (t n 2 h, Yn 21hK1 ) K 3 F (t n 3 h, Yn 31hK1 32 hK2 )
例
y (0) 1.
步长都取为 h 0.1 分别用以下两种系数：
1 1. a , 改进的Euler 法： 2 1 c1 c2 , 2 21 1. 2 积分公式： yn 1 y 0.1 y
n 2 n
yn 0.1 y
2 2 n
2
1 2. a , 3 2 1 3 c1 , c2 , 2 21 . 3 3 2 积分公式：
2 2 yn 1 yn 0.1 2 yn yn 3 0.1 yn 2
3
2
结果及比较
三阶显式Runge-Kutta方法
在推导二阶显式方法的过程中，注意到局部截断误差表达式中h3项 包含了以下表达式： Yn Ftt(t n ,Yn ) 2 FtY (t n ,Yn )Fn FYY (t n ,Yn )Fn2 FY(t n ,Yn )Ft n ,Yn ) FY(t n ,Yn )Fn (t 因此若要在局部截断误差中消去h3项，必须增加包含了以上各项的 多个方程，同时我们注意到r=2时，只有c1,2 , 1 , 21 等四个待定系数， 少于方程的数目，所以这样的系数不存在。故： r=2时Runge-Kutta 方法只能是二阶的。要得到三阶的方法，则必须有r=3。
类似前面的推导,可以导出更高阶的Runge Kutta公式. 关于Runge Kutta方法,有以下几点需要特别指出：
龙格-库塔方法-文档资料
3
c3a 2b32
c3a3 1
6
; 2
O (h4)
常见的2种三阶方法：
库塔三阶方法
h yn1yn6(k14k2k3)
k1
f(xn,yn);
k2
hh f(xn2,yn2k1)
k 3 f(x n h ,y n h k 1 2 h k 2 ) •5
四级方法：N = 4
局部截断误差 O ( h 5 )
可见误差随着 x n 的增加呈指数函数递减
当 f y 0 时，微分方程是不稳定的； 而 f y 0 时，微分方程是稳定的。
上面讨论的稳定性，与数值方法和方程中 f 有关
•21
实验方程： y y C ,R e () 0
D e f 3 对单步法 yn 1ynh(xn,yn,h )应用实验方程，
e n 1 e n h [ ( x n ,y ( x n ) , h ) ( x n ,y n , h ) ] T n 1
•15
因为单步法是 p 阶的：h0,0hh0满足|Tn1|Chp1
|e n 1| |e n| h L |e n| C h p 1|en |
其中 1hL,C hp1
•18
三、绝对稳定性 /*Absolute Stibility*/ 计算过程中产生的舍入误差对计算结果的影响
首先以Euler公式为例，来讨论一下舍入误差的传播:
yn1ynhf(xn,yn)
设实际计算得到的点 x n 的近似函数值为 yn yn n，
其中 y n 为精确值， n 为误差
yn1ynhf(xn,yn)
通过适当选取参数 1,2和 p 的值，使得公式具有 2阶精度！！
•3
由泰勒公式展开，要使公式具有 2 阶精度，只需
c3a 2b32
c3a3 1
6
; 2
O (h4)
常见的2种三阶方法：
库塔三阶方法
h yn1yn6(k14k2k3)
k1
f(xn,yn);
k2
hh f(xn2,yn2k1)
k 3 f(x n h ,y n h k 1 2 h k 2 ) •5
四级方法：N = 4
局部截断误差 O ( h 5 )
可见误差随着 x n 的增加呈指数函数递减
当 f y 0 时，微分方程是不稳定的； 而 f y 0 时，微分方程是稳定的。
上面讨论的稳定性，与数值方法和方程中 f 有关
•21
实验方程： y y C ,R e () 0
D e f 3 对单步法 yn 1ynh(xn,yn,h )应用实验方程，
e n 1 e n h [ ( x n ,y ( x n ) , h ) ( x n ,y n , h ) ] T n 1
•15
因为单步法是 p 阶的：h0,0hh0满足|Tn1|Chp1
|e n 1| |e n| h L |e n| C h p 1|en |
其中 1hL,C hp1
•18
三、绝对稳定性 /*Absolute Stibility*/ 计算过程中产生的舍入误差对计算结果的影响
首先以Euler公式为例，来讨论一下舍入误差的传播:
yn1ynhf(xn,yn)
设实际计算得到的点 x n 的近似函数值为 yn yn n，
其中 y n 为精确值， n 为误差
yn1ynhf(xn,yn)
通过适当选取参数 1,2和 p 的值，使得公式具有 2阶精度！！
•3
由泰勒公式展开，要使公式具有 2 阶精度，只需
龙格库塔法介绍
h 0.005稳定.
2) 改进欧拉法(预测 — 校正，即二阶R K法)：
yn
1
yn
h
k1 2
k2 2
,
k1
f (xn, yn ) yn,
k2 f (xn h, yn hk1) ( yn hyn )，
即yn1 故
yn
当初值准确即e0 0时，整体误差为en O(h p ).
证明
考察单步法的收敛性归结为验证增量函数(x, y,h)是否
满足Lipschitz条件.
欧拉法 : (x, y,h) f (x, y)，L L.
改进的欧拉法：
yn1
yn
h[ 2
f
(xn, yn )
f
( xn 1,
xn
f
( x,
y ( x))dx
r
h ci
i 1
f
( xn
ih,
y ( xn
ih)).
yn1 yn h(xn, yn, h),
(3.4)
其中
r
(xn, yn, h) ciki ,
(3.5)
i 1
k1 f (xn, yn )，
欧拉法r 1, p 1.改进 欧拉法r 2, p 2.
k1)
k3）
k3 f (xn h, yn hk1 2hk2 )
称为库塔三阶方法.
阶数p和段数r（计算函数值次数）的关系
r12 p1 2
3 4 5 6 7 r≥8 3 4 4 5 6 r-2
常用的经典四阶龙格 库塔方法:
第二节 龙格-库塔方法汇总.ppt
ki f ( xn aih, yn h bijk j );i 1, ..., N
.精品课件.j 1
11
(1)一级二阶的隐式中点方法：
yn1 yn hk1
k1
f ( xn
h, 2
yn
hk1 ) 2
(2)二级四阶的隐式R-K方法：
yn1
yn
h(k1 k2 ) 2
k1
1 f [xn ( 2
第二节 龙格-库塔法
基本思想：利用 f ( x, y) 在某些特殊点上的函数值的 线性组合来构造高阶单步法的平均斜率。
什么叫平均斜率？
对差商 y( xi1) y( xi ) 应用微分中值定理，有，
h
y( xi1) y( xi ) hy( xi h) 利用微分方程 y f ( x, y)，有
y( xn1 )
y(h) n1
Ch5
Step2:取
h 为步长，再从 2
x n 出发，经过两步计算得到
y( xn1)
y(h2) n1
2C( h)5 2
.精品课件.
13
y( xn1)
y(h/ 2) n1
y( xn1)
y(h) n1
1 16
记
|
y(h/ 2) n1
y(h) n1
|
16( y( xn1)
y( xi1) y( xi ) hf ( xi h, y( xi h))
这里的 f ( xi h, y( xi h)) 称为平均斜率。
.精品课件.
1
可将改进的欧拉格式改写成
yi 1
1 2
( K1
K2 ),
K1 yi hf ( xi , yi ),
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
31
hK
1
32
hK )
2
... ... K
m
f ( x h, y
i m
m1
hK
1
m2
hK ...
2
m m 1
hK
m 1
)
其中i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij ( i = 2, …, m; j = 1, …, i1 ) 均为待定系数，确 定这些系数的步骤与前面相似。
5 (h ／ 2) n 1
h y ( xn 1 ) y 2c , 2 (h ／ 2) y ( xn 1 ) yn 1 1 由上两式 . (h) y ( xn 1 ) yn 1 16
(h ／ 2) y ( xn 1 ) yn 1
1 (h ／ 2) ( h) [ yn yn 1 1 ]. 15
一般地，RK方法设近似公式为
p yn 1 yn h ci K i i 1 K1 f ( xn , yn ) i 1 K i f ( xn ai h, yn h bij K j ) (i 2, 3 , p) j 1 确定原则是使近似公式在( xn , yn )处的Taylor展开式与
HW: p.201 #6-8
变步长的龙格—库塔方法
以经典四阶龙格—库塔公式为例。从节点xn出发，以h为
(h) 步长求一近似值yn 1 , (h) 5 y ( xn 1 ) yn 1 ch ,
将步长折半，即取 值y
(h ／ 2) n 1
h 为步长从xn 跨两步到xn 1，求一近似 2
5
h , 每跨一步的截断误差是c ,因此有 2
2 1 2 2 2
来自hK1 1 hK 2 1
2 2 2 2
hK )
3
hK )
2
§2 Runge-Kutta Method
注：
龙格-库塔法的主要运算在于计算 Ki 的值，即计算 f 的
值。Butcher 于1965年给出了计算量与可达到的最高精 度阶数的关系：
Gill公式：4阶经典龙格-库塔公式的一种改进
yi 1 yi h 6 K1 2 2 K 2 2 2 K 3 K 4 K1 K2 K3 K4 f ( xi , yi ) h f ( xi h , y i 2 2 K1 ) f ( xi h 2 , yi f ( xi h , yi
公 式
Euler公式 改进Euler公式 yn 1 yn hK1 K1 f ( xn , yn ) 1 1 y y h ( K K2 ) n 1 n 1 2 2 K f (x , y ) 1 n n K 2 f ( xn h, yn hK1 )
（ 2 ）如果Δ<ε，反复将步长加倍，直到Δ>ε
为止，并以上一次步长的计算结果作为 y i 1 。
这种通过步长加倍或折半来处理步长的方法称为
变步长法。表面上看，为了选择步长，每一步都要 反复判断Δ，增加了计算工作量，但在方程的解y(x)
变化剧烈的情况下，总的计算工作量得到减少，结
果还是合算的。
高阶Runge—Kutta Method
§2 Runge-Kutta Method
最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta Method */ ： h
y i 1 yi 6 ( K 1 2 K 2 2 K 3 K 4 ) f ( xi , yi ) h f ( xi h , y K1 ) i 2 2 f ( xi h , yi h K2 ) 2 2 f ( x i h, yi hK 3 ) K1 K2 K3 K4
每步须算Ki 的个数
2
3
4
5
6
7
O(h6 )
n8
O(hn2 )
4 可达到的最高精度 O(h2 ) O(h3 ) O(h4 ) O(h ) O(h5 )
由于龙格-库塔法的导出基于泰勒展开，故精度主要受
解函数的光滑性影响。对于光滑性不太好的解，最好 采用低阶算法而将步长h 取小。
变步长的Runge—Kutta Method
y ( x )在xn处的Taylor展开式的前面项尽可能多地重合。
7.4.6 变步长的龙格-库塔法
在微分方程的数值解中，选择适当的步长是非常 重要的。单从每一步看，步长越小，截断误差就越 小；但随着步长的缩小，在一定的求解区间内所要 完成的步数就增加了。这样会引起计算量的增大， 并且会引起舍入误差的大量积累与传播。因此微分 方程数值解法也有选择步长的问题。 以经典的四阶龙格-库塔法(7.20)为例。从节点 (h) xi出发，先以h为步长求出一个近似值，记为 yi 1 ， 由于局部截断误差为 O(h 5 ) ，故有
y( xi 1 ) y
( h) i 1
ch
5
当h值不大时，式中的系数c可近似地看作为常数。
然后将步长折半,即以为 步长,从节点xi出发,跨 h ( ) 两步到节点xi+1,再求得一个近似值 yi 21 ,每跨一步的 h ,因此有 截断误差是 c
5
h 2
2
y ( xi 1 ) y ( x
h) yi( 1)
这表明以 y 作为 y( xi 1 ) 的近似值，其误差可用步 h 长折半前后两次计算结果的偏差 ( ) 来判断所选步长是否适当
(h) 2 yi y 1 i 1
当要求的数值精度为ε时：
（1）如果Δ>ε，反复将步长折半进行计算，直
y i 1 至Δ<ε为止,并取其最后一次步长的计算结果作为
§2 Runge-Kutta Method
Rn1 ch p1 y( p1) xn
Q: 由局部截断误差可以看出，步长 h 越小，局部截断 误差越小；但步长减小，在一定求解范围（区间）内 要完成的步数就增加了，步数增加会引起计算量增大， 导致舍入误差积累。因此要选取适当的步长。 选择步长时要考虑两个问题： 1.如何衡量和检验计算结果的精度？ 2.如何根据所获得的精度处理步长？
y
i 1
§2 Runge-Kutta Method
y h [ K K ... K ]
i 1 1 2 2 m m i i
K f (x , y )
1
K f ( x h, y
2 i 2 i
21
hK )
1
K f ( x h, y
3 i 3 i
h ) 2 i 1 (
这样
由此可得
2 y ( x i 1 ) y i 1 1 (h) 16 y ( x i 1 ) y i 1
h ) 2c 2 h ( )
5
y( xi 1 ) y
h ( ) 2 i 1
h ( ) 2 i 1
1 (y 15
h ( ) 2 i 1