《数值分析》第九章常微分方程的数值解同步课程

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

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

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

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

阶精度.定义第七页，共69页。Lipschitiz条件：若存在正数L,使得对一切x,y1,y2有则称f(x,y)满足Lipschitiz条件.欧拉法的总体截断误差：那么设为局部截断误差，所以第八页，共69页。第九页，共69页。特别当n=m-1时，有总体误差与h是同阶的.上式还说明，当时，有即也就是说，ym收敛到方程的准确解第十页，共69页。后退Euler公式（隐式欧拉法）（隐式欧拉公式）利用向后差商近似导数第十一页，共69页。由于未知数yn+1同时出现在等式的两边，不能直接得到，故称为隐式欧拉公式，而前者称为显式欧拉公式.一般先用显式计算一个初值，再迭代求解.隐式欧拉法的局部截断误差：即隐式欧拉公式具有1阶精度.第十二页，共69页。2梯形公式和改进Euler方法梯形公式设y=y(x)是的解,故由此得到用梯形公式近似第十三页，共69页。用yn来近似y(xn),用yn+1来近似y(xn+1),得梯形公式梯形公式是隐式的，可以用迭代法求解.第十四页，共69页。具有2阶精度.梯形公式的局部截断误差第十五页，共69页。中点欧拉公式中心差商近似导数假设,则可以导出即中点公式具有2阶精度.需要2个初值y0和y1来启动递推过程，这样的算法称为双步法，而前面的三种算法都是单步法.第十六页，共69页。方法

显式欧拉隐式欧拉梯形公式中点公式简单精度低稳定性最好精度低,计算量大精度提高计算量大精度提高,显式多一个初值,可能影响精度有没有一种方法，既有这些方法的优点，而没有它们的缺点？第十七页，共69页。改进欧拉法(1)先用显式欧拉公式作预测，算出(2)再将代入梯形公式的右边作校正,得到注：此法亦称为预测-校正法.可以证明该算法具有2阶精度，同时可以看到它是个单步递推格式，比隐式公式的迭代求解过程简单.第十八页，共69页。例用梯形公式求解初值问题(步长h=0.2)解：梯形公式为于是整理得由y(1)=y0=2依次可得y1,y2,y3,y4,y5.第十九页，共69页。例用改进欧拉法求解初值问题要求步长h=0.2，并计算y(1.2)和y(1.4)解：改进欧拉法公式为即第二十页，共69页。由y(1)=y0=1计算得第二十一页，共69页。二、Runge-Kutta方法建立高精度的单步递推格式.单步递推法的基本思想是从(xn,yn)点出发，以某一斜率沿直线达到(xn+1

,yn+1

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

吗？第二十二页，共69页。首先希望能确定系数

1、

2、p，使得到的算法格式有2阶精度，即在的前提假设下，使得

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

将K2在(xn

,yn

)点作Taylor展开1二阶Runge-Kutta方法第二十三页，共69页。(2)

将K2代入第1式，得到第二十四页，共69页。(3)将yn+1与y(xn+1)在xn点的泰勒展开作比较要求,则必须有：这里有个未知数，个方程。32所以存在无穷多个解！第二十五页，共69页。所有满足上式的统称为2阶Runge-Kutta格式.若则改进的欧拉方法若则中点公式第二十六页，共69页。2四阶Runge-Kutta方法其中

i(i=1,…,m)，

i(i=2,…,m)

和

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

)

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

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

的值,即计算f的值.Butcher于1965年给出了计算量与可达到的最高精度阶数的关系：第三十页，共69页。例用标准四阶Runge-Kutta法求初值问题在x=0.1处的近似值，取步长为h=0.1.解：所以第三十一页，共69页。那么例用标准四阶Runge-Kutta法求初值问题在x=0.4处的近似值，取步长为h=0.2.第三十二页，共69页。解：所以而所以第三十三页，共69页。

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

(同时i

)时有yi

y(xi

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

.显式单步法的共同特征是它们都是将yn加上某种形式的增量，得出yn+1,计算公式如下：增量函数第三十四页，共69页。Euler方法的增量函数改进Euler方法的增量函数

设y(x)是微分方程初值问题的准确解，定义则称为显式单步法在xn+1处的局部截断误差.第三十五页，共69页。例：考察欧拉显式格式的收敛性:解：该问题的精确解为

欧拉公式为对任意固定的x=xi=ih，有

第三十六页，共69页。

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

总有即方法收敛.第三十八页，共69页。在f(x,y)对y满足Lipschitz条件下,Euler法，改进Euler法和Runge-Kutta法的增量函数

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

欧拉隐式欧拉显式

节点xi

1.0000

2.00004.0000

8.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]上的解，分别用欧拉显、隐式格式和改进的欧拉格式计算数值解.第四十三页，共69页。

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

为复数且Re(λ)<0设在节点值yn处有扰动令那么于是第四十四页，共69页。反复应用可得为使则可得如下定义

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

的区域称为方法的绝对稳定区域,它与实轴的交称为绝对稳定区间定义下面讨论已知的几种方法的绝对稳定区间和绝对稳定区域.第四十五页，共69页。显式欧拉法:0-1-2ReImg在复平面上的绝对稳定区域是即以-1为中心，1为半径的圆域所以相应的绝对稳定区间是第四十六页，共69页。隐式欧拉法(后退欧拉法)：210ReImg在复平面上的绝对稳定区域是是以1为中心，1为半径的圆的外域所以相应的绝对稳定区间是即如果只考虑λ

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

设y(x)是微分方程初值问题的准确解，定义称为线性多步法在xn+1上的局部截断误差.若则称该方法具有p阶精度.若则称局部截断误差的主项为为误差常数.第五十三页，共69页。例设yn+1=yn-1+2hf(xn,yn)为求解常微分初值问题的线性二步法，试求该二步公式的局部截断误差主项，和精度.解:由局部截断误差的定义可知考虑在xn处的Taylor展式第五十四页，共69页。代入可得所以局部截断误差的主项为具有二阶精度.例试建立求解为微分方程初值问题具有如下形式的线性二步法，并使该方法具有二阶精度，同时求其局部截断误差的主项.解:局部截断误差为第五十五页，共69页。考虑在xn处的Taylor展式于是为使方法具有二阶精度则第五十六页，共69页。解得因此该方法为局部截断误差的主项为例试建立求解为微分方程初值问题具有如下形式的线性二步法，并使该方法具有三阶精度，同时求其局部截断误差的主项.解:局部截断误差为第五十七页，共69页。考虑在xn处的Taylor展式第五十八页，共69页。所以解得局部截断误差的主项为第五十九页，共69页。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-