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

常微分方程初值问题数值解法第一页，讲稿共八十一页哦 在高等数学中，对于常微分方程的求解，给出了一在高等数学中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等。数齐次线性方程的解法、常系数非齐次线性方程的解法等。但能求解的常微分方程仍然是有限的，大多数的常微分但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解。方程是不可能给出解析解。譬如譬如 这个一阶微分方程就不能用初等函数及其积分来表达这个一阶微分方程就不能用初等函数及其积分来表达它的解。它的解。第二页，讲稿共八十一页哦再如，方程再如，方程 的解的解 ,虽然有表可查虽然有表可查,但对于表但对于表上没有给出上没有给出 的值的值,仍需插值方法来计仍需插值方法来计算算第三页，讲稿共八十一页哦从实际问题当中归纳出来的微分方程，通常主要依靠数从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决。本章主要讨论一阶常微分方程初值问题值解法来解决。本章主要讨论一阶常微分方程初值问题(9.1)在区间在区间a x ba x b上的数值解法上的数值解法。可以证明可以证明,如果函数在带形区域如果函数在带形区域 R=axb,R=axb,-y y内连续，且关于内连续，且关于y y满足李普希兹满足李普希兹(Lipschitz)(Lipschitz)条件，即存在常数条件，即存在常数L(L(它与它与x,yx,y无关无关)使使对对R内任意两个内任意两个 都成立都成立,则方程则方程(9.1)的解的解 在在 a,b 上存在且唯一。上存在且唯一。第四页，讲稿共八十一页哦数值方法的基本思想数值方法的基本思想 对常微分方程初值问题对常微分方程初值问题(9.1)式的数值解法，就是要式的数值解法，就是要算出精确解算出精确解y(x)y(x)在区间在区间 a,b 上的一系列离散节点上的一系列离散节点 处的函数值处的函数值 的近似值的近似值。相邻两个节点的间距。相邻两个节点的间距 称为步长，步称为步长，步长可以相等，也可以不等。本章总是假定长可以相等，也可以不等。本章总是假定h为定数，称为定数，称为为定步长定步长，这时节点可表示为，这时节点可表示为数值解法需要把连续性的问题加以离散化，从而求出离散数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解。节点的数值解。第五页，讲稿共八十一页哦 对常微分方程数值解法的对常微分方程数值解法的基本出发点就是离散化。基本出发点就是离散化。其数值解法有两个基本特点，它们都采用其数值解法有两个基本特点，它们都采用“步进式步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息描述这类算法，要求给出用已知信息 计算计算 的递推公式的递推公式。建立这类递。建立这类递推公式的基本方法是在这些节点上用推公式的基本方法是在这些节点上用数值积分、数值微分、数值积分、数值微分、泰勒展开等离散化方法泰勒展开等离散化方法，对初值问题，对初值问题中的导数中的导数 进行不同的离散化处理进行不同的离散化处理。第六页，讲稿共八十一页哦对于初值问题对于初值问题的数值解法，首先要解决的问题就是的数值解法，首先要解决的问题就是如何如何对微分方程进行离对微分方程进行离散化，建立求数值解的递推公式。递推公式通常有两类，一散化，建立求数值解的递推公式。递推公式通常有两类，一类是计算类是计算yi+1时只用到时只用到xi+1,xi和和yi，即前一步的值，因此，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为有了初值以后就可以逐步往下计算，此类方法称为单单步法步法；其代表是；其代表是龙格龙格库塔库塔法。另一类是计算法。另一类是计算y yi+1i+1时，除时，除用到用到x xi+1i+1,x,xi i和和y yi i以外，还要用到以外，还要用到，即前面，即前面k步的值，此类方法称为步的值，此类方法称为多步法多步法；其代表；其代表是亚当斯法。是亚当斯法。第七页，讲稿共八十一页哦9.2简单简单的数的数值值方法与基本概念方法与基本概念9.2.1Euler公式公式 欧拉（欧拉（Euler）方法是解初值问题的最简）方法是解初值问题的最简单的数值方法。初值问题单的数值方法。初值问题的解的解y=y(x)y=y(x)代表通过点代表通过点 的一条称之为微的一条称之为微分方程的积分曲线。积分曲线上每一点分方程的积分曲线。积分曲线上每一点 的切线的斜率的切线的斜率 等于函数等于函数 在这点的在这点的值。值。第八页，讲稿共八十一页哦Euler法的求解过程是法的求解过程是:从初始点从初始点P0(即点即点(x(x0 0,y,y0 0)出发出发,作积分曲线作积分曲线y=y(x)y=y(x)在在P0点上切点上切线线 (其斜率为其斜率为 ),),与与x=xx=x1 1直线直线相交于相交于P1点点(即点即点(x(x1 1,y,y1 1),),得到得到y y1 1作为作为y(xy(x1 1)的近似值的近似值,如如上图所示。过点上图所示。过点(x(x0 0,y,y0 0),),以以f(xf(x0 0,y,y0 0)为为斜率的切线方程为斜率的切线方程为 当当 时时,得得 这样就获得了这样就获得了P P1 1点的坐标。点的坐标。第九页，讲稿共八十一页哦同样同样,过过点点P1(x x1 1,y,y1 1),),作积分曲线作积分曲线y=y(x)y=y(x)的切线的切线交直线交直线x=xx=x2 2于于P2点点,切线切线 的斜率的斜率 =直线方程为直线方程为当当 时时,得得 第十页，讲稿共八十一页哦当当 时时,得得由此获得了由此获得了P P2 2的坐标。重复以上过程的坐标。重复以上过程,就可获得一系列的点就可获得一系列的点:P P1 1,P P1 1,P Pn n。对已求得点对已求得点以以 =为斜率作直线为斜率作直线 取取第十一页，讲稿共八十一页哦 从图形上看从图形上看,就获得了一条近似于曲线就获得了一条近似于曲线y=y(x)y=y(x)的折线的折线 。这样这样,从从x x0 0逐个算出逐个算出对应的数值解对应的数值解 第十二页，讲稿共八十一页哦通常取通常取 (常数常数),),则则Euler法的计算格式法的计算格式 i=0,1,n(9.2)还可用还可用数值微分数值微分、数值积分法数值积分法和和泰勒展开法泰勒展开法推导推导EulerEuler格式。格式。以数值积分为例进行推导。以数值积分为例进行推导。将方程将方程 的两端在区间的两端在区间 上积分得，上积分得，选择不同的计算方法计算上式的积分项选择不同的计算方法计算上式的积分项 ,就会得到不同的计算公式。就会得到不同的计算公式。(9.3)第十三页，讲稿共八十一页哦 用左矩形方法计算积分项用左矩形方法计算积分项 代入代入(9.3)(9.3)式式,并用并用y yi i近似代替式中近似代替式中y(xy(xi i)即可得到向前欧即可得到向前欧拉（拉（EulerEuler）公式）公式 由于数值积分的矩形方法精度很低，所以欧由于数值积分的矩形方法精度很低，所以欧拉（拉（EulerEuler）公式当然很粗糙。）公式当然很粗糙。第十四页，讲稿共八十一页哦例例9.1用欧拉法解初值问题用欧拉法解初值问题取步长取步长h=0.2,h=0.2,计算过程保留计算过程保留4 4位小数位小数 解解:h=0.2,:h=0.2,欧拉迭代格式欧拉迭代格式 当当k=0,x1=0.2时，已知时，已知x0=0，y0=1，有，有y(0.2)y1=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(4-0.40.6144)=0.4613第十五页，讲稿共八十一页哦clear;y=1,x=0,%初始化for n=1:10 y=1.1*y-0.2*x/y,x=x+0.1,endy=1 x=0 y=1.1000 x=0.1000y=1.1918 x=0.2000y=1.2774 x=0.3000y=1.3582 x=0.4000y=1.4351 x=0.5000y=1.5090 x=0.6000y=1.5803 x=0.7000y=1.6498 x=0.8000y=1.7178 x=0.9000y=1.7848 x=1.0000第十六页，讲稿共八十一页哦9.2.2梯形公式梯形公式为了提高精度为了提高精度,对方程对方程 的两端在区间上的两端在区间上 积分得，积分得，改用梯形方法计算其积分项，即改用梯形方法计算其积分项，即 (9.4)代入代入(7.4)(7.4)式式,并用近似代替式中即可得到梯形公式并用近似代替式中即可得到梯形公式 (9.5)由由于于数数值值积积分分的的梯梯形形公公式式比比矩矩形形公公式式的的精精度度高高，因因此此梯梯形公式（形公式（9.59.5）比欧拉公式）比欧拉公式(9.2)(9.2)的精度高一个数值方法。的精度高一个数值方法。第十七页，讲稿共八十一页哦(9.5)(9.5)式的右端含有未知的式的右端含有未知的y yi+1i+1,它是一个关于它是一个关于y yi+1i+1的函数方程的函数方程,这类数值方法称为这类数值方法称为隐式方法隐式方法。相反地。相反地,欧拉欧拉法是关于法是关于y yi+1i+1的一个直接的计算公式，的一个直接的计算公式，这类数值方法这类数值方法称为称为显式方法。显式方法。第十八页，讲稿共八十一页哦9.2.3两步欧拉公式两步欧拉公式对方程对方程 的两端在区间上的两端在区间上 积分得积分得 (9.6)改用中矩形公式计算其积分项，即改用中矩形公式计算其积分项，即 代入上式代入上式,并用并用y yi i近似代替式中近似代替式中y(xy(xi i)即可得到两步欧拉即可得到两步欧拉公式公式 (9.7)第十九页，讲稿共八十一页哦 前面介绍过的数值方法前面介绍过的数值方法,无论是欧拉方法无论是欧拉方法,还是梯形方法，它们都是单步法还是梯形方法，它们都是单步法,其特点是在其特点是在计算计算y yi+1i+1时只用到前一步的信息时只用到前一步的信息y yi i;可是公式可是公式(7.7)中除了中除了y yi i外外,还用到更前一步的信息还用到更前一步的信息y yi-1i-1,即即调用了前两步的信息调用了前两步的信息,故称其为两步欧拉公式故称其为两步欧拉公式第二十页，讲稿共八十一页哦9.2.4欧拉法的局部截断误差欧拉法的局部截断误差衡衡量量求求解解公公式式好好坏坏的的一一个个主主要要标标准准是是求求解解公公式式的的精精度度,因此引入局部截断误差和阶数的概念。因此引入局部截断误差和阶数的概念。定定义义9.1在在yi准准确确的的前前提提下下,即即时时,用用数数值值方方法法计计算算yi+1的的误误差差,称称为为该该数数值值方方法法计计算算时时yi+1的局部截断误差。的局部截断误差。对于欧拉公式，假定对于欧拉公式，假定，则有，则有而将真解而将真解y(x)在在xi处按二阶泰勒展开处按二阶泰勒展开 因此有因此有 第二十一页，讲稿共八十一页哦定定义义9.2数数值值方方法法的的局局部部截截断断误误差差为为,则则称称这这种种数数值值方方法法的的阶阶数数是是P。步步长长(hN结束。结束。第二十六页，讲稿共八十一页哦（2）改改进进欧欧拉拉法法的的流流程程图图 第二十七页，讲稿共八十一页哦(3)3)程序实现程序实现(改进欧拉法计算常微改进欧拉法计算常微 分方程初值问题分方程初值问题)例例9.2 9.2 用改进欧拉法解初值问题用改进欧拉法解初值问题 区间为区间为 0,10,1,取步长取步长h=0.1h=0.1 解解:改进欧拉法的具体形式改进欧拉法的具体形式 本题的精确解为本题的精确解为 ,第二十八页，讲稿共八十一页哦clearx=0,yn=1%初始化for n=1:10yp=yn+0.1*(yn-2*x/yn);%预测x=x+0.1;yc=yn+0.1*(yp-2*x/yp);yn=(yp+yc)/2%校正end第二十九页，讲稿共八十一页哦例例9.3对初值问题对初值问题 证明用梯形公式求得的近似解为证明用梯形公式求得的近似解为并证明当步长并证明当步长h h0 0时时,y,yn n收敛于精确解收敛于精确解证明证明:解初值问题的梯形公式为解初值问题的梯形公式为 整理成显式整理成显式 反复迭代反复迭代,得到得到 第三十页，讲稿共八十一页哦由于由于 ，有，有 证毕证毕 第三十一页，讲稿共八十一页哦9.3 9.3 龙格龙格-库塔（库塔（Runge-KuttaRunge-Kutta）法）法9.3.1 9.3.1 龙格龙格-库塔库塔(Runge-Kutta)(Runge-Kutta)法的基本思想法的基本思想 Euler Euler公式可改写成公式可改写成 则则yi+1i+1的的表表达达式式y(xi+1i+1)与与的的TaylorTaylor展展开开式式的的前前两两项项完完全全相相同同,即局部截断误差为即局部截断误差为 。改进的改进的EulerEuler公式又可改写成公式又可改写成 第三十二页，讲稿共八十一页哦 上述两组公式在形式上有一个共同点上述两组公式在形式上有一个共同点:都是用都是用f(x,y)f(x,y)在在某些点上值的线性组合得出某些点上值的线性组合得出y(xy(xi+1i+1)的近似值的近似值y yi+1i+1,而且增加计而且增加计算的次数算的次数f f(x x,y y)的次数的次数,可提高截断误差的阶。如欧拉公式可提高截断误差的阶。如欧拉公式:每步计算一次每步计算一次f f(x x,y y)的值的值,为一阶方法。改进欧拉公式需为一阶方法。改进欧拉公式需计算两次计算两次f f(x x,y y)的值，它是二阶方法。它的局部截断误的值，它是二阶方法。它的局部截断误差为差为 。第三十三页，讲稿共八十一页哦 于是可考虑用函数于是可考虑用函数f(x,y)在若干点上的函数值在若干点上的函数值的线性组合来构造近似公式，构造时要求近似公式的线性组合来构造近似公式，构造时要求近似公式在在(xi i,yi i)处的处的Taylor展开式与解展开式与解y(x)在在xi i处的处的Taylor展开式的前面几项重合，从而使近似公式展开式的前面几项重合，从而使近似公式达到所需要的阶数。既避免求偏导达到所需要的阶数。既避免求偏导,又提高了计又提高了计算方法精度的阶数。或者说算方法精度的阶数。或者说,在在 这一步内多这一步内多预报几个点的斜率值，然后将预报几个点的斜率值，然后将其加权平均作为平均斜率其加权平均作为平均斜率，则可构造出更高精度，则可构造出更高精度的计算格式，这就是龙格的计算格式，这就是龙格库塔（库塔（Runge-Kutta）法的基本思想。法的基本思想。第三十四页，讲稿共八十一页哦9.3.2二阶龙格二阶龙格库塔法库塔法 在在 上取两点上取两点xi i和和 ,以该两点处的斜以该两点处的斜率值率值k1 1和和k2 2的加权平均的加权平均(或称为线性组合或称为线性组合)来求取平均斜来求取平均斜率率k*的近似值的近似值K，即，即式中式中:k1 1为为xi i点处的切线斜率值，点处的切线斜率值，k2 2为为 点处的切线斜率值点处的切线斜率值,比照改进的欧拉法比照改进的欧拉法,将将 视为视为 ，即可得，即可得 对常微分方程初值问题对常微分方程初值问题(9.1)(9.1)式的解式的解 y=y(x),),根据微分中值定根据微分中值定理，存在点理，存在点 ，使得，使得 第三十五页，讲稿共八十一页哦式中式中 K K可看作是可看作是y=y(x)y=y(x)在区间在区间 上的平均斜率。所以可上的平均斜率。所以可得计算公式为：得计算公式为：（9.14）将将y(xy(xi i)在在x=xx=xi i处进行二阶处进行二阶TaylorTaylor展开：展开：（9 9.15）也即也即 （9.13）第三十六页，讲稿共八十一页哦将将 在在x=xx=xi i处进行一阶处进行一阶TaylorTaylor展开：展开：将以上结果代入（将以上结果代入（9.149.14）得：）得：（9.16）对式对式(9.15)(9.15)和和(9.16)(9.16)进行比较系数后可知进行比较系数后可知,只要只要 （9.17）成立成立,格式格式(9.14)(9.14)的局部截断误差就等于的局部截断误差就等于有有2 2阶阶精度精度第三十七页，讲稿共八十一页哦式式(9.17)(9.17)中具有三个未知量中具有三个未知量,但只有两个方程但只有两个方程,因而有无因而有无穷多解。若取穷多解。若取 ,则则p p=1=1，这是无穷多解中的一个解，这是无穷多解中的一个解，将以上所解的值代入式将以上所解的值代入式(9.14)(9.14)并改写可得并改写可得 不难发现，上面的格式就是改进的欧拉格式。凡不难发现，上面的格式就是改进的欧拉格式。凡满足条件式（满足条件式（9.179.17）有一簇形如上式的计算格式，这）有一簇形如上式的计算格式，这些格式统称为二阶龙格些格式统称为二阶龙格库塔格式。因此改进的欧拉格库塔格式。因此改进的欧拉格式是众多的二阶龙格式是众多的二阶龙格库塔法中的一种特殊格式。库塔法中的一种特殊格式。第三十八页，讲稿共八十一页哦若取若取 ,则则 ，此时二阶龙格，此时二阶龙格-库塔库塔法的计算公式为法的计算公式为 此计算公式称为变形的二阶龙格此计算公式称为变形的二阶龙格库塔法。式中库塔法。式中 为区间为区间 的中点。的中点。第三十九页，讲稿共八十一页哦9.3.3三阶龙格三阶龙格-库塔法库塔法 为了进一步提高精度，设除为了进一步提高精度，设除 外再增加一点外再增加一点 并用三个点并用三个点 ,的斜率的斜率k1 1,k2 2,k3 3加权平均加权平均得出平均斜率得出平均斜率k*的近似值，这时计算格式具有形式的近似值，这时计算格式具有形式:（9.18）为了预报点为了预报点 的斜率值的斜率值k3 3,在区间在区间 内有两内有两个斜率值个斜率值k1 1和和k2 2可以用可以用,可将可将k1 1,k2 2加权平均得出加权平均得出 上的平均斜率上的平均斜率,从而得到从而得到 的预报值的预报值 第四十页，讲稿共八十一页哦于是可得于是可得运用运用Taylor展开方法选择参数展开方法选择参数 ,可以使格式可以使格式(9.18)的局部截断误差为的局部截断误差为 ,即具有三阶精度，这类即具有三阶精度，这类格式统称为格式统称为三阶龙格三阶龙格库塔方法库塔方法。下列是其中的一种，。下列是其中的一种，称为称为库塔（库塔（Kutta）公式。）公式。（9.19）第四十一页，讲稿共八十一页哦第四十二页，讲稿共八十一页哦9.3.4四阶龙格四阶龙格库塔法库塔法如如果果需需要要再再提提高高精精度度，用用类类似似上上述述的的处处理理方方法法，只只需需在在区区间间上上用用四四个个点点处处的的斜斜率率加加权权平平均均作作为为平平均均斜斜率率k*的的近近似似值值，构构成成一一系系列列四四阶阶龙龙格格库库塔塔公公式式。具具有有四四阶阶精度，即局部截断误差是精度，即局部截断误差是。由于推导复杂，这里从略，只介绍最常用的一种由于推导复杂，这里从略，只介绍最常用的一种四阶经四阶经典龙格典龙格库塔公式库塔公式。（9.20）第四十三页，讲稿共八十一页哦9.3.5 9.3.5 四阶龙格四阶龙格库塔法算法实现库塔法算法实现(1)(1)计算步骤计算步骤 输入输入 ,h,Nh,N 使用龙格使用龙格库塔公式（库塔公式（9.209.20）计算出）计算出y y1 1 输出输出 ，并使，并使 转到转到 直至直至n n N N 结束。结束。第四十四页，讲稿共八十一页哦(2 2)四四阶阶龙龙格格库库塔塔算算法法流流程程图图第四十五页，讲稿共八十一页哦(3)(3)程序实现程序实现(四阶龙格四阶龙格-库塔法计库塔法计 算常微分方程初值问题算常微分方程初值问题)例例9.4 9.4 取步长取步长h=0.2h=0.2，用经典格式求解初值问题，用经典格式求解初值问题 解解:由四阶龙格由四阶龙格-库塔公式可得库塔公式可得 第四十六页，讲稿共八十一页哦可同样进行其余可同样进行其余y yi i的计算。本例方程的解为的计算。本例方程的解为，从表中看到所求的数值解具有，从表中看到所求的数值解具有4位有效数字。位有效数字。龙格龙格库塔方法的库塔方法的推导基于推导基于Taylor展开方法，因展开方法，因而它要求所求的解具有较好的光滑性。如果解的光滑而它要求所求的解具有较好的光滑性。如果解的光滑性差，那么，使用四阶龙格性差，那么，使用四阶龙格库塔方法求得的数值解，库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法其精度可能反而不如改进的欧拉方法。在实际计算时，。在实际计算时，应当针对问题的具体特点选择合适的算法。应当针对问题的具体特点选择合适的算法。第四十七页，讲稿共八十一页哦9.3.6变步长的龙格变步长的龙格-库塔法库塔法在在微微分分方方程程的的数数值值解解中中，选选择择适适当当的的步步长长是是非非常常重重要要的的。单单从从每每一一步步看看，步步长长越越小小，截截断断误误差差就就越越小小；但但随随着着步步长长的的缩缩小小，在在一一定定的的求求解解区区间间内内所所要要完完成成的的步步数数就就增增加加了了。这这样样会会引引起起计计算算量量的的增增大大，并并且且会会引引起起舍舍入入误误差差的的大大量量积积累累与与传传播播。因因此此微微分分方方程程数数值值解解法法也有选择步长的问题。也有选择步长的问题。以经典的四阶龙格以经典的四阶龙格-库塔法库塔法(9.20)为例。从节点为例。从节点x xi i出出发，先以发，先以h为步长求出一个近似值，记为为步长求出一个近似值，记为 ，由于局部，由于局部截断误差为截断误差为 ，故有，故有当h h值不大时，式中的系数值不大时，式中的系数c c可近似地看作为常数。可近似地看作为常数。第四十八页，讲稿共八十一页哦然后将步长折半然后将步长折半,即以为即以为 步长步长,从节点从节点x xi i出发出发,跨两跨两步到节点步到节点x xi+1i+1,再求得一个近似值再求得一个近似值 ,每跨一步的截断误差每跨一步的截断误差是是 ,因此有因此有这样这样 由此可得由此可得 这表明以这表明以 作为作为 的近似值，其误差可用步长折的近似值，其误差可用步长折半前后两次计算结果的偏差半前后两次计算结果的偏差 来判断所选步长是否适当来判断所选步长是否适当第四十九页，讲稿共八十一页哦当要求的数值精度为当要求的数值精度为时：时：（1 1）如如果果，反反复复将将步步长长折折半半进进行行计计算算，直直至至为止为止,并取其最后一次步长的计算结果作为并取其最后一次步长的计算结果作为 （2 2）如如果果为为止止，并并以上一次步长的计算结果作为以上一次步长的计算结果作为 。这这种种通通过过步步长长加加倍倍或或折折半半来来处处理理步步长长的的方方法法称称为为变变步步长长法法。表表面面上上看看，为为了了选选择择步步长长，每每一一步步都都要要反反复复判判断断，增增加加了了计计算算工工作作量量，但但在在方方程程的的解解y(x)y(x)变变化化剧剧烈烈的的情况下，总的计算工作量得到减少，结果还是合算的。情况下，总的计算工作量得到减少，结果还是合算的。第五十页，讲稿共八十一页哦9.4算法的稳定性及收敛性算法的稳定性及收敛性9.4.1稳定性稳定性稳稳定定性性在在微微分分方方程程的的数数值值解解法法中中是是一一个个非非常常重重要要的的问问题题。因因为为微微分分方方程程初初值值问问题题的的数数值值方方法法是是用用差差分分格格式式进进行行计计算算的的，而而在在差差分分方方程程的的求求解解过过程程中中，存存在在着着各各种种计计算算误误差差，这这些些计计算算误误差差如如舍舍入入误误差差等等引引起起的的扰扰动动，在在传传播播过过程程中中，可可能能会会大大量量积积累累，对对计计算算结结果果的的准准确确性性将将产产生生影影响响。这就涉及到算法稳定性问题。这就涉及到算法稳定性问题。第五十一页，讲稿共八十一页哦 当当在在某某节节点点上上x xi i的的y yi i值值有有大大小小为为的的扰扰动动时时，如如果果在在其其后后的的各各节节点点 上上的的值值y yi i产产生生的的偏偏差差都都不不大大于于，则称这种方法是稳定的。，则称这种方法是稳定的。稳定性不仅与算法有关，而且与方程中函数稳定性不仅与算法有关，而且与方程中函数f(x,y)f(x,y)也也有关，讨论起来比较复杂。有关，讨论起来比较复杂。为简单起见，通常只针对模型为简单起见，通常只针对模型方程方程来来讨讨论论。一一般般方方程程若若局局部部线线性性化化，也也可可化化为为上上述述形形式式。模模型型方方程程相相对对比比较较简简单单，若若一一个个数数值值方方法法对对模模型型方方程程是是稳稳定定的的，并并不不能能保保证证该该方方法法对对任任何何方方程程都都稳稳定定，但但若若某某方方法法对对模型方程都不稳定，也就很难用于其他方程的求解。模型方程都不稳定，也就很难用于其他方程的求解。第五十二页，讲稿共八十一页哦先考察显式先考察显式EulerEuler方法的稳定性。模型方程方法的稳定性。模型方程的的EulerEuler公式为公式为 将上式反复递推后，可得将上式反复递推后，可得 或或式中式中第五十三页，讲稿共八十一页哦 要使要使y yi i有界，其充要条件为有界，其充要条件为 即即 由于由于 ，故有，故有 （9.269.26）可见，如欲保证算法的稳定，显式可见，如欲保证算法的稳定，显式EulerEuler格式的步长格式的步长h h的选取要受到式（的选取要受到式（9.269.26）的限制。）的限制。的绝对值越大，则的绝对值越大，则限制的限制的h h值就越小。值就越小。用隐式用隐式EulerEuler格式，对模型方程格式，对模型方程 的计算公式为，可化为的计算公式为，可化为第五十四页，讲稿共八十一页哦由于由于 ,则恒有则恒有 ,故恒有故恒有 因此，隐式因此，隐式EulerEuler格式是绝对稳定的（无条件稳定）的格式是绝对稳定的（无条件稳定）的（对任何（对任何h0h0）。）。9.4.2 9.4.2 收敛性收敛性 常常微微分分方方程程初初值值问问题题的的求求解解，是是将将微微分分方方程程转转化化为为差差分分方方程程来来求求解解，并并用用计计算算值值y yi i来来近近似似替替代代y(xy(xi i)，这这种种近近似似替替代代是是否否合合理理，还还须须看看分分割割区区间间 的的长长度度h h越越来来越越小小时时，即即 时时，是否成立。若成立，则称该方法是收敛的，否则称为不收敛。是否成立。若成立，则称该方法是收敛的，否则称为不收敛。这里仍以这里仍以EulerEuler方法为例，来分析其收敛性。方法为例，来分析其收敛性。EulerEuler格式格式第五十五页，讲稿共八十一页哦设设表示取表示取时时,按按Euler公式的计算结果公式的计算结果,即即Euler方法局部截断误差为：方法局部截断误差为：设有常数设有常数 ，则，则 （9.27）总体截断误差总体截断误差 （9.28）又又 由于由于f(x,y)f(x,y)关于关于y y满足李普希茨条件满足李普希茨条件,即即 第五十六页，讲稿共八十一页哦代入代入(9.28)上式，有上式，有再利用式（再利用式（9.27），式（），式（9.28）即即上式反复递推后，可得上式反复递推后，可得（9.29）第五十七页，讲稿共八十一页哦设设 （T T为常数）为常数）因为因为 所以所以 把上式代入式（把上式代入式（9.299.29），得），得 若不计初值误差，即若不计初值误差，即 ，则有，则有 （9.30）式式(9.30)(9.30)说明说明,当时当时 ,即即 ,所以所以EulerEuler方法是收敛的，且其收敛速度为方法是收敛的，且其收敛速度为 ，即具有，即具有一阶收敛速度。同时还说明一阶收敛速度。同时还说明EulerEuler方法的整体截断误差方法的整体截断误差为为 ，因此算法的精度为一阶。，因此算法的精度为一阶。第五十八页，讲稿共八十一页哦9.5亚当姆斯方法亚当姆斯方法9.5.1亚当姆斯格式亚当姆斯格式龙龙格格-库库塔塔方方法法是是一一类类重重要要算算法法，但但这这类类算算法法在在每每一一步步都都需需要要先先预预报报几几个个点点上上的的斜斜率率值值，计计算算量量比比较较大大。考考虑虑到到计计算算yi+1之之前前已已得得出出一一系系列列节节点点上上的的斜斜率率值值，能能否否利利用用这这些些已已知值来减少计算量呢？知值来减少计算量呢？这就是亚当姆斯（这就是亚当姆斯（Adams）方法的设计思想。）方法的设计思想。第五十九页，讲稿共八十一页哦 设用设用x xi i,x,xi+1i+1两点的斜率值加权平均作为区间两点的斜率值加权平均作为区间 上的平均斜率，有计算格式上的平均斜率，有计算格式（9.21）选取参数选取参数，使格式（，使格式（9.219.21）具有二阶精度。）具有二阶精度。第六十页，讲稿共八十一页哦将将 在在x xi i处处Taylor展开展开 代入计算格式代入计算格式(9.21)(9.21)化简化简,并假设并假设,因此有因此有 与与y(xi+1)在在xi处的处的Taylor展开式展开式相比较相比较,需取需取 才使格式才使格式(9.21)具有二阶精度。这样导出的计算格式具有二阶精度。这样导出的计算格式称之为二阶亚当姆斯格式。类似地可以导出三阶亚当姆称之为二阶亚当姆斯格式。类似地可以导出三阶亚当姆斯格式斯格式。第六十一页，讲稿共八十一页哦和和四阶亚当姆斯格式。四阶亚当姆斯格式。(9.22)这里和下面均记这里和下面均记上述几种亚当姆斯格式都是显式的，算法比较简单，上述几种亚当姆斯格式都是显式的，算法比较简单，但用节点但用节点的斜率值来预报区间的斜率值来预报区间上的上的平均斜率是个平均斜率是个外推过程，外推过程，效果不够理想。为了进一步改善效果不够理想。为了进一步改善精度，变外推为精度，变外推为内插内插，即增加节点，即增加节点xi+1的斜率值来得出的斜率值来得出上的平均斜率。譬如考察形如上的平均斜率。譬如考察形如第六十二页，讲稿共八十一页哦（9.23）的隐式格式的隐式格式,设设(9.23)右端的右端的Taylor展开有展开有 可见要使格式可见要使格式(9.23)(9.23)具有二阶精度具有二阶精度,需令需令 ,这样就可构造二阶隐式亚当姆斯格式这样就可构造二阶隐式亚当姆斯格式 其实是梯形格式。类似可导出三阶隐式亚当姆斯格式其实是梯形格式。类似可导出三阶隐式亚当姆斯格式 第六十三页，讲稿共八十一页哦和四阶隐式亚当姆斯格式和四阶隐式亚当姆斯格式（9.249.24）9.5.2亚当姆斯预报亚当姆斯预报-校正格式校正格式参照改进的欧拉格式的构造方法，以四阶亚当姆斯为例，参照改进的欧拉格式的构造方法，以四阶亚当姆斯为例，将显式（将显式（9.22）和隐式（）和隐式（9.24）相结合，用显式公式做预报，）相结合，用显式公式做预报，再用隐式公式做校正，可构成亚当姆斯预报再用隐式公式做校正，可构成亚当姆斯预报-校正格式校正格式（9.259.25）预报：预报：校正：校正：第六十四页，讲稿共八十一页哦 这种预报这种预报-校正格式是四步法，它在计算校正格式是四步法，它在计算y yi+1i+1时不但用到前一步的信息时不但用到前一步的信息 ，而且要用到，而且要用到再前面三步的信息再前面三步的信息 ，因此它不能，因此它不能自行启动。在实际计算时，可借助于某种单步自行启动。在实际计算时，可借助于某种单步法，譬如四阶龙格法，譬如四阶龙格库塔法提供开始值库塔法提供开始值 。第六十五页，讲稿共八十一页哦例例9.5取步长取步长h=0.1,h=0.1,用亚当姆斯预报用亚当姆斯预报-校正公式求解校正公式求解 初值问题初值问题的数值解。的数值解。解解:用四阶龙格用四阶龙格-库塔公式求出发值库塔公式求出发值 ,计算得：计算得：表中的表中的 ,y ,yi i和和y(xy(xi i)分别为预报值、校正值和准确解分别为预报值、校正值和准确解(),(),以比较计算结果的精度。以比较计算结果的精度。再使用亚当姆斯预报再使用亚当姆斯预报-校正公式校正公式(9.25),(9.25),第六十六页，讲稿共八十一页哦9.6一阶常微分方程组与高阶方程一阶常微分方程组与高阶方程我我们们已已介介绍绍了了一一阶阶常常微微分分方方程程初初值值问问题题的的各各种种数数值值解解法法，对对于于一一阶阶常常微微分分方方程程组组，可可类类似似得得到到各各种种解解法法，而而高高阶常微分方程可转化为一阶常微分方程组来求解。阶常微分方程可转化为一阶常微分方程组来求解。9.6.1一阶常微分方程组一阶常微分方程组对于一阶常微分方程组的初值问题对于一阶常微分方程组的初值问题（9.31）可以把单个方程可以把单个方程中的中的f和和y看作向量来处理，看作向量来处理，这样就可把前面介绍的各种差分算法推广到求一阶方程组初这样就可把前面介绍的各种差分算法推广到求一阶方程组初值问题中来。值问题中来。第六十七页，讲稿共八十一页哦设设 为节点上的近似解，为节点上的近似解，则有改进的则有改进的EulerEuler格式为格式为 预报：预报：校正：校正：（9.32）又，相应的四阶龙格又，相应的四阶龙格库塔格式（经典格式）为库塔格式（经典格式）为 （9.33）第六十八页，讲稿共八十一页哦式中式中（9.34）第六十九页，讲稿共八十一页哦把节点把节点xi上的上的yi和和zi值代入式值代入式(7.34),依次算出依次算出,然然后后把把它它们们代代入入式式(7.33),算算出节点出节点xi+1上的上的yi+1和和zi+1值。值。对于具有三个或三个以上方程的方程组的初值问题对于具有三个或三个以上方程的方程组的初值问题,也可用也可用类似方法处理类似方法处理,只是更复杂一些而已。此外只是更复杂一些而已。此外,多步方法也同样可多步方法也同样可以应用于求解方程组初值问题。以应用于求解方程组初值问题。例例7.6 7.6 用改进的用改进的EulerEuler法求解初值问题法求解初值问题 取步长取步长h=0.1h=0.1，保留六位小数。，保留六位小数。解解:改进的改进的EulerEuler法公式为法公式为预报：预报：校正：校正：将将 及及h=0.1h=0.1代入上式代入上式,得得 第七十页，讲稿共八十一页哦由初值由初值 ,计算得计算得第七十一页，讲稿共八十一页哦9.6.2高阶方程组高阶方程组 高阶微分方程高阶微分方程(或方程组或方程组)的初值问题的初值问题,原则上都原则上都可以归结为一阶方程组来求解。例如可以归结为一阶方程组来求解。例如,有二阶微分方有二阶微分方程的初值问题程的初值问题（9.359.35）在引入新的变量在引入新的变量 后后,即化为一阶方程组初值问题即化为一阶方程组初值问题:（9.369.36）式（式（9.369.36）为一个一阶方程组的初值问题，对此可用前面中）为一个一阶方程组的初值问题，对此可用前面中介绍的方法来求解。例如应用四阶龙格介绍的方法来求解。例如应用四阶龙格-库塔公式得库塔公式得 第七十二页，讲稿共八十一页哦（9.379.37）（9.389.38）第七十三页，讲稿共八十一页哦消去消去 ，上式简化为：，上式简化为：（9.399.39）（9.409.40）上述方法同样可以用来处理三阶或更高阶的微分方程上述方法同样可以用来处理三阶或更高阶的微分方程（或方程组）的初值问题（或方程组）的初值问题 第七十四页，讲稿共八十一页哦例例9.7 9.7 求解下列二阶微分方程的初值问题求解下列二阶微分方程的初值问题 取步长取步长h=0.1 h=0.1 解解:先作变换：令先作变换：令 ，代入上式，得一阶方程组，代入上式，得一阶方程组 用四阶龙格用四阶龙格-库塔方法求解库塔方法求解,按式按式(9.37)(9.37)及及(9.38)(9.38)进行计算：进行计算：取步长取步长