《常微分方程数值解法(精品)》由会员分享，可在线阅读，更多相关《常微分方程数值解法(精品)（16页珍藏版）》请在金锄头文库上搜索。
1、第九章 常微分方程的数值解法在自然科学的许多领域中，都会遇到常微分方程的求解问题。然而，我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。这些方法可以给出解的近似表达式，通常称为近似解析方法。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。利用计算机解微分方程主要使用数值方法。我们考虑一阶常微分方程初值问题在区间a, b上的解，其中f (x, y)为x, y的已知函数，y0为给定的初始值，将上述问题的精确解记为y(x)。数值方法的基本思想是：在解的存在区间上取n + 1个
2、节点这里差，i = 0,1, , n称为由xi到xi+1的步长。这些hi可以不相等，但一般取成相等的，这时。在这些节点上采用离散化方法，（通常用数值积分、微分。泰勒展开等）将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解yn作为y(xn)的近似值。这样求得的yn就是上述初值问题在节点xn上的数值解。一般说来，不同的离散化导致不同的方法。1 欧拉法与改进欧拉法1.欧拉法1对常微分方程初始问题用数值方法求解时，我们总是认为(9.1)、(9.2)的解存在且唯一。欧拉法是解初值问题的最简单的数值方法。从（9.2）式由于y (x0) = y0已给定，因而可以算出设x1 = h充分小，则近似地
3、有：(9.3)记 从而我们可以取作为y (x1)的近似值。利用y1及f (x1, y1)又可以算出y(x2)的近似值：一般地，在任意点xn+1 = (n + 1)h处y(x)的近似值由下式给出（9.4）这就是欧拉法的计算公式，h称为步长。不难看出，近似解的误差首先是由差商近似代替微商（见(9.3)）引起的，这种近似代替所产生的误差称为截断误差。还有一种误差称为舍入误差，这种误差是由于利用（9.4）进行计算时数值舍入引起的。例9.1 用欧拉法求初值问题当h = 0.02时在区间0, 0.10上的数值解。解 把代入欧拉法计算公式。就得具体计算结果如下表：nxnyny(xn)en = y(xn) -
4、 yn001.00001.0000010.020.98200.98250.000520.040.96500.96600.000530.060.94890.95030.001440.080.93360.93540.001850.100.91920.9230.0021在上表中y(xn)列，乃是初值问题（9.5）、（9.6）的真解在xn上的值。为近似值yn的误差。从表中可以看出，随着n的增大，误差也在增大，所以说，欧拉法计算简便，对一些问题有较大的使用价值，但是，它的误差较大，所得的数值解精确度不高。2改进欧拉法为了构造比较精确的数值方法，我们从另一角度重新分析一下初值问题。一般说来，一阶方程的初值
5、问题与积分方程（9.7）是等价的，当x = x1时，（9.8）要得到y(x1)的值，就必须计算出（9.8）式右端的积分。但积分式中含有未知函数，无法直接计算,只好借助于数值积分。假如用矩形法进行数值积分则因此有可见，用矩形法计算右端的积分与用欧拉法计出的结果完全相同。因此也可以说欧拉法的精度之所以很低是由于采用矩形法计算右端积分的结果。可以想象，用梯形公式计算（9.8）式右端的积分，可期望得到较高的精度。这时将这个结果代入（9.3）并将其中的y(x1)用y1近似代替，则得这里得到了一个含有y1的方程式，如果能从中解出y1，用它作为y(x1)的近似值，可以认为比用欧拉法得出的结果要好些。仿照求y
6、1的方法，可以逐个地求出y2, y3，。一般地当求出yn以后，要求yn+1，则可归结为解方程：这个方法称为梯形法则。用梯形法则求解，需要解含有yn+1的方程式，这常常很不容易。为此，在实际计算时，可将欧拉法与梯形法则相结合，计算公式为（9.9）这就是先用欧拉法由(xn, yn)得出y(xn+1)的初始近似值，然后用（9.9）中第二式进行迭代，反复改进这个近似值，直到（e为所允许的误差）为止，并把取作y(xn+1)的近似值yn+1。这个方法就叫改进欧拉方法。显然，应用改进欧拉法，如果序列收敛，它的极限便满足方程即序列的极限可取作yn+1。可以证明，如果有界，则只要h取得适当小，上述序列必定收敛。
7、这样当h取得充分小，就可保证上述迭代过程收敛到一个解。当步长h取得适当时，欧拉方法算出的值已是较好的近似，因此改进欧拉法收敛很快，通常只需二、三次迭代即可。如果迭代很多步仍不收敛，这表明表长h选的过大，应缩小步长后再计算。通常把（9.9）叫做预报校正公式，其中第一式叫预报公式，第二式叫校正公式。这个公式还可写为(9.9)3公式的截断误差现在来考察两个公式的截断误差：有多大？这里假定前一步得的结果yn=y(xn)是准确的。写出y(xn+1)的泰勒展开式为由欧拉法得两式相减得即欧拉法的截断误差为0(h2)，当h 0时它与h2是同阶无穷小量。对于改进的欧拉方法，我们以迭代一次的预报校正格式(9.9)
8、为例来说明。因为用它们代入（9.9）式第二式即得这里末把含有h的三次幂以上的项写出，因此有即迭代一次的预报校正格式（9.9）的截断误差为0(h3)。可见改进的欧拉方法比欧拉法的阶提高了。例9.2在区间0, 1.5上，取h = 0.1，求解。解：（1）用欧拉法计算公式如下：（2）用迭代一次的改进欧拉法计算公式如下：本题的精确解为，可用来检验近似解的精确程度。计算结果如下表：xn欧拉法yn迭代一次改进欧拉法yn准确解01110.11.11.0959091.0954450.21.1918181.1840961.1832160.31.2774381.2602011.2649110.41.3582131
9、.3433601.3416410.51.4351331.4161021.4142140.61.5089661.4829561.4832400.71.5803381.5525151.5491930.81.6497831.6164761.6124520.91.7177791.6781681.6733201.01.7847701.7378691.7320511.11.851181.7958221.7888541.21.9174641.8522421.8439091.31.9840461.9073231.8973671.42.0514041.9612531.9493591.52.1200522.014
10、2072.0000002 龙格库塔法由上节知道，截断误差的阶是衡量一个方法精度高低的主要依据。能否用提高截断误差阶来提高方法的精确度呢？回答是肯定的。本节介绍的泰勒级数法和龙格库塔法就是基于这种思想构造出来的。1泰勒级数法如果初值问题（9.10）的精确解y(x)在x0, x上存在k + 1阶导数且连续，那么由泰勒公式（9.11）其中截断误差为（9.12）略去截断误差，并用近似值代替真值则得（9.13）用公式(9.13)解初值问题的数值方法称为泰勒级数法，当 k = 3时，（9.13）变为（9.14）这时的截断误差是（9.15）从截断误差的表示式中看出，如果微分方程（9.10）的真解y = y(
11、x)为次数不超过3的多项式时，公式（9.14）精确地成立，因此（9.14）是3阶方法。例9.3 导出用三阶泰勒级数法解方程的计算公式解:因故而其中表示f(x, y)对x的k阶偏导数在x = xn点上的值。泰勒级数法只要初值问题的真解充分光滑，就可以获得精确度较高的数值解。但是须计算y(x)的各阶导数，这当f (x, y)的表达式复杂时是很繁琐的。因此泰勒级数法一般只用于求“表头”（即开头几点的数值解，如y1, y2, y3, , y4等）。另外，用上述级数法计算表头时，还可以得到选择步长h的信息。假定我们要求计算误差不超过e ，那么，当h满足条件（A）（B）时，应该认为是最好的。因为，当条件（
12、A）不满足时，达不到指定精确度，而当条件（B）不满足则表明h过小。能否构造一种格式，既保留泰勒级数法精确度较高的优点，又避免过多的计算f (x, y)的各阶偏导数呢？下面介绍的龙格库塔方法就能办到这一点。2龙格库塔法从理论上讲，只要函数y = y(x)在区间a, b上充分光滑，那么它的各阶导数值y(k)(xn)与函数y(x)在区间a, b上某些点的值就相互有联系，就是说，函数值可用各阶导数值近似地表示出来，反之，各阶导数值也可用函数在一些点上值的线性组合近似地表示出来。事实上，欧拉法和改进欧拉法也可以看成是导数值用函数值的线性组合表示的特例，例如，改进欧拉法可以写成（9.16）此公式也可称为二
13、阶龙格库塔公式。为了导出龙格库塔法的一般公式，我们取如下的线性组合形式，（9.17）其中（9.18）即而w1, w2, wv；b1 =0, b2, b3, bv；a21, a31, , avv-1除b1=0外均为待定系数。适当选取这些系数，使得局部截断误差的阶尽可能高即可。显然，当g = 1时，（9.12）式就是欧拉公式。下面我们先导出g =2时的公式。将k1, k2在同一点（xn, yn）泰勒展开，则有（9.19）将(9.19)代入(9.17)并与y(xn+h)在xn点的泰勒展开式：逐项比较，令h、h2项的系数相等，便得到把b2作为自由参数来确定w1和w2，如取b2 = 1，则w1 = w2 = ，a21 = 1，这时（9.17）正好就是改进的欧拉方法，截断误差的阶为0(h3)。对于g = 3的情形，我们也可以完全仿上述方法推导出三阶龙格库塔公式。这时参数满足下列条件（9.20）（9.20）比较简单的一组解为：b2 =，b3 = 1，a21 = ，a31 = -1，a32 = 2，w1 = ，w2 = ，w3 = 将它代入（9.17）得
