1、毕业论文文献综述 信息与计算科学 常微分方程初值问题的预估 -校正解法 一、前言部分 在生产实际和其他数学分支中，都会不断地遇到常微分方程，而在这些方程中，仅有很少的一部分能通过初等积分法给出通解或通积分，大多数积分必须数值计算。所以，一开始就使用数值方法求解通常更有效 1 。 解常微分方程初值问题的数值方法通常可以分为两类 2 ：（ 1）单步法，例如 Euler 方法和 Runge-Kutta 方法；（ 2）多步法，例如线性多步法。 我们 将同阶的显式公式与隐式公式相比，前者使用方便，计算量较小；而后者一般需用迭代法求解，计算量大，但其局部截断误差较小，稳定性较好。两种方法各有长处和不足。因
2、此，常常将它们配合起来使用，以发挥它们的优点，弥补各自的不足 3 。 这样将显式公式和隐式公式联合使用，前者提供预测值，而后者将预测值加以校正，使数值解更精确。由此形成的算法通常被称作预估 -校正算法（简称为 PC 算法） 原则上任一显式多步法和隐式多步法都可以搭配成预估校正算法及各种计算方案，但不是任一种方案都是可用的。 一个 好 的计算方案应该计算稳定，具有所需的精度，并且节约计算量 4 。几种常见的预估 -校正算法 5 ：（ 1） Adams 四阶预估 -校正算法； (2)Milne 方法(3)Hamming 算法。 本文综述 常微分 初值问题的数值解法及其误差估计（ 相容性、稳定性和收
3、敛性分析 ），重点介绍了 预估 -校正算法。 二、主题部分 2.1 常微分方程的起源和发展 6 许多有关微分方程的教材都会提到发现海王星的故事 。 海 王星的发现是人类智慧的结晶，也是常微分方程巨大作用的体现 ， 体现了数学演绎法的强大威力 。 1781年发现天王星后，人们注意到它所在的位置总是和万有引力定律计算出来的结果不符 。 于是有人怀疑万有引力定律的正确性 ； 但也有人认为 ， 这可能是受另外一颗尚未发现的行星吸引所致 。 当时虽有不少人相信后一种假设 ， 但缺乏去寻找这颗未知行星的办法和勇气 。 23岁的英国剑桥大学的学生亚当斯承担了这项任务 ， 他利用引力定律和对天王星的观测资料建
4、立起微分方程 ， 来求解和推算这颗未知行星的轨道 。 1843年 10月 21日，他把计算结果寄给格林威治天文台 台长艾利 ，但艾利不相信“小人物”的成果 ， 置之不理 。 两年后 ， 法国青年勒威耶也开始从事这项研究。1846年 9月 18日 ， 他把计算结果告诉了柏林天文台助理员卡勒 。 6日晚 ， 卡勒果然在勒威耶预言的位置上发现了海王星 。 对于数学 ， 特别是数学的应用 ， 微分方程所具有的重大意义主要在于 ： 很多物理与技术问题可以化归为微分方程的求解问题 。 常微分方程发展过程中所经历的四个重要时期： 1.常微分方程的经典阶段 以通解为主要研究内容； 2.常微分方程的适定性理论阶
5、段 以定解问题的适定性理论为研究内容；3.常微分方程的解析理论阶段 以 解析理论为研究内容； 4.常微分方程的定性理论阶段 以定性与稳定性理论为研究内容。 总之 ， 微分方程是一门十分有用又十分有魅力的学科 。 自 1693年微分方程概念的提出到动力系统的长足发展 ， 常微分方程经历漫长而又迅速的发展 ， 极大丰富了数学家园的内容 。随着社会技术的发展和需求 ， 微分方程会有更大的发展 ， 比如偏微分方程的迅速发展 。 可以预测 ： 随着依赖数学为基础的其它学科的发展 ， 微分方程还会继续扩展 。 2.2 常微分方程初值 问题的数值解法 7 $1.基本概念 讨论“初值问题”，通 常集中于讨论“
6、 1阶初值问题”： 00 )(),( yxy yxfy （ 2.2.1） 这是因为，任何高阶方程或方程组的初值问题，经过适当的变换可以化为 1阶方程组的初值问题，而 1阶方程组的初值问题写成向量形式，并把向量形式写作标量形式来叙述，这就是初值问题的基本形式（ 2.2.1）。 设（ 2.2.1）中的 ),( yxf 是在 RybxayxD ,|),( 上的连续函数；又 ),( yxf关于 y 满足 Lipschitz条件，即存在与 21,yy 无关的常数 0L ，使得 2121 ),(),( yyLyxfyxf , Ryy 21, 成立，则 Rybax 00 , ，初值问题（ 2.2.1）的解存
7、在、唯一、连续可微且连续依赖于初始条件 8 。 $2.常微分方程初值问题的数值解法 单步法 (1)Euler方法 Euler方法是求解初值问题 （ 2.2.1）的一种最简单、最基本的数值方法，它不一定作为一种独立的求解方法在实际中使用，但却提供了数值解法的本质思想。 利用 Taylor展开法，得近似计算公式 ),()()( 1 nnnn yxhfxyxy ， ,.)1,0( n （ 2.2.2） 从 0x 何初值 00)( yxy 出发进行迭代。（ 2.2.2）称为 Euler公式，它是一种迭代公式，也是一种差分公式；它利用前一步的数值解 ny 便可算出下一步的数值解 1ny ，所以 Eule
8、r公式也成为显式单步法。 利用向后差商，得 近似计算公式 ),()()( 111 nnnn yxhfxyxy （ 2.2.3） 待求解的 1ny 还隐含在计算公式等号右边的 f 中，故称它为隐式 Euler公式。隐式公式的计算，除非 f 比较特殊，一般只好采用迭代法求解，像求解一元非线性方程的迭代法一样，比较麻烦。 把显式（ 2.2.2）和隐式（ 2.2.3）作算术平均，则得 ),(),(2)()(111 nnnnnn yxfyxfhxyxy（ 2.2.4） 称为梯形公式。当然，它仍是隐式公式，但它的精度 提高了。 取一次迭代（即由显式 Euler公式算出 1ny 并记为 1ny ，代入梯形公
9、式右端的 1ny ，使梯形公式由隐式变为显式）可得如下的格式： ),(),(2),(1111nnnnnnnnnnyxfyxfhyyyxhfyy (2.2.5) 这称为改进的 Euler公式 (2)Runge-Kutta方法 设想要构造的目标公式的形式为 ). . .,(.),(),(). . .(11,1112122122111RRRRnRnRnnnnRRnnhkqhkqyhpxfkhkqyhpxfkyxfkkkkhyy （ 2.2.6） 其中， isii qp, 等均为待定参数， 1R ，式中用到 R 个斜率值 ),.,2,1( Riki 作加权平均以提高方法的阶，故称为 R 级的 Rung
10、e-Kutta方法。 对 1R ，式（ 2.2.6）实际就是 Euler公式。 对 4R ，推导可得含 13个未知数 11个方程的参数约束条件，从而得 4级 4阶 Runge-Kutta方法族。其中，一种最著名、最常用的称为经典 4阶 Runge-Kutta公式 ),()2,2()2,2(),()22(6342312143211hkyhxfkkhyhxfkkhyhxfkyxfkkkkkhyynnnnnnnnnn（ 2.2.7） $3.常微分方程初值问题的数值解法 多步法 设线性公式关于 ky 和 kf 为线性，考虑 l 步法，则一般形式为 (2.2.8) )()()( 10111 lklk k
11、nknn khxyhkhxyhxyT （ 2.2.9） 特别地，当 01 时，公式为显式公式，当 01 时，公式为隐式公式。 1 1101 lk knkknklkn fhyy (1) Adams 4步 4阶显式公式（也称 Adams-Bashforth方法） 取 4,4 lp ， 0,0 3211 ，得 )9375955(243211 nnnnnn ffffhyy(2.2.10) )()(7 2 02 5 1 6)5(51 hOxyhT nn (2.2.11) (2) Adams 3步 4阶隐式公式（也称 Adams-Moulton方法） 取 3,4 lp ， 0,0 211 ，得 )5199
12、(242111 nnnnnn ffffhyy(2.2.12) )()(72019 6)5(51 hOxyhT nn (2.2.13) (3) Simpson 2步 4阶隐式公式 取 2,4 lp ， 0,0 201 ，得 )4(31111 nnnnn fffhyy(2.2.14) )()(901 6)5(51 hOxyhT nn (2.2.15) (4) Milne 4步 4阶显式公式 )22(342131 nnnnn fffhyy(2.2.16) )()(4514 6)5(51 hOxyhT nn (2.2.17) (5) Hamming 3步 4阶隐式公式 )2(83)9(811121 n
13、nnnnn fffhyyy(2.2.18) )()(401 6)5(51 hOxyhT nn (2.2.19) $4.预估 -校正算法 由上述诸公式的局部截断误差可以粗略看到，隐式公式比显式公式精度略高。在处置问题数值解得稳定性分析中，我们也提到，隐式公式比显式公式稳定性好。因此，实际计算中，通常利用同阶显式公式先做“预估”，再用隐式公式做一次“校正”，甚至为了提高精度，还可利用事后误差估计“修正”。通过这些操作，可以提供相当理想的预估 -校正计算方案。 常用的一种预估 -校正方案是用 4阶 4步 Adams显式公式做预测（ Predictor） ,同阶 3步隐式公式做一次校正（ Correc
14、tor） ,这就是下列的 （ 1） Adams预估 -校正格式（称 PC格式） )519),(9(24:)9375955(24:21011132101nnnnnnnnnnnnnfffyxfhyyCffffhyyP （ 2.2.20） 实际计算中，在预估和校正之后分别紧跟一个为下一步做准备的函数值计算（ Evaluation） ),(),( 1110 110 1 nnnnnn yxffyxff 如果再考虑修正技巧，由（ 2.2.20） 第一式和第二式，以及注意到 Adams的局部截断误差式(2.2.11)和 (2.2.13)，于是有 )(7 2 019),(7 2 02 5 1 0 12 12
15、1100 11 1 nnnnnnnn yyyyyyyy 于是，可得含修正技巧的 Adams预估 -校正格式（ PMECME格式）： ),(:72019:)5199(24:),(:720251:)9375955(24:1110121211211121111110011132101nnnnnnnnnnnnnnnnnnnnnnnnnnyxffEyyyyMffffhyyCyxffEyyyyMffffhyyP另一个预估 -校正方案是将 4阶 4步 Milne显式公式做一次预测 ,同阶 3步 Hamming隐式公式做一次校正 ,并类似上述的设计，利用事后误差估计做修正，可以得下列的所谓 （ 2） Miln
16、e-Hamming预估 -校正格式（ PMECME） ),(:1209:)2(83)9(81:),(:121112:)22(34:1110121211111221111110011121301nnnnnnnnnnnnnnnnnnnnnnnnnyxffEyyyyMfffhyyyCyxffEyyyyMfffhyyP$5.收敛性与稳定性分析 10,9 考虑求解初值问题（ 2.2.1）的线性 k 步方法，即把式（ 2.2.8） 的线性多步法改写成形式略有差异的 k 步线性公式 （ 2.2.21） 记 ). .()( 122110 kkkkk aaaa 1221101 . .)( kkkklk abbb
17、b 分别称 )( 和 )( 为第一特征多项式和第二特征多项式。 如果线性 k 步公式 （ 2.2.21） 的第一特征多项式 )( 的零点的模均不超过 1，并且模等于 1 的零点为单零点，则称 k 步公式（ 2.2.21） 满足根条件。 可以验证 Adams 显式公式、 Adams 隐式公式、 Simpson 隐式公式和 Hamming 公式均满足根条件。 线性 k 步法 （ 2.2.21） 是 k （ 1k ）阶 相容的，则其收敛的充分必要条件是根条件满足。 线性 k 步法 （ 2.2.21） 稳定 的充分必要条件是它满足根条件 12,11 。 现在假设由初值和线性 k 步法求得的解表示如下：
18、 ),( 00 hyfhyiikj jnjkj jnj ， )1,.,1,0( li （ 2.2.22） 其中， )()( 000 xyhy 。 设 初值问题（ 2.2.1）中 ),( yxf 连续，且关于 y 满足 Lipschitz 条件。如果线性多步法 （ 2.2.22）的解 ny 对任意的 nhxx 0 ，有 0 ()lim nh y y x ， 且初始条件满足 kj jnjkj jnj fhy 00 00 ()lim ih hy ， )1,.2,1( li 则称线性 k 步法 （ 2.2.21） 是收敛的 13 。 把 线性 k 步法 （ 2.2.21） 用于求解试验方程 yy ，
19、)0)(Re( ， 这时 线性 k 步法简化为 kj jnjkj jnj yhy 00 (2.2.23) 记 hh ，相应于 (2.2.23)的特征方程为 0 hr hh （ 2.2.24） 对于给定的 h ，如果特征方程（ 2.2.24） 的根 i 都满足 1i ，则称线性多步法 (2.2.23)关于此 h 是绝对稳定的。 设 AZ 为 h 复平面内的一个集合，如果对全部 AZh ，线性多步法 (2.2.23)都是绝对稳定的，则称此线性多步法有绝对稳定域 AZ ； AZ 与实轴之交称为此线性多步法的绝对稳定区间。 下面引述若干在应用中可供参考的数据 14 （ 1） 步数 41k 的 Adam
20、s 显式与隐式方法的绝对稳定区间如下表： 步数 k Adams 显式方法 Adams 隐式方法 阶 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 方法组成预
21、估 -校正公式，则有较好的稳定性。 （ 3） Simpson 方法的绝对稳定域为空集。所以方法不绝对稳定。实数 对应yf，不管yf取什么符号，误差将连续增长，所以 在实际计算中不推荐 Simpson 方法。 三、总结部分 常微分方程是建立物理、生物科学以及工程学的模型时所使用的最重要的数学工具之一。本文考虑的是求解常微分方程的数值方法。首先讨论了常微分方程初值问题的理论，并到引入微分方程稳定性的概念。给出求解的最简单数值方法 Euler 方法 ，它并不是有效地数值方法，但研究它有利于引入和理解许多重要思想。随后谈论在实际计算中具有较好数值稳定性的方法 向后 Euler 方法、梯形方法。最后着重
22、于更成熟、收敛更快速的方法 Runge-Kutta 方法、 Adams-Bashforth 方法、 Adams-Moulton 方法、 Simpson 方法、 Milne方法和 Hamming 方法。由于 显式公式与隐式公式 和 各有长处和不足，故 实际计算中，通常将两者结合起来，利用同阶显式公式先做“预估”，再用隐式公式做一次“校正”，甚至为了提高精度，还可利用事后误差估计“修正”。通过这些操作，可以提供相当理想的预估 -校正计算方案，如 Adams 预估 -校正算法 和 Milne-Hamming 预估 -校正算法。 并 给出相应的误差分析，包括相容性、稳定性和收敛性分析。 上述各种预估
23、-校正算法通常需要通过编程上机计算。 MATLAB提供了多个常微分方程求解的函数 ,ode45,ode113,ode15s,ode23s,ode23t,ode23tb,其中两个为 ode23 () 和 ode45 () , 分别采用了二阶 P三阶的 RTF 方法和四 P五阶层 RKF方法 , 并采用自适应变步长的求解方法 , 即当解的变化较慢时采用较大的计算步长 , 从而使得计算速度很快 , 当方程的解变化得较快时 , 积分步长会自动地变小 , 从而使得计算的精度很高 15 。 四、参考文献 1Kendall Atkinson，韩渭敏 著，王国荣，徐兆亮，孙劼译数值分析导论 (第三版 )M北京
24、 ： 人民邮电出版社， 2009 2余德浩，汤华中微分方程数值解法 M北京 ： 科学出版社， 2003 3孙志忠，袁慰平，闻震初数值分析 (第二版 )M南京 ： 东南大学出版社， 2002 4李荣华，冯果忱微分方程数值解法 M北京 ： 高等教育出版社， 1995 5汤怀民，胡建伟微分方程数值方法 (第三版 )M南开大学出版社， 1990 6张良勇，董晓芳常微分方程的起源与发展 J高等函数学报 (自然科学版 )， 2006， 20(3)： 1-2 7郑咸义，姚迎新，雷秀仁，陆子强，黄凤辉应用数值分析 M广州 ： 华南理工大学出版社， 2008 8Arieh Lserles 著，刘晓艳，刘学深译微
25、分方程数值分析基础教程 M北京 ： 清华大学出版社， 2005 9Richard L.Burden and J.Douglas NUMERICAL ANALYSISM Wadsworth Group Brooks， 2001 10Ascher U.M. and Petzold L.R Computer methods for ordinary differential equations and differential-algebraic equationsM Society for Industrial Mathematics， 1998 11林成森数值计算方法 M北京 ： 科学出版社， 1998 12李立康，於崇华，朱政华微分方程数值解法 M上海 ： 复旦大学出版社， 1999 13关治，陆金甫数值分析基础 M北京 ： 高等教育出版社， 1998 14清华大学应用数学系现代应用数学手册 编写组现代应用数学书册（计算方法分册） M.北京 ： 北京出版社， 1990 15王沫然 ， 肖劲松 MATLAB5.X 与科学计算 M 北京 ： 清华大学出版社 ， 2000
