数值分析 第9章 常微分方程数值解法

1第9章常微分方程初值问题数值解法/*NumericalMethodforOrdinaryDifferentialEquations*/9.1引言9.2简单的数值方法9.3龙格-库塔方法9.4单步法的收敛性与稳定性9.5线性多步法9.6线性多步法的收敛性与稳定性9.7一阶方程组与刚性方程29.1引言待求解的问题.)(],,[),,(000yxybxxyxfy（1.1）（1.2）R,,.),(),(212121yyyyLyxfyxf（1.3）即存在实数，使得0L称为的利普希茨常数（简称Lips.常数）.Lf解的存在唯一性定理1设在区域上连续，关于满足利普希茨条件，fR},),({ybxayxDy/*Initial-ValueProblem*/:)(xyR],,[00ybax],[bax则对任意，常微分方程(1.1),(1.2)式当时存在唯一的连续可微解.一阶常微分方程初值问题3解对扰动的敏感性结论定理2设在区域（如定理1所定义）上连续，且关于满足利普希茨条件，设初值问题fDysxyyxfy)(),,(0的解为，则),(sxy.),(),(21210ssesxysxyxxL这个定理表明解对初值的敏感性，它与右端函数有关，当的Lips.常数比较小时，解对初值和右端函数相对不敏感，可视为好条件.若较大则可认为坏条件，即为病态问题.LffL4如果右端函数可导，由中值定理有.,,),(),(),(212121之间在yyyyyxfyxyyxy若假定在域内有界，设，则Dyyxf),(Lyyxf),(.),(),(2121yyLyxyyxy它表明满足利普希茨条件，且的大小反映了右端函数关于变化的快慢，刻画了初值问题(1.1)和(1.2)式是否是好条件.fLfy5解析解法:如何求解计算解函数y(x)在一系列节点a=x0x1…xn=b处的近似值),...,1()(nixyyii节点间距为步长，通常采用等距节点，即取hi=h(常数)。)1,...,0(1nixxhiii求解所有的常微分方程步进式：根据已知的或已求出的节点上的函数值计算当前节点上的函数值，一步一步向前推进。因此只需建立由已知的或已求出的节点上的函数值求当前节点函数值的递推公式即可。常数变易法、Lapalace变换等分离变量法、变量代换、数值解法:678一类是计算时只用到前一点的值，称为单步法.1nyny另一类是用到前面点的值，1nyk11,,,knnnyyy称为步法.k其次:要研究公式的局部截断误差和阶，数值解与精确解的误差估计及收敛性，还有递推公式的计算稳定性等问题.ny)(nxy首先:对方程离散化，建立求数值解的递推公式.),(yxfy9欧拉（Euler）方法是解初值问题的最简单的数值方法。9.2简单的数值方法9.2.1欧拉法与后退欧拉法00)(),(yxyyxfy初值问题),(00yx的解y=y(x)代表通过点的一条称之为微分方程的积分曲线。),(yx)(xy),(yxf积分曲线上每一点的切线的斜率等于函数在这点的值。321!3)(!2)()()()()(hxyhxyhxyxyhxyxynnnnnn导数用均差近似，即))(,()()()(1nnnnnxyxfxyhxyxyEuler方法的导出1然后用代替，即得ny()nyx).,(1nnnnyxhfyy（2.1）,2,1,0n10从初始点P0(即点(x0,y0))出发,作积分曲线y=y(x)在P0点上切线(其斜率为),10PP),()(000yxfxyPi+1Pny=y(x)P1PiPnPi+1P0x0x1xixi+1xnPiP1若初值已知，则依公式可逐步算出0y),,(0001yxhfyy),,(1112yxhfyy(即点(x1,y1),得到y1作为y(x1)的近似值,如上图所示。))(,(0000xxyxfyy过点(x0,y0),以f(x0,y0)为斜率的切线方程为),(0001yxhfyy获得了P1点的坐标。1xx当时,得Euler法的求解过程与x=x1直线相交于P1点).,(1nnnnyxhfyy11同样,过点P1(x1,y1),作))(,(1111xxyxfyy),(1112yxhfyy当时,得2xxPi+1Pny=y(x)P1PiPnPi+1P0x0x1xixi+1xnPiP1由此获得了P2的坐标。重复以上过程,就可获得一系列的点:P1,P1,…,Pn。对已求得点，以=为斜率作直线),(nnnyxP)(nxy),(nnyxf))(,(nnnnxxyxfyy当时,得1nxx),(1nnnnyxhfyynnyxy)(取获得了一条近似于曲线y=y(x)的折线这样,从x0逐个算出对应的数值解nxxx,,21nyyy,,21nPPPP321),,(11nnnnnnyxfxxyy交直线x=x2于P2点,),(11yxf斜率=12例1用欧拉公式求解初值问题.1)0(),10(2yxyxyy（2.2）欧拉公式为).2(1nnnnnyxyhyy取步长，计算结果见表9-1.1.0hxy21)(nxyny解7321.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计算结果对比表欧拉方法的精度很差),2(),(yxyyxf13误差的几何直观假设，即顶点落在积分曲线上，那么，按欧拉方法做出的折线便是过点的切线（图9-2）.)(nnxyynP)(xyy1nnPP)(xyynP图9-2将)(1nxy在处展开，则有nx)()(1hxyxynn)()(2)()(32hOxyhhxyxynnnEuler法的误差分析在的前提下，)(nnxyy)())(,(),(nnnnnxyxyxfyxf11)(nnyxy误差（2.3）),(22nxyh显然，这种近似有一定误差，而且步长越大，误差越大，如何估计这种误差y(xn+1)yn+1？在假设yn=y(xn)，即第n步计算是精确的前提下，考虑公式或方法本身带来的误差:Rn+1=y(xn+1)yn+1,称为局部截断误差/*localtruncationerror*/。14截断误差:实际上，y(xn)yn，yn也有误差，它对yn+1的误差也有影响，见下图。但这里不考虑此误差的影响，仅考虑方法或公式本身带来的误差，因此称为方法误差或截断误差。局部截断误差的分析：由于假设yn=y(xn)，即yn准确，因此分析局部截断误差时将y(xn+1)和yn+1都用点xn上的信息来表示，工具：Taylor展开。欧拉法的局部截断误差：Tn+1的主项/*leadingterm*/定义若某算法的局部截断误差为O(hp+1)，则称该算法有p阶精度。欧拉法具有1阶精度。)()(2321hOxyhTnn15.))(,()()(11nnxxnndttytfxyxy（2.4）对方程从到积分，得nx1nx),(yxfy左矩形公式近似))(,(nnxyxfh再以代替ny),(nxy代替也得到欧拉法)(1nxy1ny后退的Euler法如果在（2.4）中右端积分用右矩形公式))(,(11nnxyxfh),,(111nnnnyxhfyy（2.5）近似，))(,()()()(11111nnnnnnnxyxfxyxxxyxyEuler方法的导出2).,(111nnnnyxhfyy,2,1,0n),(1nnnnyxhfyy,2,1,0n显式隐式16逐步显示化隐式公式不能直接求解，一般需要用Euler显式公式得到初值，然后用Euler隐式公式迭代求解。因此隐式公式较显式公式计算复杂，但稳定性好（后面分析）。).,1,0(),,(),()(11)1(1)0(1kyxfhyyyxfhyyknnnknnnnn（2.6）迭代的收敛性由于对满足利普希茨条件),(yxfy),(),(11)(111)1(1nnknnnknyxfyxfhyy1)(1nknyyhL由此可知，只要迭代法(2.6)就收敛到解.1hL1ny迭代公式两端取极限即得.)(1)0(11nnkyyhL).,(111nnnnyxhfyy17后退欧拉法的局部误差))(,()()(1111nnnnnxyxhfxyxyT).()(232hOxyhn),,(111nnnnyxhfyy这里，是1阶方法，局部截断误差主项为.1p)(22nxyh)()(2)()(32hOxyhxyhxynnn)]()()([2hOxyhxyhnn18分别用代替得到计算公式1,nnyy),(),(1nnxyxy)],,(),([2111nnnnnnyxfyxfhyy（2.7）梯形方法隐式单步法，用迭代法求解.)],(),([2))(,(111nnnnxxyxfyxfhdttytfnn1))(,()()(1nnxxnndttytfxyxy);,()0(1nnnnyxfhyy)],(),([2)(11)1(1knnnnnknyxfyxfhyy（2.8）).,2,1,0(k迭代公式9.2.2梯形方法19)],,(),([2)(1111)1(11knnnnknnyxfyxfhyy迭代公式的收敛性)],,(),([2111nnnnnnyxfyxfhyy)],(),([2)(11)1(1knnnnnknyxfyxfhyy如果选取充分小，使得h,12hL则当时有，k1)(1nknyy迭代过程(2.8)是收敛的.)(11)1(112knnknnyyhLyy式中为关于的利普希茨常数.),(yxfyL.)2()0(111nnkyyhL20梯形法的局部误差估计)]()([2)()(111nnnnnxyxyhxyxyT).()(1243hOxyhn)],,(),([2111nnnnnnyxfyxfhyy梯形方法是2阶的，其局部误差主项为)(123nxyh)()()(!3)(2)()(432nnnnnxyhOxyhxyhxyhxy)()](2)()()([242hOxyhxyhxyxyhnnnn21)],,(),([2111nnnnnnyxfyxfhyy迭代一两次就转入下一步的计算.预测值精度可能很差改进的欧拉公式校正一次1ny迭代一次得，这个结果称校正值.)],(),([2)(11)1(1knnnnnknyxfyxfhyy迭代计算量大Step1:先用显式欧拉公式作预测，算出),(1nnnnyxfhyyStep2:再将代入隐式梯形公式1ny)],(),([2111nnnnnnyxfyxfhyy9.2.3改进欧拉公式梯形法的缺点改进22预测-校正系统：预测),,(1nnnnyxfhyy校正)].,(),([2111nnnnnnyxfyxfhyy（2.9）平均化形式),,(nnnpyxfhyy),,(1pnncyxfhyy)