2Runge-Kutta方法
龙格库塔法解轨道参数
龙格库塔法解轨道参数引言龙格库塔法（Runge-Kutta method）是一种常用的数值计算方法，用于求解常微分方程的数值解。
在天体力学中，我们经常需要通过数值方法来计算天体的轨道参数，如轨道椭圆的长短轴、离心率、倾角等。
本文将介绍龙格库塔法在解轨道参数中的应用，并详细探讨该方法的原理和实现过程。
基本原理龙格库塔法是一种迭代求解的方法，在每个时间步长内利用当前的状态来估计下一个状态。
具体而言，龙格库塔法将微分方程的求解问题转化为一个迭代的求解问题，通过逐步迭代来逼近精确解。
在解轨道参数的问题中，我们通常需要根据已知的初始条件以及天体的质量和力学模型来求解天体的轨道参数。
常用的力学模型有开普勒模型和牛顿模型。
龙格库塔法可以根据力学模型的不同进行相应的求解。
开普勒模型下的轨道参数求解步骤一：确定初始条件在使用龙格库塔法求解轨道参数之前，我们需要确定一些初始条件。
这些初始条件包括天体的质量、位置和速度。
步骤二：选择时间步长在求解过程中，我们需要选择一个合适的时间步长。
时间步长越小，计算的精度会越高，但计算的时间会增加。
步骤三：迭代求解利用龙格库塔法进行迭代求解的具体步骤如下：1.根据当前时刻的位置和速度，计算天体在该时刻的加速度。
2.根据当前时刻的位置、速度和加速度，计算下一个时刻的位置和速度。
3.更新当前时刻的位置和速度为新的位置和速度。
4.重复上述步骤，直到达到指定的终止条件。
步骤四：计算轨道参数通过迭代求解，我们可以得到天体在不同时刻的位置和速度。
根据这些位置和速度，我们可以计算出轨道参数，如离心率、倾角、长短轴等。
常用的轨道参数计算公式如下：1.离心率：e=√1+2El2μ(GM⊕)2)2.倾角：i=arccos(ℎzℎ3.长轴：a=−μT22E4.短轴：b=a√1−e2其中，E表示能量，l表示轨道角动量，μ表示标准引力参数，G表示引力常数，M⊕表示地球的质量，ℎ和ℎz分别表示轨道角动量和轨道角动量在z轴上的分量。
龙格-库塔方法(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值(即由前面(i-1)个已算出的表示，故公式是显式的.例)，而ki如当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 。
§5.3 Runge-Kutta方法ppt课件
k1
f ( xi , yi )
k2
f ( xi
ph, yi
phk1 )
yi1 yi h (1 )k1 k2 k3
k1 f ( xi , yi ) k2 f ( xi ph, yi phk1 )
为了预报点 xi的q 斜率值k3,在区间 xi内, x有iq两
h 6
(K1
2K2
2K3
K4)
K1 f ( xi, yi )
K2
f
( xi
h 2
,
yi
h 2
K1)
K3
f
( xi
h 2
,
yi
h 2
K2)
K4 f ( xi h, yi hK3 )
四阶经典的龙格-库塔格式
例 用标准4阶R-K公式求：
y( x) x y( x) 0 x 1
在y((xxi i,y1i))点y作i1TayOlo(rh展3)开fx (
x,
y)
f y ( x,
y)
f
(x,
y)
K2 f ( xi ph, yi phK1 ) f ( xi , yi ) phfx ( xi , yi ) phK1 f y ( xi , yi ) O(h2 ) y( xi ) phy( xi ) O(h2 )
5.3.2 Runge-Kutta方法的基本思想
单步递推法的基本思想是从 ( xi , yi ) 点出发，以 某一斜率沿直线达到 ( xi+1 , yi+1 ) 点。欧拉法及 其各种变形所能达到的最高精度为2阶。
矩阵特征值与特征向量的计算-Rung-Kutta方法
每步须算Ki 的个数 2 3
4
5
6
可达到的最高精度 O(h2 ) O(h3 ) O(h4 ) O(h4 ) O(h5 )
7
O(h6 )
n8
O(hn−2 )
由于龙格-库塔法的导出基于泰勒展开，故精度主要受
解函数的光滑性影响。对于光滑性不太好的解，最好 采用低阶算法而将步长h 取小。
R − K方法的主要优缺点
二级R-K方法
二级R-K方法的形式为
其局部截断误差为 将 中的各项作Taylor展开
Taylor展开有 可得 令
二级R-K能达到的 最高阶数是二阶
常用的二级二阶R-K方法
取
，得
该方法称为改进的Euler公式（梯形公式的预估校正格式）
取
，得
该方法称为中点公式
取
，得
该方法称为Heun(休恩)方法
当 为实数时，得Euler法的绝对稳定区间是
二级R-K方法的绝对稳定区间
二阶二级R-K方法的计算公式为
由此可知，二阶二级R-K方法的绝对稳定区间是 当 为实数时，得绝对稳定区间是
一些常用方法的绝对稳定区间
R-K法的绝对稳定区域
k = 4 • 3. k =3
• 2. k=2 k = 1 • 1.
题相容的充分必要条件是该单步法至少是一阶方
法。
我们本章讨论的数值方法都是与原初值问题相容的！
收敛性
定义：对任意固定的 步法产生的解 ，均有
, 若初值问题的单
则称该方法是收敛的。
我们本章讨论的数值方法都是收敛的！
收敛性判别
定理7.3：设增量函数
在区域
上连续，并对变量y和h满足Lipschitz条件。如果单步 法与微分方程初值问题相容，则单步法收敛。
2维的龙格库塔算法
2维的龙格库塔算法
二维龙格-库塔算法是一种数值计算方法，用于解决微分方程的数值解。
它使用龙格-库塔方法（Runge-Kutta methods）在二维空间中迭代地计算方程的解。
以下是二维龙格-库塔算法的基本步骤：
1.定义微分方程：给定一个微分方程 dy/dx = f(x, y)，其中 y 是未知函数，x 是
自变量。
2.初始化参数：选择初始值 x0 和 y0，以及步长 h 和迭代的最大步数 n。
3.初始化状态向量：创建一个大小为 (n+1)×2 的数组，用于存储 x 和 y 的值。
4.迭代计算：
a. 计算 k1 = h f(x0, y0)
b. 计算 k2 = h f(x0 + h/2, y0 + k1/2)
c. 计算 k3 = h f(x0 + h/2, y0 + k2/2)
d. 计算 k4 = h f(x0 + h, y0 + k3)
e. 更新 y0 = y0 + (k1 + 2k2 + 2k3 + k4) / 6
f. 更新 x0 = x0 + h
g. 如果 x0 小于或等于 n*h，则跳转到步骤 4，否则算法结束。
5. 输出结果：输出 x 和 y 的值。
以上是二维龙格-库塔算法的基本步骤，具体的实现可能会根据不同的微分方程和需求而有所不同。
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方法计算初值问题
Runge-Kutta积分方法(共24张PPT)
0.8 2.014457 1.914603 2.014459 1.68E-6 9.99E-2
0.9 2.150395 2.021086 2.150397 1.96E-6 1.29E-1
1.0 2.281716 2.117800 2.281718 2.32E-6 1.64E-1
第十二页，共二十四页。
结果(jiē 及 guǒ) 比较
将展开式代入,得到局部截断误差：
dn1 (1 c1 c2 )hFn
1 2
c22
h
2
Ft(tn
,
Yn
)
1 2
c2 21
h2
FY
(t
n
, Yn
)Fn
O(h3 ).
要使得(shǐ de)方法是二阶的，则局部截断误差应该为三阶小量，即：
1 c1 c2 0;
1 2
c22
0;
1 2
c2 21
分别用二阶四阶步长都取为四阶二阶四阶误差二阶误差0010000001000000100000000000000000001110482911024501104829160e7238e302121859712115071218597340e7709e303134014113257661340141548e7144e204146817514436711468175769e7245e205160127815635061601279995e7378e206173788016833741737881120e6545e207187624618011791876247142e6751e208201445719146032014459168e6999e209215039520210862150397196e6129e110228171621178002281718232e6164e1结果及比较结果及比较10451060101040时误差为而二阶公式相对误差仅为仍然是相当精确的结果时误差为时误差为对四阶公式rungekuttarungekuttarungekutta类似前面的推导可以导出更高阶的公式
二阶runge-kutta法介数证明
二阶runge-kutta法介数证明
二阶Runge-Kutta法是一种常用的数值解常微分方程的方法。
该方法基于以下假设：
1. 求解的常微分方程为一阶方程：dy/dt = f(t, y)。
2. 采用离散化的方式对时间轴进行分割，将时间区间[tn, tn+1]划分为n个子区间，每个子区间的长度为h = (tn+1 - tn)/n。
3. 求解的方程在每个子区间上取常值，即在子区间[tn, tn+1]上可以近似为：y(t) ≈ ytn。
4. 采用迭代的方式，首先将当前y的值初始化为y0，然后通过迭代得到下个时间点的近似值。
5. 第二次迭代的近似值将作为下个时间点的初始值。
基于以上假设，二阶Runge-Kutta法可以描述为以下步骤：
Step 1: 给定初始条件 y0 和迭代时间点 t0，设定步长 h。
Step 2: 计算 k1 = hf(t0, y0)。
Step 3: 计算 k2 = hf(t0 + h/2, y0 + k1/2)。
Step 4: 计算下个时间点的近似值 y1 = y0 + k2。
Step 5: 更新时间点和迭代值，即 t1 = t0 + h，并取 y0 = y1。
Step 6: 重复步骤2到步骤5，直到达到所需的时间点。
通过以上的步骤，二阶Runge-Kutta法可以较好地逼近求解的常微分方程。
其中，k1为在当前时间点t0处的斜率值，k2则是在下个时间点t0 + h/2处的斜率值。
该方法利用了斜率在当前时间点和下个时间点之间的平均变化率来逼近方程的解。
一高阶显式单步法的构造方法
k2
f (xn
h 2 , yn
h 2
k1
)
(
yn
h 2 k1)
(
2h
2 ) yn
k3
f (xn
h, 2
yn
h 2
k2
)
( yn
h 2 k2)
(
2h
2
3h2
4 ) yn
k4 f ( xn h, yn hk3 ) ( yn hk3 ) yn hk3
yn1
(1
h
2h2
k1 f ( xn , yn )
k2 f ( xn ph, yn phk1) f ( xn , yn ) phfx phy ' f y O(h2 ) f ( xn , yn ) phf '( xn , yn ) O(h2 )
f ( xn , yn ) phy "( xn ) O(h2 )
2!
3h3
3!
4h4
4!
) yn
特征方程旳根是 1 h 2h2 3h3 4h4
2! 3! 4!
当：
1
经典R－K格式是条件稳定旳，稳定条件是：
h 2.78
可见，显式公式一定是条件稳定旳； 绝对稳定旳公式，一定是隐式旳。
因为单步法是 p阶旳：h0 , 0 h h0 满足| Tn1 | Chp1
| En1 || En | hL | En | Ch p1 | En |
其中 1 hL, Ch p1
En O(hp )
| En | | En1 | 2 | En2 |
3 | En3 | (1 2 )
计算过程中产生旳舍入误差对计算成果旳影响
首先以Euler公式为例，来讨论一下舍入误差旳传播:
runge-kutta方法
runge-kutta方法Runge-Kutta方法是求解常微分方程组的一种常用数值方法，它是由德国数学家卡尔·龙格和马丁·库塔（Martin Kutta）发明的。
该方法是一种高阶精度的方法，可以求解一阶或高阶的常微分方程组。
$$ y_{n+1} = y_n + \sum_{i=1}^s b_i k_i $$$y_n$表示第$n$个步骤中的输出，$y_{n+1}$是下一个步骤的输出。
$k_i$表示第$i$个中间变量，$b_i$是权值，它们是根据所采用的方法计算得出的。
我们将在下文中详细介绍这些变量和权值的计算方法。
在数值计算中，我们通常采用自适应步长的方法，这需要我们给出一个误差容限值$\epsilon$，根据误差限定步长大小。
如果当前步长不满足误差要求，我们会将步长缩小再进行计算，或者选择更高阶的方法进行计算，以提高精度。
接下来我们将详细介绍如何使用Runge-Kutta方法求解常微分方程组。
1. Runge-Kutta方法的中间变量$k_i$的计算我们需要计算中间变量$k_i$的值。
它们的计算公式如下：$$ k_1 = hf(t_n, y_n) $$$$ k_2 = hf(t_n + c_2 h, y_n + a_{21}k_1) $$$$ ... $$$$ k_s = hf(t_n + c_s h, y_n + \sum_{j=1}^{s-1}a_{sj}k_j) $$$h$表示步长大小，$t_n$和$y_n$分别表示当前时间和状态，$c_i$和$a_{ij}$是与所采用的方法相关的常数。
$$ c_2 = \frac{1}{2} , a_{21} = \frac{1}{2}$$- 改进的Runge-Kutta方法（RK4）：根据具体选择的方法，我们可以计算出$k_i$的值。
$b_{ij}$也是与所采用的方法相关的常数，对于不同的方法而言，它们的计算公式不同。
根据计算出的中间变量$k_i$和权值$b_i$，我们可以进行Runge-Kutta方法的迭代，求解常微分方程组。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
h h ( x, y, h) f ( x, y ) F ( x, y ) (G f y F )( x , y ( x )) y ) 二级方法：N =2 k2 f ( x a2h, y hb21k1 ) f ( x, y) [a2hf x hb21k1 f y ]
2
h y( x ) 若取 ( x , y , h) y( x ) 2!
则
h p1 ( p ) y ( x) p!
d y( x ) f ( x , y( x )) y( x )的各阶导数的计算 dx 2 d y( x ) f ( x , y( x )) y( x ) 2 f ( x , y( x )) ...... dx
§2 龙格-库塔（Runge-Kutta）法
一、高阶显式单步法的构造方法
显式单步法的一般形式： yn1 yn h ( xn , yn , h) 构造高阶方法，即如何确定增量函数 ( x, y, h) 使得该方法的局部截断误差的阶数尽量高。
y( xn1 ) yn1 O(h
1、Taylor级数方法
四级方法：N =4
常见的2种四阶方法：
局部截断误差 O(h5 )
经典Runge-Kutta方法 Kutta四阶方法（见教材） h yn1 yn (k1 2k2 2k3 k4 ) 6 k1 f ( xn , yn )
h h k2 f ( xn , yn k1 ) 2 2 h h k3 f ( xn , yn k2 ) 2 2
1 3 2 c1 , c2 , a2 4 4 3
h 2 2h yn1 yn [ f ( xn , yn ) 3 f ( xn h, yn f ( xn , yn ))] 4 3 3
三级方法：N =3
类似于N =2的推导方法，可得到
1 c1 c2 c3 1; c2 a2 c3a3 ; 2 1 1 2 2 c2 a2 c3a3 ; c3a2b32 6 3
常见的2种三阶方法：
O(h4 )
Kutta三阶方法
Heun三阶方法（见教材）
h h k1 f ( xn , yn ); k2 f ( xn , yn k1 ) 2 2
h yn1 yn ( k1 4k2 k3 ) 6
k3 f ( xn h, yn hk1 2hk2 )
此时将步长折半一次计算，得到的为最终结果。
§3 单步法的收敛性、相容性和绝对稳定性
一、收敛性 /*Convergence*/
dy f ( x, y ) Def 1 对于初值问题 dx () 的一种 y( x0 ) y0 ; x x0
单步法 yn1 yn h ( xn , yn , h) 产生的近似解，如果
p 1
) p 为尽可能大的整数
dy f ( x, y ) 设 y( x )满足初值问题： dx ( ) y( x0 ) y0 ; x x0
h y( xn1 ) yn1 y( xn ) hy( xn ) y( xn ) 2 ! h p ( p) p 1 y ( xn ) O( h ) yn h ( xn , yn , h) p!
h h2 ( x, y, h) f F (G f y F ) 2 3!
比较关于 h 的幂次相同项的系数，得到方程组：
c1 c2 1 1 方程组有无穷多解：二级方法有无穷多种 c2 a 2 2
常见的3种二级方法：
O(h )
3
1 中点法（修正的Euler法） 取 c1 0, c2 1, a2 2 h h
j 1
N
(1)一级二阶的隐式中点方法：
yn1 yn hk1
hk1 h k1 f ( xn , yn ) 2 2
(2)二级四阶的隐式R-K方法：
h( k1 k2 ) yn1 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
三、变步长方法 基本思想：根据精度自动地选择步长 对于经典Runge-Kutta方法： n 0,1, 2,
Step1:设从
xn 出发，以 h 为步长，经过一步计算得到
y( xn1 ) y
( h) n1
Ch
5
h Step2:取 为步长，再从 xn 出发，经过两步计算得到 2
y( xn1 ) y
前式中参数的求解方法
利用二元函数的Taylor展开公式：
1 i f ( x h, y k ) (h k ) f ( x , y ) x y i 0 i !
将每个 ki 展开，然后代入 yn1 yn h ( xn , yn , h)
与Taylor级数法中的增量函数比较 h幂次相同项的系数 从而得到关于 ai , bij , ci 的方程组，解之即得。
| ( x, y1 , h) ( x, y2 , h) | L | y1 y2 | y1 , y2
p y ( x ) y O ( h ) 则该单步法是收敛的，且 n n
证明： 记 en
y( xn ) yn
由截断误差的定义
y( xn1 ) y( xn ) h ( xn , y( xn ), h) Tn1
yn y ( x n ) ， 对于任一固定的 xn x0 nh ，均有 lim h 0
则称该单步法是收敛的。
类似地可以定义隐式单步法、多步法（§4）的收敛性
Th3.1 设初值问题（*）对应的下列单步法是 p阶的，
yn1 yn h ( xn , yn , h)
且函数 满足对 y 的Lipschitz条件，即存在常数 L 0
Heun（休恩）二级方法
要求 h 项的系数尽量相同
2 2 2 2
2
cah 3 ( x , y, h) (c1 c2 ) f c2a2 hF G O( h ) 2 2 h h ( x, y, h) f F (G f y F ) 2 3! 1 1 2 得到方程组： c1 c2 1; c2 a2 ; c2 a 2 2 3
0.1 0.2 0.3 0.4 0.5 1.0954 1.1832 1.2649 1.3416 1.4142 0.6 0.7 0.8 0.9 1.0 1.4832 1.5492 1.6125 1.6733 1.7321
二、高阶和隐式Runge-Kutta方法 注：对于显式N级R-K方法，最多只能得到N级方法； 已经证明N级R-K方法的阶 p( N )具有下列关系： N
( x, y, h) f ( x, y( x))
即为Euler方法
2、Runge-Kutta方法 基本思想：利用 f ( x, y ) 在某些特殊点上的函数值的线性 组合来构造高阶单步法的增量函数 N级（阶）Runge-Kutta方法的一般形式：
yn1 yn h ( xn , yn , h)
k4 f ( xn h, yn hk3 )
2x 例2：用经典的Runge-Kutta方法 dy y y x (0,1) 求解下列初值问题 h 0.1 。 dx y ( 0) 1 解： 经典的四阶Runge-Kutta公式：
h yn1 yn (k1 2k2 2k3 k4 ) 6 h 2 ( x ) n 2 xn h 2 k1 yn ; k2 yn k1 h yn 2 yn k1 2 h 2 ( x ) n h 2( xn h) 2 k3 yn k2 ; k4 yn hk3 h 2 yn hk3 yn k 2 2
en1 en h[ ( xn , y( xn ), h) ( xn , yn , h)] Tn1
其中 ( x , y , h)
c k
i 1 i
N
i
i 1
k1 f ( x, y )
j 1
ki f ( x ai h, y h bij k j ), i 2, 3,..., N
ai bij
j 1
i 1
N =1：Euler方法
当N >1时，适当选取式中参数，使该方法的阶数尽量高
xn yn xn yn
0.1 0.2 0.3 0.4 0.5 1.0954 1.1832 1.2649 1.3416 1.4142 0.6 0.7 0.8 0.9 1.0 1.4832 1.5492 1.6125 1.6733 1.7321
同保留5位的精确值完全一致： y
2x 1
xn yn xn yn
p( N )
1,2,3,4 N
5,6,7 N-1
8,9 N-2
10,11,…
N 2
若要得到N阶以上方法，则使用N级隐式R-K方法
N级隐式R-K方法的一般形式：
yn1 yn h ci ki
i 1 N
N级隐式R-K法 可以达到2N阶
ki f ( xn ai h, yn h bij k j ); i 1,..., N
2
G( x, y( x)) f y ( x, y( x))F ( x, y( x))
2
以此类推，可以求出y( x )的各阶导数，则有
h h ( x, y, h) f ( x, y ) F ( x, y ) (G f y F )( x , y ( x )) 2! 3!
1 [(a2 h)2 f xx 2a2b21h2 k1 f xy (b21h)2 k12 f yy ] O(h3 ) 2 f ( x, y) a2h[ f x ff y ]