常微分方程数值解法讲稿.ppt

《常微分方程数值解法讲稿.ppt》由会员分享，可在线阅读，更多相关《常微分方程数值解法讲稿.ppt（48页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、常微分方程数值解法第一页，讲稿共四十八页哦 所谓数值解法,就是设法将常微分方程离散化,建立差分方程,给出解在一些离散点上的近似值.a=x0 x1x2xnxN=b其中剖分节点xn=a+nh,n=0,1,N,h称为剖分步长.数值解法就是求精确解y(x)在剖分节点xn上的近似值yny(xn),n=1,2,N.假设初值问题(8.1)的解y=y(x)唯一存在且足够光滑.对求解区域a,b做剖分 我们采用数值积分方法来建立差分公式.11.2 .2 构造数值解法的基本思想构造数值解法的基本思想 在区间xn,xn+1上对方程(8.1)做积分,则有第二页，讲稿共四十八页哦对右边的积分应用左矩形公式，则有第三页，讲

2、稿共四十八页哦梯形公式oxyab左矩形公式y=(x)右矩形公式中矩形公式第四页，讲稿共四十八页哦对右边的积分应用左矩形公式，则有因此，建立节点处近似值yn满足的差分公式称之为EulerEuler公式公式.称为梯形公式梯形公式.若对(8.2)式右边的积分应用梯形求积公式，则可导出差分公式第五页，讲稿共四十八页哦 利用Euler方法求初值问题 解解 此时的Euler公式为称为EulerEuler中点公式中点公式或称双步双步EulerEuler公式公式.若在区间xn-1,xn+1上对方程(8.1)做积分,则有对右边的积分应用中矩形求积公式，则得差分公式例例1的数值解.此问题的精确解是y(x)=x/(

3、1+x2).第六页，讲稿共四十八页哦分别取步长h=0.2,0.1,0.05,计算结果如下第七页，讲稿共四十八页哦hxnyny(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.

4、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.00227第八页，讲稿共四十八页哦Euler中点公式则不然,计算yn+1时需用到前两步的值yn,yn

5、-1,称其为两步方法两步方法,两步以上的方法统称为多步法多步法.在Euler公式和梯形公式中,为求得yn+1,只需用到前一步的值yn,这种差分方法称为单步法单步法,这是一种自开始方法.隐式公式中,每次计算yn+1都需解方程,要比显式公式需要更多的计算量,但其计算稳定性较好.在Euler公式和Euler中点公式中,需要计算的yn+1已被显式表示出来,称这类差分公式为显式公式显式公式,而梯形公式中,需要计算的yn+1隐含在等式两侧,称其为隐式公式隐式公式.第九页，讲稿共四十八页哦 从数值积分的角度来看,梯形公式计算数值解的精度要比Euler公式好,但它属于隐式公式,不便于计算.实际上,常将Eule

6、r公式与梯形公式结合使用:2 改进的改进的Euler方法和方法和Taylor展开方法展开方法 2.1 改进的改进的Euler方法方法第十页，讲稿共四十八页哦 由迭代法收敛的角度看，当 (是给定的精度要求)时,取 就可以保证迭代公式收敛,而当h很小时,收敛是很快的.而且,只要 通常采用只迭代一次的算法:称之为改进的改进的Euler方法方法.这是一种单步显式方法.改进的Euler方法也可以写成第十一页，讲稿共四十八页哦 y=y-2x/y ,0 x1的数值解,取步长h=0.1.精确解为y(x)=(1+2x)1/2.例例2 求初值问题 y(0)=1 解解 (1)利用Euler方法第十二页，讲稿共四十八

7、页哦计算结果如下:(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.4

8、142141.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阶方法阶方法.这里p为非负整数.显然,阶数越高,方法的精度越高.研究差

9、分公式阶的重要手段是Taylor展开式,一元函数和二元函数的Taylor展开式为:第十五页，讲稿共四十八页哦另外,在yn=y(xn)的条件下,考虑到y(x)=(x,y(x),则有 y(xn)=(xn,y(xn)=(xn,yn)=n y(xn)=第十六页，讲稿共四十八页哦 yn+1=yn+h(xn,yn)对Euler方法，有 =yn+(xn,yn)h+O(h2)从而有:y(xn+1)-yn+1=O(h2)所以Euler方法是一阶方法.再看改进Euler方法,因为可得第十七页，讲稿共四十八页哦所以,改进的Euler方法是二阶方法.而从而有:y(xn+1)-yn+1=O(h3)2.3 Taylor展

10、开方法展开方法 设y(x)是初值问题(8.1)的精确解,利用Taylor展开式可得第十八页，讲稿共四十八页哦称之为p阶Taylor展开方法.因此，可建立节点处近似值yn满足的差分公式其中第十九页，讲稿共四十八页哦所以,此差分公式是p阶方法.由于Taylor展开方法涉及很多复合函数(x,y(x)的导数的计算,比较繁琐,因而很少直接使用,经常用它为多步方法提供初始值.然而,Taylor展开方法给出了一种构造单步显式高阶方法的途径.Euler方法可写为 可见,公式的局部截断误差为:y(xn+1)-yn+1=O(hp+1).3 Runge-Kutta方法方法 3.1 Runge-Kutta方法的构造方

11、法的构造第二十页，讲稿共四十八页哦 构造差分公式 改进的Euler方法可写为其中i,i,ij为待定参数.若此公式的局部截断误差为第二十一页，讲稿共四十八页哦由于 yn+1=yn+h1n+h2(n+hxn+hn yn)+O(h3)O(hp+1),称公式为p p阶阶Runge-kuttaRunge-kutta方法方法,简称p p阶阶R-KR-K方法方法.对于p=2的情形,应有 =yn+h(1+2)n+h22(xn+n yn)+O(h3)所以,只要令 1+2=1,2=1/2,2=1/2 (8.4)第二十二页，讲稿共四十八页哦 一般地一般地,参数由参数由(8.4)(8.4)确定的一族差分公式确定的一族

12、差分公式(8.3)(8.3)统称为统称为二阶二阶R-KR-K方法方法.称之为中点公式中点公式,或可写为若取=1,则得1=2=1/2,=1,此时公式(8.3)就是改进的Euler公式;若取1=0,则得2=1,=1/2,公式(8.3)为 高阶R-K公式可类似推导.下面列出常用的三阶、四阶R-K公式.第二十三页，讲稿共四十八页哦 四阶标准四阶标准R-KR-K公式公式 三阶三阶R-KR-K公式公式第二十四页，讲稿共四十八页哦 解解 四阶标准R-K公式为例例3 用四阶标准R-K方法求初值问题 y=y-2x/y ,0 x1 y(0)=1的数值解,取步长h=0.2.计算结果如下:第二十五页，讲稿共四十八页哦

13、nxnyny(xn)nxnyny(xn)0120.00.20.41.001.18321.34171.001.18321.34163450.60.81.01.48331.61251.73211.48321.61251.7321 也可以构造隐式R-K方法,其一般形式为称之为p p级隐式级隐式R-KR-K方法方法,同显式R-K方法一样确定参数.如第二十六页，讲稿共四十八页哦是二级二阶隐式R-K方法,也就是梯形公式.但是p级隐式R-K方法的阶可以大于p,例如,一级隐式中点公式为或写为它是二阶方法.3.2 变步长变步长Runge-Kutta方法方法 以p阶R-K方法为例讨论.设从xn以步长h计算y(xn

14、+1)的近似值为 ,局部截断误差为其中,C是与h无关的常数.第二十七页，讲稿共四十八页哦 如果将步长减半,取h/2为步长,从xn经两步计算得到y(xn+1)的近似值记为 ,其局部截断误差为于是有从而,得到事后误差估计可见,当 成立时,可取 .否则,应将步长再次减半进行计算.第二十八页，讲稿共四十八页哦 求解初值问题的单步显式方法可一统一写为如下形式 yn+1=yn+h(xn,yn,h)(8.5)对于Euler方法,有4 单步方法的收敛性和稳定性单步方法的收敛性和稳定性 4.1 单步方法的收敛性单步方法的收敛性 y=(x,y),axb y(a)=其中(x,y,h)称为增量函数增量函数.(x,y,

15、h)=(x,y)对于改进的Euler方法,有 (x,y,h)=1/2(x,y)+(x+h,y+h(x,y)第二十九页，讲稿共四十八页哦 设y(x)是初值问题(8.1)的解,yn是单步法(8.5)产生的近似解.如果对任意固定的点xn,均有y(xn),则称单步法(8.5)是收敛的.可见,若方法(8.5)是收敛的,则当h0时,整体截断误差en=y(xn)-yn将趋于零.定理定理8.18.1 设单步方法(8.5)是p1阶方法,增量函数(x,y,h)在区域axb,-yn)的变化均不超过,则称此差分方法是绝对稳定绝对稳定的.讨论数值方法的稳定性,通常仅限于典型的试验方程 y=y 其中是复数且Re()0.在

16、复平面上,当方法稳定时要求变量h的取值范围称为方法的绝对稳定域绝对稳定域,它与实轴的交集称为绝对稳定绝对稳定区间区间.第三十三页，讲稿共四十八页哦 将Euler方法应用于方程y=y,得到 设在计算yn时产生误差n,计算值yn=yn+n,则n将对以后各节点值计算产生影响.记ym=ym+m,mn,由上式可知误差m满足方程 m=(1+h)m-1=(1+h)m-nn,mn 对隐式单步方法也可类似讨论.如将梯形公式用于方程y=y,则有 yn+1=yn+h/2(yn+yn+1)yn+1=(1+h)yn 可见,若要|m|n|,必须且只须|1+h|1,因此Euler法的绝对稳定域为|1+h|1,绝对稳定区间是

17、-2Re()h0.解出yn+1得 第三十四页，讲稿共四十八页哦类似前面分析,可知绝对稳定区域为由于Re()0,所以此不等式对任意步长h恒成立,这是隐式公式的优点.一些常用方法的绝对稳定区间为方 法方法的阶数稳 定 区 间Euler方法梯形方法改进Euler方法二阶R-K方法三阶R-K方法四阶R-K方法122234(-2 ,0)(-,0)(-2 ,0)(-2 ,0)(-2.51 ,0)(-2.78 ,0)第三十五页，讲稿共四十八页哦 解解 因y0=1,计算得y10=1024,而y(1)=9.35762310-14.例例4 考虑初值问题 y=-30y ,0 x1 y(0)=1取步长h=0.1,利用

18、Euler方法计算y10y(1).y(x)=e-30 x 这是因为h=-3不属于Euler方法的绝对稳定区间.若取h=0.01,计算得y100=3.23447710-16.若取h=0.001,计算得y1000=5.91199810-14.若取h=0.0001,计算得y10000=8.94505710-14.若取h=0.00001,计算得y100000=9.315610-14.第三十六页，讲稿共四十八页哦 单步显式方法的稳定性与步长密切相关,在一种步长下是稳定的差分公式,取大一点步长就可能是不稳定的.收敛性是反映差分公式本身的截断误差对数值解的影响;稳定性是反映计算过程中舍入误差对数值解的影响.

19、只有即收敛又稳定的差分公式才有实用价值.5 线性多步方法线性多步方法 由于在计算yn+1时,已经知道yn,yn-1,及(xn,yn),(xn-1,yn-1),利用这些值构造出精度高、计算量小的差分公式就是线性多步法.5.1 利用待定参数法构造线性多步方法利用待定参数法构造线性多步方法 r+1步线性多步方法的一般形式为第三十七页，讲稿共四十八页哦当-10时,公式为隐式公式,反之为显式公式.参数i,i的选择原则是使方法的局部截断误差为 y(xn+1)-yn+1=O(h)r+2 选取参数,0,1,2,使三步方法 yn+1=yn+h(0n+1n-1+2n-2)这里,局部截断误差是指,在yn-i=y(x

20、n-i),i=0,1,r的前提下,误差y(xn+1)-yn+1.为三阶方法.例例5 解解 设yn=y(xn),yn-1=y(xn-1),yn-2=y(xn-2),则有 第三十八页，讲稿共四十八页哦 n=(xn,y(xn)=y(xn)y(xn+1)=y(xn)+hy(xn)+1/2h2y(xn)+1/6h3y(xn)于是有若使:y(xn+1)-yn+1=O(h4),只要,0,1,2满足:n-1=(xn-1,y(xn-1)=y(xn-1)=y(xn-h)=y(xn)-hy(xn)+1/2h2y(xn)-1/6h3y(4)(xn)+O(h4)n-2=y(xn)-2hy(xn)+2h2y(xn)-4/

21、3h3y(4)(xn)+O(h4)yn+1=y(xn)+h(0+1+2)y(xn)-h2(1+22)y(xn)+h3(1/21+22)y(xn)-h4/6(1+82)y(4)(xn)+O(h5)+1/24h4y(4)(xn)+O(h5)第三十九页，讲稿共四十八页哦 =1,0+1+2=1,1+22=-1/2,1+42=1/3于是有三步三阶显式差分公式设pr(x)是函数(x,y(x)的某个r次插值多项式,则有解之得:yn+1=yn+h/12(23n-16n-1+5n-2)因为 5.2 利用数值积分构造线性多步方法利用数值积分构造线性多步方法 其中 第四十页，讲稿共四十八页哦 选取不同的插值多项式p

22、r(x),就可导出不同的差分公式.下面介绍常用的AdamsAdams公式公式.设已求得精确解y(x)在步长为h的等距节点xn-r,xn上的近似值yn-r,yn,记k=(xk,yk),利用r+1个数据(xn-r,n-r),(xn,n)构造r次Lagrange插值多项式由此,可建立差分公式 1.Adams1.Adams显式公式显式公式 其中 第四十一页，讲稿共四十八页哦由此,可建立差分公式 由于 hrj 则有 称之为r+1r+1步步AdamsAdams显式公式显式公式.第四十二页，讲稿共四十八页哦下面列出几个带有局部截断误差主项的Adams显式公式 r=0 yn+1=yn+hn+(1/2)h2y(

23、xn)2.Adams2.Adams隐式公式隐式公式 r=1 yn+1=yn+(h/2)(3n-n-1)+(5/12)h3y(xn)r=2 yn+1=yn+(h/12)(23n-16n-1+5n-2)+(3/8)h4y(4)(xn)r=3 yn+1=yn+(h/24)(55n-59n-1+37n-2-9n-3)+(251/720)h5y(5)(xn)如果利用r+1个数据(xn-r+1,n-r+1),(xn+1,n+1)构造r次Lagrange插值多项式pr(x),则可导出数值稳定性好的隐式公式,称为AdamsAdams隐式公式隐式公式,其一般形式为第四十三页，讲稿共四十八页哦其中系数为 下面列出

24、几个带有局部截断误差主项的Adams隐式公式 r=0 yn+1=yn+hn+1-(1/2)h2y(xn)r=1 yn+1=yn+(h/2)(n+n+1)-(1/12)h3y(xn)r=2 yn+1=yn+(h/12)(5n+1+8n-n-1)-(1/24)h4y(4)(xn)r=3 yn+1=yn+(h/24)(9n+1+19n-5n-1+n-2)-(19/720)h5y(5)(xn)第四十四页，讲稿共四十八页哦 3.Adams3.Adams预估预估-校正公式校正公式 由显式公式提供一个预估值，再用隐式公式校正一次，求得数值解，称为预估校正方法预估校正方法。校正 yn+1=yn+(h/24)(

25、9n+1+19n-5n-1+n-2)一般预估公式和校正公式都采用同阶公式。例如：预估 yn+1=yn+(h/24)(55n-59n-1+37n-2-9n-3)n+1=(xn+1,yn+1),n=3,4,称为四阶Adams预估校正公式.实际计算时通常用四阶单步方法(如四阶R-K公式)为它提供起始值y1,y2,y3.例例6 用四阶Adams预估校正公式求解初值问题 第四十五页，讲稿共四十八页哦 y=y-2x/y,0 x1 y(0)=1取步长h=0.1.解 用四阶R-K公式提供起始值,计算结果如下xnR-k法yn预估值yn校正值yn精确值y(xn)00.10.20.30.40.50.60.70.80

26、.9111.0954461.1832171.2649121.3415511.4140451.4830171.5489171.6121141.6729141.7315661.3416411.4142131.4832391.5491921.6124501.6733181.73204811.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.732051第四十六页，讲稿共四十八页哦练习题练习题第第250页页 习题习题88-5,8-7,8-8,8-11,8-12,8-13,8-15 第四十七页，讲稿共四十八页哦课间休息课间休息第四十八页，讲稿共四十八页哦