数值计算方法-第9章-常微分方程数值解

1,第9章常微分方程数值解法,复旦大学力学与工程科学系,9-2,很多实际问题往往归结为一组微分方程,但它们的解析解很难得到,往往采用数值方法得到近似解.本章主要讨论常微分方程(组)初值问题的数值解法.,本章的目的,概述,问题的提法,常微分方程初值问题的提法,假设连续且关于y满足李谱希兹(Lipschitz)条件:,L为给定的常数.,根据常微分方程理论,上面的方程一定存在唯一的连续可微解.,求方程的解y=y(x).,(1),复旦大学力学与工程科学系,9-3,边值问题的适定性,在许多实际问题中,是由观测得到的,存在一定的误差.如果它们的误差是微小的,那么能否保证解的误差也是微小的.,定义:假设初值问题(1)有唯一解.如果存在正常数,使得对任何,当,以及,摄动的初值问题:,有唯一解z(x),并且满足:,则称初值问题(1)是适定的.,定理:若连续,且关于y满足Lipschitz条件,则初值问题(1)是适定的.,复旦大学力学与工程科学系,9-4,离散变量法及离散误差,离散变量法,在若干个离散点:a=x0x1xn=b上,求出函数y(x)的近似值yi(i=1,2,n).,称yi为原问题的数值解(即用它来近似y(xi)的值).hi=xi+1-xi为xi到xi+1的步长.一般取步长为常数h.,要求数值解,必须对微分方程进行离散.常用方法有以下几种:,用差商近似导数,如用向前差商代替导数,用y(xn)的近似值yn代入上式,可得y(xn+1)的近似值yn+1:,原方程的近似解:,也可用其它格式的差商代替导数.,复旦大学力学与工程科学系,9-5,用数值积分方法,用左端点的矩形公式计算积分,得y(xn+1)的近似值yn+1,用梯形公式计算积分,用yn近似y(xn),yn+1近似y(xn+1),得,复旦大学力学与工程科学系,9-6,用Talor展开近似的方法,用y(xn)的近似值yn代入上式,可得y(xn+1)的近似值yn+1:,上面三种方法都能导出常微分方程数值解的公式,其中Talor展开方法还可以用来估计截断误差,所以推导时都用该方法.,令x=xn,则,若取p=1,得,复旦大学力学与工程科学系,9-7,一.欧拉(Euler)方法,Euler方法的各种形式,Euler方法,称为Euler方法,近似解是通过(x0,y0)的一条折线,每个折线段的方向与左端点处f(x)的切线方向一致.故Euler方法又称为Euler折线法.,Euler方法的几何解释,Euler方法是显式的,可直接递推求解.,复旦大学力学与工程科学系,9-8,向后Euler方法,这是隐式公式,一般用迭代法求解:,若用向后差商代替导数,即:,用y(xn)的近似值yn代入上式,可得y(xn+1)的近似值yn+1:,中心Euler方法,若用中心差商代替导数,即:,因为要用到前面两步的结果yn-1,yn,故又称为Euler两步公式.,复旦大学力学与工程科学系,9-9,Euler方法的算法描述,1.输入a,b,h,f(x,y),初值y0;2.n0,xa,yy0;3.输出n,(x,y);4.nn+1,yy+hf(x,y);5.xx+h;若x0,使得,整体截断误差:,由闭区间上连续函数的性质,Rn+1在a,b区间有界.,复旦大学力学与工程科学系,9-14,Euler方法的整体截断误差(续),反复递推得,因为在a,b中求解,故有:,利用不等式:,Euler方法的整体截断误差与h同阶.当h0时,eN0,复旦大学力学与工程科学系,9-15,二.改进的欧拉(Euler)方法,梯形公式,梯形公式的导出,梯形公式的误差,取a=xn,b=xn+1,得梯形方法的局部误差:,用梯形积分公式在xn,xn+1上求积分:,用近似值yn代替y(xn),yn+1代替y(xn+1),可得:,由定理7.2梯形积分公式的误差:,复旦大学力学与工程科学系,9-16,梯形公式的误差(续),梯形公式是隐式格式,通常采用迭代方法求解.,梯形公式的迭代格式,梯形公式的收敛性,由于f(x,y)关于y满足Lipschitz条件,若,则迭代收敛.但这样做计算量太大.,L是Lipschitz常数.,复旦大学力学与工程科学系,9-17,改进Euler法,改进Euler法的思想,称为预测-校正系统,或称改进Euler法.它属于显式单步法.,先用Euler法求yn+1的近似值;再用梯形公式迭代一次,得到yn+1.,编程采用的公式,复旦大学力学与工程科学系,9-18,改进Euler法的算法描述,输入a,b,f(x,y),区间等分数N,初值y0;,fori=1toNstep1计算各点的函数值,停机.,计算步长h=(b-a)/N;置n0,xa,yy0;输出(x,y);,3.1计算yp=y+hf(x,y);3.2置xx+h;计算yq=y+hf(x,yp);3.3置y0.5(yp+yq);3.4输出(x,y);,Endfori,改进Euler法的误差,可以证明,改进Euler法的局部截断误差为O(h3),故属于二阶方法.,复旦大学力学与工程科学系,9-19,改进Euler法的算例,从结果可以看出:精确度比Euler法明显提高.,取步长h=0.1;结果如下:,复旦大学力学与工程科学系,9-20,三.龙格-库塔(Runge-Kutta)法,Runge-Kutta法的基本思想,如何提高精度,局部截断误差由Talor展开式的余项决定,一般地,具有p阶精度,局部截断误差为:,其中,特例.p=1时,即为Euler方法.,复旦大学力学与工程科学系,9-21,障碍及对策,如,Euler法可表示为:,改进Euler法可表示为:,障碍:必须求f(x,y)的各阶偏导数,计算复杂且工作量大.,对策:用f(x,y)在某些点上函数值的线性组合,来计算近似值yn+1,并使它的Talor展开式与y(xn+1)的展开式前面若干项完全相同,从而达到一定阶数的精度这就是Runge-Kutta法(简称RK方法)的基本思想.,只计算一点的函数值,得到的是一阶方法.,须计算两点处的函数值,其Talor展开的前三项与y(xn+1)的展开式相同,得到的是二阶方法.,复旦大学力学与工程科学系,9-22,RK方法的构造,构造方法,二阶RK方法,取p=2,近似公式为:,在(xn,yn)处Talor展开:,整理得:,复旦大学力学与工程科学系,9-23,二阶RK方法（续）,理论上可以证明,无论怎样选取参数c1,c2,a2,b21,上面的公式不可能有更高的精度.即,每步计算两个函数值,只能得出二阶公式.,在xn处Talor展开,共4个未知量,有无穷多个解.代入RK近似公式,其局部截断误差都是O(h3).统称为二阶方法.,复旦大学力学与工程科学系,9-24,几个常用的二阶RK公式,改进Euler公式,取,中点公式,取,Heun二阶公式,取,复旦大学力学与工程科学系,9-25,几个常用的三阶RK公式,Kutta三阶公式,RK一般公式中,取p=3可得三阶RK公式.其中包含c1,c2,c3,a2,a3,b21,b31,b32共6个系数,也有无穷多个解.常用的三阶公式有:,Heun三阶公式,复旦大学力学与工程科学系,9-26,几个常用的四阶RK公式,经典Kutta四阶公式,RK一般公式中,取p=4可得四阶RK公式.常用的四阶公式有:,它是实际应用中最常用的单步法.,复旦大学力学与工程科学系,9-27,几个常用的四阶RK公式（续）,Gill公式,RK方法中p值与阶数的关系,p=1,2,3,4时,得到RK公式的最高阶数分别为一、二、三、四阶.但最高阶数并不一定等于p,当p=5,6,7,8,9时,RK公式的最高阶数分别为:4,5,6,6,7.,复旦大学力学与工程科学系,9-28,经典四阶Kutta方法的算法描述,输入a,b,f(x,y),区间等分数N,初值y0;,fori=1toNstep1计算各点的函数值,停机.,计算步长h=(b-a)/N;置xa,yy0;输出(x,y);,3.1计算K1=f(x,y);3.2置xx+0.5*h;计算:K2=f(x,y+0.5*h*K1);K3=f(x,y+0.5*h*K2);3.3置xa+i*h;计算:K4=f(x,y+h*K3);3.4置yy+h*K1+2(K2+K3)+K4/6;3.5输出(x,y);,Endfori,fori=1toNstep1计算各点的函数值,复旦大学力学与工程科学系,9-29,经典四阶RK公式算例,从结果可以看出:精确度比改进Euler法明显提高.,取步长h=0.1;结果如下:,复旦大学力学与工程科学系,9-30,变步长的RK方法,为何要变步长,由于y(x)是不均匀的,采用等不长计算时,有些点上精度较高,而有些点则精度较低.是否能根据精度要求,自动调整当前计算步的步长呢?,局部截断误差O(hp+1),即,c与有关,假设内变化不大,则可看作常数.,.(1),截断误差O(h/2)p+1,即,.(2),从误差的分析说起,如果分两步推进,复旦大学力学与工程科学系,9-31,变步长控制方法-1,具体实施步骤,可根据的大小来选择合适的步长.,以步长h,从xn点出发计算y(xn+h)的近似值;,.(3),以步长h/2,从xn点出发,用两步计算y(xn+h)的近似值;,如果(3)式的满足要求的精度,即,且两者相差不多,则仍以h为步长继续计算;,如果,说明步长太小.反复以2h为步长计算,直到为止.这时,前一次所用的步长h即为合适的步长,继续推进;,如果,说明步长太大.反复以h/2为步长计算,直到为止.此时所用的步长h即为合适的步长,继续推进;,复旦大学力学与工程科学系,9-32,变步长控制方法的改进,具体实施步骤的改进,即步长减半后,误差降为原来的约1/2p倍.由此可推得,若步长为,则误差降为原来的约倍.,13步与前相同;,如果,说明步长太小.将步长放大至2kh,则误差放大至.令,取满足条件的最大整数k,下一步的步长取为2kh,继续推进;,如果,说明步长太大.将步长缩小至h/2k,则误差缩小至.令,取满足条件的最小整数k,下一步的步长取为h/2k,继续推进;,复旦大学力学与工程科学系,9-33,四.线性多步法,多步法的概念,单步法:计算yn+1时,只用到前一步的信息yn;,RK方法也能提高精度,但需要计算多个点的函数值,计算量较大.,多步法:计算yn+1时,利用前几步的已有信息yn,yn-1,yn-2,从而提高计算精度多步法的基本思想.,最常用的多步法线性多步法的一般形式,其中,如果,需用到fn+1(未知),则公式为隐式的;若,则公式为显式的.,复旦大学力学与工程科学系,9-34,线性多步法的导出,基本方法,将多步法公式在xn处Talor展开,与y(xn+1)在xn处的Talor展开式相比,要求前面的若干项相同,从而确定系数i,i.,以线性两步法为例,r=1,假设f(x)充分光滑.,.(*),线性两步法的推导过程,线性两步法公式:,Talor展开:,其中,假设,都是精确的.由(*)式可得,复旦大学力学与工程科学系,9-35,线性两步法的推导过程(续-1),代入两步法公式:,复旦大学力学与工程科学系,9-36,线性两步法的推导过程(续-2),整理得,y(xn+1)的Talor展开式:,令(1)与(2)式各项相同,得:,.(2),.(1),选取使其满足(3)式前p+1个(p=1,2,3,4)方程,则可得p阶公式.,复旦大学力学与工程科学系,9-37,线性两步法的推导例子,(3)式满足前3式.,二阶公式（p=2）,共5个未知量,有无穷个解.取,满足上面三个方程.二阶公式:,四阶公式（p=4）,(3)式满足前5式.,得唯一解:,得四阶公式,复旦大学力学与工程科学系,9-38,一般形式线性多步法的系数确定,令(4)与(5)式前p+1项相等,可得,有r+1个i,r+2个i,共2r+3个未知量.,.(4),.(5),yn-i及fn-i在xn处作Talor展开,得,复旦大学力学与工程科学系,9-39,一般形式线性多步法的系数确定(续),h0项(即yn项)的系数相等,hk项(即项)的系数相等,局部截断误差,共2r+3个未知量,p+1个方程.要求:p+12r+3,由于p2r+2,即(4)式的线性多步法最多可达到2r+2阶精度.,规定:,复旦大学力学与工程科学系,9-40,常用的线性多步法,阿达姆斯(Adams)公式,线性多步公式为:,取r=3,p=42r+2;,即,局部截断误差为:,复旦大学力学与工程科学系,9-41,阿达姆斯(Adams)公式（续）,仍取r=3,p=4;并取,即,线性多步公式为:,局部截断误差为:,复旦大学力学与工程科学系,9-42,米尔恩(Milne)公式,取r=3,p=4;并取,即,线性多步公式为:,局部截断误差为:,复旦大学力学与工程科学系,9-43,海明(Hamming)公式,取r=2,p=40,使得:,整体截断误差为,设显式单步法具有p阶精度,其增量函数关于y满足Lipschitz条件,初值y0=y(x0)是精确的.则显式单步法的整体截断误差为:,利用Lipschitz条件:,L为常数.,(2),复旦大学力学与工程科学系,9-61,定理9.1证明(续),同理,因为在a,b中求解,故有:,利用不等式:,证毕#,复旦大学力学与工程科学系,9-62,推论1,连续,且关于y满足Lipschitz条件,则显式单步法是收敛的.,设显式单步法具有p阶精度,其增量函数在区域,证明:由定理9.1的结论可直接推得.,由此可得,当f(x,y)在区域D:axb,-0.故三阶RK方法的绝对稳定区间为:,可得的实根,另两个是复根.,复旦大学力学与工程科学系,9-69,四阶RK方法的稳定性,同理可推得四阶RK方法的绝对稳定区域:,四阶RK方法的绝对稳定区间为:,方程的零点:另两个是复根.,方程的零点:四个是复根,可见式恒成立.,当时,式成立.可见,复旦大学力学与工程科学系,9-70,RK方法的稳定区域图,复旦大学力学与工程科学系,9-71,一般形式方程的稳定性考虑,