世界因皒而精彩
111(4)4yn?1?yn?yn?h?yn??h2?yn???h3?ynh2624 （21）
相比较，有
c1?c2?c3?c4?1,c2a2?c3a3?c4a4?c2a22?c3a32?c4a42?11,c2a23?c3a33?c4a43?,34
1,2
又由于Runge—Kutta公式(17)中既要反映x的变化，又要反映y的变化，所以简单的选取f(x,y)为一对线性函数，这样将(13)式用于
y??xy,y?0??0 （22）
得
k1?xnyn
'''2K2?K1?a2ynh?a2K1h2
''22''2K3?K1?a3ynh?(a2b32y''xn?a3K1)h2?(a2b32K1xn?a2a3b32yn)h3?a2a3b32K1h4, ''2K4?K1?a4ynh?[(a2b42?a3b43)y''xn?a4K1]h2?2''2''[(a2b42K1?a2b32b42ynxn?b43a3K1)xn?a4(a2b42?a3b43)yn]h3, （23）
yn?1?yn?h(c1K1?c2K2?c3K3?c4K4)
'?yn(c1?c2?c3?c4)yn?(c2a2?c3a3?c4a4)y''h2
22?[(c2a22?c3a)4n'y?(c?3?c4a3a2b32''?{[c3a2a3b32?c4a4(a2b42?a3b43)]yn?2(c3a2b3?22c42a2b?ca)3b424n'4n3c4a2?b42''c4a3)bb]hn3n243yxy?x''2c4a2bb}xn32n4 y2h对（22）式求导数有
4yn????2yn??xnyn??,yn???3yn???xn2yn???2xnyn?,
将上两式代入(23)式并与Taylor公式(21)比较有：
c3a2b32?c4a2b42?c4a3b43?11,c3a2a3b32?c4a4?a2b42?a3b43??,68
- 12 -
c3a2b32?c4a2b42?c4a3b4322211?,c4a2b32b43?,1224
于是总结有如下四阶Runge—Kutta一般显格式中13个参数所满足的11个参数方程
a2?b21,a3?b31?b32,a4?b41?b42?b43,c1?c2?c3?c4?1,c2a2?c3a3?c4a4?1,2
c2a22?c3a32?c4a42?11,c2a23?c3a33?c4a43?,3411,c3a2a3b32?c4a4?a2b42?a3b43??,68
c3a2b32?c4a2b42?c4a3b43?c3a22b32?c4a22b42?c4a32b43?11,c4a2b32b43?,1224 （24）
于是经典的四阶显示Runge—Kutta公式为
?1?yn?1?yn??K1?2K2?2K3?K4?,6??K1?hf?xn,yn?,?h1???K?hfx?,y?K1?,?2n?n22????h1???K3?hf?xn?,yn?K2?,22????K?hf?x?h,y?K?.nn3?4 (25)
四阶Runge—Kutta方法中比较常用的显格式为经典的公式(25)，也是根据上述参数方程(24)来确定一般格式中的参数所得到的。因为显格式(17)共有13个参数，而总计11个参数方程，由解的判定定理知，此方程(24)有无穷多个解，从而有四阶Runge—Kutta方法中的显格式不是唯一的。我们可以尝试推出新算法，但目前为止这个算法是最实用的。
Runge—Kutta方法作为最重要的单步方法，是一类具有相当实用价值的方法。它关于初值是稳定的，其解连续地依赖于初值．它是一类便于应用的单步方法，为了计算yn?1，只用到前一步的值yn即可，因此每步的步长可以独立取定，可以按照绝对稳定性、精度等项要求随时更换。常用的Runge—Kutta方法精度较高，为了达到预定的精度，与Euler方法和梯形法相比，步长办可取得大一些，求解区间上的总步数可以少一些。但Runge—Kutta方法也有一些缺点，比如四阶Runge
- 13 -
—Kutta方法每算一步需四次计算f(x,y)的值，计算量较大(对于较复杂的f(x,y)而言)。
2、 线性多步法
在逐步推进的求解过程中，计算yn?1之前事实上已经求出了一系列的近似值
y0,y1,…,yn，如果充分利用前面多步的信息来预测yn?1，则可以期望会获得较高的
精度．这就是构造所谓线性多步法的基本思想．构造多步法的主要途径是基于数值积分方法和基于泰勒展开方法，前者可直接由方程(1)两端积分后利用插值求积公式得到．一般的线性多步法公式可表示为
yn?1???iyn?i?h??ifn?i， (26)
ii?0k?1k其中yn?i为y(xn?i)的近似，fn?i?f(xn?i,yn?i)，xn?1?x0?ih，?i,?i为常数，?0,?0不全为零，则称(26)为线性k步法，计算时需先给出前面k个近似值y0,y1,…,yk?1,再由（26)逐次求出yk,yk?1,….如果?k?0，称(26)为显式k步法，这时yn?k可直接由(26)算出；如果?k?0，则(26)称为隐式k步法，求解时与改进欧拉法相同，要用迭代法方可算出yn?k，(26)中系数?i,?i可根据方法的局部截断误差及阶确定，其定义为：设y(x)是初值问题(1),(2)的准确解。
阿当姆斯显式与隐式公式
考虑形如
yn?k?yn?k?1?h??ifn?i， (27)
i?0k的k步法，称为阿当姆斯(Admas)方法．?k?0为显式方法，亦称Adams-Bashforth公式；?k?0为隐式方法，亦称Adams-Monlton公式，直接由方程(1)两端从xn?k?1到xn?k积分求得。
例3 用四阶阿当姆斯显式和隐式方法解初值问题
?y???y?x?1, ?y(0)?1.?取步长h?0.1．
解 本题fn??yn?xn?1,xn?nh?0.1n．从四阶阿当姆斯显式公式得到
- 14 -
yn?4?yn?3??h?55fn?3?59fn?2?37fn?1?9fn?241?18.5yn?3?5.9yn?2?3.7yn?1?0.9yn?0.24n?3.24?24
对于四阶阿当姆斯隐式公式得到
h?9fn?3?19fn?2?5fn?1?fn?24yn?3?yn?2?1???0.9yn?3?22.1yn?2?0.5yn?1?0.1yn?0.24n?3?24
由此可直接解出yn?3而不用迭代，得到
yn?3?1?22.1yn?2?0.5yn?1?0.1yn?0.24n?3?． 24.9计算结果如表3，其中显式方法中的y0,y1,y2,y3及隐式方法中的y0,y1,y2均用准确解y(x)?e?x?x计算得到，对一般方程，可用四阶R-K方法计算初始近似．
- 15 -
表3 例3计算结果
阿当姆斯显式公式 阿当姆斯隐式格式 1.070032292 1.10653548 1.14881841 1.19659340 1.24933816 1.30657962 1.36788996 2.87e-006 4.82e-006 1.04081801 1.07031966 2.1 e-007 3.9 e-007 5.2 e-007 6.3 e-007 7.1 e-007 7.7 e-007 8.2 e-007 8.5 e-007 0.3 1.04081822 0.4 1.07032005 0.5 1.10653066 0.6 1.14881164 0.7 1.19658530 0.8 1.24932896 0.9 1.30656966 1.0 1.36787944 6.77 e-006 1.10653014 8.10 e-006 1.14881101 9.20 e-006 1.19658459 9.96 e-006 1.24932819 1.05 e-006 1.30656884 1.36787859 从以上例子看到同阶的阿当姆斯方法，隐式方法要比显式方法误差小，这可以从
p?1两种方法的局部截断误差主项cp?1hp?1y???xn?的系数大小得到解释，这里cp?1分别
为251/720及－19/720．
总结
我们分析了解常微分方程的数值解法，从中可知，不同的数值解法具有各自的优势，会导致不同的误差，从而使得数值方法给出不同的数值结果．以此为基础，我们可以进一步寻求及长时间研究能给出常微分方程良好数值解的数值解法．
- 16 -
参考文献:
[1]李庆扬,王能超,易大义.数值分析[M].第5版.清华大学出版社,2008. [2]A.Iserle.微分方程数值分析基础教程[M].清华大学出版社,2005. [3]王高雄,朱思铭等.常微分方程[M].第3版.高等教育出版社,2006. [4]关治,陆金甫.数值方法[M].清华大学出版社,2006. [5]黄云清等.数值计算方法[M].科学出版社,2009.
[6]胡建伟.常微分方程数值方法[M].第2版.科学出版社,2007.
[7]林立军,郭松云.常微分方程数值解法—Runge—Kutta法的历史浅析[J]．辽宁师范大学学报(自然科学版)26（2）：117—120，2003．6
[8]徐萃微,孙绳武.计算方法引论[M]．北京：高等教育出版社, 2002.
[9]李信真,车刚明,欧阳洁.计算方法[M]．西安：西北工业大学出版社,2000 [10]冯康.数值计算方法[M]．杭州：浙江大学出版社,2003．
- 17 -
