《西安石油大学现代数值计算方法第8章.pptx》由会员分享，可在线阅读，更多相关《西安石油大学现代数值计算方法第8章.pptx（70页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、本章着重讨论一阶常微分方程初值问题 的数值解法。8.0 概述 常微分方程初值问题的数值解是求上述初值问题的解y(x)在区间a，b中的点列上的近似值 .以下设不变，记为h-步长。第1页/共70页定理：如果f(x，y)满足李普希兹（Lipschitz）条件则上述微分方程有唯一解y(x)第2页/共70页 假设解y(x)在区间a,b上是存在而且唯一的，并且具有充分的光滑度，因此，要求f(x,y)也充分光滑。初值问题的解析解(理论解）用 表示,数值解法的精确解用 表示。常微分方程数值解法一般分为：（1）一步法：在计算 时，只用到 ,和 ，即前一步的值。（2）多步法：计算 时，除用到 ,和 以外，还要用
2、和 ,即前 k步的值。（3）显式格式与隐式格式。第3页/共70页8.1 欧拉法与梯形法设节点为 ，得欧拉方法计算公式为：一、欧拉(Euler)法下面通过几种常用的方法来推导该公式。1、泰勒展开法假设在 附近把y（x）做Taylor展开，有：第4页/共70页取h的线性部分,并用 表示 的近似值,得2、数值积分法 从 到 +h对等式 y(t)=f(t,y(t)进行积分得到 再利用左矩形公式，得从而得到Euler公式。由第5页/共70页3、数值微分法4、几何方法过点(xn,yn)作以f(xn,yn)为斜率的直线方程：将x=xn+1处该直线上的函数值做为y(xn+1)的近似值，则有Euler公式。这实
3、质上是在每个小区间上利用折线来代替曲线的结果，故Euler法又称Euler折线法。第6页/共70页二、梯形法在式 中,将积分用梯形公式来代替，则有 从而得到梯形公式：梯形方法关于yn+1是隐式的，而Euler方法是显式的。一般情形下不容易从上式解出yn+1，因而可将上式与Euler公式联合使用，即第7页/共70页使用上式时，先用第一式算出xn+1处yn+1的初始近似再用第二式反复迭代，得到数列用来控制迭代次数，这里为允许误差。把满足误差要求的可以证明,当f(x,y)满足Lipschitz条件,即：(L为Lipschitz常数)时,上述数列收敛。作为y(xn+1)的近似值yn+1.类似地可以得出
4、yn+2,yn+3,第8页/共70页证明：由和有：反复使用不等式有：第9页/共70页 实用中,在h 取得较小时,用梯形公式计算,第二式只迭代一次就结束,得到Euler预估-校正格式:第一式称为预估公式，第二式称为校正公式。三、Euler预估-校正格式第10页/共70页四、方法的误差估计、收敛性和稳定性 定义1：为 某一数值方法在xn处的整体截断误差(不考虑舍入误差的影响)。定义2：对单步法，在 的假设下，称为在 处的局部截断误差。（P232定义1）Remark1:Euler法的局部截断误差为(由泰勒余项):Remark2:梯形方法的局部截断误差为（由梯形积分）第11页/共70页用泰勒展开法推导
5、Euler预估校正格式的局部截断误差改写Euler预估校正公式为：在的假定下，第12页/共70页而第13页/共70页而因此有故Euler预估校正方法为的局部截断误差阶为O(h3)。第14页/共70页定义3：若一个方法的局部截断误差为 ，则称该方法为p阶方法,或称该方法具有p阶精度。P232定义2 截断误差Remark：Euler方法是一阶方法，梯形法和Euler预估校正法是二阶方法。第15页/共70页整体截断误差与局部截断误差的关系且局部截断误差有界：则Euler法的整体截断误差n满足估计式：其中L为李普希兹常数，b-a为求解区间长度，定理：如果f(x，y)满足李普希兹（Lipschitz）条
6、件第16页/共70页收敛性与稳定性收敛性定义：如果某一数值方法对于任意固定的xn=x0+nh，当h0(同时n )时有yn y(xn)，则称该方法收敛。定义 用一个数值方法，求解微分方程初值问题时，对给定的步长h0，若在计算 时引入误差 （也称扰动），但由此引起计算后面的 时的误差按绝对值均不增加，则称这个数值方法是稳定的。Remark：该定理表明，整体截断误差比局部截断误差低一阶。对其它方法，也有类似的结论。稳定性定义第17页/共70页 稳定性Remark：由于稳定性问题比较复杂，通常的做法是将满足李普希兹条件的微分方程模型化。设f y=常数，此时微分方程为线性方程 y=y。为保证微分方程的稳
7、定性，假定0。讨论某方法的稳定性，就是讨论该方法对模型方程的稳定性。第18页/共70页 稳定性结论Euler法的稳定性条件是：梯形法是绝对稳定的。Euler预估校正格式的稳定性条件是：对非线性方程，应视，此时将是变化的。的变化将引起h的变化，属于绝对稳定区域，则认为对如果步长h固定，此时，若此方程而言，方法是稳定的。第19页/共70页8.2 泰勒展开法与龙格-库塔（RungeKutta）方法问题：利用泰勒展开法推导高阶单步的求解常微分方程初值问题的数值方法。从提高截断误差阶的阶数入手。第20页/共70页 假定初值问题的解y（x）及函数f（x，y）是充分光滑的，则:当n 充分小时，略去余项 ，则
8、有p阶计算公式 一、Taylor 方法 第21页/共70页其中，第22页/共70页 上式称为p阶Taylor方法。特别地，当p1时，就是Euler公式。当p2时，得二阶Taylor方法：当Taylor方法的阶数p取的较大时，需计算f(x，y)的高阶导数值，计算量较大。特别当f(x，y)较复杂时，y(x)的高阶导数会很复杂。因此Taylor方法很少单独使用，但可以用它来启发思路。第23页/共70页二、RungeKutta 方法 基本思想：用不同点的函数值作线性组合，构造近似公式，把近似公式和解的Taylor展开比较，使前面的若干项吻合，从而使近似公式达到一定的阶数。一般的显式R-K方法，可以写成
9、 第24页/共70页其中，为常数，选取这些常数的原则是，要求第一式的右端在 处泰勒展开后，按h 的幂次重新整理，得到 与微分方程的解的Taylor展开式 有尽可能多的项重合，即要求 第25页/共70页 上述公式叫做N级的Runge-Kutta方法，其局部截断误差为其中表示显然，Euler法是一级一阶R-K方法。下面以二级R-K公式为例，来说明R-K方法的推导过程。第26页/共70页二阶龙格-库塔公式 适当选择，p，使yn+1具有2阶精度注意到第27页/共70页将 在 处展开，有 第28页/共70页而y(xn+1)在xn处的Taylor展式为：将k1，k2表示式代入第29页/共70页Euler预
10、估-校正格式 若取第30页/共70页Remark1：我们可以构造无穷多个二级R-K方法，这些方法的截断误差均为O(h3),即都是二阶方法。其中二阶Heun方法是截断误差项数最少，且允许f 任意变化的情况下截断误差最小的二阶方法。Remark2：二级R-K方法不可能达到三阶Remark3：同样可构造其他阶的R-K方法，它们都有无穷多组解，且三级R-K方法阶数不超过3，四级R-K方法阶数不超过4。Remark4：更高阶的方法由于计算量较大，一般不再采用。第31页/共70页常用的三阶R-K公式（具有三阶精度）第32页/共70页标准（经典）四阶R-K公式（有四阶精度）第33页/共70页关于R-K方法计
11、算量的讨论 二阶R-K方法需计算两个函数值，四阶R-K方法需计算四个函数值，但精度要比二阶方法高出两阶。因此，要达到同样的精度，用低阶方法需步长取得比较小，但若用高阶方法则可以将步长取得大一些，从而降低计算量。四阶经典RK方法的稳定性条件是关于R-K方法稳定性的讨论二阶RK方法的稳定性条件是 第34页/共70页 线性多步法的基本思想：如果充分利用前面多步的信息来预测yn+k，则可以期望获得较高的精度。8.3 线性多步法 前面的RK方法是增加一些非节点处的函数值的计算来提高单步法的精度的，这样使计算量增加了许多。本节介绍多步法，是在不过分增加计算量的情况下取得较高的计算精度。线性多步法公式的构造
12、一般用两种方法，即Taylor展开法与数值积分法。第35页/共70页线性多步法的一般形式式中都为实数，且。当10时上式为隐式方法，当10时，上式为显示方法。由于求yn+1用到前面yn,yn-1,yn-r等r1个值，且关于yn-j和fn-j(j=0,1,2,r)都是线性的，因此称上式为线性r1步方法。第36页/共70页 一、用数值积分方法构造线性多步法将 方程两端从 积分得（1）对 取等距插值节点 ，对应的函数值为 如果k取不同的值，以及F(x)取不同的插值多项式近似，由上式就可以推导出不同的线性多步公式。第37页/共70页其插值余项为：1.阿达姆斯(Adams)外插公式 在(1)式中取k=0，
13、并选择xn,xn-1,xn-2,xn-3作为插值节点，作函数F(x)的三次插值多项式：第38页/共70页把F(x)=L3(x)+R3(x)代入（1）式，有 略去上式右端第三项，得对于上式积分部分用变量代换x=xn+th，并注意到第39页/共70页则第40页/共70页从而得到线性四步Adams显式公式：其局部截断误差就是数值积分的误差第41页/共70页因(x-xn)(x-xn-1)(x-xn-2)(x-xn-3)在xn,xn+1上不变号，并设F(4)(x)在xn,xn+1上连续，利用积分中值定理，存在n xn,xn+1，使得 因为插值多项式L3(x)是在xn3,xn上作出的，而积分区间为xn,x
14、n+1，故上式称为Adams外插公式。第42页/共70页2.阿达姆斯(Adams)内插公式 若在(1)式中取k=2，并选择xn1,xn,xn-1,xn-2作为插值节点，作函数F(x)的三次插值多项式。类似于上面的外插公式，有该公式也称为Adams内插公式，为三步隐式方法。第43页/共70页3.阿达姆斯(Adams)预估-校正公式 由于Adams内插公式是隐式方法，故用它做计算需使用迭代法。通常把Adams外插公式与内插公式结合起来使用，先由前者提供初值，再由后者进行修正，即第44页/共70页当在求解区域内成立时，迭代收敛。若上式中的第二式只迭代一次，便得到Adams预估-校正格式。第45页/共
15、70页 Taylor展开法更具一般性，不仅可以构造用数值积分法得出的数值方法，而且还可导出积分法得不到的方法。它比积分法更加灵活。下面仅举一例说明如何用这种方法构造线性多步法。二、用Taylor展开法构造线性多步公式 首先以xn-1,xn,xn+1为节点，构造形如的公式。假设上式右边第46页/共70页将函数在x=xn处展开，有:代入给定公式并按h的幂次整理得到下式：第47页/共70页第48页/共70页将上式与比较，选择系数i(i=0,1)和i(i=-1,0,1)使两式中关于h的同次幂的系数有尽可能多的项相等。故有：第49页/共70页求解上述方程组，得出0，1，1，0，1。所得到的算式的局部截断
16、误差为O(h5)。Reamrk：我们也可以只要求前面几个方程组成立，如要求前面4个方程组成立时，所得算式的局部截断误差为O(h4)。如令00，带入上式的前4个方程，解得1=1，1=1=1/3，0=4/3,于是得到计算公式为：第50页/共70页 此时上式中第5式也恰巧成立。可以得到上式得截断误差为：称上式为辛浦生（Simpson）公式，它可由数值积分方法而得到。我们也可以用类似的方法构造其它的线性多步法，如前面的Adams公式等。第51页/共70页三、出发值的计算 使用线性k步法求解初值问题时，需要知道k个出发值y0,y1,yk-1才能进行计算。然而初值问题只提供一个yn,还有k1个出发值，需要
17、通过别的方法计算出来。常用的方法是一步方法。由于初值对于确定微分方程的解有重要作用，因而在求解数值解时，对出发值的精度也必须有相应的要求。为了保证多步方法的精确度，用于计算出发值的一步方法的阶数至少不低于多步方法的阶。第52页/共70页 理论上讲，可用Taylor展开法和Runge-Kutta方法，计算出发值。但由于Taylor展开法要计算高阶导数值，故最常用的方法还是选择与多步法同阶的Runge-Kutta方法。一旦出发值计算出来，线性多步法的计算量（特别是显式公式）就会很小，因为每次只须计算一次f值。Remark：有关线性多步法的整体截断误差、收敛性及数值稳定性的讨论可参考有关文献。第53
18、页/共70页8.4数值算例求解常微分方程初值问题在0,1上的解，取步长h0.1。计算结果如下图所示：第54页/共70页第55页/共70页第56页/共70页第57页/共70页第58页/共70页第59页/共70页第60页/共70页返回第61页/共70页第8章例题例1：用欧拉预校方法求解初值问题（要求取步长h=0.2，计算y(1.2)和y(1.4)的近似值，小数点后保留5位小数）：解:欧拉预校格式为:第62页/共70页由y(1)=y0=1计算得:于是有:第63页/共70页例2:用二阶泰勒展开法求初值问题:解:二阶泰勒展开为:因为:代入上式并略去高阶项o(h3)，则得求解公式为：由y(1)=y0=1计
19、算得:求x=1.5的近似值(取步长为0.25,小数点后保留5位)第64页/共70页例3（作业5）证明求解初值问题的如下单步方法是二阶方法。在的假定下，证明：第65页/共70页因此有故所给方法是二阶方法。而对直接在处泰勒展开有：证毕第66页/共70页例4 设求解常微分方程初值问题的如下线性二步格式：其中：，试确定参数，使该格式为三阶格式。解：为考虑局部截断误差，设于是所给格式可以写为：（1）分别将在xn处泰勒展开，有:第67页/共70页代入（1）式并按h的幂次整理后有：（2）比较（2）（3）式h幂次相同的项并令其系数相等有：解得：而由泰勒展开有：（3）第68页/共70页此时（3）式与（2）相减有：即所给格式为三阶格式，具体为：#第69页/共70页感谢您的观看！第70页/共70页