常微分方程数值解法 6VIP

第四章,常微分方程数值解法,本章研究常微分方程初值问题(1)的数值解法.这里假定f(x,y)满足解的存在唯一性定理及相当光滑等条件，从而(1)有唯一解y(x).,本章的数值解法，它不是求方程的解y(x)的解析表达式或近似表达式，而是直接求一系列离散点上的解值y(xi)的近似值yi.利用计算机解微分方程主要使用数值方法.,取一系列点x0,x1,xn,y(x0)=y0,y(x1)y1,y(xn)yn,y0,y1,yn,称为数值解h=xnxn-1称为步长,数值解法,利用离散化方法将问题(1)转化为关于离散量的相应问题，相应问题的解yn作为y(xn)的近似结果.常用离散化方法：数值微分，数值积分，Taylor展开等.问题(1)可以有很多种数值解法，可以建立许多数值公式.,Euler折线法与改进的Euler法,问题(1)等价于若已知解y(x)在点xn上的函数值y(xn),则有用不同的近似公式计算此式中的定积分值,就得到解初值问题的不同数值解法.,Euler法,用矩形求积公式计算有y(x1)y(x0)+hf(x0,y(x0)=y0+hf(x0,y0)=y1y(x2)y(x1)+hf(x1,y(x1)y1+hf(x1,y1)=y2y(xn)y(xn-1)+hf(xn-1,y(xn-1)yn-1+hf(xn-1,yn-1)=yny(xn+1)y(xn)+hf(xn,y(xn)yn+hf(xn,yn)=yn+1yn+1=yn+hf(xn,yn)（n=0,1,2,）(2)称此方法为Euler折线法或矩形法,改进的Euler法（梯形法）,用梯形求积公式计算积分得(3)这个方法称为改进的Euler法或梯形法运用它常采用下面的迭代格式,下面讨论yn+1(k)的收敛性.由于f(x,y)满足解的存在唯一性定理的条件，从而只要hL/21,yn+1(k)就收敛，设，当k时，等式两边取极限有因此只要h足够小，就能保证yn+1(k)收敛于(3)中的yn+1,Euler预报校正公式,若改进的Euler方法只迭代一次，便得Euler预报-校正公式,Euler预报校正公式,Euler预报-校正公式也可写成形式,局部截断误差设yn=y(xn)，则称y(xn+1)-yn+1为方法的局部截断误差.方法的阶数若数值方法的局部截断误差为O(hp+1)，则称这种方法为p阶方法，这里p为非负整数.,由泰勒展开式知设yn=y(xn)设对于Euler法yn+1=yn+hf(xn,yn)=y(xn)+hf(xn,y(xn)=y(xn)+hy(xn)y(xn+1)-yn+1=O(h2)当h0时，y(xn+1)-yn+1是h2的同阶无穷小,Euler公式yn+1=yn+hf(xn,yn)（n=0,1,2,）它是显格式，该方法的局部截断误差为O(h2)，是一阶方法.,对于Euler预报-校正公式，利用泰勒展开式k1=hf(xn,yn)=hf(xn,y(xn)=hy(xn),y(xn+1)-yn+1=O(h3)当h0时，y(xn+1)-yn+1是h3的同阶无穷小.,Euler预报校正公式,局部截断误差为O(h3)，是二阶方法.,单步法：若求yn+1，只需利用它前一步的信息yn，则称这种方法为单步法。它由y0出发，可求得y1，y2，y3多步法：若求yn+1，需利用它前面至少两个点的息，则称这种方法为多步法.Euler折线法，Euler预报校正公式，Runge-Kutta方法是单步法;Adams外推法，Adams内插法是多步法.,Runge-Kutta方法,Taylor展开式理论上讲，只要解y(x)充分光滑，通过保留Taylor展开式的若干项就可得到任意阶的近似公式，但计算y(x)的各阶导数很麻烦。可间接利用这种思想,Euler法也可写成形式其局部截断误差为O(h2)，是一阶方法.每步计算f的值一次,Euler预报-校正公式也可写成形式局部截断误差为O(h3)，是二阶方法每步计算f的值二次可以通过增加计算f的值的次数，提高公式的阶数（精度）,以f在不同点上的函数值的线性组合来代替yn+1yn，其中有一些可待定选取的待定参数，通过Taylor展开式确定这些待定参数使建立的数值方法按要求达到一定的阶数，这种思想就是Runge-Kutta方法的思想,二阶Runge-Kutta公式,一般形式其中R1,R2,a,b为待定常数其局部截断误差为O(h3)，是二阶方法每步计算f的值二次,设yn=y(xn)k1=hf(xn,yn)=hf(xn,y(xn)把k2中f在(xn,y(xn)处泰勒展开,再将k1,k2代入yn+1中，将其与y(xn+1)泰勒展开式比较，要使y(xn+1)-yn+1=O(h3)，含h0,h1,h2的项相同即有,个未知数，个方程满足条件的解不止一组取就是Euler预报-校正公式,三阶Runge-Kutta公式,一般形式其中R1,R2,R3,a2,a3,b21,b31,b32为待定常数其局部截断误差为O(h4)，是三阶方法每步计算f的值三次,四阶Runge-Kutta公式,一般形式其中有13个待定常数局部截断误差为O(h)，是四阶方法每步计算f的值四次,设yn=y(xn)把k2,k3,k4中f在(xn,y(xn)处泰勒展开后，将k1,k2,k3,k4代入yn+1中，再将yn+1按h的幂重新整理后与y(xn+1)泰勒展开式比较，要使y(xn+1)-yn+1=O(h5)，含h0,h1,h2,h3,h4的项相同从而确定各参数,13个未知数，11个方程满足条件的解不止一组最常用的是,标准四阶Runge-Kutta公式,（n=0,1,2,）局部截断误差为O(h5)，是四阶方法.,例,用Euler方法求解问题取h=0.1解设f(x,y)=-y+x+1，x0=0,y0=1,xn=x0+nh=0.1n（n=0,1,10）Euler格式为由y0=1出发，按上面公式的计算结果并与精确解y(x)进行比较，如表所示,xnyny(xn+1)|y(xn+1)yn+1|01.0000001.00000000.11.0000001.0048370.0048370.21.0100001.0187310.0087310.31.0290001.0408180.0118180.41.0561001.0703200.0142200.51.0904901.1065310.0160410.61.1314411.1488120.0173710.71.1782971.1965850.0182880.81.2304671.2493290.0188620.91.2874201.3065700.0191501.01.3486781.3678790.019201,利用Euler方法计算积分在点x=0.5,1,1.5,2处的数值解.,解令则有等价的问题对此一阶常微分方程初值问题，取步长h=0.5，设,由Euler格式,从y0=0出发计算的数值解如表xnyn0.50.5000001.01.1420131.52.5011542.07.245022,用改进的Euler法（梯形公式）解初值问题,取步长h=0.2，小数点值至少保留5位.解设f(x,y)=8-3y，x0=1,y0=2,xn=x0+nh=1+0.2n（n=0,1,5）梯形公式为于是,整理得显格式由y0=2出发，计算结果如表所示xnynxnyn121.62.562591.22.307691.82.610621.42.473372.02.63649,用Euler预报校正格式求解初值问题,要求取步长h=0.2，计算y(1.2)及y(1.4)的近似值，小数点后至少保留5位.解设f(x,y)=-y-y2sinx,x0=1,y0=1,xn=x0+nh=1+0.2n,Euler预报校正格式为,于是有由y0=1计算得,取步长h=0.4，写出用标准四阶Runge-Kutta方法求解初值问题,的计算公式，并计算y(1.8)的近似值，小数点后至少保留6位.解设f(x,y)=xsin(x+y),x0=1,y0=0,xn=x0+nh=1+0.4n，（n=0,1,20）标准四阶Runge-Kutta公式为,代入f(x,y)=xsin(x+y)有,由y0=0计算得y(1.4)y1=0.460389y(1.8)y2=0.911704,例,初值问题有精确解若xn=nh，yn是用Euler法得到的y(x)在x=xn处的近似值.证明：,Euler格式为yn+1=yn+hf(xn,yn)（n=0,1,2,）代入f(x,y)=ax+b，则yn+1=yn+h(axn+b)由y0=0得,因xn=nh，于是所以,求系数a，b，c和d使数值计算公式有y(xn+1)yn+1=O(h5),解设则都在xn处作Taylor展开到y(4)(x)，即,再按h作升幂排列，即而,要使y(xn+1)yn+1=O(h5)，比较h的同次幂系数，必有,Adams显式（外推）公式,取4个等距步长h的插值点xn-3,xn-2,xn-1,xn,用f(x,y(x)在这四个节点上的函数值作f(x,y(x)的三次插值多项式3(x),则用yn-i代替y(xn-i),fn-i=f(xn-i,yn-i),得四阶Adams显式(外推)公式其局部截断误差为,它是一个四阶方法.,Ada