1、第十章常微分方程数值解法-10.2 Euler 方法方法第十章第十章 常微分方法的数值解法常微分方法的数值解法10.2.2 局部误差和方法的阶局部误差和方法的阶10.2.1 Euler 方法及其有关的方法方法及其有关的方法第十章常微分方程数值解法第第10章章 常微分方法的数值解法常微分方法的数值解法教学目的教学目的 1. 1. 掌握解常微分方程的单步法:掌握解常微分方程的单步法:EulerEuler方法、方法、TaylorTaylor方法和方法和Runge-KuttaRunge-Kutta方法;方法;2. 2. 掌握解常微分方程的多步法:掌握解常微分方程的多步法:AdamsAdams步法、步法
2、、SimpsonSimpson方法和方法和MilneMilne方法等;方法等;3. 3. 了解单步法的收敛性、相容性与稳定性;多步了解单步法的收敛性、相容性与稳定性;多步法的稳定性。法的稳定性。教学重点及难点教学重点及难点 重点重点是解常微分方程的单步法:是解常微分方程的单步法:EulerEuler方法、方法、TaylorTaylor方法和方法和Runge-KuttaRunge-Kutta方法和解常微分方程的方法和解常微分方程的多步法:多步法:AdamsAdams步法、步法、SimpsonSimpson方法和方法和MilneMilne方法等;方法等;难点难点是理解单步法的收敛性、相容性与稳定性
3、及多是理解单步法的收敛性、相容性与稳定性及多步法的稳定性。步法的稳定性。第十章常微分方程数值解法1基本概念基本概念 包含自变量、未知函数及未知函数的导数或微分的方包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程。在微分方程中程称为微分方程。在微分方程中, , 自变量的个数只有一个自变量的个数只有一个, , 称为称为常微分方程常微分方程. .。自变量的个数为两个或两个以上的微分。自变量的个数为两个或两个以上的微分方程叫方程叫偏微分方程偏微分方程。微分方程中出现的未知函数最高阶导。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。如果未知函数数的阶数称为微分方程的阶数。如果未
4、知函数 y 及其各阶及其各阶导数导数都是一次的都是一次的, ,则称它是线性的则称它是线性的, ,否则称为非线性的。否则称为非线性的。 )(,nyyy 第十章常微分方程数值解法 在高等数学中,对于常微分方程的求解,给出了一些典在高等数学中,对于常微分方程的求解,给出了一些典型方程求解析解的基本方法,如可分离变量法、常系数齐次型方程求解析解的基本方法,如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。但能求线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是有限的,大多数的常微分方程是不可解的常微分方程仍然是有限的,大多数的常微分方程是不可能给出解析解。能
5、给出解析解。 譬如譬如 22yxy 这个一阶微分方程就不能用初等函数及其积分来表达这个一阶微分方程就不能用初等函数及其积分来表达它的解。它的解。 第十章常微分方程数值解法再如,方程再如,方程 1)0(yyy的解的解 , ,虽然有表可查虽然有表可查, ,但对于表上没有给出的但对于表上没有给出的 值值, ,仍需插值方法来计算仍需插值方法来计算. .xey xe 大多数情况下,常微分方程只能用近似方法求解。大多数情况下,常微分方程只能用近似方法求解。这种近似解法可分为两大类:这种近似解法可分为两大类:一类是一类是近似解析法近似解析法,如级,如级数解法、逐次逼近法等;数解法、逐次逼近法等;另一类是另一
6、类是数值解法数值解法,它给出方,它给出方程在一些离散点上的近似值。程在一些离散点上的近似值。第十章常微分方程数值解法 在具体求解微分方程时,需具备某种定解条件,微分方在具体求解微分方程时,需具备某种定解条件,微分方程和定解条件合在一起组成程和定解条件合在一起组成定解问题定解问题。定解条件有两种:一种是给出积分曲线在初始点的状态,定解条件有两种:一种是给出积分曲线在初始点的状态,称为称为初始条件初始条件,相应的定解问题称为,相应的定解问题称为初值问题初值问题。另一类是给出积分曲线首尾两端的状态,称为另一类是给出积分曲线首尾两端的状态,称为边界条件边界条件,相应的定解问题称为相应的定解问题称为边值
7、问题边值问题。第十章常微分方程数值解法其中其中 x 是质量,是质量,m是离开平衡位是离开平衡位o的距离,的距离,t为时间,为时间,c为弹簧系数。为弹簧系数。 例如例如:弹簧一质量系统的振动问题:弹簧一质量系统的振动问题经一定的简化后可用一个二阶常微经一定的简化后可用一个二阶常微分方程分方程oxmcdtxd22来描述。来描述。mxxoc 000000( )( ):t tttx txdxx txx tdt当弹簧在振动开始时刻时的初始位置和初始速度都确定时,弹簧振动规律唯一确定。这就可以写成一个初值问题2200000()()dxcxd tmxtxxtx第十章常微分方程数值解法我们现在讨论常微分方程的
8、数值解法。先从最简单的我们现在讨论常微分方程的数值解法。先从最简单的一阶一阶常微分方程的初值问题常微分方程的初值问题出发出发开始讨论。开始讨论。问题问题: : 考虑一阶常微分方程初值的问题:考虑一阶常微分方程初值的问题: 00)(),(yxyyxfy求函数求函数0( ),y xaxxb满足满足其中其中),(yxf为已知函数为已知函数, ,0y是已知值是已知值. .(1 .1)(1 .2 ): - ( , ) ( ),( , )( )-2y2yxy2y4xy4f x y4xxy 13dy2yf x y4dxxy 13 例例 方方程程令令:且且给给出出初初值值就就得得到到一一阶阶常常微微分分方方程
9、程的的初初值值问问题题:第十章常微分方程数值解法由常微分方程理论可知:只要上式中的函数由常微分方程理论可知:只要上式中的函数f (x, y)在区域在区域G=axb,y内连续,且关于内连续,且关于 y 满足满足Lipschitz条件,条件,即存在与即存在与 x, y 无关的常数无关的常数L,使,使1212( ,)( ,) (1.3)f x yf x yL yy12,(1.1)(1.2)()Gyyyy x对内 任 意 两 个都 成 立 , 则 方 程的 解必 定 存 在 且 唯 一 。Lipschitz条件条件的几何解释的几何解释: 即一条曲线上任意两点连线的斜率即一条曲线上任意两点连线的斜率的绝
10、对值都有小于某一个数的绝对值都有小于某一个数.下面分析均假定满足上述条件。下面分析均假定满足上述条件。第十章常微分方程数值解法 初值问题初值问题(8.1)(8.1)的数值解法是通过微分方程离散化而给的数值解法是通过微分方程离散化而给出解在出解在某某些些离散点上离散点上( (节点上节点上) )的近似值,为了讨论问题方便的近似值,为了讨论问题方便, ,引入以下引入以下概念。概念。在在,ba上引入节点上引入节点,0nkkx ,210bxxxxan 1 kkkxxh常用等步长常用等步长: :,nabh 则有则有khaxk )., 0(nk (1.1)(1.1)的准确解记为的准确解记为y( (x),),
11、.),(kkkfyxf 记记求初值问题数值解的方法是求初值问题数值解的方法是步进法步进法. .即逐个节点计算即逐个节点计算, ,由由)(kiyi .1 ky计计算算称为步长。称为步长。), 2 , 1(nk )(kxy的近似解的近似解记记为为,ky步进法步进法单步法单步法: :多步法多步法: : 仅仅由仅仅由ky计算计算.1 ky计算计算)1, 1 , 0( ljyjk由由.1 ky共用到共用到l个值个值.即即,11 lkkkyyy称为称为l步法。步法。 第十章常微分方程数值解法( (2)2)理论分析理论分析: :单步法比单步法比1 l的多步法容易分析的多步法容易分析( (稳定性稳定性).).
12、( (3)3)选步长方面选步长方面: : 单步法容易改变步长单步法容易改变步长. .( (4)4)精度精度: :多步法精度高一些多步法精度高一些. .单步法与多步法的区别单步法与多步法的区别: :( (1)1)计算方面计算方面: :l步方法只用于步方法只用于lkyk ,的计算的计算, ,110, lyyy的的计算计算要用其它方法。要用其它方法。), 2 , 1( ,kpyxpipi计算计算yi+1时只用到时只用到xi+1, xi 和和yi,即前一步的值,因此有了初值以即前一步的值,因此有了初值以后就可以逐步往下计算,此类方法称为后就可以逐步往下计算,此类方法称为单步单步法法. 其代表是其代表是
13、龙龙格格库塔法库塔法。另一类是计算另一类是计算yi+1时,除用到时,除用到xi+1 , xi 和和 yi 以外以外,还要用还要用到到 ,即前面,即前面k步的值,此类方法步的值,此类方法称为称为多步法多步法;其代表是其代表是亚当斯法亚当斯法。第十章常微分方程数值解法单步法与多步法又都有单步法与多步法又都有显式方法显式方法和和隐式方法隐式方法之分之分. .计算计算公式公式依依次可写成次可写成: :显式单步法显式单步法: :),(1hyxhyykkkk (1.4)隐式单步法隐式单步法: :),(11hyyxhyykkkkk (1.5)该式右端项含有该式右端项含有,1 ky因此若求因此若求,1 ky需
14、要解方程需要解方程(1.5)(1.5)。注注: :迭代法迭代法显式多步法显式多步法: :),(111hyyyxhyylkkkkkk (1.6)隐式多步法隐式多步法: :),(111hyyyxhyylkkkkkk (1.7)线性多步法线性多步法: :ikliiikliikfhyy 11101 (1.8)1, lk1, lk其中其中101, 110, ll 是独立于是独立于k和和f 的常数。的常数。则则若若, 01 第十章常微分方程数值解法注注: :(1.8)(1.8)关于关于ikikfy ,都是线性的。都是线性的。(1.8)(1.8)是显式是显式 ( (右端不含有右端不含有 ) );1 ky,
15、01 若若则则(1.8)(1.8)是隐式的。是隐式的。数值解法要讨论的内容数值解法要讨论的内容: :方法构造方法构造误差分析误差分析稳定性稳定性 显式与隐式两类方法各有特点,使用显式算法远比隐式算显式与隐式两类方法各有特点,使用显式算法远比隐式算法方便,但考虑数值稳定性等因素,人们常选用隐式算法。法方便,但考虑数值稳定性等因素,人们常选用隐式算法。第十章常微分方程数值解法微分方程离散化常用方法1(,)1.()() nnnnxynnAy xy xdydxxx用差商代替微商11. : ( , )(0 ,1,)nnnnxxxxBdydxf x y dxndx数值积分用数值积分方法离散化第十章常微分方
16、程数值解法11111 ,(), (), ( , ) (, ()(0,)1,nnnnnnnxnnxnnyyy xy xyyyf x y dxhf xy xn用代替对右端积分采用数值积分公式(如左矩形公式等),从而得到关于的方程11 ( , )(0,1,)nnnnxxxxdydxf x y dxndx第十章常微分方程数值解法22C ( ) ()()()()2 ()(, ()()2nnnnnnnnnxy xTaylorhy xhy xhy xy xhy xhf xy xy x在 附近的展开:11 ( ) () ( ,) 0, 1, 2, nnnnnnnhyy xy xyyhf x ynTaylor取
17、的线性部分, 且得的近似值:展开法不仅可得到求数值解的公式,且容易估计截断误差。第十章常微分方程数值解法欧拉欧拉(Euler)方法是解初值问题的最简单的数值方法。对方法是解初值问题的最简单的数值方法。对初值问题(初值问题(1.11.1)(1.2)(1.2)的解的解y=y(x)代表通过点代表通过点 的一条称之为微分方程的积的一条称之为微分方程的积分曲线。积分曲线上每一点分曲线。积分曲线上每一点 的切线的斜率的切线的斜率 等等于函数于函数 在这点的值。在这点的值。 00)(),(yxyyxfy),(00yx),(yx)(xy),(yxf2 欧拉(欧拉(EulerEuler)法)法0,axxb第十章
18、常微分方程数值解法Euler法的求解过程是法的求解过程是:从初始点从初始点P0(即点即点(x0,y0)出发出发,作积分曲线作积分曲线y=y(x)在在P0点上切线点上切线 (其斜其斜率为率为 ),与与x=x1直线直线10PP),()(000yxfxy相交于相交于P1点点( (即点即点( (x1 1, ,y1 1),),得到得到y1 1作为作为y( (x1 1) )的近似值的近似值, ,如上图如上图所示。过点所示。过点( (x0 0, ,y0 0),),以以f( (x0 0, ,y0 0) )为为斜率的切线方程为斜率的切线方程为 当当 时时, ,得得 )(,(0000 xxyxfyy1xx )(,
19、(010001xxyxfyy这样就获得了这样就获得了P1点的坐标。点的坐标。 yxo),(000yxp),(111yxp)(xy0 xa 1x2xbxn 1 1o o 显式显式 EulerEuler方法方法第十章常微分方程数值解法同样同样, 过点过点P1(x1,y1),作积分曲线作积分曲线y=y(x)的切线交直线的切线交直线x=x2于于P2点点,切线切线 的斜率的斜率 =直线方程为直线方程为21PP)(1xy),(11yxf)(,(1111xxyxfyy)(,(121112xxyxfyy当当 时时, ,得得 2xx yxo),(000yxp),(111yxp)(xy),(222yxp0 xa
20、1x2xbxn 第十章常微分方程数值解法当当 时时, ,得得由此获得了由此获得了P P2 2的坐标。重复以上过程的坐标。重复以上过程, ,就可获得一系列的就可获得一系列的点点: :P P1 1, ,P P1 1, , ,P Pn n。对已求得点对已求得点以以 = = 为斜率作直线为斜率作直线 ),(nnnyxP)(nxy),(nnyxf)(,(nnnnxxyxfyy1nxx11(,)()nnnnnnyyf xyxxnnyxy)(取取yxo),(000yxp),(111yxp)(xy),(222yxp),(nnnyxp0 xa 1x2xbxn 第十章常微分方程数值解法 从图形上看从图形上看, ,
21、就获得了一条近似于曲线就获得了一条近似于曲线y=y(x)的折线的折线 。这样这样, ,从从x0 0逐个算出逐个算出对应的数值解对应的数值解 nxxx,21nyyy,21nPPPP321yxo),(000yxp),(111yxp)(xy),(222yxp),(nnnyxp0 xa 1x2xbxn 第十章常微分方程数值解法通常取通常取 ( (常数常数),),则则EulerEuler法的计算格式法的计算格式 hhxxiii1)(),(001xyyyxhfyyiiii i=0,1,n ( 2.1 ) 此外,还可用数值微分、数值积分法和泰勒展开此外,还可用数值微分、数值积分法和泰勒展开法推导法推导Eul
22、erEuler格式:格式:欧拉方法又称为欧拉方法又称为折线法折线法。第十章常微分方程数值解法以以TaylorTaylor展开法为例进行推导展开法为例进行推导将将)(1 kxy在在kxx 点进行点进行TaylorTaylor展开展开,! 2)()(,()()(121 kkkkkkkkkkxxhyxyxfhxyxy 忽略忽略2kh这一高阶项这一高阶项,)(kxy 由于由于因此因此分别用分别用,1 kkyy),(kkkyxff 近似近似),(),(1 kkxyxy),(,(kkxyxf得得,1kkkkfhyy 结合初值条件结合初值条件00()y xy即得即得(2.1).00()yy x,1kkkkf
23、hyy 1, 1 , 0 nk(2.1) )(,()(kkkxyxfxy ),()(yxfdxdyxy (2.2)第十章常微分方程数值解法以以向前差分近似微分法向前差分近似微分法为例进行推导为例进行推导向前差商向前差商kkkhxyxy)()(1 近似近似),(kxy 得得)(,()()(1kkkkkxyxfhxyxy 将近似号改为等号,将近似号改为等号,),(,1kkkkkyxffyy 用用近似近似),(),(1 kkxyxy),(,(kkxyxf并结合初始条件即得并结合初始条件即得(2.1)。(2.3)第十章常微分方程数值解法以以左矩形数值积分法左矩形数值积分法为例进行推导为例进行推导将将(
24、1.1)两端从两端从kx到到1 kx积分积分,得得( , ),(1.1)dyf x ydx 1)(,()()(1kkxxkkdxxyxfxyxy数值积分采用左矩形公式数值积分采用左矩形公式 , 即即)(,()(,(1kkkxxxyxfhdxxyxfkk 由初始条件亦得由初始条件亦得(2.1).oxy)(,(xyxf)(,(kkxyxfkx1 kxkh用用),(,1kkkkkyxffyy 近似近似),(kxy),(,(kkxyxf得得,1kkkkfhyy ),(1 kxy,1kkkkfhyy 1, 1 , 0 nk(2.1) 00()yy x从从x0处的初值处的初值y0 开开始,按(始,按(2.
25、1)可逐)可逐步计算以后各点上步计算以后各点上的值。称(的值。称(2.1)式)式为为显式显式Euler法。法。(2.4)第十章常微分方程数值解法例例8.18.1 用欧拉法解初值问题用欧拉法解初值问题 1)0()6 . 00(2yxxyyy取步长取步长h=0.2 ,h=0.2 ,计算过程保留计算过程保留4 4位小数位小数 解解: h=0.2, x1,2,3=0.2, 0.4,0.6, 欧拉迭代格式欧拉迭代格式 2),(xyyyxf21(,)nnnnnnnnyyhf xyyhyhx y0.2(4)(0,1,2)nnnyx yn当 k=0, x1=0.2时,已知x0=0,y0=1,有 y(0.2)y
26、1=0.21(401)0.8当 k=1, x2=0.4时,已知x1 =0.2, y1 =0.8,有 y(0.4) y2 =0.20.8(40.20.8)0.6144当 k=2, x3 =0.6时,已知x2 =0.4, y2 =0.6144,有 y(0.6) y3=0.20.6144(40.40.6144)=0.4613 xn0 0.20.40.6yn1 0.80.6411 0.4613第十章常微分方程数值解法,! 2)()(,()()(12111 kkkkkkkkkkxxhyxyxfhxyxy 将将)(kxy在在1 kxx点进行点进行TaylorTaylor展开展开忽略忽略2kh项项, ,用用
27、),(,1111 kkkkkyxffyy分别近似分别近似),(),(1 kkxyxy),(,(11 kkxyxf可得隐式可得隐式/ /后退欧拉方法后退欧拉方法: :),(111 kkkkkyxfhyy1, 1 , 0 nk(2.5)2 2o o 隐式隐式/ /后退后退 EulerEuler方法方法)(,()()(1kkkkkxyxfhxyxy )(1 kxy在在kxx 点展开点展开显式显式EulerEuler方法方法,11 kkkkfhyy1, 1 , 0 nk,1 kkkxx ,! 2)(2kkhy 第十章常微分方程数值解法说明说明: :(1)(1)隐式隐式 EulerEuler方法也可用方
28、法也可用kkkhxyxy )()(1近似微分近似微分,)(11 kkxyyxxdxdy向后向后差商差商, ,即即( (自己推自己推) )(2)(2)隐式隐式 EulerEuler方法方法(2.5)(2.5)是关于是关于1 ky需要解方程需要解方程. .的方程,的方程,1ky或者用或者用右矩数值右矩数值求积公式来建立求积公式来建立. .若求若求第十章常微分方程数值解法 隐隐式式算法(算法(2.52.5)常用迭代法来实现,而迭代过程实质上)常用迭代法来实现,而迭代过程实质上是逐步显式化。是逐步显式化。 设用显式Euler格式算出 作迭代初值 ,以此代入(2.5)右端,使之转化为显示,直接计算得:
29、,再用 代入(2.5)右端又有:),()0(1iiiiyxhfyy)0(1iy(1)(0)111(,)iiiiyyhf xy)1(1iy),()1(11)2(1iiiiyxhfyy若迭代过程收敛,则极限值 必为隐式方程的解,从而可获得后退Euler方法的解。如此反复进行可得序列 。)(1hiy, 2 , 1 , 0),()(11)1(1kyxhfyykiiiki1)(1limikikyy第十章常微分方程数值解法如果将(如果将(2.1)和()和(2.5)两式作算术平均,就得)两式作算术平均,就得梯形公式。梯形公式。111()()012iiiiiihyyf xyf xy,i, , 。(2.6)梯形
30、公式也是隐式公式。梯形公式也是隐式公式。 显式与隐式两类方法各有特点,使用显式算法远比隐式算法显式与隐式两类方法各有特点,使用显式算法远比隐式算法方便,但考虑数值稳定性等因素,人们常选用隐式算法。方便,但考虑数值稳定性等因素,人们常选用隐式算法。3 3o o 梯形公式梯形公式),(111 kkkkkyxfhyy1, 1 , 0 nk(2.5),1kkkkfhyy 1, 1 , 0 nk(2.1)第十章常微分方程数值解法此外,还用数值积分推导出梯形方法:此外,还用数值积分推导出梯形方法:将将(1.1)两端从两端从kx到到1 kx积分积分,得得( , ),(1.1)dyf x ydx11()()(
31、 , ( ) (2.7)kkxkkxy xy xf x y x dx用数值积分的用数值积分的梯形方法梯形方法计算其积分项,即计算其积分项,即 1111, ( )( , ( )(, ()2iixiiiiiixxxf x y x dxf x y xf xy x代入代入(2.7)(2.7)式式, ,并用并用yi近似代替式中近似代替式中y(xi), 即可得到即可得到梯形公式梯形公式 ),(),(2111iiiiiiyxfyxfhyy( 2.6 )( 2.6 ) 第十章常微分方程数值解法 由于数值积分的梯形公式比矩形公式的精度高,因此梯由于数值积分的梯形公式比矩形公式的精度高,因此梯形公式(形公式(2.
32、62.6)比欧拉公式)比欧拉公式( 2.1 )( 2.1 )的精度高一个数值方法。的精度高一个数值方法。 从几何上看,梯形公式是取从几何上看,梯形公式是取xi , xi+1区间两端点处斜率区间两端点处斜率的平均斜率。的平均斜率。1iixxy=f(x)xyoAB2Q1QP Euler方法是过Q1(xi ,yi)点以斜率f (xi ,yi)引直线交x=xi+1的点A。后退Euler方法是以点Q2(xi ,yi)处的斜率f (xi+1 ,yi+1)为斜率,从x=xi+1点引直线交Q1于另一点B。 A、B两点都偏离点Q2点,梯形方法就是取A、B两点的中点P作为Q2近似,上图表明梯形方法确实改善了精度。
33、第十章常微分方程数值解法例例8.2 取取h=0.1,用,用Euler方法、隐式方法、隐式Euler方法和梯形方法解方法和梯形方法解。,1)0(1yyxy解解 本题有本题有 如果用如果用Euler方法,由方法,由(2.1)并代入)并代入h=0.1得得 。,11)(0yyxyxf10.10.90.1iiiyxy同理,用隐式同理,用隐式Euler方法方法, 并可整理成并可整理成显式表达显式表达,有,有111(0.10.1)1.1iiiyxy。,1kkkkfhyy ,11 kkkkfhyy第十章常微分方程数值解法以上三种方法及准确解以上三种方法及准确解 的数值结果如表的数值结果如表8-1所示。所示。x
34、exxy)(用梯形公式并可整理成用梯形公式并可整理成显式表达显式表达,有,有。)105.095.01 .0(05.111nnnyxy111(,)(,)2iiiiiihyyf xyf xy第十章常微分方程数值解法表表8-1 Euler方法方法 隐式隐式Euler方法方法 梯形法梯形法 准确解准确解nx 0 1 1 1 10.1 1.000000 1.009091 1.004762 1.0048370.2 1.010000 1.026446 1.018549 1.018731 0.3 1.029000 1.051315 1.040633 1.040818 0.4 1.056100 1.083013
35、 1.070096 1.070320 0.5 1.090490 1.120921 1.106278 1.106531 从表中看到,在从表中看到,在 处,处,Euler方法和隐式方法和隐式Euler方法的误差方法的误差5 . 0nxnnyxy)(2104 . 12106 . 14105 . 2分别是分别是 和和 ,而梯形方法的误差,而梯形方法的误差却是却是第十章常微分方程数值解法例例8.38.3 对初值问题对初值问题 1)0(0yyy证明用梯形公式求得的近似解为证明用梯形公式求得的近似解为 2,2nnhyhnhx 并证明当步长并证明当步长h0时时,yn收敛于精确解收敛于精确解证明证明: 解初值问
36、题的梯形公式为解初值问题的梯形公式为 xe),(),(nnnnnnyxfyxfhyyyyxf),( 211nnnnyyhyy 整理成整理成显式显式 nnyhhy反复迭代反复迭代, ,得到得到 yhhyhhyhhyhhynnnnn.10ynnhhy22 第十章常微分方程数值解法由于由于 ,有,有 nhx xxxxhxhhhxhnhhhhhyeee2121lim22limlim222222000nnhhy22xnhy elim0 证毕证毕 第十章常微分方程数值解法对方程对方程 的两端在区间上的两端在区间上 积分得积分得 ),(yxfy 11,iixx1111()(),( )iixiixy xy x
37、fx y xdx改用改用中矩形公式中矩形公式计算其积分项,即计算其积分项,即 )(,)(,1111iiiixxxyxfxxdxxyxfii代入上式代入上式, ,并用并用yi近似代替式中近似代替式中y(xi)即可得到两步欧拉公式即可得到两步欧拉公式 ),(211iiiiyxhfyy ( 2.8 ) 4 4o o 两步欧拉公式两步欧拉公式第十章常微分方程数值解法 前面介绍过的数值方法前面介绍过的数值方法,无论是欧拉方法无论是欧拉方法, 还是梯形方法,还是梯形方法,它们都是单步法它们都是单步法, 其特点是在计算其特点是在计算yi+1时只用到前一步的信时只用到前一步的信息息yi; 可是公式可是公式(8
38、.7)中除了中除了yi外外, 还用到更前一步的信息还用到更前一步的信息yi-1, 即调用了前两步的信息即调用了前两步的信息, 故称其为两步欧拉公式故称其为两步欧拉公式 。第十章常微分方程数值解法 对于不同的方法,计算值对于不同的方法,计算值yn与准确解与准确解y(xn) 的误的误差各不相同。所以有必要讨论方法的截断误差。我差各不相同。所以有必要讨论方法的截断误差。我们称们称en= yn -y (xn)为某一方法在为某一方法在xn点的点的整体截断误差整体截断误差。显然,显然, en不单与不单与xn这步的计算有关,它与以前各步这步的计算有关,它与以前各步的计算也有关,所以误差被称为整体的。分析和估
39、的计算也有关,所以误差被称为整体的。分析和估计整体截断误差计整体截断误差en是复杂的。为此,我们是复杂的。为此,我们假设假设xn处处的的yn没有误差,即没有误差,即yn=y (xn) ,考虑从考虑从xn到到xn+1这一步这一步的误差,这就是如下的的误差,这就是如下的局部误差局部误差的概念。的概念。10.2.2 局部误差和方法的阶局部误差和方法的阶第十章常微分方程数值解法 衡量求解公式好坏的一个主要标准是求解公式的精衡量求解公式好坏的一个主要标准是求解公式的精度度, , 因此引入局部截断误差和阶数的概念。因此引入局部截断误差和阶数的概念。 定义定义8.2 数值方法的局部截断误差为数值方法的局部截
40、断误差为O(hp+1), 则则称这种数值称这种数值方法的阶数是方法的阶数是P。步长。步长(h N 结束。结束。 10,xx)(21),(),(11cpipiiciiipyyyyxhfyyyxhfyy11, yx0101,yyxx第十章常微分方程数值解法(2)改进欧拉法的流程图)改进欧拉法的流程图 开 始 输 入x0, y0,h , N 1 n x0 + h x1 y0+ h f( x0,y0 ) yp y0+ h f( x1,yp) yc ( yp+ yc) /2 y1 输 出x1, y1 n + 1 n n = N ? x1 x0 y1 y0 结 束 n y 第十章常微分方程数值解法 欧拉公
41、式欧拉公式和和改进欧拉公式分改进欧拉公式分iiipipicicpiyyyyyyyyyyyy9 .0)(1 .091.0)(1 .0905.0211及及两种格式的计算结果分别列表如下两种格式的计算结果分别列表如下别为:别为:000,1xy解:取步长解:取步长h=0.1, iiiyyyyi9 . 0)( 1 . 01ydxdyy1)0( 1 , 0 x并比较两法所得计算结果并比较两法所得计算结果的精度。的精度。例例6:试分别用:试分别用欧拉欧拉公式和改进的欧拉公式求解:公式和改进的欧拉公式求解:第十章常微分方程数值解法 由表可见,与精确解 相比,改进的Euler公式的精度较Euler公式有明显的提
42、高。 xey ixEuler公式改进的Euler公式精确解iyiy iiyxy iiyxy ixy01110.10.20.310.90000000.81000000.72900000.34867843108374. 43107308. 821018182. 121092010. 10.90500000.81902500.74121760.36854100.90483740.81873080.74081820.3678794410626. 1410942. 2410994. 3410616. 6End 第十章常微分方程数值解法EulerEuler 上面介绍的改进公式,其预测公式(公式)的精度差,
43、与校正公式(梯形公式)不匹配。 下面下面再再看两步看两步Euler公式(公式(2.8),除了给出初值外,还),除了给出初值外,还需要借助其它单步法(如需要借助其它单步法(如Euler公式,后退公式,后退Euler公式及梯形公式及梯形公式等)再提供一个启动值公式等)再提供一个启动值y1, 然后才能启动计算公式依次计然后才能启动计算公式依次计算算y2, y3 , , ),(211iiiiyxhfyy ( 2.8 ) 用两步用两步Euler公式与梯形公式相匹配,又可得到下面预测公式与梯形公式相匹配,又可得到下面预测-校正系统:校正系统: 第十章常微分方程数值解法 用两步用两步Euler公式与梯形公式
44、相匹配,又可得到下面预测公式与梯形公式相匹配,又可得到下面预测-校正系统:校正系统: 校正:校正:),(),(2111iiiiiiyxfyxfhyyiiiiyxhfyy,211(2.18) 预测:预测: 与改进的与改进的Euler公式相比较易见(公式相比较易见(2.18)-(2.19)的一个)的一个突出特点是它的预测公式与校正公式具有同阶精度。据此可突出特点是它的预测公式与校正公式具有同阶精度。据此可以比较方便的估计截断误差,并基于这种估计,可以提供一以比较方便的估计截断误差,并基于这种估计,可以提供一种提高精度的简易方法。种提高精度的简易方法。 (2.19)第十章常微分方程数值解法311()
45、( )(2.20)3iiihy xyyx311()( )(2.21)12iiihy xyyx 若(若(2.18)预测公式中的)预测公式中的yi和和yi-1都是准确的。即都是准确的。即 yi=y(xi) , yi-1=y(xi-1)则由两步则由两步Euler公式的截断误差公式公式的截断误差公式(2.12)知:知:再由梯形公式截断误差公式再由梯形公式截断误差公式(2.11)知:知: (2.19)校正公式的截)校正公式的截断误差为:断误差为:校正:校正:),(),(2111iiiiiiyxfyxfhyyiiiiyxhfyy,211 预测:预测:(2.18)(2.19)第十章常微分方程数值解法11)(
46、iiyxy41)()(1111iiiiyxyyxy 比较(比较(2.20)()(2.21)可见,校正值的误差)可见,校正值的误差y(xi+1) -yi+1 大约大约只是预测值只是预测值 的误差的的误差的1/4(符号相反符号相反).即即 由此导出误差的事后估计由此导出误差的事后估计:111111114()()5(2.22)1()()5iiiiiiiiy xyyyy xyyy 第十章常微分方程数值解法 1预测: 2改进: 3计算: *112iiihyyp)(5411iiiicppm),(11*1iiimxfm可以期望利用这样估计出的误差作为计算结果的一种 补偿,有可能使精度得到改善. 以pi和 c
47、i分别表示第i步的预测值和校正值,按估计式(2.22), 和 分别作为pi+1和 ci+1的改进值,在校正值ci+1尚未求出之前,可用上一步的偏差pi- ci替代 pi+1- ci+1 来改进预测值pi+1.这样设计的计算方案有如下六个环节:)(54111iiicpp)(51111iiicpc第十章常微分方程数值解法4校正: 5改进:6计算:)(2*11iiiiymhyc111151iiiicpcy),(11*1iiiyxfy运用上述方案计算 时,要用到先一步的信息 , , 和 更前一步的 。因此启动算法之步必须给出启动值 和 。 可用其它单步法(例如改进的Euler方法)来计算, 则一般取为0。计算结果表明,这种简单的处理方法通常可以获得令人满意的效果。1iy1iy1y11cp 1y11cp *iiiiyypc