十常微分方程的解法.doc

179-第十五章第十五章 常微分方程的解法常微分方程的解法建立微分方程只是解决问题的第一步，通常需要求出方程的解来说明实际现象， 并加以检验如果能得到解析形式的解固然是便于分析和应用的，但是我们知道，只 有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样 的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如，于是对于用微分方程解决实际问题来说，数值解法就22xydxdy是一个十分重要的手段§1 常微分方程的离散化 下面主要讨论一阶常微分方程的初值问题，其一般形式是(1) 0)(),(yaybxayxfdxdy在下面的讨论中，我们总假定函数连续，且关于满足李普希兹(Lipschitz)条),(yxfy 件，即存在常数，使得L||| ),(),(|yyLyxfyxf 这样，由常微分方程理论知，初值问题(1)的解必定存在唯一 所谓数值解法，就是求问题（1）的解在若干点)(xybxxxxaNL210 处的近似值的方法，称为问题（1）的数值解，),, 2 , 1(NnynL),, 2 , 1(NnynL称为由到的步长。

今后如无特别说明，我们总取步长为常量nnnxxh1nx1nxh建立数值解法，首先要将微分方程离散化，一般采用以下几种方法： （i）用差商近似导数若用向前差商代替代入（1）中的微分方程，则得hxyxynn)()(1)( 'nxy), 1 , 0())(,()()(1Lnxyxfhxyxynnnn化简得 ))(,()()(1nnnnxyxhfxyxy 如果用的近似值代入上式右端，所得结果作为的近似值，记为，)(nxyny)(1nxy1ny则有), 1 , 0(),(1Lnyxhfyynnnn （2） 这样，问题（1）的近似解可通过求解下述问题  )(), 1 , 0(),(01 ayynyxhfyynnnnL-180-（3） 得到，按式（3）由初值可逐次算出式（3）是个离散化的问题，称为0yL,,21yy差分方程初值问题 需要说明的是，用不同的差商近似导数，将得到不同的计算公式 （ii）用数值积分方法 将问题（1）的解表成积分形式，用数值积分方法离散化例如，对微分方程两 端积分，得(4)1), 1 , 0())(,()()(1nnxxnnndxxyxfxyxyL右边的积分用矩形公式或梯形公式计算。

（iii）Taylor 多项式近似 将函数在处展开，取一次 Taylor 多项式近似，则得)(xynx ))(,()()( ')()(1nnnnnnxyxhfxyxhyxyxy 再将的近似值代入上式右端，所得结果作为的近似值，得到离)(nxyny)(1nxy1ny散化的计算公式),(1nnnnyxhfyy 以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式 的计算公式其中的 Taylor 展开法，不仅可以得到求数值解的公式，而且容易估计截 断误差§2 欧拉（Euler）方法 2.1 Euler 方法 Euler 方法就是用差分方程初值问题（3）的解来近似微分方程初值问题（1）的 解，即由公式（3）依次算出的近似值这组公式求问题（1）)(nxy), 2 , 1(Lnyn的数值解称为向前 Euler 公式 如果在微分方程离散化时，用向后差商代替导数，即，则得计算公式hxyxyxynn n)()()( '1 1   )(), 1 , 0(),(0111 ayynyxhfyynnnnL（5） 用这组公式求问题（1）的数值解称为向后 Euler 公式。

向后 Euler 法与 Euler 法形式上相似，但实际计算时却复杂得多向前 Euler 公式 是显式的，可直接求解向后 Euler 公式的右端含有，因此是隐式公式，一般要1ny用迭代法求解，迭代公式通常为  ), 2 , 1 , 0(),(),()( 11)1( 1)0( 1 Lkyxhfyyyxhfyyk nnnk nnnnn（6） 2.2 Euler 方法的误差估计对于向前 Euler 公式（3）我们看到，当时公式右端的都是近似的，L, 2 , 1nny-181-所以用它计算的会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的1ny所谓局部截断误差 假定用（3）式时右端的没有误差，即，那么由此算出ny)(nnxyy ))(,()(1nnnnxyxhfxyy （7） 局部截断误差指的是，按（7）式计算由到这一步的计算值与精确值nx1nx1ny之差为了估计它，由 Taylor 展开得到的精确值是)(1nxy11)(nnyxy)(1nxy)()( ' '2)( ')()(321hOxyhxhyxyxynnnn（8） （7） 、 （8）两式相减（注意到）得),('yxfy )()()( ' '2)(23211hOhOxyhyxynnn（9）即局部截断误差是阶的，而数值算法的精度定义为：2h若一种算法的局部截断误差为，则称该算法具有阶精度阶精度。

)(1phOp 显然越大，方法的精度越高式（9）说明，向前 Euler 方法是一阶方法，因此p 它的精度不高§3 改进的 Euler 方法 3.1 梯形公式 利用数值积分方法将微分方程离散化时，若用梯形公式计算式(4)中之右端积分， 即))](,())(,([2))(,(111  nnnnxxxyxfxyxfhdxxyxfnn 并用代替，则得计算公式1,nnyy)(),(1nnxyxy)],(),([2111nnnnnnyxfyxfhyy这就是求解初值问题（1）的梯形公式 直观上容易看出，用梯形公式计算数值积分要比矩形公式好梯形公式为二阶方 法 梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为 (10)    ), 2 , 1 , 0( )],(),([2),()( 11)1( 1)0( 1Lkyxfyxfhyyyxhfyyk nnnnnk nnnnn由于函数关于满足 Lipschitz 条件，容易看出),(yxfy||2||)1( 1)( 1)( 1)1( 1  k nk nk nk nyyhLyy-182-其中为 Lipschitz 常数。

因此，当时，迭代收敛但这样做计算量较大L120hL如果实际计算时精度要求不太高，用公式（10）求解时，每步可以只迭代一次，由此 导出一种新的方法—改进 Euler 法 3.2 改进 Euler 法 按式（5）计算问题（1）的数值解时，如果每步只迭代一次，相当于将 Euler 公 式与梯形公式结合使用：先用 Euler 公式求的一个初步近似值，称为预测值，1ny1ny然后用梯形公式校正求得近似值，即1ny校正预测] ),(),([2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy（11） 式（11）称为由 Euler 公式和梯形公式得到的预测—校正系统，也叫改进 Euler 法 为便于编制程序上机，式（11）常改写成)(21),(),(1qpnpnnqnnnpyyyyhxhfyyyxhfyy（12） 改进 Euler 法是二阶方法§4 龙格—库塔（Runge—Kutta）方法 回到 Euler 方法的基本思想—用差商代替导数—上来实际上，按照微分中值定 理应有10),( ')()(1hxyhxyxynnn注意到方程就有),('yxfy ))(,()()(1hxyhxhfxyxynnnn （13）不妨记，称为区间上的平均斜率。

可见给出一))(,(hxyhxfKnn],[1nnxx种斜率， （13）式就对应地导出一种算法K 向前 Euler 公式简单地取为，精度自然很低改进的 Euler 公式可理),(nnyxfK解为取，的平均值，其中，这种处K),(nnyxf),(11nnyxf),(1nnnnyxhfyy 理提高了精度 如上分析启示我们，在区间内多取几个点，将它们的斜率加权平均作为],[1nnxx，就有可能构造出精度更高的计算公式这就是龙格—库塔方法的基本思想K 4.1 首先不妨在区间内仍取 2 个点，仿照（13）式用以下形式试一下],[1nnxx-183-1,0),,(),()(12122111hkyhxfkyxfkkkhyynnnnnn（14） 其中为待定系数，看看如何确定它们使（14）式的精度尽量高为此我们,,,21 分析局部截断误差，因为，所以（14）可以化为11)(nnyxy)(nnxyy )())(,( ))(,())(,( ))(,()( '))(,()()(2 112122111hOxyxfhkxyxhfxyxfhkxyhxfkxyxyxfkkkhxyynnynnxnnnnnnnnn（15） 其中在点作了 Taylor 展开。

（15）式又可表为2k))(,(nnxyx)()()( ')()(32 2211hOfffhxhyxyyyxnnn注意到)()( ' '2)( ')()(321hOxyhxhyxyxynnnn中，，可见为使误差，只须令fy 'yxfffy' ')()(3 11hOyxynn，， 1212121（16） 待定系数满足（16）的（15）式称为 2 阶龙格—库塔公式由于（16）式有 4 个未知数而只有 3 个方程，所以解不唯一不难发现，若令 ，，即21211为改进的 Euler 公式可以证明，在内只取 2 点的龙格—库塔公式精度最高],[1nnxx为 2 阶 4.2 4 阶龙格—库塔公式 要进一步提高精度，必须取更多的点，如取 4 点构造如下形式的公式：),(),(),(),()(3625143423122311121443322111hkhkhkyhxfkhkhkyhxfkhkyhxfkyxfkkkkkhyynnnnnnnnnn（17） 其中待定系数共 13 个，经过与推导 2 阶龙格—库塔公式类似、但更复杂的iii,,-184-计算，得到使局部误差的 11 个方程。

取既满足这些方程、又)()(5 11hOyxynn 较简单的一组，可得iii,,),()2,2()2,2(),()22(6342 31 2143211hkyhxfkhkyhxfkhkyhxfkyxfkkkkkhynnnnnnnnn（18） 这就是常用的 4 阶龙格—库塔方法（简称 RK 方法） §5 线性多步法 以上所介绍的各种数值解法都是单步法，这是因为它们在计算时，都只用到1ny前一步的值，单步法的一般形式是ny) 1,, 1 , 0(),,(1NnhyxhyynnnnL（19） 其中称为增量函数，例如 Euler 方法的增量函数为，改进 Euler 法),,(hyx),(yxf 的增量函数为))],(,(),([21),,(yxhfyhxfyxfhyx如何通过较多地利用前面的已知信息，如，来构造高精度的算法rnnnyyy,,,1L计算，这就是多。