科学计算—理论、方法及其基于MATLAB 的程序实现与分析微分方程（组）数值解法§1 常微分方程初值问题的数值解法微分方程（组）是科学研究和工程应用中最常用的数学模型之一。
如揭示质点运动规律的Newton 第二定律：()()()⎪⎩⎪⎨⎧'='==000022x t x x t x t F dt xd m （１） 和刻画回路电流或电压变化规律的基尔霍夫回路定律等，但是，只有一些简单的和特殊的常微分方程及常微分方程组，可以求得用公式给出的所谓“解析解”或“公式解”，如一阶线性微分方程的初值问题：()()00y y t f ay dtdy=+= （２） 的解为：()()()τττd f e y e t y tt a at⎰-+=00 （３）但是，绝大多数在实际中遇到的常微分方程和常微分方程组得不到“解析解”，因此，基于如下的事实：１、绝大多数的常微分方程和常微分方程组得不到（有限形式的）解析解；２、实际应用中往往只需要知道常微分方程（组）的解在（人们所关心的）某些点处的函数值（可以是满足一定精度要求的近似值）；如果只需要常微分方程（组）的解在某些点处的函数值，则没有必要非得通过求得公式解，然后再计算出函数值不可，事实上，我们可以采用下面将介绍的常微分方程（组）的初值问题的数值解法，就可以达到这一目的。
一般的一阶常微分方程（组）的初值问题是指如下的一阶常微分方程（组）的定解问题：()()000,y t y t t t y t F dtdyf=≤≤= （７）其中()()()()⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=t y t y t y t y n 21 （８）()()()()⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=y t f y t f y t f y t F n ,,,,21（９） 常微分方程（组）的初值问题通常是对一动态过程（动态系统、动力系统）演化规律的描述，求解常微分方程（组）的初值问题就是要了解和掌握动态过程演化规律。
§１.1 常微分方程（组）的Cauch 问题数值解法概论假设要求在点（时刻）k t ，n k ,,2,1 =处初值问题（７）的解()00,,y t t y y =的（近似）值,如果已求得k t 时刻的值()k t y 或它的近似值k y (如0t 时刻的值()0t y )，那么将式（７）的两端在区间[]1,+k k t t 上积分()()()()dt t y t f t y t y dt dtdyk kk kt t k k t t ⎰⎰++=-=+11,1 (10)可得()()()()dt t y t f t y t y k kt t k k ⎰++=+1,1 (11)或()()()dt t y t f y y t y k kt t k k k ⎰++=≈++1,11 (12)显然,为了利用式（11）或（12）求得()1+k t y 的精确值（近似值1+k y ），必须计算右端的积分，这是问题的关键也是难点所在，如前所述，一般得不到精确的公式解，因此需要采用数值积分的方法求其近似解, 可以说,不同的式值积分方法将给出不同的Cauch 问题的数值解法。
§1.2 最简单的数值解法——Euler 方法假设要求在点（时刻）kh t t k +=0，nt t h f 0-=，n k ,,2,1 =处初值问题（７）的解()t y y =的近似值。
首先对式（７）的两端积分，得()()()()dt t y t f t y t y dt dtdyk k k kt t k k t t ⎰⎰++=-=+11,1 （13） 对于式（13）的右边，如果用积分下限k t t =处的函数值()()k k t y t f ,代替被积函数作积分（从几何上的角度看，是用矩形面积代替曲边梯形面 积），则有()()()()()[]k k t tk k t y t hf dt t y t f t y t y k k,,11≈=-⎰++ （14）进而得到下式给出的递推算法—Euler 方法()()001,,2,1,y t y nk y t hf y y k k k k ==+=+ （15）例１ 用Euler 方法解如下初值问题，取3.0=h ，()1030sin 23=≤≤-=y t t y y dtdy解：由（15）得()1010,,2,1sin 6.03.131==-=+y k t y y y kk k k结果如下：如果取1.0=h ，其结果如下图所示：Euler_Method§1.3 改进的Euler 方法对于（15）的右边，如果被积函数用积分限k t t =和1+=k t t 处的函数值的算术平均值代替（几何上，是用梯形面积代替曲边梯形面积），则有()()()()()()()()[]111,,2,1++++≈=-⎰+k k k k t t k k t y t f t y t f hdt t y t f t y t y k k（16） 进而得到下式给出的递推算法：()()[]()00111,,2,1,,2y t y nk y t f y t f hy y k k k k k k ==++=+++ （17）通常算法（17）比Euler 方法（15）的精度高，但是，按算法（17）求1+k y 时要解（非线性）方程（组），这是算法（17）不如Euler 方法的方面，为了１） 尽可能地保持算法（17）精度高的优点； ２） 尽可能地利用Euler 方法计算简单的长处； 人们采取了如下的称之为改进的Euler 方法的折衷方案：预测 ()k k k k y t hf y y ,01+=+ （18） 修正 ()()[]n k y t f y t f h y y k k k k k k ,,2,1,,2111 =++=+++ （19）()00y t y =例２ Euler 方法与改进的Euler 方法的比较 下图是当3.0=h 时比较的结果：open Improved_Euler_Method.m§1.4 Euler 方法和改进的Euler 方法的误差分析由Taylor 公式()()()()()()()()()()22111,!21h O h t y t f t y t t t y t t t y t y t y k k k k k k k k k k k ++=+-''+-'+=+++（19） 说明Euler 方法的截断误差是()2h O ，类似地，由()()()()()()321!21,h O h t y h t y t f t y t y k k k k k +''++=+ （20） ()()()()()()321111!21,h O h t y h t y t f t y t y k k k k k +''+-=++++ （21） 以及()()()()()()h O t y t y h t y t y t y k k k k k +''=''⇒+'''+''=''++11 (22)让式（20）的两端减式（21）的两端，可得()()()()()()[]()()()()()()()[]()31132111,,241,,2h O t y t f t y t f ht y h O h O h t y t f t y t f ht y t y k k k k k k k k k k k +++=++++=+++++ （23）从上述推导Euler 方法、改进的Euler 方法的过程以及例１、例２容易看出，改进的Euler 方法Euler 方法的精度高，其原因在于： 1 在推导Euler 方法时，我们是用待求解函数()t y y =在一点处的变化 率()()k k t y t f ,代替()t y y =在区间],[1+k k t t 上的平均变化率：()()()()t y t hf dt t y t f k kt t ,,1=⎰+ （24）2 而在推导改进的Euler 方法时，我们是用待求解函数()t y y =在两点处变化率的平均值()()()()[]11,,21+++k k k k t y t f t y t f 代替()t y y =在区间],[1+k k t t 上的平均变化率；显然，通常()()()()[]11,,21+++k k k k t y t f t y t f 比()()k k t y t f ,更接近于()t y y =在区间],[1+k k t t 上的平均变化率。
由此启发人们：适当地选取区间],[1+k k t t 上函数()t y y =若干点处的变化率，用它们加权平均值代替()t y y =在区间],[1+k k t t 上的平均变化率，近似解的精度应更高。
下面将要介绍的Runge —Kutta 法就是基于上述想法得到的。
§2 Runge —Kutta 法Runge —Kutta 法是按选取区间],[1+k k t t 上函数()t y y =变化率()y t f ,的个数的多少和截断误差的阶数()m h O 来区分的一系列方法，如 １ 二阶的Runge —Kutta 法（改进的Euler 方法）()()()⎪⎪⎩⎪⎪⎨⎧++=+==++21111212,,K K hy y hK y t f K y t f K k k k k k k （25） ２ 三阶的Runge —Kutta 法()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧+++=++=⎪⎭⎫ ⎝⎛++==+32112312146,2,2,K K K h y y hK y h t f K K h y h t f K y t f K k k k k k k k k （26） ３ 四阶的Runge —Kutta 法 1) 古典形式()()()⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧++++=++=⎪⎭⎫ ⎝⎛++=⎪⎭⎫ ⎝⎛++==+432113423121226,2,22,2,K K K K hy y hK y h t f K K h y h t f K K h y h t f K y t f K kk k k k k k k k k （27） 2) Gill 公式(具有减小舍入误差的优点)()()()()⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎨⎧+++-++=⎪⎪⎭⎫ ⎝⎛++-+=⎪⎪⎭⎫ ⎝⎛-+-++=⎪⎭⎫ ⎝⎛++==+432113242131212222622222,222212,22,2,K K K K hy y hK hK y h t f K hK hK y h t f K K h y h t f K y t f K k k k k k k k k k k （28） 4 Runge —Kutta 法的一般形式()∑=++=+=ni i i k k k k k K a h y h y t h y y 11,,φ (29)其中()h y t k k ,,φ称为增量函数(Increment Function)以及()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧++++=+++=++==-----11,111,1122212123111121,,,,n n n n k n k n k k k k k k hK q hK q y h p t f K hK q hK q y h p t f K hK q y h p t f K y t f K （30）需要特别指出的是：在确定（29）、（30）中的参数i a ，i p 和ij q 时，应该使（29）的右端和适当阶的（如n 阶）Taylor 展式一致，这样，至少对于底阶的R-K 法来说，参与加权的斜率个数n 与方法的阶数是一致的。