chap常微分方程的数值解法PPT课件

1、1第第7章章 常微分方程（组）的常微分方程（组）的数值解法数值解法刘东毅刘东毅天津大学理学院数学系天津大学理学院数学系第1页/共56页2第7章 常微分方程（组）的数值解法主要目标： 掌握常微分方程初值问题数值解法的基本理论掌握计算机上的常用算法主要内容：初值问题计算格式的建立Runge-Kutta 方法一阶常微分方程组与高阶方程的数值解法第2页/共56页3第7章 常微分方程（组）的数值解法 在科学研究和工程实践中会遇到很多微分方程，虽然从理论上可以证明其解的存在性，但其解的解析表达式往往是很难求解的，或者即使可以写出来，但也难于计算，此时，只能借助数值解来解决问题. 常微分方程（组）定解问题是

2、自然科学和工程技术领域中常见的数学模型. 本章介绍求解此类问题的基本理论和数值解法。 第3页/共56页4yR定义定义7.0.1 若存在常数若存在常数 L 0, 使得对一切的使得对一切的 xa , b 及及 y , , 均有均有 |( , )( , )|,f x yf x yL yy则称则称 f (x, y) 在在 D 上关于上关于 y 满足满足 Lipschitz 条件条件, 其中其中 L 称为称为 Lipschitz 常数常数. 我们首先考虑一阶常微分方程我们首先考虑一阶常微分方程初值问题初值问题0( , ),( ).yf x yaxby ay 其中其中 f (x , y) 是区域是区域(

3、, )|,Dx yaxbyR上的实值函数上的实值函数.(7.0.1)第4页/共56页5 我们首先给出常微分方程初值问题解的存在惟一性定理。定理 假设 f (x, y)C(D), ,且关于 y 满足 Lipschitz 条件, 则一阶常微分方程初值问题(7.0.1) 存在唯一解. 下面在此前提下, 我们讨论上述初值问题 (7.0.1) 的数值解法。0( , ),( ).yf x yaxby ay 第5页/共56页6012,Naxxxxb然后在节点上建立逼近于原初值问题的计算格式然后在节点上建立逼近于原初值问题的计算格式 (或差分格式或差分格式), 由此计算出由此计算出原问题原问题的解的解 y(

4、x ) 在节点在节点 x1 , x2 , . . . , xN 处的近似值：处的近似值： y1 , y2 , . . ., yN， 称它们为称它们为常微分方程初常微分方程初值问题的数值解值问题的数值解. 相邻两个节点的距离相邻两个节点的距离 hn= xn+1 - xn 称为称为步长步长, 通常取定步长通常取定步长h 0, 即即节点节点 xn = x0 + nh , n = 0,1, , N. 其基本思想是在区间 a , b 上引入一系列节点第6页/共56页7 初值问题计算格式的建立1. 数值微分方法数值微分方法在等距节点下讨论问题在等距节点下讨论问题. 利用两点数值微分公式利用两点数值微分公式

5、.1 计算格式的建立计算格式的建立将上式代入初值问题将上式代入初值问题 (7.0.1),)()(0yayx,yfy1()()(, ()(),2nnnnny xy xhf xy xyh有有11()()()(),()2nnnnnnny xy xhy xyxxh(7.1.1) 第7页/共56页8略去余项, 并以数值解 yn , yn+1 代替 y (xn) 及 y (xn+1), 则得差分方程差分方程上式称为上式称为 Euler 公式公式. 利用此式可由利用此式可由初值初值 y0 出发按出发按“步进式步进式” 方法方法, 逐步逐步求得数值解求得数值解y1 , y2 , . . . , yN .)2

6、. 1 . 7(),(1nnnnyxhfyy由于计算由于计算yn+1时时, 只只用到它前一步的结用到它前一步的结果果yn , 这类公式称为这类公式称为单步法单步法. 又因为其关又因为其关于于yn+1是显式形式是显式形式, 故称该故称该Euler公式为公式为显格式显格式.第8页/共56页9如果利用下列数值微分公式111()()()(),()2nnnnnnny xy xhy xyxxh由由 类似的可导出类似的可导出)()(xx,yfxy),(111nnnnyxhfyy上述公式称为上述公式称为后退的后退的 Euler 公式公式, 此公式为单步法公式此公式为单步法公式. 又因为它关于又因为它关于 yn

7、+1 成成隐式形式隐式形式, 所以该公式为所以该公式为隐式公式隐式公式，简称，简称隐格式隐格式.第9页/共56页102111()()()(),()26nnnnnnny xy xhy xyxxh类似地，可导出类似地，可导出上述公式称为上述公式称为Euler两步法两步法公式公式. 这因为，当计算这因为，当计算 yn+1 时时, 要用到要用到 yn -1 与与 yn . 显然它显然它也是显格式也是显格式.如果利用下列三点数值微分公式112).(,nnnnyyhf xy(7.1.3) 第10页/共56页11设设 y (x)C2a , b, 由由 Taylor 公式有公式有由于由于 故上式即为故上式即为

8、 ( )( , ( ) ,y xf x y x略去余项略去余项, 并以并以 yn , yn+1 代替代替 y (xn) 及及 y (xn+1), 得到的差分方程正是得到的差分方程正是Euler 公式公式.211()()()()().2nnnnnnnhy xy xhy xyxx(7.1.4) )5 . 1 . 7()(! 2)(,()()(21nnnnnyhxyxhfxyxy 2. Taylor 展开展开法法第11页/共56页123. 数值积分方法数值积分方法对对 ，在区间，在区间 xn , xn+1 上积分，得上积分，得( )( , )y xf x y11( )( , ( ),nnnnxxxx

9、y x dxf x y x dx则有则有11()()( , ( ).nnxnnxy xy xf x y x dx对上式中的积分采用不同的数值积分公式可得到不同的差分方程对上式中的积分采用不同的数值积分公式可得到不同的差分方程. 例如例如, 对上式对上式的积分采用左矩形公式的积分采用左矩形公式, 可得到可得到 Euler 公式公式. 第12页/共56页1311131()() (, ()(, ()2(, ()().12nnnnnnnnnnnhy xy xf xy xf xy xhfyxx若对此式的积分采用梯形公式, 11()()( , ( ).nnxnnxy xy xf x y x dx则有则有若

10、略去余项若略去余项, 以以 yn , yn+1 代替代替 y (xn) 及及 y (xn+1), 得到的差分方程得到的差分方程111 (,)(,2.)nnnnnnhyyf xyf xy第13页/共56页14上式称为上式称为梯形公式梯形公式. 由于它关于由于它关于 yn+1 成隐式形式成隐式形式, 故其为故其为隐格式隐格式. 隐格式求解比较困难隐格式求解比较困难, 当当 yn 已知时已知时, 要求要求 yn+1 ，需解关于需解关于 yn+1 的非线性方程的非线性方程. 在实际应用时在实际应用时, 上式常与上式常与 Euler 公式联合使用公式联合使用, 构成如下计算格式构成如下计算格式:111

11、(,)(,2.)nnnnnnhyyf xyf xy(7.1.6)(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk第14页/共56页15隐式梯形公式的迭代格式(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk(7.1.7)由由上式可以得到一个序列上式可以得到一个序列: , k = 0,1, , 关于此序列的收敛性关于此序列的收敛性, 有如下的定理有如下的定理.( )1kny第15页/共56页16 设 f (x , y) 在区域 D 上关于 y

12、满足 Lipschitz 条件, 即|( , )( , )|()|.f x yf x yLyy其中其中 L 为为 Lipschitz 常数常数, 当步长当步长 时时, 对任意的初值对任意的初值 按格式（按格式（7.1.7）2hL(0)1ny(0)1(1)( )111(,), (,)(,)(0,1,)2nnnnkknnnnnnyyh f xyhyyf xyf xyk生成的序列生成的序列 收敛于梯形公式（收敛于梯形公式（7.1.6）( )1kny111 (,)(,)2nnnnnnhyyf xyf xy1ny的解的解 .第16页/共56页17为了减少计算量为了减少计算量, 可采用可采用预测预测- -

13、校正格式校正格式. 方法是先用方法是先用Euler公式求得一个初公式求得一个初始近似值始近似值 称为称为预测值预测值, 再把再把 带入梯形公式右端计算一次求得带入梯形公式右端计算一次求得 yn+1 称之为称之为校正校正值值, 即即1ny1ny预测预测:校正校正:1(,),nnnnyyhf xy111 (,)(,).2nnnnnnhyyf xyf xy上式称为上式称为预测预测 - - 校正公式校正公式或或改进的改进的 EulerEuler 公式公式. 上式也可写成如下形式上式也可写成如下形式:11 (,)(,(,)2.nnnnnnnnhyyf xyf xyhf xy第17页/共56页18例：利用

14、Euler公式与改进的Euler公式求解初值问题（步长）. 1)0(, 10,2yxyxyy解：由步长解：由步长h=0.1，知节点知节点 设数值解为设数值解为 利用利用Euler公式得公式得,1 . 00nnhxxn.10, 2 , 1 , 0n,nynnnnnnnnnnnyxyyxyyyxhfyy2 . 01 . 121 . 0),(1第18页/共56页19计算结果见下表(见书P227表7.1) 此初值问题的解析解为此初值问题的解析解为 , 从上表可以看出从上表可以看出, 数值解数值解 yn与解析与解析解解 y(xn) 比较比较, yn精度较差精度较差. xy21第19页/共56页20解此问

15、题的改进的Euler公式为 ),(1nnnnyxhfyy),(),(2111nnnnnnyxfyxfhyy.1 . 011. 0105. 12221 . 0,2 . 01 . 121 . 01111111nnnnnnnnnnnnnnnnnnnnnyxyxyyxyyxyyyyxyyxyyy同同Euler公式比较公式比较, 改进的改进的Euler法显然精度提高了法显然精度提高了.由于误差大小是评价计算格式优劣的重要依据由于误差大小是评价计算格式优劣的重要依据, 故需要给出有关误差的概念故需要给出有关误差的概念. 计算结计算结果见下果见下表表(见书见书P228表表7.2)第20页/共56页21 截断

16、误差与方法的精度定义定义 7 称误差称误差 en+1 = y ( xn+1 ) - - yn+1为为数值方法在点数值方法在点 xn+1 的截断误差的截断误差, 又称又称整体截断误差整体截断误差. 设设 yk= y ( xk ) (k = 0,1,. . . , n)，则则 为数值方法在点为数值方法在点 xn+1 的的局部截断误差局部截断误差.111()nnny xy第21页/共56页22整体截断误差 en+1 是在没有引进舍入误差的情况下, 纯粹因为不准确的计算格式造成的, 故又称为方法误差.它不仅与 x = xn+1 这一步的计算有关, 而且和 xn , xn-1 ,. . . , x1 这

17、几步的计算都有关系. 局部截断误差是假设局部截断误差是假设 xn 之前各数值解没有误差之前各数值解没有误差, 仅由仅由 xn 到到 xn+1 这一步这一步计算由计算计算由计算格式引起的误差格式引起的误差.第22页/共56页23如Euler公式1(,)nnnnyyhf xy在点在点 xn+1 的整体截断误差的整体截断误差 en+1 = y (xn+1)- - yn+122()(, ()()(,)2(, ()(,)()2nnnnnnnnnnnnnhy xhf xy xyyhf xyheh f xy xf xyy局部截断误差局部截断误差2111()()2nnnnhy xyy第23页/共56页24定义

18、定义7.1.2 若某数值方法的局部截若某数值方法的局部截断断误差为误差为 则称该方法具有则称该方法具有 P 阶精度阶精度, 或称其为或称其为 P 阶方法阶方法.11(),pnO h可以证明可以证明:Euler 方法的局部截断误差方法的局部截断误差 其具有其具有一阶一阶精度精度. 梯形方法的局部截断误差梯形方法的局部截断误差 其具有其具有二阶二阶精度精度.改进的改进的 Euler 方法方法的局部截断误差的局部截断误差 具有具有二阶二阶精度精度.21(),nO h31(),nO h31(),nO h第24页/共56页257.2 Runge-Kutta 方法 继续讨论前面的继续讨论前面的 Taylo

19、r 展开展开法。法。设设 y (x)C2a , b, 由由 Taylor 公式有公式有211()()()()().2nnnnnnnhy xy xhy xyxx 由由 故上式即为故上式即为 ( )( , ( ) ,y xf x y x)(! 2)(,()()(21nnnnnyhxyxhfxyxy 略去余项略去余项, 并以并以 yn , yn+1 代替代替 y (xn) 及及 y (xn+1), 得到得到Euler 公式公式.第25页/共56页26进一步假设设 y (x)Cp+1a , b,由 Taylor 公式有) 1 . 2 . 7(,)(!)(! 2)()()()(21nnppnnnnRxy

20、phxyhxyhxyxy )2 . 2 . 7.(),()()!1(11)1(1nnnpnppnxxhOyphR其中其中 由由 故式故式(7.2.1)即为即为 ( )( , ( ) ,y xf x y x第26页/共56页277.2 Runge-Kutta 方法略去余项略去余项, 并以数值解并以数值解 yn , yn+1 代替代替 (7.2.3) 中的中的 解析解解析解y (xn) 及及 y (xn+1), 可得到一个可得到一个差分方程，即差分方程，即nnppnnnnRxyphxyhxyhxyxy )(!)(! 2)()()()(21) 3 . 2 . 7(,)(,(!)(,(! 2)(,()

21、()()1(21nnnppnnnnnnRxyxfphxyxfhxyxhfxyxy., )()(,()!1(11)(1nnnpnnppnxxhOyfphR其中余项可写成其中余项可写成)(,(),()()(xyxfdxdyxfkkk注：这里注：这里第27页/共56页28在(7.2.3)中略去余项,用yn , yn+1 代替 y (xn) 及 y (xn+1)4 . 2 . 7(.! 2! 2)1(1)1(21pnpnnnpnpnnnnfphfhfhyfphfhfhyy, ),(nnnyxff . ) 1, 2 , 1(),()()(piyxffnniin其中其中称称 (7.2.4) 式为求解常微分

22、方程初值问题数值解式为求解常微分方程初值问题数值解Taylor的格式的格式 .21(1)()()(,()(,()2!(,(),(7.2.3)!nnnnnnppnnnhy xy xhf xy xfxy xhfxy xRp第28页/共56页29由于局部截断误差 ),(O)(1111pnnnhyxy可知它是一个可知它是一个 p 阶方法。当阶方法。当p=1时，上式正是时，上式正是Euler 公式。但当公式。但当 p 2 时，需要计算时，需要计算f (x, y(x) ) 的高阶导数，特别是对于复杂函数的高阶导数，特别是对于复杂函数 f (x, y(x) 的求导，这无疑是大大增加计的求导，这无疑是大大增加

23、计算量，这是它最大的缺点。因此高阶的算量，这是它最大的缺点。因此高阶的Taylor方法是不实用的。方法是不实用的。 德国数学家德国数学家C.Runge 及提出了一种改进策略，得到了至今还被作为高精度的单及提出了一种改进策略，得到了至今还被作为高精度的单步法广泛使用龙格步法广泛使用龙格-库塔法库塔法 ( Runge- Kutta method )。 第29页/共56页307.2.1 Runge-Kutta方法的基本思想.! 2)1(11pnpnnnnfphfhfhyyRunge-Kutta方法是方法是利用利用 f 在在某些点处函数值的线性组合替代（某些点处函数值的线性组合替代（7.2.4）步长）

24、步长 h 后面括号中的因子来构造差分方程后面括号中的因子来构造差分方程, 从而避免了高阶导数的计算从而避免了高阶导数的计算, 这就是这就是 Runge-Kutta 方法的基本思想方法的基本思想. )1(1! 2pnpnnfphfhf用用f 在在某些点处函数值的线性组合替代这一部分某些点处函数值的线性组合替代这一部分第30页/共56页31其一般形式为:11111,(,),(,)(2,3, )rnniiinniininijjjyyhkkf xykf xh yhkir其中其中 r 是上式中调用是上式中调用 f 的个数的个数, r 称为级数，称为级数， 为待定参数为待定参数, 适当确定这些参数适当确定

25、这些参数, 可使之具有尽可能高的精度可使之具有尽可能高的精度. 如局部截断如局部截断误差满足误差满足,iii j . )(O)(1111rnnnhyxy第31页/共56页32 二阶 Runge-Kutta 方法考虑考虑 r = 2 的情况的情况, 此时有此时有11 122121(),(,),(,).nnnnnnyyhkkkf xykf xh yh k利用二元函数的利用二元函数的 一阶一阶Taylor 公式，即全微分公式公式，即全微分公式 希望适当选择参数希望适当选择参数 使上式的局部截断误差为使上式的局部截断误差为12, , 3111()()nnny xyO h即为二阶方法即为二阶方法.,)(

26、)()(,()(,(),(),(22nnnnnynnnxnnyyxxOyyyxfxxyxfyxfyxf下面将下面将yn+1与与y(xn+1)作比较作比较第32页/共56页332222( ,)(,)(,)()(,)()()()() ()() ()()()nnxnnnynnnnndefnxnnynnnnf x yf xyfxyxxfxyyyOxxyyffxxfyyOxxyy. )O()()()(212hkfhfhfknynxn从而有, )(22111kkhyynn将上式代入),(12khyhxfknn再由第33页/共56页34得到. )()()()()O()()()()(3222121211221

27、11hOfffhfhyhkfhfhfhkhykkhyynynnxnnnynxnnnn在下面要将yn+1与y(xn+1)作比较，使它们的局部截断误差满足. )O()()()(212hkfhfhfknynxn, )(O)(3111hyxynnn为此考虑y(xn+1)。第34页/共56页35再根据 y(xn+1) 在点 xn的一元 3 阶 Taylor 展开式. )()()(! 2)(),(! 2)()(! 2)()()(323),(2321hOfffhhfyhOyxfhhfyhOxyhxyhxyxynynnxnnyxnnnnnnnn )()()()(322211hOfffhfhyynynnxnnn

28、由刚才已得到的, )(O)(3111hyxynnn让它们满足第35页/共56页36即由122211212左式含有四个未知元三个左式含有四个未知元三个方程方程, 因此因此解不唯一解不唯一. 参数参数满足左式的满足左式的一族公式一族公式统称统称二阶二阶 Runge-Kutta 公式公式. 可得参数应满足下列方程组:. )()()(! 2)(321hOfffhhfyxynynnxnnn)()()()(322211hOfffhfhyynynnxnnn第36页/共56页37 取上式称为上式称为中点公式中点公式 .1210,1,212,nnyyhk1(,),nnkf xy21(,).22nnhhkf xy

29、k 取上式称为上式称为 Heun 公式公式 .112(3),4nnhyykk1(,),nnkf xy2122(,).33nnkf xh yhk12132,443 可见可见, 二阶二阶 Runge-Kutta 公式公式, 每计算一步需要每计算一步需要 两次调用两次调用 f 的函数值的函数值. 取121,1,2112121(),2(,),(,).nnnnnnhyykkkf xykf xh yhk得这正是这正是改进的改进的 Euler 公式公式. 第37页/共56页38 四阶Runge-Kutta方法当当 r = 4 时时, 类似地可导出四阶类似地可导出四阶 Runge-Kutta 公式公式, 这种

30、公式也有一族这种公式也有一族, 其中常用地有其中常用地有:标准标准 (经典经典) 的的 Runge-Kutta 方法方法11234(22),6nnhyykkkk1(,),nnkf xy21(,),22nnhhkf xyk32(,),22nnhhkf xyk43(,).nnkf xh yk第38页/共56页3911234(22)(22),6nnhyykkkk1(,),nnkf xy21(,),22nnhhkf xyk3122122(,),222nnhkf xyhkhk423222(,).22nnkf xh yhkhk Gill 公式公式 Gill 公式是标准的公式是标准的 Runge-Kutta

31、 公式的改进形式公式的改进形式, 这种算法可节省存储单元这种算法可节省存储单元, 并能控并能控制舍入误差的增长制舍入误差的增长.四阶四阶 Runge-Kutta 公式公式, 每一步计算需四次调用每一步计算需四次调用 f 的函数值的函数值, 计算量较大计算量较大, 但其局部截断但其局部截断误差可达误差可达 O(h5), 精度较高精度较高.第39页/共56页40例例7.2.1 用标准的用标准的 四阶四阶Rung-Kutta 法解初值问题法解初值问题,取取步长步长h=0.2.2, (0)1,(01).xyyyxy 解解: 解此问题的计算公式为解此问题的计算公式为112340.2(22),6nnyyk

32、kkk12,nnnxkyy21122,22()nnnhxhkykhyk32222,22()nnnhxhkykhyk4332,2()nnnxhhkykyhkxn0.20.40.60.81.0ynyn- y (xn)1.183 21.341 71.483 3 1.612 5 1.732 10.000 00.000 00.000 0 0.000 10.000 1计算结果如下计算结果如下:显然在计算量大致相同的情显然在计算量大致相同的情况下况下, 标准的标准的 Runge-Kutta方方法比改进的法比改进的 Euler 方法精确方法精确度更高度更高. （参见（参见p227和和p228的表）的表）第40

33、页/共56页41一阶常微分方程组与高阶方程初值问题的数值解法 7.5.1 一阶常微分方程组初值问题 )y,y,y(x,fdxdy)y,y,y(x,fdxdy)y,y,y(x,fdxdym21mmm2122m2111( )( )( )1122mmy a = s ,ya = s , ya = s ,xa,b第41页/共56页42写成向量形式： 0( ,)( ) (7.5.2)dYF x YdxY a =Y,)()()(Yn21n21xyxyxyyyy.s ,s ,sYTn210,)()()()(21212211mmmm,y,yx,yf,y,yx,yf,y,yx,yfx,YF其中其中注意注意: 在形

34、式上在形式上 (7.5.2) 与与 (7.0.1) 一样一样, 所以可以把求解常微分方程初值问题的各所以可以把求解常微分方程初值问题的各种数值方法推广到方程组上来种数值方法推广到方程组上来. 第42页/共56页43利用向量值函数的微积分理论, 很容易推导出一阶常微分方程组初值问题的数值解法. )3 . 5 . 7(),(1nnnnYxhFYY如如Euler公式公式,),(T,2, 1nmnnnyyyY其中其中 .),(),(),(),(, 2, 1, 2, 12, 2, 11nmnnnmnmnnnnmnnnnnyyyxfyyyxfyyyxfYxF第43页/共56页44(7.5.3)的分量形式为

35、 )4 . 5 . 7(. ), 2, 1(),(, 2, 1,1,miyyyxhfyynmnnninini),(),(),(, 2, 1, 2, 12, 2, 11, 2, 11,1, 21, 1nmnnnmnmnnnnmnnnnmnnnmnnyyyxfyyyxfyyyxfhyyyyyy或或 第44页/共56页45四阶标准的Runge-Kutta公式 )5 . 5 . 7()22(643211KKKKhYYnn),(1nnYxFK )2,2(12KhYhxFKnn)2,2(23KhYhxFKnn),(34hKYhxFKnn设设 则则 (7.5.5)的分量形式为的分量形式为 ,)4, 3, 2

36、, 1(),(T,2, 1jkkkKjmjjj第45页/共56页46四阶标准的Runge-Kutta公式的分量形式 , )2,2,2,2(),2,2,2,2(),2,2,2,2(),(),22(63 ,3 , 23 , 14,2,2, 22, 13 ,1 ,1 , 21 , 12, 2, 11 ,4,3 ,2,1 ,1,inmininniiinmininniiinmininniinmnnniiiiiininikhykhykhyhxfkkhykhykhyhxfkkhykhykhyhxfkyyyxfkkkkkhyy., 2, 1mi其中其中第46页/共56页47例：试写出用中点公式解下列初值问题的

37、计算公式. 3)0(, 1)0(,35,643zyzyzzyxy解解:令令 ,zyY则则 .nnnzyY.35643),(zyzyxYxF,)2 , 1(, 2, 1jkkKjjj再取再取由向量形式的中点公式由向量形式的中点公式第47页/共56页48中点公式的向量形式 )2,2(),(12121KhYhxFKYxFKhKYYnnnnnn上述中点公式的分量计算形式为上述中点公式的分量计算形式为2, 22, 111kkhzyzynnnnnnnnnzyzyxkk356431 , 21 , 1)2( 3)2( 5)2(6)2(4)2( 31 , 21 , 11 , 21 , 12, 22, 1khzkhykhzkhyhxkknnnnn分量计算形式为分量