变步长的龙格库塔法
- 格式：ppt
- 大小：542.50 KB
- 文档页数：8
下载文档原格式
合集下载
龙格库塔积分算法
龙格库塔法龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。
这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。
由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。
龙格库塔法是一种在工程上应用广泛的高精度单步算法，可以应用在物理、工程、控制、动力学中，如模糊控制、弹道分析以及分析光纤特性等，在系统仿真中得到广泛应用。
龙格库塔法源自于相应的泰勒级数方法，在每一插值节点用泰勒级数展开,其截断误差阶数也是，根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。
其中4 阶龙格库塔法是最常用的一种方法。
因为它相当精确、稳定、容易编程。
在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。
如果需要较高的精度, 可采取减小步长的方法即可。
4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。
1、初值问题对于一阶常微分方程的初值问题根据常微分方程的理论可知，此初值问题的解在区间[a，b]上存在，且唯一。
2、离散化取步长h=(b-a)/n，将区间[a , b]分成n个子区间：a=<=b在其中任意两点的曲线段上，根据积分中值定理，一段光滑曲线上至少有一点，它的斜率与整段曲线的平均斜率相同，得=y’() (0<<1)其中，=可以将上式改写成y()=y()+h*K (2.1)其中K为平均斜率，K=f()公式（2.1）表明，如果能够确定平均斜率K，就可以根据（2.1）式得到y()的值。
欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K，从而求得y()的近似值。
3、Euler法欧拉法虽然精度低，但它是最简单的一种显式单步法，也是龙格库塔法的基础。
首先，令、为y() 及y()的近似值，并且令平均斜率K=f()，即以点的斜率作为平均斜率K，便得到欧拉公式=+h* f() （3.1）4、改进的欧拉法此种方法是取、两点的斜率的平均值作为平均斜率K，即K= ,其中、均为y()以及y()的近似值，就得到改进后的欧拉公式（4.1）其中、分别为、两点的斜率值，即= ,=在上面的（4.1）式中，k2是未知的，采用一种叫预报法的方法来求解。
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 )
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
常用的经典四阶龙格 库塔方法:
变步长四阶龙格库塔方法在水文水力计算中的应用
变步长四阶龙格库塔方法在水文水力计算中的应用
梁建;李宝春;祁涛;卢俊
【期刊名称】《水资源与水工程学报》
【年(卷),期】2015(0)3
【摘要】将变步长四阶龙格库塔方法用于水库调洪演算与水面线的计算中。
变步长龙格库塔方法,可以有效地避免峰值的错过。
以安徽省小水库除险加固工程为例,对该水库采用传统试算方法和变步长龙格库塔方法进行调洪演算和水库溢洪道水面线计算,分别得到了下泄洪水过程线、水位变化过程曲线和溢洪道水面线。
结果表明:本方法易于程序设计,精度满足工程要求,可以用于水文水力计算中微分方程的求解。
【总页数】5页(P161-164)
【关键词】水库调洪;水面线计算;水文水力计算;变步长龙格库塔
【作者】梁建;李宝春;祁涛;卢俊
【作者单位】安徽省.水利部淮委水利科学研究院
【正文语种】中文
【中图分类】TV131
【相关文献】
1.变步长龙格库塔方法在安徽中小水库调洪中的应用 [J], 梁建;祁涛;卢俊
2.四阶龙格—库塔法在电机数字仿真中的适用范围与最佳步长的确定 [J], 邬学义
3.四阶龙格-库塔方法的程序设计与应用 [J], 罗丽珍; 吴庆军
4.四阶龙格库塔法在高士沟水库调洪计算中的应用 [J], 陈猛
5.溢洪道泄槽水面线计算的四阶龙格—库塔方法 [J], 范如琴
因版权原因，仅展示原文概要，查看原文内容请购买。
第二节 龙格-库塔方法
h( k1 + k2 ) yn +1 = yn + 2 hk1 1 1 3 3 k1 = f [ xn + ( + )h, yn + )hk2 ] +( + 2 6 4 4 6 hk2 1 3 1 3 k2 = f [ xn + ( − )h, yn + ( − )hk1 + ] 2 6 4 6 4
1 库塔三阶方法 库塔三阶方法
h h k1 = f ( xn , yn ); k2 = f ( xn + , yn + k1 ) 2 2
h yn +1 = yn + ( k1 + 4k2 + k3 ) 6
k3 = f ( xn + h, yn − hk1 + 2hk2 )
四级方法： 四级方法：N = 4 常见的2种四阶方法： 常见的 种四阶方法： 种四阶方法 经典龙格 库塔方法 龙格1经典龙格-库塔方法
1 中点法（修正的Euler法） 取 c1 = 0, c2 = 1, a2 = 1 中点法（修正的 法 2 h h
y n +1 = y n + hf ( x n +
常见的3种二级方法： 常见的 种二级方法： 种二级方法
二阶龙格库塔方法 2 二阶龙格库塔方法
h yn +1 = yn + [ f ( xn , yn ) + f ( xn + h, yn + hf ( xn , yn ))] 2
1 取 c1 = c2 = , a2 = 1 2
2
, yn +
2
f ( x n , y n ))
类似于N 的推导方法， 的推导方法 三级方法： 三级方法：N = 3 类似于 = 2的推导方法，可得到 1 c1 + c 2 + c 3 = 1; c 2 a 2 + c 3 a 3 = ; 2 O ( h4 ) 1 1 2 2 c 2 a 2 + c 3 a 3 = ; c 3 a 2 b32 = 6 3 常见的2种三阶方法： 常见的 种三阶方法： 种三阶方法
变步长 ode23tb 算法
变步长ode23tb算法是一种用于解非刚性常微分方程初值问题的数值算法。
该算法由MathWorks公司的工程师开发，其特点是能够自适应调整步长以提高求解精度和效率。
ode23tb算法是ode23算法的改进版本，通过引入更高阶的解析技术和自适应的步长控制，使得求解过程更加可靠和高效。
二、算法原理ode23tb算法基于龙格-库塔(Runge-Kutta)方法，通过描述性的一阶创建步骤和更高阶的增补步骤来近似解微分方程。
算法会根据当前的解的精度和误差情况，自适应地调整步长以确保求解的精度和效率。
ode23tb算法还引入了自动时步的策略，以处理不同变化速度的微分方程。
三、算法特点1. 自适应调整步长：ode23tb算法能够根据当前解的情况实时调整步长，避免不必要的计算和提高求解效率。
2. 高阶解析技术：通过引入更高阶的解析技术，ode23tb算法能够更准确地近似微分方程的解，在不增加过多计算成本的情况下提高求解精度。
3. 自动时步策略：算法能够根据微分方程的变化速度自动调整时步，适应不同变化速度的微分方程的求解。
作为一种高效可靠的数值求解算法，ode23tb算法在科学计算、工程领域和其他需要求解微分方程的问题中得到了广泛的应用。
例如在控制系统的建模和仿真中，ode23tb算法能够准确地求解系统的状态方程，为系统性能分析提供有力的支持。
在生物医学领域，ode23tb算法被用于建立生物模型和仿真生物系统的动态行为。
ode23tb算法还被应用于气象学、天文学、材料科学等领域，为复杂的微分方程问题提供了强大的求解能力。
五、算法发展随着科学技术的发展和对求解微分方程精度和效率要求的不断提高，ode23tb算法也在不断地进行改进和优化。
未来的发展方向主要包括进一步提高算法的自适应调整能力，完善算法的并行计算和局部精度控制，以满足各种复杂微分方程求解问题的需求。
六、总结ode23tb算法作为一种自适应调整步长的高效数值求解算法，在科学计算和工程实践中发挥着重要作用。
(完整版)第二节龙格-库塔方法
en1 en h[( xn , y( xn ), h) ( xn , yn , h)] Tn1
因为单步法是 p阶的：h0 , 0 h h0 满足| Tn1 | Chp1
| en1 || en | hL | en | Ch p1 | en |
其中 1 hL, Ch p1
en O(hp )
即取
K * 1K1 2 K2 yi1 yi h(1K1 2 K2 )
其中1 和 2 是待定常数。若取 K1 f ( xi , yi ) ，则
问题在于如何确定 xi p 处的斜率 K2 和常数 1 和 2 。
仿照改进的欧拉方法，用欧拉方法预测 y( xi p ) 的值，
yi p yi phK1
k2
f ( xn
h 2
,
yn
h 2 k1)
hh k3 f ( xn 2 , yn 2 k2 )
k4 f ( xn h, yn hk3 )
例2：用经典的龙格-库塔方法
求解下列初值问题 h 0.1
dy 。 dx
y
2x y
x (0,1)
解：经典的四阶龙格-库塔公式： y(0) 1
E(h) 1 h
绝对稳定域： 1 h 1
当 R 时， 1 h 1 2 h 0
绝对稳定区间：(2, 0)
❖经典的R-K公式：yn1
yn
h 6
(k1
2k2
2k3
k4 )
k1 f ( xn , yn ) yn
k2
f
(
yn1 yn hf ( xn , yn )
n1 yn1 yn1
n1 n h[ f ( xn , yn ) f ( xn , yn )] [1 hf y ( xn , )]n
因为单步法是 p阶的：h0 , 0 h h0 满足| Tn1 | Chp1
| en1 || en | hL | en | Ch p1 | en |
其中 1 hL, Ch p1
en O(hp )
即取
K * 1K1 2 K2 yi1 yi h(1K1 2 K2 )
其中1 和 2 是待定常数。若取 K1 f ( xi , yi ) ，则
问题在于如何确定 xi p 处的斜率 K2 和常数 1 和 2 。
仿照改进的欧拉方法，用欧拉方法预测 y( xi p ) 的值，
yi p yi phK1
k2
f ( xn
h 2
,
yn
h 2 k1)
hh k3 f ( xn 2 , yn 2 k2 )
k4 f ( xn h, yn hk3 )
例2：用经典的龙格-库塔方法
求解下列初值问题 h 0.1
dy 。 dx
y
2x y
x (0,1)
解：经典的四阶龙格-库塔公式： y(0) 1
E(h) 1 h
绝对稳定域： 1 h 1
当 R 时， 1 h 1 2 h 0
绝对稳定区间：(2, 0)
❖经典的R-K公式：yn1
yn
h 6
(k1
2k2
2k3
k4 )
k1 f ( xn , yn ) yn
k2
f
(
yn1 yn hf ( xn , yn )
n1 yn1 yn1
n1 n h[ f ( xn , yn ) f ( xn , yn )] [1 hf y ( xn , )]n
龙格库塔法的基本思想
龙格库塔法的基本思想
龙格库塔法，又称拉格朗日-龙格-库塔(Rallgendre-Kutta)法，是由德国数学家拉格朗日于1867年和德国数学家库塔于1901年分别提出的数值解法，是一种用于解决初值问题的数值分析方法，是飞行器运动动力学中用于近似解微分方程的方法之一。
龙格库塔法的基本思想是：将积分分成多个小步，将每一步步长求解结果串联起来，从而求出整个时间轴上每一个点的解。
为了满足既定的精度要求，需要不断减小步长，从而增加划分的步数。
如果步长继续减小，微分方程求解速度会变慢很多，故使用龙格库塔方法可较快求解方程。
龙格库塔法的实际应用中，通常是把初值问题用一系列积分步骤替代，每个积分步骤用四阶龙格库塔公式来模拟，每一步的计算结果都有效，最后得到的结果都是基于此，因此更加稳定。
此外，每个积分步骤可以由四步单步计算得出，因此实现和追踪也很容易。
总而言之，龙格库塔法比其他方法计算更快、更准确，可以用于解决复杂的初值问题，近年来大受工业界的欢迎，形成了一种重要的数值解法。
变步长四阶龙格库塔法原理
具体地说，将区分以下两种情况处理：
4
1. 对于给定的精度 ，如果 ，反复将步长折半 进行计算，直至 为止.
h ) 2 n 1 (
这时取最终得到的 y
作为结果；
2. 如果 ，反复将步长加倍，直到 为止， 这时再将步长折半一次，就得到所要的结果. 这种通过加倍或折半处理步长的方法称为变步长方法.
n 1
2
因此有
y ( xn 1 ) y
h ) 2 n 1 (
h 2c , 2
5
（3.15）
比较（3.14）式和（3.15）式我们看到，步长折半后， 1 误差大约减少到 ， 16
3
即有
y ( xn 1 ) y y ( xn 1 ) y
h ) 2 n 1 (h) n 1 (
表面上看，为了选择步长，每一步的计算量增加了，
但总体考虑往往是合算的.
事后估计式
y ( xn 1 ) y
h ) 2 n 1 ( ( ) 1 ( 2 [ yn 1 ynh ) ]. 1 15 h
这样，可以通过检查步长，折半前后两次计算结果的偏差
y
h ) 2 n 1 ( ( ynh ) 1
来判定所选的步长是否合适.
1
2° 如何依据所获得的精度处理步长？ 考察经典的四阶龙格-库塔公式
h yn 1 yn 6 ( K1 2 K 2 2 K 3 K 4 ), K1 f ( xn , yn ), h h K 2 f ( xn , yn K1 ), 2 2 K f ( x h , y h K ), n n 2 3 2 2 K 4 f ( xn h, yn hK 3 ).
第3章03龙格-库塔方法
提供 y(xn p ) 的预报值 yn p ，有 yn p yn phk1 ，因此 k2 f (xn p , yn p ) f (xn ph, yn phk1)
这样有计算公式：
yn1 yn h[(1 )k1 k2 )
k1 f (xn , yn ), k2 f (xn ph, yn phk1)
y(h) n1
1 16
y( xn1)
(h)
y2 n1
1 15
(h)
y2 n1
y(h) n1
误差事后估计
可以通过检查步长折半前后两次计算结果的偏差
y y ( h ) 2
(h)
n1
n1
来判断所选取的步长是否合适.
1) 对于给定的精度ε,如果δ> ε,则反复将步长折半 进行计算,直到δ< ε为止,这时取步长折半后的 新值作为结果;
k4
k3 f ( xn 2h, yn 21hk1 22hk2 ) f ( xn 3h, yn 31hk1 32hk2 33hk3)
四阶龙格－库塔公式每步要四次计算函数值，具有四阶精 度，局部截断误差是O(h5) .
四阶龙格－库塔格式
下列经典格式是其中常用的一种:
yn1 k2
yk
1.000000 1.005000 1.019025 1.041218 1.070800 1.107076 1.149404 1.197210 1.249975 1.307228 1.368514
y( xk ) yk
0.0 1.6×10－4 2.9×10－4 4.0×10－4 4.8×10－4 5.5×10－4 5.9×10－4 6.2×10－4 6.5×10－4 6.6×10－4 6.6×10－4
调用 2 次函数。有 2 阶精度。
这样有计算公式：
yn1 yn h[(1 )k1 k2 )
k1 f (xn , yn ), k2 f (xn ph, yn phk1)
y(h) n1
1 16
y( xn1)
(h)
y2 n1
1 15
(h)
y2 n1
y(h) n1
误差事后估计
可以通过检查步长折半前后两次计算结果的偏差
y y ( h ) 2
(h)
n1
n1
来判断所选取的步长是否合适.
1) 对于给定的精度ε,如果δ> ε,则反复将步长折半 进行计算,直到δ< ε为止,这时取步长折半后的 新值作为结果;
k4
k3 f ( xn 2h, yn 21hk1 22hk2 ) f ( xn 3h, yn 31hk1 32hk2 33hk3)
四阶龙格－库塔公式每步要四次计算函数值，具有四阶精 度，局部截断误差是O(h5) .
四阶龙格－库塔格式
下列经典格式是其中常用的一种:
yn1 k2
yk
1.000000 1.005000 1.019025 1.041218 1.070800 1.107076 1.149404 1.197210 1.249975 1.307228 1.368514
y( xk ) yk
0.0 1.6×10－4 2.9×10－4 4.0×10－4 4.8×10－4 5.5×10－4 5.9×10－4 6.2×10－4 6.5×10－4 6.6×10－4 6.6×10－4
调用 2 次函数。有 2 阶精度。
matlab ppt_13_龙格-库塔法
第十三讲 1. 基本思想
龙格――库塔法
1.1 数值法解常微分方程 已知 y = f (x, y ), (a ≤ x ≤ b) y (x 0 ) = y 0 取步长为h = (b − a)/n，将区间分成n个子区间， a = x0 < x 1 · · · < x i < · · · < x n = b 求离散点上的y (xi), y (xi+1) − y (xi) 由中值定理， = y (xi + θh) h 得 y (xi+1) = y (xi) + hkave 平均斜率为 kave = f (xi + θh, y (xi + θh)) 。 1.2 欧拉公式 取kave = f (xi, yi)，得欧拉公式 yi+1 = yi + hf (xi, yi) 取kave = (K1 + K2)/2，则得改进的欧拉公式 h yi+1 = yi + (K1 + K2) 2 其中K1和K2分别为xi 和xi+1 两点的斜率值。 其中K2 的计算为：先用欧拉法得到y (xi+1)的预报值 y i+1 = yi + hK1 = yi + hf (xi, yi) 再用预报值计算xi+1处的斜率值 K2 = f (xi+1, y i+1) = f (xi+1, yi + hf (xi, yi)) 可得 h yi+1 = yi + (K1 + K2) 2 h = yi + [f (xi, yi) + f (xi+1, yi + hf (xi, yi))] 2 (0 < θ < 1)
龙格――库塔法
1.1 数值法解常微分方程 已知 y = f (x, y ), (a ≤ x ≤ b) y (x 0 ) = y 0 取步长为h = (b − a)/n，将区间分成n个子区间， a = x0 < x 1 · · · < x i < · · · < x n = b 求离散点上的y (xi), y (xi+1) − y (xi) 由中值定理， = y (xi + θh) h 得 y (xi+1) = y (xi) + hkave 平均斜率为 kave = f (xi + θh, y (xi + θh)) 。 1.2 欧拉公式 取kave = f (xi, yi)，得欧拉公式 yi+1 = yi + hf (xi, yi) 取kave = (K1 + K2)/2，则得改进的欧拉公式 h yi+1 = yi + (K1 + K2) 2 其中K1和K2分别为xi 和xi+1 两点的斜率值。 其中K2 的计算为：先用欧拉法得到y (xi+1)的预报值 y i+1 = yi + hK1 = yi + hf (xi, yi) 再用预报值计算xi+1处的斜率值 K2 = f (xi+1, y i+1) = f (xi+1, yi + hf (xi, yi)) 可得 h yi+1 = yi + (K1 + K2) 2 h = yi + [f (xi, yi) + f (xi+1, yi + hf (xi, yi))] 2 (0 < θ < 1)
4.2 龙格库塔法
重复积分过程，直到相邻两次的积分值R2n 与Rn 满足 下列关系：
| R2n Rn | | R2n | 1
R2 n Rn | | R2 n
其中 为允许的误差限。
| R2n | 1
例1：利用龙贝格法计算单摆的振荡周期 T l /2 d T 4 1 g 0 2 0 [1 sin ( )sin 2 ]2 2 g 9.8m/s2 ， 其中 l 5m ， 0 20/ ， 105 program main external f call romb(0.,1.5708,1e-5,r2,f) write(*,"(3x,'B=',e12.6)") r2 End function f(x) f=4*sqrt(5/9.8)/(1-sin(10/3.14159)**2*sin(x)**2)**(1/2) end
subroutine romb(a,b,ep,r2,f) h=b-a t1=h/2.*(f(a)+f(b)) n=1 5 s=0.0 do k=0,n-1 s=s+f(a+(k+0.5)*h) end do t2=t1/2.+h/2.*s s2=t2+(t2-t1)/3. if(n/=1)goto 20 15 n=n+n h=h/2 t1=t2 s1=s2 goto 5 20 c2=s2+(s2-s1)/15. if(n/=2)goto 40
共有 n 1个分点 xk a kh, h b a , k 0,1,2,, n 用Tn 表示简单梯形法所求得的积分值。 1、考察一个子段 [ xk , xk 1 ], 其中点为 xk 1
2
n
1 ( xk xk 1 ), 2
变步长的龙格库塔法
然后将步长折半,即以为 h
2
步长,从节点xi出发,跨
两步到节点xi+1,再求得一个近似值
y
h ()
2 i1
,每跨一步的
截断误差是 c h 5 ,因此有
y2 ( xi1)y(xi( h 21))( h )2ch 25
这样
y(xi1) yi21 1
y(xi1)
y(h) i1
16
由此可得
(h)
y(xi1)yi 21
并且会引起舍入误差的大量积累与传播。因此微分
方程数值解法也有选择步长的问题。
以经典的四阶龙格-库塔法(7.20)为例。从节点
xi出发，先以h为步长求出一个近似值，记为
y (h) i1
，
由于局部截断误差为 O(h5 ) ，故有
y(xi1)yi( h1 ) ch 5
当h值不大时，式中的系数c可近似地看作为常数。
m
i
m
m1
1
m2
2
m m 1
m 1
其中i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij
( i = 2, …, m; j = 1, …, i1 ) 均为待定系数，确 定这些系数的步骤与前面相似。
§2 Runge-Kutta Method
➢ 最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta
（2）如果Δ<ε，反复将步长加倍，直到Δ>ε
为止，并以上一次步长的计算结果y i作1 为
。
这种通过步长加倍或折半来处理步长的方法称为
变步长法。表面上看，为了选择步长，每一步都要
反复判断Δ，增加了计算工作量，但在方程的解y(x)
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
变步长的龙格—库塔方法
以经典四阶龙格—库塔公式为例。从节点xn出发，以h为
步长求一近似值yn(h1) ,
y( xn1 )
y(h) n 1
ch5 ,
将步长折半，即取
h 2
为步长从xn
跨两步到xn
，求一近似
1
值y ( h／2) n 1
,
每跨一步的截断误差是c
h 2
5
,因此有
y( xn1 )
y ( h／2) n 1
2c
h 2
5
,
由上两式
y( xn1 )
y ( h／2) n 1
1.
y( xn1 )
y(h) n 1
16
y( xn1 )
y ( h／2) n 1
1 15
[
y ( h／2) n 1
y(h) n 1
].
公式
Euler公式 改进Euler公式
yn1 K1
yn f ( xn ,
hK1 yn )
完成的步数就增加了。这样会引起计算量的增大，
并且会引起舍入误差的大量积累与传播。因此微分
方程数值解法也有选择步长的问题。
以经典的四阶龙格-库塔法(7.20)为例。从节点
xi出发，先以h为步长求出一个近似值，记为
y (h) i 1
，
由于局部截断误差为 O(h5 ) ，故有
y( xi1 )
y (h) i 1
Q: 由局部截断误差可以看出，步长 h 越小，局部截断 误差越小；但步长减小，在一定求解范围（区间）内 要完成的步数就增加了，步数增加会引起计算量增大， 导致舍入误差积累。因此要选取适当的步长。
选择步长时要考虑两个问题： 1.如何衡量和检验计算结果的精度？ 2.如何根据所获得的精度处理步长？
HW: p.201 #6-8
（1）如果Δ>ε，反复将步长折半进行计算，直
至Δ<ε为止,并取其最后一次步长的计算结果作为y i 1
（2）如果Δ<ε，反复将步长加倍，直到Δ>ε
为止，并以上一次步长的计算结果yi作1 为
。
这种通过步长加倍或折半来处理步长的方法称为
变步长法。表面上看，为了选择步长，每一步都要
反复判断Δ，增加了计算工作量，但在方程的解y(x)
j 1
(i 2, 3L , p)
确定原则是使近似公式在( xn , yn )处的Taylor展开式与
y( x)在xn处的Taylor展开式的前面项尽可能多地重合。
7.4.6 变步长的龙格-库塔法
在微分方程的数值解中，选择适当的步长是非常
重要的。单从每一步看，步长越小，截断误差就越
小；但随着步长的缩小，在一定的求解区间内所要
7
O(h6 )
n8
O(hn2 )
由于龙格-库塔法的导出基于泰勒展开，故精度主要受
解函数的光滑性影响。对于光滑性不太好的解，最好 采用低阶算法而将步长h 取小。
➢ 变步长的Runge—Kutta Method
§2 Runge-Kutta Method
Rn1 ch p1 y( p1) xn
§2 Runge-Kutta Method
➢ 最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta
Method */ ：
yi 1
yi
h 6
(K1
2K2
2K3
K4)
K1 f ( xi, yi )
K2
f
( xi
h2 ,
yi
h 2
K1)
K3
f
( xi
h2 ,
yi
h 2
K2)
K4 f ( xi h, yi hK3 )
➢ Gill公式：4阶经典龙格-库塔公式的一种改进
yi1
yi
h 6
K1
2
2 K2 2
2 K3 K4
K1 f (xi , yi )K2f( xih 2
,
yi
h 2
K1)
K3
f
( xi
h 2
,
yi
2 1 2
hK1
1
2 2
hK2 )
yn
1
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
Ki
f ( xn ai h, yn h bij K j )
1
y( x i1 )
y (h) i 1
16
由此可得
y(xi1 )
(h)
y2 i 1
1 15
(
(h)
y2 i 1
y
(h) i 1
)
h
这表明以
()
y2 i 1
作为
y(xi1) 的近似值，其误差可用步
长折半前后两次计算结果的偏差
来判断所选步长是否适当
y y ( h ) 2 i 1
(h) i 1
当要求的数值精度为ε时：
i
21
1
K f ( x h, y hK hK )
3
i
3
i
31
1
32
2
... ...
K f ( x h, y hK hK ... hK )
m
i
m
m1
1
m2
2
m m 1
m 1
其中i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij
( i = 2, …, m; j = 1, …, i1 ) 均为待定系数，确 定这些系数的步骤与前面相似。
变化剧烈的情况下，总的计算工作量得到减少，结
果还是合算的。
➢ 高阶Runge—Kutta Method
§2 Runge-Kutta Method
y y h[ K K ... K ]
i 1
i
11
22
mm
K f(x , y )
1
ii
K f ( x h, y hK )
2
i
2
K4
f (xi h ,
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 )
ch5
当h值不大时，式中的系数c可近似地看作为常数。
然后将步长折半,即以为 h 2
步长,从节点xi出发,跨
两步到节点xi+1,再求得一个近似值
h ()
y2 i 1
,每跨一步的
截断误差是 c h 5 ,因此有
2
y( xi1 )
h
()
y(
xi
2 1
)
h
2c
h 5 2
()
这样
y( x i1 )
y2 i 1
以经典四阶龙格—库塔公式为例。从节点xn出发，以h为
步长求一近似值yn(h1) ,
y( xn1 )
y(h) n 1
ch5 ,
将步长折半，即取
h 2
为步长从xn
跨两步到xn
，求一近似
1
值y ( h／2) n 1
,
每跨一步的截断误差是c
h 2
5
,因此有
y( xn1 )
y ( h／2) n 1
2c
h 2
5
,
由上两式
y( xn1 )
y ( h／2) n 1
1.
y( xn1 )
y(h) n 1
16
y( xn1 )
y ( h／2) n 1
1 15
[
y ( h／2) n 1
y(h) n 1
].
公式
Euler公式 改进Euler公式
yn1 K1
yn f ( xn ,
hK1 yn )
完成的步数就增加了。这样会引起计算量的增大，
并且会引起舍入误差的大量积累与传播。因此微分
方程数值解法也有选择步长的问题。
以经典的四阶龙格-库塔法(7.20)为例。从节点
xi出发，先以h为步长求出一个近似值，记为
y (h) i 1
，
由于局部截断误差为 O(h5 ) ，故有
y( xi1 )
y (h) i 1
Q: 由局部截断误差可以看出，步长 h 越小，局部截断 误差越小；但步长减小，在一定求解范围（区间）内 要完成的步数就增加了，步数增加会引起计算量增大， 导致舍入误差积累。因此要选取适当的步长。
选择步长时要考虑两个问题： 1.如何衡量和检验计算结果的精度？ 2.如何根据所获得的精度处理步长？
HW: p.201 #6-8
（1）如果Δ>ε，反复将步长折半进行计算，直
至Δ<ε为止,并取其最后一次步长的计算结果作为y i 1
（2）如果Δ<ε，反复将步长加倍，直到Δ>ε
为止，并以上一次步长的计算结果yi作1 为
。
这种通过步长加倍或折半来处理步长的方法称为
变步长法。表面上看，为了选择步长，每一步都要
反复判断Δ，增加了计算工作量，但在方程的解y(x)
j 1
(i 2, 3L , p)
确定原则是使近似公式在( xn , yn )处的Taylor展开式与
y( x)在xn处的Taylor展开式的前面项尽可能多地重合。
7.4.6 变步长的龙格-库塔法
在微分方程的数值解中，选择适当的步长是非常
重要的。单从每一步看，步长越小，截断误差就越
小；但随着步长的缩小，在一定的求解区间内所要
7
O(h6 )
n8
O(hn2 )
由于龙格-库塔法的导出基于泰勒展开，故精度主要受
解函数的光滑性影响。对于光滑性不太好的解，最好 采用低阶算法而将步长h 取小。
➢ 变步长的Runge—Kutta Method
§2 Runge-Kutta Method
Rn1 ch p1 y( p1) xn
§2 Runge-Kutta Method
➢ 最常用为四级4阶经典龙格-库塔法 /* Classical Runge-Kutta
Method */ ：
yi 1
yi
h 6
(K1
2K2
2K3
K4)
K1 f ( xi, yi )
K2
f
( xi
h2 ,
yi
h 2
K1)
K3
f
( xi
h2 ,
yi
h 2
K2)
K4 f ( xi h, yi hK3 )
➢ Gill公式：4阶经典龙格-库塔公式的一种改进
yi1
yi
h 6
K1
2
2 K2 2
2 K3 K4
K1 f (xi , yi )K2f( xih 2
,
yi
h 2
K1)
K3
f
( xi
h 2
,
yi
2 1 2
hK1
1
2 2
hK2 )
yn
1
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
Ki
f ( xn ai h, yn h bij K j )
1
y( x i1 )
y (h) i 1
16
由此可得
y(xi1 )
(h)
y2 i 1
1 15
(
(h)
y2 i 1
y
(h) i 1
)
h
这表明以
()
y2 i 1
作为
y(xi1) 的近似值，其误差可用步
长折半前后两次计算结果的偏差
来判断所选步长是否适当
y y ( h ) 2 i 1
(h) i 1
当要求的数值精度为ε时：
i
21
1
K f ( x h, y hK hK )
3
i
3
i
31
1
32
2
... ...
K f ( x h, y hK hK ... hK )
m
i
m
m1
1
m2
2
m m 1
m 1
其中i ( i = 1, …, m )，i ( i = 2, …, m ) 和 ij
( i = 2, …, m; j = 1, …, i1 ) 均为待定系数，确 定这些系数的步骤与前面相似。
变化剧烈的情况下，总的计算工作量得到减少，结
果还是合算的。
➢ 高阶Runge—Kutta Method
§2 Runge-Kutta Method
y y h[ K K ... K ]
i 1
i
11
22
mm
K f(x , y )
1
ii
K f ( x h, y hK )
2
i
2
K4
f (xi h ,
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 )
ch5
当h值不大时，式中的系数c可近似地看作为常数。
然后将步长折半,即以为 h 2
步长,从节点xi出发,跨
两步到节点xi+1,再求得一个近似值
h ()
y2 i 1
,每跨一步的
截断误差是 c h 5 ,因此有
2
y( xi1 )
h
()
y(
xi
2 1
)
h
2c
h 5 2
()
这样
y( x i1 )
y2 i 1