（凝聚态物理专业论文）关于线性obrechkoff多步方法求解非线性振动方程的研究.pdf

上海大学硕士学位论文 v 摘要 现代科学、技术、工程中大量数学模型都可用微分方程描述，但只有少数问 题有解析解。整个2 0 世纪，关于如何得到周期二阶微分方程数值解的研究得到 了长足发展。特别是，随着上世纪4 0 年代电子计算机的发明，更多新思想与新 方法得到实现。近几十年，一些数学软件，例如m a t h e m a t i c a ，m a t l a b 帮助研究 者做出了更深入的分析、实现了更复杂的方法。 本硕士论文开展的是：周期二阶微分方程的数值方法的研究。 通过改进原有的数值方法，发展了一种新的高精度、高效率的p 稳定 o b r e c h k o f f 线性多步方法。通过高阶微商的使用来提高方法的精度，并在计算 过程中采用n e w t o n 线性化大大简化高次代数方程组的求解。利用这种方法可以 对常见的周期初值问题进行求解数值求解，例如s t i e f e l - b e t t i s 问题和d u f f i n g 方程，都能得到高精度并且稳定的数值解。为了显示这类新方法在精度和效率方 面和以前的方法相比，具有压倒性的优势，本文中选用了大量的数值例子进行对 比。从计算结果来看，新方法的精度要比以前的方法高好几个数量级。 本论文是基于作者硕士期间所做工作和发表论文所作深入和全面的总结。 关键词：数值差分方法，p 稳定，精度，o b r e c h k o f f 方法，n e w t o n 线性化 上海大学硕士学位论文 v i a b s t r a c t i nv a r i o u sa r e a so fm o d e ms c i e n c e s ，e n g i n e e r i n ga n dt e c h n o l o g ya p p l i c a t i o n s ， m a n yp r a c t i c a lp r o b l e m sl e a dt ot h em o d e l s ，w h i c hc a nb ed e s c r i b e db ys e c o n d o r d e r o r d i n a r yd i f f e r e n t i a le q u a t i o n s h o w e v e r , o n l yq u i t eaf e wo ft h e mh a v et h ee x a c t l y a n a l y t i c a ls o l u t i o n s t h r o u g ht h ew h o l ec e n t u r yo f2 0 “，h o wt oo b t a i nt h en u m e r i c a l s o l u t i o nf o rap e r i o d i co d e sh a sa t t a i n e db r i l l i a n tp r o g r e s s p a r t i c u l a r l y , w i t ht h e i n v e n t i o no fe l e c t r o n i cc o m p u t e rs i n c et h e1 9 4 0 s ，m o r ea n dm o r en e wi d e a sa n d a p p r o a c h e sh a v eb e e nr e a l i z e d i nr e c e n td e c a d e s ，s o m en e wd e v e l o p e ds o f t w a r ef o r m a t h e m a t i c ss u c ha sm a t h e m a t i c a ，m a t l a bh a v eo p e n e dan e w e p o c h , w h i c ha l l o w su s t om a k ed e e pt h i n k i n g ，a n a l y s i s ，p r o g r a m m ea n di m p l e m e n tc o m p l e xf l g o r i t h mi np c e n v i r o m e n t t h i sd i s s e r t a t i o ni sa c c o m p l i s h e do nm yw o r ka n dt h ep u b l i s h e dp a p e r s ，w h i c h i n c l u d et h en u m e r i c a lm e t h o d sd e v e l o p e df o rt h ep e r i o d i cs e c o n d - o r d e ro r d i n a r y d i f f e r e n t i a le q u a t i o nw i t ht h ei n i t i a l - v a l u e b ya d d i n gt h eh i g h o r d e rd e r i v a t i v e sa n da p p l y i n gt h en e w t o nl i n e a r i z a t i o nt o t h er e s u l t a n th i 曲o r d e ra l g e b r a i ce q u a t i o n s ，ih a v eg e n e r a l i z e dt h et r a d i t i o n a lp - s t a b l e l i n e a rs y m m e t r i c a lm u l t i s t e pm e t h o d st ot h ep - s t a b l el i n e a rs y m m e t r i co b r e c h k o f f m u l t i s t e pm e t h o d s , w h i c hh a v ev e r yh i g hl o c a lt r u n c a t i o ng l t o r t h en u m e r i c a l e x p e r i m e n t sf o r t h ew e l l k n o w np e r i o d i cp r o b l e m s ，s u c ha st h es t i e f e l - b e t t i sp r o b l e m a n dd u f f m ge q u a t i o nd e m o n s t r a t et h a to u rn e wm e t h o dh a st h ea d v a n t a g eo v a l t h e p r e v i o u sm e t h o d si ns t a b i l i t y ，a c c u r a c ya n de f f i c i e n c y a l lt h ew o r ka c c o m p l i s h e dd u r i n g o - 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 n t h i sd i s s e r t a t i o n k e yw o r d s ：n u m e r i c a ld i f f e r e n c em e t h o d ，p - s t a b l e ，a c c u r a c y , o b r e c h k o f f m e t h o d , n e w t o nl i n e a r i z a t i o n 上海大学硕士学位论文 原创性声明 本人声明：所呈交的论文是本人在导师指导下进行的研究工作。 除了文中特别加以标注和致谢的地方外，论文中不包含其他人已发表 或撰写过的研究成果。参与同一工作的其他同志对本研究所做的任何 贡献均已在论文中作了明确的说明并表示了谢意。 签名：至鳌! 玺e t 期：竺盟! ：12 本论文使用授权说明 本人完全了解上海大学有关保留、使用学位论文的规定，即：学 校有权保留论文及送交论文复印件，允许论文被查阅和借阅；学校可 以公布论文的全部或部分内容。 ( 保密的论文在解密后应遵守此规定) 签名：逖导师签名：礁墨茸丞日期：迎妇：盘、! 上海大学硕士学位论文 。 第一章绪论 现代科学、技术、工程中大量数学模型都可用微分方程来描述，很多现代自 然科学本身就是建立在微分方程的基础上。从微积分理论形成以来，微分方程就 成为人们描述、解释或预见各种自然现象的基本手段并不断取得了显著的成效。 遗憾的是，绝大多数微分方程定解问题的解是不能用解析形式表示的，这就产生 了理论进一步发展的要求。通过大量数学模型建立了反映客观现象的各种微分方 程，但人们又无法得到这些方程的解，就不能准确和定量解释或预见客观过程 【l 】 随着计算机的发展和出现，为了解决上述矛盾的一门新的学科一微分方程的 数值方法得到了前所未有的发展和应用。今天，计算机的飞速发展正把数值方法 推向人类科学活动的前沿，使它上升为一种重要的科学方法，特别是对科学的定 量化起了重要作用。 微分方程数值方法的主要任务就是求解各种各样的微分方程，特别是一些大 规模、非线性、几何非规则性的方程。例如由物理学中的一些基本原理，质量守 恒，能量守恒，动量守恒，角动量守恒，电荷守恒，引力规律，牛顿第二定律， 电磁规律等得到的数学模型主要是用微分方程一常微分方程和偏微分方程来描 述的。因此掌握和应用微分方程数值方法的已不仅限于计算科学专业，大量从事 力学、物理学、天文学、航空、航天、土木、机械、电子、地质勘测等领域的技 术人员，甚至连金融工程、风险投资等非科技领域都已把这门学科作为自己领域 的一种重要研究手段和方法。 值得注意的是，人类的计算能力决不是单方面地依赖于电子计算机，还取决 于计算方法的稳定性，精度和效率。数值方法的发展对于提高计算能力与研制新 一代计算机是同样重要的。在今天，即使计算机发展速度之快让人预料不及，但 实际应用需求增长更快，大规模科学与工程计算仍得不到满足，需依靠发展新的 数值方法来完成这些规模庞大的计算。 过去几十年里，数值方法的研究已有长足的发展。从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 o r m e r - c o w e l l 方法。另外1 9 世 上海大学硕士学位论文 2 纪末b a s h f o r t h 和a d a m s 首次提出线性多步方法以及r u n g e k u t t a 提出单步方法， 求解常微分方程的数值方法就成为充满活力的研究领域。随着2 0 世纪中期电子 计算机的问世，大型、高速计算机的出现，利用数值方法解决实际问题的需求变 得日益迫切。在各个自然科学和社会科学研究领域，有大量的问题依赖于数值方 法，并且所要解决的问题越来越复杂，原有的数值方法已无能为力。计算机的发 展和计算方法的发展推动了计算机语言的发展，特别是革命性的数学软件的问 世，为数值方法的发展铺平了道路，使数值方法向着高精度、高效率、高稳定性 的道路发展。 本文主要研究具有周期性质的二阶常微分初值问题 y ”( 石) = f ( z ，y ) ，y ( x o ) = y o ，y ( x o ) = y o ，( 1 1 ) 的数值解法，其中f ( x ，y ) 可含有y 的非线性项。 方程( 1 1 ) 在科学和工程领域都有广泛的应用，一直以来，求解方程( 1 1 ) 的数 值方法吸引着大家的注意睁1 5 ，3 5 5 7 ，c h a w l a 和n e t a 发展了两步四阶的p 稳定 方法【1 6 】，j a i n 1 7 改进了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 方法 2 2 - 2 4 ，s i m o s 发展了三角拟合的p 稳定o b r e c h k o f f 方法【2 0 】。 w a n g 研究发现一阶公式在o b r e c h k o f f 方法中的重要性，通过高精度一阶公式的 引入大大地提高了o b r e c h k o f f 方法的精度【1 8 】，并且还发展了另一种新的三角拟 合的方法【1 9 】，简化了s i m o s 的p 稳定o b r e c h k o f f 方法的系数【2 0 】，提高了小步 长时候方法的稳定性【5 l 】，同时又构建了其他o b r e c h k o f f 多步方法来提高效率 【2 9 【3 4 】。 下面我们介绍几种物理学中常见的二阶非线性常微分方程【2 】。 1 d u f f i n g 方程 d u f f i n g 方程的形式为 工+ 2 ，+ 国孑x + 占，孑x = 4c oso t ( 1 2 ) 它是描述在双势阱中受迫振动系统的一种动力学方程，具有非线性行为，并 且只有在一定条件下有周期解。物理学家对这样的系统有极大的兴趣，因为在一 定的初值条件和参数下会产生不可预测的动力学行为一一混沌现象。因此如何有 效地得到的数值解也是计算科学工作者感兴趣的问题。 上海大学硕士学位论文 与线性谐振子相比，它除了有线性恢复力一矿x 和与速度成正比的阻尼力 一2 ；以外，还有一个由双势阱产生的非线性恢复力一s 廊，3 ，其中，风为正常 数，h t l 。当s ，o 时，它表征非线性恢复力大于线性恢复力，称为硬非线性； 当sc0 时，它表征非线性恢复力小于线性恢复力，称为软非线性。一c o sq f 为强 迫振荡，其中_ q 为常数。 无阻尼且无外驱动力的d u f f i n g 方程，通常描写弹性体的非线性振动，形式 如下： 簟+ f v o z x + f 卢孑工3 = 0 ( 1 3 ) 无阻尼且有强迫振动项的d u f f i n g 方程， j + 磊工+ 占席x 3 = a c o s c o o t ( 1 4 ) 这种方程在一般情况有周期解，并且解的周期与外驱动力的频率相同。经过许多 人的研究，发现这个方程有解析解，并且得到精度很高的解析解。所以在数值研 究方面，很多人喜欢用这个方程来检验、比较新的数值方法的稳定性、精度和效 率。本文也用该方程来检验我们发展的新方法。 2 v a n d e r p o l 方程 为了分析三极管的振荡电路的参数对稳定性的影响，v a n d c r p o i 提出了一个 数学模型，可用微分方程来描述。该方程称为v a nd c rp o l 方程。 这个方程也可用力学模型来解释：设质点的位移为x ，其加速度和速度分别 为叠和l ，若它受到一个线性恢复力面z ( 嘞为线性振荡圆频率) ，同时还 受到一个非线性阻尼力- 2 ( 言一1 ；，其中，。称为阻尼系数，是一个正常数。 这样，表征质点运动的方程为 号一- 卜币一 s ， 显然，质点的运动是与阻尼项z ( 孝 一 系数的符号决定的。当 i 叫c 时表征负阻尼，当h ，时表征正阻尼。也就是说，当位移增加并且超过 上海大学硕士学位论文 4 临界值时，阻尼机制的性质将由负变为正。因此在质点的位移小于临界值，就会 产生自发振荡并且维持这种状态，振幅不断增加；直到达到临界值，电路变成衰 减振荡，振幅不断减少，直到小于临界值。v a nd e rp o l 方程在一般情况下都有 周期解。这个方程用来分析真空电子管振荡电路非常有效。 3 单摆运动方程 无阻尼的单摆运动方程 设有竖直平面内运动的一单摆，摆球为单位质量，摆长为i 。若摆球由初始 位置岛开始摆动，摆动的角度为0 ，则在重力作用下的单摆运动方程如下： 扩+ ；s i np ) = 0 ( 1 6 ) 其中j 表示角加速度。 有阻尼的单摆运动方程 如果考虑单摆运动受到的阻尼力，且阻尼的大小与角速度成正比，这样，单 摆运动方程可以写为 + 2 百+ 曲；s i n ( 0 ) = 0 ( 1 7 ) 其中，0 称为阻尼系数，_ 2 称为阻尼力 在本文中我们主要把两类方程作为例子，一类是类似检测方程的 州x ) - - - - j y ( ) 吲习，另一类是非衰减的d u f f i n g 方程，但我们的方法对普遍的二 阶非线性方程都适用。 本文结构安排如下：下一章，介绍解微分方程最基本的数值方法一一差分方 程，常见的结构、l a m b e r t 稳定性理论以及p 稳定方法。第三章，给出我们发展 的p 稳定高精度方法。第四章，通过用新方法得到的数值结果，展示新方法的优 越性。最后给出总结。 上海大学硕士学位论文 第二章用线性差分方法解周期问题 2 1 差分方程 解微分方程度基本方法是采用差分方程。将一个微积分问题离散化，转换成 代数问题即差分方程，就可达到解微分方程的目的。主要过程为： 1 ) 将物理问题的连续区域进行离散化； 2 ) 利用差分公式计算每个格点上的各阶导数： 3 ) 将各阶导数的差分表达式代入微分方程从而得到最后的迭代公式； 4 ) 利用迭代公式进行计算。 d 似) f 卜叫户：= ：才一 图2 - 1 一l 解的区域d ( x ) 与有限差分格点 一维物理问题连续区域d ( x ) 的离散化如图2 1 1 所示。首先将区域分成间 距为h 的离散格点，然后用离散格点的有限差分方程近似表示连续区域的微分方 程。图2 一l - l 中的函数y ( x ) 在格点，z 处可近似表示为 y ( ) = 虬( 2 1 1 ) 利用差分方法对微分方程的各阶导数进行替换，函数以功在矗点的泰勒级 数可写为： y n + l 咒+ 儿| j i + 丢以”厅2 + + 去( 埘) 矿+ r ”1 ， ( 2 1 2 a ) ：儿一只t 厅+ yn 矗2 + + 生孚以佃) 矿+ 月一i ( 2 1 2 b ) zl! 把( 2 1 2 b ) 与( 2 1 2 a ) 相减，忽略址的平方项及高阶项，得到一阶微分的中心差商 为： 驴鼍产 ( 2 1 3 ) 把( 2 1 2 b ) 和( 2 1 2 a ) 相加，忽略，的三次方项及高阶项，得到二阶微分的中心差 上海大学硕士学位论文 商为： 以。= 丛等立 ( 2 1 4 ) 根据( 2 1 2 a ) 至( 2 1 4 ) 可构造出如下微分方程形式： y 。= f ( x ，y ) ( 2 1 5 a ) 的差分格式 ) ，肿1 + ) o 一2 = 炉z 。( 2 1 5 b ) 在上面的差商表示中，如果高阶微商可以展开到第m + l 项时，其截断误差为 旷| - 丽1 ) ，1 恸“ ( 2 1 6 ) 在大多数情况下，截断误差随h 哼0 变为零，阶数越高则误差越小，收敛 越快。 在实际计算中产生误差是不可避免的，一方面计算机字长限制带来误差，称 为舍入误差；另外数值方法或所用的差分方程本身不精确带来的误差，称为截断 误差，最后初始值不精确也可能引入误差。这样经过逐次计算，这些误差就会积 累并传播下去，对以后的计算结果都精度产生影响，有可能使后面的结果越来越 差。 圭堂奎堂堡主堂鱼堡茎 1 2 2 常见的数值方法 1 n u m e m v 方法 n u m e r o v 方法是由原苏联天文学家b v n u m e r o v 在1 9 3 3 年求解天体轨道时 提出的数值方法。在许多数值问题中，该方法要比简单的r u n g e k u t t a 方法精 度高一些。下面对该方法作些简单介绍。设在区间【口，6 】上求解二阶微分方程， 以h 为步长均匀划分区间【a ，b 】，在第玎个格点上， 靠= x + n h ，栉= o ，1 ， = y ( x d ( 2 2 1 ) 甄俐分别珂y m - l , 进，仃黍勒级数展开，得到 - = y ( x n + h ) = y 。+ b y ；+ 譬y ：+ i h 3y ，卜 + 百h 4y h 岳j ，辨而h 6 ) + d ( ，) ， 琢，叫矿j j ) = y n - - b y ；+ 譬小譬y ：3 ) 渊) + 万h 4y 卜盖5z 5 ，+ 面h 6 + 口( 厅，) 然后g t ( 2 2 2 ) 式中的两式相加，得到 + 一轨= 扳+ h 1 2 4 ，v 。( 4 ) + 兰妒+ 耐) ( 2 2 3 ) 由( 2 1 5 b ) 可知，( 2 2 3 ) 的左边是y 的二阶中心差商。对( 2 2 3 ) 两边对石微商两次， 可得如下的y 。二阶中心导数， 2 战+ 。+ 矗。一2 以) = 矗4 蟛+ 等磴+ 。( h 8 ) ( 2 2 4 ) 用( 2 2 3 ) 式和( 2 2 4 ) 式消去) ，项，并利用( 1 1 ) 式，将b 州) 写成，进行整 理之后可得n u m e r o v 方法的差分公式 y n + l + 一1 一巩= 西h 2 ( 矗l + 1 0 厶+ 丘1 ) ( 2 2 s ) 上海大学硕士学位论文 l t e ( h ) 一赤矿y 6 ) ( nv n u m e r o v 方法是经典的两步方法，递推公式非常简洁，只要知道了前两个节 点的数值，就可以利用( 2 2 5 ) 式求该点的数值，从而一步一步向前推进，求得后 面每一个节点的数值。该方法很实用，但是有一定的局限性，就是精度还不够高。 2 s t s r m e r - c o w e l l 方法 s t 6 r m e r - c o w e l l 方法是天体物理学家非常熟悉的数值方法，m i l a n i 等人在 1 9 8 6 1 9 8 8 年曾用该方法对太阳系外的行星轨道进行数值积分，得到长时间( 1 0 8 年) 的轨道跟踪数据【5 8 】。 下面对该方法进行简单介绍：先对( 1 1 ) 式积分两次，可得到 y ( x + ) 一y ( x ) = 砂x x ) + e + 6 ( z + h - t ) f ( t ，y ( t ) ) d t ( 2 2 7 ) 这个公式可以看成带有积分余项的泰勒公式，但我们不希望，出现，为此，我 们将( 2 2 7 ) 式中的h 换成h 便可得到 y ( x - h ) - y ( x ) = 一砂( x ) + e 瑚( x 一而一t ) f ( t ，y ( t ) ) d t ( 2 2 8 ) 将( 2 2 7 ) 式和( 2 2 8 ) 式相加，并进行化简整理可得 灭x + 协一2 贝x ) + 贝x 一= f 6 + 而一f ) 盯p ) + f ( 2 x f ) ) 毋 ( 2 2 9 ) 在( 2 2 9 ) 式中将朋换成在节点h ，x 。一口上的q 次插值多项式。对于不同的 x ，珂，q 可有不同的选法，因此可得到许多不同的数值方法。例如在( 2 - 2 9 ) 式中 令x = 靠，工+ h = x l ，q 0 便可得到s t 6 r m e r 公式 j ，州一2 y 。+ y 川：| 1 1 2 盯。v m 厶 ( 2 2 1 0 ) 坍= 0 在( 2 2 9 ) 式中令x = 稚l ，x + h = 而，q 2 便可得到c o w e l l 公式 y 。一2 y 州+ y m = 办2 q 盯肼v 坍厶 ( 2 2 1 1 ) 小= 0 其中，( 2 2 1 0 ) 式为显式，( 2 2 1 1 ) 式为隐式，这里的系数o 我们不再详细讨论， 详见【3 】。 s t o r m e r - c o w e l l 方法是在上世纪6 0 年代最为流行的数值方法。但这个方法 上海大学硕士学位论文 在积( 1 1 ) 时遇到了克星，数值解是不稳定的。较为有名的例子是计算天体问题的 两体轨道，发现当s t 6 r m e r - c o w e l l 的多步方法的步数超过2 以后，在数值计算中 就会表现出不稳定性【2 7 】。 3 高阶微商方法 1 9 4 2 年，o b r e c h k o f f 【6 】提出了利用高阶微商构建差分数值方法求解常微 分方程。最早的o b r e c h k o f f 方法可表示为 k， k 哆) - - y h 岛，o 州，嚷= 1 ( 2 2 1 2 ) j = o i = o j = o k 取不同值对应不同的数值方法，k = 1 时为单步方法，当k 2 对应为多步 方法。根据l a m b e r t 2 5 ，误差系数随着，增加而减小的速度远比随k 变化的速度 要快，也就是说增加高阶微商比增加步数更能提高精度。相对于多步方法来说， 单步方法的稳定性更好控制。其次，单步o b r e c h k o f f 结构简单，需初值较少，方 便求解初值问题。 下面以单步o b r e c h k o f f 方法为例进行讨论： 摊 艇+ l 00 00 o ：_ 呻 00 图2 - 2 1 单步o b r e c h k o f f 方法格点图 当，= 4 时，o b r e c h k o f f 方法的形式为： 胁争吨一多 峭办( y 协+ 争+ y 协一争鹕序( y 瓢+ 7 h 叫w 一争 怕牙护 + 争+ 尸 一争+ 层扩 + 争尹 毫) ) = q 上海大学硕士学位论文 1 0 其中我们用到y ( x + 匐。和y ( x 一尝 叠以，通过在石点的t a y i o r 级数 展开，可求得( 2 2 1 3 ) 的系数为： 届，反，屈，尼 = 一圭，西3 ，一面1 ，而1 ， ( 2 2 1 4 ) 相应的截断误差为 三t e ( 耻轰惫硎o ) ( 2 2 1 5 ) a ( 2 2 1 5 ) 可以看出，单步的o b r e c h k o f r 方法具有很小的误差系数互；孑i 话丽。 但从数值方法的应用历史来看，这种数值方法并没有像其他一步r u n g e - - k u t t a 方法那样得到广泛应用。一个重要原因是包含高阶微商使计算大大复杂了。 对于一个数值方法，除了精度以外，稳定性是必须考虑的首要因素。在实际 计算中误差是不可避免的，由于计算机字长的限制、初值精确不够、截断误差的 存在，误差会逐渐积累和传播造成不稳定。另外对有些问题，例如周期问题，刚 性问题，计算方法没有足够的稳定区间，这样的方法就根本不能使用。因此在采 用一个数值方法时，首先要搞清这种方法的稳定性，是不是适合我们的具体问题。 上海大学硕士学位论文 2 3 稳定性理论 解特殊的二阶常微分方程( 不包含一阶微商项) 的初值问题时，如果采用 s t s r m e r - c o w e l l 方法时，当步数超过2 时数值解都会出现所谓轨道不稳定的问 题。例如对于描述在圆的轨道上做规则运动的检验问题，无论选取什么步长，用 s t s r m e r - c o w e l l 方法得出的数值解都向内盘旋转移动虽然也有各种各样的改进 的s t o r m e r - c o w e l l 方法能够克服这一问题。但是它们都需要一个预先知道的频 率。1 9 7 6 年l a m b e r t 和w a t s o n 2 7 首先认识到对于周期问题，数值方法存在稳定 性问题，提出了对称多步方法的稳定区域的概念，并指出造成s t i r m e r - c o w e l l 方法不稳定的原因是该方法的稳定区域等于零。对于已知频率的周期问题，一个 数值方法的稳定性质可由步长与频率乘积来描述，即所属的稳定区间的宽度。下 面引用l a m b e r t 的论文【2 7 】对稳定性理论作简单介绍。 对于特殊的二阶常微方程y 。( 功= f ( x , y ) ，l a m b e r t 考虑采用如下的对称多 步方法 ff q ，= h 2 尼，舯， k - 2 ( 2 3 1 ) i = 0f 1 0 其中z = ( 毛，) ，并且系数有如下对称性质 = 口f 。属= 反f ， f = o ，1 置 ( 2 3 2 ) 方法( 2 3 1 ) 相应的特征多项式是： r 户( f ) = q f ，盯( f ) = 尼f ， f c j = ol i o 我们也称方法( 2 3 1 ) 为方法( p ，仃) 对于任意充分光滑的函数z ( x ) 和充分小的h ， 若有 if o t l z ( x + i h ) - h 2 e f l f ( x + i h ) = ( 2 h 9 “z “( x ) + o ( 向”3 ) ， ( 2 3 4 ) 则称方法( p ，仃) 为尸阶，其误差系数为c 0 2 a 如果特征多项式函数( p ，盯) 同时满足以下条件 上海大学硕士学位论文 k 0 ) a k = l ，l a ol + l p oi o ，i p , l o ； j ，o ( 2 ) p ，仃无公共因子； ( 3 ) p ( 1 ) = p x l ) = o ，p ”( 1 ) = 2 a ( 1 ) ，即( p ，盯) 的阶数至少大于一； 则称方法( p ，盯) 为零稳定，即特征多项式p 所有的根都满足l t 阵l 。 f 根据递推性质，表示为数值结果为y 。= 以等，因此方法( 2 3 1 ) 的数值性 质可由一组特征多项式函数( p ，盯) 来表示。 以上讨论的是在h 争0 时的渐进稳定性，而在实际的计算过程中h 总要取一 定大小的值，为了保证方法的稳定性一绝对稳定性，步长h 的取值要受到一定的 限制。为了找到方法的稳定性条件，下面考虑把测试方程 y 。= 一( - 0 2 y ，出，y r 代入( 2 3 1 ) 式，可以得到 k ( q + 日2 乃) + ，= o ，h = h c o j = 0 再把测试方程的解 o ) = 一”代入上式，得到 k ( 口，+ 日2 局) 一州“埘= 0 ， = o 约去一一，并设e = ，可得到特征方程为 ( 口j + h2 岛) ，l0 ， ，- 0 ( 2 3 6 ) ( 2 。3 7 ) 该式的一般解为，= 4 誓，j = 1 2 丘，这里是特征方程( 2 3 8 ) 式的根， s = l 在文献 2 7 】中，l a m b e r t 详细证明了，如果存在一个区间( 0 ，风2 ) ，在该区间内的 h 使( 2 3 8 ) 式的根满足如下条件： 上海大学硕士学位论文 ，i = 肿，匕= f 州研， 悱1 ，s ：3 ，”毛r 那么，该方法的稳定区间就是( o ，风2 ) 。( 2 3 9 ) 式就是该线性多步方法的稳定性条 件。如果一个方法的稳定区间为日2 ( o ，呦，则该方法称为p 稳定的方法。 根据h e n r i c i 5 6 对稳定问题的进一步讨论，发现如果数值方法的特征多项式 函数p 有两个根为i f l 2 | _ 1 ，其它根为i 幺喀1 ，则( 岛盯) 仍为收敛方法。并且 特征多项式函数户的根可表示为，s = l 2 ，k ，其中自= 乞= + l 为主根， ，s = 3 ，4 ，k ，为伪根。 在具体分析时可以用更为简便的方法，考虑多步方法 q ) o ，= 2 层厶， ( 2 3 1 0 ) 将测试方程y _ - - c 0 2 y 代入( 2 3 1 0 ) ，直接令) 争州7 ”，如“争名，其 中i = 0 1 k ，即可得到特征多项式 r f 2 ( 2 ；h 2 ) = ( q + 日2 尼m = o 1 5 0 其中h = 缈 。 解特征多项式n ( ；h 2 ) ，同样可以得到满足条件( 2 3 9 ) 的稳定区间 ( 0 ，砜2 ) 。 自从1 9 7 6 年，l a m b e r t 2 7 提出了针对周期问题的对称多步方法的p 稳定性， c a s h 【3 0 1 ( 1 9 8 1 ) ，c h a w l a 【3 1 1 ( 1 9 8 1 ) 和h a i r e r 【3 2 ( 1 9 7 9 ) 发展了一些较高精度的p 稳定方法之后，c h a w l a & n e t a 1 6 3 3 5 4 】，a n a n t h a k r i s h n a i a h 2 2 2 3 】，s i m o s 1 2 0 1 1 4 5 - 4 7 1 1 5 3 ，又发展出了一系列p 稳定方法。 其中较为典型的是a n a n t h a k r i s h n a i a h 的两个p 稳定方法【2 3 】和s i m o s 发展的 利用三角拟和技术与高阶微商构建的两步p 稳定o b r e e h k o f f 方法【2 0 】。下面我们 将对这几个方法做简单介绍。 上海大学硕士学位论文 1 4 2 4a n a n t h a k r i s h n a i a h 的p 稳定方法 1 9 8 7 年，a n a n t h a k r i s h n a i a h 2 3 根据前面所述的l a m b e r t 的稳定性理论【2 7 】 提出了一种p 稳定的o b r e c h k o f f 方法来求解二阶微分初值问题，包括两个方法， 我们分别称之为2 1 1 方法和2 1 2 方法。 主结构如下： m 嗍咖协一2 ) 伪= 耖( f l i o ( y ( 2 0 ( x + h ) + 尹纠+ 确尹嘲( 2 4 1 ) 将( 2 4 1 ) 代入测试方程y ”= 一m 2 y ，得到特征多项式 p ( 五) = 彳( h ) 旯2 2 b ( h ) a + a ( h ) = 0 ( 2 4 2 ) 其中， 彳( 日) = l + ( - 0 ”1 风h “， ： ( 2 4 3 ) m 、 b ( h ) = 1 - ( 一1 ) 尾h “ 这里h = 蛔，在本文后面的内容也将这样表示。如果( 2 4 3 ) 式满足， i 堕盟l 1 ， l 爿( h ) l ( 2 4 4 ) 则( 2 4 2 ) 式的根是复数，且模为1 。 设 2 = p + 徊m ，则特征方程为： - 2 c o s 8 ( h ) 月+ 1 = 0 ( 2 4 5 ) 这里 e o s o ( h ) = b ( h ) a ( h ) ( 2 4 6 ) 根据相位滞后的定