精选数值分析第九章讲义

（优选）数值分析第九章现在是1页\一共有68页\编辑于星期日第九章常微分方程的数值解一、Euler方法三、单步法的收敛性和稳定性二、Runge-Kutta方法四、线性多步法现在是2页\一共有68页\编辑于星期日很多科学技术和工程问题常用常微分方程的形式建立数学模型.但是对于绝大多数的微分方程问题，很难或者根本不可能得到它的解析解.本章重点考察一阶方程的初值问题的数值解法，就是寻求解y(x)在一系列离散点处的近似值的方法.相邻两个节点间的距离称为步长.现在是3页\一共有68页\编辑于星期日一、Euler方法1欧拉公式由初值条件表示积分曲线从出发，并在处的切线斜率为因此可以设想积分曲线在x=x0附近可以用切线近似的代替曲线.切线方程为当x=x1时，代入有这样得到y(x1)的近似值y1的方法.现在是4页\一共有68页\编辑于星期日重复上述方法，当x=x2时依次可以计算出x3,x4,…处的近似值y3,y4,…由此得到Euler公式：由于用折线近似代替方程的解析解，所以Euler方法也称为Euler折线法.例用Euler法计算初值问题的解在x=0.3时的近似值，取步长h=0.1.现在是5页\一共有68页\编辑于星期日解：Euler公式的截断误差局部截断误差：一步Euler公式产生的误差;总体截断误差：Euler公式的累积总误差;现在是6页\一共有68页\编辑于星期日

在假设yn=y(xn)，即第i步计算是精确的前提下，考虑的截断误差Rn=y(xn+1)

yn+1称为局部截断误差.定义欧拉法的局部截断误差：所以欧拉法具有1阶精度.

若某算法的局部截断误差为O(hp+1)，则称该算法有p

阶精度.定义现在是7页\一共有68页\编辑于星期日Lipschitiz条件：若存在正数L,使得对一切x,y1,y2有则称f(x,y)满足Lipschitiz条件.欧拉法的总体截断误差：那么设为局部截断误差，所以现在是8页\一共有68页\编辑于星期日现在是9页\一共有68页\编辑于星期日特别当n=m-1时，有总体误差与h是同阶的.上式还说明，当时，有即也就是说，ym收敛到方程的准确解现在是10页\一共有68页\编辑于星期日后退Euler公式（隐式欧拉法）（隐式欧拉公式）利用向后差商近似导数现在是11页\一共有68页\编辑于星期日由于未知数yn+1同时出现在等式的两边，不能直接得到，故称为隐式欧拉公式，而前者称为显式欧拉公式.一般先用显式计算一个初值，再迭代求解.隐式欧拉法的局部截断误差：即隐式欧拉公式具有1阶精度.现在是12页\一共有68页\编辑于星期日2梯形公式和改进Euler方法梯形公式设y=y(x)是的解,故由此得到用梯形公式近似现在是13页\一共有68页\编辑于星期日用yn来近似y(xn),用yn+1来近似y(xn+1),得梯形公式梯形公式是隐式的，可以用迭代法求解.现在是14页\一共有68页\编辑于星期日具有2阶精度.梯形公式的局部截断误差现在是15页\一共有68页\编辑于星期日中点欧拉公式中心差商近似导数假设,则可以导出即中点公式具有2阶精度.需要2个初值y0和y1来启动递推过程，这样的算法称为双步法，而前面的三种算法都是单步法.现在是16页\一共有68页\编辑于星期日方法显式欧拉隐式欧拉梯形公式中点公式简单精度低稳定性最好精度低,计算量大精度提高计算量大精度提高,显式多一个初值,可能影响精度有没有一种方法，既有这些方法的优点，而没有它们的缺点？现在是17页\一共有68页\编辑于星期日改进欧拉法(1)先用显式欧拉公式作预测，算出(2)再将代入梯形公式的右边作校正,得到注：此法亦称为预测-校正法.可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单.现在是18页\一共有68页\编辑于星期日例用梯形公式求解初值问题(步长h=0.2)解：梯形公式为于是整理得由y(1)=y0=2依次可得y1,y2,y3,y4,y5.现在是19页\一共有68页\编辑于星期日例用改进欧拉法求解初值问题要求步长h=0.2，并计算y(1.2)和y(1.4)解：改进欧拉法公式为即现在是20页\一共有68页\编辑于星期日由y(1)=y0=1计算得现在是21页\一共有68页\编辑于星期日二、Runge-Kutta方法建立高精度的单步递推格式.单步递推法的基本思想是从(xn,yn)点出发，以某一斜率沿直线达到(xn+1

,yn+1

)点.欧拉法及其各种变形所能达到的最高精度为2阶.考察改进的欧拉法，可以将其改写为：斜率一定取K1K2的平均值吗？步长一定是一个h

吗？现在是22页\一共有68页\编辑于星期日首先希望能确定系数1、2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得

将改进欧拉法推广为：(1)

将K2在(xn

,yn

)点作Taylor展开1二阶Runge-Kutta方法现在是23页\一共有68页\编辑于星期日(2)

将K2代入第1式，得到现在是24页\一共有68页\编辑于星期日(3)将yn+1与y(xn+1)在xn点的泰勒展开作比较要求,则必须有：这里有个未知数，个方程。32所以存在无穷多个解！现在是25页\一共有68页\编辑于星期日所有满足上式的统称为2阶Runge-Kutta格式.若则改进的欧拉方法若则中点公式现在是26页\一共有68页\编辑于星期日2四阶Runge-Kutta方法其中i(i=1,…,m)，i(i=2,…,m)

和ij(i=2,…,m;j=1,…,i1

)

均为待定系数，确定这些系数的步骤与前面相似.现在是27页\一共有68页\编辑于星期日由于方程的个数少于未知量的个数，所以方程有无穷多个解，可以根据情况得到几种常用的解，即得到相应的四阶公式.最常用为四阶经典龙格-库塔法也称为标准四阶龙格-库塔公式现在是28页\一共有68页\编辑于星期日Gill公式现在是29页\一共有68页\编辑于星期日753可达到的最高精度642每步须算Ki的个数(2)龙格-库塔法的导出基于泰勒展开，故精度主要受解函数的光滑性影响.对于光滑性不太好的解,最好采用低阶算法而将步长h

取小.注：(1)龙格-库塔法的主要运算在于计算Ki

的值,即计算f的值.Butcher于1965年给出了计算量与可达到的最高精度阶数的关系：现在是30页\一共有68页\编辑于星期日例用标准四阶Runge-Kutta法求初值问题在x=0.1处的近似值，取步长为h=0.1.解：所以现在是31页\一共有68页\编辑于星期日那么例用标准四阶Runge-Kutta法求初值问题在x=0.4处的近似值，取步长为h=0.2.现在是32页\一共有68页\编辑于星期日解：所以而所以现在是33页\一共有68页\编辑于星期日

若某算法对于任意固定的x=xi=x0+ih,当h0

(同时i)时有yi

y(xi

)，则称该算法是收敛的.定义1单步法的收敛性三、单步法的收敛性和稳定性单步法是在计算yn+1时只用到前一步的信息yn

.显式单步法的共同特征是它们都是将yn加上某种形式的增量，得出yn+1,计算公式如下：增量函数现在是34页\一共有68页\编辑于星期日Euler方法的增量函数改进Euler方法的增量函数

设y(x)是微分方程初值问题的准确解，定义则称为显式单步法在xn+1处的局部截断误差.现在是35页\一共有68页\编辑于星期日例：考察欧拉显式格式的收敛性:解：该问题的精确解为

欧拉公式为对任意固定的x=xi=ih，有现在是36页\一共有68页\编辑于星期日

设y(x)是微分方程初值问题的准确解，定义若存在最大整数p，使显式单步法的局部截断误差满足则称该方法具有p阶精度，或称为p方法.Tn+1按h展开的第一项，又称为主项.若局部截断误差的展开式写成则称为局部截断误差的主项现在是37页\一共有68页\编辑于星期日单步法的收敛定理设单步法具有p阶精度其增量函数关于y满足Lipschitz条件即存在常数L，使对任何的及任意的x有又设初值y0是准确的，即则总体截断误差是p阶的，也就是特别的当时,不论n为何值,

总有即方法收敛.现在是38页\一共有68页\编辑于星期日在f(x,y)对y满足Lipschitz条件下,Euler法，改进Euler法和Runge-Kutta法的增量函数

都对y满足Lipschitz条件，所以上述结论对这些方法都成立.例设是求解微分方程的单步法，试求其局部截断误差的主项，并说出它具有几阶精度.解:现在是39页\一共有68页\编辑于星期日考虑在xn处的Taylor展式所以该方法的局部截断误差的主项是具有一阶精度.现在是40页\一共有68页\编辑于星期日例设试求出它具有几阶精度.解:考虑在xn处的Taylor展式现在是41页\一共有68页\编辑于星期日所以该方法的局部截断误差的主项是具有二阶精度.2单步法的稳定性收敛性是在假定每一步计算都准确的前提下,讨论步长时，方法的总体截断误差是否趋于零的问题.稳定性是讨论舍入误差的积累能否对计算结果有严重的影响.现在是42页\一共有68页\编辑于星期日例：考察初值问题在区间0.00.10.20.30.40.5精确解改进欧拉法

欧拉隐式欧拉显式

节点xi

1.00002.00004.00008.00001.6000101

3.2000101

1.00002.5000101

6.25001021.56251023.90631039.76561041.00002.50006.25001.56261013.90631019.76561011.00004.97871022.47881031.23411046.14421063.0590107[0,0.5]上的解，分别用欧拉显、隐式格式和改进的欧拉格式计算数值解.现在是43页\一共有68页\编辑于星期日

若一种数值解法仅在节点值yn上有大小为δ的扰动，于以后各节点值ym(m>n)上，仅由δ所引起的扰动都不超过δ时，称该方法是稳定的.定义一般分析时为简单起见，只考虑试验方程λ

为复数且Re(λ)<0设在节点值yn处有扰动令那么于是现在是44页\一共有68页\编辑于星期日反复应用可得为使则可得如下定义

若中，则称单步法是绝对稳定的.在复平面上，λh满足

的区域称为方法的绝对稳定区域,它与实轴的交称为绝对稳定区间定义下面讨论已知的几种方法的绝对稳定区间和绝对稳定区域.现在是45页\一共有68页\编辑于星期日显式欧拉法:0-1-2ReImg在复平面上的绝对稳定区域是即以-1为中心，1为半径的圆域所以相应的绝对稳定区间是现在是46页\一共有68页\编辑于星期日隐式欧拉法(后退欧拉法)：210ReImg在复平面上的绝对稳定区域是是以1为中心，1为半径的圆的外域所以相应的绝对稳定区间是即如果只考虑λ

<0的实数，则相应的绝对稳定区间对于任意的h都成立，所以是无条件稳定的现在是47页\一共有68页\编辑于星期日梯形公式：在复平面上的绝对稳定区域是是复平面的左半平面即也是无条件稳定的相应的绝对稳定区间是现在是48页\一共有68页\编辑于星期日龙格-库塔法：而显式1~4阶方法的绝对稳定区域为k=1k=2k=3k=4-1-2-3---123ReImg现在是49页\一共有68页\编辑于星期日例设是求解微分方程的单步法，分析它的稳定性.解:所以绝对稳定区域是即为复平面的左半平面.在实数域上是无条件稳定的.现在是50页\一共有68页\编辑于星期日解：将代入得即例讨论求解初值问题的求解公式：的稳定性.(λ>0为实数)所以绝对稳定区域是所以因此是条件稳定的.现在是51页\一共有68页\编辑于星期日四、线性多步法在逐步推进的求解过程中，计算yn+1之前已经求出了一系列的近似值y0,y1,…,yn,如果充分利用前面信息来预测yn+1，则可期望会获得较高的精度，这就是线性多步法的基本思想.1线性多步法的一般公式最常用的线性多步法公式为其中为常数，yn-k为y(xn-k)的近似值fn-k=f(xn-k,yn-k)现在是52页\一共有68页\编辑于星期日特别的当时,上式为显式,否则是隐式.

设y(x)是微分方程初值问题的准确解，定义称为线性多步法在xn+1上的局部截断误差.若则称该方法具有p阶精度.若则称局部截断误差的主项为为误差常数.现在是53页\一共有68页\编辑于星期日例设yn+1=yn-1+2hf(xn,yn)为求解常微分初值问题的线性二步法，试求该二步公式的局部截断误差主项，和精度.解:由局部截断误差的定义可知考虑在xn处的Taylor展式现在是54页\一共有68页\编辑于星期日代入可得所以局部截断误差的主项为具有二阶精度.例试建立求解为微分方程初值问题具有如下形式的线性二步法，并使该方法具有二阶精度，同时求其局部截断误差的主项.解:局部截断误差为现在是55页\一共有68页\编辑于星期日考虑在xn处的Taylor展式于是为使方法具有二阶精度则现在是56页\一共有68页\编辑于星期日解得因此该方法为局部截断误差的主项为例试建立求解为微分方程初值问题具有如下形式的线性二步法，并使该方法具有三阶精度，同时求其局部截断误差的主项.解:局部截断误差为现在是57页\一共有68页\编辑于星期日考虑在xn处的Taylor展式现在是58页\一共有68页\编辑于星期日所以解得局部截断误差的主项为现在是59页\一共有68页\编辑于星期日2Adams外推公式考虑用r+1个点(xn-k,f(xn-k,yn-k))构造一个r次多项式来近似的被积函数f(x,y(x)),这里用yn-k作为y(xn-k)的近似值，令fn-k=f(xn-k,yn-