山东大学数值分析课第9章 微分方程初值问题的数值解法-1

第九章微分方程初值问题的数值解法内容提纲引言Euler法及其改进Runge-Kutta方法线性多步法误差分析数值解法的收敛性、相容性和稳定性边值问题数值解法简介引言初值问题的数值解法:求初值问题的解在一系列节点的值y(xn)的近似值yn的方法.本章数值解法的特点:都是采用“步进式”,即求解过程顺着节点排列的次序一步步向前推进.基本知识:(1)定理1:如果函数f(x,y)在区域上连续,且关于y满足Lipschitz条件常微分方程初值问题:00(,),[,]()dyfxyxabdxyxy求未知函数y=y(x).{(,)|,}DxyaxbyR121212|(,)(,)|||(,),(,),0fxyfxyLyyxyDxyDL此时Lipschitz条件显然成立.故常用在D上连续有界来代替f(x,y)关于y满足Lipschitz条件.注:如无特别说明,总假设(1)的解存在唯一且足够光滑.在f(x,y)对变量y可微的情形下,若偏导数连续有界,则可取L为除了要保证(1)有唯一解外,还需保证微分方程本身是稳定的,即(1)的解连续依赖于初始值和函数f(x,y).也就是说,当初始值y0及函数f(x,y)有微小变化时,只能引起解的微小变化.(其中L称为Lipschitz常数),则对任何,初值问题(1)在[a,b]上存在唯一连续可微解y=y(x).00(,)xyDfy(,)(,)max||xyDfxyLyfy定理2:如果函数f(x,y)在区域上关于y满足Lipschitz条件,则(1)是稳定的.{(,)|,}DxyaxbyR单步迭代:计算yn+1时仅用yn;初值问题(1)与下列积分方程的解等价:00()(,())xxyxyftytdt初值问题的数值解就是求一系列节点上函数y=y(x)的近似值.称为步长.一般取等步长h.01naxxxb01,,,,nyyy1iiixxx多步迭代:计算yn+1时除用yn外,还要用到yn1,yn2,…;k步迭代要用到yn1,yn2,…,ynk+1.显式单步迭代:1(,,)nnnnyyhxyh隐式单步迭代:11(,,,)nnnnnyyhxyyh(2)一、Euler方法及其改进将[a,b]n等分,记,(0,1,,)kbahxakhknn微分法:111100()()()()()(,)0,1,,1()kkkkkkkkkkkyxyxyxyxyxxxhyyhfxyknyxy积分法:11(1)()()(,())kkxkkxyxyxftytdt积分项利用矩形公式计算11(,())(,)()()(,)kkxkkkkkkxftytdthfxyyxyxhfxy1.显式Euler方法(★)Taylor公式推导:2111()()()(),2(,)0,1,,1kkkkkkkkkkkhyxyxhyxyxxyyhfxyknEuler公式几何意义:P1P2Pk也称折线法P0xy2.梯形法称之为梯形公式.这是一个隐式公式,通常用迭代法求解.具体做法:取先用Euler法求出初值,即,将其代入梯形公式的右端,使之转化为显式公式,即注:当f(x,y)关于y满足Lipschitz条件且步长h满足直至满足:若采用梯形公式计算(★)中的积分项,则有111111()()[(,())(,())]2[(,)(,)]2kkkkkkkkkkkkhyxyxfxyxfxyxhyyfxyfxy(0)1ky(0)1(,)kkkkyyhfxy(1)()111[(,)(,)]2llkkkkkkhyyfxyfxy(1)()11||llkkyy(1)11lkkyy类似地,可得23,,kkyy(☆)112hL时,迭代格式(☆)收敛.3.改进的Euler方法把Euler法作为预报(称为预估公式),把隐式的梯形公式作为校正(称为校正公式),则得改进的Euler方法:1111(,)[(,)(,)]2kkkkkkkkkkyyhfxyhyyfxyfxy或11[(,)(,(,))]2kkkkkkkkhyyfxyfxyhfxy也称为预估-校正法.有时为了方便,预估-校正格式也写成下面形式:11122121[](,)(,)kkkkkkyyKKKhfxyKhfxhyhK二、单步法的局部截断误差及精度Def1:先假设,再估计误差这种误差称为单步迭代法在xk+1处的局部截断误差.()kkyxy,1()[()(,(),)]khkkkkRyxyxhxyxhDef2:若某种数值方法的局部截断误差为,则称该数值方法的精度为P阶的.1()pOh注:通常情况下,P越大,h越小,则截断误差越小,数值方法越精确.设10.Euler方法是一阶方法.()kkyyx所以Euler方法为一阶方法.而211211()()()(),,2()()(,())()kkkkkkkkkkkyxyxhyxhyxxyxyxhfxyxOh122,11(,),()()()2kkkkkhkkkyyhfxyhRyxyyOh20.梯形法是二阶方法.2312311123111()()()()()2()()()()()22[()()][()()][()()]()2kkkkkkkkkkkkkkhyxyxhyxyxOhhyxyxhyxyxOhhyxyxhyxyxyxyxOhTaylor展开将代入上式,得1()()()kkyxyxOh311()()[()()]()2kkkkhyxyxyxyxOh11131131111()()[()()]2[(,)(,)]()2[(,())(,)]()2kkkkkkkkkkkkkkhyxyyxyxyxyhfxyfxyOhhfxyxfxyOh11111(,)11(,())(,)|(())kkkkkxkkffxyxfxyyxyy而代入上式得:13(,)11(1|)(())()2kxkkhfyxyOhy当h充分小时,若,则可选取h,使得||1fLy故梯形法的精度为2.同样可以证明改进的Euler法也是二阶方法.梯形法的局部截断误差为:3,111()(),12khkkkkkkhRyxyfxx1(,)|||12kxhfy从而111(,)(,)33,11(,)11|21|2()(1|)()()2kkkxxkhkkxhfhfyyhfRyxyOhOhy例1:1)0(20,2yxxyydxdy取步长h=2/10,2/20,2/30,2/40,分别用欧拉法、改进的欧拉法和梯形法求解.解:记f(x,y)=y－xy2,xk=kh(k=0,1,2,···,n)(1).Euler法:yk+1=yk+h(yk－xkyk2)(k=0,1,···,n)y0=1当h=2/10时,n=10.由Euler公式可得:k01234yk+11.21.38241.5061.535041.46503k56789yk+11.328771.170771.021130.891690.783788(2).改进的Euler法:2212210[()()2(())],0,1,,1kkkkkkkkkkkkkkhyyyxyyhyxyxyhyxyknyk01234yk+11.19121.343841.423481.419051.3473k56789yk+11.237261.114240.9941510.8847510.788666(3).梯形法(计算过程略)221111[()],0,1,,2kkkkkkkkhyyyxyyxykn(1)2()()21111[()()],0,1,2,2lllkkkkkkkkhyyyxyyxyln10203040h0.20.10.06670.05误差0.10590.05210.03420.0256Euler法误差:改进的Euler法误差:n10203040h0.20.10.06670.05误差0.01230.00260.00115.9612e-004-101234500.511.5预-校方法,h=0.2时误差最大值:0.0123-101234500.511.52欧拉方法,h=0.2时误差最大值:0.1059xexxy211)(解析解:三、Runge-Kutta方法1、Taylor级数法设初值问题有解y(x),由Tayler公式得:00(,),(),yfxyyxyaxb2()11()()()()()()2!!pppkkkkkhhyxyxhyxyxyxOhp令当时,有.此时①为p阶Taylor方法.p=1时即为Euler公式.称之为Taylor级数法.其中例2:取步长h=0.1,用一阶、二阶和四阶Taylor方法求解下列初值问题①2()12!!ppkkkkkhhyyhyyyp()()(),0,1,2,,iikkyyxip()()()iikkyyx111()()pkkyxyOh解:(1)一阶Taylor法21,0.2(0)1yyxy210.1kkkyyyk01234yk+11.11.2211.370081.557791.80046(2)二阶Taylor法232231()220.10.122!kkkkyyyyyyyyyk01234yk+11.111.246891.421751.652631.97088(3)四阶Taylor法324(4)3523423451(2)6624240.10.10.10.126242!3!4!kkkkkkyyyyyyyyyyyyyyyk01234yk+11.11111.249961.428481.666441.99942记由得称为[xk,xk+1]上的平均斜率.故2、Runge-Kutta方法11()()()()kkkkyxyxyxxh1()()()()(,())kkkyxyxhyyxhfy*(,())Kfy只要对K*提供不同的算法,就会得出不同的计算公式.如取则得改进的Euler公式,它是利用xk,xk+1两点的斜率值K1,K2的算术平均值作为K*,精度比Euler法高.则得Euler公式;取*1()()kkyxyxhK*(,)kkKfxy*1212111(),(,),(,)2kkkkKKKKfxyKfxyhKRunge-Kutta法的基本思想:设法在[xk,xk+1]内多预报几个点的斜率,再将它们的加权平均值作为平均斜率K*一般显式Runge-Kutta公式为:1111(,1,2,,rkkiiiiikikijjjyyhckkfxhyhkir其中为待定参数,且.称为r级Runge-Kutta方法计算公式.,,iiijc10②注:式中待定参数的确定:先将②式右端在(xk,yk)处展成h的幂级数(即将yk+1展成h的幂级数);再将y(xk+1)作Taylor级数展开;最后比较两式中hk(k=0,1,2,…)的系数,以确定出所有待定参数.即可得p个方程,从而确定出待定参数.代入表达式即可得到计算公式.如果要求两个表达式的前p+1项完全重合,即局部