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

第5章 常微分方程数值解法,5.1 引言 包含自变量、未知函数及未知函数的导数或微分的方程称为微分方程.在微分方程中, 自变量的个数只有一个, 称为常微分方程.自变量的个数为两个或两个以上的微分方程叫偏微分方程.微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数.如果未知函数y及其各阶导数 都是一次的,则称它是线性的,否则称为非线性的.,在常微分方程中，对于常微分方程的求解，给出了一些典型方程求解析解的基本方法，如可分离变量法、常系数齐次线性方程的解法、常系数非齐次线性方程的解法等.但能求解的常微分方程仍然是有限的，大多数的常微分方程是不可能给出解析解. 譬如,这个一阶微分方程就不能用初等函数及其积分来表达它的解.,再如，方程,的解 ,虽然有表可查,但对于表上没有给出 的值,仍需插值方法来计算,从实际问题当中归纳出来的微分方程，通常主要依靠数值解法来解决.本章主要讨论一阶常微分方程初值问题,( 5.1 ),在区间axb上的数值解法.,数值方法的基本思想 对常微分方程初值问题(5.1)式的数值解法，就是要算出精确解y(x)在区间a,b上的一系列离散节点 处的函数值 的近似值,相邻两个节点的间距 称为步长，步 长可以相等，也可以不等.本章总是假定h为定数，称为定步长，这时节点可表示为,数值解法需要把连续性的问题加以离散化，从而求出离散节点的数值解.,对常微分方程数值解法的基本出发点就是离散化.其数值解法有个基本特点，它们都采用“步进式”，即求解过程顺着节点排列的次序一步一步地向前推进，描述这类算法，要求给出用已知信息 计算 的递推公式. 建立这类递推公式的基本方法是在这些节点上用数值积分、数值微分、泰勒展开等离散化方法，对初值问题 中的导数 进行不同的离散化处理.,差分格式,对于初值问题 的数值解法，首先要解决的问题就是如何对微分方程进行离散化，建立求数值解的递推公式. 递推公式通常有两类，一类是计算yn+1时只用到xn+1, xn和yn，即前一步的值，因此有了初值以后就可以逐步往下计算，此类方法称为单步法. 另一类是计算yn+1时，除用到xn+1,xn和yn以外，还要用到 即前面k步的值，此类方法称为多步法.,5.2 Euler方法 5.2.1 Euler格式 Euler方法是解初值问题的最简单的数值方法. 初值问题 的解y=y(x)称为它的积分曲线.积分曲线上每一点 处的切线的斜率 等于函数 在这点的值.,即有,在点xn处为,取,则有,设y(xn)的近似值yn已知，用它带入上式右端进行计算，并取计算结果yn+1作为y(xn+1)的近似值,则可得,( 5.2 ),Euler格式,Euler法的求解过程是:从初始点P0(即点(x0,y0)出发, 作积分曲线y=y(x)在P0点上切线 (其斜率为 ),与x=x1直线,相交于P1点(即点(x1,y1),得到y1作为y(x1)的近似值,如上图所示.过点(x0,y0),以f(x0,y0)为斜率的切线方程为 当 时,得,这样就获得了P1点的坐标.,同样, 过点P1(x1,y1),以斜率 作直线交x=x2于P2点,线段 所在的直线方程为,当 时,得,由此获得了P2的坐标.,当 时,得,重复以上过程,就可获得一系列的点:P1,P1,Pn.对已求得点 以 为斜率作直线,取,从图形上看,就获得了一条近似于曲线y=y(x) 的折线,这样,从x0逐个算出 对应的数值解,Euler格式算法实现 （1）计算步骤 输入 并计算出 使用Euler格式进行计算 输出 ，并使 转到 直至n N 结束.,（2）Euler格式的流程图,例5.1 用Euler格式解初值问题,取步长h=0.1.,Euler格式的具体形式为,计算公式的精度,常以Taylor展开为工具来分析计算公式的精度.,为简化分析，假定yn是准确的，即在 的前提下估计误差,局部截断误差,Euler格式的局部截断误差,由,从而,对,在点xn+1处为,取,则有,设y(xn)的近似值yn已知，用它带入上式右端进行计算，并取计算结果yn+1作为y(xn+1)的近似值,则可得,( 5.3 ),后退的Euler格式,5.2.2 后退的Euler格式,隐式方程(5.3) 常用迭代法求解,用它带入(5.3)的右端，使之转化为显式，直接计算得,然后再用 带入(5.3)的右端，又有,如此反复进行迭代，得,后退的Euler格式的局部截断误差(假设 ),（在xn点Taylor展开）,从而有,移项整理得,可得,后退的Euler格式的局部截断误差,Euler格式的局部截断误差,5.2.3 梯形格式 比较Euler格式和后退的Euler格式的误差公式，可以看到，如果将这两种方法进行算术平均，即可消除误差的主要部分 而获得更高的精度.即,( 5.4 ),梯形方法,梯形方法的局部截断误差(假设 ),（在xn点Taylor展开）,从而有,移项整理得,可得,梯形方法(5.4)是隐式的，可用迭代法求解,梯形法的迭代公式为,注：迭代过程是收敛的,( 5.5 ),5.2.4 改进的欧拉公式 显式Euler公式计算工作量小,但精度低.梯形公式虽提高了精度,但为隐式公式,需用迭代法求解,计算工作量大.综合Euler公式和梯形公式便可得到改进的Euler公式. 先用欧拉公式(5.2)求出一个初步的近似值 称为预测值, 它的精度不高, 再用梯形公式(5.4)对它校正一次,即迭代一次,求得yn+1,称为校正值, 这种预测-校正方法称为改进的Euler格式：,预测,校正,它可以表示为嵌套形式.,或者表示成下列平均化形式,( 5.6 ),改进Euler格式算法实现 （1）计算步骤 输入 并计算出 使用改进Euler格式进行计算 输出 ，并使 转到 直至n N 结束.,（2）改进Euler格式的流程图,例5.2 用改进的Euler方法解初值问题,取步长h=0.1.,解: 改进Euler方法的具体形式,5.2.5 Euler两步格式 改进的Euler格式中,下面构造与梯形方法在精度上相匹配的显式方法,在点xn处为,则有,设y(xn-1),y(xn) 的近似值yn-1,yn已知，用它带入上式右端进行计算，并取计算结果yn+1作为y(xn+1)的近似值,则可得,( 5.7 ),Euler两步格式,取,局部截断误差（设 ）,前面介绍过的数值方法,无论是Euler方法,后退的Euler方法, 梯形方法，还是改进的Euler方法,它们都是单步法,其特点是在计算yn+1时只用到前一步的信息yn 可是公式(5.7)中除了yn外,还用到更前一步的信息yn-1,即调用了前两步的信息,故称其为Euler两步法,( 5.8 ),由前面知对于(5.8),( 5.9 ),对于(5.9),假定 则,从而有,可得事后估计式,可以期望，利用这样估计出的误差作为计算结果的一种补偿，有可能使精度得到改善,设以pn和cn分别表示第n步的预估值和校正值,预测,改进,校正,改进,设以pn和cn分别表示第n步的预估值和校正值,预测,改进,校正,改进,计算,计算,注：算法开始之前必须先给出y0,y1和p1-c1,5.3 Runge-Kutta方法,定义5.1：如果一种方法的局部截断误差为 则称该方法具有p阶精度.,5.3.2 Runge-Kutta方法的基本思想,则,由微分中值定理有,y(x)=f(x,y(x),设,称K*为区间xn,xn+1上的平均斜率,注:只要对K*提供一种算法，就相应的导出一种计算格式,Euler格式 可改写成,改进的Euler格式,可改写成,这个处理过程提示我们，如果设法在(xn,xn+1)内多预测几个点的斜率值，然后将它们加权平均作为平均斜率K*，则有可能构造出具有更高精度的计算格式,Runge-Kutta方法的基本思想,5.3.3 二阶Runge-Kutta方法 在 上取两点xn和 以这两点处的斜率值K1和K2的加权平均(或称为线性组合)来求取平均斜率K*的近似值K，即,式中:K1为xn点处的切线斜率值， K2为 点处的切线斜率值,比照改进的Euler法,将 视为 ，即可得,这样设计出的计算格式具有如下形式：,上式含有三个待定系数 下面通过适当选取这些系数的值，使得格式具有二阶精度.,( 5.10 ),将以上结果代入,得,设 则,其中,（在(xn ,yn)Taylor展开）,比较系数后可知,只要,成立,格式(5.10)的局部截断误差就等于,二阶Runge-Kutta格式,( 5.11 ),(1)若取 则 这是无穷多解中的一个解，将以上所解的值代入式(5.10)并改写可得,(2)若取 则 此时二阶Runge-Kutta法的计算公式为,此计算公式称为变形的Euler格式.,5.3.4 三阶Runge-Kutta方法,为了进一步提高精度，设除 外再增加一点,并用三个点 的斜率K1,K2,K3加权平均 得出平均斜率K*的近似值，这时计算格式具有形式:,为了预报点 的斜率值K3,在区间 内有两 个斜率值K1和K2可以用,可将K1,K2加权平均得出 上的平均斜率,从而得到 的预报值,这样设计出的计算格式具有如下形式：,可以通过适当选取 这些系数的值，使得格式具有三阶精度.,( 5.12 ),从而,和前面的推导类似，只要,成立,格式(5.12)的局部截断误差就等于,三阶Runge-Kutta格式,( 5.13 ),下列是其中的一种，称为Kutta公式.,5.3.5 四阶Runge-Kutta方法,如果需要再提高精度，用类似上述的处理方法，只需在区间 上用四个点处的斜率加权平均作为平均斜率K*的近似值，构成一系列四阶Runge-Kutta塔公式. 具有四阶精度，即局部截断误差是 .,龙格库塔方法的推导基于Taylor展开方法，因而它要求所求的解具有较好的光滑性.如果解的光滑性差，那么，使用四阶龙格库塔方法求得的数值解，其精度可能反而不如改进的欧拉方法.在实际计算时，应当针对问题的具体特点选择合适的算法.,9.3.5 四阶龙格库塔法算法实现 (1) 计算步骤 输入 ,h,N 使用龙格库塔公式（9.20）计算出y1 输出 ，并使 转到 直至n N 结束.,(2) 四阶龙格库塔算法流程图,程序实现(四阶龙格-库塔法计 算常微分方程初值问题),例5.4 取步长h=0.2，用经典格式求解初值问题,解: 由四阶龙格-库塔公式可得,可同样进行其余yi的计算.本例方程的解为 ，从表中看到所求的数值解具有4位有效数字.,5.3.6 变步长的龙格-库塔法 在微分方程的数值解中，选择适当的步长是非常重要的.单从每一步看，步长越小，截断误差就越小；但随着步长的缩小，在一定的求解区间内所要完成的步数就增加了.这样会引起计算量的增大，并且会引起舍入误差的大量积累与传播.因此微分方程数值解法也有选择步长的问题. 以经典的四阶龙格-库塔法(5.14)为例.从节点xn出发，先以h为步长求出一个近似值，记为 ，由于局部截断误差为 ，故有,当h值不大时，式中的系数c可近似地看作为常数.,然后将步长折半,即以 为步长,从节点xn出发,跨两步到节点xn+1,再求得一个近似值 ,每跨一步的截断误差是 ,因此有,这样,由此可得,这表明以 作为 的近似值，其误差可用步长折半前后两次计算结果的偏差,来判断所选步长是否适当,当要求的数值精度为时：,（1）如果，反复将步长折半进行计算，直至为止，并以上一次步长的计算结果作为 . 这种通过步长加倍或折半来处理步长的方法称为变步长法.表面上看，为了选择步长，每一步都要反复判断，增加了计算工作量，但在方程的解y(x)变化剧烈的情况下，总的计算工作量得到减少，结果还是合算的.,Euler方法,后退的Euler方法,改进的Euler方法,显式单步法,隐式单步法,显式单步法,一阶常微分方程初值问题,O(h2),O(h2),O(h3),二阶R-K方法,三阶R-K方法,显式单步法,显式单步法,Runge-Kutta方法,O(h3),O(h4),经典的四阶R-K方法,显式单步法,O(h5),5.4 单步法的收敛性及稳定性,5.4.1 单步法的收敛性,常微分方程,差分方程,是否合理？,数值解法的基本思想,定义5.2：若一种数值方法对于任意固定的xn=x0+nh，有 则称该方法是收敛的.,例5.5.就简单的初值问题 考查Euler方法的收敛性.,求解方程的Euler格式为,从而,解：该问题的准确解为：,得,由,收敛！,下面考查一般的单步法.,显式单步法的计算公式均可写为,其中 称为增量函数.,注：不同的显式单步法对应不同的增量函数.,Euler格式,改进的Euler格式,( 5.15 ),单步法的收敛性定理,定理5.1.假设单步法(5.15)满足以下条件,(2)增量函数 关于y满足Lipschitz条件,(1)具有p阶精度,(即局部截断误差为O(hp+1),(3)初值y0是准确的，,即y0 =y(x0),则其整体截断误差,注：从而判断单步法的收敛性，归结为验证增量函数 是否满足Lipschitz条件.,(1)对于Euler方法，由于 故当 满足Lipschitz条件时它是收敛的.,(2)对于改进的Euler方法，,假定f(x,y)关于y满足Lipschitz条件,记其Lipschitz常数为L,则,设限定 上式表明 关于y满足Lipschitz 条件,其Lipschitz常数为,收敛！,(3)类似可以验证Runge-Kutta方法的收敛性.,练习：验证二阶Runge-Kutta方法的收敛性,设限定 上式表明 关于y满足Lipschitz 条件,其Lipschitz常数为,收敛！,5.4.2 单步法的稳定性 稳定性在微分方程的数值解法中是一个非常重要的问题.因为微分方程初值问题的数值方法是用差分格式进行计算的，而在差分方程的求解过程中，存在着各种计算误差，这些计算误差如舍入误差等引起的扰动，在传播过程中，可能会大量积累，对计算结果的准确性将产生影响.这就涉及到算法稳定性问题.,例5.6.考察初值问题 在区间0, 0.5 上的解.分别用Euler方法、后退的Euler方法和改进的欧拉方法计算数值解.,1.0000 2.0000 4.0000 8.0000 16.0000 32.0000,1.0000 2.5000101 6.2500102 1.5625102 3.9063103 9.7656104,1.0000 2.5000 6.2500 15.625 39.0625 97.6563,1.0000 4.9787102 2.4788103 1.2341104 6.1442106 3.0590107,Whats wrong?!,定义5.3.若一种数值方法在节点值yn上产生大小为的扰动,如果在其后的各节点值ym(mn)上产 生的偏差均不超过，则称该方法是稳定的.,稳定性不仅与算法有关，而且与方程中函数f(x,y)有关，讨论起来比较复杂.为简单起见，通常只针对模型方程,为保证微分方程本身的稳定性，假定,一般方程若局部线性化，也可化为上述形式. 模型方程相对比较简单，若一个数值方法对模型方程是稳定的，并不能保证该方法对任何方程都稳定. 但若某方法对模型方程都不稳定，也就很难用于其他方程的求解.,可见,要保证误差是不增长的,只要选取h充分小,使得,Euler方法的稳定性.,模型方程 的Euler公式为,设在节点值yn上有一扰动值 它的传播使节点值yn+1 产生大小为 的扰动值.,假设用 按照Euler格式得出 的计算过程不再有新的误差，则扰动值满足,条件稳定,假设用 按照后