《微分方程初值问题的数值解法讲稿.ppt》由会员分享,可在线阅读,更多相关《微分方程初值问题的数值解法讲稿.ppt(52页珍藏版)》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、关于微分方程初值问题的数值解法第一页,讲稿共五十二页哦 包含自变量、未知函数及未知函数的导数或微分的方程称为包含自变量、未知函数及未知函数的导数或微分的方程称为微分微分方程方程。在微分方程中在微分方程中,自变量的个数只有一个自变量的个数只有一个,称为常微分方程称为常微分方程。自变量的个数为两个或两个以上的微分方程叫自变量的个数为两个或两个以上的微分方程叫偏微分方程偏微分方程。微分方程中出现的未知函数最高阶导数的阶数称为微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数微分方程的阶数。如果未知函数如果未知函数y y及其各阶导数及其各阶导数)(,nyyy 8.18.1 引引 言言8.1.1
2、8.1.1 微分方程知识回顾微分方程知识回顾都是一次的都是一次的,则称它是则称它是线性线性的的,否则称为否则称为非线性非线性的。的。第二页,讲稿共五十二页哦 在高等数学中,对于常微分方程的求解,给出了一些典型在高等数学中,对于常微分方程的求解,给出了一些典型方程求解析解的基本方法:方程求解析解的基本方法:一阶:一阶:可分离变量法、齐次方程、一阶线性方程以及伯努利方可分离变量法、齐次方程、一阶线性方程以及伯努利方程程 高阶:高阶:可降阶方程、可降阶方程、常系数齐次线性方程的解法、常系数非齐次线常系数齐次线性方程的解法、常系数非齐次线性方程的解法性方程的解法 但能求解的常微分方程仍然是有限的,大多
3、数的常微分方程是求但能求解的常微分方程仍然是有限的,大多数的常微分方程是求不出解析解的。不出解析解的。22yxy 这个一阶微分方程就不能用初等函数及其积分来表这个一阶微分方程就不能用初等函数及其积分来表达它的解。达它的解。8.1.18.1.1 微分方程知识回顾微分方程知识回顾 例如例如 第三页,讲稿共五十二页哦 从实际问题当中归纳出来的微分方程,通常主要依靠数值解法来从实际问题当中归纳出来的微分方程,通常主要依靠数值解法来解决。本章主要讨论解决。本章主要讨论一阶常微分方程初值问题一阶常微分方程初值问题 00)(),(yxyyxfy(8.1(8.1)在区间在区间 可以证明可以证明,如果函数在带形
4、区域如果函数在带形区域 R=axb,-R=axb,-y y 内连续,内连续,且关于且关于y y满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件,即存在常数条件,即存在常数L(L(它与它与x,yx,y无关无关)使使 2121),(),(yyLyxfyxf 对对R R内任意两个内任意两个 都成立都成立,则方程则方程(8.1)(8.1)的解的解 在在 a,b 上存在且唯一。上存在且唯一。21,yy)(xyy 上的数值解法上的数值解法。axb第四页,讲稿共五十二页哦 常微分方程初值问题常微分方程初值问题(8.1)(8.1)式的数值解法,首先要算出精确解式的数值解法,首先要算出精
5、确解y(x)y(x)在区间在区间 a,ba,b 上的一系列离散节点上的一系列离散节点 处的函数值处的函数值 的近似值的近似值:相邻两个节点的间距相邻两个节点的间距 称为步长,本章总称为步长,本章总是假定是假定h h为定数,称为定步长,这时节点可表示为为定数,称为定步长,这时节点可表示为bxxxxann110)(,),(),(10nxyxyxynyyy,10iixxh1niihxxi,2,1 ,08.1.28.1.2 数值方法的基本思想数值方法的基本思想 1 1、数值解法需要把连续性的问题加以离散化,从而求出、数值解法需要把连续性的问题加以离散化,从而求出离散节点处的数值解。离散节点处的数值解。
6、第五页,讲稿共五十二页哦 描述这类算法,要求给出用已知信息描述这类算法,要求给出用已知信息 计算计算 的递推公式。的递推公式。建立这类递推公式的基本方法是在这些节点上用数值建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法,对初值问题积分、数值微分、泰勒展开等离散化方法,对初值问题中的导数中的导数 进行不同的离散化处理进行不同的离散化处理。021,yyyyiii1iyy00)(),(yxyyxfy 2 2、数值解法的基本特点是采用、数值解法的基本特点是采用“步进式步进式”:即求解过程按照递:即求解过程按照递推公式顺着节点排列的次序一步一步地向前推进。推公式顺着节
7、点排列的次序一步一步地向前推进。8.1.28.1.2 数值方法的基本思想数值方法的基本思想第六页,讲稿共五十二页哦 递推公式通常有两类,递推公式通常有两类,一类是计算一类是计算y yi+1i+1时只用到时只用到x xi+1i+1,x,xi i 和和y yi i,即前一步的值,因此有,即前一步的值,因此有了初值以后就可以逐步往下计算,此类方法称为了初值以后就可以逐步往下计算,此类方法称为单步法,单步法,其代表是其代表是龙格龙格库塔法。库塔法。另一类是计算另一类是计算y yi+1i+1时,除用到时,除用到x xi+1i+1,x,xi i和和y yi i以外,还要用到以外,还要用到 ,即前面,即前面
8、k k步的值,此类方法称为多步的值,此类方法称为多步法,步法,其代表其代表是亚当斯法。是亚当斯法。),2,1(,kpyxpipi8.1.28.1.2 数值方法的基本思想数值方法的基本思想第七页,讲稿共五十二页哦一、一、EulerEuler方法及其改进方法及其改进 将将 a,b n 等分等分,记记,(0,1,)kbahxakhknn微分法微分法:111100()()()()()(,)0,1,1()kkkkkkkkkkky xy xy xy xy xxxhyyh f xykny xy积分法积分法:11(1)()()(,()kkxkkxy xy xf t y tdt积分项利用矩形公式计算积分项利用矩
9、形公式计算 11(,()(,)()()(,)kkxkkkkkkxf t y t dth f xyy xy xh f xy1.1.显式显式EulerEuler方法方法()第八页,讲稿共五十二页哦TaylorTaylor公式推导公式推导:2111()()()(),2(,)0,1,1kkkkkkkkkkkhy xy xhy xyxxyyh f xykn第九页,讲稿共五十二页哦Oyxy=y(x)(x1,y1)p1p0 x0 x1x2xixi+1xn-1xnpipi+1pn-1pnp1p2pipi+1pn-1pn切线切线p0p1的斜率为的斜率为f(x0,y0)p2(x2,y2)欧拉公式的几何意义:欧拉公
10、式的几何意义:第十页,讲稿共五十二页哦EulerEuler法的求解过程是法的求解过程是:从初始点从初始点P P0 0(即点即点(x(x0 0,y,y0 0)出发出发,作积分曲线作积分曲线y=y(x)y=y(x)在在P P0 0点上切线点上切线 (其斜率为其斜率为 ),),与直线与直线x=xx=x1 1相相交于交于P P1 1点点(即点即点(x(x1 1,y,y1 1),),得到得到y y1 1作为作为y(xy(x1 1)的近似值)的近似值)10PP),()(000yxfxy)(,(0000 xxyxfyy)(,(010001xxyxfyy这样就获得了这样就获得了P P1 1点的坐标。点的坐标。
11、当当 时时,得得1xx 重复以上过程重复以上过程,就可获得一系列的点就可获得一系列的点:p:p1 1,p,p2 2,p,pn n,相应的可求出相应的可求出y y1 1,y,y2 2,y,yn n,取取nnyxy)(第十一页,讲稿共五十二页哦从图形上看从图形上看,就获得了一条近似于曲线就获得了一条近似于曲线y=y(x)y=y(x)的折线的折线nPPPP321通常取通常取 (常数常数),),则则EulerEuler法的计算格式法的计算格式 hhxxiii1)(),(001xyyyxhfyyiiii i=0,1,n (7.2)第十二页,讲稿共五十二页哦2.2.梯形法梯形法 称之为梯形公式称之为梯形公
12、式.这是一个隐式公式这是一个隐式公式,通常用迭代法求解通常用迭代法求解.具体做法具体做法:取取 先用先用EulerEuler法求出初值法求出初值 ,即即 ,将其代入梯形将其代入梯形公式的右端公式的右端,使之转化为显式公式使之转化为显式公式,即即 注注:当当 f(x,y)关于关于y满足满足Lipschitz条件且步长条件且步长h 满足满足 直至满足直至满足:若采用梯形公式计算若采用梯形公式计算()中的积分项中的积分项,则有则有111111()()(,()(,()2(,)(,)2kkkkkkkkkkkkhy xy xf xy xf xy xhyyf xyf xy(0)1ky(0)1(,)kkkky
13、yh f xy(1)()111(,)(,)2llkkkkkkhyyf xyf xy(1)()11|llkkyy(1)11lkkyy类似地类似地,可得可得23,kkyy()第十三页,讲稿共五十二页哦112h L 时时,迭代格式迭代格式()收敛收敛.3.3.改进的改进的EulerEuler方法方法 把把EulerEuler法作为预报法作为预报(称为预估公式称为预估公式),),把隐式的梯形公式把隐式的梯形公式作为校正作为校正(称为校正公式称为校正公式),),则得改进的则得改进的EulerEuler方法方法:1111(,)(,)(,)2kkkkkkkkkkyyh f xyhyyf xyf xy或或11
14、(,)(,(,)2kkkkkkkkhyyf xyf xyh f xy也称为预估也称为预估-校正法校正法.第十四页,讲稿共五十二页哦11212112(,)(,)kkkkkkyyKKKh fxyKh fxh yhK有时为了方便有时为了方便,预估预估-校正格式也写成下面形式校正格式也写成下面形式:第十五页,讲稿共五十二页哦 改进的欧拉公式比欧拉公式精度高的原因是:改进欧拉公式改进的欧拉公式比欧拉公式精度高的原因是:改进欧拉公式用梯形面积代替曲边梯形面积,而欧拉公式用矩形面积代替曲边用梯形面积代替曲边梯形面积,而欧拉公式用矩形面积代替曲边梯形面积。梯形面积。数值积分的梯形公式比矩形公式的精度高。数值积
15、分的梯形公式比矩形公式的精度高。第十六页,讲稿共五十二页哦二、单步法的局部截断误差及精度二、单步法的局部截断误差及精度 Def 1:先假设先假设 ,再估计误差再估计误差这种误差称为单步迭代法在这种误差称为单步迭代法在 xk+1处的局部截断误差处的局部截断误差.()kky xy,11()k hkkRy xyDef 2:若某种数值方法的局部截断误差为若某种数值方法的局部截断误差为 ,则称该数值方法的则称该数值方法的精度为精度为P 阶的阶的.1()ph注注:通常情况下通常情况下,P 越大越大,h 越小越小,则截断误差越小则截断误差越小,数值方法越精确数值方法越精确.第十七页,讲稿共五十二页哦所以所以
16、EulerEuler方法为一阶方法方法为一阶方法.而而 211211()()()(),2()()(,()()kkkkkkkkkkky xy xh yxh yxxy xy xh fxy xO h122,11(,),()()()2kkkkk hkkkyyh fxyhRy xyyO h设设 1 10 0.Euler.Euler方法是一阶方法方法是一阶方法.()kkyy x第十八页,讲稿共五十二页哦2 20 0.梯形法是二阶方法梯形法是二阶方法.2312311123111()()()()()2()()()()()22()()()()()()()2kkkkkkkkkkkkkkhy xy xhy xyxO
17、 hhy xy xhy xyxO hhy xy xh y xy xyxyxO hTaylorTaylor展开展开 第十九页,讲稿共五十二页哦将将 代入上式代入上式,得得 1()()()kkyxyxO h311()()()()()2kkkkhy xy xy xy xO h11131131111()()()()2(,)(,)()2(,()(,)()2kkkkkkkkkkkkkkhy xyy xy xy xyhf xyf xyO hhf xy xf xyO h11111(,)11(,()(,)|()kkkkkxkkff xy xf xyy xyy而而代入上式得代入上式得:13(,)11(1|)()(
18、)2kxkkhfy xyO hy当当h充分小时充分小时,若若 ,则可选取则可选取 h,使得使得|1fLy第二十页,讲稿共五十二页哦故梯形法的精度为故梯形法的精度为2.同样可以证明同样可以证明改进的改进的EulerEuler法也是二阶方法法也是二阶方法.梯形法的梯形法的局部截断误差局部截断误差为为:3,111()(),12k hkkkkkkhRy xyfxx 1(,)|12kxhfy从而从而111(,)(,)33,11(,)11|21|2()(1|)()()2kkkxxk hkkxhfhfyyhfRy xyO hO hy第二十一页,讲稿共五十二页哦例例1:1)0(20,2yxxyydxdy取步长
19、取步长 h=2/10,2/20,2/30,2/40,分别用欧拉法、改进的欧拉法和梯形法求分别用欧拉法、改进的欧拉法和梯形法求解解.解解:记记 f(x,y)=y x y2,xk=k h (k=0,1,2,n)(1).Euler法法:yk+1=yk+h(yk xk yk2)(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法法:221
20、2210()()2(),0,1,1kkkkkkkkkkkkkkhyyyx yyh yx yxyh yx yknyk01234yk+11.19121.343841.423481.419051.3473k56789yk+11.237261.114240.9941510.8847510.788666(3).梯形法梯形法(计算过程略计算过程略)221111(),0,1,2kkkkkkkkhyyyx yyxykn(1)2()()21111()(),0,1,2,2lllkkkkkkkkhyyyx yyxyl第二十三页,讲稿共五十二页哦 n 10 20 30 40 h 0.2 0.1 0.0667 0.05
21、误差误差 0.1059 0.0521 0.0342 0.0256Euler法误差法误差:改进的改进的Euler法误差法误差:n 10 20 30 40 h 0.2 0.1 0.0667 0.05误差误差 0.0123 0.0026 0.0011 5.9612e-004第二十四页,讲稿共五十二页哦-101234500.511.5预预-校方法校方法,h=0.2时时误差最大值误差最大值:0.0123-101234500.511.52欧拉方法欧拉方法,h=0.2时时误差最大值误差最大值:0.1059xexxy 211)(解析解解析解:第二十五页,讲稿共五十二页哦三、三、Runge-Kutta 方法方法
22、1 1、Taylor 级数级数法法 设初值问题设初值问题 有解有解 y(x),由由Tayler公式得公式得:00(,),(),yf x yy xyaxb2()11()()()()()()2!pppkkkkkhhy xy xhy xyxyxO hp令令当当 时时,有有 .此时为此时为 p 阶阶Taylor方法方法.p=1时即为时即为Euler公式公式.称之为称之为Taylor级数法级数法.其中其中例例2:取步长取步长 h=0.1,用一阶、二阶和四阶用一阶、二阶和四阶Taylor方法求解下列初值问方法求解下列初值问题题2()12!ppkkkkkhhyyh yyyp()()(),0,1,2,iikk
23、yyxip()()()iikkyyx111()()pkky xyO h第二十六页,讲稿共五十二页哦解解:(1)一阶一阶Taylor法法21,0.2(0)1yyxy210.1kkkyyyk01234yk+11.11.2211.370081.557791.80046(2)二阶二阶Taylor法法232231()220.10.122!kkkkyyyyyyyyyk01234yk+11.111.246891.421751.652631.97088第二十七页,讲稿共五十二页哦(3)四阶四阶Taylor法法324(4)3523423451(2)6624240.10.10.10.126242!3!4!kkkk
24、kkyyy yyyy yyyyyyyyk01234yk+11.11111.249961.428481.666441.99942第二十八页,讲稿共五十二页哦记记由由得得称为称为xk,xk+1上的平均斜率上的平均斜率.故故2 2、Runge-Kutta方法方法11()()()()kkkky xy xyxxh1()()()()(,()kkky xy xhyy xh fy*(,()Kfy只要对只要对K*提供不同的算法提供不同的算法,就会得出不同的计算公式就会得出不同的计算公式.如取如取则得改进的则得改进的Euler公式公式,它是利用它是利用xk,xk+1两点的斜率值两点的斜率值K1,K2 的算术平均的
25、算术平均值作为值作为K*,精度比精度比Euler法高法高.则得则得Euler公式公式;取取*1()()kky xy xhK*(,)kkKf xy*1212111(),(,),(,)2kkkkKKKKf xyKf xyhK第二十九页,讲稿共五十二页哦Runge-Kutta法的法的基本思想基本思想:设法在设法在xk,xk+1内多预报几个点的斜率内多预报几个点的斜率,再将它们的加权平均值作为再将它们的加权平均值作为平均斜率平均斜率K*一般显式一般显式Runge-Kutta公式公式为为:1111.,1,2,rkkiiiikkikiki jjjyyhc kkhfx yfxh yhkir其中其中 为待定参
26、数为待定参数,且且 .称为称为r 级级Runge-Kutta方法计算公式方法计算公式.,iii jc 10第三十页,讲稿共五十二页哦即可得即可得 p 个方程个方程,从而确定出待定参数从而确定出待定参数.代入表达式即可得到计算公式代入表达式即可得到计算公式.如如果要求两个表达式的前果要求两个表达式的前p+1项完全重合项完全重合,即局部截断误差达到即局部截断误差达到 ,则则称式为称式为 p 阶阶 r 级级的的Runge-Kutta方法方法.常用的是常用的是 r=2,3,4 级的级的R-K方法方法,且适当选取参数使得且适当选取参数使得 p=r.如要求如要求:2311232()11112!3!()()
27、()()()()2!kkpppkkkkkyyrhr hr hhhy xy xhy xyxyxO hp()12(),(),()pkkpkryxryxryx1()pO h注注:式中待定参数的确定式中待定参数的确定:先将式右端在先将式右端在(xk,yk)处展成处展成h的幂级数的幂级数(即将即将 yk+1 展成展成 h 的幂级数的幂级数);再将再将 y(xk+1)作作Taylor 级数展开级数展开;最后比较两式中最后比较两式中hk(k=0,1,2,)的系数的系数,以确定出所有待定参数以确定出所有待定参数.第三十一页,讲稿共五十二页哦Runge-Kutta方法的推导方法的推导(以以r=2为例为例):当当
28、r=2 时时11 122122211()(,)(,)kkkkkkyyh c kc kkf xykf xh yhk1122211(,)(,)kkkkkkyyhc f xyhc f xh yhk记记(,),(,),(,)kkxxkkyykkff xyffxyffxy第三十二页,讲稿共五十二页哦则则1222112221 111 12223122221 1,(,)(,)()(),()()()()kkkkxykkkxykfkf xh yhkf xyhfk fO hyyc kc khyccf hcfk fhO h又又(,),kkxyyf xyfyff f23123()()()()()2()()2kkkkk
29、xyhy xy xh yxyxO hhyh fff fO h第三十三页,讲稿共五十二页哦这是一个四个参数三个方程的非线性方程组这是一个四个参数三个方程的非线性方程组.它有一个自由度它有一个自由度.称满足上称满足上述方程组的一族公式为二级二阶述方程组的一族公式为二级二阶Runge-Kutta方法方法.为使局部截断误差为为使局部截断误差为 ,比较上述两式右端同次幂系数比较上述两式右端同次幂系数,应应取取3()O h122222111/21/2cccc第三十四页,讲稿共五十二页哦(1)常用的二阶常用的二阶Runge-Kutta方法方法:预估预估-校正算法校正算法(2)122211/2,1cc1121
30、21()/2(,)(,)kkkkkkyyh kkkfxykfxh yhk122210,1,1cc1212122(,)(,)kkkkhhkkyyhkkfxykfxyk中间点方法中间点方法 第三十五页,讲稿共五十二页哦注注:二级二级Runge-Kutta方法的精度最高是二阶的方法的精度最高是二阶的,不可能达到三阶不可能达到三阶.要提高要提高计算方法的阶计算方法的阶,就必须增加预报点就必须增加预报点.常用的三阶常用的三阶Runge-Kutta方法方法(r=3):(1)Heun(休恩休恩)方法方法 (3)122211/4,3/4,2/3cc1121222133(3)/4(,)(,)kkkkkkyyh
31、kkkfxykfxh yhk11231112122312(4)/6(,)(,)(,2)kkkkkkkkyyh kkkkfxykfxh yhkkfxh yhkhk三阶三阶Kutta方法方法 第三十六页,讲稿共五十二页哦(1)三阶三阶Heun方法方法 标准标准(经典经典)四阶四阶Runge-Kutta方法方法 (2)1131112133223233(3)/4(,)(,)(,)kkkkkkkkyyh kkkfxykfxh yhkkfxh yhk常用的四阶常用的四阶Runge-Kutta方法方法(r=4):11234111212211322243(22)/6(,)(,)(,)(,)kkkkkkkkkk
32、yyh kkkkkfxykfxh yhkkfxh yhkkfxh yhk第三十七页,讲稿共五十二页哦(2)称为称为Gill(吉尔吉尔)方法方法 1123411121222122131222222242322(22)(22)/6(,)(,)(,)(,)kkkkkkkkkkyyh kkkkkfxykfxh yhkkfxh yhkhkkfxh yhkhk注注:从理论上讲从理论上讲,可以构造任意高阶的计算方法可以构造任意高阶的计算方法.但事实上但事实上,精度的阶数精度的阶数与预报点的个数之间并非等量关系与预报点的个数之间并非等量关系.预报点的个数预报点的个数 r123456789r 10精度的阶数精度
33、的阶数123445667 r-2一般情况下一般情况下,四阶四阶Runge-Kutta方法已可满足精度要求方法已可满足精度要求.第三十八页,讲稿共五十二页哦例例3:用经典用经典Runge-Kutta方法求解下列初值问题方法求解下列初值问题(取取 h=0.1)2,01(0)1yxyxy 解解:0010.1,0,1,0.1(0,1,9)kkhxyxxk标准标准Runge-Kutta公式为公式为:1123412132430.1(22)/622(0.05)0.052(0.05)0.052(0.1)0.1kkkkkkkkkkyykkkkkxykxykkxykkxyk计算结果见下表计算结果见下表.为比较在相
34、同计算量条件下近似解的精度为比较在相同计算量条件下近似解的精度,表中列表中列出了出了Euler法法(h=0.025)和改进的和改进的Euler法法(h=0.05)在相应节点上的计算结在相应节点上的计算结果果.第三十九页,讲稿共五十二页哦xiEuler法法h=0.025改进改进Euler法法h=0.05经典经典R-K法法h=0.1准确解准确解0.11.1114391.1153801.1155121.1155130.21.2552091.2639141.2642081.2642080.31.4346671.4490891.4495761.4495760.41.6535171.6747561.675
35、4731.6754740.51.9158491.9451711.9461621.9461640.62.2261782.2650402.2663542.2663560.72.5894852.6395612.6412552.6412580.83.0112713.0744793.0766193.0766230.93.4976063.5761443.5788043.5788091.04.0551924.1515734.1548394.154845注注:用表中每种方法计算用表中每种方法计算 yi 都需要计算四次都需要计算四次 f 的值的值,即它们的计算量基本即它们的计算量基本相等相等.第四十页,讲稿共五
36、十二页哦四、四、单步单步法的进一步讨论法的进一步讨论收敛性、相容性与稳定性收敛性、相容性与稳定性注注:由定义可知由定义可知,数值方法的收敛性并不涉及计算过程的舍入误差数值方法的收敛性并不涉及计算过程的舍入误差,只与方法只与方法的截断误差有关的截断误差有关.若格式收敛若格式收敛,则整体截断误差必趋于零则整体截断误差必趋于零.Def:(整体截断误差整体截断误差)称称 为某一数值方法在点为某一数值方法在点 xk 处的整体截断误差处的整体截断误差.它不仅与它不仅与 xk 有关有关,也也与与xk-1,xk-2,x1,x0 有关有关.()kkkey xy则称该单步法收敛则称该单步法收敛.Def:对满足解存
37、在唯一性条件的初值问题对满足解存在唯一性条件的初值问题(1),如果一个显式单步如果一个显式单步法法(3)产生的近似解对于任一固定的产生的近似解对于任一固定的 ,均有均有 00,),xx bxxnh0lim()nhyy x1.收敛性收敛性第四十一页,讲稿共五十二页哦由于由于 ,且且 关于关于 y 满足满足Lipschitz条件条件,得得 则存在常数则存在常数 c 0 使得使得 且单步法中函数且单步法中函数 关于关于 y 满足满足Lipschitz条件条件,则则 定理定理1:若初值问题的一个单步法的局部截断误差为若初值问题的一个单步法的局部截断误差为 记记 证证:由局部截断误差的定义知由局部截断误
38、差的定义知1,()(1)pn hRO hp(,)x y h111()()pnnney xyO h1,1()()(,(),)()pn hnnnnRy xy xhxy xhO h()(,(),)nnnyy xhxy xh11|()|pny xych1(,)nnnnyyhxyh(,)x y h1|()|(,(),)(,)|(1)|()|nnnnnnnnnyyy xyhxy xhxyhhLy xy第四十二页,讲稿共五十二页哦故故 从而有从而有 故故 若若 y(x0)=y0,则则e0=0,由不等式由不等式 得得 111111|()|()|(1)|nnnnnpney xyy xyyychhLe111112
39、112+10+11+10|(1)(1)|1(1)(1)|1(1)(1)(1)(1)|(1)1(1)|(1)1ppnnpnpnnnpnechhLchhLechhLhLechhLhLhLhLehLchhLehL210112hLhLhLhLe!0(1)nnhLhLe+(+1)1(+1)11|1(1)1nhLppnhLnecechhehLL第四十三页,讲稿共五十二页哦设单步法为设单步法为 注注:定理表明定理表明,数值方法的整体截断误差比局部截断误差低一阶数值方法的整体截断误差比局部截断误差低一阶.收敛的方法收敛的方法至少是一阶方法至少是一阶方法.在该定义条件下在该定义条件下,Euler方法是一阶的方法
40、是一阶的,预估预估-校正方法校正方法是二阶是二阶.当当f(x,y)关于关于 y 也满足也满足Lipschitz条件条件,r 级级Runge-Kutta方法中的方法中的 关于关于 y 也满足也满足Lipschitz条件条件,故定理中的条件得到满足故定理中的条件得到满足,解的收敛性解的收敛性得到保证得到保证.由于由于R n,h0(h0),且且 xn为任意点为任意点,故该式相当于用近似方程故该式相当于用近似方程 当当x=xn+1固定时固定时,所以有所以有 10(+1)nnhxxba+=-()11|1pb-a Lpncehec hL2.相容性相容性第四十四页,讲稿共五十二页哦通过在通过在 x=xn 处
41、求解近似方程而获得原方程的近似解处求解近似方程而获得原方程的近似解.因此因此,必须要求当必须要求当h0 时时,近似方程应逼近于原方程近似方程应逼近于原方程.来代替来代替 因此因此,要使要使 h0 时时,近似方程的极限状态为原微分方程近似方程的极限状态为原微分方程,需且只需下列需且只需下列极限成立极限成立:由于由于 由于假设由于假设 是连续函数是连续函数,故上式可表示为故上式可表示为 Def:如果当如果当h0时时,近似方程逼近微分方程近似方程逼近微分方程,则称数值公式则称数值公式与原微分方程相容与原微分方程相容.相容的充要条件相容的充要条件:第四十五页,讲稿共五十二页哦事实上事实上:Remark
42、:可以证明若单步法的阶大于或等于可以证明若单步法的阶大于或等于1,则单步法与微分方程相容则单步法与微分方程相容;反之反之,如果单步法与微分方程相容如果单步法与微分方程相容,且且 关于关于 y 满足满足Lipschitz条件条件,则单步则单步法至少为一阶方法法至少为一阶方法.(h0)(1)若单步法的阶大于或等于若单步法的阶大于或等于1,由由 知知 即单步法与微分方程相容即单步法与微分方程相容.故有故有 (2)如果单步法与微分方程相容如果单步法与微分方程相容,且且 关于关于y 满足满足Lipschitz条件条件,则则第四十六页,讲稿共五十二页哦关于关于 单步法的收敛性以及收敛性定理都是在计算过程中
43、无任何舍入误差的前单步法的收敛性以及收敛性定理都是在计算过程中无任何舍入误差的前提条件下建立的提条件下建立的,但在实际计算时通常会有舍入误差及其积累但在实际计算时通常会有舍入误差及其积累,数值求解数值求解微分方程的过程是一个递推公式微分方程的过程是一个递推公式,必须考必须考 即与微分方程相容的单步法至少为一阶方法即与微分方程相容的单步法至少为一阶方法.Remark:在定理条件下在定理条件下,Euler方法、预估方法、预估-校正方法以及校正方法以及Runge-Kutta方方法都与原微分方程相容法都与原微分方程相容.中连续中连续,且关于变量且关于变量y 满足满足Lipschitz条件条件,则单步法
44、收敛的充要条件为相容性则单步法收敛的充要条件为相容性条件成立条件成立.Th1.设增量函数设增量函数 在区域在区域 3.稳定性稳定性第四十七页,讲稿共五十二页哦 如果数值方法在计算过程中舍入误差的积累越来越大如果数值方法在计算过程中舍入误差的积累越来越大,得不到有效控制得不到有效控制,则称其是不稳定的则称其是不稳定的;反之如果计算结果对初始数据的误差及计算过程中反之如果计算结果对初始数据的误差及计算过程中的误差不敏感的误差不敏感,即舍入误差不增长即舍入误差不增长,则称相应的算法是稳定的则称相应的算法是稳定的.数值方法的数值方法的稳定性有各种定义稳定性有各种定义,这里仅考虑绝对稳定性概念这里仅考虑
45、绝对稳定性概念.虑误差积累能否得到控制虑误差积累能否得到控制.Remark:从上面的定义可知从上面的定义可知,单步法是绝对稳定的单步法是绝对稳定的,与模型方程与模型方程 设某数值方法在节点设某数值方法在节点 xn 处对初值问题的数值解为处对初值问题的数值解为 yn,实际计算实际计算得到的近似解为得到的近似解为 ,称为第称为第 n 步数值解得扰动步数值解得扰动.=-Def:若某种数值方法在计算若某种数值方法在计算 yn 时有扰动时有扰动 但在计算后面的但在计算后面的 ym(m n)由由 引起的误差引起的误差 满足满足 则称该数值方法则称该数值方法是绝对稳定的是绝对稳定的.设设 f(x,y)关于关
46、于 y 满足满足Lipschitz条件条件,下面仅对典型方程下面仅对典型方程(模型方程模型方程)进行讨论进行讨论.(其中其中为复常数为复常数,且且 )第四十八页,讲稿共五十二页哦(1)Euler显式公式显式公式:几个常用公式的稳定性几个常用公式的稳定性要求误差不增加要求误差不增加,即即 必须必须中的复数中的复数以及所用步长以及所用步长 h 有关有关.若对复平面上的某个区域若对复平面上的某个区域G,当当 时单步法绝对收敛时单步法绝对收敛,则称则称G为单步法的绝对收敛域为单步法的绝对收敛域,G与实与实轴的交集称为绝对稳定区间轴的交集称为绝对稳定区间.显然绝对稳定域越大显然绝对稳定域越大,数值方法的
47、绝数值方法的绝对稳定性越好对稳定性越好.将将Euler显式公式用于模型方程显式公式用于模型方程 ,则有则有设设 yn 有误差有误差 ,参与运算的量为参与运算的量为 由此引起的由此引起的 yn+1 有误差有误差 ,则实际得到近似则实际得到近似 yn+1的量为的量为 ,即即故有故有第四十九页,讲稿共五十二页哦故故当当为实数时为实数时,得绝对稳定区间为得绝对稳定区间为(-2,0).即即 是保持绝对稳定性对步长是保持绝对稳定性对步长 h 所加的限制所加的限制.因此因此Euler法的法的绝对稳定域为绝对稳定域为(Euler法是条件稳定的法是条件稳定的):(2)梯形公式梯形公式:将梯形公式用于模型方程将梯形公式用于模型方程 ,则有则有由于由于设设第五十页,讲稿共五十二页哦当当x=Re(h)0时时,上式右端总小于上式右端总小于1.故梯形法的绝对稳定区域为故梯形法的绝对稳定区域为:Re(h)0,即左半平面即左半平面.(3)四阶经典四阶经典R-K方法方法:四阶经典四阶经典R-K方法用于模型方程方法用于模型方程 ,则有则有第五十一页,讲稿共五十二页哦感谢大家观看感谢大家观看第五十二页,讲稿共五十二页哦