东北大学数值分析第八章 常微分方程数值解法

第8章常微分方程数值解法§1.1为什么要研究数值解法一阶常微分方程初值问题的一般形式为y=(x,y),axb§1引言(8.1)y(a)=其中(x,y)是已知函数,为给定的初值.如果函数(x,y)在区域axb,-y上连续且关于y满足Lipschitz条件其中L0为Lipschitz常数,则初值问题(8.1)有唯一解.yxyyLyxfyxf,,),(),(所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值.a=x0x1x2…xn…xN=b其中剖分节点xn=a+nh,n=0,1,…,N,h称为剖分步长.数值解法就是求精确解y(x)在剖分节点xn上的近似值yny(xn),n=1,2,…,N.假设初值问题(8.1)的解y=y(x)唯一存在且足够光滑.对求解区域[a,b]做剖分我们采用数值积分方法来建立差分公式.§1.2构造数值解法的基本思想在区间[xn,xn+1]上对方程(8.1)做积分,则有对右边的积分应用左矩形公式，则有)2.8())(,()()(11nnxxnndxxyxfxyxy梯形公式oxyab左矩形公式y=(x)babfafabdxxf)]()([2)(baafabdxxf)()()(右矩形公式babfabdxxf)()()(中矩形公式babafabdxxf)2()()(对右边的积分应用左矩形公式，则有)2.8())(,()()(11nnxxnndxxyxfxyxy因此，建立节点处近似值yn满足的差分公式称之为Euler公式.称为梯形公式.))(,()()(1nnnnxyxhfxyxy),(1nnnnyxhfyy1,2,1,0,0Nny若对(8.2)式右边的积分应用梯形求积公式，则可导出差分公式1,2,1,0,0Nny)],(),([2111nnnnnnyxfyxfhyy利用Euler方法求初值问题解此时的Euler公式为称为Euler中点公式或称双步Euler公式.若在区间[xn-1,xn+1]上对方程(8.1)做积分,则有11))(,()()(11nnxxnndxxyxfxyxy对右边的积分应用中矩形求积公式，则得差分公式),(211nnnnyxhfyy1,2,1,0,0Nny例120,21122xyxy0)0(y的数值解.此问题的精确解是y(x)=x/(1+x2).分别取步长h=0.2,0.1,0.05,计算结果如下)211(221nnnnyxhyy2,1,0,00nyhxnyny(xn)y(xn)-ynh=0.20.000.400.801.201.602.000.000000.376310.542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.201.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.491800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.00227Euler中点公式则不然,计算yn+1时需用到前两步的值yn,yn-1,称其为两步方法,两步以上的方法统称为多步法.在Euler公式和梯形公式中,为求得yn+1,只需用到前一步的值yn,这种差分方法称为单步法,这是一种自开始方法.隐式公式中,每次计算yn+1都需解方程,要比显式公式需要更多的计算量,但其计算稳定性较好.在Euler公式和Euler中点公式中,需要计算的yn+1已被显式表示出来,称这类差分公式为显式公式,而梯形公式中,需要计算的yn+1隐含在等式两侧,称其为隐式公式.从数值积分的角度来看,梯形公式计算数值解的精度要比Euler公式好,但它属于隐式公式,不便于计算.实际上,常将Euler公式与梯形公式结合使用:§2改进的Euler方法和Taylor展开方法§2.1改进的Euler方法)],(),([2111nnnnnnyxfyxfhyy1,2,1,0,0Nny),(]0[1nnnnyxhfyy)],(),([2][11]1[1knnnnnknyxfyxfhyy1,2,1,0,0Nny由迭代法收敛的角度看，当(是给定的精度要求)时,取就可以保证迭代公式收敛,而当h很小时,收敛是很快的.而且,只要||][1]1[1knknyy.]1[11knnyy,12Lyfh),(1nnnnyxhfyy)],(),([2111nnnnnnyxfyxfhyy1,2,1,0,0Nny通常采用只迭代一次的算法:称之为改进的Euler方法.这是一种单步显式方法.改进的Euler方法也可以写成)(2211KKhyynn),(1nnyxfK1,2,1,0,0Nnyy=y-2x/y,0x1的数值解,取步长h=0.1.[精确解为y(x)=(1+2x)1/2.]),(12hKyhxfKnn例2求初值问题y(0)=1解(1)利用Euler方法nnnnyxyy/2.01.119,2,1,0,10ny)(05.0211KKyynnnnnyxyK/219,2,1,0,10ny计算结果如下:1121.0)1.0(21.0KyxKyKnnn(2)利用改进Euler方法nxnEuler方法yn改进Euler法yn精确解y(xn)01234567891000.10.20.30.40.50.60.70.80.9111.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.7177791.78477011.0959091.1840961.2662011.3433601.4164021.4859561.5525151.6164761.6781681.73786911.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.732051在节点xn+1的误差y(xn+1)-yn+1,不仅与yn+1这一步计算有关,而且与前n步计算值yn,yn-1,…,y1都有关.为了简化误差的分析,着重研究进行一步计算时产生的误差.即假设yn=y(xn),求误差y(xn+1)-yn+1,这时的误差称为局部截断误差,它可以反映出差分公式的精度.§2.2差分公式的误差分析如果单步差分公式的局部截断误差为O(hp+1),则称该公式为p阶方法.这里p为非负整数.显然,阶数越高,方法的精度越高.研究差分公式阶的重要手段是Taylor展开式,一元函数和二元函数的Taylor展开式为:另外,在yn=y(xn)的条件下,考虑到y(x)=(x,y(x)),则有321!3)(!2)()()()()(hxyhxyhxyxyhxyxynnnnnn2222222),(),(2),(!21),(),(),(),(kyyxfhkyxyxfhxyxfkyyxfhxyxfyxfkyhxfnnnnnnnnnnnnnny(xn)=(xn,y(xn))=(xn,yn)=ny(xn)=nnnnnfyfxfxyxfdxd))](,([nnnnnnnnnfyfyfxffyffyxfxf2222222)(2][)(nnnnfyfxfdxdxyyn+1=yn+h(xn,yn)对Euler方法，有21!2)()()()()(hxyhxyxyhxyxynnnnn=yn+(xn,yn)h+O(h2)从而有:y(xn+1)-yn+1=O(h2)所以Euler方法是一阶方法.再看改进Euler方法,因为),(12hKyhxfKnn1hKyfhxffnnn)(221321222122222hOKhyfKhyxfhxfnnn可得所以,改进的Euler方法是二阶方法.而nnnnnnfyxfhhfyy221)(244222223hOfyffyxfxfhnnnnn)(!3)(2)()()()(4321hOhxyhxyhxyxyxynnnnnnnnnnfyfxfhhfy22)()(26422222223hOfyfyfxffyffyxfxfhnnnnnnnnn从而有:y(xn+1)-yn+1=O(h3)§2.3Taylor展开方法设y(x)是初值问题(8.1)的精确解,利用Taylor展开式可得称之为p阶Taylor展开方法.………………1)1()(21)!1()(!)(!2)()()()(pppnpnnnnhpyhPxyhxyhxyxyxy)())(,(!))(,(!2))(,()(1)1()1(2pnnppnnnnnhOxyxfPhxyxfhxyxhfxy因此，可建立节点处近似值yn满足的差分公式),(!),(!2),()1()1(21nnppnnnnnnyxfPhyxfhyxhfyy1,2,1,0,0Nny),(),(),(),()1(yxfyyxfxyxfyxffyfyfxffyffyxfxfyxf2222222)2()(2),(其中所以,此差分公式是p阶方法.由于Taylor展开方法涉及很多复合函数(x,y(x))的导数的计算,比较繁琐,因而很少直接使用,经常用它为多步方法提供初始值.然而,Taylor展开方法给出了一种构造单步显式高阶方法的途径.Euler方法可写为可见,公式的局部截断误差为:y(xn+1)-yn+1=O(hp+1).§3Runge-Kutta方法§3.1Runge-Kutta方法的构造hKyynn1),(nnyxfK构造差分公式…………改进的Euler方法可写为)(2211KKhyynn),(1nnyxfK),(12hKyhxfKnn)(22111ppnnKKKhyy),(1nnyxfK),(12122KhyhxfKnn),(11piipinPnPKhyhxfK其中i,i,ij为待定参数.若此公式的局部截断误差为由于yn+1=yn+h1n+h2(n+hxn+hnyn)+O(h3)O(hp+1),称公式为p阶Runge-kutta方法,简称p阶R-K方法.对于p=2的情形,应有)(22111KKhyynn),(1n