西安石油大学现代数值计算方法第8章.pptx

常微分方程初值问题的数值解法常微分方程初值问题的数值解法8.8.1 1 欧拉法与梯形法欧拉法与梯形法8.8.2 2 泰勒展开法与龙格泰勒展开法与龙格-库塔库塔（RungeRungeKuttaKutta）方法）方法8.3 8.3 线性多步法线性多步法第八章第八章8.08.0 概述概述8.4 8.4 数值算例数值算例本章着重讨论一阶常微分方程初值问题本章着重讨论一阶常微分方程初值问题 的数值解法。的数值解法。8.08.0 概述概述 常微分方程初值问题的数值解是求上述初常微分方程初值问题的数值解是求上述初值问题的解值问题的解y y(x x)在区间在区间 a a，b b 中的点中的点列列上的近似值上的近似值 .以下设以下设不变，记为不变，记为h-h-步长步长 。定理定理：如果：如果f f(x x，y y)满足李普希兹（满足李普希兹（LipschitzLipschitz）条件）条件则上述微分方程有唯一解则上述微分方程有唯一解y(x)假设解假设解y y(x x)在区间在区间 a a,b b 上是存在而且唯一的，上是存在而且唯一的，并且具有充分的光滑度，因此，要求并且具有充分的光滑度，因此，要求f f(x,yx,y)也充也充分光滑。初值问题的解析解分光滑。初值问题的解析解(理论解）用理论解）用 表示表示,数值解法的精确解用数值解法的精确解用 表示。表示。常微分方程数值解法一般分为：常微分方程数值解法一般分为：（1）一步法：在计算 时，只用到 ,和 ，即前一步的值。（2）多步法：计算 时，除用到 ,和 以外，还要用 和 ,即前 k步的值。（3）显式格式与隐式格式。8.8.1 1 欧拉法与梯形法欧拉法与梯形法设节点为设节点为 ，得欧拉方法计算公得欧拉方法计算公式为：式为：一、欧拉一、欧拉(Euler)(Euler)法法下面通过几种常用的方法来推导该公式。下面通过几种常用的方法来推导该公式。1 1、泰勒展开法、泰勒展开法假设在假设在 附近把附近把y y（x x）做做TaylorTaylor展开，有：展开，有：取取h h的线性部分的线性部分,并用并用 表示表示 的近似值的近似值,得得2 2、数值积分法、数值积分法 从 到 +h对等式 y(t)=f(t,y(t)进行积分得到 再利用左矩形公式，得再利用左矩形公式，得从而得到从而得到EulerEuler公式。公式。由由3 3、数值微分法、数值微分法4 4、几何方法、几何方法过点过点(x xn n,y yn n)作以作以f f(x xn n ,y yn n)为斜率的直线方程：为斜率的直线方程：将将x x=x xn+n+1 1处该直线上的函数值做为处该直线上的函数值做为y y(x xn n+1+1)的近似值，的近似值，则有则有EulerEuler公式。这实质上是在每个小区间上利用折公式。这实质上是在每个小区间上利用折线来代替曲线的结果，故线来代替曲线的结果，故EulerEuler法法又称又称EulerEuler折线法折线法。二、梯形法二、梯形法在式在式 中中,将积分用将积分用梯形公式来代替，则有梯形公式来代替，则有 从而得到梯形公式：从而得到梯形公式：梯形方法关于梯形方法关于y yn n+1+1是隐式的，而是隐式的，而EulerEuler方法是显方法是显式的。一般情形下不容易从上式解出式的。一般情形下不容易从上式解出y yn n+1+1，因而可将，因而可将上式与上式与EulerEuler公式联合使用，即公式联合使用，即使用上式时，先用第一式算出使用上式时，先用第一式算出x xn n+1+1处处y yn n+1+1的初始近似的初始近似再用第二式反复迭代，得到数列再用第二式反复迭代，得到数列用用来控制迭代次数，这里来控制迭代次数，这里 为允许误差。把满足误差要求的为允许误差。把满足误差要求的可以证明可以证明,当当f f(x,yx,y)满足满足LipschitzLipschitz条件条件,即：即：(L L为为LipschitzLipschitz常数常数)时时,上述数列收敛。上述数列收敛。作为作为y y(x xn n+1+1)的近似值的近似值y yn n+1+1.类似地可以得出类似地可以得出y yn+2n+2,y,yn+3n+3,证明：证明：由由和和有：有：反复使用不等式有：反复使用不等式有：实用中实用中,在在h h 取得较小时取得较小时,用梯形公式计算用梯形公式计算,第第二式只迭代一次就结束二式只迭代一次就结束,得到得到EulerEuler预估预估-校正格式校正格式:第一式称为第一式称为预估公式预估公式，第二式称为，第二式称为校正公式校正公式。三、三、EulerEuler预估预估-校正格式校正格式四、方法的误差估计、收敛性和稳定性四、方法的误差估计、收敛性和稳定性 定义1：为 某一数值方法在xn处的整体截断误差(不考虑舍入误差的影响)。定义2：对单步法，在 的假设下，称为在 处的局部截断误差。（P232定义1）Remark1:Euler法的局部截断误差为法的局部截断误差为(由泰勒余项由泰勒余项):Remark2:梯形方法的局部截断误差为（由梯形积分）用泰勒展开法推导用泰勒展开法推导EulerEuler预估校正预估校正格式的局部截断误差格式的局部截断误差改写改写EulerEuler预估校正公式为：预估校正公式为：在在的假定下，的假定下，而而而而因此有因此有故故E Euleruler预估校正方法为的局部截断误差阶为预估校正方法为的局部截断误差阶为O O(h h3 3)。定义定义3 3：若一个方法的局部截断误差为：若一个方法的局部截断误差为 ，则则称该方法为称该方法为p p阶方法阶方法,或称该方法具有或称该方法具有p p阶精度阶精度。P232P232定义定义2 2 截断误差截断误差RemarkRemark：EulerEuler方法是一阶方法，梯形法和方法是一阶方法，梯形法和EulerEuler预估校正法是二阶方法。预估校正法是二阶方法。整体截断误差与局部截断误差的关系整体截断误差与局部截断误差的关系且局部截断误差有界：且局部截断误差有界：则则EulerEuler法的整体截断误差法的整体截断误差 n n满足估计式：满足估计式：其中其中L L为李普希兹常数，为李普希兹常数，b-ab-a为求解区间长度，为求解区间长度，定理定理：如果：如果f f(x x，y y)满足李普希兹（满足李普希兹（LipschitzLipschitz）条件）条件收敛性与稳定性收敛性与稳定性收敛性定义收敛性定义：如果某一数值方法对于任意固定的：如果某一数值方法对于任意固定的x xn n=x x0 0+nhnh，当，当h h0(0(同时同时n n )时有时有y yn n y y(x xn n)，则，则称该方法称该方法收敛收敛。定义定义 用一个数值方法，求解微分方程初值问题时，对给定的步长h0，若在计算 时引入误差 （也称扰动），但由此引起计算后面的 时的误差按绝对值均不增加，则称这个数值方法是稳定的稳定的。RemarkRemark：该定理表明，整体截断误差比局部截断：该定理表明，整体截断误差比局部截断误差低一阶。对其它方法，也有类似的结论。误差低一阶。对其它方法，也有类似的结论。稳定性定义稳定性定义 稳定性稳定性RemarkRemark：由于稳定性问题比较复杂，通常的：由于稳定性问题比较复杂，通常的做法是将满足李普希兹条件的微分方程模型做法是将满足李普希兹条件的微分方程模型化。设化。设 f f y y=常数，此时微分方程为线性常数，此时微分方程为线性方程方程 y y=y y。为保证微分方程的稳定性，。为保证微分方程的稳定性，假定假定 00。讨论某方法的稳定性，就是讨论。讨论某方法的稳定性，就是讨论该方法对模型方程的稳定性。该方法对模型方程的稳定性。稳定性结论稳定性结论EulerEuler法的稳定性条件是：法的稳定性条件是：梯形法梯形法是绝对稳定的。是绝对稳定的。EulerEuler预估校正格式的稳定性条件是：预估校正格式的稳定性条件是：对对非线性方程非线性方程，应视，应视，此时，此时 将是变化的。将是变化的。的变化将引起的变化将引起 h h的变化，的变化，属于绝对稳定区域，则认为对属于绝对稳定区域，则认为对如果步长如果步长h h固定，固定，此时，若此时，若此方程而言，方法是稳定的。此方程而言，方法是稳定的。8.8.2 2 泰勒展开法与龙格泰勒展开法与龙格-库库塔（塔（RungeRungeKuttaKutta）方法）方法问题问题：利用泰勒展开法推导高阶单步的：利用泰勒展开法推导高阶单步的求解常微分方程初值问题的数值方法。求解常微分方程初值问题的数值方法。从提高截断误差阶的阶数入手。从提高截断误差阶的阶数入手。假定初值问题的解假定初值问题的解y y（x x）及函数及函数f f（x x，y y）是充分光滑的，则是充分光滑的，则:当当n n 充分小时，略去余项充分小时，略去余项 ，则有，则有p p阶计算公式阶计算公式 一、一、Taylor Taylor 方法方法 其中，其中，上式称为上式称为p p阶阶TaylorTaylor方法方法。特别地，当。特别地，当p p1 1时，就是时，就是EulerEuler公式。当公式。当p p2 2时，得二时，得二阶阶TaylorTaylor方法：方法：当当TaylorTaylor方法的阶数方法的阶数p p取的较大时，需计取的较大时，需计算算f f(x x，y y)的高阶导数值，计算量较大。特别的高阶导数值，计算量较大。特别当当f f(x x，y y)较复杂时，较复杂时，y y(x)(x)的高阶导数会很复的高阶导数会很复杂。因此杂。因此TaylorTaylor方法很少单独使用，但可以方法很少单独使用，但可以用它来启发思路。用它来启发思路。二、二、RungeRungeKuttaKutta 方法方法 基本思想基本思想：用不同点的函数值作线性组合，构：用不同点的函数值作线性组合，构造近似公式，把近似公式和解的造近似公式，把近似公式和解的TaylorTaylor展开比较，展开比较，使前面的若干项吻合，从而使近似公式达到一定使前面的若干项吻合，从而使近似公式达到一定的阶数。一般的显式的阶数。一般的显式R-KR-K方法，可以写成方法，可以写成 其中，其中，为常数，选取这些常数的原则是，为常数，选取这些常数的原则是，要求第一式的右端在要求第一式的右端在 处泰勒展开后，按处泰勒展开后，按h 的的幂次重新整理，得到幂次重新整理，得到 与微分方程的解的与微分方程的解的TaylorTaylor展开式展开式 有尽可能多的项重合，即要求有尽可能多的项重合，即要求 上述公式叫做上述公式叫做N N级的级的Runge-KuttaRunge-Kutta方法，其局部方法，其局部截断误差为截断误差为其中其中表示表示显然，显然，EulerEuler法是一级一阶法是一级一阶R-KR-K方法。方法。下面以二级下面以二级R-KR-K公式为例，来说明公式为例，来说明R-KR-K方法的推导方法的推导过程。过程。二阶龙格-库塔公式 适当选择适当选择，p，使，使yn+1具有具有2阶精度阶精度注意到注意到将将 在在 处展开，有处展开，有 而而y y(x xn n+1+1)在在x xn n处的处的TaylorTaylor展式为：展式为：将将k1，k2表示式代入表示式代入EulerEuler预估预估-校正格式校正格式 若取若取Remark1Remark1：我们可以构造无穷多个二级我们可以构造无穷多个二级R-KR-K方方法，这些方法的截断误差均为法，这些方法的截断误差均为O O(h h3 3),),即都是二即都是二阶方法。其中二阶阶方法。其中二阶HeunHeun方法是截断误差项数最方法是截断误差项数最少，且允许少，且允许f f 任意变化的情况下截断误差最小任意变化的情况下截断误差最小的二阶方法。的二阶方法。Remark2Remark2：二级二级R-KR-K方法不可能达到三阶方法不可能达到三阶Remark3Remark3：同同样样可可构构造造其其他他阶阶的的R-KR-K方方法法，它它们们都都有有无无穷穷多多组组解解，且且三三级级R-KR-K方方法法阶阶数数不不超超过过3 3，四级，四级R-KR-K方法阶数不超过方法阶数不超过4 4。Remark4Remark4：更更高高阶阶的的方方法法由由于于计计算算量量较较大大，一一般不再采用。般不再采用。常用的三阶常用的三阶R-KR-K公式（具有三阶精度）公式（具有三阶精度）标准（经典）四阶标准（经典）四阶R-KR-K公式（有四阶精度）公式（有四阶精度）关于关于R-KR-K方法计算量的讨论方法计算量的讨论 二阶二阶R-KR-K方法需计算两个函数值，四阶方法需计算两个函数值，四阶R-KR-K方法需方法需计算四个函数值，但精度要比二阶方法高出两阶。因计算四个函数值，但精度要比二阶方法高出两阶。因此，要达到同样的精度，用低阶方法需步长取得比较此，要达到同样的精度，用低阶方法需步长取得比较小，但若用高阶方法则可以将步长取得大一些，从而小，但若用高阶方法则可以将步长取得大一些，从而降低计算量。降低计算量。四阶经典四阶经典R RK K方法的稳定性条件是方法的稳定性条件是关于关于R-KR-K方法稳定性的讨论方法稳定性的讨论二阶二阶R RK K方法的稳定性条件是方法的稳定性条件是 线性多步法的线性多步法的基本思想基本思想：如果充分利用：如果充分利用前面多步的信息来预测前面多步的信息来预测y yn n+k k，则可以期望获，则可以期望获得较高的精度。得较高的精度。8.3 8.3 线性多步法线性多步法 前面的前面的R RK K方法是增加一些非节点处的函方法是增加一些非节点处的函数值的计算来提高单步法的精度的，这样使计数值的计算来提高单步法的精度的，这样使计算量增加了许多。本节介绍多步法，是在不过算量增加了许多。本节介绍多步法，是在不过分增加计算量的情况下取得较高的计算精度。分增加计算量的情况下取得较高的计算精度。线性多步法公式的线性多步法公式的构造构造一般用两种方法，一般用两种方法，即即TaylorTaylor展开法与数值积分法。展开法与数值积分法。线性多步法的一般形式线性多步法的一般形式式中式中都为实数，且都为实数，且。当。当 1 1 0 0时上式为隐式方法，当时上式为隐式方法，当 1 10 0时，上式为显示方法。由于求时，上式为显示方法。由于求y yn n+1+1用到前面用到前面y yn n,y yn n-1 1,y yn n-r-r等等r r1 1个值，且关于个值，且关于y yn n-j j和和f fn n-j j(j j=0,1,2,=0,1,2,r r)都是线性的，因此称上式为线性都是线性的，因此称上式为线性r r1 1步方法。步方法。一、用数值积分方法构造线性多步法用数值积分方法构造线性多步法将将 方程两端从方程两端从 积分得积分得 （1）对对 取等距插值节点取等距插值节点 ，对应对应的函数值为的函数值为 如果如果k k取取不同的值，以及不同的值，以及F F(x x)取不同的插值多项式近似，由取不同的插值多项式近似，由上式就可以推导出不同的线性多步公式。上式就可以推导出不同的线性多步公式。其插值余项为：其插值余项为：1.1.阿达姆斯阿达姆斯(Adams)(Adams)外插公式外插公式 在在(1)(1)式中取式中取k k=0=0，并选择，并选择x xn n,x xn-n-1 1,x xn n-2-2,x xn n-3-3作为作为插值节点，作函数插值节点，作函数F F(x x)的三次插值多项式：的三次插值多项式：把把F F(x x)=)=L L3 3(x x)+)+R R3 3(x x)代入（代入（1 1）式，有）式，有 略去上式右端第三项，得略去上式右端第三项，得对于上式积分部分用变量代换对于上式积分部分用变量代换x=xx=xn n+th+th，并注意到，并注意到则则从而得到线性四步从而得到线性四步AdamsAdams显式公式：显式公式：其局部截断误差就是数值积分的误差其局部截断误差就是数值积分的误差因因(x x-x xn n)()(x x-x xn n-1-1)()(x x-x xn n-2-2)()(x x-x xn n-3-3)在在 x xn n,x xn n+1+1 上不上不变号，并设变号，并设F F(4)(4)(x x)在在 x xn n,x xn n+1+1 上连续，利用积分中上连续，利用积分中值定理，存在值定理，存在 n n x xn n,x xn n+1+1，使得，使得 因为插值多项式因为插值多项式L L3 3(x x)是在是在 x xn n3 3,x xn n 上作出的，上作出的，而积分区间为而积分区间为 x xn n,x xn n+1+1，故上式称为，故上式称为AdamsAdams外插公外插公式。式。2.2.阿达姆斯阿达姆斯(Adams)(Adams)内插公式内插公式 若在若在(1)(1)式中取式中取k k=2=2，并选择，并选择x xn n1 1,x xn n,x xn n-1-1,x xn n-2-2作为插值节点，作函数作为插值节点，作函数F F(x x)的三次插值多项式。的三次插值多项式。类似于上面的外插公式，有类似于上面的外插公式，有该公式也称为该公式也称为AdamsAdams内插公式，为三步隐式方法。内插公式，为三步隐式方法。3.3.阿达姆斯阿达姆斯(Adams)(Adams)预估预估-校正公式校正公式 由于由于AdamsAdams内插公式是隐式方法，故用它做计内插公式是隐式方法，故用它做计算需使用迭代法。通常把算需使用迭代法。通常把AdamsAdams外插公式与内插公外插公式与内插公式结合起来使用，先由前者提供初值，再由后者式结合起来使用，先由前者提供初值，再由后者进行修正，即进行修正，即当在求解区域内成立时，迭代收敛。在求解区域内成立时，迭代收敛。若上式中的第二式只迭代一次，便得到若上式中的第二式只迭代一次，便得到AdamsAdams预预估估-校正格式。校正格式。TaylorTaylor展展开开法法更更具具一一般般性性，不不仅仅可可以以构构造造用用数数值值积积分分法法得得出出的的数数值值方方法法，而而且且还还可可导导出出积积分分法法得得不不到到的的方方法法。它它比比积积分分法法更更加加灵灵活活。下下面面仅仅举举一一例例说明如何用这种方法构造线性多步法。说明如何用这种方法构造线性多步法。二、用二、用TaylorTaylor展开法构造线性多步公式展开法构造线性多步公式 首先以首先以x xn n-1-1,x xn n,x xn n+1+1为节点，构造形如为节点，构造形如的公式。假设上式右边的公式。假设上式右边将函数将函数在在x x=x xn n处展开，有处展开，有:代入给定公式并按代入给定公式并按h的幂次整理得到下式：的幂次整理得到下式：将上式与将上式与比较，选择系数比较，选择系数 i i(i i=0,1)=0,1)和和 i i(i i=-1,0,1)=-1,0,1)使两式中关使两式中关于于h h的同次幂的系数有尽可能多的项相等。故有：的同次幂的系数有尽可能多的项相等。故有：求解上述方程组，得出求解上述方程组，得出 0 0，1 1，1 1，0 0 ，1 1。所。所得到的算式的局部截断误差为得到的算式的局部截断误差为O O(h h5 5)。ReamrkReamrk：我们也可以只要求前面几个方程组成立，：我们也可以只要求前面几个方程组成立，如要求前面如要求前面4 4个方程组成立时，所得算式的局部截断个方程组成立时，所得算式的局部截断误差为误差为O O(h h4 4)。如令。如令 0 00 0，带入上式的前，带入上式的前4 4个方程，个方程，解得解得 1 1=1=1，1 1=1 1=1/3=1/3 ，0 0=4/3,=4/3,于是得到计算公式于是得到计算公式为：为：此时上式中第此时上式中第5 5式也恰巧成立。可以得到上式式也恰巧成立。可以得到上式得截断误差为：得截断误差为：称上式为称上式为辛浦生（辛浦生（SimpsonSimpson）公式）公式，它可由数值积，它可由数值积分方法而得到。分方法而得到。我们也可以用类似的方法构造其它的线性多步法，我们也可以用类似的方法构造其它的线性多步法，如前面的如前面的AdamsAdams公式等。公式等。三、出发值的计算三、出发值的计算 使用线性使用线性k k步法求解初值问题时，需要知道步法求解初值问题时，需要知道k k个个出发值出发值y y0 0,y y1 1,y yk k-1-1才能进行计算。然而初值问题才能进行计算。然而初值问题只提供一个只提供一个y yn n,还有还有k k1 1个出发值，需要通过别的方个出发值，需要通过别的方法计算出来。常用的方法是一步方法。由于初值对法计算出来。常用的方法是一步方法。由于初值对于确定微分方程的解有重要作用，因而在求解数值于确定微分方程的解有重要作用，因而在求解数值解时，对出发值的精度也必须有相应的要求。为了解时，对出发值的精度也必须有相应的要求。为了保证多步方法的精确度，用于计算出发值的一步方保证多步方法的精确度，用于计算出发值的一步方法的阶数至少不低于多步方法的阶。法的阶数至少不低于多步方法的阶。理论上讲，可用理论上讲，可用TaylorTaylor展开法和展开法和Runge-Runge-KuttaKutta方法，计算出发值。但由于方法，计算出发值。但由于TaylorTaylor展开展开法要计算高阶导数值，故最常用的方法还是选法要计算高阶导数值，故最常用的方法还是选择与多步法同阶的择与多步法同阶的Runge-KuttaRunge-Kutta方法。一旦出方法。一旦出发值计算出来，线性多步法的计算量（特别是发值计算出来，线性多步法的计算量（特别是显式公式）就会很小，因为每次只须计算一次显式公式）就会很小，因为每次只须计算一次f f值。值。RemarkRemark：有关线性多步法的整体截断误差、：有关线性多步法的整体截断误差、收敛性及数值稳定性的讨论可参考有关文献。收敛性及数值稳定性的讨论可参考有关文献。8.48.4数值算例数值算例求解常微分方程初值问题求解常微分方程初值问题在在0,10,1上的解，取步长上的解，取步长h h0.10.1。计算结果如下图所示：计算结果如下图所示：返回第第8 8章例题章例题例例1 1：用欧拉预：用欧拉预校方法求解初值问题（要求取步长校方法求解初值问题（要求取步长h=0.2h=0.2，计算，计算y(1.2)y(1.2)和和y(1.4)y(1.4)的近似值，小数点后保的近似值，小数点后保留留5 5位小数）：位小数）：解解:欧拉预欧拉预校格式为校格式为:由由y(1)=y0=1计算得计算得:于是有于是有:例例2:用二阶泰勒展开法求初值问题用二阶泰勒展开法求初值问题:解解:二阶泰勒展开为二阶泰勒展开为:因为因为:代入上式并略去高阶项代入上式并略去高阶项o(ho(h3 3)，则得求解公式为：，则得求解公式为：由由y(1)=y0=1计算得计算得:求求x=1.5的近似值的近似值(取步长为取步长为0.25,小数点后保留小数点后保留5位位)例例3 3（作业（作业5 5）证明求解初值问题）证明求解初值问题的如下单步方法是二阶方法。的如下单步方法是二阶方法。在在的假定下，的假定下，证明：证明：因此有因此有故故所给方法所给方法是二阶方法。是二阶方法。而对而对直接在直接在处泰勒展开有：处泰勒展开有：证毕证毕例例4 4 设求解常微分方程初值问题设求解常微分方程初值问题的如下线性二步格式：的如下线性二步格式：其中：其中：，试确定参数，试确定参数，使该格式为三阶格式。，使该格式为三阶格式。解：为考虑局部截断误差，设解：为考虑局部截断误差，设于是所给格式可以写为：于是所给格式可以写为：（1）分别将分别将在在x xn n处泰勒展开，有处泰勒展开，有:代入（代入（1 1）式并按）式并按h h的幂次整理后有：的幂次整理后有：（2）比较（比较（2 2）（）（3 3）式）式h h幂次相同的项并令其系数相等有：幂次相同的项并令其系数相等有：解得：解得：而由泰勒展开有：而由泰勒展开有：（3）此时（此时（3 3）式与（）式与（2 2）相减有：）相减有：即所给格式为三阶格式，具体为：即所给格式为三阶格式，具体为：#