1、 本科 毕业论文 （ 设计 ） （ 20 届） 常微分方程初值问题的预估 -校正解法 所在学院 专业班级 信息与计算科学 学生姓名 学号 指导教师 职称 完成日期 年 月 I 摘要： 求解常微分方程初值问题的有限差分法一般分为显式方法和隐式方法。显式方法计算简单 ，而隐式方法稳定性好、精度高，故 实际计算中，通常将两者结合起来，利用同阶显式公式先做“预估”，再用隐式公式做一次“校正”，即所谓的预估 -校正算法。本文综述了预估 -校正算法的理论基础、算法构造以及算法实现，并通过数值例子验证理论分析。 关键词： 有限差分法； 预估 -校正算法 ； 算例 II Predictor - Correct
2、or Algorithms for Initial-Value Problems for Ordinary Differential Equations Abstract: Finite difference methods for solving Initial-value problems for ordinary differential equations can be divided into explicit methods and implicit methods generally. For the explicit methods, it can be calculated
3、easily than the corresponding implicit methods, but the implicit methods has good stability and higher degree of accuracy, so in practical calculation, we usually combine them together, using explicit method as a “predictor“, then implicit method as a “corrector“, it is called Predictor - Corrector
4、Algorithms. In this paper, we summarize the theory basis, algorithms construction and implementation of Predictor - Corrector Algorithms, then several numerical examples are provided to confirm the theoretical analysis. Key words: finite difference methods ； predictor - corrector algorithms；numerica
5、l examples III 目 录 1 绪 论 . 错误 !未定义书签。 1.1 问题的背景 . 错误 !未定义书签。 1.2 有限差分法的简介 . 错误 !未定义书签。 2 常微分方程初值问题的预估 -校正算法 . 错误 !未定义书签。 2.1 预估 -校正算法的若干准备 . 错误 !未定义书签。 2.1.1 常微分方程初值问题的基本概念 . 错误 !未定义书签。 2.1.2 常微分方程初值问题的数值解法 单步法 . 3 2.1.3 常微分方程初值问题的数值解法 多步法 . 8 2.2 若干个预估 -校正算法 . 11 2.2.1 预估 -校正算法的构造 . 11 2.2.2 预估 -校算
6、法的实现 . 错误 !未定义书签。 3 数值算例 . 21 3.3.1 算例 1. 21 3.3.2 算例 2. 22 4 结 论 . 错误 !未定义书签。 致 谢 . 错误 !未定义书签。 参考文献 . 错误 !未定义书签。 附 录 . 28 1 1 绪论 1.1 问题的背景 许多有关微分方程的教材都会提到发现海王星的故事 。 海王星的发现是人类智慧的结晶，也是常微分方程巨大作用的体现 ， 体现了数学演绎法的强大威力 。 1781年发现天王星后，人们注意到它所在的位置总是和万有引力定律计算出来的结果不符 。 于是有人怀疑万有引力定律的正确性 ； 但也有人认为 ， 这可能是受另外一颗尚未发现的
7、行星 吸引所致 。 当时虽有不少人相信后一种假设 ， 但缺乏去寻找这颗未知行星的办法和勇气 。 23岁的英国剑桥大学的学生亚当斯挑下了这项任务 ， 他利用引力定律和对天王星的观测资料建立起微分方程 ， 来求解和推算这颗未知行星的轨道 。 1843年 10月 21日，他把计算结果寄给格林威治天文台台长艾利 ， 但艾利不相信“小人物”的成果 ，对其 置之不理 。 两年后 ， 法国青年勒威耶也开始从事这项研究。 1846年 9月 18日 ， 他把计算结果告诉了柏林天文台助理员卡勒 。 6日晚 ， 卡勒果然在勒威耶预言的位置上发现了海王星 1。 对于数学 ， 特别是数学的应用 ， 微分方 程所具有的重
8、大意义主要在于 ： 很多物理与技术问题可以转化为微分方程的求解问题 ，而这些问题 仅有很少的一部分能通过初等积分法给出通解或通积分，大多数积分必须进行数值计算。所以，一开始就使用数值方法求解通常更加有效。常用的解常微分方程初值问题的数值方法有有限差分法和有限元法，其中有限差分法可以分为两类：（ 1）单步法，例如 Euler方法和 Runge-Kutta方法；（ 2）多步法，例如线性多步法 2。我们将同阶的显式公式与隐式公式相比较，发现前者使用方便，计算量较小；而后者一般需用迭代法求解，计算量大，但其局部截断 误差较小，稳定性较好。两种方法各有其长处和不足。因此，常常将它们配合起来使用，以发挥它
9、们的优点，弥补各自的不足 3。而将显式公式和隐式公式联合起来使用，前者提供预测值，后者将预测值加以校正，可以使数值解更精确。这种算法通常被称作预估 -校正算法（简称为 PC算法）。 1.2 有限差分法的简介 45 有限差分法是 微分方程和积分方程数值解的方法 。 因其比较简单，所以在实际应用的数值方法中占很大比例。由于数字电子计算机只能存储有限数据和作有限次运算，所以任何一种适用于计算机的解题方法，都必须把连续问题离散化， 最终化成有限形式的线性代数方程组。 有限差分 法将求解域划分为差分网格，用有限网格节点代替连续的求解域。 对其用 Taylor级数展开，把控制方程中的导数用网格节点上的函数
10、值的差商 来 代替，从而建立以网格节点上的值为未知数的 线性 代数方程组。该方法是一种直接将微分问题变为代数问题的近似数值解法， 对于矩形区域，有限差分数值解是一种行之有效的方法，由于差分法 数学概念直观，表达简单， 因而成为 发展较早且较成熟的数值方法 。对于微分方程定解问题中的每一个差分方法，基本上按照（ 1）差分2 格式的建立；（ 2）差分格式的求解；（ 3）算例； （ 4）差分格式的可解性、收敛性和稳定性等 4个方面展开。 本文将从以下两个方面进行展开：首先，在第二部分， 讨论了常微分方程初值问题的理论，并引入常微分方程稳定性的概念。然后给出求解的最简单数值方法 Euler方法 ，实际
11、计算中具有较好数值稳定性的方法 向后 Euler 方法、梯形方法，以及更成熟、收敛更快速的方法 Runge-Kutta 方法、 Adams-Bashforth 方法、 Adams-Moulton 方法、 Simpson 方法、 Milne 方法和Hamming方法。实际计算中，通常将两者结合起来，利用同阶显式公 式先做“预估”，再用隐式公式做一次“校正”，甚至为了提高精度，还可利用事后误差估计“修正”。通过这些操作，可以得到相当理想的预估 -校正计算方案，如 Adams预估 -校正算法 和 Milne-Hamming预估 -校正算法 ；最后，在第三部分，数值算例中，我们通过几个算例验证了我们的
12、理论分析。 3 2 常微分方程初值问题的预估 -校正算法 2.1 预估 -校正算法的若干准备 2.1.1 常微分方程初值问题的基本概念 讨论“初值问题”，通常集中于讨论“ 1阶初值问题”： 00 )(),( yxy yxfy （ 2.1.1-1） 这是因为，任何高阶方程或方程组的初值问题，经过适当的变换，都可以化为 1阶方程组的初值问题，而 1阶方程组的初值问题写成向量形式，并把向量形式用标量形式来叙述，就是初值问题的基本形式（ 2.1.1-1）。 设（ 2.1.1-1）中的 ),( yxf 是在 RybxayxD ,|),( 上的连续函数；且),( yxf 关于 y 满足 Lipschitz
13、条件，即存在与 21,yy 无关的常数 0L ，使得 2121 ),(),( yyLyxfyxf , Ryy 21, 成立，则 Rybax 00 , ，初值问题（ 2.1.1-1）的解存在、唯一、连续可微且连续依赖于初始条件 6。 2.1.2 常微分方程初值问题的数值解法 单步法 7 (1)Euler方法 Euler方法是求解初值问题 （ 2.1.1-1）的一种最简单、最基本的数值方法，它不一定作为一种独立的求解方法在实际中使用，但却提供了数值解法的本质思想。 利用 Taylor展开法，设（ 2.1.1-1）的解 )(xy 充分光滑，在已取定的 ,.)1,0( bxn 处Taylor展开有 )
14、(!2 )()()()( 321 hOhxyhxyxyxy nnnn 略去 2h 以上的项，并分别用 ny ， 1ny 近似表示 )(nxy ， )( 1nxy ，且根据方程可知),()(,()( nnnnn yxfxyxfxy ，于是 得近似计算公式 ),()()( 1 nnnn yxhfxyxy ， ,.)1,0( n （ 2.1.2-1） 这就使得从初值 0x 和 )( 00 xyy 出发，用式（ 2.1.2-1）可得 )(1xy 的近似值 1y ，再由 1x ，1y ，用式（ 2.1.2-1）又可得 )(2xy 的近似值 2y ，如此继续，即可得初值问题 （ 2.1.1-1）4 的 E
15、uler方法，式 （ 2.1.2-1）称为 Euler公式，它是一种迭代公式，也是一 种差分公式；它利用前一步的数值解 ny 算出下一步的数值解 1ny ，所以 Euler公式也称为显式单步法。 下面继续推导求解初值问题 （ 2.1.1-1）的其他方法。现在考虑 利用向前差商 )(,()()()( 1 nnnnn xyxfxyh xyxy 和向后差商 )(,()()()( 1111 nnnnn xyxfxyh xyxy ， 整理并分别以 1ny ， ny 近似表示 )( 1nxy ， )(nxy ，且令近似式为等式，即可得 近似计算公式 ),(1 nnnn yxhfyy 和 ),( 111 n
16、nnn yxhfyy （ 2.1.2-2） 显然，前者仍是 Euler公式（ 2.1.2-1），而后者已有了本质的不同，待求解的 1ny 还隐含在计算公式等号右边的 f 中，故称它为隐式 Euler公式。对应地，称公式（ 2.1.2-1）为显式 Euler公式。隐式公式的计算，除非 f 比较特殊，一般只好采用迭代法求解，像求解一元非线性方程的迭代法一样，比较麻烦。所以我们换为另一种思路： 把显式（ 2.1.2-1）和隐式（ 2.1.2-2）作算术平均，则得 ),(),(2)()(111 nnnnnn yxfyxfhxyxy（ 2.1.2-3） 这称为梯形公式。当然，它仍是隐式公式，但它的精度提
17、高了。下面我们以此 为例，来看看对隐式公式如何做迭代求解。对每一点 nx 以显式公式（ 2.1.2-1）提供迭代初值，对梯形公式构造迭代格式 , . . .1,0, . . . ;1,0),(),(2),()(11)1(1)0(1knyxfyxfhyyyxhfyyknnnnnknnnnn5 这时，如果迭代序列 )1(1kny 当 k 时收敛于 1ny ，便可取适当的 )1(1kny 作为 1ny 的近似值。这里容易证明，迭代 对每一个 n ，当 k 时是收敛的。事实上，因为 ),( yxf 关于 y 满足Lipschitz条件，因此 |2|),(),(|2| )( 11)( 1111)1( 1
18、1 knnknnnnknn yyLhyxfyxfhyy 可见，只要选取 h ，使得 12 Lh ，则有 | )( 11)1( 11 knnknn yyyy 它表示每多迭代一次，误差总是小一点，故 当 k 时，便有 1)( 1 nkn yy ，即迭代收敛。 既然迭代格式是收敛的，则取一次迭代（即由显式 Euler公式算出 1ny 并记为 1ny ，代入梯形公式右端的 1ny ，使梯形公式由隐式变为显式）可得如下的格式： ),(),(2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy (2.1.2-4) 这称为改进的 Euler公式（其实质是为了使用较高精度的梯形公式，借助显式 E
19、uler公式使梯形公式显式化）。 (2)Runge-Kutta方法 设想要构造的目标公式的形式为 ). . .,(.),(),(). . .(11,1112122122111RRRRnRnRnnnnRRnnhkqhkqyhpxfkhkqyhpxfkyxfkkkkhyy （ 2.1.2-5） 其中， isii qp, 等均为待定参数， 1R ，式中用到 R 个斜率值 ),.,2,1( Riki 作加权平均以提高方法的阶，故称为 R 级的 Runge-Kutta方法。 对 1R ，式（ 2.1.2-5）实际就是 Euler公式。 6 对 2R ，有 ),(),()(12122111phkyphxf
20、kyxfkkkhyynnnnnn ， （ 2.1.2-6） 其中， )10(, 21 pp 为 3个待定系数。通过确定这 3个参数 ，使得式（ 2.1.2-6）具有 2阶精度。 为此，按局部截断误差的要求，分别计算出 )( 1nxy 和 1ny ，并考虑 )( nn xyy 。现设初值问题（ 2.1.1-1）的 ),( yxf 和解 )(xy 充分光滑，于是，一方面，由 Taylor展开有 )()(2)()()( 321 hOxyhxyhxyxy nnnn )()(,(2)(,()( 32 hOxyxfhxyxhfxy xnnnnn )()(2)()( 3)(2)( hOfffhfhxy ny
21、xnn （ 2.1.2-7） 其中以下标“ )(n ”表示 f 及其导数在 )(,( nn xyx 处取值。另一方面，当解为式（ 2.1.2-6）时，由式（ 2.1.2-6）中的 1k 和 2k 作二元 Taylor展开，并以下标“ n ”表示 f 及其导数在),( nn yx 处取值，有 ),(,(),( 211 nnnnnnnn yxphfyphxfyxfhyy )()()()()()( 221 hOfphfphffhfhy nnynxnnn )()()()( 32221 hOfffphfhy nyxnn （ 2.1.2-8） 由 于 考 虑 的 是 局 部 截 断 误 差 ， 应 有 nn yxy )( ， 故 有 nn ff )()( )( ，)()( nyx fff = nyx fff )( 。于是，由 式 （ 2.1.2-7）与式（ 2.1.2-8）两边相减可知，要使得式（ 2.1.2-6）具有二阶精度，需且只需 2112)(222122221phphhhh （ 2.1.2-9）
