151第八章 常微分方程数值解在工程和科学技术的实际问题中，常常需求解常微分方程。

但由常微分方程理论可知，常微分方程中往往只有少数较简单和典型的方程可求出其解析解。

在大多数情况下，常微分方程只能用近似法求解。

这种近似解法可分为两大类：一类是近似解析法，如级数解法、逐次逼近法等；另一类则是数值解法，它给出方程在一些离散点上的近似解。

本章主要讨论一阶常微分方程的初值问题：()()⎪⎩⎪⎨⎧==0,y a y y x f dx dyb x a ≤≤ （8.1） 从理论上讲，只要方程中的()y x f ,连续且关于y 满足李普希兹（Lipschitz ）条件，即存在常数L ，使()()2121,,y y L y x f y x f -≤-则常微分方程存在唯一解)(x y y =。

所谓微分方程数值解，就是求微分方程的解()x y 在一系列离散节点 b x x x x a n n =<<<<=-110处()i x y 的近似值i y ),,1,0(n i =. 相邻的两个节点之间的距离i i i x x h -=+1称为由i x 到1+i x 的步长，通常取为常数h 。

求数值解，首先应将微分方程离散化，常用的方法有： （1） 用差商代替微商 若用向前差商代替微商，即()()()()()i i i i i x y x f x y hx y x y ,1='≈-+ )1,,1,0(-=n i代入（8.1）中的微分方程，则得()1+i x y ()()()i i i x y x hf x y ,+≈152 记)(i x y 的近似值i y ，则由上式右端可计算出)(1+i x y 的近似值，即()i i i i y x hf y y ,1+=+ )1,,1,0(-=n i （8.2）（2） 数值积分法 利用数值积分法左矩形公式()()i i x y x y -+1=()()()i i x x y x hf dx x y x f i i,,1≈⎰+可得同样算法 ()i i i i y x hf y y ,1+=+（3） 用泰勒（Taylor ）公式将函数)(x y 在i x 处展开，取一次Taylor 多项式近似，则得()()h x y x y i i +=+1()()i i x y h x y '+≈()()()i i i x y x hf x y ,+=从而也得到离散化得计算公式 ()i i i i y x hf y y ,1+=+§1 欧拉（Euler ）方法1.1欧拉方法对一阶微分方程（8.1），把区间[]b a ,作n 等分：b x x x x a n n =<<<<=-110 ， 则分点为 ih a x i +=， nab h -=),2,1(n i = 由以上讨论可知，无论用一阶向前差商，还是用数值积分法左矩形公式，或者用泰勒公式取前两项都可得到同样的离散化计算公式()i i i i y x hf y y ,1+=+并将初值条件代入，则得到数值算法：()()⎩⎨⎧=+=+a y y y x hf y y i i i i 01, ),2,1(n i = （8.3） 称其为欧拉方法。

几何上欧拉方法就是用一条折线近似表示曲线()x y y =（如图8-1）。

因此欧拉方法又称为欧拉折线方法。

153y图8-1 欧拉方法 1.2欧拉方法的误差估计定义1 假设)(i i x y y =为准确值，考虑计算一步所产生得误差，即用某种数值算法计算)(1+i x y 1+i y 所产生的误差()111+++-=i i i y x y R ，称为该数值算法的局部截断误差。

定义2 考虑用某种数值算法计算时，因前面的计算不准确而引起的准确解)(i x y 与数值解i y 的误差，()i i i y x y e -=称为该数值算法的整体截断误差。

设函数),(y x f 充分光滑，问题（8.1）的解()x y 在],[b a 上有二阶连续导数，由泰勒公式有()()h x y x y i i +=+1=()()()i i i y h x y h x y ξ''+'+221=()()i i i i y h y x hf y ξ''++221,所以154 ()111+++-=i i i y x y R =)(212i y h ξ''，1+<<i i i x x ξ （8.4） 定义3 如果一数值解法的局部截断误差为)(1+p h O ，则称该算法为p 阶算法。

当h 充分小时，由（8.4）知欧拉方法的局部截断误差为)(2h O ，因此欧拉方法是一个 一阶方法，计算结果的精度较差。

1.3 改进的欧拉方法由微分方程数值解的三种基本构造方法知，若取不同的差商（如向后差商），不同的数 值积分公式（如梯形公式），以及泰勒公式取前三项、四项等可得不同的算法。

如果用向后差商近似代替导数，则有()()()()()i i i i i x y x f x y hx y x y ,1='≈-- ),,1(n i =即 ()111,)()(++++≈i i i i y x hf x y x y )1,,1,0(-=n i所以有 ()111,++++=i i i i y x hf y y )1,,1,0(-=n i （8.5） （8.5）式称为隐式欧拉公式。

如果用梯形公式计算积分：()()()()()()[]11,,2,1+++≈⎰+i i i i x x x y x f x y x f hdx x y x f i i()()[]111,,2+++++=i i i i i i y x f y x f hy y （8.6） 且 ()111+++-=i i i y x y R = ()ξy h '''-3121（8.7）由于此方程为1+i y 的隐式方程，不易求解。

一般将其与欧拉方法联合使用，可得算法()()()()()()[]⎪⎩⎪⎨⎧++=+=+++++k i i i i i k i i i i i y x f y x f h y y y x hf y y 111101,,2, （8.8）)1,2,1,0;,2,1,0(-==n i k按式（8.8）计算问题（8.1）的数值解时，如果每步只迭代一次，相当于将欧拉公式与梯形公式结合使用，即在实际计算中，当h 比较小时，常取一次迭代后的近似值()11+i y 为1+i y ，155于是有改进的欧拉方法⎪⎩⎪⎨⎧++=+=++++)]~,(),([2),(~1111i i i i i i i i i i y x f y x f h y y y x hf y y )1,2,1,0(-=n i （8.9） 例1 用欧拉方法和改进的欧拉方法求微分方程()[]7.0,0,10322∈⎪⎩⎪⎨⎧==x y xydx dy 的数值解（取h=0.1）。

解 由欧拉方法（8.3），得数值计算公式1+i y =i y ＋0.1×232i i y x 计算结果如表8-1由改进的欧拉方法（8.8），得数值计算公式⎪⎩⎪⎨⎧++=+=++++]~3232[05.0321.0~2112121i i i i i i i i i i y x y x y y y x y y 计算结果如表8-2表8-1x i 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 y i 1.0000 1.0069 1.0208 1.0391 1.0628 1.0923 1.1269 1.1643 误差 0.0000 0.0037 0.0077 0.0993 0.0120 0.0151 0.0189 0.0222表8-2x i 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 y i 1.0000 1.0033 1.0132 1.0292 1.0506 1.0773 1.1079 1.1422 误差 0.0000 0.0000 0.0000 0.0002 0.0020 0.0101 0.0122 0.0053例2 用欧拉法、改进欧拉法求微分方程数值解（h=0.1）。

156 ⎪⎩⎪⎨⎧=-='1)0(y y x y y解 由欧拉方法（8.3），得数值计算公式1+i y =i y ＋0.1⎪⎪⎭⎫⎝⎛-i i i y x y 由改进的欧拉方法（8.8），得数值计算公式⎪⎪⎩⎪⎪⎨⎧-+-+=-+=+++++)]~~()[(05.0)(1.0~11111i i i i i i i i i i i i i y xy y x y y y y x y y y计算结果如下§2 龙格－库塔（Runge －Kutta ）方法2.1 泰勒展开法由于欧拉方法为一阶方法，为了提高算法的阶，有必要讨论更高阶的方法。

在泰勒展 开式中取更多的项，如取p ＋1项可得p 阶算法。

()p i p i i i i y p h y h y h y y !!221++''+'+=+157()()()i p p i yp h R ξ111!1++++= 其中()k i y ）可用复合函数求导法则计算。

如p=2时得二阶泰勒方法()221(,),[]2!2i i i i i i i i i x y x y h h y y hy y y hf x y f f f +'''''=++=+++⋅2.2 龙格－库塔法为了避免计算高阶导数，龙格－库塔方法利用()y x f ,某些点处的值的线性组合构造计算公式，使其按泰勒公式展开后与初值问题解的泰勒展开式比较，有尽可能多的项相同。

龙格－库塔法的一般形式为：()()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧++=++==++++=+h y h x f K h y h x f K y x f K K a K a K a h y y m i m i m i i i m m i i μλμλ,,,222122111 （8.10） 下面以二阶龙格－库塔法为例说明龙格－库塔法的构造过程。

二阶龙格－库塔公式为（8.11）其中221,,λa a 为待定参数。

事实上，将K 2在()i i y x ,处按泰勒公式展开，则()i i i i i y x hf a y hK a hK a y y ,122111+=++=+＋()()()()()],,,,[2222h O y x f yy x hf y x f x hy x f h a i i i i i i i i +∂∂+∂∂+λλ ()()()⎪⎩⎪⎨⎧++==++=+1222122111,,hK y h x f K y x f K K a K a h y y i i i i i i λλ158 ()()()()()()322221,,,,h O y x f y y x f y x f x h a y x hf a a y i i i i i i i i i +⎥⎦⎤⎢⎣⎡∂∂+∂∂+++λ=()()()()()322221h O x y h a x y h a a x y i i i +''+'++λ 另一方面()()h x y x y i i +=+1=()()()()3221h O x y h x y h x y i i i +''+'+于是为使局部截断误差的阶尽可能高，应使⎪⎩⎪⎨⎧==+2112221λa a a 方程组有无穷多组解，取定参数则得到许多具体的二阶龙格－库塔公式。