常微分方程初值问题数值解法

1第9章常微分方程初值问题数值解法9.1引言本章要着重考察的一阶方程的初值问题.)(),,(00yxyyxfy（1.1）（1.2）只要函数适当光滑——譬如关于满足利普希茨(Lipschitz）条件),(yxfy.),(),(yyLyxfyxf（1.3）理论上就可以保证初值问题（1.1），（1.2）的解存在并且唯一.)(xyy2所谓数值解法，就是寻求解在一系列离散节点)(xy121nnxxxx上的近似值.相邻两个节点的间距称为步长.,,,,,121nnyyyynnnxxh1如不特别说明，总是假定为定数，这时节点为.),2,1(ihhi),2,1,0(0inhxxn初值问题（1.1），（1.2）的数值解法的基本特点是采取“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进.描述这类算法，只要给出用已知信息计算的递推公式.,,,21nnnyyy1ny首先，要对方程（1.1）离散化，建立求数值解的递3推公式.一类是计算时只用到前一点的值，称为单步法.另一类是用到前面点的值，称为步法.1nyny1nyk11,,,knnnyyyk其次，要研究公式的局部截断误差和阶，数值解与精确解的误差估计及收敛性，还有递推公式的计算稳定性等问题.ny)(nxy49.2简单的数值方法与基本概念9.2.1欧拉法与后退欧拉法在平面上，微分方程（1.1）的解称作它的积分曲线.xy)(xyy积分曲线上一点的切线斜率等于函数的值.),(yx),(yxf如果按函数在平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.),(yxfxy基于上述几何解释，从初始点出发，先依方向场在该点的方向推进到上一点，然后再从依方向场的方向推进到上一点，循此前进做出),(000yxP1xx1P1P2xx2P5一条折线（图9-1).210PPP一般地，设已做出该折线的顶点，过依方向场的方向再推进到，显然两个顶点的坐标有关系),(nnnyxPnP),(111nnnyxP1,nnPP图9-16),,(11nnnnnnyxfxxyy即).,(1nnnnyxhfyy（2.1）这就是著名的欧拉（Euler）公式.若初值已知，则依公式（2.1）可逐步算出0y),,(),,(11120001yxhfyyyxhfyy例1求解初值问题.1)0(),10(2yxyxyy（2.2）7解欧拉公式的具体形式为).2(1nnnnnyxyhyy取步长，计算结果见表9-1.1.0h7321.17848.10.14142.14351.15.06733.17178.19.03416.13582.14.06125.16498.18.02649.12774.13.05492.15803.17.01832.11918.12.04832.15090.16.00954.11000.11.0)()(19nnnnnnxyyxxyyx计算结果对比表8初值问题（2.2）有解，按这个解析式子算出的准确值同近似值一起列在表9-1中，两者相比较可以看出欧拉方法的精度很差.xy21)(nxyny还可以通过几何直观来考察欧拉方法的精度.假设，即顶点落在积分曲线上，那么，按欧拉方法做出的折线便是过点的切线（图9-2）.)(nnxyynP)(xyy1nnPP)(xyynP图9-29从图形上看，这样定出的顶点显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的.1nP为了分析计算公式的精度，通常可用泰勒展开将在处展开，则有)(1nxynx).,()(2)()()()(121nnnnnnnnxxyhhxyxyhxyxy在的前提下，.于是可得欧拉法（2.1）的公式误差)(nnxyy)())(,(),(nnnnnxyxyxfyxf),(2)(2)(2211nnnnxyhyhyxy（2.3）称为此方法的局部截断误差.10如果对方程（1.1）从到积分，得nx1nx.))(,()()(11nnxxnndttytfxyxy（2.4）右端积分用左矩形公式近似，再以代替代替也得到（2.1），局部截断误差也是（2.3）.))(,(nnxyxfhny1),(nnyxy)(1nxy如果在（2.4）中右端积分用右矩形公式近似，则得另一个公式))(,(11nnxyxfh),,(111nnnnyxhfyy（2.5）称为后退的欧拉法.后退的欧拉公式与欧拉公式有着本质的区别，后者是关于的一个直接的计算公式，这类公式称作是显式的；1ny11然而公式（2.5）的右端含有未知的，它实际上是关于的一个函数方程，这类公式称作是隐式的.1ny1ny隐式方程（2.5）通常用迭代法求解，而迭代过程的实质是逐步显示化.设用欧拉公式),()0(1nnnnyxhfyy给出迭代初值，用它代入（2.5）式的右端，使之转化为显式，直接计算得)0(1ny),,()0(11)1(1nnnnyxhfyy然后再用代入（2.5）式，又有)1(1ny).,()1(11)2(1nnnnyxhfyy12如此反复进行，得).,1,0(),,()(11)1(1kyxhfyyknnnkn（2.6）由于对满足利普希茨条件（1.3）.由（2.6）减（2.5）得),(yxfy.),(),(1)(111)(111)1(1nknnnknnnknyyhLyxfyxfhyy由此可知，只要迭代法（2.6）就收敛到解.关于后退欧拉方法的公式误差，从积分公式看到它与欧拉法是相似的.1hL1ny139.2.2梯形方法为得到比欧拉法精度高的计算公式，在等式（2.4）右端积分中若用梯形求积公式近似，并用代替代替，则得ny1),(nnyxy)(1nxy)],,(),([2111nnnnnnyxfyxfhyy（2.7）称为梯形方法.梯形方法是隐式单步法，可用迭代法求解.同后退的欧拉方法一样，仍用欧拉方法提供迭代初值，则梯形法的迭代公式为14).,2,1,0()],(),([2);,()(11)1(1)0(1kyxfyxfhyyyxfhyyknnnnnknnnnn（2.8）为了分析迭代过程的收敛性，将（2.7）式与（2.8）式相减，得)],,(),([2)(1111)1(11knnnnknnyxfyxfhyy于是有,2)(11)1(11knnknnyyhLyy式中为关于的利普希茨常数.),(yxfyL15如果选取充分小，使得h,12hL则当时有，这说明迭代过程（2.8）是收敛的.k1)(1nknyy169.2.3单步法的局部截断误差与阶初值问题(1.1)，(1.2)的单步法可用一般形式表示为),,,,(11hyyxhyynnnnn（2.9）其中多元函数与有关，当含有时，方法是隐式的，若不含则为显式方法，所以显式单步法可表示为),(yxf1ny1ny),,,(1hyxhyynnnn（2.10）称为增量函数，例如对欧拉法（2.1）有),,(hyx).,(),,(yxfhyx它的局部截断误差已由（2.3）给出.17对一般显式单步法则可如下定义.定义1设是初值问题(1.1)，(1.2)的准确解，称)(xyy)),(,()()(11hxyxhxyxyTnnnnn（2.11）为显式单步法（2.10）的局部截断误差.之所以称为局部的，是假设在前各步没有误差.1nTnx当时，计算一步，则有)(nnxyy.)),(,()()()],,([)()(11111nnnnnnnnnnnThxyxhxyxyhyxhyxyyxy18根据定义，显然欧拉法的局部截断误差),()(2)()()())(,()()(3211hOxyhxyhxyhxyxyxhfxyxyTnnnnnnnnn即为（2.3）的结果.这里称为局部截断误差主项.显然.)(22nxyh)(21hOTn的公式误差.所以，局部截断误差可理解为用方法（2.10）计算一步的误差，也即公式（2.10）中用准确解代替数值解产生)(xy19定义2设是初值问题(1.1)，(1.2)的准确解，若存在最大整数使显式单步法(2.10)的局部截断误差满足)(xyp),(),,()()(11pnhOhyxhxyhxyT（2.12）则称方法（2.10）具有阶精度.p若将（2.12）展开式写成),())(,(211ppnnnhOhxyxT则称为局部截断误差主项.1))(,(pnnhxyx以上定义对隐式单步法（2.9）也是适用的.20对后退欧拉法（2.5）其局部截断误差为).()(2)]()()([)()(2)())(,()()(322321111hOxyhhOxyhxyhhOxyhxyhxyxhfxyxyTnnnnnnnnnn这里，是1阶方法，局部截断误差主项为.1p)(22nxyh21对梯形法（2.7）有).()(12)()](2)()()([2)(!3)(2)()]()([2)()(434232111hOxyhhOxyhxyhxyxyhxyhxyhxyhxyxyhxyxyTnnnnnnnnnnnnn所以梯形方法（2.7）是二阶的，其局部误差主项为.)(123nxyh229.2.4改进的欧拉公式梯形方法虽然提高了精度，但其算法复杂，在应用迭代公式（2.9）进行实际计算时，每迭代一次，都要重新计算函数的值，而迭代又要反复进行若干次，计算量很大，而且往往难以预测.),(yxf为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法.具体地，先用欧拉公式求得一个初步的近似值,称之为预测值，预测值的精度可能很差，再用梯形公式（2.7）将它校正一次，即按（2.8）式迭代一次得，这个结果称校正值，而这样建立的预测-校正系统通常称为改进的欧拉公式：1ny1ny23预测),,(1nnnnyxhfyy校正)].,(),([2111nnnnnnyxfyxfhyy（2.13）或表为下列平均化形式).(21),,(),,(11cpnpnncnnnpyyyyxhfyyyxhfyy24例2用改进的欧拉方法求解初值问题（2.2）.解改进的欧拉公式为).(21),2(),2(11cpnpnpncnnnnpyyyyxyhyyyxyhyy仍取，计算结果见表9-2.同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度.1.0h257321.17379.10.14142.14164.15.06733.16782.19.03416.13434.14.06125.16153.18.02649.12662.13.05492.15525.17.01832.11841.12.04832.14860.16.00954.10959.11.0)()(29nnnnnnxyyxxyyx计算结果对比表269.3龙格-库塔方法9.3.1显式龙格-库塔法的一般形式上节给出了显式单步法的表达式),,,(1hyxhyynnnn其局部截断误差为),(),,()()(11pnhOhyxhxyhxyT对欧拉法，即方法为阶，若用改进欧拉法，它可表为)(21hOTn1p))].,(,(),([21nnnnnnnnyxhfyhxfyxfhyy（3.1）27此时增量函数))].,(,(),([21),,(nnnnnnnnyxhfyh