第八章--常微分方程数值解法课件

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,第8章 常微分方程数值解法,1,.1 为什么要研究数值解法,一阶常微分方程初值问题的一般形式为,y,=(x,y),axb,1 引言,(,8.1),y,(a)=,其中(,x,y),是已知函数,为给定的初值.,如果函数(,x,y),在区域,axb,-y0,为,Lipschitz,常数,则初值问题(8.1)有唯一解.,所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值.,a=x,0,x,1,x,2,x,n,x,N,=b,其中剖分节点,x,n,=a+,nh,n=0,1,N,h,称为剖分步长.数值解法就是求精确解,y(x),在剖分节点,x,n,上的近似值,y,n,y(,x,n,),n=1,2,N.,假设初值问题(8.1)的解,y=y(x),唯一存在且足够光滑.对求解区域,a,b,做剖分,我们采用数值积分方法来建立差分公式.,1,.2 构造数值解法的基本思想,在,区间,x,n,x,n,+1,上对方程(8.1)做积分,则有,对,右边的积分应用左矩形公式，则有,梯形公式,o,x,y,a,b,左,矩形公式,y=,(x),右矩形公式,中矩形公式,对,右边的积分应用左矩形公式，则有,因此，建立节点处近似值,y,n,满足的差分公式,称之为,Euler,公式,.,称为,梯形公式,.,若对(8.2)式右边的积分应用梯形求积公式，则可导出差分公式,利用,Euler,方法求初值问题,解,此时的,Euler,公式为,称为,Euler,中点公式,或称,双步,Euler,公式,.,若在区间,x,n,-1,x,n,+1,上对方程(8.1)做积分,则有,对,右边的积分应用中矩形求积公式，则得差分公式,例1,的,数值解.此问题的精确解是,y(x)=x/(1+x,2,).,分别取步长,h=0.2,0.1,0.05,计算结果如下,h,x,n,y,n,y(,x,n,),y(,x,n,)-,y,n,h=0.2,0.00,0.40,0.80,1.20,1.60,2.00,0.00000,0.37631,0.54228,0.52709,0.46632,0.40682,0.00000,0.34483,0.48780,0.49180,0.44944,0.40000,0.00000,-0.03148,-0.05448,-0.03529,-0.01689,-0.00682,h=0.1,0.00,0.40,0.80,1.20,1.60,2.00,0.00000,0.36085,0.51371,0.50961,0.45872,0.40419,0.00000,0.34483,0.48780,0.49180,0.44944,0.40000,0.00000,-0.01603,-0.02590,-0.01781,-0.00928,-0.00419,h=0.05,0.00,0.40,0.80,1.20,1.60,2.00,0.00000,0.35287,0.50049,0.50073,0.45425,0.40227,0.00000,0.34483,0.48780,0.49180,0.44944,0.40000,0.00000,-0.00804,-0.01268,-0.00892,-0.00481,-0.00227,Euler,中点公式则不然,计算,y,n,+1,时需用到前两步的值,y,n,y,n,-1,称其为,两步方法,两步以上的方法统称为,多步法,.,在,Euler,公式和梯形公式中,为求得,y,n,+1,只需用到前一步的值,y,n,这种差分方法称为,单步法,这是一种自开始方法.,隐式公式中,每次计算,y,n,+1,都需解方程,要比显式公式需要更多的计算量,但其计算稳定性较好.,在,Euler,公式和,Euler,中点公式中,需要计算的,y,n,+1,已被显式表示出来,称这类差分公式为,显式公式,而梯形公式中,需要计算的,y,n,+1,隐含在等式两侧,称其为,隐式公式,.,从数值积分的角度来看,梯形公式,计算数值解的精度要比,Euler,公式好,但它属于隐式公式,不便于计算.,实际上,常将,Euler,公式与梯形公式结合使用:,2,改进的,Euler,方法和,Taylor,展开方法,2.1,改进的,Euler,方法,由迭代法收敛的角度看，当,(,是,给定的精度要求)时,取,就可以保证迭代公式收敛,而当,h,很小时,收敛是很快的.,而且,只要,通常采用只迭代一次的算法:,称之为,改进的,Euler,方法,.,这是一种单步显式方法.,改进的,Euler,方法也可以写成,y,=y-2x/y ,0 x1,的数值解,取步长,h=0.1.,精确解为,y(x)=(1+2x),1/2,.,例2,求初值问题,y,(0)=1,解,(1)利用,Euler,方法,计算结果如下:,(2)利用改进,Euler,方法,n,x,n,Euler,方法,y,n,改进,Euler,法,y,n,精确解,y(,x,n,),0,1,2,3,4,5,6,7,8,9,10,0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1,1.1,1.191818,1.277438,1.358213,1.435133,1.508966,1.580338,1.649783,1.717779,1.784770,1,1.095909,1.184096,1.266201,1.343360,1.416402,1.485956,1.552515,1.616476,1.678168,1.737869,1,1.095445,1.183216,1.264991,1.341641,1.414214,1.483240,1.549193,1.612452,1.673320,1.732051,在节点,x,n,+1,的误差,y(,x,n,+1,)-,y,n,+1,不仅与,y,n,+1,这一步计算有关,而且与前,n,步计算值,y,n,y,n,-1,y,1,都有关.,为了简化误差的分析,着重研究进行一步计算时产生的误差.即假设,y,n,=y(,x,n,),求误差,y(,x,n,+1,)-,y,n,+1,这时的误差称为,局部截断误差,它可以反映出差分公式的精度.,2.2,差分公式的误差分析,如果单步差分公式的局部截断误差为,O(h,p+1,),则称该公式为,p,阶方法,.这里,p,为非负整数.显然,阶数越高,方法的精度越高.,研究差分公式阶的重要手段是,Taylor,展开式,一元函数和二元函数的,Taylor,展开式为:,另外,在,y,n,=y(,x,n,),的条件下,考虑到,y(x)=(x,y(x),则有,y(,x,n,)=(,x,n,y(,x,n,)=(,x,n,y,n,)=,n,y(,x,n,)=,y,n,+1,=,y,n,+h(,x,n,y,n,),对,Euler,方法，有,=,y,n,+(,x,n,y,n,)h+O(h,2,),从而有:,y(,x,n,+1,)-,y,n,+1,=O(h,2,),所以,Euler,方法是一阶方法.,再看,改进,Euler,方法,因为,可得,所以,改进的,Euler,方法是二阶方法.,而,从而有:,y(,x,n,+1,)-,y,n,+1,=O(h,3,),2.3,Taylor,展开方法,设,y(x),是,初值问题(8.1)的精确解,利用,Taylor,展开式可得,称之为,p,阶,Taylor,展开方法.,因此，可建立节点处近似值,y,n,满足的差分公式,其中,所以,此差分公式是,p,阶方法.,由于,Taylor,展开方法涉及很多复合函数(,x,y(x),的导数的计算,比较繁琐,因而很少直接使用,经常用它为多步方法提供初始值.然而,Taylor,展开方法给出了一种构造单步显式高阶方法的途径.,Euler,方法可写为,可见,公式的局部截断误差为:,y(,x,n,+1,)-,y,n,+1,=O(h,p+1,).,3,Runge,-,Kutta,方法,3.1,Runge,-,Kutta,方法的构造,构造差分公式,改进的,Euler,方法可写为,其中,i,i,ij,为待定参数.,若,此公式的局部截断误差为,由于,y,n,+1,=,y,n,+h,1,n,+h,2,(,n,+h,xn,+h,n,yn,)+O(h,3,),O(,h,p+1,),称公式为,p,阶,Runge,-,kutta,方法,简称,p,阶,R-K,方法,.,对于,p=2,的情形,应有,=,y,n,+h(,1,+,2,),n,+h,2,2,(,xn,+,n,yn,)+O(h,3,),所以,只要令,1,+,2,=1,2,=1/2,2,=1/2 (8.4),一般地,参数由(8.4)确定的一族差分公式(8.3)统称为二阶,R-K,方法.,称之为,中点公式,或可写为,若取,=1,则得,1,=,2,=1/2,=1,此时公式(8.3)就是改进的,Euler,公式;,若取,1,=0,则得,2,=1,=1/2,公式(8.3)为,高阶,R-K,公式可类似推导.,下面列出常用的三阶、四阶,R-K,公式.,四阶,标准,R-K,公式,三阶,R-K,公式,解,四,阶,标准,R-K,公式为,例3,用四阶标准,R-K,方法求初值问题,y,=y-2x/y ,0 x1,y,(0)=1,的数值解,取步长,h=0.2.,计算结果如下:,n,x,n,y,n,y(,x,n,),n,x,n,y,n,y(,x,n,),0,1,2,0.0,0.2,0.4,1.00,1.1832,1.3417,1.00,1.1832,1.3416,3,4,5,0.6,0.8,1.0,1.4833,1.6125,1.7321,1.4832,1.6125,1.7321,也可以构造隐式,R-K,方法,其一般形式为,称之为,p,级隐式,R-K,方法,同显式,R-K,方法一样确定参数.如,是二级二阶隐式,R-K,方法,也就是梯形公式.但是,p,级隐式,R-K,方法的阶可以大于,p,例如,一级隐式中点公式为,或写为,它是二阶,方法.,3.2,变步长,Runge,-,Kutta,方法,以,p,阶,R-K,方法为例讨论.设从,x,n,以步长,h,计算,y(,x,n,+1,),的近似值为,局部截断误差为,其中,C,是与,h,无关的常数.,如果将步长减半,取,h/2,为步长,从,x,n,经两步计算得到,y(,x,n,+1,),的近似值记为,其,局部截断误差为,于是有,从而,得到事后误差估计,可见,当,成立时,可取,.否则,应将步长再次减半进行计算.,求解初值问题,的,单步显式方法可一统一写为如下形式,y,n,+1,=,y,n,+h(,x,n,y,n,h)(8.5),对于,Euler,方法,有,4,单步方法的收敛性和稳定性,4.1,单步方法的收敛性,y,=(x,y),axb,y,(a)=,其中(,x,y,h),称为,增量函数,.,(,x,y,h)=(x,y),对于改进的,Euler,方法,有,(,x,y,h)=1/2(x,y)+(x+h,y+h(x,y),设,y(x),是初值问题(8.1)的解,y,n,是单步法(8.5)产生的近似解.如果对任意固定的点,x,n,均有,y(,x,n,),则称,单步法(8.5)是收敛的,.,可见,若方法(8.5)是收敛的,则当,h0,时,整体截断误差,e,n,=y(,x,n,)-,y,n,将趋于零.,定理8.1,设单步方法(8.5)是,p1,阶方法,增量函数(,x,y,h),在,区域,axb,-yn),的变化均不超过,则称此差分方法是,绝对稳定,的,.,讨论数值方法的稳定性,通常仅限于典型的试验方程,y=y,其中是复数且,Re()0.,在复平面上,当方法稳定时要求变量,h,的取值范围称为方法的,绝对稳定域,它与实轴的交集称为,绝对稳定区间,.,将,Euler,方法应用于方程,y=y,得到,设在,计算,y,n,时产生误差,n,计算值,y,n,=,y,n,+,n,则,n,将对以后各节点值计算产生影响