Euler方法有各种格式，但其精度最高不超过2阶，一般难以满足实际计算的精度要求。因此，有必要构造精度更高的数值计算公式求解微分方程。Runge-Kutta方法就是一种高精度的经典的解常微分方程的单步方法。
1. Runge-Kutta 方法的基本思想
对于函数 y = y ( x ) y=y(x) y=y(x)，设 y n = y ( x n ) y_n=y(x_n) yn?=y(xn?)，将 y ( x n + 1 ) y(x_{n+1}) y(xn+1?)在点 x n x_n xn?处展开为Taylor级数：
y ( x n + 1 ) = y ( x n ) + h y ′ ( ? ) x n < ? < x n + 1 y(x_{n+1})=y(x_n)+hy'(\epsilon) \quad x_n<\epsilon<x_{n+1} y(xn+1?)=y(xn?)+hy′(?)xn?<?<xn+1?
利用 y ′ ( x n ) = f ( x n , y n ) y'(x_n)=f(x_n,y_n) y′(xn?)=f(xn?,yn?)可得：
y ( x n + 1 ) = y ( x n ) + h f ( ? , y ( ? ) ) y(x_{n+1})=y(x_n)+hf(\epsilon,y(\epsilon)) y(xn+1?)=y(xn?)+hf(?,y(?))
其中， f ( ? , y ( ? ) ) f(\epsilon,y(\epsilon)) f(?,y(?))为区间 ( x n , x n + 1 ) (x_n,x_{n+1}) (xn?,xn+1?)上的平均斜率，记为 k ? k^* k?。因此，只要能设法提供一种算法求得 k ? k^* k?，就可相应地导出一种计算格式。
回顾一下前面所讲的Euler方法可知：
若取 x n x_n xn?点处的斜率 k 1 = f ( x n , y n ) k_1=f(x_n,y_n) k1?=f(xn?,yn?)作为 k ? k^* k?，则可得到显示Euler格式，但只有一阶精度；
若取 x n x_n xn?和 x n + 1 x_{n+1} xn+1?两点处的斜率 k 1 = f ( x n , y n ) k_1=f(x_n,y_n) k1?=f(xn?,yn?)和 k 2 = f ( x n + 1 , y n + h k 1 ) k_2=f(x_{n+1},y_{n}+hk_1) k2?=f(xn+1?,yn?+hk1?)的算术平均值作为平均斜率 k o k^o ko，即 k o = ( k 1 + k 2 ) / 2 k^o=(k_1+k_2)/2 ko=(k1?+k2?)/2，则得到改进的Euler格式，则有二阶精度。
由此推而广之，如果设法在区间 ( x n , x n + 1 ) (x_n,x_{n+1}) (xn?,xn+1?)内取若干个点的斜率，然后将它们加权平均作为平均斜率 k o k^o ko，则可以构造一种具有更高精度的计算格式。这就是Runge-Kutta方法的基本思想。
2. 二阶Runge-Kutta格式
在改进的Euler格式的基础上加以修改，即在区间 [ x n , x n + 1 ] [x_n,x_{n+1}] [xn?,xn+1?]内，除 x n x_n xn?外再取区间中任一点p，则有：
x n + p = x n + p h 0 < p ≤ 1 x_{n+p}=x_n+ph \quad 0<p\leq 1 xn+p?=xn?+ph0<p≤1
取 x n x_n xn?和 x n + p x_{n+p} xn+p?两点处的斜率 k 1 k_1 k1?和 k 2 k_2 k2?加权平均作为平均斜率 k o k^o ko，即：
k o = ( 1 ? λ ) k 1 + λ k 2 λ 为 待 定 系 数 k^o=(1-\lambda)k_1+\lambda k_2 \quad \lambda 为待定系数 ko=(1?λ)k1?+λk2?λ为待定系数
这样构造出的计算格式为：
y n + 1 = y n + h [ ( 1 ? λ ) k 1 + λ k 2 ] (1a) y_{n+1}=y_n+h[(1-\lambda)k_1+\lambda k_2] \tag{1a} yn+1?=yn?+h[(1?λ)k1?+λk2?](1a)
其中：
k 1 = f ( x n , y n ) ( 1 b ) k 2 = f ( x n + p , y n + p h k 1 ) ( 1 c ) k_1=f(x_n,y_n) \quad (1b)\\ k_2=f(x_{n+p},y_n+phk_1) \quad (1c) k1?=f(xn?,yn?)(1b)k2?=f(xn+p?,yn?+phk1?)(1c)
（1）式称为二阶的Runge-Kutta格式，其中 p , λ p,\lambda p,λ为待定参数。在确定 p , λ p,\lambda p,λ的同时，使二阶Runge-Kutta格式具有较高的精度，具体推导如下：
根据假设 y n = y ( x n ) y_n=y(x_n) yn?=y(xn?)，有：
k 1 = f ( x n , y n ) = y ′ ( x n ) k 2 = f ( x n + p , y n + p h k 1 ) = f ( x n + p h , y n + p h k 1 ) k_1=f(x_n,y_n)=y'(x_n) \\ k_2=f(x_{n+p},y_n+phk_1)=f(x_n+ph,y_n+phk_1) k1?=f(xn?,yn?)=y′(xn?)k2?=f(xn+p?,yn?+phk1?)=f(xn?+ph,yn?+phk1?)
将 k 2 k_2 k2?右端在点 ( x n , y n ) (x_n,y_n) (xn?,yn?)处展开为Taylor级数：
k 2 = f ( x n , y n ) + p h ? f x ( x n , y n ) + p h y ′ ( x n ) ? f y ( x n , y n ) + O ( h 2 ) k_2=f(x_n,y_n)+ph·f_x(x_n,y_n)+phy'(x_n)·f_y(x_n,y_n)+O(h^2) k2?=f(xn?,yn?)+ph?fx?(xn?,yn?)+phy′(xn?)?fy?(xn?,yn?)+O(h2)
而
p h ? f x ( x n , y n ) + p h ? y ′ ( x n ) ? f y ( x n , y n ) = p h ? f ′ ( x n , y n ) = p h ? y ′ ′ ( x n ) ph·f_x(x_n,y_n)+ph· y'(x_n)·f_y(x_n,y_n)=ph·f'(x_n,y_n)=ph·y''(x_n) ph?fx?(xn?,yn?)+ph?y′(xn?)?fy?(xn?,yn?)=ph?f′(xn?,yn?)=ph?y′′(xn?)
故
k 2 = f ( x n , y n ) + p h y ′ ′ ( x n ) + O ( h 2 ) = y ′ ( x n ) + p h y ′ ′ ( x n ) + O ( h 2 ) k_2=f(x_n,y_n)+phy''(x_n)+O(h^2)=y'(x_n)+phy''(x_n)+O(h^2) k2?=f(xn?,yn?)+phy′′(xn?)+O(h2)=y′(xn?)+phy′′(xn?)+O(h2)
将 k 1 k_1 k1?和 k 2 k_2 k2?代入 y n + 1 = y n + h [ ( 1 ? λ ) k 1 + λ k 2 ] y_{n+1}=y_n+h[(1-\lambda)k_1+\lambda k_2] yn+1?=yn?+h[(1?λ)k1?+λk2?]中，有：
y n + 1 = y n + h [ ( 1 ? λ ) k 1 + λ k 2 ] = y ( x n ) + h y ′ ( x n ) + λ p h 2 y ′ ′ ( x n ) + λ O ( h 3 ) y_{n+1}=y_n+h[(1-\lambda)k_1+\lambda k_2] \\ =y(x_n)+hy'(x_n)+\lambda ph^2y''(x_n)+\lambda O(h^3) yn+1?=yn?+h[(1?λ)k1?+λk2?]=y(xn?)+hy′(xn?)+λph2y′′(xn?)+λO(h3)
再将 y ( x n + 1 ) y(x_{n+1}) y(xn+1?)在点 ( x n , y n ) (x_n,y_n) (xn?,yn?)处展开为Taylor级数：
y ( x n + 1 ) = y ( x n + h ) = y ( x n ) + h y ′ ( x n ) + 1 2 h 2 y ′ ′ ( x n ) + O ( h 3 ) y(x_{n+1})=y(x_n+h)=y(x_n)+hy'(x_n)+\frac{1}{2}h^2y''(x_n)+O(h^3) y(xn+1?)=y(xn?+h)=y(xn?)+hy′(xn?)+21?h2y′′(xn?)+O(h3)
比较两式，得：
y ( x n + 1 ) ? y n + 1 = ( 1 2 h 2 ? λ p h 2 ) y ′ ′ ( x n ) ? λ O ( h 3 ) y(x_{n+1})-y_{n+1}=(\frac{1}{2}h^2-\lambda ph^2)y''(x_n)-\lambda O(h^3) y(xn+1?)?yn+1?=(21?h2?λph2)y′′(xn?)?λO(h3)
因此，当
λ p = 1 2 \lambda p=\frac{1}{2} λp=21?
时，二阶Runge-Kutta格式具有为2阶精度。可见二阶Runge-Kutta格式有很多具体形式：
若取 λ = 1 2 , p = 1 \lambda=\frac{1}{2},p=1 λ=21?,p=1，则可得改进的Euler格式；
若取 λ = 3 4 , p = 2 3 \lambda=\frac{3}{4},p=\frac{2}{3} λ=43?,p=32?，则可得Heun（休恩）格式；
若取 λ = 1 , p = 1 2 \lambda=1,p=\frac{1}{2} λ=1,p=21?，则可得变形的Euler格式，即中点格式：
{ y n + 1 = y n + h k 2 k 1 = f ( x n , y n ) k 2 = f ( x n + 1 2 , y n + h 2 k 1 ) \begin{cases} y_{n+1}=y_n+hk_2 \\ k_1 = f(x_n,y_n) \\ k_2 = f(x_{n+\frac{1}{2}},y_n+\frac{h}{2}k_1) \end{cases} ??????yn+1?=yn?+hk2?k1?=f(xn?,yn?)k2?=f(xn+21??,yn?+2h?k1?)?
3. 四阶Runge-Kutta格式
在区间 [ x n , x n + 1 ] [x_n,x_{n+1}] [xn?,xn+1?]内，使用两个不同的点可以构造出二阶Runge-Kutta格式。依此规律，在区间 [ x n , x n + 1 ] [x_n,x_{n+1}] [xn?,xn+1?]内，取3个不同的点可以构造出三阶Runge-Kutta格式；取4个不同的点可以构造出四阶Runge-Kutta格式。对此可以加以证明。在实际中，应用最广泛的是四阶经典的Runge-Kutta格式：
{ y n + 1 = y n + h 6 ( k 1 + 2 k 2 + 2 k 3 + k 4 ) k 1 = f ( x n , y n ) k 2 = f ( x n + 1 2 , y n + h 2 k 1 ) k 3 = f ( x n + 1 2 , y n + h 2 k 2 ) k 4 = f ( x n + 1 , y n + h k 3 ) \begin{cases} y_{n+1}=y_n+\frac{h}{6}(k_1+2k_2+2k_3+k_4) \\ k_1=f(x_n,y_n) \\ k_2=f(x_{n+\frac{1}{2}},y_n+\frac{h}{2}k_1) \\ k_3=f(x_{n+\frac{1}{2}},y_n+\frac{h}{2}k_2) \\ k_4=f(x_{n+1},y_n+hk_3) \end{cases} ????????????????yn+1?=yn?+6h?(k1?+2k2?+2k3?+k4?)k1?=f(xn?,yn?)k2?=f(xn+21??,yn?+2h?k1?)k3?=f(xn+21??,yn?+2h?k2?)k4?=f(xn+1?,yn?+hk3?)?
四阶Runge-Kutta方法的优点如下：
（1）它是一种高精度的单步法，可达四阶精度 O ( h 5 ) O(h^5) O(h5)；
（2）数值稳定性较好；
（3）只需知道一阶导数，无需明确定义或计算其他高阶导数；
（4）只需给出 y n y_n yn?就能计算出 y n + 1 y_{n+1} yn+1?，所以能够自启动；
（5）编程容易。
四阶Runge-Kutta法除经典的格式外，还有其他的格式。例如：
（1）Kutta格式
{ y n + 1 = y n + h 8 ( k 1 + 3 k 2 + 3 k 3 + k 4 ) k 1 = f ( x n , y n ) k 2 = f ( x n + 1 3 , y n + h 3 k 1 ) k 3 = f ( x n + 2 3 , y n ? h 3 k 2 + h k 2 ) k 4 = f ( x n + 1 , y n + h k 1 ? h k 2 + h k 3 ) \begin{cases} y_{n+1}=y_n+\frac{h}{8}(k_1+3k_2+3k_3+k_4) \\ k_1=f(x_n,y_n) \\ k_2=f(x_{n+\frac{1}{3}},y_n+\frac{h}{3}k_1) \\ k_3=f(x_n+\frac{2}{3},y_n-\frac{h}{3}k_2+hk_2) \\ k_4=f(x_{n+1},y_n+hk_1-hk_2+hk_3) \end{cases} ????????????????yn+1?=yn?+8h?(k1?+3k2?+3k3?+k4?)k1?=f(xn?,yn?)k2?=f(xn+31??,yn?+3h?k1?)k3?=f(xn?+32?,yn??3h?k2?+hk2?)k4?=f(xn+1?,yn?+hk1??hk2?+hk3?)?
（2）Gill格式：
{ y n + 1 = y n + h 6 ( k 1 + ( 2 ? 2 ) k 2 + ( 2 + 2 ) k 3 + k 4 ) k 1 = f ( x n , y n ) k 2 = f ( x n + 1 2 , y n + 1 2 h k 1 ) k 3 = f ( x n + 1 2 , y n + 2 ? 1 2 h k 1 + 2 ? 2 2 h k 2 ) k 4 = f ( x n + 1 , y n ? 2 2 h k 2 + 2 + 2 2 h k 3 ) \begin{cases} y_{n+1}=y_n+\frac{h}{6}(k_1+(2-\sqrt{2})k_2+(2+\sqrt{2})k_3+k_4) \\ k_1=f(x_n,y_n) \\ k_2=f(x_{n+\frac{1}{2}},y_n+\frac{1}{2}hk_1) \\ k_3=f(x_{n+\frac{1}{2}},y_n+\frac{\sqrt{2}-1}{2}hk_1+\frac{2-\sqrt{2}}{2}hk_2) \\ k_4=f(x_{n+1},y_n-\frac{\sqrt{2}}{2}hk_2+\frac{2+\sqrt{2}}{2}hk_3) \end{cases} ????????????????yn+1?=yn?+6h?(k1?+(2?2?)k2?+(2+2?)k3?+k4?)k1?=f(xn?,yn?)k2?=f(xn+21??,yn?+21?hk1?)k3?=f(xn+21??,yn?+22??1?hk1?+22?2??hk2?)k4?=f(xn+1?,yn??22??hk2?+22+2??hk3?)?