（2.2） q3 ?cph[fx(xn,yn)?fy(xn,yn)f(xn,yn)]?O(h).p22yn?1?y(xn)?(c1?c2)hf(xn,yn)再将y(xn?1)在x?xn作Taylor展开
h2y(xn?1)?y(xn)?hy?(xn)?[fx(xn,yn)?fy(xn,yn)f(xn,yn)]?O(h3). (2.3)
2!若要求局部截断误差达到O(h3)，则通过比较上面两式知，参数c1,c2,p,q必须满足
11c1?c2?1,c2p?,c2q?
221这是四个未知量，三个方程式的方程组，因此解不惟一。当我们选取c1?c2?,p?q?1时，
2(2.2)和(2,.3)的前三项完全一致。 将参数值代入(2.1)便得到二阶Runge-Kutta方法 (简记为R-K方法)的一种迭代方法
hh?y?y?k?k2,n1?n?122? (2.4) ?k1?f(xn,yn),?k?f(x?h,y?hk),nn1?2?它的局部截断误差为O(h2)，(2.4)实际上就是改进的Euler方法(1.15).
若取c1?0,c2?1,p?q?1,则得到二阶R-K方法的又一种迭代公式 2??yn?1?yn?hk2,? ?k1?f(xn,yn),?hh?k2?f(xn?,yn?k1),?22或写成
h1yn?1?yn?hf(xn?,yn?hf(xn,yn)). (2.5)
22公式(2.4)、(2.5)都是显式单步二阶公式.
可以证明了不论这四个参数如何选择，都不能使局部截断误差达到O(h). 这说明在计算两次函数值的情况下，局部截断误差的阶最高为O(h)，要再提高阶就必须增加计算函数值的次数。
34二、三阶及四阶Runge-Kutta方法
为了提高方法的精度,考虑每步计算三次函数f(x,y)值。根据两次计算函数值的做法，很自然地取yn?1的形式为
9
?yn?1?yn?h(c1k1?c2k2?c3k3),?k?f(x,y),?1nn ??k2?f(xn?a1h,yn?b21hk1),??k3?f(xn?a2h,yn?b31hk1?b32hk2).适当选择参数c1、c2、c3、a1、a2、a3、b21、b31、b32，使截断误差y(xn?1)?yn?1的阶尽可能
高。由于在推导公式时只考虑局部截断误差，故设yn?y(xn)。类似前面二阶R-K公式的推导方法，将yn?1在(xn,yn)处作Taylor展开，然后再将y(xn?1)在x?xn处作Taylor展开（这里不详细写出展开式），只要两个展开式的前四项相同，便有y(xn?1)?yn?1?0(h4)，而要两个展开式的前四项相同，参数必须满足
c1?c2?c3?1,a1?b21,a3?b31?b32,
111 (2.6) 22a1c2?a2c3?，a1c2?a2c3?,a1c3b32?.236这是8个未知数6个方程的方程组，解不是惟一的，可以得到很多公式。当我们取方程组的一组解为a1?11141,a2?1,b21?,b31??1,b32?2,c1?,c2?,c3?时，就得到一22666种常用的三阶Runge-Kutta公式
h?y?y?(k1?4k2?k3),n?n?16??k1?f(xn,yn), (2.7) ??k?f(x?h,y?hk),nn1?222?k?f(x?h,y?hk?2hk).nn12?3从上面导出公式(2.7)的过程知，它的局部截断误差为O(h).
如果每步计算四次函数f(x,y)的值，完全类似的，可以导出局部截断误差为O(h5)的四阶Runge-Kutta公式。详细推导过程这里略去，只给出结果
41?y?y?(k1?2k2?2k3?k4),n?n?16??k1?f(xn,yn),?hh? (2.8) ?k2?f(xn?,yn?k1),22?1h?k?f(x?h,y?k2),nn?322???k4?f(xn?h,yn?hk3).公式(2.8)就是一种常用的四阶Runge-Kutta公式，也称为标准的(或经典的)四阶R-K方法。
三阶、四阶Runge-Kutta公式都是单步显式公式。
10
例2 试用Euler方法、改进的Euler方法及四阶经典R-K方法在不同步长下计算初值问题
?dy???y(1?xy),0?x?1, ?dx??y(0)?1在0.2、0.4、0.8、1.0处的近似值，并比较它们的数值结果.
解 对上述三种方法,每执行一步所需计算f(x,y)??y(1?xy)的次数分别为1、2、4。为了公正起见，上述三种方法的步长之此应为1:2:4。因此，在用Euler方法、改进的Euler方法及四阶经典R-K方法计算0。2、0。4、0。8、1。0处的近似值时，它们的步长应分别取为0。05、0。1、0。2，以使三种方法的计算量大致相等。
Euler方法的计算格式为
yn?1?yn?0.05?[yn(1?xnyn)].
改进的Eluer方法的计算格式为
??yp?yn?0.1?[yn(1?xnyn)],???yc?yn?0.1?[yp(1?xn?1yp)], ??yn?1?1(yp?yc).?2?四阶经典R-K方法的计算格式为
0.2?y?y??(k1?2k2?2k3?k4),n?1n?6??k1??yn(1?xnyn),?0.20.20.2?k??(y??k)[1?(x?)(y??k1)], ?2n1nn222?0.20.20.2?k??(y??k)[1?(x?)(y??k2)],n2nn?3222???k4??(yn?0.2?k3)[1?(xn?0.2)(yn?0.2?k3)]初始值均为y0?y(0)?1，将计算结果列于表9.2.
表9.2 Euler方法 （步长h=0.05） 改进的Euler方法 （步长h=0.1） 四阶经典R-K方法 （步长h=0.2） 准确解 xn 0.2 0.4 0.6 yn 0.8031866 0.6271777 0.4825586 yn 0.8052632 0.6325651 0.4905510 yn 0.8046363 0.6314653 0.4891979 y(xn) 0.8046311 0.6314529 0.4891800 11
0.8 1.0 0.3693036 0.2827482 0.3786397 0.2923593 0.3772249 0.2910086 0.3772045 0.2909884 从表9.2可以看出，在计算量大致相等的情况下，Euler方法计算的结果只有2位有效数字，改进的Euler方法计算的结果有3位有效数字，而四阶经典R-K方法计算的结果却有5位有效数字，这与理论分析是一致的。例1和例2的计算结果说明，在解决实际问题时，选择恰当的算法是非常必要的。
需要指出的是Runge-Kutta方法的基于Taylor展开法，因而要求解具有足够的光滑性。如果解的光滑性差，使用四阶Runge-Kutta方法求得数值解的精度，可能不如改进的Euler方法精度高。因此，在实际计算时，要根据具体问题的特性，选择合适的算法。
三、变步长的Runge-Kutta方法
上面导出的Runge-Kutta方法都是定步长的，单从每一步来看，步长h越小，局部截断误差也越少，但随着步长的减小，在一定范围内要进行的步数就会增加，而步数增加不仅增加计算量，还有可能相起舍入误差的积累过大。由于Runge-Kutta方法是单步法，每一步计算步长都是独立的，所以步长的选择具有较大的灵活性。因此根据实际问题的具体情况合理选择每一步的步长是非常有意义的。正如第五章建立变长的Simpson求积公式那样，我们来建立变步长的Runge-Kutta公式。
以四阶标准Runge-Kutta公式为例进行说明，从基点x0出发，先选一个步长h，利用公式
(h)(2.8)求出的近似值记为y1, 由于公式的局部截断误差为O(h5), 所以有
y(x1)?y1(k)?ch2, (2.9)
x?x1h其中c为常数，然后将步长h进行折半, 即取步长为, 由基点x0出发, 到0算一步, 由
22h()x0?x1的x1再算一步, 将求得的近似值记为y12, 因此有 2y(x1)?y1?h2?由(2.9)及(2.10)得
?h??2c??. (2.10)
?2?5y(x1)?y1??1 ?y(x1)?y1(h)16.h2一般地, 从xn出发, 照上面的做法也可得到
??y(xn?1)?yn1?1?, (h)y(xn?1)?yn?116h2或写成
12