常微分方程初值问题的数值解法.doc

第八章 常微分方程数值解法常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。第一节 欧拉法求解常微分方程初值问题 （1） 的数值解，就是寻求准确解在一系列离散节点上的近似值 称为问题的数值解，数值解所满足的离散方程统称为差分格式，称为步长，实用中常取定步长。显然，只有当初值问题(1)的解存在且唯一时，使用数值解法才有意义，这一前提条件由下 面定理保证。定理 设函数在区域上 连续，且在区域D内满足李普希兹(Lipschitz)条件，即存在正数L，使得对于R 内任意两点与,恒有则初值问题( 1)的 解存在并且唯一。一 欧拉法(欧拉折线法)若将函数在点处的导数用两点式代替， 即，再用近似地代替,则初值问题(1)变为 (2)(2)式就是著名的欧拉(Euler)公式。以上方法称 为欧 拉法或欧拉折线法。欧拉公式有明显的几何意义。从几何上看，求解初值问题(1)就是平面上求一条通过点的曲线，并使曲线上任意一点处的切线斜率为。 欧拉公式的几何意义就是从点出发作一斜率为的直线交直线于点，点的纵坐标就是的近似值；再从点作一斜率为的直线交直线于点, 点的纵坐标就是的近似值；如此继续进行，得一条折线。该折线就是解的近似图形，如图8-1。 欧拉法的其它几种解释：1 假设在附近展开成泰勒级数取的线性部分，并用作为的近似值，得 2 对方程两边从到积分，得 (3) 用 矩形公式计算上式右侧积分，即 并用作为的近似值，得故欧拉法也称为矩形法。二 改进的欧拉法(梯形法) 欧拉法形式简单，低，特别当曲线y=y(x)计算方便，但精度比较的曲率较大时，欧拉法的效果更差。为了达到较高精度的计算公式，对欧拉法进行改进，用梯形公式计算(3)式 右侧积分，即并用y作为y(x)的近似值，得到改进 的欧拉公式 (4)上述方法称为改进的欧拉法，也称为梯形法。不难发现，欧拉公式 是 关于y的显式，即只要已知y,经过一次计算便可得y的值，而改进的 欧拉公式是以y的隐式方程给出，不能直接得到y。隐式方程(4)通常用 迭代法求解，而迭代过程的实质是逐步显式化。先用欧拉 公式给出y的迭代初值，然后 再用改进的欧拉公式(4)进行迭代，即有 (5)迭代过程进行到连续两次迭代结果之差的绝对值小于给定的精度即为止，这时取然后再转入下一步计算。下面讨论是否收敛；若收敛，它的极限是否满足(4)式。假设f(x,y)满足李普希兹条件 则由此可见，只要(这里只要步长h足够小即可)，当k时，有，所以收敛。又因为f(x ,y)对y连续，当k时，对等式两端取极限， 得因 此，只要步长h足够小，就可保证收敛到满足(4)式的。三 预估一校正法改进的欧拉公式在实际计算时要进行多次 迭代，因而计算量较大。在实用上，对于改进的欧拉公式(5)只迭代一次，即先用欧拉公式算出y的预估值y，再用改进的欧拉公式(4)进行一次迭代得到 校正值y，即 (6)预估校正公式也常写成下列形式： (7)四 公式的截断误差定义 若某种微分方程数值解 公式的截断误差是O(h)，则称这种方法是k阶方法。为了简化分析，在进行误差分析时，我们假设前一步的结果是准确的，即在 y=y(x)的前提下，考虑用y作为y(x)的近似值而产生的截断误差，这种误差称为局部截断误差。由泰勒公式对于欧拉公式，有 于是 则欧拉公式的截断误差为O(),所以 欧拉法是一阶方法。对于预估校正公式，有而 于是因此所以y(x)- y= O()则预估校正公式的截断误差为O(),也即预估校正法是二阶方法。 可以证明，改进的欧拉公式与预估校正公式的截断误差相同，均为O()。这里略去证明。例1 在区间0,1上以h=0.1为步长，分别用欧拉法与预估 校正法求初值问题的数值解解 该方程为贝努利方程，其精确解为欧拉公式的具体形式为预估校正公式的具体形式为取步长,计算结果见下表：表8-1表8-1 欧拉法预估-矫正法准确解近似解与准确解比较，欧拉法的结果大致只有两位有效数字，而预估校正法的结果则 有3位有效数字。第二节 龙格库塔法前面讨论的欧拉法与改进的欧拉法都是一步 法，即计算y时，只用到前一步值。龙格库塔(Runge-Kut ta)法(简称为R-K方法)是一类高精度的一步法，这类方法与泰勒级数法有着密 切的关系。一 泰勒级数法设有初值问题由 泰勒展开式 （1） 若令 (2)则即公式(2)为k阶方法。 从理论上讲，只要解y(x)有任意阶导数，泰勒展开方法就可以构造任意阶求y的公 式，但由于计算这些导数是非常复杂的。如所以这种方法实际上不能用来解初值问题。二 龙格 库塔方法(R-K方法)R-K方法不是通过求导数的方法构造近似公式，而是通过计算 不同点上的函数值，并对这些函数值作线性组合，构造近似公式，再把近似公式与解的泰勒展开 式进行比较，使前面的若干项相同，从而使近似公式达到一定的阶数。我们先分析欧拉法 与预估校正法。对于欧拉法 每步计算f的值一次，其截断误差为O()。对于预估校正法。 每步计算的值两次，其截断误差为O()。下面对预估校正法进行改进，将该公式写成更一般的形式 (3)其中为待定常数。选择这 些常数的原则是在的前提下，使)的阶尽量 高。为此，作泰勒展开 其中都是在处的函数值。将代入得与泰勒展开式(2)进行比较，要使得只要 四个参数满足 (4)若,即得预估校正公式。满足( 4)式的可以有各种不同的取法，但不管如何取法，都要计算两次的值(即 计算f在两个不同点的函数值)，截断误差都是。满足条件(4)的一族公式(3)统称为 二阶龙格库塔公式。人们容易想到，如果不增加计算函数值的次数，能否适当地选择这 四个参数，使近似公式的局部截断误差的阶再提高，比如达到。为此，把多展 开 一项，有 （4）所以 而在的泰勒展开式由于展开式的项中是无法通过选择参数来消去的，所以不论四个参 数如何选择，都不能使局部截断误差达到。要想提 高近似公式的阶，只能继续增加计算f的值的次数。如果每步计算三次的值，可将公式写 成下列形式： (5)与二阶龙格库塔公式的讨论方法类似，要使，只需8个参数满足方程 组包含6个方程，8个未知量，其解不唯一。满足条件(6)的一族公式(5)统称为三阶龙格库 塔公式。一个比较简单的三阶龙格龙塔公式是 截断误差为的四阶龙格库塔公式 是常用的公式，每步都要计算四次f的值。它的一般形式是 (6)(6)式中13个待定常数需满足下列11个方程的方程组最常用的四阶龙格库 塔公式是标准四阶 龙格库塔公式：和 吉尔公式标准四阶龙格库塔公式手算时采用表8-2所示的表格计算。例 用标准四阶龙格库塔方法求解初值问题解 计算过程和结果 如表8-3所示 因此有 表8-3对该例，用几种不同的一步法计算的结果如表8-4。由表8-4可见，虽然四阶龙格库塔方法每步要计算四次f的值，但以=0.2为 步长的计算结果就有5位有效数字，而欧拉法与预估校正法以=0.1为步长的计算结果 才具有2位与3位有效数字。如果步长也取0.2，则结果的精度会更低。表8-4欧拉法预估-校正法四阶法准确解第三节亚当斯方法求解初值问题的一步法在计算时只用到前面一步的结果，所以要提高精度时，需要增加中间函数值的计算，这就加大了计算量。如果在计算yn+1时，不仅用到xn上的近似值yn，还用到前面若干节点xn-1 ,xn-2 ，上的近似值yn-1 ,yn-2 ,，称这种方法为多步法。亚当斯(Adams)方法是多步法中的一种。我们知道，求解初值问题 (8.3.1)等价于求解积分方程 (8.3.2)选用不同的数值方法计算(8.3.2)式右端的积分项，就会导出不同的计算公式。例如，用矩形法计算积分项代入(8.3.2)式得离散化即得欧拉公式，其截断误差为O(h2)。为了提高精度，改用梯形法计算积分项代入(8.3.2)式得 离散化得到改进的欧拉公式，其截断误差是O(h3)。为此，我们想到，基于插值原理可以建立一系列的数值积分方法，运用这些方法可以导出求初值问题的一系列计算公式。一般地，若用插值多项式k(x)，代替f(x,y(x),用作为的近似值，离散化后得公式 (8.3.3)8.3.1亚当斯显式公式设初值问题(8.3.1)的解y(x)在xn ,xn-1 ,xn-k上各点的近似值yn ,yn-1 ,yn-k都已计算出来，用点(xn ,fn),(xn-1 ,fn-1),，(xn-k ,fn-k)作后插公式其中x=xn+th，将k(x)代入(8.3.3)式，得亚当斯显式公式 其中它的前几个值见表8-5 表8-5j 0 1 2 3 aj 1 1/2 5/12 3/8当k=0时，(8.3.4)式为欧拉公式。当k=1时，(8.3.4)式为 当 k=3时，(8.3.4)式为 (8.3.5)利用差分与函数值之间的关系，(8.3.5)可写成下列形式 (8.3.6)假设(8.3.6) 式右端的值精确，即yn-i=y(xn-i)，i=0,1,2,3.那么fn-i=f(xn-i ,yn-i)=f(xn-i ,y(xn-i)=y(xn-i)代入(8.3.6)式 将上式右端各项在点xn处展开，得另一方 面，对于准确解y(xn+1)，有于是所以公式(8.3.6)的局部截断误差为O(h5)，它是四阶方法，称(8.3.6)式为四阶亚当斯显式公式。8.3.2亚当斯隐式公式亚当斯显式方法选取xn ,xn-1 ,xn-k 作为插值点，这时，用插值函数k(x)在求积区间xn ,xn+1 上逼近函数f(x,y(x)实际上是一个外推过程，效果不太理想。为了提高精度，我们变外推为内插，改用xn+1 ,xn ,xn-k+1为插值节点建立后插公式k(x),然后重复前一段的推导过程，得到亚当斯隐式公式其中 它的前几个值见表8-6。 表8-6 j 0 1 2 3 bj 1 -1/2 -1/12 -1/24当k=0时，(8.3.7)式为 当k=1时，(8.3.7)式为梯形公式当 k=3时，(8.3.7)式为 即 (8.3.8)类似于亚当斯显式公式，可以导出隐式公式(8.3.8)的局部截断误差为所以隐式公式(8.3.8)是四阶方法。8.3.3亚当斯预估一校正公式 把亚当斯显式公式与隐式公式联立使用，构造亚当斯预估一校正公式。以四阶亚当斯方法为例有下列公式： (8.3.9)迭代过程进行到连续 两次迭代结果之差的绝对值小于给定的精度为止。这时取，然后转入下一步计算。只迭代1次的公式 (8.3.10)称为亚当斯预估一校正公式。第1式称为预估公式，第2式称为校正公式。其局部截断误差为下面我们具体讨论误差估计问题。设 由(8.3.10)式求得准确解 y(xn+1)的近似解为yn+1，由于上面二式相减得 所以于是有若 则可持续求解yn+2,，否则,应缩小步长h重新计算。与同阶的龙格库塔方法相比较，亚当斯方法计算量小，公式简单，程序易于实现。但是由于亚当斯方法在计算yn+1时，不仅用到前面一步的信息，而且还要用到更前面几步的信息，因此它不能自动开始。开始的前几个值必须借助于一步法获得。例 分别用四阶亚当斯显式公式与预估一校正公式求解初值问题 解用第二节例1，按标准四阶龙格库塔方法求得的结果y1,y2,y3作为开始值，然后用显式公式(8.3.6)与预估一校正方法(8.3.10) 进行计算，计算结果如表8-7所示。表8-7 xn ynR-K方法 显式方法预估计-校正法 准确解00.10.20.30.40.50.60.70.80.91.001.095445531.183216751.264912231.341551761.414046421.483018191.548918881.612116341.672917041.731569761.3416411361.414213831.483239831.549193381.612451541.673320001.7320507211.095445121.183215961.264911061.341640791.414213561.483239701.549193341.612351551.673320051.73205081近似解与准确解进行比较，显式方法的结果有4位有效数字，而预估一校正法的结果则有7位有效数字。 第四节线性多步法线性多步法的一般公式为= +h() (8.4.1) 当-1 =0时，公式(8.4.1)是显式的；当-1 0时 ，公式(8.4.1)是隐式的。“线性”是指是, , , ,的线性组合，“多步”是指计算yn+1时 ，不仅用到前一步的结果yn,而且用到更前面几步的结果。所有形如(8.4.1)式的公式都 可以利用泰勒展开的方法构造出来。作为例子，我们推导形如 = +h() (8.4.2)的数值计算公式。假设 =y(), i=n, n-1, n-2 =f(,)i=n,n-1,n-2,n-3 将上述各项在xn处泰勒展开 Yn-i=y(xn-ih)=y(xn)-ih+-+, i=1,2= =-+-+, i=1,2,3代入(8.4.2)式 得=()y()+()+()+()+()+. (8.4.3)将y()在处泰勒展开(8.4.4)比较(8.4.3)式与(8.4.4)式，让h的同次幂的系数相等，得 (8.4.5) 如果(8.4.5)的前7个方程求出未知数,则其局部截断误差为.当然，也可以从前k个方程解出未知数(k7)，其局部截断误差为.比如从前5 个方程解出7个未知数，因此，有2个自由未知数，若令,则从(8.4.5)的前5 个方程解得 对应的公式是它正好是四阶阿达姆斯显式公式。类似地， 若设计算公式为= +h() 为使这类公式为四阶方法，系数应满足 (8.4.7) 特别地，取，则从(8.4.7)式解得 对应的公式是它正好是四 阶亚当斯隐式公式。第五节方程组与高阶方程的数值解法前面几节介绍了一阶微分方程初值问题的各种数值解法，这些解法同样适用 于一阶微分方程组与高阶方程。为了避免书写的复杂，下面仅就两个未知函数的方程组和二阶方程为例说明各种方法的计算公式。8.5.1一阶方程组设有一阶微分方程组初值问题 (8.5.1)1欧拉法的计算公式 (8.5.2)2改进的欧拉法的计算公式对n=0,1,2,，计算 (8.5.3) 当连续两次迭代结果之差的绝对值小于给定的精度,即 时，取,，然后转入下一步运算。3四阶标准龙格库塔公式 (8.5.4)4四阶亚当斯显式公式 (8.5.5) 其中 需要注意，前几步的值由其它单步法求得。8.5.2二阶方程设有二阶微分方程初值问题 (8.5.6)令z=y，则(8.5.6)化为一阶微分方程初值问题 (8.5.7)令(8.5.7)式中,于是应用四阶龙格库塔公式(8.5.4),有 则又因为 所以于是得初值问题(8.5.6)的计算公式为 第六节 边值问题的数值解法在具体求解微分方程时，需要附加某种定解条件。微分方程与定解条件一起组成定解问题。对于高阶微分方程，定解条件通常有两种给法，第一种是给出积分 曲线在初始时刻的性态，称这类条件为初始条件，相对应的定解问题称为初值问题；第二种是给出积分曲线在首末两端的性态，称这类条件为边值条件，相对应的定解问题为边值问题 。本节以二阶微分方程y=f(x,y,y),axb为例讨论边值问题，其边值条件 可分为三类：第一边值条件y(a)=a, y(b)=第二边值条件 y(a)=a, y(b)= 第三边值条件y(a)-a0y(a)=a , y(b)- 0y(b)= 其中 a00,00,a0 +00.本节介绍解线性方程边值问题的差分方法以及 适用于非线性方程边值问题的试射法。8.6.1解线性方程边值问题的差分方法二阶线性方程的一般形式为 y+p(x)y+q(x)y=f(x), axb (8.6.1)首先将区间a,b进行等距划分，分点 , i=0,1,2,n其中一般地称其次，在各个节点上，将y,y用差商近似表示。这里要求有相同阶数的截断误差，以保证精度协调。对于内部节点，二阶导数用二阶中心差商表示，得 i=1,2,n-1一 阶导数用一阶中心差商表示，得 i=1,2,n-1 假设，则于是得方程(8.6.1)的差分方程 + i=1,2,n-1 (8.6.2)其中将(8.6.2)整理，可以写成下列形式 i=1,2,n-1 (8.6.3)其中 (8.6.3)是含有n+1个未知数yi (i=0,1,2,n)的线性方程组，方程的个数为n-1。要使方程组(8.6.3)有唯一解，还需要有两上边值条件补充两个方程。对于第一边值条件，直接就得到两个方程 于是得到第一边值问题的差分方程组 (8.6.4)这个方程组是三对角方程组，可以利用追赶法求解。 对于第二及第三 边值条件，由于条件中包含了导数，所以边值条件也必须用差商来近似表示。为使截断误差达到，需要用到牛顿等距插值公式，得 将这两个近似公式代入边值条件中，就得到两个方程，再与(8.6.3)联立，就可得到对应的差分方程组， 用追赶法解出yi(i=0,1,n)。对于边值问题的收敛性，即考虑当h0时，差分 方