1、毕业论文 开题报告 信息与计算科学 常微分方程初值问题的预估 -校正解法 一、选题的背景、意义 1.1 历史背景 许多有关微分方程的教材都会提到发现海王星的故事 。 海王星的发现是人类智慧的结晶，也是常微分方程巨大作用的体现 ， 体现了数学演绎法的强大威力 。 1781年发现天王星后，人们注意到它所在的位置总是和万有引力定律计算出来的结果不符 。 于是有人怀疑万有引力定律的正确性 ； 但也有人认为 ， 这可能是受另外一颗尚未发现的行星吸引所致 。 当时虽有不少人相信后一种假设 ， 但缺乏去寻找这颗未知行星的办法和勇气 。 23岁的英国剑桥大学的学生亚当斯承担了这项任务 ， 他利用引力定律和对天
2、王星的观测资料建立起微分方程 ， 来求解和推算这颗未知行星的轨道 。 1843年 10月 21日，他把计算结果寄给格林威治天文台台长艾利 ，但艾利不相信“小人物”的成果 ， 置之不理 。 两年后 ， 法国青年勒威耶也开始从事这项研究。1846年 9月 18日 ， 他把计算结果告诉了柏林天文台助理员卡勒 。 6日晚 ， 卡勒果然在勒威耶预言的位置上发现了海王星 1 。 对于数学 ， 特别是数学的应用 ， 微分方程所具有的重大意义主要在于 ： 很多物理与技术问题可以化归为微分方程的 求解问题 ，而这些问题 仅有很少的一部分能通过初等积分法给出通解或通积分，大多数积分必须数值计算。所以，一开始就使用
3、数值方法求解通常更有效 2 。解常微分方程初值问题的数值方法通常可以分为两类：（ 1）单步法，例如 Euler方法和 Runge-Kutta方法；（ 2）多步法，例如线性多步法 3 。我们将同阶的显式公式与隐式公式相比，前者使用方便，计算量较小；而后者一般需用迭代法求解，计算量大，但其局部截断误差较小，稳定性较好。两种方法各 有长处和不足。因此，常常将它们配合起来使用，以发挥它们的优点，弥补各自的不足 4 。这样将显式公式和隐式公式联合使用，前者提供预测值，而后者将预测值加以校正，使数值解更精确。由此形成的算法通常被称作预估 -校正算法（简称为 PC算法） 1.2 国内外研究现状 在生产实际和
4、其他数学分支中，都会不断地遇到常微分方程。 常微分方程在发展过程中经历了四个重要时期： 1.常微分方程的经典阶段 以通解为主要研究内容； 2.常微分方程的适定性理论阶段 以定解问题的适定性理论为研究内容； 3.常微分方程的解析理论阶段 以解析理论为研究内容； 4.常微分方程的定性理论阶段 以定性与稳定性理论为研究内容。 解常微分方程初值问题的数值方法也经历了一些发展：从最简单的单步法，例如 Euler方法和 Runge-Kutta方法，到相对复杂的多步法，例如线性多步法，最后我们将显式公式和隐式公式联合使用，前者提供预测值，而后者将预测值加以校正，形成了预估 -校正算法。 现在 预估校正算法
5、的运用非常广泛，国内外对其的研究 表明： 原则上任一显式多步法和隐式多步法都可以搭配成预估校正算法及各种计算方案，但不是任一种方 案都是可用的。一个 好 的计算方案应该计算稳定，具有所需的精度，并且节约计算量 5 。几种常见的预估 -校正算法 6 ：（ 1） Adams 四阶预估 -校正算法； (2)Milne 方法 (3)Hamming 算法。 1. 3 发展趋势 常微分方程是一门十分有用又十分有魅力的学科 。 自 1693年微分方程概念的提出到动力系统的长足发展 ， 常微分方程经历漫长而又迅速的发展 ， 极大丰富了数学家园的内容 。 随着社会技术的发展和需求 ， 微分方程会有更大的发展 。
6、 由于在实践的 运用中，常微分方程往往结构复杂，要给出一般方程的表达式是十分困难的，甚至还存在大量的不能用初等函数来表达的方程，此时，就只能用求其近似的数值解，若采用常规的人工推导、求解无意是效率非常低下的，但如果借助计算机来进行辅助，则将大大提高求解的速度和精度 7 。 所以，利用计算机运算速度快的优点，将欧拉公式、经典龙格 -库塔方法等进行计算机程序算法分析与实现，将成为研究求解常微分方程初值问题的趋势。 二、研究的基本内容与拟解决的主要问题 2.1 基本内容 8 $1.基本概念 讨论“初值问题”，通常集中于讨论“ 1阶初值问题”： 00 )(),( yxy yxfy （ 2.1.1） 这
7、是因为，任何高阶方程或方程组的初值问题，经过适当的变换可以化为 1阶方程组的初值问题，而 1阶方程组的初值问题写成向量形式，并把向量形式写作标量形式来叙述，这就是初值问题的基本形式（ 2.1.1）。 设（ 2.2.1）中的 ),( yxf 是在 RybxayxD ,|),( 上的连续函数；又 ),( yxf关于 y 满足 Lipschitz条件，即存在与 21,yy 无关的常数 0L ，使得 2121 ),(),( yyLyxfyxf , Ryy 21, 成立，则 Rybax 00 , ，初值问题（ 2.2.1） 的解存在、唯一、连续可微且连续依赖于初始条件 9 。 $2.常微分方程初值问题的
8、数值解法 单步法 (1)Euler方法 Euler方法是求解初值问题 （ 2.1.1）的一种最简单、最基本的数值方法，它不一定作为一种独立的求解方法在实际中使用，但却提供了数值解法的本质思想。 利用 Taylor展开法，得近似计算公式 ),()()( 1 nnnn yxhfxyxy ， ,.)1,0( n （ 2.1.2） 从 0x 何初值 00)( yxy 出发进行迭代。（ 2.1.2）称为 Euler公式，它是一种迭代公式，也是一种差分公式；它利用前一步的数值解 ny 便可算出下一步的数值解 1ny ，所以 Euler公式也成为显式单步法。 利用向后差商，得 近似计算公式 ),()()(
9、111 nnnn yxhfxyxy （ 2.1.3） 待求解的 1ny 还隐含在计算公式等号右边的 f 中，故称它为隐式 Euler公式。隐式公式的计算，除非 f 比较特殊，一般只好采用迭代法求解，像求解一元非线性方程的迭代法一样，比较麻烦。 把显式（ 2.1.2）和隐式（ 2.1.3）作算术平均，则得 ),(),(2)()(111 nnnnnn yxfyxfhxyxy（ 2.1.4） 称为梯形公式。当然，它仍是隐式公式，但它的精度提高了。 取一次迭代（即由显式 Euler公式算出 1ny 并记为 1ny ，代入梯形公式右端的 1ny ，使梯形公式由隐式变为显式）可得如下的格式： ),(),(
10、2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy (2.1.5) 这称为改进的 Euler公式 (2)Runge-Kutta方法 设想要构造的目标公式的形式为 ). . .,(.),(),(). . .(11,1112122122111RRRRnRnRnnnnRRnnhkqhkqyhpxfkhkqyhpxfkyxfkkkkhyy （ 2.1.6） 其中， isii qp, 等均为待定参数， 1R ，式中用到 R 个斜率值 ),.,2,1( Riki 作加权平均以提高方法的阶，故称为 R 级的 Runge-Kutta方法。 对 1R ，式（ 2.1.6）实际就是 Euler公式
11、。 对 4R ，推导可得含 13个未知数 11个方程的参数约束条件，从而得 4级 4阶 Runge-Kutta方法族。其中，一种最著名、最常用的称为经典 4阶 Runge-Kutta公式 ),()2,2()2,2(),()22(6342312143211hkyhxfkkhyhxfkkhyhxfkyxfkkkkkhyynnnnnnnnnn（ 2.1.7） $3.常微分方程初值问题的数值解法 多步法 设线性公式关于 ky 和 kf 为线性，考虑 l 步法，则一般形式为 (2.1.8) )()()( 10111 lklk knknn khxyhkhxyhxyT （ 2.1.9) 特别地，当 01 时
12、，公式为显式公式，当 01 时，公式为隐式公式。 (1) Adams 4步 4阶显式公式（也称 Adams-Bashforth方法） 取 4,4 lp ， 0,0 3211 ，得 )9375955(243211 nnnnnn ffffhyy(2.1.10) )()(7 2 02 5 1 6)5(51 hOxyhT nn (2.1.11) (2) Adams 3步 4阶隐式公式（也称 Adams-Moulton方法） 取 3,4 lp ， 0,0 211 ，得 )5199(242111 nnnnnn ffffhyy(2.1.12) )()(72019 6)5(51 hOxyhT nn (2.1.
13、13) (3) Simpson 2步 4阶隐式公式 取 2,4 lp ， 0,0 201 ，得 )4(31111 nnnnn fffhyy(2.1.14) 1 1101 lk knkknklkn fhyy )()(901 6)5(51 hOxyhT nn (2.1.15) (4) Milne 4步 4阶显式公式 )22(342131 nnnnn fffhyy(2.1.16) )()(4514 6)5(51 hOxyhT nn (2.1.17) (5) Hamming 3步 4阶隐式公式 )2(83)9(811121 nnnnnn fffhyyy(2.1.18) )()(401 6)5(51 h
14、OxyhT nn (2.1.19) 2.2 主要内容 8 $1.预估 -校正算法 由上述诸公式的局部截断误差可以粗略看到，隐式公式比显式公式精度略高。在处置问题数值解得稳定性分析中，我们也提到，隐式公式比显式公式稳定性好。因此，实际计算中，通常利用同阶显式公式先做“预估”，再用隐式公式做一次“校正”，甚至为了 提高精度，还可利用事后误差估计“修正”。通过这些操作，可以提供相当理想的预估 -校正计算方案。 常用的一种预估 -校正方案是用 4阶 4步 Adams显式公式做预测（ Predictor） ,同阶 3步隐式公式做一次校正（ Corrector） ,这就是下列的 （ 1） Adams预估
15、-校正格式（称 PC格式） )519),(9(24:)9375955(24:21011132101nnnnnnnnnnnnnfffyxfhyyCffffhyyP （ 2.2.1） 实际计算中，在预估和校正之后分别紧跟一个为下一步做准备的函数值计算（ Evaluation） ),(),( 1110 110 1 nnnnnn yxffyxff 如果再考虑修正技巧，由（ 2.2.1） 第一式和第二式，以及注意到 Adams的局部截断误差式(2.1.11)和 (2.1.13)，于是有 )(7 2 019),(7 2 02 5 1 0 12 12 1100 11 1 nnnnnnnn yyyyyyyy
16、于是，可得含修正技巧的 Adams预估 -校正格式（ PMECME格式）： ),(:72019:)5199(24:),(:720251:)9375955(24:1110121211211121111110011132101nnnnnnnnnnnnnnnnnnnnnnnnnnyxffEyyyyMffffhyyCyxffEyyyyMffffhyyP另一个预估 -校正方案是将 4阶 4步 Milne显式公式做一次预测 ,同阶 3步 Hamming隐式公式做一次校正 ,并类似上述的设计，利用事后误差估计做修正，可以得下列的所谓 （ 2） Milne-Hamming预估 -校正格式（ PMECME） )
17、,(:1209:)2(83)9(81:),(:121112:)22(34:1110121211111221111110011121301nnnnnnnnnnnnnnnnnnnnnnnnnyxffEyyyyMfffhyyyCyxffEyyyyMfffhyyP$2.收敛性与稳定性分析 11,10 考虑求解初值问题（ 2.1.1）的线性 k 步方法，即把式（ 2.1.8） 的线性多步法改写成形式略有差异的 k 步线性公式 （ 2.2.2） 记 ). .()( 122110 kkkkk aaaa 1221101 . .)( kkkklk abbbb 分别称 )( 和 )( 为第一特征多项式和第二特征多
18、项式。 如果线性 k 步公式 （ 2.2.2） 的第一特征多项式 )( 的零点的模均不超过 1，并且模等 kj jnjkj jnj fhy 00 于 1 的零点为单零点，则称 k 步公式（ 2.2.2） 满足根条件。 可以验证 Adams 显式公式、 Adams 隐式公式、 Simpson 隐式公式和 Hamming 公式均满足根条件。 线性 k 步法 （ 2.2.2） 是 k （ 1k ）阶 相容的，则其收敛的充分必要条件是根条件满足。 线性 k 步法 （ 2.2.2） 稳定 的充分必要条件是它满足根条件 13,12 。 现在假设由初值和线性 k 步法求得的解表示如下： ),( 00 hyf
19、hyiikj jnjkj jnj ， )1,.,1,0( li （ 2.2.3） 其中， )()( 000 xyhy 。 设 初值问题（ 2.1.1）中 ),( yxf 连续，且关于 y 满足 Lipschitz 条件。如果线性多步法 （ 2.2.3）的解 ny 对任意的 nhxx 0 ，有 0 ()lim nh y y x ， 且初始条件满足 00 ()lim ih hy ， )1,.2,1( li 则称线性 k 步法 （ 2.2.2） 是收敛的 14 。 把 线性 k 步法 （ 2.2.2） 用于求解试验方程 yy ， )0)(Re( ， 这时 线性 k 步法简化为 kj jnjkj jn
20、j yhy 00 (2.2.4) 记 hh ，相应于 (2.2.4)的特征方程为 0 hr hh （ 2.2.5） 对于给定的 h ，如果特征方程（ 2.2.5） 的根 i 都满足 1i ，则称线性多步法 (2.2.4)关于此 h 是绝对稳定的。 设 AZ 为 h 复平面内的一个集合，如果对全部 AZh ，线性多步法 (2.2.4)都是绝对稳定的，则称此线性多步法有绝对稳定域 AZ ； AZ 与实轴之交称为此线性多步法的绝对稳定区间。 下面引述若干在应用中可供参考的数据 15 （ 1） 步数 41k 的 Adams 显式与隐式方法的绝对稳定区间如下表： 步数 k Adams 显式方法 Adam
21、s 隐式方法 阶 p 绝对稳定区间 阶 p 绝对稳定区间 1 1 0,2 2 0, 2 2 0,1 3 0,6 3 3 0,116 4 0,3 4 4 0,103 5 0,4990 其中， Adams 显式方法在 1k 时就是单步 Euler 方法，所得绝对 稳定区间与单步法的结论一致； Adams 隐式方法在 1k 时就是梯形公式，所得结论也一致；此外，同阶的 Adams 方法中，隐式公式的绝对稳定性比显示公式的绝对稳定性好。 （ 2） Milne 方法对任何 它都不是决对稳定的。但与 Hamming 方法组成预估 -校正公式，则有较好的稳定性。 （ 3） Simpson 方法的绝对稳定域为
22、空集。所以方法不绝对稳定。实数 对应yf，不管yf取什么符号，误差将连续增长，所以在实际计算中不推荐 Simpson 方法。 三、研究的方法与技术路线、研究难点，预期达到的目标 3.1 研究的方法和路线 在自然科学和经济的许多领域，常常会遇到常微分方程的初值问题。本文综述 常微分初值问题 的数值解法及其误差估计，重点介绍了 预估 -校正算法 （该算法是结合显示公式和隐式公式两者优点发展起来的一种数值算法） 。 具体描述构造 给定问题的预估 校正 格式，给出这 个 格式的 相容性、稳定性和收敛性分析 。 最后通过计算机编程， 给出数值例子，验证理论分析。 3.2 研究的难点 原则上任一显式多步法
23、和隐式多步法都可以搭配成预估校正算法及各种计算方案，但不是任一种方案都是可用的。一个 好 的计算方案应该计算稳定，具有所需的精度，并且节约计算量。 论文的难点在于 怎样更好地建立预估 校正格式，使其在不过分增加计算工作量的前提下精度得到提高。 3.3 预期目标 通过对常微分方程初值问题及其预估 校正算法的研究，能够 构造 给定问题的预估 校正算法 ， 给出相应的误差分析，包括相容性、稳定性和 收敛性 分析，并 通过程序（ C 语言、FORTRAN 语言或 MATLAB 语言）编程，给出数值例子，验证理论分析。 四、论文详细工作进度和安排 第 7 学期第 9 周（ 2010 年 11 月 5 号）至第 7 学期第 19 周（ 2011 年 1 月 10 号） 完成毕业论文文献检索、文献综述、外文翻译及开题报告。 第 7 学期第 19 周（ 2011 年 1 月 10 号）至第 8 学期第 3 周（ 2011 年 3 月 11 号） 完成毕业论文的数据收集、论文初稿。 第 8 学期第 3 周（ 2011 年 3 月 11 号）至第 8 学期第 11 周（ 2011 年 5 月 3 号） 1、 进入实习单位进行毕业实习，对论文进行修改； 2、 第 11 周（ 2011 年 5 月 3 号）前必须返校，完成毕业实习 返校，并递交毕业实习报告，进一步完善毕业论文；
