chap常微分方程的数值解法实用

1第7章 常微分方程（组）的数值解法 在科学研究和工程实践中会遇到很多微分方程，虽然从理论上可以证明其解的存在性，但其解的解析表达式往往是很难求解的，或者即使可以写出来，但也难于计算，此时，只能借助数值解来解决问题.常微分方程（组）定解问题是自然科学和工程技术领域中常见的数学模型.本章介绍求解此类问题的基本理论和数值解法。第1页/共55页2定义7.0.1 若存在常数 L 0,使得对一切的 xa,b 及 y,均有 则称 f(x,y)在 D 上关于 y 满足 Lipschitz 条件,其中 L 称为 Lipschitz 常数.我们首先考虑一阶常微分方程初值问题其中 f(x,y)是区域上的实值函数.(7.0.1)第2页/共55页3 我们首先给出常微分方程初值问题解的存在惟一性定理。定理 假设 f(x,y)C(D),且关于 y 满足 Lipschitz 条件,则一阶常微分方程初值问题(7.0.1)存在唯一解.下面在此前提下,我们讨论上述初值问题(7.0.1)的数值解法。第3页/共55页4然后在节点上建立逼近于原初值问题的计算格式(或差分格式),由此计算出原问题的解 y(x)在节点 x1,x2,.,xN 处的近似值：y1,y2,.,yN，称它们为常微分方程初值问题的数值解.相邻两个节点的距离 hn=xn+1-xn 称为步长,通常取定步长h 0,即节点 xn=x0+nh,n=0,1,N.其基本思想是在区间 a,b 上引入一系列节点第4页/共55页57.1 初值问题计算格式的建立1.数值微分方法数值微分方法在等距节点下讨论问题.利用两点数值微分公式7.1.1 计算格式的建立将上式代入初值问题(7.0.1)有(7.1.1)第5页/共55页6略去余项,并以数值解 yn,yn+1 代替 y(xn)及 y(xn+1),则得差分方程差分方程上式称为 Euler 公式.利用此式可由初值 y0 出发按“步进式”方法,逐步求得数值解y1,y2,.,yN.由于计算yn+1时,只用到它前一步的结果yn,这类公式称为单步法.又因为其关于yn+1是显式形式,故称该Euler公式为显格式.第6页/共55页7如果利用下列数值微分公式由 类似的可导出上述公式称为后退的 Euler 公式,此公式为单步法公式.又因为它关于 yn+1 成隐式形式,所以该公式为隐式公式，简称隐格式.第7页/共55页8类似地，可导出上述公式称为Euler两步法公式.这因为，当计算 yn+1 时,要用到 yn-1 与 yn.显然它也是显格式.如果利用下列三点数值微分公式(7.1.3)第8页/共55页9设 y(x)C2a,b,由 Taylor 公式有由于 故上式即为 略去余项,并以 yn,yn+1 代替 y(xn)及 y(xn+1),得到的差分方程正是Euler 公式.(7.1.4)2.Taylor 展开展开法法第9页/共55页103.数值积分方法数值积分方法对 ，在区间 xn,xn+1 上积分，得则有对上式中的积分采用不同的数值积分公式可得到不同的差分方程.例如,对上式的积分采用左矩形公式,可得到 Euler 公式.第10页/共55页11若对此式的积分采用梯形公式,则有若略去余项,以 yn,yn+1 代替 y(xn)及 y(xn+1),得到的差分方程第11页/共55页12上式称为梯形公式.由于它关于 yn+1 成隐式形式,故其为隐格式.隐格式求解比较困难,当 yn 已知时,要求 yn+1，需解关于 yn+1 的非线性方程.在实际应用时,上式常与 Euler 公式联合使用,构成如下计算格式:(7.1.6)第12页/共55页13隐式梯形公式的迭代格式(7.1.7)由上式可以得到一个序列:,k=0,1,关于此序列的收敛性,有如下的定理.第13页/共55页14定理 7.1.1 设 f(x,y)在区域 D 上关于 y 满足 Lipschitz 条件,即其中 L 为 Lipschitz 常数,当步长 时,对任意的初值 按格式（7.1.7）生成的序列 收敛于梯形公式（7.1.6）的解 .第14页/共55页15为了减少计算量,可采用预测-校正格式.方法是先用Euler公式求得一个初始近似值 称为预测值,再把 带入梯形公式右端计算一次求得 yn+1 称之为校正值,即预测:校正:上式称为预测-校正公式或改进的 EulerEuler 公式.上式也可写成如下形式:第15页/共55页16例7.1.1：利用Euler公式与改进的Euler公式求解初值问题（步长h=0.1）解：由步长h=0.1，知节点 设数值解为 利用Euler公式得第16页/共55页17计算结果见下表(见书P227表7.1)此初值问题的解析解为 ,从上表可以看出,数值解 yn与解析解 y(xn)比较,yn精度较差.第17页/共55页18解此问题的改进的Euler公式为 同Euler公式比较,改进的Euler法显然精度提高了.由于误差大小是评价计算格式优劣的重要依据,故需要给出有关误差的概念.计算结果见下表(见书P228表7.2)第18页/共55页197.1.2 截断误差与方法的精度定义 7.1.1 称误差 en+1=y(xn+1)-yn+1为数值方法在点 xn+1 的截断误差,又称整体截断误差.设 yk=y(xk)(k=0,1,.,n)，则 为数值方法在点 xn+1 的局部截断误差.第19页/共55页20整体截断误差 en+1 是在没有引进舍入误差的情况下,纯粹因为不准确的计算格式造成的,故又称为方法误差.它不仅与 x =xn+1 这一步的计算有关,而且和 xn,xn-1,.,x1 这几步的计算都有关系.局部截断误差是假设 xn 之前各数值解没有误差,仅由 xn 到 xn+1 这一步计算由计算格式引起的误差.第20页/共55页21如Euler公式在点 xn+1 的整体截断误差 en+1=y(xn+1)-yn+1局部截断误差第21页/共55页22定义7.1.2 若某数值方法的局部截断误差为 则称该方法具有 P 阶精度,或称其为 P 阶方法.可以证明:Euler 方法的局部截断误差 其具有一阶精度.梯形方法的局部截断误差 其具有二阶精度.改进的 Euler 方法的局部截断误差 具有二阶精度.第22页/共55页237.2 Runge-Kutta 方法 继续讨论前面的 Taylor 展开法。设 y(x)C2a,b,由 Taylor 公式有 由 故上式即为 略去余项,并以 yn,yn+1 代替 y(xn)及 y(xn+1),得到Euler 公式.第23页/共55页24进一步假设设 y(x)Cp+1a,b,由 Taylor 公式有其中 由 故式(7.2.1)即为 第24页/共55页257.2 Runge-Kutta 方法略去余项,并以数值解 yn,yn+1 代替(7.2.3)中的 解析解y(xn)及 y(xn+1),可得到一个差分方程，即其中余项可写成注：这里第25页/共55页26在(7.2.3)中略去余项,用yn,yn+1 代替 y(xn)及 y(xn+1)其中称(7.2.4)式为求解常微分方程初值问题数值解Taylor的格式.第26页/共55页27由于局部截断误差 可知它是一个 p 阶方法。当p=1时，上式正是Euler 公式。但当 p 2 时，需要计算f(x,y(x)的高阶导数，特别是对于复杂函数 f(x,y(x)的求导，这无疑是大大增加计算量，这是它最大的缺点。因此高阶的Taylor方法是不实用的。德国数学家C.Runge 及M.W.Kutta提出了一种改进策略，得到了至今还被作为高精度的单步法广泛使用龙格-库塔法(Runge-Kutta method)。第27页/共55页287.2.1 Runge-Kutta方法的基本思想Runge-Kutta方法是利用 f 在某些点处函数值的线性组合替代（7.2.4）步长 h 后面括号中的因子来构造差分方程,从而避免了高阶导数的计算,这就是 Runge-Kutta 方法的基本思想.用f 在某些点处函数值的线性组合替代这一部分第28页/共55页29其一般形式为:其中 r 是上式中调用 f 的个数,r 称为级数，为待定参数,适当确定这些参数,可使之具有尽可能高的精度.如局部截断误差满足第29页/共55页30 7.2.2 二阶 Runge-Kutta 方法考虑 r=2 的情况,此时有利用二元函数的 一阶Taylor 公式，即全微分公式 希望适当选择参数 使上式的局部截断误差为即为二阶方法.,下面将yn+1与y(xn+1)作比较第30页/共55页31从而有将上式代入再由第31页/共55页32得到在下面要将yn+1与y(xn+1)作比较，使它们的局部截断误差满足为此考虑y(xn+1)。第32页/共55页33再根据 y(xn+1)在点 xn的一元 3 阶 Taylor 展开式由刚才已得到的让它们满足第33页/共55页34即由左式含有四个未知元三个方程,因此解不唯一.参数满足左式的一族公式统称二阶 Runge-Kutta 公式.可得参数应满足下列方程组:第34页/共55页35 取上式称为中点公式.取上式称为 Heun 公式.可见,二阶 Runge-Kutta 公式,每计算一步需要 两次调用 f 的函数值.取得这正是改进的 Euler 公式.第35页/共55页36 7.2.3 四阶Runge-Kutta方法当 r=4 时,类似地可导出四阶 Runge-Kutta 公式,这种公式也有一族,其中常用地有:标准(经典)的 Runge-Kutta 方法第36页/共55页37 Gill 公式 Gill 公式是标准的 Runge-Kutta 公式的改进形式,这种算法可节省存储单元,并能控制舍入误差的增长.四阶 Runge-Kutta 公式,每一步计算需四次调用 f 的函数值,计算量较大,但其局部截断误差可达 O(h5),精度较高.第37页/共55页38例7.2.1 用标准的 四阶Rung-Kutta 法解初值问题,取步长h=0.2.解:解此问题的计算公式为xn0.20.40.60.81.0ynyn-y(xn)1.183 2 1.341 7 1.483 3 1.612 5 1.732 10.000 0 0.000 0 0.000 0 0.000 10.000 1计算结果如下:显然在计算量大致相同的情况下,标准的 Runge-Kutta方法比改进的 Euler 方法精确度更高.（参见p227和p228的表）第38页/共55页397.5一阶常微分方程组与高阶方程初值问题的数值解法 7.5.1 一阶常微分方程组初值问题 第39页/共55页40写成向量形式：其中注意:在形式上(7.5.2)与(7.0.1)一样,所以可以把求解常微分方程初值问题的各种数值方法推广到方程组上来.第40页/共55页41利用向量值函数的微积分理论,很容易推导出一阶常微分方程组初值问题的数值解法.如Euler公式其中 第41页/共55页42(7.5.3)的分量形式为 或 第42页/共55页43四阶标准的Runge-Kutta公式 设 则(7.5.5)的分量形式为 第43页/共55页44四阶标准的Runge-Kutta公式的分量形式 其中第44页/共55页45例7.5.1：试写出用中点公式解下列初值问题的计算公式.解:令 则 再取由向量形式的中点公式第45页/共55页46中点公式的向量形式 上述中点公式的分量计算形式为分量计算形式为第46页/共55页47整理得分量计算格式第47页/共55页487.5.2 高阶常微分方程初值问题的数值解法高阶常微分方程初值问题的一般形式为 引入新变量则(7.5.7)可化为如下一阶常微分方程组初值问题.第48页/共55页49等价的一阶常微分方程组初值问题：于是，可采用前面7.5.1所介绍的方法求此问题的数值解.第49页/共55页50例7.5.2写出用标准四阶Runge-Kutta公式求解的计算公式.解:令 ，将此问题转化为一阶常微分方程组初值问题 再令 则 第50页/共55页51求此问题数值解的标准四阶Runge-Kutta公式为第51页/共55页527.5.3 高阶常微分方程组初值问题的数值解法 对于高阶常微分方程组的初值问题,考察下面的二阶常微分方程组初值问题:令 将此问题转化为一阶常微分方程组初值问题 第52页/共55页53等价的一阶常微分方程组初值问题第53页/共55页54本章小结基本知识建立初值问题数值解的计算格式的方法数值微分方法、Taylor 展开法和数值积分方法截断误差与方法的精度Runge-Kutta方法基本思想一阶、二阶、四阶Runge-Kutta方法一阶常微分方程组数值解法高阶方程（组）的数值解法第54页/共55页7:ODE55感谢您的欣赏！第55页/共55页