数值计算方法第七章

第七章常微分方程初值问题数值解法7.1引言科学技术中常常需要求解常微分方程的定解问题．这类问题最简单的形式，是本章将要着重考察的一阶方程的初值问题'00,yfxyyxy1.11.2我们只要，只要函数,fxy适当光滑—譬如关于y满足利普希茨Lipschitz条件,,fxyfxyLyy.1.3理论上就可以保证初值问题1.1，1.2的解yyx存在并且唯一.虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法.所谓数值解法，就是寻求解yx在一系列离散节点121nnxxxx上的近似值121,,,,nnyyyy.相邻两个节点的间距1nnnhxx成为步长.今后如不特别说明，总是假定1,2,ihhi为定数，这时节点为0,0,1,2,.nxxnhn初值问题1.1，1.2的数值解法有个基本特点，它们都采取“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进．描述这类算法，只要给出用已知信息,12,,nnnyyy计算1ny的递推公式.首先，要对方程1.1离散化，建立求数值解的递推公式．一类是计算1ny时只用到前一点的值ny，称为单步法．另一类是用到1ny前面k点的值11,,,nnnkyyy称为k步法．其次，要研究公式的局部截断误差和阶，数值解ny与精确解nyx的误差估计及收敛性，还有递推公式的计算稳定性等问题．7.2简单的数值方法与基本概念7.2.1欧拉法与后退欧拉法我们知道，在xy平面上，微分方程1.1的解yyx称为它的积分曲线．积分曲线上一点,xy的切线斜率等于函数,fxy的值．如果按函数,fxy在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点方向相一致．基于上述几何解释，我们从初始点000,pxy出发，先依方向场在该点的方向推进到1xx上一点1p，然后再从1p依方向场的方向推进到2xx上一点2p，循此前进做出一条折线012ppp（图91）．一般地，设已做出该折线的顶点np，过,nnnpxy依方向场的方向再推进到111,nnnpxy，显然两个定点np，1np的坐标有关系11,,_nnnnnnyyfxyxx即1,nnnnyyhfxy2.1这就是著名的欧拉Euler公式．若初值0y已知，则依公式2.1可逐步算出10002111,,,,yyhfxyyyhfxy例1求解初值问题'201xyyyy01x，2.2解为便于进行比较，本章将用多种数值方法求解上述初值问题．这里先用欧拉方法，欧拉公式的具体形式为12.nnnnnxyyhyy取步长0.1,h计算结果见表71表71计算结果对比nxnynyxnxnynyx0.10.20.30.40.51.10001.19181.27741.35821.43511.09541.18321.26491.34161.41420.60.70.80.91.01.50901.58031.64981.71781.78481.48321.54921.61251.67331.7321初值问题2.2有解12yx，按这个解析式子算出的准确值nyx同近似值ny一起列在表71中，两者相比较可以看出欧拉方法的精度很差．还可以通过几何直观来考察欧拉方法的精度．假设nnyyx，即顶点np落在积分曲线yyx上，那么，按欧拉方法做出的切线np1np便是yyx过点np的切线．从图形上看，这样定出顶点1np显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的．为了分析计算公式的精度，通常可采用泰勒展开将1nyx在nx处展开，则有1nnyxyxh2''',2nnnhyxyxhy1,.nnnxx在nnyyx的前提下，',,.nnnnnfxyfxyxyx于是可得欧拉法2.1的公式误差22''''1122nnnnhhyxyyyx，2.3称为此方法的局部截断误差．如果对方程1.1从nx到1nx积分，得11,nnnnxyxyxftytdtx．2.4右端积分用左矩形公式,nnhfxyx近似，再以ny代替nyx，1ny代替1nyx也得到2.1，局部截断误差也是2.3．如果在2.4中右端积分用右矩形公式11,nnhfxyx近似，则得另一个式111,nnnnyyhfxyx2.5称为后退得欧拉法．后退的欧拉公式与欧拉公式有着本质的区别，后者是关于1ny的一个直接的计算公式，这类公式称作是显式的；然而公式2.5的右端含有未知的1ny，它实际上是关于1ny的一个函数方程，这类公式称作是隐式的．显式与隐式两类方法各有特点．考虑到数值稳定性等其他因素，人们有时需要选用隐式方法，但使用显式算法远比隐式方便．隐式方程2.5通常用迭代法求解，而迭代过程的实质是逐步显示化．设用欧拉公式01,nnnnyyhfxy给出迭代初值01ny，用它代入2.5式的右端，使之转化为显示，直接计算得10111,nnnnyyhfxy，然后再用11ny代入2.5式，又有21111,nnnnyyhfxy．如此反复进行，得1111,kknnnnyyhfxy0,1,k．2.6由于,fxy对y满足利普希茨Lipschitz条件1.3．由2.6减2.5得1111111,,kknnnnnnyyhfxyfxy11.knnhLyy由此可知，只要1hL迭代法2.6就收敛到解1ny．关于后退欧拉法的公式误差，从积分公式看到它与欧拉法是相似的．7.2.2梯形方法为得到比欧拉法精度高的计算公式，在等式2.4右端积分中若用梯形求积公式近似，若用ny代替nyx，1ny代替1nyx，则得111,,2nnnnnnhyyfxyfxy2.7称为梯形方法．梯形方法是隐式单步法，可用迭代法求解．同后退的欧拉法一样，仍用欧拉方法提供迭代初值，则梯形法的迭代公式为011111,;,,;20,1,2,.nnnnkknnnnnnyyhfxyhyyfxyfxyk2.8为了分析迭代过程的收敛性，将2.7式与2.8式相减，得1111111,,2kknnnnnnhyyfxyfxy，于是有11111,2kknnnnhLyyyy式中L为,fxy关于y的利普希茨常数．如果选取h充分小，使得12hL，则当k时有1kny1ny，这说明迭代过程2.8是收敛的．7.2.3单步法的局部截断误差与阶初值问题1.1,1.2的单步法可用一般形式表示为11,,,,nnnnnyyhxyyh2.9其中多元函数与,fxy有关，当含有1ny时，方法是隐式的，若不含1ny则为显式方法，所以显式单步法可表示为1,,,nnnnyyhxyh2.10,,xyh称为增量函数，例如对欧拉法2.1有,,xyh,fxy．它的局部截断误差已由2.3给出，对一般显式单步法则可如下定义．定义1设yx是初值问题1.1,1.2的准确解，称11,,nnnnnTyxyxhxyxh2.11为显式单步法2.10的局部截断误差．1nT之所以称为局部的，是假设在nx前各步没有误差．当nynyx时，计算一步，则有111,,nnnnnnyxyyxyhxyh11,,nnnnnyxyxhxyxhT．所以，局部截断误差可理解为用方法2.10计算一步的误差，也即公式2.10中用准确解yx代替数值解产生的公式误差．根据定义，显然欧拉法的局部截断误差11,nnnnnTyxyxhxyx'1nnnyxyxhyx2''3,2nhyxh即为2.3的结果．这里2''2nhyx称为局部截断误差主项．显然21nTh．一般情形的定义如下，定义2设yx是初值问题1.1,1.2的准确解，若存在最大整数p使显式单步法2.10的局部截断误差满足11,,pnTyxhyxhxyhh，2.12则称方法2.10具有p阶精度．若将2.12展开式写成121,ppnnnTxyxhh，则1,pnnxyxh称为局部截断误差主项．以上定义对隐式单步法2.9也是适用的．例如，对后退欧拉法2.5其局部截断误差为1111,nnnnnTyxyxhfxyx2'''32nnhhyxyxh'''2nnhyxhyxh2''3.2nhyxh这里1p，是1阶方法，局部截断误差主项为2''2nhyx．同样对梯形法2.7有''1112nnnnnhTyxyxyxyx23''''''23!nnnhhhyxyxyx2'''''''422nnnnhhyxyxhyxyxh3'''4.12nhyxh所以梯形方法2.7是二阶的，其局部误差主项为3'''12nhyx．7.2.4改进的欧拉公式我们看到，梯形方法虽然提高了精度，但其算法复杂，在应用迭代公式2.9进行实际计算时，每迭代一次，都要重新计算函数,fxy的值，而迭代又要反复进行若干次，计算量很大，而且往往难以预测．为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法．具体地说，我们先用欧拉公式求得一个初步的近似值1ny，称之为预测值，预测值1ny的精度可能很差，再用梯形公式2.7将它校正一次，即按2.8式迭代一次得1ny，这个结果称校正值，而这样建立的预测-校正系统通常称为改进的欧拉公式：预测1ny,,nnnyhfxy校正1ny11,,.2nnnnnhyfxyfxy2.13或表为下列平均化形式11,,,,1.2pnnncnnpnpcyyhfxyyyhfxyyyy例2用改进的欧拉方法求解初值问题2.2．解改进的欧拉公式为112,2,1.2npnnnncnppnpcxyyhfyyxyyhfyyyyy仍取0.1h，计算结果见表92．同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度．表72计算结果对比nxnynyxnxnynyx0.10.20.30.40.51.09591.18411.26621.34241.41641.09541.18321.26491.34161.41420.60.70.80.91.01.48601.55251.61531.67821.73791.48321.54921.61651.67331.73217.3龙格-库塔方法7.3.1显式龙格-库塔法的一般形式上