数值分析第四版第九章常微分方程初值问题数值解法

第九章常微分方程初值问题的数值解法,科学研究和工程技术中的问题往往归结为求某个常微分方程的定解问题.常微分方程的理论指出，很多方程的定解虽然存在，但可能十分复杂难于计算，也可能不能用简单的初等函数表示,因此常求其能满足精度要求的近似解.,1引言,常微分方程的数值解法常用来求近似解，由于它提供的算法能通过计算机便捷地实现，因此近年来得到迅速的发展和广泛的应用常微分方程数值解法的特点是：对求解区间进行剖分，然后把微分方程离散成在节点上的近似公式或近似方程，最后结合定解条件求出近似解.因此数值解法得到的近似解是一个离散的函数表.,1,设方程问题的解y(x)的存在，令a=x0x1xn其中hk=xk+1-xk,如是等距节点xk=a+kh,h称为步长。y(x)的解析表达式不容易得到或根本无法得到，我们用数值方法求得y(x)在每个节点xk上y(xk)的近似值，用yk表示，即yky(xk)，这样y0,y1,.,yn称为微分方程的数值解。,如方法中要计算yn+1时只用到前一点的值yn则称为单步法.又如需要前k点值yn,yn-1,yn-k+1来计算yn+1称其为k步法.其次我们还得研究算法的局部误差和阶,数值解与精确解间的误差及收敛性,还有数值算法的稳定性.,微分方程的数值解：,2,基本知识复习一阶常微方程的初值问题：,其中f(x,y)是已知函数，(2)是定解条件，常微分方程的理论指出：当f(x,y)定义在区域G=(axb,y)，若存在正的常数L，使：f(x,y1)f(x,y2)Ly1-y2|,对任意的xa,b和任意的y1,y2成立，则称f(x,y)对y满足利普希茨(Lipschitz)条件,L称利普希茨常数.若f(x,y)在区域G连续，关于y满足利普希茨条件，则上述一阶常微分方程的初值问题的解存在且唯一.我们以下的讨论，都在满足上述条件下进行.,3,一阶常微分方程组和高阶常微分方程,一阶常微分方程组常表述为：,写成向量形式为,4,还有高阶常微分方程定解问题，如二阶定解问题：,可化作一阶常微方程组的定解问题：,本章着重讨论一阶常微方程初值问题的各种数值解法.这些解法都可以写成向量形式而用于一阶常微分方程组的初值问题.这样，也就解决了高阶方程的定解问题.,5,返回,2简单的数值方法与基本概念,对问题,最简单而直观的方法是欧拉方法.欧拉方法在精度要求不高时，仍不失为一实用方法.同时通过欧拉方法的讨论，容易弄清常微方程初值问题数值解法的一些基本概念和构造方法的思路.,欧拉方法(Euler)的导出,取步长为h，节点xn=x0+nh,n=0,1,2,，其中x0=a又设y(x)为上述问题的解.,对上式两端积分，积分区间取xn,xn+1,6,记g(x)=f(x,y(x)对上述(*)积分用左矩形公式,得到：,舍去Rn，并令yn=y(xn),得到：,这就是欧拉公式.,7,欧拉公式还可由泰勒展开得到.假定y(x)二阶连续可导，把y(xn+1)在xn点展开：,舍去h2项，并令yn=y(xn)，注意到y(xn)=f(xn,yn),也得到,欧拉公式有明显的几何意义，如图9-1(P337)，过点(x0,y0)的曲线是解y(x).欧拉方法是在(x0,y0)作y(x)的切线，它与直线x=x1交于(x1,y1)，过(x1,y1)作过此点的积分曲线的切线，又与x=x2交于(x2,y2)，如此下去，得到一条折线，欧拉方法就是用这条折线近似地代替曲线y(x)，故欧拉方法有时也称欧拉折线法.,8,推出：,不难看出，(2.5)式中公式两端都含有yn+1,一般情况不能由yn的值计算yn+1,而需要解方程，故称为欧拉隐式公式,隐式方法一般用迭代方式求解:设由欧拉方法得:,代如(2,5)式得:,如此有:,对(*)式的右端用右矩形公式代替积分，可得到,欧拉隐式公式和欧拉中点公式,(后退的欧拉公式),9,推出：,称为欧拉中点公式(二步法).容易看出,中点公式计算yn+1时,不仅需要yn的值，还需要yn-1的值.,又如把积分区间改为xn-1,xn+1，并对右端用中矩形公式代替积分，可得到：,利用利普希茨条件:由(2.5)减(2.6)式得:,由此只要hLn)上产生的偏差均不超过n,即mn则称该方法是稳定的.否则称此方法是不稳定的.,39,定义:单步法(4.1)用于试验方程(4.8),如果得到的解yn+1=E(h)yn,满足E(h)1,则称方法(4.1)是绝对稳定的.在=h的复平面上,使E(h)1的变量围成的区域,称为绝对稳定区域,它与实轴的交称为绝对稳定区间.,单步法(4.1)用于试验方程(4.8)得到的解形式yn+1=E(h)yn故得扰动方程n+1=E(h)n,Euler法的绝对稳定区域,从而,绝对稳定区域越大,h可选大些,方法的适应性就越强.,后退的Euler方法稳定性,得,只要R(h)0,则有E(h)1,故后退的Euler方法是无条件稳定性的,但收敛性要求0h1,h仍受限制.,41,梯形公式的稳定性,所以梯形公式是无条件绝对收敛的.,42,四阶R-K方法的绝对稳定区域,43,所以,所以绝对稳定区域为:,44,返回,45,5线性多步法,单步法只利用前一步的结果,只要给出初值,就能开始计算但也因为它只利用前一步的值,为了提高精度就要计算一些非结点处的函数值,增加了计算量.R-K法就是通过这一途径提高精度的.下面介绍的线性多步法,在求yn+1时,不仅用到yn的值,还用到前若干步的yn-1,yn-k的值,这些值都是已知的,因此可在计算量增加不多的情况下提高精度.形成线性多步法.,用待定系数法构造线性多步法,线性k步法的一般形式是：,46,其中i,i都是常数,0+00,yn+i为y(xn+i)近似值,fn+j=f(xn+j,yn+j),xn+i=x0+ih,由(5.1)可看出要计算yn+k,要利用它前面的k个值yn,yn+1,yn+k-1,又因为(5.1)关于yn+j和fn+j都是线性组合,所以这一类方法都称为线性k步法.如果k=0,则称显式k步法;如果k0,则称隐式k步法,定义:设y(x)是初值问题的准确解,线性多步法在xn+k上的局部误差为,若,称(5.1)是p阶方法.p1,则称(5.1)方法与初值问题(1.1)是相容的.,把y(xn+ih)和y(xn+ih)作Taylor展开：,代入(5.1)得,其中,满足,47,q=2,3,若选择j,j,使C0=C1=Cp=0,Cp+10,则就是一种p阶的线性k步方法,Tn+k=Cp+1hp+1y(p+1)(xn)+O(hp+2)称第一项为局部截断误差的主项,Cp+1称为误差常数.,由相容性定义,p1,即有:C0=C1=0得:,所以(5.1)与微分(1.1)相容的充要条件是(5.6)成立.,48,当k=1时,若1=1则由(5.6)得:0=0=1,此时即为欧拉法,由(5.4)可得C2=1/20所以方法为1阶精度.局部误差为:,对k=1时,若10此时方法为隐式公式.为了确定0,0,1可由C0=C1=C2=0解得0=1,0=1=1/2,得梯形法,由(5.4)得C3=-1/120故p=2,为二阶方法,对k2,的多步法可由(5.4)确定j,j,且由(5.5)可得局部误差.,49,阿当姆斯显式与隐式公式,形如:,的k步法,称为阿当姆斯(Adams)方法.k=0为显式方法,k0为隐式方法,称为阿当姆斯显式与隐式公式.也称Adams-Bashforth公式和Adams-Monlton公式.这可由方程(1.1)两端积分求得.,对(1.1)两端在xn+k-1,xn+k上积分,50,用节点xn+k,xn+k-1,xn,或xn+k-1,xn-1的F(x)的k次插值多项式(x)代替F(x)求积分,即得线性多步公式.,选择节点xn+3,xn+2,xn+1,xn作插值节点作函数f(x,y(x)的三次插值多项式：,可得,就是四步4阶阿达姆斯(显式)公式,其局部截断误差为,51,也可由(5.4)中C0=C1=Cp=0推出.对比(5.1)与(5.7)式,得:此时0=1=k-2=0,k-1=1.由(5.4)显有C0=0.令:C1=Ck+1=0则可得0,1,.k(若k=0,则令C0=C1=Ck=0来求得0,1,.k-1),例如k=3时,由C1=C2=C3=C4=0,根据(5.4)得:,52,若3=0(显式),则由前三个方程解得:0=5/121=-16/122=23/12可得k=3时的阿当姆斯显式公式,再由(5.4)可得:C4=3/8,所以(5.8)是三阶方法.且,若30(隐式),则由四个方程解得:0=1/241=-5/242=19/243=3/8可得k=3时的阿当姆斯隐式公式,53,(5.9)是四阶方法.局部误差是,类似可以求得常用的阿当姆斯显式或隐式公式,kp阿当姆斯显式公式(P364表9-5)Cp+1,54,kp阿当姆斯隐式公式(P364表9-6)Cp+1,55,米尔尼(Milne)方法与辛普森(Simpson)方法,考虑另一个k=4的显式公式,其中i是待定常数,可根据使公式的阶尽可能的高这一条件来确定其数值.,根据(5.4)得,C0=0令C1=C2=C3=C4=0,得:,可得:0=01=8/32=-4/33=8/3,56,57,故所求的公式为,称为米尔尼(Milne)公式,为4阶方法,它的局部截断误差为：,米尔尼(Milne)方法也可用积分得:,类似还可以构造二步四阶的辛普森(Simpson)公式：,局部截断误差为,右端积分用辛普森求积公式可得辛普森方法:,汉明方法辛普森公式是二步法中阶数最高的,但它的稳定性较差.可考虑另一类三步法公式,58,0,1,2,1,2,3为常数,为使公式为四阶的,若取1,=1,则可的辛普森公式,若取1,=0,利用泰勒展开由(5.4),令C0=C1=C2=C3=C4=0得,解之得：0=-1/82,=9/81=-3/82=6/83=3/8,得称为汉明(Hamming)公式,59,60,局部截断误差为：,它是一个四阶隐式三步法.,且可知C5=-1/40,线性多步法要用到二个以上的初始值才能开始计算,不能“自启动”.一般应选择同阶或高阶的单步法算出前几个值,再使用线性多步法.由于线性多步法每一步实际上只计算一次函数值,又能获得较高的精度,所以应用比较广泛.,61,例隐式三步法,解:k=3,0=1=0,2=-1,3=1,解之,0=1/24,1=-5/24,2=19/24,3=9/24.,公式为：,62,63,局部截断误差为,它是一个线性三步四阶隐式公式,称为四阶阿达姆斯内插公式,应用十分普遍.,此外常用的还有线性四步四阶的显式阿达姆斯公式,常称为阿达姆斯外推公式：,64,预估校正法,对于隐式线性多步法,计算时要进行迭代,计算量大,通常采用显式方法先给出yn+k的初始近似,记为ypn+k称为P(Predictor)预估,再计算fn+k的函数值E(Evaluation),继而以隐式公式计算yn+k称为C(Corrector)校正及M(Modifier)改进,阿达姆斯四阶预测-校正公式的PECE模式,65,阿达姆斯公式的PMECME模式,由四阶阿达姆斯公式的截断误差,对于PECE有,二式相减,66,由此得事后误差估计：,于是加上两个改进步(M),最后再求一次函数值(E)供下一次预测用,就形成阿达姆斯法的PMECME模式,67,修正米尔尼-汉明预测-校正法(PMECME),利用米尔尼公式(5.11)进行预测,再用哈明公式(5.15)进行校正,中间利用两公式局部截断误差的主项系数的比例进行改进两次,就得到PMECME模式,此模式称修正米尔尼-汉明预测-校正法,有时侯并不一定要按系数公式(5.4)确定多步法(5.1)的系数i,i,可直接用泰勒展开进行构造.,例解初值问题的显式二步法,试确定系数i,i使方法阶数尽量的高,并求局部截断误差,解:,68,69,为使方法阶数尽量高,令,解得:,此时为三阶公式,且所求局部误差为:,例证明存在的一个值,使线性多步法是四阶的,解:只需证明局部误差Tn+1=O(h5)由泰勒展开可得:,70,71,当=9时,Tn+1=O(h5).故方法是四阶的.,72,返回,6一阶常微分方程组和高阶方程(P373),一阶常微分方