数值分析第九章常微分方程数值解法

1,第九章常微分方程数值解法,1Euler方法2改进的Euler方法3Runge-Kutta法4线性多步法,2,近似求解,3,一阶常微分方程初值问题,为保证上述问题解的存在唯一性,对函数f做如下假设:f关于y的偏导数有界,即存在常数L0,使得,本章中的L均指此数.,常微分方程解的存在唯一性,4,数值解法:就是寻求解y(x)在一系列离散点,上的近似值,相邻两个节点的间距称为步长.本章均假设,此时节点,常微分方程数值解法,5,初值问题的各种差分方法都采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进.,只要给出从已知信息计算的递推公式，此类计算格式统称为差分格式.,常微分方程差分格式,6,数值求解一阶常微分方程初值问题,难点:如何离散y?,常见离散方法,差商近似导数数值积分方法Taylor展开方法,7,结果有,设用y(xn)的近似值yn代入上式右端，记所求结果为yn+1，这样导出的计算公式,将方程y=f(x,y)中的y在xn点用向前差商近似,差商近似导数(Euler公式),Euler公式,单步法,显式公式,8,Euler公式计算总结,将区间a,bN等分,分点xn=x0+nh(n=1,2,N)步长,用Euler公式,求得yn(n=1,2,N.),9,结果有,由此得差分格式计算公式,将方程y=f(x,y)中的y在xn+1点用向后差商近似,向后Euler公式,单步法,隐式公式,差商近似导数(向后Euler公式),10,设将方程y=f(x,y)的两端从xn到xn+1求积分,得,用不同的数值积分方法近似上式右端积分,可以得到计算y(xn+1)的不同的差分格式.,用左矩形公式计算积分项,再离散化，得Euler公式,数值积分方法(Euler公式),11,用右矩形公式计算积分项,再离散化，得向后Euler公式,数值积分方法(向后Euler公式),12,用梯形法计算积分项,再离散化，即可得如下计算公式,梯形公式,单步法,隐式公式,数值积分方法(梯形公式),13,Taylor展开方法(Euler公式),由此亦可得Euler公式,将y(xn+1)在xn点作Taylor展开,得,14,1Euler方法,用差分格式,来近似微分方程初值问题,15,例用Euler公式求解初值问题,解,取步长h=0.1,计算结果见下表,计算公式为,16,y(xn)1.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321,xn0.10.20.30.40.50.60.70.80.91.0,yn1.10001.19181.27741.35821.43511.50901.58031.64981.71781.7848,两者比较可以看出Euler法的精度较差.,17,h,h,h,Straightlineapproximation,Euler法的几何意义,x0,x1,x2,x3,y0,18,向后Euler法,这是一个隐式公式.已知yn,上式是关于yn+1的非线性方程,不能从中直接求出yn+1,需要用非线性方程求根的方法来求解,通常用迭代法.,19,已知yn求yn+1的迭代法,则当hL1时,迭代函数,向后Euler法计算量大,但数值稳定性好!,20,Euler法的局部截断误差,局部截断误差:用Euler法计算一步所产生的误差,即假设第n步的计算是精确的yn=y(xn),y(xn+1)与yn+1的差,即,这种误差称为局部截断误差,简称截断误差.,21,当yn=y(xn)时,局部截断误差主项,Euler法的局部截断误差,22,向后Euler法的局部截断误差,定义其局部截断误差为,向后Euler法的计算公式,23,向后Euler法的局部截断误差,局部截断误差主项,24,在不考虑舍入误差的情况下,称y(xn+1)与yn+1之差为整体截断误差,记为,令,则,Euler法的整体截断误差,25,Euler法的整体截断误差,26,注意e0=0,故,27,所以,这表明,Euler方法的整体截断误差与h同阶.当h0时,eN0.,如果某种数值方法的局部截断误差为O(hp+1),则称它的精度是p阶的,或称之为p阶方法.由此可知Euler方法仅为一阶方法.,28,梯形公式,2改进的Euler法,从直观上看,用梯形公式计算数值积分要比矩形公式精度高.相应地求解常微分方程的梯形公式的精度也要比Euler公式精度高.,29,梯形公式的局部截断误差,局部截断误差主项,故梯形公式是二阶方法,30,则当hL2时,梯形公式是一个隐式公式.已知yn,用迭代法计算yn+1,用梯形法计算,从yn到yn+1要进行迭代计算,计算量较大.实际计算常把Euler法和梯形法结合起来,这就是改进的Euler法.,31,先用Euler法求得一个初步的近似值，记为，称之为预测值，然后用它替代梯形法右端的yn+1再直接计算fn+1,得到校正值yn+1,这样建立的预测校正系统称为改进的Euler法：,预测,校正,改进的Euler法,32,实践表明,改进Euler法的精度明显高于Euler法.利用Taylor展开可以证明,改进Euler的局部截断误差为O(h3),故它是二阶方法.,它有下列平均化形式：,33,例用改进Euler法求解,解,取步长h=0.1,计算结果见下表,改进Euler法公式为,34,y(xn)1.09541.18321.26491.34161.41421.48321.54921.61251.67331.7321,xn0.10.20.30.40.50.60.70.80.91.0,yn1.09591.18411.26621.34341.41641.48601.55251.61651.67821.7379,同Euler法的计算结果比较,改进Euler法明显改善了精度.,35,改进Euler法的几何意义,(K1+K2)/2,K2,K1,xn,xn+1,36,考察差商根据微分中值定理，存在点，利用所给方程y=f(x,y)得,称K*=f(,y()为区间xn,xn+1上的平均斜率,这样只要对平均斜率K*提供一种算法，相应地便导出一种计算格式.,3Runge-Kutta法,Runge-Kutta法的基本思想,37,梯形法：,改进Euler法：,其中,Euler法：,向后Euler法：,38,RungeKutta法的基本思想就是设法在xn,xn+1内多预报几个点的斜率值，然后把它们加权平均作为平均斜率，以期望构造出更高精度的计算格式.,Runge-Kutta法的基本思想,39,考察区间xn,xn+1内一点xn+p=xn+ph,0n)上产生的偏差值均不超过,则称该方法是绝对稳定的.,用某种数值方法求解上述模型方程，对固定的h,使得该方法绝对稳定的h和的全体称为该方法的绝对稳定区域.,88,Euler方法的绝对稳定区域,模型问题,计算格式,设实际计算yn时产生误差n,它使节点值yn+1产生误差n+1,则,从而,为使误差不扩大,需,89,向后Euler方法的绝对稳定区域,模型问题,计算格式,设实际计算yn时产生误差n,它使节点值yn+1产生误差n+1,则,向后Euler方法的绝对稳定区域,90,梯形公式的绝对稳定区域,模型问题,计算格式,91,类似于Euler方法的讨论,得梯形法的绝对稳定区域为,即,梯形公式的绝对稳定区域,92,RK方法的绝对稳定区域,二阶RK方法的绝对稳定区域,对于模型问题,改进Euler法,93,当是实数时,2.51h0,当是实数时,2.78h0,二阶RK方法的绝对稳定区域,当是实数时,2h0,三阶和四阶RK方法的绝对稳定区域分别为,94,例y=20y(0x1),y(0)=1,分别取h=0.1及h=0.2用经典的四阶RK方法计算.,h=0.24.9825.0125.0625.03125.0,xn0.20.40.60.81.0,h=0.10.931010.121010.141020.151030.17104,解,用四解RK方法计算,其误差见下表,本例中=20,h分别为2及4,前者在绝对稳定区间内,后者则不在.从以上结果看到,如果步长h不满足绝对稳定条件,误差增长很快.,95,前面研究了单个方程y=f的数值解法,只要把y和f理解为向量,所提供的各种算法即可推广应用到一阶方程组的情形.,6一阶微分方程组与高阶微分方程的数值解法,一阶方程组的数值解法,96,以yn,zn表示节点xn上y(x),z(x)的近似解，则其改进的Euler法具有形式,例对于方程组,令,原方程组变为,97,预报,即,98,校正,即,99,关于高阶微分方程（或方程组）的初值问题，原则上总可以归结为一阶微分方程组来求解.,引进新变量z=y，即可化为一阶方程组的初值问题,从而可以运用前面介绍的算法来求解.,高阶微分方程的数值解法,例对于下列二阶方程的初值问题,100,三阶Adams显式公