第七章-常微分方程的数值解法 1

科学计算与数学建模,中南大学数学科学与计算技术学院,常微分方程的数值解法,第七章常微分方程的数值解法,科学和工程技术中常常要求解常微分方程。根据实际背景不同,所遇到的问题可分为两类:一类是常微分方程初值问题；一类是常微分方程边值问题。一般地,要找出这两类问题的解析解往往非常困难,甚至是不能的。本章将介绍它们的数值解法。所谓数值解法,就是在没有办法知道未知函数的解析表达式的情况下,我们近似计算未知函数在其定义域中的某些离散点上的函数值。当然,如果这些离散点在函数的定义域内的发布很密,且相应点的函数值的计算又非常准确,那么就意味着基本上找到了微分方程的解:微分方程的数值解。,函数是事物的内部联系在数量方面的反映,如何寻找变量之间的函数关系,在实际应用中具有重要意义。在许多实际问题中，往往不能直接找出变量之间的函数关系,但是根据问题所提供的情况，有时可以列出含有要找的函数及其导数的关系式。这就是所谓的微分方程，从而得出微分方程模型。,7.1常微分方程模型的举例,例1物体冷却过程的数学模型,将物体放置于空气中，在时刻时，测量的它的温度为,10分钟后测量得温度为。我们要求此物体的温度和时间的关系,并计算20分钟后物体的温度。这里我们假定空气温度保持为。,解:为了解决上述问题，需要了解有关热力学的一些基本规律。例如，热量总是从温度高的物体向温度低的物体传导的；在一定的温度范围内，一个物体的温度变化速度与这个物体的温度和其所在的介质温度的差值成正比。这是已为实验证实了的牛顿（Newton）冷却定律。设物体在时刻的温度为，则温度的速度以来示。,注意到热量总是从温度高的物体向温度低的物体传导，因而。所以温度差恒正；又因物体将随时间而逐渐冷却,故温度变化速度恒负。故有：(7.1.1)这里是比例常数。方程（7.1.1）就是物体冷却过程的数学模型，它含有未知函数及它的一阶导数,这样的方程称为一阶微分方程。为了解出物体的温度和时间的关系,我们要从方程（7.1.1）中解出。,注意到常数,且,可将(7.1.1)改写成（7.1.2）这样和就被分离开了。两边积分，得到（7.1.3）这里的c是任意常数。上式可写成,令,则有（7.1.4）再根据初始条件可得,于是（7.1.5）如果k的数值确定了,(7.1.5)就完全决定了温度和时间的关系。根据条件时，得到由此得到从而（7.1.6）20分钟后物体的温度就是,从方程（7.1.6）还可得到，当时，,这可解释为：经过一段时间后物体的温度和空气的温度没有什么差别了。微分方程的“解”可以用图形来表示。这往往给我们一个简明直观的了解。从例1中可以大体上看出用微分方程解决实际问题的基本步骤:（1）建立起实际问题的数学模型,也就是建立反映这个实际问题的微分方程；（2）求解这个微分方程；（3）用所得的数学结果解决实际问题,从而预测到某些物理过程的特定性质,以便达到能动改造世界,解决实际问题的目的。在找到了变量之间所要满足的微分方程后,还需要找出代表所考虑的问题的初始状态的条件,这就是所谓的初始条件。求一个微分方程满足一定的初始条件的解的问题，称为微分方程的初值问题。,微分方程初值问题的适定性对于一个微分方程的初值问题，通常要讨论如下三个问题：（1）解的存在性，即初值问题是否有解？（2）解的唯一性，即处置问题的解是否只有一个？（3）解的稳定性，即处置问题的解关于初值、参数等的连续性、可导性等等。以上三个问题也叫做微分方程的适定性。当一个微分方程的定解问题的解是存在、唯一且解关于初值、参数等是稳定的时，就说这个问题是适定的。否则，就说是不适定的。微分方程的初值问题解的适定性，具有重要的实际意义。微分方程模型通常是用来描述确定性的模型。对一个有实际问题所建立的微分方程模型,如果其初值问题的解不存在,或解不唯一,这样的模型本身就不是合理的,是没有实际意义的。因为在一定的条件下物理现象到最后总会有确定的结果,这反映在模型上,就是定解问题有唯一解.而解的稳定性更是具有重要的实际运用背景。,由于在测量初始条件的值和测量方程中各项系数（或参数）等的值时，不可避免地出现测量误差，致使我们得到的微分方程模型，通常只能是近似地描述所讨论的实际问题。自然会问：当测量的数据出现“小”的变动时，相应模型的“解”是否也只有“小”的变动？如果回答是肯定的,我们就说这个模型的解（在某种意义下）是稳定的,否则,就说这个模型是不稳定的。当然，只有“稳定的”解才具有可靠性，只有“稳定的”解才会有使用价值。相反，“不稳定的”解是不会有任何使用价值的。因为初值、参数等的微小误差或干扰将导致“差之毫厘，谬以千里”的严重后果。不过并不是所有的微分方程模型都可用上述方法求解出来。1638年莱布尼茨向全世界提出如下的求解方法.该问题1886年才被数学家刘维尔证明没有解析解。只能借助于本章要介绍的数值解法。,7.2初值问题数值解法的推导方式及常用解法,我们先考虑常微分方程初值问题。它的数学模型是：求，使之满足（7.2.1）其中是已知常数,是已知函数,且满足条件:（7.2.2）（7.2.2）式中的是不依赖的常数,称为李普希兹(Lipschitz)常数。由常微分方程理论可知,上述问题存在唯一解。我们的目标就是计算区间上等距点处该未知函数的函数值其中为此，可以用下列方式设计计算法。,7.2.1数值微分法,在处，微分方程（7.2.1）成为（7.2.3）用数值求导数的公式,分别代替方程（7.2.3）的左边，可得如下由或计算的显示或隐式公式:（7.2.4）,略去含有的误差项并整理,可得如下由或计算的近似计算公式:（7.2.5）其中,分别表示和的近似值。,显然,是已知的函数值,利用公式（7.2.5）就可以逐步计算其它所有节点的近似值。公式（7.2.5）确定的三种方法叫做欧拉法（Euler）、后退欧拉法和中点法。三种方法的不同之处在于：欧拉法可以直接由计算；后退欧拉法在由计算时，则需要求解关于的方程;中点法则必须由两个值来计算。因此，通常称后退欧拉法为隐式方法；称中点法为二步法；称欧拉法为单步的显示方法。,一般地，如果计算时，用到了前m(m1)步的值该方法就称为多步法,否则称为单步法.在推导上述公式时，略去了误差项。我们用表示公式的局部截断误差,之所以称局部截断误差,是因为它是在假设前面若干步计算没有误差,即的条件下近似值与准确值的误差为如果公式的局部截断误差是的同阶无穷小(时),就称这个公式或相应的数值解法是p阶方法.由(7.11)易知,欧拉法、后退欧拉法和中点法分别是1、1、2阶方法。,显然，寻找微分方程的数值解时，计算的每一步都会有局部截断误差，因此的实际截断误差的截断误差都有关系,我们把由此产生的误差称为整体截断误差.例7.2.1分别用欧拉法、后退欧拉法和中点法求解如下初值问题：取步长h=0.02.,解因为步长h=0.02，所以各节点因为利用欧拉法的计算公式可取利用后退欧拉法的计算公式,可解得的显示表达于是按中点法的计算公式,需要知道两个初值.在此，我们利用后退欧拉法计算的结果再依次计算例7.2.1的解析解是用它计算各节点的函数值可得把上述三种方法计算的结果同准确对照,可以看出它们确实都是准确值的近似值,只是误差不一样.欧拉法的误差偏小,后退欧拉法偏大,中点法最小.这跟它们的局部截断误差的符号、阶数、大小是一致的。,7.2.2数值积分,设计初值问题数值解法的第二种方式是数值积分法。它的做法是：对原方程（7.2.1）两边在区间上计算定积分得(7.2.6)再把右边的定积分用数值积分公式计算,从而可以得到关于的方程.最简单的数值积分方法是使用梯形求积公式,即,略去误差项,并代入(7.2.6)式得(7.2.7)其中基于公式(7.2.7)的数值解法称为求解初值问题的梯形法.从上述推导过程可以看出,梯形法的局部截断误差是(7.2.8)从式中用到因此,梯形法是二阶、单步、隐式法。,例7.2.2用梯形法求解如下初值问题：取步长.解因为步长,所以各节点为由梯形法的计算公式:可得所以,根据可依次计算,从例7.2.1和例7.2.2可以看出,利用隐式法求解微分方程都需要求解关于的线性或非线性方程.对线性方程的情形，容易得到的解析表达式，但对线性方程的情形，一般很难得到的解析表达式。因此，只有用第六章介绍的非线性方程的数值解法计算，例如对梯形法来说（参看公式（7.2.7），可构造迭代公式（7.15）其中迭代函数具有压缩性质。这是因为其中L为李普希兹常数，这样，h充分小时，,当然，我们也可以构造的牛顿迭代公式,但不管哪一种迭代公式都需要给定初始值在设计微分方程数值解法时,我们常常先用某个显示法,如欧拉法、中点法等，给出的初始值，又称之为的预测值，再利用迭代法计算下去。实践证明，利用这种方式给出的初始值，迭代法收敛很快。如果求解非线性方程时只迭代一次，那么就得到了一类重要的数值方法，即预测-校正法，它是指出为了计算，先用显示法作预测，再用迭代公式作一次校正。,例7.2.3用欧拉法作预测，用迭代公式（7.2.9）作一次校正,则得到如下预测-校正法:(7.2.10)公式(7.2.10)也可以写成如下显示公式:(7.2.11)上式又称为改进的欧拉公式.,例7.2.4用改进的欧拉法求解如下初值问题:取步长h=0.1.解因为步长h=0.1，所以各节点为因为改进的欧拉法公式可写成,所以根据可依次计算所以例7.2.3的解析解是用它计算各节点处的函数值可得可见改进的欧拉法确实给出了的较好的近似解。,需要指出，利用数值积分设计初值问题的数值解法时，也可以对原问题（7.7）两边在区间上计算定积分，即得（7.2.12）再把右边的定积分用数值积分公式，如辛普森求积公式计算.类似于前面的讨论过程，可以推导出一些初值问题的数值解法及局部截断误差。读者自己可以试一试。,7.2.3泰勒级数法与龙格-库塔法,既然设计微分方程的数值解法就是要推导近似计算的公式，那么利用泰勒公式也是一种好方法。这是因为：（7.2.13）其中是介于与之间的常数，是未知函数在点的k阶导数，但它的值可以利用微分方程本身计算：,（7.2.14）,略去（7.2.14）中的含的项，则得到如下求解初值问题的p阶泰勒级数法：,（7.2.15）,P阶泰勒级数的局部截断误差上述泰勒级数法需要计算的高阶偏导数，计算量较大，仅适用于一些简单的方程。但仿照公式（7.2.15）的形式，可以设计更实用的龙格-库塔（Runge-Kutta）方法,简称R-K方法.它不计算的高阶偏导数，代之以计算在区间上的若干函数值。它的一般形式是：其中和都是待定参数。,参数的选择规则当然是保证所设计出来的方法是尽可能的高阶方法。一般地，如果m给定，则上述龙格-库塔方法至少是m阶方法。下面以二阶龙格-库塔方法的推导为例，说明上述参数的确定方法。二阶龙格-库塔方法的形式如下：（7.2.16）其局部截断误差是（7.2.17）,根据泰勒公式知：,这样，可以写成,为了使二阶龙格-库塔方法成为尽可能高阶的方法，中h和前面的系数应该为零，即要求参数满足条件：又因为的任意性，所以观察的表达式可知，一般不能要求的系数为零。这就是说，m=2时，龙格-库塔方法只能是2阶方法。由于二阶龙格-库塔方法包含四个参数，但只需要满足上述三个条件，所以龙格-库塔方法不止一个求解公式。,理论上说，任选其中一个参数，就能得到一种龙格-库塔公式。容易验证，如果固定则对应的龙格-库塔方法是改进的欧拉法。类似于二阶龙格-库塔方法的推导，可以得到三阶、四阶、五阶、.，龙格-库塔方法。但随着阶数增高，使用公式在区间上要计算的函数值的个数也增多，因而，计算量增大，所以，实际中最常用的龙格-库塔方法是四阶或五阶龙格-库塔方法.,下面列出三个常用公式：,（1）,（2）,（3）,称（1）为标准（或经典）的龙格-库塔公式，称（3）为英格兰（England）公式。为了分析龙格-库塔方法的计算量和计算精度，我们把标准的龙格-库塔方法同欧拉法及改进的欧拉法进行比较。在相同步长的情况下，欧拉法每一步只计算一个函数值，改进的欧拉法每步计算两个函数值，标准龙格-库塔方法每步需计算四个函数值，就是说，标准龙格-库塔方法的计算量差不多是欧拉法的四倍，是改进的欧拉法的二倍。为了比较它们的计算精度，我们可以将欧拉法的步长取为4h。这样，如果用三种方法求解同一个初值问题，则它们的计算量相当。在计算量相当的条件下，比较它们的计算结果，就能够看出它们的精度差异。,例7.2.5分别用欧拉法、改进的欧拉法和标准龙格-库塔方法求解如下初值问题：以比较三种方法的计算精度。解：为了比较三种方法的计算精度，我们将欧拉法的步长取为0.025，改进的欧拉法的步长取为0.05，龙格-库塔方法的步长取为0.1，则三种方法的变量每增加0.1时，都需要计算4个函数值。现将计算结果列于表7.2.1。,表7.2.1,从计算结果可以看出,标准龙格-库塔方法比另外两种方法的精度好很多.在x=0.5处,三种方法的误差分别是,7.3求解初值问题的线性多步法,7.3.1阿达姆斯方法上一节介绍了数值积分法可用来设计初值问题的数值解法.下面用这个方法推导求解初值问题的一类重要方法:阿达姆斯(Admas)方法.首先,对原方程两边在区间上计算定积分,得(7.3.1)假设前r+1(r0)个节点处的函数值已知,即的近似值已算出，从而函数值也已知这样就可以利用r+1组数据点，构造被积函数的r此插值多项式。,其中是拉格朗日插值基函数。利用插值多项式，就可以近似计算（7.3.1）式右端的定积分,从而可得（7.3.2）记则（7.3.2）式可以简单地写成：（7.3.3）式就是求解初值问题的阿达姆斯显式计算公式，它属于多步法。又因为是关于的线性表达式，所以说它是线性多步法。在上述阿达姆斯显式公式的推导中，选用了，作为插值节点，但构造出来的插值多项式是代替区间上的未知函数，因此属于“外插”。,第二章已指出,外插的误差较大。为了弥补这个缺陷,可以将上述推导过程中的节点改为,。这时,公式(7.3.3)相应地变为（7.3.4）注意到（7.3.4）式右端第二项含有（可能是非线性表达式）,所以（7.3.4）式属于隐式公式,称为阿达姆斯隐式公式.对指定的，表7.3.1、表7.3.2分别给出了阿达姆斯显式公式和隐式公式的系数值。,表7.3.1,表7.3.2,表7.3.3,利用插值多项式的余项,可以求出阿达姆斯方法的局部截断误差,对指定的r，表7.3.3列出了它们的局部截断误差主项.,7.3.2米尔尼方法和哈明方法,前面介绍的阿达姆斯显式或隐式公式,以及其它所有公式都可以归纳成如下形式:(7.3.5)当时,(7.3.5)是显式公式;当时,(7.3.5)是隐式公式.因为与之间是线性关系，所以当时，(7.3.5)是线性多步法.下面将用上一节介绍的泰勒级数法构造线性多步法,保证它具有尽量高阶的精度.,首先,公式(7.3.5)的局部截断误差(7.3.6)利用原微分方程得(7.3.7)因为有泰勒展开式,所以(7.3.8)要使线性多步法具有尽量高阶的精度,分析上述的表达可知,只要令,的系数为零,从而得到如下关于,的线性方程组:,(7.3.9)这时,线性多步法的局部截断误差为(7.3.10)值得指出的是，当时，（7.3.8）式要稍加改变，而（7.3.10）式仍然成立.,接下来,我们指定，看一看公式（7.3.5）的几个具体形式。由（7.3.9）式知这时的线性方程组由5个方程,9个未知量,因此有4个自由变量.如果取，则得到容易看出,这时的线性多步法就是的阿达姆斯显式公式，其局部截断误差为（7.3.12）,如果取,则能够导出四阶隐式阿达姆斯公式:(7.35)它的余项是(7.36)如果取则得到这种线性多步法就是著名的四阶显式米尔尼（Milne）公式，它的余项是(7.37),如果则得到著名的四阶隐式哈明（Hmaming）公式：(7.3.15)它的余项(7.3.16)如果则得到隐式辛普森公式:(7.3.17)它的余项是(7.3.18),在求解标准的初值问题（7.2.1）时，多步法需要单步法首先为其提供更多初值，隐式法需要显式法为其提供迭代初值，因此，实际应用中常常将这些方法配合着使用.例如,首先利用四阶的单步法,如标准的龙格-库塔方法为四阶阿达姆斯显式公式提供初值;再利用四阶阿达姆斯显式公式为四阶阿达姆斯隐式公式提供迭代初值.如果四阶阿达姆斯隐式公式只迭代一次,则得到阿达姆斯预测-校正公式:(7.3.19),如果首先利用四阶的单步法，如标准的龙格-库塔方法为米尔尼显式公式提供四个初值，再利用米尔尼显式公式为哈明隐式公式提供的预测值，则得到米尔尼-哈明预测校正公式：（7.3.20）如此配合使用，除了能够克服各种方法单独使用的缺陷以外，还能够及时估计每步的局部截断误差，以米尔尼-哈明预测校正公式为例，根据米尔尼显式公式和哈明隐式公式的局部截断误差可知：,（7.3.21）假定，将以上两式相减得：即,将之代入（7.3.21）式得实用误差估计式：(7.3.22)类似于逐次线性插值和龙贝格求积算法的思想，利用得到的实用误差估计式，我们就可以进一步修正预测值和校正值，从而得到精确度更好的预测校正公式：,(7.3.23)其中是预测值的修正值，为校正值，为的修正值，并作为的近似值，该公式称为米尔尼-哈明预测-校正公式。它仅当时才能用。时，除已知外，都需要用四阶的单步法，如标准的龙格-库塔方法提供，时，=0.此外，每次计算得到的相比，修正的预测-校正公式优点是它可以通过计算较少的函数值达到较高的计算精度，它也可以随时估计局部截断误差。,例7.3.1用修正的米尔尼-哈明预测-校正法求解如下初值问题：步长。要求用标准的方法提供更多初值。解:利用标准的公式：,可依次计算：(1)当时，,(2)当i=1时,(3)当i=2时,接下来，利用修正的米尔尼-哈明预测-校正法公式依次计算,(1)当i=3时,(2)当i=4时,因为是的修正值，其误差限应该不超过的误差限，利用实用误差估计式知，的误差约为所以的误差限为,即准确到了小数点后7位数.,7.4常微分方程边值问题的数值解法,常微分方程边值问题的数学模型：求一个函数，使之满足（7.4.1）当关于是线性函数时，问题（7.4.1）称为线性两点边值问题。常微分方程边值问题的基本数值解法分为两类，一类是将它转化成初值问题来求解，第二类方法是利用数值微商的方法将它转化成线性或非线性方程组。,7.4.1试射法,试射法的基本步骤如下：首先，将边值问题转化成如下形式初值问题，（7.4.2）即依据边值条件寻求与它等价的初始条件：其次，令，从而使问题（7.4.2）转化为如下一阶微分方程组：,最后，令其中，这样方程(7.4.2）就具有形式（7.4.3）显然，（7.4.3）与标准的常微分方程初值问题（7.2.1）具有相同形式，因而所有求解初值问题（7.2.1）的方法都可以用来求解（7.4.3）。所不同的是，只需要把原公式中等分别换成向量和例如，欧拉公式的向量形式是,或用分量形式表示成：因此，利用试射法求解边值问题的关键是如何把边值条件转化为等价的初值条件，即确定。具体方法如下：1.凭经验提供的两个预测值，并按这两个斜率值“试射”。所谓试射，就是按上述试射法的基本步骤分别求解对应的初值问题。假设他们的解为,计算以得到两个近似值；,2.如果均不满足预定的精度，就用线性插值方法校正，即选择新的斜率值3.用试射，又会得到对应的初值问题的解和的近似值。如果不满足预定的精度，回到计算过程2，即利用和选择新的斜率值。重复上述计算过程，直到找到合适的斜率值的近似值，显然在该初值条件下得到的初值问题的解也是原问题（7.4.1）解。,例7.4.1用试射法求解如下问题解:假设则对应一阶微分方程组（初值问题）为：（7.4.5）,选取，用标准龙格库塔方法求解上述方程组，求得；再选取，用标准龙格库塔方法求解上述方程组，求得。由作线性插值，计算m的新的近似值：并由此