第10章-常微分方程数值解VIP

1、常微分方程与解,为n阶常微分方程。,如果函数在区间a,b内n阶可导，称方程,为方程满足定解条件的解。,第10章常微分方程的数值解,10.1引言,科学研究和工程技术中的问题往往归结为求某个常微分方程的定解问题.常微分方程的理论指出，除少数简单情况能获得初值问题的初等解（用初等函数表示的解）外，绝大多数情况下是求不出初等解的.有些初值问题即便有初等解，也往往由于形式过于复杂而不便处理。常微分方程的数值解法常用来求近似解，由于它提供的算法能通过计算机便捷地实现，因此近年来得到迅速的发展和广泛的应用。,10.2初值问题解法的基本概念,科学技术中常常需要求解常微分方程的定解问题.这类问题最简单的形式，是本章将着重考察的一阶方程的初值问题,我们知道，只有f(x,y)适当光滑譬如关于y满足利普希茨(Lipschitz)条件,理论上就可以保证初值问题的解yf(x)存在并且唯一.我们以下的讨论，都在满足上述条件下进行。,虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法.,所谓数值解法,就是寻求解y(x)在一系列离散节点,上的近似值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时只用到xn+1，xn和yn，即前一步的值。因此，有了初值以后就可以逐步往下计算，其代表是龙格库塔法称为单步法.另一类是用到yn+1前面k点的值yn,yn-1,yn-k+1，称为多步法.其次，要研究公式的局部截断误差和阶，数值解yn与精确解y(xn)的误差估计及收敛性，还有递推公式的计算稳定性等问题.,数值解的思想,（1）将连续变量离散为,（2）用代数的方法求出解函数在点的近似值,如果找不到解函数数学界还关注：解的存在性解的唯一性解的光滑性解的振动性解的周期性解的稳定性解的混沌性,第一步：连续变量离散化,第二步：用直线步进,1、Euler格式,10.3简单单步法,10.3.1欧拉(Euler)方法,过做以为切线斜率的方程,当,时，得,，取,当,时，得,，取,过,做以,为切线斜率的方程,一般地，过,做以,为切线斜率的方程,当,时，得,，取,例1用欧拉公式求解初值问题,解取步长h=0.1，欧拉公式的具体形式为,其中xn=nh=0.1n(n=0,1,10),已知y0=1,由此式可得,依次计算下去，部分计算结果见下表.,与准确解相比，可看出欧拉公式的计算结果精度很差.,欧拉公式具有明显的几何意义,就是用折线近似代替方程的解曲线，因而常称公式(3.1)为欧拉折线法.,还可以通过几何直观来考察欧拉方法的精度.假设yn=y(xn),即顶点Pn落在积分曲线y=y(x)上，那么，,按欧拉方法做出的折线PnPn+1便是y=y(x)过点Pn的切线.从图形上看,这样定出的顶点Pn+1显著地偏离了原来的积分曲线，可见欧拉方法是相当粗糙的.,18世纪最杰出的数学家之一，13岁时入读巴塞尔大学，15岁大学毕业，16岁获得硕士学位。1727年-1741年（20岁-34岁）在彼得堡科学院从事研究工作，在分析学、数论、力学方面均有出色成就，并应俄国政府要求，解决了不少地图学、造船业等实际问题。24岁晋升物理学教授。1735年（28岁）右眼失明。,1741年-1766（34岁-59岁）任德国科学院物理数学所所长，任职25年。在行星运动、刚体运动、热力学、弹道学、人口学、微分方程、曲面微分几何等研究领域均有开创性的工作。1766年应沙皇礼聘重回彼得堡，在1771年（64岁）左眼失明。Euler是数学史上最多产的数学家，平均以每年800页的速度写出创造性论文。他去世后，人们用35年整理出他的研究成果74卷。,14,方法一化导数为差商的方法,由于在逐步求解的过程中，y(xn)的准确值无法求解出来，因此用其近似值代替。为避免混淆，以下学习简记：,y(xn)：待求函数y(x)在xn处的精确函数值yn：待求函数y(x)在xn处的近似函数值,欧拉(Euler)方法（几种推导法）,15,代入初值问题表达式可得：,根据y0可以一步步计算出函数y=y(x)在x1,x2,x3x4,上的近似值y1,y2,y3,y4,常微分方程数值解是一组离散的函数值数据，它的精确表达式很难求解得到，但可以进行插值计算后用插值函数逼近y(x),为了分析计算公式的精度，通常可用泰勒展开将y(xn+1)在xn处展开，则有,在yn=y(xn)的前提下，f(xn,yn)=f(xn,y(xn)=y(xn).于是可得欧拉法(3.1)的公式误差为,称为此方法的局部截断误差.,方法二泰勒级数展开法,17,方法三数值积分法,同样以近似值yn代替精确值y(xn)可得：,将微分方程y=f(x,y)在区间xn,xn+1上积分：,18,2.隐式欧拉法(后退),在数值积分法推导中，积分的近似值取为积分区间宽度与右端点处的函数值乘积，即：,这样便得到了隐式欧拉法：,（3.3）,隐式欧拉公式与欧拉公式有着本质的区别,后者是关于yn+1的一个直接计算公式，这类公式称作是显式的；前者公式的右端含有未知的yn+1，它实际上是关于yn+1的一个函数方程,这类方程称作是隐式的.,显式与隐式两类方法各有特点，考虑到数值稳定性等其他因素，人们有时需要选用隐式方法，但使用显式算法远比隐式方便.,隐式方程通常用迭代法求解，而迭代过程的实质是逐步显式化.,设用欧拉公式,给出迭代初值，用它代入(3.1)式的右端，使之转化为显式，直接计算得,然后再用代入(3.1)式，又有,如此反复进行，得,由于f(x,y)对y满足Lipschitz条件(2.1).由(3.4)减(3.3)得,由此可知，只要hLn)上产生的扰动为，如果：,定义：设在节点xn处用数值算法得到的理想数值解为yn，而实际计算得到的近似解为，称差值：,为第n步的数值解的扰动。,则称该数值方法是稳定的。,下面以欧拉法为例考察计算稳定性.,例4用欧拉公式求解初值问题,解用欧拉法解方程y=-100y得,其准确解是一个按指数曲线衰减很快的函数.,若取步长h=0.025，则欧拉公式的具体形式为,计算结果见表,明显计算过程不稳定,但取h=0.005,yn+1=-1.5yn,则计算过程稳定.,对后退的欧拉公式，取h=0.025时，则计算公式为yn+1=-(1/3.5)yn.计算结果见表,这时计算过程是稳定的.,例题表明稳定性不但与方法有关，也与步长h有关，当然与方程中的f(x,y)有关.为了只考察数值方法本身，通常只检验数值方法用于解模型方程的稳定性，模型方程为,其中为一直实数或复数(Re()0)，这个方程分析较简单，对一般方程可以通过局部线性化化为这种形式。,53,定义6单步法（4.2）用于解模型方程（4.4），若得到的解，满足，则称方法（4.1）是绝对稳定的.,在的平面上，使的变量围成的区域，称为绝对稳定域，,它与实轴的交称为绝对稳定区间.,54,欧拉法：,考察模型方程：,即：,假设在节点值yn上有扰动n，在节点值yn+1上有扰动n+1，且n+1仅由n引起（即：计算过程中不再引起新的误差）,55,欧拉法稳定,即：,欧拉法稳定的条件：,针对模型方程：的显式欧拉法：,化简得：,56,隐式欧拉法：,考察模型方程：,即：,化简为：,假设yn上有扰动，则yn+1的扰动为：,隐式欧拉法稳定,，上式均成立，所以：,隐式欧拉法稳定是恒稳定的,57,故,对有，,故绝对稳定域为的左半平面，,梯形法的稳定性,绝对稳定区间为，即时梯形法均是稳定的.,隐式欧拉法与梯形方法的绝对稳定域均为在具体计算中步长的选择只需考虑计算精度及迭代收敛性要求而不必考虑稳定性，具有这种特点的方法特别重要.,8.3龙格-库塔(Runge-Kutta)法,欧拉方法是显式的一步法,使用方便,但精度较低.本节将构造出高精度的显式一步法：龙格-库塔法,简称R-K法.,8.3.1二阶R-K法,欧拉法的公式为：yi+1=yi+hf(xi,yi)i=0,1,2,n-1决定其精度的是函数f(xi,yi).如能改进这个函数,就可能提高公式的精度.为此把公式改写成：,yi+1=yi+h(xi,yi,h)i=0,1,2,n-1(8.10),选择函数(xi,yi,h),一种方法是用若干个点的函数值的线性组合代替(xi,yi,h),如：,其中cj,aj,bjl是待定参数,aj和bjl满足,以上方法称为p级R-K法,选择cj,aj和bjl,可能使以上方法为p阶方法.显然欧拉法就是一阶R-K法.,二级R-K法的形式是:,此时,由二元函数的泰勒展开：,其中所有的偏导数都是它们在点(xi,y(xi)的值,下同,又由于：,所以,代入(8.11),代入(8.10),而Taylor展开式,二式相减,得局部截断误差,只要c1,c2,a2满足以上方程,就得到一个二阶的R-K法.这是一个不定方程,有无穷多解.比如：,63,结束,(1)取c1=c2=1/2,a2=1得,这实际上是(8.9)公式,即梯形公式的预估-校正公式只迭代一次的形式,通常称为改进的欧拉法.,(2)取c1=0,c2=1,a2=1/2得:,这公式又称中点公式.我们还可以构造其他的二阶R-K法.,8.3.2四阶R-K法,用类似的方法可以确定三级和四级R-K法的参数,构造出三阶和四阶的R-K法.但最常用的是四阶R-K法,四阶R-K法也不只一个,下面给出的是最常用的四阶经典的R-K公式：,65,结束,例3用经典的四阶R-K法计算例2题目,取步长为0.2,且与准确值比较.,计算结果列入表8-3：可见即使用h=0.2计算,也比一阶和二阶方法精度好得多,66,结束,8.4线性多步法,单步法只利用前一步的结果,只要给出初值,就能开始计算但也因为它只利用前一步的值,为了提高精度就要计算一些非结点处的函数值,增加了计算量.R-K法就是通过这一途径提高精度的.下面介绍的线性多步法,在求yi+1时,不仅用到yi的值,还用到前若干步的yi-1,yi-k的值,这些值都是已知的,因此可在计算量增加不多的情况下提高精度.,8.4.1用待定系数法构造线性多步法,线性多步法的一般形式是：,或写为：,其中j,j(j=0,1,k)都是实常数,且k0,0+00,fi+j=f(xi+j,yi+j),j=0,1,k,由(8.15)可看出要计算yi+k,要利用它前面的k个值yi,yi+1,yi+k-1,又因为(8.15)关于yi+j和fi+j都是线性组合,所以这一类方法都称为线性k步法.欧拉法,隐式欧拉法和梯形法都是线性一步法,欧拉中点公式是线性二步法.,把y(x+jh)和y(x+jh)作Taylor展开：,68,结束,其中,69,结束,若选择j,j,使C0=C1=Cp=0,Cp+10,则,将x=xi代入上式,设k=1,并注意到y(xi+jh)=f(xi+jh,y(xi+jh),可推出,即,70,结束,设yi=y(xi),yi+1=y(xi+1),yi+k-1=y(xi+k-1),记左端为yi+k,并舍去最后两项,(8.17)变为：,就是一种p阶的线性k步方法,Cp+1hp+1y(p+1)(xi)称为局部截断误差的主项.当k=0时,是显式方法,当k0时是隐式方法.,下面构造几个实用的线性多步法公式,例4形如：yi+4=-0yi+h(1fi+1+2fi+2+3fi+3)的线性4步法公式,试确定0,1,2,3并求其局部截断误差主项.,71,结束,解:由(8.15)知1=2=3=0,4=1,0=4=0,因为有4个待定系数,由(8.16)写出前4个方程：,解之,得0=-1,1=3=8/3,2=-4/3,故所求的公式为,72,结束,将0,1,2,3代入C4,再求C5,公式(8.19)称为米尔尼(Milne)公式,它的局部截断误差为：,它也可以写成,73,结束,例5试确定下列公式的系数和局部截断误差yi+1=ayi+byi-2+h(cfi+1+dfi+efi-1),解按(8.15)和(8.16)可知k=3,0=-b,1=0,2=-a,3=1,0=0,1=e,2=d,3=c,有五个未知参数,写出前五个方程：,74,结束,解之得：a=9/8,b=-1/8,c=3/8,d=6/8,e=-3/8代入C5,此公式称为哈明(Hamming)公式,写为：,局部截断误差记为：,它是一个四阶隐式三步法.,在(8.15)中若0=1=k-2=0,k=1,则称此类方法为阿达姆斯(Adams)类方法.当k=0时,称显式,否则为隐式.这类方法也可以用以上待定系数法求出这类方法的通式为：,75,结束,例6求(8.24)中的隐式三步法解:k=3,0=1=0,2=-1,3=1,76,结束,解之,0=1/24,1=-5/24,2=19/24,3=9/24.,公式为：,77,结束,局部截断误差为,它是一个线性三步四阶隐式公式,称为四阶阿达姆斯内插公式,应用十分普遍.,此外常用的还有线性四步四阶的显式阿达姆斯公式,常称为阿达姆斯外推公式：,78,结束,8.4.2用数值积分法构造线性多步法公式,在(8.15)式中,若令k=1,0,1,k-1中只有一个不为零,这样的线性多步法公式也可用数值积分方法构造.如阿达姆斯外推公式.,对(8.1)两端在xi,xi+1上积分,选择节点xi-3,xi-2,xi-1,xi作插值节点作函数f(x,y(x)的三次插值多项式：,79,结束,且插值余项,把f(x,y(x)=P3(x)+R3(x)代入(9.29)得：,经计算积分后得,80,结束,所以得,就是阿达姆斯外推公式,其局部截断误差为,用积分中值定理依然可以推出:,81,结束,用类似的方法还可以构造二步四阶的辛普生(Simpson)公式：,它的局部截断误差为,线性多步法要用到二个以上的初始值才能开始计算,不能“自启动”.一般应选择同阶或高阶的单步法算出前几个值,再使用线性多步法.由于线性多步法每一步实际上只计算一次函数值,又能获得较高的精度,所以应用比较广泛.,82,结束,8.5预估校正法,在8.2中介绍过梯形公式的预估校正法,本节中介绍几种高阶的实用的预估校正法.,介绍几个英文缩写：P(Predictor)预估C(Corrector)校正E(Evaluation)求函数值M(Modifier)改进,83,结束,8.5.1阿达姆斯公式的PEC模式,它是用阿达姆斯外推公式作预测,用阿达姆斯内插公式作一次校正,例7用Adams外推公式和AdamsPEC方法计算,取步长h=0.1,前三步用四阶R-K法求值,计算值保留七位有效数字.,84,结束,解:计算结果列表如下(略),8.5.2阿达姆斯公式的PMECME模式,令Pi和Ci分别为(8.34)中的第i步的预测值和校正值：,二式相减,设y(5)(x)变化不大,都约等于y(5)(),85,结束,由此得：,于是在(8.34)式中加上两个改进步(M),最后再求一次函数值(E)供下一次预测用,就形成阿达姆斯法的PMECME模式,86,结束,开始求M4时,无C3和P3,认为C3-P3=0,8.5.3哈明(Hamming)法PMECME模式,利用米尔尼公式(8.20)进行预测,再用哈明公式(8.22)进行校正,中间利用两公式局部截断误差的主项系数的比例进行改进两次,就得到另一种PMECME模式,此模式称哈明法,哈明法为：,开始求M4时,无C3和P3,认为P3-C3=0,87,结束,8.6一阶常微分方程组和高阶方程,本章8.1.2已提到,一阶