常微分方程数值解法ppt课件.ppt

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

1、第七章第七章 常微分方程初值问题数值解法常微分方程初值问题数值解法数值分析数值分析18:09:1518:09:1518:09:152本章内容本章内容18:09:1618:09:1618:09:163简单的数值方法与基本概念 科学技术中常常需要求解常微分方程的定解问科学技术中常常需要求解常微分方程的定解问题题. 这类问题最简单的形式，是本章将着重考察的这类问题最简单的形式，是本章将着重考察的一一阶方程的初值问题阶方程的初值问题00( , ),().yf x yy xy 我们知道，只有我们知道，只有f(x, y)适当光滑适当光滑譬如关于譬如关于y满足满足利普希茨利普希茨(Lipschitz)条件条

2、件理论上就可以保证初值问题的解理论上就可以保证初值问题的解yf(x)存在并且唯一存在并且唯一.),(),(yyLyxfyxf 18:09:1618:09:1618:09:164 虽然求解常微分方程有各种各样的解析方法，但虽然求解常微分方程有各种各样的解析方法，但解析方法只能用来求解一些特殊类型的方程，实际问解析方法只能用来求解一些特殊类型的方程，实际问题中归结出来的微分方程主要靠数值解法题中归结出来的微分方程主要靠数值解法. 所谓所谓数值解法数值解法, 就是寻求解就是寻求解y(x)在一系列离散节点在一系列离散节点 121nnxxxx上的近似值上的近似值 y1,y2,yn,yn+1,. 相邻两个

3、节点的间距相邻两个节点的间距hn=xn+1- -xn称为称为步长步长. 今后如不特别说明，总是假定今后如不特别说明，总是假定 hi=h(i=1,2,)为为定数定数, 这时节点为这时节点为xn=x0+nh(i=0,1,2,) (等距节点等距节点).18:09:1618:09:1618:09:165 初值问题的初值问题的数值解法数值解法有个有个基本特点基本特点，他们都采，他们都采取取“步进式步进式”，即求解过程顺着节点排列的次序一步，即求解过程顺着节点排列的次序一步一步地向前推进一步地向前推进. 描述这类算法，只要给出用已知信描述这类算法，只要给出用已知信息息yn,yn- -1,yn- -2,计算

4、计算yn+1的的递推公式递推公式. 首先，要对微分方程离散化，建立求解数值解的首先，要对微分方程离散化，建立求解数值解的递推公式递推公式. 一类是计算一类是计算yn+1时时只用到前一点的值只用到前一点的值yn，称称为为单步法单步法. 另一类是用到另一类是用到yn+1前面前面 k 点的值点的值yn,yn-1, yn-k+1，称为称为k步法步法. 其次，要研究公式的其次，要研究公式的局部截断误差局部截断误差和和阶阶，数值解，数值解yn与与精确解精确解y(xn)的的误差估计误差估计及及收敛性收敛性，还有递推公式的还有递推公式的计算稳定性计算稳定性等问题等问题.18:09:1718:09:1718:0

5、9:176一阶常微分方程一阶常微分方程 00( , )()yf x yy xy 的解的解( )yy x是通过点是通过点00(,)xy的一条曲线的一条曲线( )yy x，称之为微分，称之为微分方程的积分曲线。积分曲线上每一点方程的积分曲线。积分曲线上每一点( , )x y的切线斜率的切线斜率( )y x等于等于函数函数( , )f x y在这点的值。在这点的值。 1 1 几何推导几何推导 从初始点从初始点000(,)P xy出发，做切线出发，做切线0()y x，与，与1xx交于交于111( ,)P x y点，用点，用1y作为曲线作为曲线( )y x上的点上的点11( , ()x y x的纵的纵坐

6、标坐标1()y x的近似值。 再从的近似值。 再从1P做切线做切线1()y x， 与， 与2xx交于交于222(,)P xy点，用点，用2y作为曲线作为曲线( )y x上的点上的点22(, ()xy x的纵坐标的纵坐标2()y x的近似的近似值。这样下去便可作出一条折线值。这样下去便可作出一条折线012P PP 。设已作。设已作出折线的顶点出折线的顶点为为nP，再从，再从nP做切线做切线()ny x，推进到，推进到111(,)nnnPxy。 欧拉法18:09:1718:09:1718:09:177过过00(,)xy做以做以000()(,)y xf xy为切线斜率的方程为切线斜率的方程 0000

7、(,)()yyf xyxx 当当1xx时，得时，得100010(,)()yyf xyxx，取，取11()y xy。 过过11(,)x y做以做以111( )( ,)y xf x y为切线斜率的方程为切线斜率的方程 1111( ,)()yyf x yxx 当当2xx时，得时，得211121( ,)()yyf x yxx，取，取22()y xy。 一般地，过一般地，过(,)nnxy做以做以()(,)nnny xf xy为切线斜率的方程为切线斜率的方程 (,)()nnnnyyf xyxx 当当1nxx时，得时，得11(,)()nnnnnnyyf xyxx，取，取1()nny xy。 从从0 x出发逐

8、个算出出发逐个算出12,nx xx，对应的数值解，对应的数值解12,nyyy。 一般取一般取1nnxxh，得欧拉公式，得欧拉公式 1(,)nnnnyyhf xy 18:09:1718:09:1718:09:17818:09:1718:09:1718:09:1792 2 欧拉法的数学推导欧拉法的数学推导: :泰勒展开法泰勒展开法 将将1()ny x在在nx处做泰勒展开处做泰勒展开 21()()()()()2!nnnnnhy xy xhy xhy xy 当当h充分小时，忽略高次项得充分小时，忽略高次项得 22()()2!nhyO h 因此，有欧拉公式因此，有欧拉公式 1(,)nnnnyyhf xy

9、 18:09:1718:09:1718:09:17102 2 欧拉法数学推导欧拉法数学推导 ：数值微分（用差商代替导数：数值微分（用差商代替导数 ）nxxxx,210 设设 等距，步长等距，步长1,0,1,nnhxxn yxfhxyxyhxyhxyhxyhxyxy,)()()( 令令x=xn , x+h=xn+1 , y(xn)yn ，y(xn+1 ) yn+1 ,初值问题离散化为初值问题离散化为 100(,) ,0,1,2,()nnnnyyh f xyny xy(欧拉公式欧拉公式) 18:09:1718:09:1718:09:17112 2 欧拉法数学推导欧拉法数学推导: : 数值积分数值积

10、分 将方程将方程( )( , )y xf x y两端从两端从nx到到1nx积分，有积分，有 11( )d( , ( )dnnnnxxxxy xxf x y xx 11()()( , ( )dnnxnnxy xy xf x y xx 算出积分项，可得算出积分项，可得1()ny x。利用左矩形公式。利用左矩形公式 1( , ( )d( , ( )nnxxf x y xxhf x y x 代入，并离散化，有欧拉公式代入，并离散化，有欧拉公式 1(,)nnnnyyhf xy 例例7-1 用欧拉公式求解初值问题用欧拉公式求解初值问题2(01),(0)1.xyyxyy 解解 取步长取步长h=0.1，欧拉公

11、式的具体形式为欧拉公式的具体形式为)2(1nnnnnyxyhyy 其中其中xn=nh=0.1n (n=0,1,10), 已知已知y0 =1, 由此式可得由此式可得191818. 1)1 . 12 . 01 . 1 ( 1 . 01 . 1)2(1 . 11 . 01)2(1111200001 yxyhyyyxyhyy18:09:1718:09:1718:09:1713依次计算下去，依次计算下去，部分计算结果部分计算结果见下表见下表. xy21 与准确解与准确解 相比，可看出欧拉公式的计算结相比，可看出欧拉公式的计算结果精度很差果精度很差. xn 欧拉公式数值解yn准确解y(xn) 误差 0.2

12、 0.4 0.6 0.8 1.0 1.191818 1.358213 1.508966 1.649783 1.784770 1.183216 1.341641 1.483240 1.612452 1.732051 0.008602 0.016572 0.025726 0.037331 0.05271918:09:1718:09:1718:09:1714 局部截断误差和阶：局部截断误差和阶： 数值公式的精度数值公式的精度定义定义 局部截断误差：假设第局部截断误差：假设第n n步是准确的，即步是准确的，即y(xn )=yn，将，将y(xn+1 ) - yn+1定义为数值方法的局部截断定义为数值方法

13、的局部截断误差误差。 由于实际上由于实际上y yn n不是准确值，因此它的误差会传播下去。不是准确值，因此它的误差会传播下去。实际计算时，每一步都可能产生舍入误差。实际计算时，每一步都可能产生舍入误差。定义定义 若局部截断误差为若局部截断误差为O(hO(hp+1p+1), p), p为正整数，则称数为正整数，则称数值公式是值公式是p p阶公式。阶公式。 18:09:1718:09:1718:09:1715 欧拉公式的截断误差是欧拉公式的截断误差是O(h2)，公式是公式是1 阶的。阶的。1(,)()()nnnnnnyyh f xyy xh y x211()() () ( )2nnny xy xy

14、xhyh二阶泰勒公式二阶泰勒公式 两式相减，由设两式相减，由设 yn=y(xn ) ，有有 22112nnhy xyyO h欧拉公式的局部截断误差和阶欧拉公式的局部截断误差和阶18:09:1718:09:1718:09:1716隐式（后退）欧拉公式：取隐式（后退）欧拉公式：取1()ny x的向后差商的向后差商 111() ()()nnny xy xy xh 替代替代111()(,)nnny xf xy中的导数项，并离散化，有隐式欧拉公中的导数项，并离散化，有隐式欧拉公式式 111(,)nnnnyyhf xy 隐式（后退）欧拉公式的局部截断误差：假设隐式（后退）欧拉公式的局部截断误差：假设()n

15、nyy x，则则 211()()2nnnhy xyyx 隐式（后退）欧拉公式隐式（后退）欧拉公式18:09:1718:09:1718:09:1717为了提高精度，改用中心差商为了提高精度，改用中心差商 111 ()()2nny xy xh 替代替代()(,)nnny xf xy中的导数项，并离散化，有两步欧拉公式中的导数项，并离散化，有两步欧拉公式 112(,)nnnnyyhf xy 两步欧拉公式是两步法，要用前两步的值。两步欧拉公式是两步法，要用前两步的值。 两步欧拉公式的局部截断误差两步欧拉公式的局部截断误差 231()()()()()( )2!3!nnnnnhhy xy xhy xhy

16、xy xy 231()()()()()( )2!3!nnnnnhhy xy xhy xhy xy xy 上二式相减，可得上二式相减，可得. . 两步欧拉公式两步欧拉公式18:09:1718:09:1718:09:1718311()()2()( )3nnnhy xy xhy xy 311()()2()( )3nnnhy xy xhy xy 设设()nnyy x，11()nnyy x前两步准确。则上式成为前两步准确。则上式成为 311()2(,)( )3nnnnhy xyhf xyy 与与 112(,)nnnnyyhf xy 相比较，因此，有相比较，因此，有 311()( )3nnhy xyy 两

17、步欧拉公式的局部截断误差是两步欧拉公式的局部截断误差是3()O h，是二阶方法。，是二阶方法。 18:09:1718:09:1718:09:1719梯形法梯形法对微分方程对微分方程y=f(x,y) 两边求两边求x xn n到到x xn+1 n+1 的定积分，有的定积分，有11()()( , ( )dnnxnnxy xy xf x y xx利用梯形公式计算积分，有利用梯形公式计算积分，有 1111( , ( )d (, ()(, ()2nnxnnnnnnxxxf x y xxf xy xf xy x将将y(xn ) 、y(xn+1 )分别用分别用yn、yn+1 代替，构造数值公式代替，构造数值公

18、式11100 (,)(,) ,0,1,2,2()nnnnnnhyyf xyf xynyy x18:09:1718:09:1718:09:1720改进的欧拉法改进的欧拉法( (预报预报- -校正公式校正公式) ) 欧拉法欧拉法 ，显式，计算量小，精度低。，显式，计算量小，精度低。梯形公式梯形公式 是隐式公式是隐式公式 ,计计算量大，精度高。算量大，精度高。 实际计算时，将二者综合之，先用欧拉公式计算出实际计算时，将二者综合之，先用欧拉公式计算出y yn+1n+1作作为初始值为初始值, ,初始值精度不高，取作预报值，代入梯形公式，初始值精度不高，取作预报值，代入梯形公式，得到校正值得到校正值y y

19、n+1n+1。写成预报写成预报- -校正公式校正公式 1111(,) (,)(,)2nnnnnnnnnnyyh f xyhyyf xyf xy1(,)nnnnyyh f xy111 (,)(,)2nnnnnnhyyf xyf xy18:09:1718:09:1718:09:1721预报预报- -校正公式校正公式又常常写成一步嵌套显式形式又常常写成一步嵌套显式形式 或写成平均化形式或写成平均化形式11 (,)(,(,)2nnnnnnnnhyyf xyf xyh f xy11(,)(,)1()2pnnncnnpnpcyyh f xyyyh f xyyyy预报预报- -校正公式的局部截断误差校正公式

20、的局部截断误差y(xn+1)- yn+1=O(h3)1111(,) (,)(,)2nnnnnnnnnnyyh f xyhyyf xyf xy18:09:1718:09:1718:09:172218:09:1818:09:1818:09:1823例 7-2 用改进欧拉法解初值问题2(0)1xyyyy 其中0,1x。 解解 改进欧拉法公式112()2()1()2npnnnncnppnpcxyyh yyxyyh yyyyy 0,1x， 0100022 0()10.1(1)1.11xyyh yy 取步长h=0.1，有n=0 18:09:1818:09:1818:09:1824龙格龙格- -库塔库塔（R

21、unge-KuttaRunge-Kutta) )法法 基本思想基本思想 二阶龙格二阶龙格- -库塔法库塔法 经典龙格经典龙格- -库塔法库塔法18:09:1818:09:1818:09:1825基本思想基本思想: : 一阶常微分方程初值问题一阶常微分方程初值问题00( , )()yf x yy xy 欧拉公式欧拉公式1(,)nnnnyyhf xy，写成，写成111(,)nnnnyyhkkf xy，调用，调用 1 1 次函数。次函数。有有 1 1 阶精度。阶精度。 改进欧拉法，预报改进欧拉法，预报- -校正公式校正公式n+111n+1(,) (,)(,)2nnnnnnnnyyh f xyhyyf

22、 xyf xy 写成写成112121()2(,),(,)nnnnnnhyykkkf xykf xh yhk调用调用 2 2 次函数。 有次函数。 有 2 2 阶精度。阶精度。 龙格-库塔（Runge-Kutta)法18:09:1818:09:1818:09:1826将改进欧拉法（预报将改进欧拉法（预报- -校正公式）改写成平均化的形式校正公式）改写成平均化的形式 1121211()2(,),(,)nnnnnnhyykkkf xykf xyhk 用用 nx与与 1nx两个点的斜率值两个点的斜率值 1k与与 2k取算术平均作为平取算术平均作为平均斜率均斜率*121()2kkk，而，而 1nx处的斜

23、率值处的斜率值 2k 则通过已知信息则通过已知信息 ny 来预测。来预测。 利用利用( , )f x y在某些点的值的线性组合构造公式，得到在某些点的值的线性组合构造公式，得到1()ny x的近似值的近似值1ny， 增加调用， 增加调用( , )f x y的次数， 可提高精度的阶的次数， 可提高精度的阶数。数。 设法在设法在1,nnxx这一步内多预报几个点的斜率，然后将其加这一步内多预报几个点的斜率，然后将其加权平均，作为平权平均，作为平均斜率，则可构造出具有更高精度的计算公式。均斜率，则可构造出具有更高精度的计算公式。 18:09:1818:09:1818:09:1827二阶龙格二阶龙格-库

24、塔法库塔法 考察考察1,nnxx内一点内一点npnxxph，其中，其中01p 用用,nnpxx点的斜率点的斜率12,k k加权平均，得平均斜率，即加权平均，得平均斜率，即 112(1)nnyyhkk 式中式中为待定系数。取为待定系数。取1(,)nnkf xy，2(,)npnpkf xy，先用欧拉法，先用欧拉法提供提供()npy x的预报值的预报值npy，有，有1npnyyphk，因此，因此 21(,)(,)npnpnnkf xyf xph yphk 这样有计算公式这样有计算公式 18:09:1818:09:1818:09:1828112121(1)(,),(,)nnnnnnyyhkkkf xy

25、kf xph yphk 其中含有两个待定参数，选取合适的值，使公式有高的精度。仍假定其中含有两个待定参数，选取合适的值，使公式有高的精度。仍假定()nnyy x，则有，则有1(,)()nnnkf xyy x，21(,)nnkf xph yphk 2(,)(,)(,)(,)()nnxnnnnynnf xyph fxyf xyfxyO h 2()()()nny xphy xO h 这时，有这时，有231()()()()nnnnyy xhy xph y xO h 和和1()ny x的泰的泰勒展开式勒展开式2311()()()()()2nnnny xy xhy xh y xO h比较，从而有构造二阶龙

26、格比较，从而有构造二阶龙格- -库塔法所需条件库塔法所需条件 12p 18:09:1818:09:1818:09:1829此此为为不不定定方方程程，有有无无穷穷多多个个解解。 （1）取取11,2p，得得 112121()2(,)(,)nnnnnnhyykkkf xykf xh yhk 这这就就是是改改进进欧欧拉拉公公式式。 （2）取取1,12p，得得 12121(,)(,)22nnnnnnyyhkkf xyhhkf xyk 该该公公式式称称为为中中点点公公式式。 18:09:1818:09:1818:09:1830经典（四阶）龙格经典（四阶）龙格- -库塔法库塔法 仿照上述的讨论，可导出四阶经

27、典龙格库塔公式仿照上述的讨论，可导出四阶经典龙格库塔公式112341213243(22)6(,),(,)22(,),(,)22nnnnnnnnnnhyykkkkhhkf xykf xykhhkf xykkf xh yhk经典经典龙格库塔公式龙格库塔公式每步要四次计算函数值，具有每步要四次计算函数值，具有四阶精度，局部截断误差是四阶精度，局部截断误差是O(h5) . 18:09:1818:09:1818:09:183118:09:1818:09:1818:09:1832 然而值得指出的是，龙格然而值得指出的是，龙格- -库塔方法的推导基于库塔方法的推导基于泰勒展开泰勒展开方法，因而它要求所求的方

28、法，因而它要求所求的解具有较好的光解具有较好的光滑性质滑性质. 反之，如果解的光滑性差，那么，使用龙反之，如果解的光滑性差，那么，使用龙格格- -库塔方法求得的数值解，其精度可能反而不如改库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法进的欧拉方法. 实际计算时，我们应当实际计算时，我们应当针对问题的具针对问题的具体特点选择合适的算法体特点选择合适的算法.18:09:1818:09:1818:09:183318:09:1818:09:1818:09:1834线性多步法线性多步法kiinikiiniknfhyy010 (n=0,1, )其中其中, xn+i= x0+(n+i)h, fn+i

29、 = f(xn+i , yn+i)局部载断误差局部载断误差kiinikiiniknknxyhxyxyT010)()()( Adamas显格式显格式: yn+2=yn+1+h(3fn+1- fn)/2 yn+3=yn+2+h(23fn+2- 16fn+1 +5fn)/1218:09:1818:09:1818:09:1835y = f (x, y) 21)(,()()(12nnxxnndxxyxfxyxy 2121)(,()(nnnnxxxxdxxyxfdxxy在在区间区间xn , xn+1上插值上插值f(x)(xn+1x)fn + (xxn)fn+1 / h 21)()(111nnxxnnnndxfxxfxxh2/ 323211122nnnnffhfhfhh 二阶二阶Adamas显格式显格式: yn+2=yn+1+h(3fn+1- fn)/2