（凝聚态物理专业论文）求解一维schrodinger方程的两步obrechkoff方法.pdf

原创性声明 本人声明：所呈交的论文是本人在导师指导下进行的研究工作。 除了文中特别加以标注和致谢的地方外，论文中不包含其他人已发表 或撰写过的研究成果。参与同一工作的其他同志对本研究所做的任何 贡献均已在论文中作了明确的说明并表示了谢意。 签名：日期： 本论文使用授权说明 本人完全了解上海大学有关保留、使用学位论文的规定，即：学 校有权保留论文及送交论文复印件，允许论文被查阅和借阅；学校可 以公布论文的全部或部分内容。 ( 保密的论文在解密后应遵守此规定) 签名：鱼塑盗导师签名：1 喇日期： 上海大学硕士学位论文6 摘要 各种科学领域问题都可用微分方程描述，例如在理论物理、粒子物理、物理 化学、量子物理和分子物理，这些问题都可用s c h r & t i n g e r 方程来研究，然而只 有少数的微分方程存在解析解。因此近似方法才是解这些微分方程的主要手段， 才能分析和研究这些问题的本质。随着上世纪4 0 年代电子计算机的发明和发展， 一门新的学科数值计算诞生并迅速发展很多寻找数值解的实用思想与方法得到 实现。特别是近几十年，一些数学软件，m a t h e m a t i c a ，m a t l a b 的出现可使研究者 用计算机更有效求解微分方程得到更精确的解。本硕士论文主要研究一种新的求 解一维s c h r 6 d i n g e r 方程的高阶微商方法，或称o b r e c h k o f f 方法。由于篇幅的限 制，下面只讨论o b r e c h k o f f 两步方法。研究表明加入奇数阶微商，可以使只含有 偶数阶微商的差分方程的截断误差大大减少。用这样的方法解一维s c h r 6 d i n g e r 方程可使精度和效率大大提高。 论文有四章，第一章介绍解s c h r o d i n g e r 方程的背景材料。 第二章，介绍差分方法的基本知识和几种解常微分方程的数值解法，包括 e u l e r 方法、r u n g e k u t t a 方法、n u m e r o v 方法、s t 6 r m e r - c o w e l l 方法和o b r e c h k o f f 方法。 第三章，介绍o b r e c h k o f f 新两步方法，这是本论文的主要内容。 第四章，给出新方法求解一维s c h r 6 d i n g e r 方程数值结果，包括熟知的 w o o d s s a x o n 势、m o r s e 势和p 6 s c h l t e l l e r 势。最后给出结论。 本论文是作者在硕士期间所做的工作和所发表的论文基础上做的深入和全 面的总结。 关键词：一维s c h r f d i n g e r 方程，数值差分方法，p 稳定，o b r e c h k o f f 方法 上海大学硕士学位论文 7 a b s t r a c t m a n yp h y s i c a lp r o b l e m so c c u r r e di nv a r i o u ss c i e n t i f i ca r e a sc a n b ed e s c r i b e db y d i f f e r e n t i a l e q u a t i o n s ，s u c ha s i nt h e o r e t i c a lp h y s i c s ，n u c l e a rp h y s i c s ，p h y s i c a l c h e m i s t r y , q u a n t u mc h e m i s t r ya n d m o l e c u l a rp h y s i c s ，w h i c hc a nb es t u d i e db y s c h r 6 d i n g e re q u a t i o n h o w e v e r , o n l yaf e wo ft h e s ed i f f e r e n t i a le q u a t i o n sc a nb e s o l v e de x a c t l y a p p r o x i m a t em e t h o d sa r et h em a i nm e a n sf o rs o l v i n g ，a n a l y z i n ga n d u n d e r s t a n d i n gp h y s i c a lp r o b l e m s w i t ht h ei n v e n t i o na n dd e v e l o p m e n to fe l e c t r o n i c c o m p u t e rs i n c et h e19 4 0 s ，t h en u m e r i c a lc o m p u t a t i o na san e ws c i e n c eh a sg r e wu p a n dp l a y e dm o r ea n dm o r ei m p o r t a n tr o l ei nt h ev a r i o u ss c i e n t i f i cf i e l d s p a r t i c u l a r l y i nr e c e n td e c a d e s ，s o m ef a m o u sm a t h e m a t i c a ls o f t w a r eh a sb e e nd e v e l o p e d ，s u c h 弱 m a t h e m a t i c aa n dm a t l a b ，w h i c hm a k et h es c i e n t i s t st ob ea b l et ou s eac o m p u t e rt o s o l v et h ed i f f e r e n t i a l e q u a t i o n sa c c u r a t e l y a n d e f f i c i e n t l y t h i s d i s s e r t a t i o n c o n c e n t r a t e st oan e wk i n do ft h eh i g l l - d e r i v a t i v em e t h o d sf o rs o l v et h ed i f f e r e n t i a l e q u a t i o n s ，o rt h eo b r e c h k o f fm e t h o d s ，a n do n l yt h et w o - s t e po b r e c h k o f f m e t h o df o ra o n e - d i m e n s i o n a ls c h r 6 d i n g e re q u a t i o nw i l lb ed i s c u s s e d i ti sf o u n dt h a tt h el o c a l t r u n c a t i o ne r r o rc a nb eg r e a t l yr e d u c e d ，i ft h eo d d - o r d e rd e r i v a t i v e sa r ea d d e dt ot h e t r a d i o n a ld i f f e r e n c ee q u a t i o nw h i c hc o n t a i n so n l yt h ee v e n o r d e r b yu s i n gt h en e w m e t h o d ，o n ec a no b t a i nt h en u m e r i c a ls o l u t i o n so fao n e - d i m e n s i o n a ls c h r 6 d i n g e r e q u a t i o nm o r ea c c u r a t e l ya n de f f i c i e n t l y t h i sp a p e rc o n s i s t so ff o u rc h a p t e r s t h ef i r s to n eo u t l i n e st h eh i s t o r yo ft h e s c h r 6 d i n g e re q u a t i o n c h a p t e r2 ，t h eb a s i c so ft h en u m e r i c a lm e t h o d sf o r t h eo r d i n a r yd i f f e r e n t i a l e q u a t i o n si sd i s c u s s e d ，w h i c hi n c l u d et h ee u l e rm e t h o d ，t h er u n g e - k u t t am e t h o d ，t h e n u m e r o vm e t h o d ，t h es t 6 r m e r - c o w e l lm e t h o da n dt h eo b r e c h k o f fm e t h o d i nc h a p t e r3 ，an e wk i n do ft w o - s t e po b r e c h k o f fm e t h o d si sp r e s e n t e d ，w h i c hi s t h em a i nt o p i co ft h ed i s s e r t a t i o n i n c h a p t e r4 ，t h r e en u m e r i c a le x a m p l e sf o rs o l v i n g t h eo n e d i m e n s i o n a l s c h r 6 d i n g e re q u a t i o nw i t ht h ew e l l - k n o ww o o d s s a x o np o t e n t i a l ，m o r s ep o t e n t i a l a n dp 6 s c h l - t e l l e rp o t e n t i a la r ep r e s e n t e d ，a n dt h ec o n c l u s i o ni sd r a w n 上海大学硕士学位论文 8 a l lt h ew o r ka c c o m p l i s h e dd u r i n gg r a d u a t ep e r i o dh a sb e e nw e l lp r e s e n t e di nt h i s k e yw o r d s ：o n e - d i m e n s i o n a ls c h r f i d i n g e re q u a t i o n , n u m e r i c a l d i f f e r e n c e m e t h o d ，p - s t a b l e ，o b r e c h k o f fm e t h o d 上海大学硕士学位论文 1 0 序言 原子的微观结构问题是原子、分子研究领域里最重要、最基本的问题。随着 科学的发展，人们发现物质宏观特性和微观结构之间存在深刻的联系，并且认识 到原子结构是研究复杂体系( 包括分子和固体系统) 的基础。早在量子力学发展 初期，就有不少科学家着手原子结构的研究和具体计算。s c h r 6 d i n g e r 用原子波 动方程求出了氢原子的解析解，但多电子原子的结构问题却无法用解析推导的办 法给以解决，而只能采用计算物理的手段处理【l 】。 应用s c h r 6 d i n g e r 方程研究原子结构时，一般认为原子核是静置于坐标原点 的质点。对于氢原子来说，核外只有一个电子，因而势函数具有最简单的形式， 此时氢原子的s c h r 6 d i n g e r 方程能够得到严格的解析解。除此之外，其他原子都 不能求出严格的解析解。随着计算机的出现和发展，微分方程的数值方法得到了 前所未有的发展和应用。今天，计算机的飞速发展正把数值方法推向人类科学活 动的前沿，使它上升为一种重要的科学方法。值得注意的是，人类的计算能力决 不是单方面依赖电子计算机，还取决于数值方法的稳定性、精度和效率【2 】。本 文主要阐述针对s c h r 6 d i n g e r 方程提出了一种新的o b r e c h k o f f 两步方法。 在叙述此工作以前，先对以前的方法作简单的介绍。这些方法包括e u l e r 方 法、r u n g e k u t t a 方法、n u m e r o v 方法、s t 6 r m e r - c o w e l l 方法和o b r e c h k o f f 方法， 以及对一些有重要影响的工作作简单讨论，包括：e x p o n e n t i a lf i t t i n g 方法 【5 8 - 【5 9 】，e m b e d d e 方法【6l 】，e i g h t - o r d e rf o r m u l a 6 2 ，o b r e c h k o f f 方法【18 】- 【2 0 】， 通过这些方法的介绍我们可以看到，计算常微分方程的数值方法是如何从低精 度、低效率，到高精度、高效率不断向前发展的。 另外，数值方法的研究已有长足的发展。从1 7 世纪开始，科学家就开始对 物理问题采用数学方法进行处理，特别是对天体问题的二阶常微分方程的成功积 分。有许多天体物理的数值技术和数据可以追溯到高斯的工作，后来经过 s t 6 r m e r ( 1 9 0 7 ) 和c o w e l l ( 1 9 1 0 ) 整理和发展，称为s t 6 r m e r - c o w e l l 方法，即线性多 步方法。r u n g e k u t t a 方法的问世又给常微分方程数值解法的研究方向增添了新 的活力。随着2 0 世纪中期电子计算机的问世，大型、高速计算机的出现，利用 上海大学硕士学位论文 数值方法解决实际问题的需求变得日益迫切。在各个自然科学和社会科学研究领 域，有大量的问题依赖于数值方法，并且所要解决的问题越来越复杂，原有的数 值方法已无能为力。计算机的发展和计算方法的发展推动了计算机语言的发展， 特别是革命性的数学软件的问世，为数值方法的发展铺平了道路，使数值方法向 着高精度、高效率、高稳定性的道路发展。 在过去的几十年里，数值方法求解微分方程一直吸引着大家的注意 【3 1 3 ，1 8 - 3 0 ，c h a w l a 和n e t a 发展了两步四阶的p 稳定方法 3 1 】，j a i n 3 2 改进了 s t i e f f e l b e t t i s 方法， a n a n t h a k r i s h n a i a h 发展了p 稳定o b r e c h k o f f 方法 3 3 3 5 】， s i m o s 发展了三角拟合的p 稳定o b r e c h k o f f 方法 3 6 】。w a n g 研究发现一阶公式在 o b r e c h k o f f 方法中的重要性，通过高精度一阶公式的引入大大地提高了 o b r e c h k o f f 方法的精度 3 7 】，并且还发展了另一种新的三角拟合的方法【3 8 】，简 化了s i m o s 的p 稳定o b r e c h k o f f 方法的系数【3 6 】，提高了小步长时方法的稳定性， 同时又构建了其他o b r e c h k o f f 多步方法来提高效率 3 9 】【4 0 】。随着数值方法精度 的提高，我们对数值方法的要求也越来越高，因此我们迫切需要发展并改进原有 的数值方法来获得更高精度的数值解。我在原有o b r e c h k o f f 方法的基础上进行 改进，发展了一种新的o b r e c h k o f f 方法，对于一般的常微分方程，可以获得很 高的精度，为数值方法的发展提供了很好的判断标准，促进了数值方法的发展。 在本文中我们主要阐述新的o b r e c h k o f f 两步方法。 本文的结构安排如下：第一章介绍s c h r 6 d i n g e r 方程历史；第二章主要介绍 线性差分方法、常微分方程的几种数值解法以及稳定性理论；第三章详细介绍新 的o b r e c h k o f f 两步方法，它具有连续高阶微商的特点；第四章举例说明新方法在 精度上的优越性。最后对本人的工作给出结论。 上海大学硕士学位论文1 2 第1 章一维s c h r i i d i n g e r 方程 本章主要介绍一维s c h r 歌l i n g e r 方程，以及对微分方程的数值解法的发展做 简要概括。 1 1 s c h r o d i n g e r 方程 s c h r 6 d i n g e r 在1 9 2 6 年提出的波动方程，它成功地解决了量子态怎样随时间 演化以及在各种具体情况下如何求出波函数的问题。通常我们用波函数( 芦，) 来 描述一个微观粒子的量子态，当l f ，( 尹，) 确定以后，粒子的任何一个力学量的平均 值以及它取各种可能测值的概率都能够完全确定。 对于经典自由粒子来说，它的能量动量关系是 e ：卫，( 1 1 1 ) 2m 。 、 m 是粒子质量。而对于量子力学中自由粒子，它们具有波粒二象性，因此我们对 于经典的能量动量关系作如下替换： 砌昙 ( 1 1 1 12 ) 西( 2 ) 多_ p = - i 壳v ， 然后作用于波函数上，就得到以下方程 f 嗉听力= 一v 2 v ( 矾 ( 1 1 3 ) 进一步考虑在势场矿( 芦) 中运动的粒子，按照经典粒子的能量关系式 _ e ：二+ y ， 2m 对于上式作替换( 1 1 3 ) ，然后作用于波函数上，即得 f 壳昙孵，f ) = 【一：2 - 二所v 2 + v l 孵力 ( 1 1 4 ) ( 1 1 5 ) 上海大学硕士学位论文 1 3 这就是s c h r o d i n g e r 方程，它揭示了原子世界中物质运动的基本规律。在1 9 2 6 年以后的几年中，原子结构和在原子水平上的物质结构，以及究竟是什么东西决 定物质的物理和化学性质这些古老而基本的问题一个接一个地迅速得到解决。例 如，物体为什么有绝缘体、半导体和导体之分。原则上的原子结构、原子辐射和 口衰变现象等很多曾是一团迷的物理现象都搞清楚了，所以说s c h r 6 d i n g e r 方程 是量子力学最重要最基本的方程，其地位与n e w t o n 方程在经典力学中的地位相 当 6 3 】。 由于( 1 1 5 ) 式是含时s c h r 6 d i n g e r 方程，并且只含波函数对时间的一次微商， 当给定体系的初态波函数杪( 尹，0 ) 后，求解方程，原则上即可确定以后任何时刻 f 0 的波函数少( 尹，f ) ，即波函数表示的是量子态随时间演化的规律。一般来讲， 这个初值问题的求解是比较困难的，因此我们主要考虑常见的而且极其重要的情 况，即势场不显含，( 在经典力学中，这种势场中的粒子是机械能守恒的) 。此时， s c h r 6 d i n g e r 方程存在下列形式的特解，即i f ，( 尹，f ) 可以分离变量， y ( 尹，f ) = 缈( 芦) 厂( f ) ( 1 1 6 ) 代入( 1 1 5 ) 式，分离变数后，得 而ih 石a f = 赤卜砑h 2v 2 川硼孵) = e ( 1 1 7 ) e 是既不依赖于t ，也不依赖于尹的常数。这样， 丢l i l 儿) = 一i i e ， ( 1 1 8 ) 所以， 厂( f ) 口e x p - ie t h 】 ( 1 1 9 ) 因此，方程的解可表示为 缈( 尹，f ) = i u e f f ) e x p ( 一ie t h ) ， ( 1 1 1 0 ) ( 1 1 1 0 ) 式所描述的量子态即为定态，其中( 芦) 是满足下列方程 【一互h 脚2 v 2 + 矿扩) 】沙( 尹) = e y ( 尹)( 1 1 1 1 ) 的解，式( 1 1 1 1 ) 为不含时间的s c h r 6 d i n g e r 方程。 上海大学硕士学位论文1 4 ( 1 1 1 1 ) 式就是粒子的能量本征方程，( 尹) 是能量本征函数，相应的e 值 、 是体系的能量本征值。在求解方程时，并非对于一切e 值所得出的解y 扩) 都满 足物理上的要求，这些要求需要根据波函数的统计诠释或者具体物理情况。例如， 对于束缚态，就要求波函数少( 尹) 在无穷远处的值必须趋于零，这样，我们求解 方程所得到的e 值就是方程的本征值。下面我们简单了解一下一维运动粒子的能 量本征值和本征函数。 一维问题在数学上处理起来比较简单，较容易得出严格的结果，从而能够对 结果进行细致的讨论，量子力学体系的许多特征，都可以在一维问题中展示出来， 而且还是处理各种复杂问题的基础，因此，我们考虑一维定态s c h r s d i n g e r 方程， 其形式如下 f 喏咐) _ 【一嘉等川瑚嘲) ( 1 1 1 2 ) 因为是定态，既具有一定能量，因此波函数形式为 ( x ，f ) = ( x ) e x p ( 一iet h ) ，( 1 1 1 3 ) 我们把( 1 1 1 3 ) 式代入( 1 1 1 2 ) 式，y ( x ) 满足 【一妥等川瑚吣郴吣) ( 1 1 1 4 ) 或， 参吣) + 等旷俐吣) - 0 ( 1 1 1 5 ) 为计算方便，令竽：1 ，简化后可以得到 厅 少飞x ) = 厂( x ) ( x ) ， ( 1 1 1 6 ) 其中厂( x ) = 矿( x ) 一e 。( 1 1 1 6 ) 式就是本文中所要求解的方程形式，在不同的边界 条件和不同的y ( x ) 有不同的本征值解。 上海大学硕士学位论文 1 2 微分方程的数值解法的发展 各种物理学领域都涉及到微分方程问题，例如理论物理、粒子物理、物理化 学、量子物理和分子物理都能演绎成一维s e h r o d i n g e r 方程，然而其中只有少数 部分有解析解。最早解方程( 1 1 1 6 ) 式的数值方法是n u m e r o v 方法，尽管这个方 法只有四阶精度，但仍被广泛运用，因为它计算简便，而且有比较大的稳定区域 ( o ，6 ) 。另一方面，很多年来提高这个方法并寻求一种更准确、更快、更可靠的 算法来获得数值解已经成为最为活跃的领域d 8 - 【2 0 【4 4 】 5 5 】- 6 2 】。我们以著名的 w o o d s - s a x o n 势为例，数值方法有很多种【1 8 】- 【2 0 】 5 8 】一【6 2 】。另外，我们发现具 有挑战性的是计算能量本征值的误差总是随着能级的增加而增加，因此一直在促 使我们发展一种新的数值方法。例如，e x p o n e n t i a lf i t t i n g 5 8 - 【5 9 】，e m b e d d e d 方 法【61 】，e i g h t - o r d e rf o r m u l a 6 2 ，o b r e c h k o f f 方法 18 - 2 0 1 。从这些方法的不断发 展中可以看出，在计算常微分方程时，数值方法从低精度、低效率，到高精度、 高效率向前发展。 下一章我们将介绍数值方法中最基本的线性差分方法、稳定性理论以及常微 分方程的几种数值解法。 上海大学硕士学位论文1 6 第2 章线性差分方法及稳定性理论 本章首先介绍一下解微分方程的最基本的方法一差分方法，然后详细介绍 常微分方程的几种数值解法，包括e u l e r 方法、r u n g e - k u t t a 方法、n u m e r o v 方法、 s t s r m e r - c o w e l l 方法和o b r e c h k o f f 方法。 2 1 差分方法 差分方法是解微分方程的一种主要的计算方法。它的基本思想是： 1 )把连续的定解区域用由有限个离散点构成的网格来代替，这些离散 点称作网格的格( 节) 点； 2 )把在连续定解区域上定义的连续变量函数用在格网上定义的离散变 量函数来近似； 3 )把原方程和定解条件中的微商用差商来近似，积分用积分和来近似； 4 )原方程和定解条件近似地代之以代数方程组，解此代数方程组就得 到原问题的近似解。 此外，为了保证计算过程的可行及计算结果的正确，还须从理论上研究差分 方程组的数值性质，包括解的存在性、唯一性、稳定性和收敛性。稳定性就是指 计算过程中舍入误差的积累应保持有界。收敛性就是指当网格无限加密时，差分 解应收敛到原问题的解。 下面我们以最常见的二阶微分方程初值问题为例来引出线性差分格式，微分 方程形式如下 y ”( x ) = f ( x ，y ) ( 2 1 1 ) 方程( 2 1 1 ) 是一个一维物理问题的微分方程，首先我们将方程离散化后，求 找其解y = j ，( x ) 在一系列离散格点x o 五 x 2 r u n g e - k u t t a 方法 由于欧拉方法精度太低没有实用价值，如欲再进一步提高解的精度，可用一 种高精度的单步法叫u n g e - k u t t a 方法，简称r - k 法。它采用了间接使用泰勒 级数法的技术，设计思想就是利用微分中值定理设法在【而，矗+ 。】内多预报几个点 的斜率值，然后把它们加权平均作为平均斜率，以期望构造出更高精度的计算格 式。其中四阶r u n g e k u t t a 方法被称为经典单步方法。 其具体形式如下： 乃“= 咒+ 尝( k + 2 g + 2 ( 3 + k 4 ) o k i = 厂( 薯，儿) k = 厂( 薯+ 虿h ，只+ 宝k ) ( 2 2 4 ) 玛= 厂( 五+ 尝，乃+ 尝叠) 该方法有以下一些性质： 1 ) 精度为四阶、每递推一步需计算导数4 次； 2 ) 局域误差为o ( 6 t 5 ) ，全局误差为o ( a t 4 ) ，有相容性； 3 ) 稳定区域为( 口j 1 2 ) 2 7 8 5 ； 4 ) 由于有相容性和足够宽的稳定区间，所以是收敛的。 但这个单步方法精度较低，实用价值不大。 n u m e r o v 方法 n u m e r o v 方法是由原苏联天文学家b vn u m e r o v 在1 9 3 3 年求解天体轨道时 提出的数值方法。在许多数值问题中，该方法要比简单的r u n g e k u t t a 方法精度 上海大学硕士学位论文 高一些。下面对该方法作简单介绍。设在区间【口，b 】上求解二阶微分方程，以h 为 步长均匀划分区间 口，纠，在第1 1 个格点上， 毛2 工+ ，n 、h ，丹= 0 ，1 ，2 ( 2 2 5 ) = 夕( 矗) 、7 y n + l - y ( 毛+ j 1 ) = 以+ j i + 丁h 2 以、i h 3 以( 3 ) 二一。：y。矗一j，：一6h h y i ：二以伪 c 2 2 6 ， 虬一l = y ( 矗一) = 一 + 虬”一以( 3 + 丢胆盖脾篆 m 7 ) + 。+ 儿- i - 2 = 乃2 以, 十西h 4 此h ，+ 羔俩，+ 。( 办8 ) ( 2 2 7 ) 对( 2 2 6 ) 式中进行两次微商计算，把得到的式子相加可以消去y 4 ，然后利用 + 。+ 以一i 一2 以= 箬( z 钉+ 1 0 五+ z 一。) ， ( 2 2 8 ) l t e ( h ) = 一赤办6 y 6 ( x ) ( 2 2 9 ) 上海大学硕士学位论文 2 1 下面对该方法作简单介绍：先对( 2 1 1 ) 积分，可得到 雎+ y ( x + 办) 一y ( x ) = y 协) + 【( x + h - t ) f ( t ，y ( t ) ) d t ( 2 2 1 0 ) 这个公式可以看成带有积分余项的泰勒公式，但我们不希望少( x ) 出现，为此， 我们将( 2 2 1 0 ) 中的h 换成一 便可得到 少( x 一 ) 一y ( x ) 一办y ( x ) + r 一( x h t ) f ( t ，y ( t ) ) d t ( 2 2 1 1 ) 将( 2 2 1o ) 式和( 2 2 11 ) 式相加，并进行化简整理可得 f x + h y ( x + 办) 一2 y ( x ) + y ( x - h ) = l ( x + - t ) ( f ( t ) + f ( 2 x - t ) ) d t ( 2 2 1 2 ) d 在( 2 2 1 2 ) 式中将厂( f ) 换成在y s ? a 吒，吒一。上的g 次插值多项式。对于不同的 x ，刀，g 可有不同的选法，因此可得到许多不同的数值方法。例如在( 2 2 1 2 ) 式中令 x = ，x + h - - x 州，q 0 便可得到s t 6 r m e r 公式 卅一2 坛+ 一，：办：兰v 朋z m = o 在( 2 2 1 2 ) 式中令x = - l ，x + h = x 。，q 2 便可得到c o w e l l 公式 q 咒- 2 y 一l + 一2 = 乃2 v ”五， m = o ( 2 2 1 3 ) ( 2 2 1 4 ) 其中( 2 2 1 3 ) 式为显式，( 2 2 1 4 ) 式为隐式，这里的系数吒我们不再讨论，详见文 献 4 3 】。 s t 6 r m e r - c o w e l l 方法是上世纪6 0 年代最流行的数值方法。但是这个方法在 求解( 2 1 1 ) 时遇到了困难，表现为数值解不稳定。 o b r e c h k o f f 方法 通常需要求解的常微分方程为y lf ( x ，y ) 和y = f ( x ，y ) ，其中f ( x ，y ) 为y 的线性或非线性函数。对于更高阶的常微分方程，可通过降阶办法变为一阶 与两阶的微分方程。1 9 4 2 年，o b r e c h k o f f 5 提出了利用高阶微商构建差分数值方 法求解常微分方程。由于高阶微商的引进大大改善了数值方法的精度，但高阶微 商的计算繁复及早期计算机的限制，没得到太大的应用。 上海大学硕士学位论文 最早的o b r e c h k o f f 方法可表示为 哆此+ ，= z a , j y n 州，- - - 1 ， ( 2 2 1 5 ) j = o i = 0 o k 取不同值对应于不同的数值方法，k = l 时为单步方法，当k 2 对应为多步方 法。根据l a m b e r t 4 6 ，误差系数随，的增加而减小的速度远比随k 变化的速度要 快，也就是说增加高阶微商比增加步数更能提高精度。相对于多步方法来说，单 步方法的稳定性更好控制。其次，单步o b r e c h k o f f 结构简单，需初值较少，方便 求解初值问题。 下面主要讨论单步o b r e c h k o f f 方法： 只+ l 一只+ 圭c ( f ) 办t ( 儿+ i ( f ) 一( 一1 ) ，虼) ：0 ( 2 2 1 6 ) j = i 图2 2 1 单步o b r e c h k o f f 方法格点图 当歹= 4 时，o b r e c h k o f f 方法的形式为： 上海大学硕士学位论文 少( x + 兰) 一y ( x i h ) + 届办( y c x + 鲁，+ y ( 石一兰，) + f 1 2 h 2 ( 少飞x + 宝) 一y ”( x j h ) ) ( 2 2 1 7 ) + f 1 3 h 3 ( y o ) c x + - h zy ( 3 ) c x 一鲁，) + f 1 4 h 4 ( y ( 4 c x + 考，一y ( 4 c x 一耋，) = 。， 其中y ( x + 尝) 口以+ 。和y ( x 一兰) 口此，通过在x i 鬟jt a y - 。r 级数展开，可求 1气11 屈，尾，8 3 ，屈) 2 吒，云，一亩，志) ， ( 2 2 1 8 ) 相应的截断误差为 l t e ( h ) ：盟1 - o ( h 1 。) ( 2 2 1 9 ) 2 5 4 0 1 6 0 0 7、7 从( 2 2 1 9 ) 可以看出，单步的o b r e c h k o f f 方法具有很小的误差系数 2 5 4 0 l 1 6 0 0 。从单步o b r e c h k o f r 方法( 2 2 1 7 ) n - ，单步方法结构简单、稳定 性好，误差系数相对较小，可以通过利用高阶微商得到较高的精度。 对于一个数值方法，除了精度以外，稳定性是必须考虑的首要因素。在实 际计算中误差是不可避免的。由于计算机字长的限制、初值精确不够、截断误差 的存在，误差会逐渐积累和传播造成不稳定。另外对有些问题，例如周期问题， 计算方法没有足够的稳定区间，这样的方法就根本不能使用。因此首先讨论数值 方法的稳定性是必要的。 2 3 稳定性理论 稳定性在微分方程的数值解法中是一个极其重要的问题。因为微分方程数值 解法是通过离散化将微分方程转化为差分方程来求解的，而在差分方程的求解过 程中，存在着各种计算误差，这些计算误差如舍入误差引起的扰动，在误差传播 上海大学硕士学位论文 过程中，可能会大量积累，以至于“淹没 了方程的真解。例如初值问题【l 】。 下面我们通过对求解二阶常微分方程的初值问题来讨论数值方法的稳定性 和引入相应的理论，二阶常微分方程的初值问题如下 y ”( x ) = f ( x ，y ) ，y ( x o ) = y o ， y ( ) = y 1 0 ( 2 3 1 ) 早在上世纪6 0 年代，当时许多数值方法很多都是不对称的，其中线性多步 方法中精度较高且最为流行的当属s t s r m e r - c o w e l l 方法。但这个方法在积( 2 3 1 ) 时遇到了困难，其数值解是不稳定的。较为有名的例子是计算天体问题的两体轨 道，发现当s t s r m e r - c o w e l l 的多步方法的步数超过2 以后，在数值计算中就会表 现出不稳定【4 4 】。1 9 7 6 年l a m b e r t 和w a t s o n 4 4 首先认识到在求解周期问题时， 数值方法存在稳定性问题，提出了的对称多步方法的稳定区域的概念。并指出造 成s t s r m e r - c o w e l l 方法不稳定的原因是该方法的稳定区域等于零。对于已知频率 的周期问题，一个数值方法的稳定性质可由步长与频率乘积来描述，即所属的稳 定区间的宽度。下面引用l a m b e r t 的文章【4 4 】对稳定性理论作简单介绍。 对于( 2 3 1 ) 式，l a m b e r t 提出了如下对称的多步方法 kx q ，t _ h 2 层z “， k 2 ( 2 3 2 ) t = 0t = 0 系数有如下对称性质 = 一，层= 厥中 i = o ，1 k 则方法( 2 3 2 ) 的数值性质可由一组特征多项式函数( p ，仃) 来表示 ki p = f ，仃= 屈f ， f c ( 2 3 3 ) ( 2 3 4 ) 我们称相应的多步方法( 2 3 2 ) 为方法( p ，仃) 。对于充分光滑的任意连续函数 z ( x ) ，如果方法( p ，仃) 有 kk z + f 厅) 一 2 f l , z ( x + i h ) = c p + 2 h ，+ 2 z ，+ 2 o ) + d ( 厅肿3 ) ， ( 2 3 5 ) j = 0j = o 则称方法( p ，仃) 为p 阶，其误差系数为q + 2 。 如果特征多项式函数( p ，盯) 同时满足以下条件 上海大学硕士学位论文 七 ( 1 ) a k = 1 ，l i + l p oi o ，i 屈i o ； ( 2 ) p ，仃无公共因子； ( 3 ) 户( 1 ) = p 。( 1 ) = o ，p ”( 1 ) = 2 0 - ( 1 ) ，即( p ，盯) 的阶数至少大于一； 则称方法( p ，o r ) 为零稳定，即特征多项式p 所有的根都满足i 缶l 1 。 根据h e n r i c i 4 7 ，发现方法( p ，仃) 是收敛的，并且如果数值方法的特征多 项式函数p 有两个根为lf ，。：| - 1 ，其它根为i 缶i 1 。并且特征多项式函数p 的 根可表示为，s = 1 ，2 ，k ，其中缶= 乞= + l 为主根，乒，s = 3 ，4 ，k ，为伪根。 考虑如下测试方程 其通解为 y ”( x ) = 一彩2 少( x ) ，国，y ，x r( 2 3 6 ) y ( x ) = ac o s ( o x + b s i no ) x ，或少( x ) = a e 蛔善+ b e 叫口工( 2 3 7 ) ta ：b ：o ，对于所有的彳和b ，其周期为丝。 彩 将方法( p ，仃) 应用于测试方程( 2 3 6 ) ，得到 足 ( q + 日2 屈) ，= o 其中日= c a b 。 ( 2 3 8 ) 根据递推性质，多步方法的数值结果可表示为以= 4 。其中， r a s = l ，2 ，七) 表示方程不同的根。因此( 2 3 8 ) 相应的特征多项式( 尸，仃) 可表示 为 q ( ，；日2 ) = p ( ，) + 日2 仃( ，) ( 2 3 9 ) 如果日2 ( o ，h 0 2 ) ，特征多项式q ( ，；日2 ) = 0 的解满足： - - e 一，r 2 = e 一坩，i ，：i 1 ，i = 3 ，4 ，k ， ( 2 3 1 0 ) 上海大学硕士学位论文 其中口( 日) 为实函数，则区间( 0 ，h 0 2 ) 称为方法( 2 3 2 ) 的稳, i g l x n 。 因此根据l 锄b e r t 和w a s t o n 4 4 ，在区间( o ，h 0 2 ) 内，方法( 2 3 2 ) 具有周期 性质，适合于求解周期性问题。 我们在分析具体问题时可以用更为简便的方法，考虑多步方法 将测试方程y 。= 缈2 y 代入方法，直接令咒卅寸e ( n + i ) l a 9 he 如 一五，i = 0 ，1 ，k 得 到特征多项式 r q ( 五；日2 ) = ( + 日2 屈) ， ( 2 3 1 2 ) 其中h = c o h 。 解特征多项式q ( a ；日2 ) = 0 ，同样可以得到上式的解旯满足条件( 2 3 1 0 ) 的稳定 区间，即( 0 ，风2 ) 。 如果一个数值方法的稳定区间为( 0 ，h 2 ) 为( 0 ，0 0 ) ，则该方法称为p 稳定的方 法。 下面用这个理论分析单步o b r e c h k o f f 方法的稳定性。可以证明单步 o b r e c h k o f f 方法具有p 稳定性质。 将单步o b r e c h k o f f7 h - 法( 2 2 1 7 ) ，代入测试方程少什= 一( 1 9 2 y ，可得到相应的 特征多项式 p ( z ) = ( 彳( 日) + ，曰( h ) ) 名一( 彳( 日) 一f b ( h ) ) = 0 ， ( 2 3 1 3 ) 其中系数为 彳( 日) = 1 一h 2 屈+ h 4 以， b 悼、) = hp i h ：p 3 q 。3 1 4 ) 特征多项式( 2 3 1 3 ) 的根可表示为 a = e t a ( m ( 2 3 1 5 ) 满足p 稳定条件，即单步o b r e c h