米尔尼方法与辛普森方法

单击此处编辑母版标题样式,单击此处编辑母版文本样式,第二级,第三级,第四级,第五级,*,9.5.3 米尔尼方法与辛普森方法,考虑与(5.7)不同的另一个 的显式公式,其中 为待定常数，可根据使公式的阶尽可,能高这一条件来确定其数值.,由（5.4）可知 ，再令 得,到,1,解此方程组得,于是得到四步显式公式,（5.11）,称为,米尔尼,(Milne）,方法,.,由于 ，故方法为4阶，其局部截断误差为,（5.12）,2,米尔尼方法也可以通过方程（1.1）两端积分,得到.,若将方程（1.1）从 到 积分，可得,右端积分通过辛普森求积公式就有,（5.13）,称为,辛普森方法,. 它是隐式二步四阶方法，其局部截断误,差为,（5.14）,3,9.5.4 汉明方法,辛普森公式是二步方法中阶数最高的，但它的稳定性,较差，为了改善稳定性，考察另一类三步法公式,其中系数 及 为常数.,如果希望导出的公式是四阶的，则系数中至少有一个,自由参数.,若取 ，则可得到辛普森公式.,若取 ，仍利用泰勒展开，由（5.4），令,则可得到,4,解此方程组得,于是有,（5.15）,5,称为,汉明,(Hamming),方法,.,由于 ，故方法是四阶的，且局部截断,误差为,（5.16）,6,9.5.5 预测-校正方法,对于隐式的线性多步法，计算时要进行迭代，计算量,较大.,为了避免进行迭代，通常采用显式公式给出 的,一个初始近似，记为 ，称为,预测,(predictor)，接着计算,的值(evaluation)，再用隐式公式计算 ，称为,校正,（corrector）.,在（2.13）中用欧拉法做预测，再用梯形法校正，得,到改进欧拉法，它就是一个二阶预测-校正方法.,一般情况下，预测公式与校正公式都取同阶的显式方,法与隐式方法相匹配.,例如用四阶的阿当姆斯显式方法做预测，再用四阶阿,当姆斯隐式公式做校正，得到以下格式：,7,预测P：,求值E：,校正C：,求值E：,此公式称为,阿当姆斯四阶预测-校正格式,（PECE).,依据四阶阿当姆斯公式的截断误差，对于PECE的预,测步P有,对校正步C有,8,两式相减得,于是有下列事后误差估计,容易看出,9,（5.17）,比 更好.,但在 的表达式中 是未知的，因此计算时用,上一步代替，从而构造一种修正预测-校正格式 (PMECME ):,P：,M：,E：,10,C：,M：,E：,注意：在PMECME格式中已将（5.17）的 及,分别改为 及 .,利用米尔尼公式（5.11）和汉明公式（5.15）相匹配，,并利用截断误差（5.12），（5.16）改进计算结果，可类,似地建立四阶修正米尔尼-汉明预测-校正格式（,PMECME,）：,11,P：,M：,E：,C：,M：,E：,12,例7,将例6的初值问题用修正的米尔尼-汉明预测-校,正公式计算 及 ，初值 仍用已算出的精,确解，即 ，,给出计算结果及误差.,解,根据修正的米尔尼-汉明预测-校正公式可得,其中,13,14,误差,从结果看，此方法的误差比四阶阿当姆斯隐式法和四阶汉,明方法小，这与理论分析一致.,15,9.5.6 构造多步法公式的注记和例,前面已指出构造多步法公式有基于数值积分和泰勒展,开两种途径，只对能将微分方程（1.1）转化为等价的积分,方程的情形方可利用数值积分方法建立多步法公式，它是,有局限性的.,即前种途径只对部分方法适用.而用泰勒展开则可构造,任意多步法公式，其做法是根据多步法公式的形式，直接,在 处做泰勒展开即可.,不必套用系数公式（5.4）确定多步法（5.1）的系数,及 ，因为多步法公式不一定如（5.1）,的形式.,另外，套用公式容易记错.,16,例8,解初值问题,用显式二步法,其中,试确定参数 使方法阶数尽可能高,并求局部,截断误差.,解,本题仍根据局部截断误差定义，用泰勒展开确定,参数满足的方程.,17,由于,18,为求参数 使方法阶数尽量高，可令,即得方程组,19,解得 ，此时公式为三阶，,而且,即为所求局部截断误差.,而所得二步法为,20,例9,证明存在 的一个值，使线性多步法,是四阶的.,证明,只要证明局部截断误差 ，则方法,为四阶.,仍用泰勒展开，由于,21,22,当 时， ，故方法是四阶的.,23,9.6 方程组和高阶方程,9.6.1 一阶方程组,前面研究了单个方程 的数值解法，只要把,和 理解为向量，那么，所提供的各种计算公式即可应用,到一阶方程组的情形.,考察一阶方程组,的初值问题，初始条件给为,若采用向量的记号，记,24,则上述方程组的初值问题可表示为,（6.1）,求解这一初值问题的四阶龙格-库塔公式为,式中,25,或表示为,其中,26,这里 是第 个因变量 在节点 的近似值.,考察两个方程的特殊情形：,27,这时四阶龙格-库塔公式具有形式,（6.2）,其中,（6.3）,28,（6.3）,这是一步法，利用节点 上的值 ，由（6.3）,式顺序计算 ，然后代入（6.2）,式即可求得节点 上的 .,29,9.6.2 化高阶方程为一阶方程组,高阶微分方程（或方程组）的初值问题，原则上总可,以归结为一阶方程组来求解.,例如，考察下列 阶微分方程,（6.4）,初始条件为,（6.5）,只要引进新的变量,30,即可将 阶方程（6.4）化为如下的一阶方程组：,（6.6）,初始条件（6.5）则相应地化为,（6.7）,初值问题（6.4），（6.5）和（6.6），（6.7）是彼此等,价的.,31,特别地，对于下列二阶方程的初值问题：,引进新的变量 ，即可化为下列一阶方程组的初值,问题：,32,针对这个问题应用四阶龙格-库塔公式（6.2），有,由（6.3）式可得,33,如果消去 ，则上述格式可表示为,这里,34,9.6.3 刚性方程组,在求解方程组（6.1）时，经常出现解的分量数量级差,别很大的情形，这给数值求解带来很大困难，这种问题称,为,刚性,(stiff),问题,.,考察以下例子.,给定系统,（6.8）,它可用解析方法求出准确解，方程右端系数矩阵,35,的特征值为,方程的准确解为,当 时， 称为稳态解，,中均含有快变分量 及慢变分量 .,对应于 的快速衰减的分量在 秒时已衰,36,减到 ，称 为,时间常数,.,当 时快变分量即可被忽略，而对应于 的,慢变分量，它的时间常数 ，它要计,算到 时，才能衰减到 ，也就是说,解 必须计算到 才能达到稳态解.,它表明方程（6.8）的解分量变化速度相差很大，是,一个刚性方程组.,如果用四阶龙格-库塔法求解, 步长选取要满足, 即 ，才能使,计算稳定.,37,而要计算到稳态解至少需要算到 ，则需计算,14 388步.,这种用小步长计算长区间的现象是刚性方程数值求解,出现的困难，它是系统本身病态性质引起的.,对一般的线性系统,（6.9）,其中,若 的特征值,相应的特征向量为 ，则方程组（6.9）的,通解为,38,（6.10）,其中 为任意常数，可由初始条件 确定，,为特解.,假定 的实部 , 则当 时，,为稳态解.,定义8,若线性系统（6.9）中 的特征值 满足条件,，且,则称系统（6.9）为,刚性方程,，称 为,刚性比,.,39,刚性比 时， 为病态矩阵，故刚性方程也称,病态方程.,通常 就认为是刚性的. 越大病态越严重.,方程（6.8）的刚性比 ，故它是刚性的.,对一般非线性方程组（6.1），可类似定义8，将 在,点 处线性展开，记 .,假定 的特征值为 , 于是由定,义8可知，当 满足条件 ,且,40,则称系统（6.1）是刚性的， 称为方程（6.1）的局部,刚性比.,求刚性方程数值解时，若用步长受限制的方法就将出,现小步长计算大区间的问题，因此最好使用对步长 不加,限制的方法，如前面已介绍的欧拉后退法及梯形法，即,A-稳定的方法，,所谓A-稳定就是指数值方法的绝对稳定域包含了,平面的左半平面.,这种方法当然对步长 没有限制，但A-稳定方法要求,太苛刻，Dahlquist已证明所有显式方法都不是A-稳定的，,而隐式的A-稳定多步法阶数最高为2，且以梯形法误差常,数为最小.,这就表明本章所介绍的方法中能用于解刚性方程的方,法很少.,41,通常求解刚性方程的高阶线性多步法是吉尔（Gear),方法，还有隐式龙格-库塔法（见文献21），这些方法,都有现成的数学软件可供使用.,42,