常微分方程数值解法的误差分析.doc

淮北师范大学 2013届学士学位论文 常微分方程数值解法的误差分析学院、专业 数学科学学院 数学与应用数学研 究 方 向 计算数学 学 生 姓 名 李 娜 学 号 20091101070 指导教师姓名 陈 昊 指导教师职称 讲 师 年 月 日常微分方程数值解法的误差分析李 娜（淮北师范大学数学科学学院，淮北，235000）摘 要 自然界与工程技术中的很多现象，往往归结为常微分方程定解问题。许多偏微分方程问题也可以化为常微分方程问题来近似求解。因此，研究常微分方程的数值解法是有实际应用意义的。数值解法是一种离散化的数学方法，可以求出函数的精确解在自变量一系列离散点处的近似值。随着计算机计算能力的增强以及数值计算方法的发展，常微分方程的数值求解方法越来越多，比较成熟的有Euler法、后退Euler法、梯形方法、RungeKutta方法、投影法和多步法，等等本文将对这些解的误差进行分析,以求能够得到求解常微分数值解的精度更好的方法。关键词： 常微分方程, 数值解法, 单步法, 线性多步法, 局部截断误差 Error Analysis of Numerical Method for Solving theOrdinary Differential EquationLi Na(School of Mathematical Science, Huaibei Normal University, Huaibei, 235000)AbstractIn nature and engineering have many phenomena , definite solution of the problem often boils down to ordinary differential equations. So study the numerical solution of ordinary differential equations is practical significance. The numerical method is a discrete mathematical methods, and exact solution of the function can be obtained in the approximation of a series of discrete points of the argument.With the enhanced computing power and the development of numerical methods,ordinary differential equations have more and more numerical solution, there are some mature methods. Such as Euler method, backward Euler method, trapezoidal method, Runge-Kutta method, projection method and multi-step method and so on. Therefore, numerical solution of differential equation is of great practical significance. Through this paper, error of these solutions will be analyzed in order to get a the accuracy better way to solve the numerical solution of ordinary differential.Keywords: Ordinary differential equations, numerical solution methods, single step methods, linear multi-step methods, local truncation error 目 录引 言1一、常微分方程11、定义12、常微分方程初值问题描述23、数值解法的基本思想与途径24、数值解的分类35、问题(1)解的存在惟一性定理4二、几种常用的数值解法及其误差分析41、单步法4(一)、欧拉法5(二)、向后EuIer方法6(三)、法7(四)、改进欧拉法7(五)RungeKutta方法92、 线性多步法14总结16参考文献:17引 言自然界中很多事物的运动规律可用微分方程来刻画。常微分方程是研究自然科学和社会科学中的事物、物体和现象运动、演化和变化规律的最为基本的数学理论和方法。物理、化学、生物、工程、航空航天、医学、经济和金融领域中的许多原理和规律都可以描述成适当的常微分方程，因此，常微分方程的理论和方法不仅广泛应用于自然科学，而且越来越多的应用于社会科学的各个领域。由于该问题比较复杂且涉及的面广，使得有些问题的解析解很难求出，而对于一些典型的微分方程(如线性方程、某些特殊的一阶非线性方程等)可以运用基本方法求出其解析解，并在理论上可以根据初值问题的条件把其中的任意常数完全确定下来。然而，在生产实际和科学研究中所遇到的微分方程往往很复杂，在很多情况下都不可能给出解的解析表达式，有时即使能求出形式的解，也往往因计算量太大而不实用，而且高次代数方程求根也并不容易，所以用求解析解的方法来计算微分方程的数值解往往是不适宜的。从实际意义来讲我们更关心的是某些特定的自变量在某一个定义范围内的一系列离散点上的近似值。本文研究的主要是针对常微分方程各种数值解法的误差进行分析。一、常微分方程- 16 -1、定义首先，我们在这部分给出所需的一些基本概念和基本知识。我们已经知道微分方程就是联系着自变量、未知函数以及其导数的关系式。如果在微分方程中，自变量的个数只有一个，我们称这种微分方程为常微分方程。方程 就是常微分方程的例子，这里y是未知函数，t是自变量。微分方程中出现的未知函数最高阶导数的阶数称为微分方程的阶数。2、常微分方程初值问题描述在自然科学和经济的许多领域中，常常会遇到一阶常微分方程的初值问题 (1)这里 是充分光滑，即关于或满足李普希茨条件的二元函数，是给定的初值，称为初始条件。3、数值解法的基本思想与途径一阶微分方程的初值问题(1)的解是区间上的连续变量的函数，因而问题(1)实际上是一个连续性的问题，求这个问题的数值解，就是要求在区间上的若干个离散点处的函数近似值，例如：，然后计算出解的近似值.一般常取为等距离的点，即或 称为步长。建立数值方法的第1步，就是把连续性问题(1)通过一定的方法化为在给定的个点上的近似的差分方程的初值问题，称这个过程为离散化。常用离散化的方法如下：(一)用差商替代导数在点处的导数可以近似地表示成差商从而把初值问题(1)化为差分问题 (2)其中表示解在点处的近似解，即。当然，用差商来近似地表示导数，方法不是唯一的，这里所用的是所谓的向前差商。(二)Taylor展开法在一点(例如点)的附近，的同次数的近似多项式中的Taylor多项式 为最好。其中为一正整数。通过微分方程，便可以逐次把各阶导数在处的值表示出来。(三)数值积分法对微分方程在区间上求积分，得 于是，初值问题(1)便可以近似地化为这样，关于上式右端的积分，可以用数值积分方法计算其近似值。4、数值解的分类常微分方程初值问题的数值解法一般分为两大类：单步法：所谓单步法是指这类方法在计算时，只用到前一步的值然后逐步往下计算。这个算法的代表是龙格-库塔算法，简称RK方法。四阶显示Runge-Kutta方法是求解普通常微分方程初值问题数值解法中的重要方法，而隐式Runge-Kutta公式是求解刚性常微分方程初值问题的重要方法。多步法：这类方法在计算时，除了用到前一步的值之外，还要用到这前面步的值，这个算法的代表就是阿达姆斯(Adams)方法。5、问题(1)解的存在惟一性定理一个常微分方程是不是有特解呢?如果有，又有几个呢?这是微分方程论中一个基本的问题，数学家把它归纳成基本定理，叫做存在和唯一性定理。因为如果没有解，而我们要去求解，那是没有意义的；如果有解而又不是唯一的，那又不好确定。因此，存在和唯一性定理对于微分方程的求解是十分重要的。这个重要的存在和唯一性就是下面列出的著名的存在惟一性定理。定理 如果在带形区域中连续，且关于y满足Lipchiz条件：即存在正常数L，使得 对所有的以及任何都成立，那么初值问题(1)存在惟一的连续可微解。二、几种常用的数值解法及其误差分析 首先来了解下初值问题数值解中最简单的一种方法欧拉法。1、单步法Euler折线法发生的历史背景。在微分方程研究之初，瑞士数学家LEuler(1707.4-1783.9)做出了开创性的工作。他和其他一些数学家在解决力学、物理学问题的过程中创立了微分方程这门学科。在常微分方程方面，Euler在1743年发表的论文中，用代换扩给出了任意阶常系数线性微分方程的古典解法，最早引入了“通解”和“特解的概念。1768年，Euler在其有关月球运行理论的著作中，创立了广泛用于求初值问题(1)的数值解的方法，次年又把它推广到二阶方程。欧拉的想法如下：选择步长，然后在情况下用解函数的切线 代替解函数。这样对于点就可以得到 在点重复如上的程序再次计算新的方向，就会得到所谓的递推公式：，这就是Euler方法。由此，再通过连接所有这些切线得到的函数被称为Euler折线。如果我们令，这些折线就会越来越接近解函数。 (一)、欧拉法Euler方法是最简单的一步法，它是一阶的，精度较差，但公式很简单，即 (3)Euler方法的几何意义在数值计算思想中已经体现出来了，实际上就是用过已知点的折线来近似代替过此点的积分曲线。因此，这种方法又称为折线法。在Euler法中，数值解的误差首先是由差商代替导数引起的，这种近似替代所产生的误差称为截断误差。另外，计算过程中还会由于数值的舍入产生另一种误差舍入误差。显然只有当初产生的误差在以后各步的计算中不会无限制扩大时，即当初始误差充分小时，以后各步的误差也可以充分小，Euler法才具有实用价值。收敛性、截断误差估计与稳定性闷题是常微分方程各种数值解法研究中必须考虑的基本问题。显然这些问题在Euler法中是得到验证的，详见下面例子分析。 在平面上，微分方程(1)的解y=y(x)称作它的积分曲线。积分曲线上一点 线斜率等于函数的值，如果按函数在平面上建立一个方向场，那么，积分曲线上每一点的切线方向均与方向场在该点的方向相一致，基于上述几何解释，从初始点出发，先依方向场在该点的方向推进到上一点，然后再从依方向场的方向推进到上一点，循此前进推出一条折线，一般地，设已做出该折线的顶点，过依方向场的方向再推进到，显然两个顶点，的坐标有关系即 (4)这就是著名的欧拉(Euler)公式。若初值已知，则依公式(4)可逐步算出例1 求解初值问题 (5)解 欧拉公式的具体形式为取步长，计算结果如下表：表1 计算结果对比0.10.20.30.40.51.10001.19181.27741.35821.43511.09541.18321.26491.34161.41420.60.70.80.91.01.50901.58031.64981.71781.78481.48321.54921.61251.67331.7321初值问题(2.1.2)有解，按这个解析式子算出的准确值同近似值一起列在表1， 两者相比较可以看出欧拉方法的精度很差。(二)、向后EuIer方法向后Euler方法和Euler方法差不多，只是把用 去代替，这时计算公式为 （6）向后Euler方法的总体截断误差也是一阶的，因此向后Euler方法是收敛的。这里需要指出它与Euler方法的一个很大不同之处，Euler方法是显式方法，即由明显地表示出来了，而向后Euler方法是隐式方法，计算时要解隐式方程(6)。通常解此方程用迭代法。因此计算较为麻烦，但比显式Euler方法精度要高。(三)、法将Euler方法公式与向后Euler方法公式作加权平均，得到如下公式： （7）称式(7)为初值问题(1)的法公式。其中如为初值条件。在此法中当时，即 （8）此时的法称为梯形公式法。梯形公式也隐式格式，用起来要进行迭代，其计算公式为 （9）这里在应用本迭代法时，是先用Euler方法求初值的近似值即： 然后将替代梯形公式(8)中的得到的式(9)。式(9)又称为预测校正公式。换言之，由Euler方法给出预测值，再用梯形法予以校正。很显然，当步长办取得适当小时，由Euler方法算出的值已是较好的近似。格式(9)收敛很快，通常只需一两次迭代即可满足精度要求，若需多次迭代，则应缩小步长后再行计算。梯形公式法比用向后Euler方法的迭代步长可以放宽一倍，它的总体截断误差为，比Euler方法高一阶。但它每积分一步要计算二次函数值，这说明的精度的提高是以增加计算量为代价的。(四)、改进欧拉法为得到比欧拉法精度高的计算公式，在等式如果对方程(1)从到积分,得 (10)右端积分中若用梯形求积公式近似，并用代替，代替，则得 (11)称为改进欧拉法改进欧拉方法是隐式单步法，可用迭代法求解用欧拉方法提供迭代初值，则改进欧拉法的迭代公式为 (12)为了分析迭代过程的收敛性，将(11)式与(10)相减，得 ，于是有 ，式中为对y满足Lipschitz常数，如果选取h充分小，使得，则当时有，这说明迭代过程(12)是收敛的例2用改进的欧拉方法求解初值问题(1)解改进的欧拉公式为仍取，计算结果见下表同例1中欧拉法的计算结果比较，改进欧拉法明显改善了精度表2 计算结果对比0.10.20.30.40.51.09591.18411.26621.34341.41641.09541.18321.26491.34161.41420.60.70.80.91.01.48601.55251.61531.67821.73791.48321.54921.61251.67331.7321(五)RungeKutta方法德国数学家CDTRunge(18561927)是数值方法发展史上具有里程碑作用的人物。1895年，他在Hanover发表了关于微分方程数值解法的经典论文常微分方程数值解法。此文成为常微分方程RungeKutta方法的发端。此后，Runge结合教学活动积极投身于发展一般的数值分析特别是各种实际应用中的RungeKutta方法(严格来说，此方法在Kutta作出工作后才能称作RungeKutta方法)。RungeKutta方法是一种特殊的单步方法，事实上，这个方法可以看作在上取若干条积分曲线的若干个点的切线斜率，再进行一次(或多次)算术(或加权)平均后产生的新斜率，再按这个斜率从出发，以直线带曲线向前推进一步的过程。与Taylor展示法相比，RungeKutta方法不用增加微商的次数就可以得到较高的阶。Runge一Kutta方法除了在微分方程求解中扮演的传统角色外，人们发现相关类型的初值问题可以用Runge一Kutta方法或适合更一般问题的RungeKutta方法求解，比如RungeKutta方法被应用到了Hamilton系统中。前面提到的几种数值解法的精度是很低的，下面给出高阶一步法RungeKutta方法。它是最常用的一种数值解法，因为它相当精确、稳定、容易编程。RungeKutta方法至今仍然得到广泛地应用。、二级二阶RungeKutta方法由RungeKutta方法的思想，我们得到二阶RungeKutta公式为 （13） (14)RungeKutta公式是在计算两次函数值的情况下，局部截断误差的阶最高是3，式(13)是允许函数任意变化情况下截断误差最小的二阶方法。要再提高阶就必须增加计算函数值的次数。上述式(14)又称为欧拉预估校正公式。、三级三阶RungeKuuta方法两个常用的三阶RungeKuuta方法分别为： （15） （16）这两个常用方法在解决实际问题中能够达到较低的精度要求。但是要更高精度要求的，我们必须了解更高阶的方法一四级四阶RungeKuuta方法。、 四级四阶RungeKutta方法这种方法在解决实际问题中常用。在这里，我们将公式进行导出，以熟悉其方法的实用过程。由于RungeKutta公式对初值问题(1)中的一般都适用，则它必然对特殊的，或也适用。从而可以定出特定的参数来，因此，这里采用的是一种待定系数法，来导出四级四阶显式RungeKutta公式，过程十分简单明了。为了计算简单，令 近似值为，这样求解问题(1)的四阶显示Runge一Kutta公式为 （17） 这里要求 （18）由于RungeKutta公式的思想与数值积分 （19）相似。所以希望四阶显示RungeKutta公式对 （20）是精确的。这样，将（17）式应用于（19）式，就有 与Taylor公式 （21）相比较，有 又由于RungeKutta公式(17)中既要反映的变化，又要反映的变化，所以简单的选取为一对线性函数，这样将(13)式用于 （22）得 （23） 对（22）式求导数有 将上两式代入(23)式并与Taylor公式(21)比较有： 于是总结有如下四阶RungeKutta一般显格式中13个参数所满足的11个参数方程 （24）于是经典的四阶显示RungeKutta公式为 (25)四阶RungeKutta方法中比较常用的显格式为经典的公式(25)，也是根据上述参数方程(24)来确定一般格式中的参数所得到的。因为显格式(17)共有13个参数，而总计11个参数方程，由解的判定定理知，此方程(24)有无穷多个解，从而有四阶RungeKutta方法中的显格式不是唯一的。我们可以尝试推出新算法，但目前为止这个算法是最实用的。RungeKutta方法作为最重要的单步方法，是一类具有相当实用价值的方法。它关于初值是稳定的，其解连续地依赖于初值它是一类便于应用的单步方法，为了计算，只用到前一步的值即可，因此每步的步长可以独立取定，可以按照绝对稳定性、精度等项要求随时更换。常用的RungeKutta方法精度较高，为了达到预定的精度，与Euler方法和梯形法相比，步长办可取得大一些，求解区间上的总步数可以少一些。但RungeKutta方法也有一些缺点，比如四阶RungeKutta方法每算一步需四次计算的值，计算量较大(对于较复杂的而言)。 2、 线性多步法在逐步推进的求解过程中，计算之前事实上已经求出了一系列的近似值，如果充分利用前面多步的信息来预测，则可以期望会获得较高的精度这就是构造所谓线性多步法的基本思想构造多步法的主要途径是基于数值积分方法和基于泰勒展开方法，前者可直接由方程(1)两端积分后利用插值求积公式得到一般的线性多步法公式可表示为 ， (26)其中为的近似，为常数，不全为零，则称(26)为线性步法，计算时需先给出前面个近似值再由（26)逐次求出如果，称(26)为显式步法，这时可直接由(26)算出；如果，则(26)称为隐式步法，求解时与改进欧拉法相同，要用迭代法方可算出，(26)中系数可根据方法的局部截断误差及阶确定，其定义为：设是初值问题(1),(2)的准确解。阿当姆斯显式与隐式公式考虑形如，(27)的步法，称为阿当姆斯(Admas)方法为显式方法，亦称Adams-Bashforth公式；为隐式方法，亦称Adams-Monlton公式，直接由方程(1)两端从到积分求得。例3用四阶阿当姆斯显式和隐式方法解初值问题取步长解本题从四阶阿当姆斯显式公式得到对于四阶阿当姆斯隐式公式得到由此可直接解出而不用迭代，得到计算结果如表3，其中显式方法中的及隐式方法中的均用准确解计算得到，对一般方程，可用四阶R-K方法计算初始近似表3例3计算结果阿当姆斯显式公式阿当姆斯隐式格式0.30.40.50.60.70.