龙格-库塔方法（Runge-Kutta）

龙格-库塔⽅法（Runge-Kutta）
龙格-库塔⽅法（Runge-Kutta）
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个，所应满⾜的⽅程为
这是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+1
K1 = 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)
end
format 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 。

于是，r=4，p=4的13个参数(c4不能为0)引出了多种⽅案和挑战，如：
参数优化使阶数增加到5阶，得到四阶五阶R-K ⽅法，matlab 中有程序ode45；
四级四阶R-K ⽅法的步长⾃动选取；结合新算法的应⽤算法构造；适应于新的领域实现求解； ……经典的四阶R-K ⽅法是：
(3.2.12)
其中也需满⾜12341c c c c +++=，这⾥为(1/6 2/6 2/6 1/6).它的局部截断误差
，故p=4，这是最常⽤的四阶R-K ⽅法，数学库中都有⽤此⽅法求解
初值问题的软件.这种⽅法的优点是精度较⾼，缺点是每步要算4个右端函数值，计算量较⼤.
例3.3 ⽤经典四阶R-K ⽅法解例1.1的初值问题，
仍取h=0.1，计算到
，并与改进Euler 法、梯形法在
处⽐较其误
差⼤⼩.
解⽤四阶R-K ⽅法公式(3.2.12)，此处
，于是当n=0时
于是，按公式
(3.2.12)可算出
此⽅法误差：
改进Euler 法误差：
梯形法误差：
可见四阶R-K ⽅法的精度⽐⼆阶⽅法⾼得多.⽤四阶R-K ⽅法求解初值问题(1.2.1)精度较⾼，但要从理论上给出误差的估计式则⽐较困难.那么应
如何判断计算结果的精度以及如何选择合适的步长h ？通常是通过不同步长在计算机上的计算结果近似估计.设
在处的值
,当
时，
的近似为()
1h n y ，于是由四阶R-K ⽅法有
若以为步长，计算两步到
，则有
于是得
即
或
(3.2.13) 它给出了误差的近似估计.如果(ε为给定精度)，则认为
以为步长的计算结果满⾜精度要求，若，则还可放⼤步
长.因此(3.2.13)提供了⾃动选择步长的⽅法.
附：经典的4级(也是4阶)R-K ⽅法的程序：
function y = DELGKT4_Rungkuta(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+1
K1 = 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/2 y(i-1)+K2*h/2]);
K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]);
y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6;
end
format short;
讲解：
求初值问题（1.2.1）的单步法主要是指Runge-Kutta法，本节主要讨论显式R－K⽅法，建⽴具体的计算公式使⽤的是Taylor 展开，形如（3.2.4）的显式R－K⽅法，当r＝1时就是Euler法，因此只要讨论的计算公式，在r
确定后如何推导公式都是⼀样的，只是r越⼤计算越复杂，为了掌握了解公式来源，只要以r＝2为例推导计算公式即可.因此本节重点就是⽤Taylor展开求出
r＝2的显式R-K⽅法的计算公式，由于⽅法的局部截断误差为（3.2.6），的右端有的项，要对它做Taylor展开，就要⽤到⼆元函数
的Taylor展开，按照⼆元函数Taylor级数
（3.2.14）
将它⽤到（3.2.6）的的展开式中，即可得到按升幂整理出的结果，对
r＝2的公式只能得到＝2阶的公式，即，于是2级R-K⽅法（3.2.5）的系数必须满⾜（3.2.7）给出的⽅程，它的解由（3.2.8）给出，
只要，求出的公式都是r=2的2阶R-K ⽅法.⽽常⽤的就是
得到的
改进Euler 法（3.1.11）和
得到的中点公式（3.2.9）.
我们知道，显式单步法的数值公式为：
如果取定右边的
则称该数值公式为显式Runge-Kutta ⽅法．（⼀）显式Runge-Kutta ⽅法 1. ⼆阶⼆级R-K ⽅法
11122()n n y y h 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 =++．四参数满⾜
举例：中点公式四个参数
212211
1, 0, 2c c a b ====
，即
1,(,)22n n n n n n h h y y hf x y f x y +??
=+++??
，
程序如下
function y = DELGKT2_mid(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+1
v1 = Funval(f,varvec,[x(i-1) y(i-1)]); t = y(i-1) + h*v1/2;
v2 = Funval(f,varvec,[x(i-1)+h/2 t]); y(i) = y(i-1)+h*v2; end
format short;
休恩（suen ）法 c1=1/4, c2=3/4, a2=b21=2/3
function y = DELGKT2_suen(f, h,a,b,y0,varvec) format long;
N = uint16((b-a)/h); y = zeros(N+1,1); y(1) = y0; x = a:h:b;
var = findsym(f); for i=2:N+1
K1 = Funval(f,varvec,[x(i-1) y(i-1)]);
K2 = Funval(f,varvec,[x(i-1)+2*h/3 y(i-1)+K1*2*h/3]); y(i) = y(i-1)+h*(K1+3*K2)/4; end
format short;
改进的Euler 法 c1=1/2, c2=1/2, a2=b21=1
function y = DEModifEuler(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+1
v1 = Funval(f,varvec,[x(i-1) y(i-1)]); t = y(i-1) + h*v1;
v2 = Funval(f,varvec,[x(i) t]); y(i) = y(i-1)+h*(v1+v2)/2; end
format short;
实验对⽐参数取法不同得到不同的⼆阶R-K ⽅法精确程度各异：这个已经可以做(同学们已做？！)
注意：⼆级R-K ⽅法只能达到⼆阶，⽽不可能达到三阶.因为r=2只有4个参数，要达到p=3则在(3.2.6)的展开式中要增加3项，即增加三个⽅程,加上(3.2.7)的三个⽅程，共计六个⽅程求4个待定参数，验证得出是⽆解的.
2．三阶三级R-K ⽅法
其中1(,)n n k f x y =，2221(,(,))n n n n k f x a h y b hf x y =++，⽽
共计8个参数待设定。

Taylor 展开分析局部截断误差，使得精度达到三级，即
误差为
4()O h 。

于是，r=3，p=3的8个参数所满⾜的关系式(c3不能为0)为
Kutta 法
function y = DELGKT3_kutta(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+1
K1 = 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; end
format short;
Suen 法
function y = DELGKT3_suen(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+1
K1 = Funval(f,varvec,[x(i-1) y(i-1)]);
K2 = Funval(f,varvec,[x(i-1)+h/3 y(i-1)+K1*h/3]); K3 = Funval(f,varvec,[x(i-1)+2*h/3 y(i-1)+K2*2*h/3]); y(i) = y(i-1)+h* (K1+3*K3)/4; end
format short;
3. 四阶四级R-K ⽅法
举例：经典四级Rung-Kutta法
程序
function y = DELGKT4_Rungkutta(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+1
K1 = 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/2 y(i-1)+K2*h/2]);
K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]);
y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6;
end
format short;
基尔(jer、Gear)法
function y = DELGKT4_jer(f, h,a,b,y0,varvec)
format long;
N = (b-a)/h;
C2 = sqrt(2);
y = zeros(N+1,1);
y(1) = y0;
x = a:h:b;
var = findsym(f);
for i=2:N+1
K1 = 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/2 y(i-1)+K1*h*(C2-1)/2+K2*h*(2-C2)/2]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)-K2*h*C2/2+K3*h* (2+C2)/2]);
y(i) = y(i-1)+h*(K1+(2-C2)*K2+(2+C2)*K3+K4)/6;
end
format short;
变型Rung-Kutta法
function y = DELGKT4_qt(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+1
K1 = Funval(f,varvec,[x(i-1) y(i-1)]);
K2 = Funval(f,varvec,[x(i-1)+h/3 y(i-1)+K1*h/3]);
K3 = Funval(f,varvec,[x(i-1)+2*h/3 y(i-1)-K1*h/3+K2*h]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*(K1-K2+K3)]); y(i) = y(i-1)+h* (K1+3*K2+3*K3+K4)/8; end
format short;
说明：
四阶四级R-K ⽅法的步长⾃动选取问题
⽤四阶R-K ⽅法求解初值问题(1.2.1)精度已经⾜够⾼了，但要从理论上给出误差
的估计式则⽐较困难.那么应如何判断计算结果的精度以及如何
选择合适的步长h ？通常是通过不同步长在计算机上的计算结果近似估计.设
()n y x 在n x 处的值()n n y y x =,当1n n x x h +=+时，1()n y x +的近似为()1h n y +，于是由四阶R-K ⽅法有
若以为步长，计算两步到，则有
于是得
即
或
(3.2.13)
它给出了误差的近似估计.如果 (ε为给定精度，如
1E-3)，则认为以为步长的计算结果满⾜精度要求，若，
则还可放⼤步长.因此(3.2.13)提供了⾃动选择步长的⽅法.
（⼆）隐式(半隐式)R-K ⽅法举例
⼀种半隐式举例，该格式可以通过除法化为显式单步法格式，受到除法复杂度限制，容易造成耗时，不稳，相对精度降低。

罗赛布诺克半隐式法
11122
1
1112
2121121[1(,)](,)
[1(,)](,)
n n y n n n n y n n n n y y w k w k k h ha f x y f x y k h ha f x c h y c k f x d h y d k +--=++??=-??=-++++? function y = DELSBRK(f, h,a,b,y0,varvec)
format long;
a1 = 1.40824829;a2 = 0.59175171; c1 = 0.17378667;c2 = c1;
w1 = -0.41315432;w2 = 1.41315432; N = (b-a)/h;
y = zeros(N+1,1); y(1) = y0; x = a:h:b;
var = findsym(f);
dy = diff(f, varvec(2)); for i=2:N+1
f1 = Funval(f,varvec,[x(i-1) y(i-1)]); dy1 = Funval(dy,varvec,[x(i-1) y(i-1)]); k1 = h*f1/(1-h*a1*dy1);
dy2 = Funval(dy,varvec,[x(i-1)+c1*h y(i-1)+c2*k1]); f2 = Funval(f,varvec,[x(i-1)+c1*h y(i-1)+c2*k1]); k2 = h*f2/(1-h*a2*dy2); y(i) = y(i-1)+w1*k1+w2*k2; end
format short;
3.3 单步法的收敛性与绝对稳定性
3.3.1 单步法的收敛性
定义 3.1 设y(x)是初值问题(1.2.1)
的精确解，
是单步法(3.2.2)在
处产⽣的近似解，若
则称⽅法(3.2.2)产⽣的数值解收敛于.
实际上，定义中是⼀固定点，当h→0时n→∞，n不是固定的.因
显然⽅法收敛，则在固定点处的整体误差，当p≥1
时.
下⾯定理给出⽅法(3.2.2)收敛的条件.
定理3.1设初值问题(1.2.1)的单步法(3.2.2)是p阶⽅法(p≥1)，且函数
对y满⾜Lipschitz条件，即存在常数L＞0，使对，均有
则⽅法(3.2.2)收敛，且.
定理证明略.
3.3.2 绝对稳定性
⽤单步法(3.2.2)求数值解，由于原始数据及计算过程舍⼊误
差影响，实际得到的不是⽽是，其中是误差，再计算下⼀步得到
以Euler法为例，若令,则
(3.3.1)
如果，则从计算到误差不增长，它是稳定的.但如果条件不满⾜就不稳定.
例3.4y′=-100y,y(0)=1,精确解为，⽤Euler法求解得
若取h=0.025，则,当，⽽
，显然计算是不稳定的.
如果⽤后退Euler法(3.1.5)解此例，仍取h=0.025,则
,即
显然当，计算是稳定的.
由此看到稳定性与⽅法有关，也与有关，在此例中.在研究⽅法的稳定性时，通常不必对⼀般的f(x,y)进⾏讨论，⽽只针对模型⽅程
(3.3.2)
这⾥可能为复数.规定是因为时微分⽅程(3.3.2)本⾝是不稳定的，⽽讨论数值⽅法(3.2.2)的稳定性，必须在微分⽅程本⾝稳定的前提下进⾏.另⼀⽅⾯，对初值问题(1.2.1)，若将f(x,y)在处线性展开，可得
于是⽅程(1.2.1)可近似表⽰为
它表明⽤模型⽅程(3.3.2)是合理的，⾄于模型⽅程(3.3.2)中所以⽤复数λ是因为初值问题(1.2.1)如果是⽅程组，即，则是(m×m)阶矩阵，其特征值可能是复数.当然对单个⽅程，λ就是实数，此时只要规定＜0
即可.
⽤单步法(3.2.2)解模型⽅程(3.3.2)可得到
(3.3.3)
其中依赖所选⽅法，如⽤Euler法，则
(3.3.4)
此时由(3.3.1)看到误差⽅程也为，与(3.3.4)是⼀样的.因此对⼀般单步法(3.2.2)误差⽅程也与(3.3.3)⼀致.下⾯再考虑⼆阶R-K⽅法有
对四阶R-K⽅法，可得
定义3.2 将单步法(3.2.2)⽤于解模型⽅程(3.3.2)，若得到(3.3.3)中的
则称⽅法是绝对稳定的.在复平⾯上复变量满⾜的区域，称为⽅法(3.2.2)的绝对稳定域，它与实轴的交点称为绝对稳定区间.
例如对Euler法，在复平⾯上是以(-1，0)为圆⼼，以1为半径的单位圆域内部，当为实数时，则得绝对稳定区间为，
在例3.4中时⽅法稳定，⽽例中取
因＜0，故有
.
h=0.025故不稳定.
对后退Euler法(3.1.5)，
因＜0，故，其绝对稳定域是以(1，0)为圆⼼的单位圆外部，绝对稳定区间为，即对任何h＞0⽅法都是绝对稳定的.
⼆阶R-K⽅法的绝对稳定区间为.
三阶R-K⽅法的绝对稳定区间为.
四阶R-K⽅法的绝对稳定区间为.
例3.5⽤经典四阶R-K⽅法计算初值问题
步长取h=0.1及0.2,给出计算误差并分析其稳定性.
解本题直接按R-K⽅法(3.2.12)的公式计算.因精确解为，其计算误差如表所⽰.
从计算结果看到，h=0.2时误差很⼤，这是由于在λ=-20,h=0.2时λh=-4，⽽四阶R-K⽅法的绝对稳定区间为［-2.785,0］,故
h=0.2时计算不稳定，误差很⼤.⽽h=0.1时=-2，其值在绝对稳定区间［-2.785,0］内，计算稳定，故结果是可靠的.
讲解：
由于微分⽅程初值问题数值解公式求出的解是⼀个逐次递推的过程，因此原始数据误差及计算过程舍⼊误差对解的影响就是数值⽅法绝对稳定性研究的问题，如果由计算误差不增长，⽅法就是绝对稳定的.为使问题得到简化通常就是将⽅法⽤于解模型⽅程（3.3.2），对于单步法得到的差分⽅程为，由于模型⽅程的，代⼊Euler法
，得，对⼆阶R-K⽅法，例如，⽤改进Euler法
于是
对三阶R-K⽅法有
对四阶R-K⽅法有
只要⽅法，就是绝对稳定的，这时的值当n增⼤式是减少的，故计算稳定.这时舍⼊误差影响可忽略不计，⽽当，则增⼤，⽅法不稳定，计算结果是不可靠的.因此⽤显式单步法必须使，也就是步长选择要满⾜这⼀要求.
对于隐式的梯形公式
将模型⽅程
，即
代⼊得
于是
注意，于是有
，对成⽴.
这就表明对任意步长h ，梯形法都是绝对稳定的.
3.4 练习题
2. 对于⼀阶微分⽅程初值问题3sin cos ,22
()0,
2dy
y x x dx y πππ?=+≤≤=??，⽤Euler 法，
改进Euler 法，⼆阶R-K ⽅法求解，并作图⽐较。

3. 对于⼀阶微分⽅程初值问题1,0 1.5
(0)1dy
y x x dx y ?=-++≤≤=?，⽤Euler 法，改
进Euler 法，⼆阶R-K 法，三阶R-K 法求解，并与真解作图⽐较，列出误
差对⽐表格。

4. 对于⼀阶微分⽅程初值问题1ln(1),01
(0)1dy
x x dx y ?=++≤≤=?，⽤三阶、四阶
R-K 法与真解作图对⽐。

5. 对于⼀阶微分⽅程初值问题2,01
(0)1dy
xy x dx y ?=≤≤=?，R-K 法的半隐式、四阶
与真解作图对⽐。

6. 对⼀阶微分⽅程初值问题8,24
(2)1dy
y x dx y ?=-≤≤=?进⾏稳定性分析：
a. 得出⼆阶R-K 法中步长的稳定区间(域)；
b. ⾃选两个步长(h1为区间内的数(如0.1), h2为区间外的数(如0.3))得到两个不同数值解和精确解作图⽐较，并列出误差表格。