《计算方法》PPT课件.ppt

第六章 常微分方程初值问题的数值解法,许多科学技术问题，例如天文学中的星体运动，空间技术中的物体飞行，自动控制中的系统分析，力学中的振动，工程问题中的电路分析等，都可归结为常微分方程的初值问题。,所谓初值问题，是函数及其必要的导数在区间的起始点为已知的一类问题，一般形式为：,我们先介绍简单的一阶问题：,常微分方程 求解求什么？应求一满足初值问题的解函数y = y(x) ，如对下列微分方程：,由常微分方程的理论可知：上述问题的解唯一存在。,高等数学中，微分方程求解，如对一阶微分方程：y =f(x,y)是求解解函数y = y(x) ，使满足上述方程。但能够求出,所谓常微分方程初值问题的数值解法，就是求它的解 ，在一系列节点 上的近似值 ，即：,称为步长，一般总取 为常数。,准确的解析函数y(x)的微分方程是很少的，高数中研究微分方程的求解，是分门别类讨论，对不同类型的微分方程，求解方法不一样，因此，要求解微分方程，首先必须认清类型。,由于在实际问题和科学研究中遇到的微分方程往往很复杂，绝大多数很难，甚至不可能求出解析函数y(x)，因此只能考虑求其数值解。,本章重点讨论一阶微分方程，在此基础上介绍一阶微分方程组与高阶微分方程的数值解法。,6.1 欧拉（Euler）法,以Euler法及其改进方法为例,说明常微分方程初值问题数值解法的一般概念，Euler法很简单，准确度也不高，介绍此方法的目的，是由于对它的分析讨论能够比较清楚地显示出方法的一些特点，而这些特点及基本方法反映了其它方法的特点。,Euler法用于求解一阶微分方程初值问题：,6.1.1 欧拉折线法,从几何意义出发：一阶常微分方程： 及初始条件： 的解解函数y=y (x) 在xoy平面上是一族解曲线， 而初值问题则是其中一条积分曲线，假定y = y(x)的曲线如下图所示，从给定的点P0(x0, y0) 出发，以P0为切点，作切线，切线斜率为曲线y(x)的切线斜率y =f (x0,y0)，因此可得切线：,切线与直线 相交：,取 作为 的近似值。,再过点 ，作斜率为 的直线方程：,它与直线 相交：,取 作为 的近似值。,它与直线 相交：,依次下去，就可以求得点 ，通过这点 作斜率为 的直线方程：,由此便得到原初值问题的解 在节点 处的近似值 。,我们把这种求解初值问题的方法，称为欧拉折线法。,可以看出，欧拉折线法的计算公式为：,通过上述计算过程可以看出，欧拉折线法的基本思想是用一系列直线所组成的折线去近似地代替曲线，并用这些直线交点处的纵坐标 作为精确解 的近似值。,欧拉折线法的截断误差：,假设第 n 步求得的 是精确的，即 ，则在第n+1 步，把 称为欧拉折线法的截断误差。,又,由泰勒展开式可知：,所以,这说明欧拉折线法的精度是很差的。,6.1.2 欧拉方法的改进,要求常微分方程 及初始条件： 的解，可以通过积分方法求得：,. 令 ，代入上面积分方程，得：,. 令 ，代入原积分方程，得：,在微分方程两边，从 到 对 求积分：,要求一个函数的数值积分 ，通过第五章的学习，我们已学会了很多方法：梯形公式、辛甫生公式、柯特斯公式以及龙贝格公式。而选用不同的数值积分方法，便导出不同的计算公式：,. 根据递推，一般有：,要求得 的近似值，只要用数值方法求出 的积分近似值就可以了。,.采用矩形公式:以左矩形面积代替曲边梯形面积,如下图所示，亦即以：,.采用矩形公式:以右矩形面积代替曲边梯形面积,如下图所示，亦即以：,.采用梯形公式:以梯形面积代替曲边梯形面积,如下图所示，亦即以：,当 时，由此便得到微分方程初值问题的一系列近似值 ，利用上述公式求解一阶微分方程初值问题的方法称为梯形法则。,代入递推公式，得：,若用 和 分别近似 和 ，则有计算公式：,梯形公式看作是以（xn,yn）(xn+1,yn+1)构造的插值多项式代替被积函数得到的，而Euler公式则是以左端点函数值近似被积函数而得到，还可以用多个点做插值多项式近似被积函数构造另一些精度更高的求解微分方程的数值公式，梯形公式比Euler公式更准确一些，误差更小。,注1 ：Euler公式为显式，右矩形，梯形公式为隐式；,注2：前面已有Euler法的局部截断误差：,显然， Euler 法为1阶方法, 且步长h越小，阶数P 越高，局部截断误差越小，当然计算精度越高；,注3 ：梯形法是几阶？梯形法精度比Euler法高，阶数肯定比Euler法高，其实我们可以利用数值积分公式的误差估计式，因为我们是用梯形数值积分公式计算,因此,由积分中梯形公式的误差知此时的局部截断误差为：,梯形法为2阶方法 ！,注4 ： Euler法与梯形法在计算上有根本区别， Euler公式中的yn为已知值或已算出的值，由yn可直接求出yn+1，这种方法通常称为显式方法，而在梯形法则中， yn+1隐含在函数f(xn+1 , yn+1)中，必须通过解方程才能求出来，因此这种方法称为隐式方法。,. 改进的欧拉方法,梯形法则的精确度要比欧拉公式的精确度要高，但却是隐式格式，而欧拉公式的精确度虽然较低，但却是显式格式。,在数值计算中，将Euler法与梯形公式合起来使用，即先使用Euler公式,由（xn,yn）算出yn+1，记为yn+1(0)，称为预测值，然后用梯形公式去提高精度，将yn+1(0) 校正为较准确的值：,在梯形法则中， ， 隐含在 中，必须通过解方程才能求出来，而在欧拉公式 中 中，已知 ，就可以直接求出 。,当 时，终止迭代，从而求得 ，把 作为 的近似值，这种方法称为改进的欧拉方法。,. 预报-校正公式,改进的欧拉方法，计算量太大，若步长 较小，一般把迭代一次的结果取作 ，于是得到：,上式就称为预报-校正公式，其中第一式为预报公式，第二式为校正公式。为上机计算编程方便，常将上式改写为:,. 预报-校正公式的截断误差,假设 ，由于：,将 和 代入公式 ，并整理：,由此可见，预报-校正公式比欧拉公式的截断误差高了一阶，是二阶方法。,所以预报-校正公式的截断误差为：,由于数值积分的梯形公式的截断误差公式为：,所以梯形公式的截断误差为：,所以梯形法则的截断误差亦为：,所以预报-校正公式的截断误差与梯形法则是同阶的，它们的计算精度是一样的，而预报-校正公式是显式的，便于计算。,到目前为止，对于常微分方程初值问题的数值解法，我们已经学了下面几种方法：,.欧拉折线法,截断误差为：,.梯形法则,截断误差为：,.改进的欧拉方法,预报-校正公式的截断误差都为：,解：用欧拉公式计算，计算公式如下：,. 预报-校正公式,用预报-校正公式计算，计算公式如下：,0,0,0,0,1,0,0.09606,0.09620,2,0.19211,0.37186,0.37288,例,这些结果在下表中，可见计算结果的精度，Euler法与准确值相比较Euler法结果偏小；梯形法则的结果精度较高。,6.2 龙格库塔法(Runge-kutta ),欧拉方法和预报-校正公式虽然算法简单，但精度比较低，不能满足工程计算的要求，因此，在工程上一般不用，必须寻求精度较高的其它方法。,6.2.1 龙格-库塔方法的基本思想：,因此, 只要对平均斜率k*提供一种算法，由公式 y(xn+1) = y(xn) + h f (xn +h, y(xn +h)便相应地得到一种微分方程的数值计算公式。,用这个观点来研究欧拉公式与预报-校正公式，可以发现欧拉公式由于仅取xn一个点的斜率值f (xn,yn)作为平均斜率 k* 的近似值，因此精度很低。而预报-校正公式却是利用了 xn与xn+1两个点的斜率值k1 = f (xn,yn)与k2=f (xn+1,yn+hk1)取算术平均作为平均斜率k*的近似值。,其中k2是通过已知信息yn利用欧拉公式求得的。,预报-校正公式比欧拉公式精度高的原因，也就在于确定平均斜率时，多取了一个点的斜率值。因此它启发我们，如果设法在xi,xi+1上多预报几个点的斜率值，然后将它们加权平均作为k*的近似值，则有可能构造出更高精度的计算公式，这是龙格-库塔方法的基本思想。,6.2.2 二阶龙格-库塔公式,首先推广预报-校正公式，考察区间xn,xn+1内任一点：,我们希望用xn和x n+1两个点的斜率值k1和k2加权平均作为 平均斜率k*的近似值 ：,然后再用预测值yn+l通过计算f 产生斜率值k2=f (xn+l ,yn+l)，这样设计出的计算公式具有形式 ：,其中1，2 为待定常数，同改进欧拉公式一样，这里仍取：,问题在于怎样预测xn+l处的斜率值k2。 仿照预报-校正公式，先用欧拉公式提供y (xn+l)的预测值 :,现在仍假定yn=y(xn)，即yn是准确的，将y(xn+1)与 yn+1都在 xi处作泰勒展开：,公式中含有三个待定参数1, 2和l，我们希望适当选取这些参数值，使得上述公式具有二阶精度，亦即使：,代入原公式，得：,比较上式与泰勒展开式，要使公式具有二阶精度，只有满足下列条件：,这里一共有三个待定参数，但只需满足两个条件，因此有一个自由度，于是满足上述条件的参数不止一组，而是一族，相应的公式求解yn+1的公式也有一族，这些公式统称为二阶龙格-库塔公式，简称二阶R-K公式 。,特别，当l=1即xn+l=xn+1时， 1= 2=1/2，二阶R-K公式就是预报-校正公式。 如果取l=1/2，则1=0, 2=1，这时二阶R-K公式称为变形的欧拉公式，其形式：,从表面上看，变形的欧拉公式仅含一个斜率值k2，但k2是通过k1计算出来的，因此每完成一步，仍然需要两次计算函数 f 的值，工作量和预报-校正公式相同。,综上所述，构造二阶R-K公式主要由以下几步产生：,在区间xn,xn+1上取二点，预报相应点的斜率值； 2) 对此两斜率值加权平均作为平均斜率值的近似值； 3) 将yn+1与y(xn+1)都在xi处作泰勒展开，为使公式达到二阶精度，比较相应系数，建立有关参数所应满足的方程组； 4) 解此方程组得一族二阶R-K公式。,基于二阶龙格-库塔公式设计思想，仿照预报-校正公式，我们可构造如下形式的计算公式：,6.2.3 高阶龙格-库塔公式,其中 ， ， 都是待定常数。上述公式称为 m 阶显式龙格-库塔公式。,.当m = 2时二阶的龙格-库塔公式。,.若取m= 4时，类似于二阶龙格-库塔法的推导过程，经过繁锁的推导，便得到实际中最常用的四阶龙格-库塔公式：,一般称它为标准四阶R-K公式，其截断误差为：,.将经典的四阶龙格-库塔公式加以改正，便得到基尔公式。,解：用四阶的龙格-库塔公式，计算公式如下：,. 龙格-库塔法存在的问题,利用龙格-库塔法每求解一次 ,需要反复计算4次 的值，计算量比较大。,在计算过程中，为达到精度要求，一般要减少步长 ，而步长 越小，截断误差也越小，但在同一区间 内，计算 的点数就越多，这样，不仅增大了计算量，而且还有可能导致舍入误差的严重积累。,且截断误差不容易估计。,龙格-库塔方法的推导是在用泰勒展开方法的基础上进行的，因而它要求所求的解具有较好的光滑性质。假若解的光滑性差，那么使用四阶R-K公式求得的数值解，其精度可能反而不如预报-校正公式。在实际计算时，应当针对问题的具体特点，选择合适的算法。,6.2.4 变步长的龙格-库塔法,当用数值方法解微分方程时，单从一步来看，步长越小，局部截断误差就越小。但从整体来看，步长越小，在求解范围内所要完成的步数就会越多。这一方面增加了计算量；另一方面也增大了舍入误差的积累，因此有个如何合理选择步长的问题。在解决这个问题时，需要考虑两个问题： 怎样衡量和检验计算结果的精度？ 如何依据获得的精度信息及时处理步长。,以四阶标准R-K公式为例，从节点xn出发，先以步长h求得xn+1处解的近似值，记作yn+1h，由于公式的局部截断误差为O(h5)，故有：,当 不大时， 可以看成是常数，然后将步长折半，即取 为步长，从 出发，经过计算求得 的近似值为 ，由于每跨一步的局部截断误差为c(h/2)5，因此有：,两式相除有：,整理后得：,上式说明以 作为 的近似值，其误差可用步长分半前后两次计算结果之差来表示。,3. 以h作基本步长，当 时，反复将步长折半，每折半一 次，由xn+1的步数增加一倍，直至,2.如果 ，又不需要节点xn+h处的数值解值，则可将步长加倍，只要保证所求节点处的y值符合精度要求即可。,这种方法，我们称为变步长的四阶-龙格库塔法。,例.,解 取基本步长h=0.1，节点xk=kh， =106，按变步长四阶标准R-K法进行计算， 其结果见下表，显然，尽管每一步的步长均只折半一次，但解的精度可达106，与准确解：,比较，说明以上分析是正确的。,6.3 一阶微分方程组的初值问题,下面以两个方程的情形为例，介绍一阶方程组的数值解法。设初值问题：,若采用向量记号，记：,则初值问题可表示为：,它与常微分方程初值问题形式上完全一致，前面介绍的解初值问题的各种数值方法均可用于解一阶方程组初值问题。,下面就以四阶的龙格-库塔法为例，求解上述一阶微分方程组一阶方程组初值问题。,设有一阶微分方程组的初值问题：,下面同样以四阶的龙格-库塔法为例，求解上述一阶微分方程组。,则初值问题可