龙格库塔方法ppt课件
合集下载
计算方法-常微分方程初值问题数值解法-Euler公式-龙格-库塔法市公开课获奖课件省名师示范课获奖课
当 x x1 时，得 y1 y0 f(x0 , y0 )(x1 x0 )
这么就取得了P1点旳坐标： (x1, y1) 。将y1作为y(x1)旳 近似值(想象(x1, y1) 在积分曲线y=y(x)上)
过点P1(x1,y1)，作积分曲线y=y(x)旳切线交直线x=x2于
P2点。注意切线 P1P2 旳斜率(近似)为 y(x1 ) f(x1 , y1 )
第9章 常微分方程初值问题数值解法
§9.1 引言
➢ 包括自变量、未知函数及未知函数旳导数旳方程称 为微分方程。
➢ 自变量个数只有一种旳微分方程称为常微分方 程。
微分方程中出现旳未知函数最高阶导数旳阶数 称为微分方程旳阶数。
假如未知函数y及其各阶导数
y, y, … , y(n)
都是一次旳,则称其为线性旳,不然称为非线性旳。
xi1 xi1 f[xi , y(xi )]
代入上式,并用yi近似替代式中y(xi)即可得到 两步欧拉公式
yi1 yi1 2hf(xi , yi ) ( 9.7 )
【注】欧拉措施和梯形措施，都是单步法,其特点是 在计算yi+1时只用到前一步旳信息yi; 而两步欧拉公式 (9.7)中除了yi外,还用到更前一步旳 信息yi-1,即调用了前两步旳信息。
当 x xi1 时,得
yi1 yi f(xi , yi )(xi1 xi )
这么,从x0逐一算出 x1 , x2 , … xn
相应旳数值解
y1 , y 2 , … yn
就取得了一系列旳点: P1, P1,…,Pn。 从图形上看,就取得了一条近似于曲线y=y(x)
旳折线 P1P2P3 … Pn 。
xi x0 ih, i 1,2, … , n
数值解法需要把连续性旳问题加以离散化，从 而求出离散节点旳数值解。
这么就取得了P1点旳坐标： (x1, y1) 。将y1作为y(x1)旳 近似值(想象(x1, y1) 在积分曲线y=y(x)上)
过点P1(x1,y1)，作积分曲线y=y(x)旳切线交直线x=x2于
P2点。注意切线 P1P2 旳斜率(近似)为 y(x1 ) f(x1 , y1 )
第9章 常微分方程初值问题数值解法
§9.1 引言
➢ 包括自变量、未知函数及未知函数旳导数旳方程称 为微分方程。
➢ 自变量个数只有一种旳微分方程称为常微分方 程。
微分方程中出现旳未知函数最高阶导数旳阶数 称为微分方程旳阶数。
假如未知函数y及其各阶导数
y, y, … , y(n)
都是一次旳,则称其为线性旳,不然称为非线性旳。
xi1 xi1 f[xi , y(xi )]
代入上式,并用yi近似替代式中y(xi)即可得到 两步欧拉公式
yi1 yi1 2hf(xi , yi ) ( 9.7 )
【注】欧拉措施和梯形措施，都是单步法,其特点是 在计算yi+1时只用到前一步旳信息yi; 而两步欧拉公式 (9.7)中除了yi外,还用到更前一步旳 信息yi-1,即调用了前两步旳信息。
当 x xi1 时,得
yi1 yi f(xi , yi )(xi1 xi )
这么,从x0逐一算出 x1 , x2 , … xn
相应旳数值解
y1 , y 2 , … yn
就取得了一系列旳点: P1, P1,…,Pn。 从图形上看,就取得了一条近似于曲线y=y(x)
旳折线 P1P2P3 … Pn 。
xi x0 ih, i 1,2, … , n
数值解法需要把连续性旳问题加以离散化，从 而求出离散节点旳数值解。
龙格库塔
其中，
xn p xn ph, 0 q 1
xnq xn qh, p q 1
于是就可以构造所谓的三阶龙格－库塔格式 ，下列库塔 格式是其中的一种：
yn
1
yn
K1
f
(xn ,
h 6 (K1 yn )
4K2
K3)
K2 f (xn p , yn hK1)
K
3
f ( xn1, yn
h(K1 2K2 ))
3.3 龙格-库塔讲义
——《数值分析简明教程》
主讲：王老师
黄淮学院数学科学系
2011-3-16
暮春三月，草长莺飞
1
❖ 考察差商 y(xn1) y(xn ) ，根据微分中值定理存在
点 ，使得
h
1.龙格 库塔y' 的f 设计思想
从而利用所给方程 得 y(xn1) y(xn ) y'( )
h
yxn＋1 yxn hf , y ，xn xn1
y y(x)
K*
xn
2011-3-16
暮春三月，草长莺飞
x n 1
x
3
2.二阶龙格-库塔方法
随意考察区间 xn, xn1 内 一点
xn＋p xn ph, 0 p 1
用两个点 xn , xn p 的斜率的K1 , K2的加权平均 代替平均斜率 K ,于是我们就得到如下计算 格式:
yn1 K1 f
则
f (Kxn ,2yn )
如下图
y
y y(x)
K1
yn1 yn hf (xn , yn )
xn
xn 1
x
2011-3-16
暮春三月，草长莺飞
5
适当选取它们的值，就可使上述格式有较高
xn p xn ph, 0 q 1
xnq xn qh, p q 1
于是就可以构造所谓的三阶龙格－库塔格式 ，下列库塔 格式是其中的一种：
yn
1
yn
K1
f
(xn ,
h 6 (K1 yn )
4K2
K3)
K2 f (xn p , yn hK1)
K
3
f ( xn1, yn
h(K1 2K2 ))
3.3 龙格-库塔讲义
——《数值分析简明教程》
主讲：王老师
黄淮学院数学科学系
2011-3-16
暮春三月，草长莺飞
1
❖ 考察差商 y(xn1) y(xn ) ，根据微分中值定理存在
点 ，使得
h
1.龙格 库塔y' 的f 设计思想
从而利用所给方程 得 y(xn1) y(xn ) y'( )
h
yxn＋1 yxn hf , y ，xn xn1
y y(x)
K*
xn
2011-3-16
暮春三月，草长莺飞
x n 1
x
3
2.二阶龙格-库塔方法
随意考察区间 xn, xn1 内 一点
xn＋p xn ph, 0 p 1
用两个点 xn , xn p 的斜率的K1 , K2的加权平均 代替平均斜率 K ,于是我们就得到如下计算 格式:
yn1 K1 f
则
f (Kxn ,2yn )
如下图
y
y y(x)
K1
yn1 yn hf (xn , yn )
xn
xn 1
x
2011-3-16
暮春三月，草长莺飞
5
适当选取它们的值，就可使上述格式有较高
第四讲龙格-库塔方法
第四讲龙格-库塔⽅法龙格-库塔⽅法3.2 Runge-Kutta法3.2.1 显式Runge-Kutta法的⼀般形式上节已给出与初值问题(1.2.1)等价的积分形式(3.2.1)只要对右端积分⽤不同的数值求积公式近似就可得到不同的求解初值问题(1.2.1)的数值⽅法，若⽤显式单步法(3.2.2)当，即数值求积⽤左矩形公式，它就是Euler法(3.1.2)，⽅法只有⼀阶精度，若取(3.2.3)就是改进Euler法，这时数值求积公式是梯形公式的⼀种近似，计算时要⽤⼆个右端函数f的值，但⽅法是⼆阶精度的.若要得到更⾼阶的公式，则求积分时必须⽤更多的f值，根据数值积分公式，可将(3.2.1)右端积分表⽰为注意，右端f中还不能直接得到，需要像改进Euler法(3.1.11)⼀样，⽤前⾯已算得的f值表⽰为(3.2.3)，⼀般情况可将(3.2.2)的表⽰为(3.2.4)其中这⾥均为待定常数，公式(3.2.2)，(3.2.4)称为r级的显式Runge-Kutta法，简称R-K⽅法.它每步计算r个f值(即)，⽽k由前⾯(i-1)个已算出的表⽰，故公式是显式的.例i如当r=2时，公式可表⽰为(3.2.5) 其中.改进Euler 法(3.1.11)就是⼀个⼆级显式R-K ⽅法.参数取不同的值，可得到不同公式.3.2.2 ⼆、三级显式R-K ⽅法对r=2的显式R-K ⽅法(3.2.5)，要求选择参数，使公式的精度阶p 尽量⾼，由局部截断误差定义11122211()()[(,())(,)]n n n n n n n T y x y x h c f x y x c f x a h y b hk ++=--+++ (3.2.6)令，对(3.2.6)式在处按Taylor 公式展开，由于将上述结果代⼊(3.2.6)得要使公式(3.2.5)具有的阶p=2,即，必须(3.2.7)即由此三式求的解不唯⼀.因r=2，由（3.2.5）式可知，于是有解(3.2.8)它表明使(3.2.5)具有⼆阶的⽅法很多，只要都可得到⼆阶精度R-K⽅法.若取，则，则得改进Euler法(3.1.11)，若取，则得，此时(3.2.5)为(3.2.9)其中称为中点公式.改进的Euler法(3.1.11)及中点公式(3.2.9)是两个常⽤的⼆级R-K⽅法，注意⼆级R-K⽅法只能达到⼆阶，⽽不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个⽅程,加上(3.2.7)的三个⽅程，共计六个⽅程求4个待定参数，验证得出是⽆解的.当然r=2,p=2的R-K⽅法(3.2.5)当取其他数时，也可得到其他公式，但系数较复杂，⼀般不再给出.对r=3的情形，要计算三个k值，即其中将按⼆元函数在处按Taylor公式展开，然后代⼊局部截断误差表达式，可得可得三阶⽅法，其系数共有8个，所应满⾜的⽅程为(3.2.10)这是8个未知数6个⽅程的⽅程组，解也是不唯⼀的，通常.⼀种常见的三级三阶R-K⽅法是下⾯的三级Kutta⽅法：(3.2.11)附：R-K 的三级Kutta ⽅法程序如下function y = DELGKT3_kuta(f, h,a,b,y0,varvec) format long; N = (b-a)/h;y = zeros(N+1,1); y(1) = y0; x = a:h:b;var = findsym(f); for i=2:N+1K1 = Funval(f,varvec,[x(i-1) y(i-1)]);K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h y(i-1)-h*K1+K2*2*h]);y(i) = y(i-1)+h*(K1+4*K2+K3)/6; %满⾜c1+c2+c3=1,(1/6 4/6 1/6)endformat short; 3.2.3 四阶R-K ⽅法及步长的⾃动选择利⽤⼆元函数Taylor 展开式可以确定(3.2.4)中r=4,p=4的R-K ⽅法，其迭代公式为111223344()n n y y h c k c k c k c k +=++++其中1(,)n n k f x y =，2221(,(,))n n n n k f x a h y b hf x y =++，⽽33311322(,)n n k f x a h y b hk b hk =+++ 44411422433(,)n n k f x a h y b hk b hk b hk =++++共计13个参数待定，Taylor 展开分析局部截断误差，使得精度达到四阶，即误差为5()O h 。
Runge-Kutta法
3½×RK· ¨
1.422
¸ ĽøEuler· ¨
1.42
1.418 0.2965 0.297 0.2975 0.298 0.2985 0.299 0.2995 0.3 0.3005 17
§ 7.2 Runge-Kutta法
龙格-库塔(Runge-Kutta)方法简称R-K法，是一种应用较广的 高精度的单步法。
所谓单步法就是在计算yi时只用到前一步信息yi-1的方法。
本节介绍R-K法的构造原理、常用公式。
1
一、Runge-Kutta方法的构造原理 对于常微分方程的初值问题
y f ( x , y ) y( a ) y0 a xb
另一方面
h2 y( xi 1 ) y( xi ) hy( xi ) y( xi ) o(h3 ) 2
9
式（6）的局部截断误差:
en (h) y ( xi 1 ) yi 1 1 h 1 (1 2 ) y( xi ) h2 22 y( xi ) o(h3 ) 2
1 2 1 1 2 2 2
有无穷多组解，从而可以得到许多具体的二阶R -K公式 如：
1 2 1，2 1
1 1 0， 2 1， 2 2
10
用类似的方法也可以构造三阶R-K公式
h yi yi 1 ( K1 4 K2 K3 ) 6 K1 f ( xi 1, yi 1 ) h h K 2 f ( xi 1 , yi 1 K1 ) 2 2 K3 f ( xi 1 h, yi 1 h(2 K2 K1 ))
因而方法(8)有4阶精度
12
例1. 使用高阶R-K方法计算初值问题
(完整版)第二节龙格-库塔方法
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
高等流体力学龙格库塔法求解实际问题汇报PPT
结论：水平方向速度u随y的变化曲线，当y很小时速度梯 度很大，随着y的增加速度梯度逐渐减小，水平方向速度 u逐渐趋于主流速度U=25m/s。 对于平板的不同位置（取x=0.01m,x=0.5m,x=2m,x=5m 处），水平方向速度u趋于主流速度U的快慢也不尽相同， 随着X的增加u的变化逐渐减慢。
由Y方向的速度公式 v Ue [ f '() f ()]，以及 y 2 x
k2 2
，yn +
l2 2
，zn +
m2 2
，tn +
t ） 2
l4 =tf（2 xn +k3，yn +l3，zn +m3，tn +t）
%求解F得到A x=0;y=0;z=1;h=0.001;t0=0;tn=6;j=1;t=t0:h:tn;N=(tn-t0)/h; x1=zeros(N,1);y1=zeros(N,1);z1=zeros(N,1); for i=1:N+1
引入流函数，将未知 函数的数量降为1个， 并做无因次相似性变 换，变为常微分方程， 定解条件相应变换。
将初边值问题转换为 初值问题
4.通过 编程进行流场分析
我们小组主要对 X Y 方向的速度、边界层厚度、排挤厚度、 以及板面切应力进行了求解分析。
程序实现
f =x，f =y，f =z f =y，f =z，f =-xz f 0 A
%求边界层厚度 板面切应力 排挤厚度 动量切应力 x=0:0.001:5; Ue=25;n=0.00002593; p=1.205; m1=5.138*sqrt(2*n/p/ Ue *x); m2=sqrt(2*n/p/25*x).*(sum(1-y1))*h; m3=sqrt(2*n/p/15*x).*(sum(y1.*(1-y1)))*h m4=0.1656*x.^(-0.5); figure(5);plot(x,m1); title(' 边 界 层 厚 度 δ 随 x 的 变 化 曲 线 '); xlabel('x(m)'); ylabel('y(m)--边界层厚度'); grid on;
龙格-库塔方法---数学课的讲解共34页
龙格-库塔方法---数学课的讲解
16、人民应该为法律而战斗，就像为 了城墙 而战斗 一样。 ——赫 拉克利 特 17、人类对于不公正的行为加以指责 ，并非 因为他 们愿意 做出这 种行为 ，而是 惟恐自 己会成 为这种 行为的 牺牲者 。—— 柏拉图 18、制定法律法令，就是为了不让强 者做什 么事都 横行霸 道。— —奥维 德 19、法律是社会的习惯和思想的结晶 。—— 托·伍·威尔逊 20、人们嘴上挂着的法律，其真实含 义是财 富。— —爱献 生
21、要知道对好事的称颂过于夸大，也会招来人们的反感轻蔑和嫉妒。——培根 22、业精于勤，荒于嬉；行成于思，毁于随。——韩愈
23、一切节省，归根到底都归结为时间的节省。——马克思 24、意志命运往往背道而驰，决心到最后会全部推倒。——莎士比亚
25、学
16、人民应该为法律而战斗，就像为 了城墙 而战斗 一样。 ——赫 拉克利 特 17、人类对于不公正的行为加以指责 ，并非 因为他 们愿意 做出这 种行为 ，而是 惟恐自 己会成 为这种 行为的 牺牲者 。—— 柏拉图 18、制定法律法令，就是为了不让强 者做什 么事都 横行霸 道。— —奥维 德 19、法律是社会的习惯和思想的结晶 。—— 托·伍·威尔逊 20、人们嘴上挂着的法律，其真实含 义是财 富。— —爱献 生
21、要知道对好事的称颂过于夸大，也会招来人们的反感轻蔑和嫉妒。——培根 22、业精于勤，荒于嬉；行成于思，毁于随。——韩愈
23、一切节省，归根到底都归结为时间的节省。——马克思 24、意志命运往往背道而驰，决心到最后会全部推倒。——莎士比亚
25、学
第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 阶精度。
《微分方程的数值解法maab四阶龙格—库塔法》PPT模板课件
k1 f (tn , yn )
k2
f (tn
1 2
h,
yn
h 2 k1)
k3
f (tn
1 2
h,
yn
h 2
k2
)
k 4 f (tn h , y n hk 3 )
四 阶 Runge-Kutta 法计算流程图
开始
h 初始条件：t
迭代次数：
；y
0
N
0
积分步长：
tn t0
for i = 1 : N
解析解： x x x1 3 2(((ttt))) 0 .0 8 1 1 2 P 8 k 0siw n t) (2 .6 3 0 3 3 P k 0siw n t) (0 .2 12 2 2 P k 0siw n t)(
第一个质量的位移响应时程
各种solver 解算指令的特点
解法指令 解题类 型
特点
ode45 非刚性 采用4、5阶Runge－Kutta法
适合场合 大多数场合的首选算法
ode23 非刚性 采用Adams算法
较低精度（10－3）场合
ode113 非刚性
ode23t ode15s
适度刚 性
刚性
ode23s 刚性 ode23tb 刚性
y2
(0)
(2.2) (2.3)
例：著名的Van der Pol方程
y (y21)y y0
令
y1y,y2y
Y
y y
1 2
降为一阶 Yyy12(y12y12)y2y1
初始条件
Y0
y1(0) y2(0)
y10 y20
3. 根据式(2.2)编写计算导数的M函数文件ODE文件
第3讲(龙格-库塔方法)
h f yi [ f ( x i , yi ) f ( xi , yi ) ( x i , yi ) h 2 x f ( xi , yi ) h f ( xi , yi ) O( h2 )] y
易见，它与二阶泰勒级数方法仅相差 O( h3 )!
这一分析给我们提供了一个重要信息，那就是 我们所遇到的泰勒级数方法中求导数的困难是可以 克服的，改进的欧拉方法就没有用到导数，而是借 助于函数在某些点处的值 (复合函数的思想）。
又 y( x ) df ( x , y( x )) f x f y y f x f y f dx
故二阶泰勒级数方法为 h2 yi 1 yi h f ( xi , yi ) ( f x ( xi , yi ) f ( xi , yi ) f y ( xi , yi )) 2! 更高阶方法更复杂，主要是求导复杂！
yi 1
h2 hk ( k ) yi h y y yi i i 2! k!
这样的数值方法称为k 阶泰勒级数方法。
yi 1
h2 hk ( k ) yi h y y yi i i 2! k!
泰勒级数方法也是单步法，且其局部截断误差为
h2 hk ( k ) LTE y( xi 1 ) y( xi ) hy( xi ) y( xi ) y ( x i ) 2! k!
第二节 龙格-库塔方法
（Runge-Kutta）
根据局部截断误差与整体误差的关系可知， 局部截断误差的阶是衡量一个方法优劣的重要依据。 考虑用提高局部截断误差的阶来提高数值方法的 精度。 泰勒级数法 龙格―库塔方法
一、泰勒级数方法
d y f ( x, y ), x I 如果初值问题 d x 的精确解 y(x) 在 I y( x ) y 0 0
易见，它与二阶泰勒级数方法仅相差 O( h3 )!
这一分析给我们提供了一个重要信息，那就是 我们所遇到的泰勒级数方法中求导数的困难是可以 克服的，改进的欧拉方法就没有用到导数，而是借 助于函数在某些点处的值 (复合函数的思想）。
又 y( x ) df ( x , y( x )) f x f y y f x f y f dx
故二阶泰勒级数方法为 h2 yi 1 yi h f ( xi , yi ) ( f x ( xi , yi ) f ( xi , yi ) f y ( xi , yi )) 2! 更高阶方法更复杂，主要是求导复杂！
yi 1
h2 hk ( k ) yi h y y yi i i 2! k!
这样的数值方法称为k 阶泰勒级数方法。
yi 1
h2 hk ( k ) yi h y y yi i i 2! k!
泰勒级数方法也是单步法，且其局部截断误差为
h2 hk ( k ) LTE y( xi 1 ) y( xi ) hy( xi ) y( xi ) y ( x i ) 2! k!
第二节 龙格-库塔方法
（Runge-Kutta）
根据局部截断误差与整体误差的关系可知， 局部截断误差的阶是衡量一个方法优劣的重要依据。 考虑用提高局部截断误差的阶来提高数值方法的 精度。 泰勒级数法 龙格―库塔方法
一、泰勒级数方法
d y f ( x, y ), x I 如果初值问题 d x 的精确解 y(x) 在 I y( x ) y 0 0
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)
常微分方程龙格-库塔方法
得中点公式 （8.2.4）
h h yn1 yn hf xn ，yn f xn，yn ， 2 2
取c=2/3得Heun公式
h 2 2 yn1 yn f xn，yn 3 f xn h，yn hf xn，yn 4 3 3
L
i 1
它与显式R-K公式的区别在于：显式公式中，对系数 aij 求和的上限是 i 1 ，从
aij构成的矩阵是一个严格下三角阵。而在隐式公式中，对系数 aij 求和的上 2， ，L 限是L，从而 aij 构成的矩阵是方阵，需要用迭代法求出近似斜率K i i 1，
而 推导隐式公式的思路和方法与显式公式法类似。
第八章常微分方程数值解法
8.2 Runge-Kutta 方法
8.2.1 Runge-Kutta方法的基本思想
8.2.2 几类显式Runge-Kutta方法
第八章常微分方程数值解法
8.2.1 Runge-Kutta方法的基本思想
显式Euler方法是最简单的单步法，它是一阶的，它可以看作Talylor展开后
f
，而不是高阶导数，将它们作线性
组合得到平均斜率，将其与解的Taylor展开相比较，使前面若干项吻合，从而 得到具有一定阶的方法。这就是 yn1 yn h i K i， Runge-Kutta方法的基本思想，其一般形式 i 1 为 (8.2.1) K1 f xn，yn ，
i 1 Ki f x c h , y c h a K 3， ，L， n i n i ij j ，i 2， j 1 L
到 K i* 。参数i，ci 和 aij 待定，确定它们的原则和方法是：将（8.2.2）式中
* 的 y xn1 在 xn 处作Taylor展开，将 K i 在 xn，y xn 处作二元Taylor展开，
就得到四阶龙格库塔法27页PPT
就得到四阶龙格库塔法
36、如果我们国家的法律中只有某种 神灵， 而不是 殚精竭 虑将神 灵揉进 宪法， 总体上 来说， 法律就 会更好 。—— 马克·吐 温 37、纲纪废弃之日，便是暴政兴起之 时。— —威·皮 物特
38、若是没有公众舆论的支持，法律 是丝毫 没有力 量的。 ——菲 力普斯 39、一个判例造出另一个判例，它们 迅速累 聚，进 而变成 法律。 ——朱 尼厄斯
40、人类法律，事物有规律，这是不 容忽视 的。— —爱献 生
66、节制使快乐增加并使享受加强。 ——德 谟克利 特 67、今天应做的事没有做，明天再早也 是耽误 了。——裴斯 泰洛齐 68、决定一个人的一生，以及整个命运 的，只 是一瞬 之间。 ——歌 德 69、懒人无法享受休息之乐。——拉布 克 70、浪费时间是一桩
36、如果我们国家的法律中只有某种 神灵， 而不是 殚精竭 虑将神 灵揉进 宪法， 总体上 来说， 法律就 会更好 。—— 马克·吐 温 37、纲纪废弃之日，便是暴政兴起之 时。— —威·皮 物特
38、若是没有公众舆论的支持，法律 是丝毫 没有力 量的。 ——菲 力普斯 39、一个判例造出另一个判例，它们 迅速累 聚，进 而变成 法律。 ——朱 尼厄斯
40、人类法律，事物有规律，这是不 容忽视 的。— —爱献 生
66、节制使快乐增加并使享受加强。 ——德 谟克利 特 67、今天应做的事没有做，明天再早也 是耽误 了。——裴斯 泰洛齐 68、决定一个人的一生，以及整个命运 的，只 是一瞬 之间。 ——歌 德 69、懒人无法享受休息之乐。——拉布 克 70、浪费时间是一桩
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
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);
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为迭代步数累 加值
5
误差分析：
四阶R-K方法的每一步需要计算四 次函数值f，可以证明其局部截断误差 为O(h5).
注意上述公式对于标量或者向 量函数(y可以是向量)都适用。
6
注意的问题
Runge-Kutta法的主要运算在于计算 Ki 的值，即计 f 的值。计算量与可达到的最高精度阶数的关系：
每步须算Ki 的 个数
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方法精度更高
2
四阶RungeKutta方法
3
4
这样，下一个值(yn+1)由现在的值(yn)加上时间间隔 (h)和一个估算的斜率的乘积 决定。该斜率是以下斜率的加权平均： k1是时间段开始时的斜率； k2是时间段中点的斜率，通过欧拉法采用斜率 k1 来 决定 y在点 tn + h/2的值； k3也是中点的斜率，但是这次采用斜率 k2决定 y值； k4是时间段终点的斜率，其 y值用 k3 决定。 当四个斜率取平均时，中点的斜率有更大的权值：
Runge-Kuttua方法和matlab原 理
1
龙格－库塔法（Runge-Kutta） 数值分析中，龙格－库塔法（Runge-Kutta） 是用于模拟常微分方程的解的重要 的一类隐式或显式迭代法。这些技术由数学家卡 尔·龙格和马丁·威尔海姆·库塔于1900年左右发 明。
经典四阶龙格库塔法 龙格库塔法的家族中的一个成员如此常用，以至 于经常被称为“RK4”或者就是 “龙格库塔法”。
实施方案
• 以经典四阶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)
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
14
五、变步长Runge-Kutta 方法 • 从每一步看，步长越小，截断误差
越小；但随着步长的缩小，在一定 求解范围内所要完成的步数就会增 加，步数的增加不但引起计算量的 增大，而且可能导致舍入误差的严 重积累，因此需要选择步长
•如何衡量和检验计算结果的精度 •如何依据所判定的精度来处理步长
15
每跨一步截断误差是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
3
4
5
6
7 n8
可达到的最高 精度
O(h2) O(h3)
O(h4 )
O(h4 )
O(h5 )
O(h6 )
O(hn2 )
R-K(高阶)方法不唯一,选择不同的 参数能得到
R不-K同方的法R的-K推公导式是基于Taylor展开 法，因而要求
解具有较好的光滑性，如果光滑性 较差精度可
能不如改进Euler方法,最好采用低 7
解 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)|
四阶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
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
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);
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为迭代步数累 加值
5
误差分析：
四阶R-K方法的每一步需要计算四 次函数值f，可以证明其局部截断误差 为O(h5).
注意上述公式对于标量或者向 量函数(y可以是向量)都适用。
6
注意的问题
Runge-Kutta法的主要运算在于计算 Ki 的值，即计 f 的值。计算量与可达到的最高精度阶数的关系：
每步须算Ki 的 个数
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方法精度更高
2
四阶RungeKutta方法
3
4
这样，下一个值(yn+1)由现在的值(yn)加上时间间隔 (h)和一个估算的斜率的乘积 决定。该斜率是以下斜率的加权平均： k1是时间段开始时的斜率； k2是时间段中点的斜率，通过欧拉法采用斜率 k1 来 决定 y在点 tn + h/2的值； k3也是中点的斜率，但是这次采用斜率 k2决定 y值； k4是时间段终点的斜率，其 y值用 k3 决定。 当四个斜率取平均时，中点的斜率有更大的权值：
Runge-Kuttua方法和matlab原 理
1
龙格－库塔法（Runge-Kutta） 数值分析中，龙格－库塔法（Runge-Kutta） 是用于模拟常微分方程的解的重要 的一类隐式或显式迭代法。这些技术由数学家卡 尔·龙格和马丁·威尔海姆·库塔于1900年左右发 明。
经典四阶龙格库塔法 龙格库塔法的家族中的一个成员如此常用，以至 于经常被称为“RK4”或者就是 “龙格库塔法”。
实施方案
• 以经典四阶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)
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
14
五、变步长Runge-Kutta 方法 • 从每一步看，步长越小，截断误差
越小；但随着步长的缩小，在一定 求解范围内所要完成的步数就会增 加，步数的增加不但引起计算量的 增大，而且可能导致舍入误差的严 重积累，因此需要选择步长
•如何衡量和检验计算结果的精度 •如何依据所判定的精度来处理步长
15
每跨一步截断误差是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
3
4
5
6
7 n8
可达到的最高 精度
O(h2) O(h3)
O(h4 )
O(h4 )
O(h5 )
O(h6 )
O(hn2 )
R-K(高阶)方法不唯一,选择不同的 参数能得到
R不-K同方的法R的-K推公导式是基于Taylor展开 法，因而要求
解具有较好的光滑性，如果光滑性 较差精度可
能不如改进Euler方法,最好采用低 7
解 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)|
四阶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
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