数值分析--第9章_常微分方程初值问题数值解法

上页下页第9章常微分方程初值问题数值解法•9.1引言•9.2简单的数值方法与基本概念•9.3龙格-库塔方法•9.4单步法的收敛性与稳定性•9.5线性多步法•9.6方程组和高阶方程上页下页9.1引言科学技术中常常需要求解常微分方程的定解问题.这类问题最简单的形式，是本章将着重考察的一阶方程的初值问题)2.1(.)()1.1(),,(00yxyyxfy我们知道，只有f(x,y)适当光滑—譬如关于y满足利普希茨(Lipschitz)条件理论上就可以保证初值问题的解y＝f(x)存在并且唯一..),(),(yyLyxfyxf上页下页虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法.所谓数值解法,就是寻求解y(x)在一系列离散节点121nnxxxx上的近似值y1,y2,,yn,yn+1,.相邻两个节点的间距hn=xn+1-xn称为步长.今后如不特别说明，总是假定hi=h(i=1,2,)为定数,这时节点为xn=x0+nh(i=0,1,2,)(等距节点).上页下页初值问题的数值解法有个基本特点，他们都采取“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进.描述这类算法，只要给出用已知信息yn,yn-1,yn-2,计算yn+1的递推公式.首先，要对微分方程离散化，建立求解数值解的递推公式.一类是计算yn+1时只用到前一点的值yn，称为单步法.另一类是用到yn+1前面k点的值yn,yn-1,,yn-k+1，称为k步法.其次，要研究公式的局部截断误差和阶，数值解yn与精确解y(xn)的误差估计及收敛性，还有递推公式的计算稳定性等问题.上页下页9.2简单的数值方法与基本概念9.2.1欧拉法与后退欧拉法我们知道，在xy平面上，微分方程(1.1)式的解y=f(x)称作它的积分曲线，积分曲线上一点(x,y)的切线斜率等于函数f(x,y)的值.如果按f(x,y)在xy平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致.基于上述几何解释，我们从初始点P0(x0,y0)出发,先依方向场在该点的方向推进到x=x1上一点P1，然后再从P1点依方向场在该点的方向推进到x=x2上一点P2,循环前进做出一条折线P0P1P2.上页下页一般地，设已做出该折线的顶点Pn,过Pn(xn,yn)依方向场的方向再推进到Pn+1(xn+1,yn+1)，显然两个顶点Pn,Pn+1的坐标有关系斜率),,(111nnnnnnnnyxfhyyxxyy)1.2(),(1nnnnyxhfyy这就是著名的(显式)欧拉(Euler)公式.若初值y0已知，则依公式(2.1)可逐次逐步算出各点数值解.即),,(0001yxhfyy),,(1112yxhfyy上页下页例1用欧拉公式求解初值问题)2.2(.1)0(),10(2yxyxyy解取步长h=0.1，欧拉公式的具体形式为)2(1nnnnnyxyhyy其中xn=nh=0.1n(n=0,1,,10),已知y0=1,由此式可得191818.1)1.12.01.1(1.01.1)2(1.11.01)2(1111200001yxyhyyyxyhyy上页下页依次计算下去，部分计算结果见下表.xy21与准确解相比，可看出欧拉公式的计算结果精度很差.xn欧拉公式数值解yn准确解y(xn)误差0.20.40.60.81.01.1918181.3582131.5089661.6497831.7847701.1832161.3416411.4832401.6124521.7320510.0086020.0165720.0257260.0373310.052719上页下页欧拉公式具有明显的几何意义,就是用折线近似代替方程的解曲线，因而常称公式(2.1)为欧拉折线法.()yyxxynx1nxnp1np1npx还可以通过几何直观来考察欧拉方法的精度.假设yn=y(xn),即顶点Pn落在积分曲线y=y(x)上，那么，按欧拉方法做出的折线PnPn+1便是y=y(x)过点Pn的切线.从图形上看,这样定出的顶点Pn+1显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的.上页下页为了分析计算公式的精度，通常可用泰勒展开将y(xn+1)在xn处展开，则有).,()(2),()()(2)()()()(1221nnnnnnnnnnnnxxyhyxhfxyyhxyhxyhxyxy在yn=y(xn)的前提下，f(xn,yn)=f(xn,y(xn))=y(xn).于是可得欧拉法(2.1)的公式误差为)3.2(),(2)(2)(2211nnnnxyhyhyxy称为此方法的局部截断误差.上页下页如果对方程(1.1)从xn到xn+1积分，得)4.2(.))(,()()(11nnxxnndttytfxyxy右端积分用左矩形公式hf(xn,y(xn))近似，再以yn代替y(xn)，yn+1代替y(xn+1)也得到欧拉公式(2.1)，局部截断误差也是(2.3).称为(隐式)后退的欧拉公式.如果右端积分用右矩形公式hf(xn+1,y(xn+1))近似，则得到另一个公式)5.2(),,(111nnnnyxhfyy上页下页后退的欧拉公式与欧拉公式有着本质的区别,后者是关于yn+1的一个直接计算公式，这类公式称作是显式的；前者公式的右端含有未知的yn+1，它实际上是关于yn+1的一个函数方程,这类方程称作是隐式的.显式与隐式两类方法各有特点，考了到数值稳定性等其他因素，人们有时需要选用隐式方法，但使用显式算法远比隐式方便.隐式方程通常用迭代法求解，而迭代过程的实质是逐步显式化.上页下页设用欧拉公式),()0(1nnnnyxhfyy给出迭代初值，用它代入(2.5)式的右端，使之转化为显式，直接计算得)0(1ny),,()0(11)1(1nnnnyxhfyy然后再用代入(2.5)式，又有)1(1ny).,()1(11)2(1nnnnyxhfyy如此反复进行，得)6.2().,1,0(),()(11)1(1kyxhfyyknnnkn上页下页由于f(x,y)对y满足Lipschitz条件(1.3).由(2.6)减(2.5)得),(),(11)(111)1(1nnknnnknyxfyxfhyy.1)(1nknyyhL由此可知，只要hL1，迭代法(2.6)就收敛到解.关于后退欧拉方法的公式误差，从积分公式看到它与欧拉法是相似的.上页下页9.2.2梯形方法为得到比欧拉法精度高的计算公式，在等式(2.4)右端积分用梯形求积公式近似,并用yn代替y(xn),yn+1代替y(xn+1)，则得)7.2(,),(),(2111nnnnnnyxfyxfhyy称为矩形方法.矩形方法是隐式单步法，用迭代法求解，同后退的欧拉方法一样，仍用欧拉法提供迭代初值，则矩形迭代公式为上页下页)8.2().,1,0(),(),(2);,()(11)1()0(11kyxfyxfhyyyxhfyyknnnnnknnnnn为了分析迭代过程的收敛性,将(2.7)与(2.8)相减,得)],,(),([2)(1111)1(11knnnnknnyxfyxfhyy于是有,2)(11)1(11knnknnyyhLyy使得,12hL则当k→∞时有,这说明迭代过程(2.8)是收敛的.1)(1nknyy上页下页9.2.3单步法的局部截断误差与阶初值问题(1.1),(1.2)的单步法可用一般形式表示为)9.2(,),,,(11hyyxhyynnnnn其中多元函数与f(x,y)有关，当含有yn+1时，方法是隐式的，若不含yn+1则为显式方法，所以显式单步法可表示为(x,y,h)称为增量函数，例如对欧拉法(2.1)有)10.2(,),,(1hyxhyynnnn.),(),,(yxfhyxnn它的局部截断误差(2.3)已由给出,对一般显式单步法则可如下定义.上页下页定义1设y(x)是初值问题(1.1),(1.2)的准确解,称)11.2().),(,(11hxyxhyyTnnnnn为显式单步法(2.10)的局部截断误差.Tn+1之所以称为局部的，是假设在xn前各步没有误差.当yn=y(xn)时，计算一步，则有.)),(,()()()],,([)()(11111nnnnnnnnnnnThxyxhxyxyhyxhyxyyxy所以，局部截断误差可理解为用方法(2.10)计算一步的误差，也即公式(2.10)中用准确解y(x)代替数值解产生的公式误差.根据定义,显然欧拉法的局部截断误差为上页下页).()(2)()()())(,()()(3211hOxyhxyhxyhxyxyxhfxyxyTnnnnnnnnn即为(2.3)的结果.这里称为局部截断误差主项.显然Tn+1=O(h2).一般情形的定义如下)(22nxyh定义2设y(x)是初值问题的准确解，若存在最大整数p使显式单步法(2.10)的局部截断误差满足)12.2().(),,()()(11pnhOhyxhxyhxyT则称方法(2.10)具有p阶精度.上页下页若将(2.10)展开式写成).()(2)]()()([)()(2)())(,()()(322321111hOxyhhOxyhxyhhOxyhxyhxyxhfxyxyTnnnnnnnnnn).())(,(211ppnnnhOhxyxT则称为局部截断误差主项.1))(,(pnnhxyx以上定义对隐式单步法(2.9)也是适用的.例如，对后退欧拉法(2.5)其局部截断误差为这里p=1是1阶方法，局部截断误差主项为).(22nxyh上页下页同样对梯形法(2.7)有).()(12)()](2)()()([2)(!3)(2)()]()([2)()(434232111hOxyhhOxyhxyhxyxyhxyhxyhxyhxyxyhxyxyTnnnnnnnnnnnnn所以梯形方法(2.7)是二阶的.其局部截断误差主项为).(124nxyh上页下页9.2.4改进的欧拉公式我们看到，梯形方法虽然提高了精度，但其算法复杂，在应用迭代公式(2.9)进行实际计算时，每迭代一次，都要重新计算函数f(x,y)的值，而迭代又要反复进行若干次，计算量很大，而且往往难以预测.为了控制计算量，通常只迭代一两次就转入下一步的计算，这就简化了算法.具体地说，我们先用欧拉公式求得一个初步的近似值，称之为预测值，此预测值的精度可能很差，再用梯形公式(2.7)将它校正一次，即按(2.8)式迭代一次，这个结果称之为校正值.1ny1ny上页下页这样建立的预测—校正系统通常称为改进的欧拉公式:或表为下列平均化形式).(21),,(),,(11cpnpnncnnnpyyyyxhfyyyxhfyy(2.13),),(1nnnnyxhfyy预测.)],(),([2111nnnnnnyxfyxfhyy校正上页下页例2用改进的欧拉法解例1中的初值问题(2.2).解仍取步长h=0.1,改进的欧拉公式为).(21),2(),2(11cpnpnpncnnnnpyyyyxyhyyyxyhyy部分计算结果见下表xnyn误差xnyn误差00.20.41.1.1840961.34336