（计算数学专业论文）高振荡微分方程的对称数值解法.pdf

中文摘要 摘要：高振荡微分方程是指其解具有高振荡性的一类微分方程，在分子动力学、 天体力学、量子化学以及原子物理等方面有着广泛的应用。因此，研究其数值解 法具有重要意义。然而对于高振荡微分方程，一般的数值解法难以给出好的计算 结果。 设计数值计算格式的一个基本想法是数值解法保持原问题的基本特征。对于 h a m i l t o n 方程，h a m i l t o n 函数( 通常表示能量) 为守恒量。按照上述原则，数值 解法也应较好地保持这一性质。 本文系统地介绍了h a m i l t o n 方程的性质、辛几何算法及对称、组合、分裂数 ，a rr 、 值方法。主要讨论了形如戈+ q 2 x = g ( 工) lg ( x ) = 一半l 的一类高振荡微分方程。此 o x 类方程可以写成h a m i l t o n 方程，相应的h a m i l t o n 函数为守恒量。f p u 问题可以表 示成该类方程的形式。针对这类方程，本文给出了两个新的对称数值格式。以f p u 问题为例进行了数值实验，实验结果显示，这两个解法具有较好的能量保守性。 关键词：高振荡微分方程；h a m i l t o n 方程；辛几何算法；脉冲法；对称数值解法； m o d u l a t e df o u r i e r 展开方法 分类号：0 2 4 1 8 1 a bs t r a c t a b s t r a c t ：h i 龇y - o s c i l l a t o r yd i f f e r e n t i a le q u a t i o n sa l e ak i n do fe q u a t i o n s 、1 1 0 s e s 0 1 u t i o i l sa r eh i g h l y - o s c i l l a t o r y , w h i c ha r ee x t e n s i v e l ya p p l i e d i nm o l e c u l a rd y n a n u c s ， c e l 础a lm e c h a i l i c s ，q u a n t u mc h e m i s t r y , a t o m i cp h y s i c sa n d s oo n t h e r e f o r e ，i ti ss 1 班 。f i c a n tt os n m yi t sn u m e r i c a lm e t h o d s f o rh i g h l y - o s c i l l a t o r y d i f f e r e n t i a le q u a t i o n s ， h o w e v 既i t i sh a r dt og i v eg o o dc o m p u t a t i o n a lr e s u l t sw i t hg e n e r a ln u m e r i c a l m e m o d s ab a s i ci d e ab 抛l dt h ed e s i g no ft h en u m e r i c a ls c h e m e si s t h a tt h e yc a l ll y r e s e e m ep r o p e n i e so ft h eo r i g i n a lp r o b l e m sa sm u c ha sp o s s i b l e h a m i l t o n i a n f u n c t l o n s “柚i c hu s u a l l ym e a ne n e r g y ) a r ec o n s e r v a t i v eq u a n t i t i e s o fh a m i l t o n i a ne q u a t i o n 3 a c e o r d i n gt 0m ea b o v ep r i n c i p l e ，w ea l s oe x p e c tt h a tn u m e r i c a lm e t h o d sc o u l dk e e p t h i sp r o p e r t yb e t t e r i i lt h i sp 印w es y s t e m a t i c a l l yi n t r o d u c et h ep r o p e r t i e so f h a m i l t o n i a ne q u a t i o 粥， 哪n p l e c t i cg e o m e t r i ca l g o r i t h m sa n ds y m m e t r i c ，c o m p o s i t i o na n ds p l l t t i i l gn 啪c n c a j m e t l 】【o d s ，w em a i n l yd i s c u s sak i n do fh i g h l y - o s c i l l a t o r yd i f f e r e n t i a le q u a t i o n sw m c h t a k et l l ef o m 珏出习fg ( 曲；一掣1 t h i sk i n do fe q u a t i o n s c a nb ew n t t e na s h 鼬i l t o n i a nc q u a t i o i l s ，a n dt h ec o r r e s p o n d i n gh a m i l t o n i a nf u n c t i o n s a r ec o n s e n r a t l v e q u a n t i t i e s t h ef p up r o b l e m sc a n b ee x p r e s s e da st h ef o r mo ft h i sk i n do fe q u a t i o n s w e 咖t w on e ws y m m e t r i cn u m e r i c a l s c h e m e sf o rt h ee q u a t i o n 8 t h em 姐唧砒 e x p e m 饥tr e s u l t sf o rf p up r o b l e m ss h o w t h a tt w os y m m e t r i cn u m e n c a lm 甜l o d sh a v e b e t t e rb e h a v i o ro fe n e r g yc o n s e r v a t i o n k e yw o r d s ：h i g h l y o s c i l l a t o r y d i f f e r e n t i a le q u a t i o n s ；h a m i l t o n i a ne q u a t i o n s ； s 卿l c c t i cg e o m e t r i ca l g o r i t h m s ；i m p u l s em e t h o d ；s y m m e t r i c h u m e r i c a lm 劬。d ； m o d u l a t e df o u r i e re x p a n s i o nm e t h o d c l a s s n o ：0 2 4 1 8 1 学位论文版权使用授权书 本学位论文作者完全了解北京交通大学有关保留、使用学位论文的规定。特 授权北京交通大学可以将学位论文的全部或部分内容编入有关数据库进行检索， 并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学校向国 家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 学位论文作者签名：豫铀 签字日期：矽略年多月弓1 7 t 导师签名：起千福 签字日期：矽口缉月弓日 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的研 究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表或 撰写过的研究成果，也不包含为获得北京交通大学或其他教育机构的学位或证书 而使用过的材料。与我同工作的同志对本研究所做的任何贡献均已在论文中作 了明确的说明并表示了谢意。 学位论文作者签名：陈销签字日期：刀口8 年乡月多日 3 5 致谢 本论文工作的顺利完成，首先要感谢我的导师赵平福老师对我的悉心指导和 亲切关怀，也离不开赵老师对我的严格要求和敦敦教诲。无论是在基础课学习阶 段，还是在论文的选题、研究及写作阶段，赵老师始终都给予我许多的支持和帮 助，在此特向赵老师表示深深的敬意和感激。 两年来，赵老师严谨的治学态度和科学的工作方法使我从中学到了很多东西， 这将是我今后学习生活和工作中宝贵的财富，在此再次向赵老师表示由衷的感谢! 感谢研究生期间给予我传道授业解惑和生活学习上关心帮助我的所有老师， 同时感谢所有一路走来、互相勉励的同学和朋友们，正是你们增添了我动力和不 断前进的信念。 感谢我学习生活了6 年的交大母校，在这里我度过了人生最宝贵的一段岁月， 是知行校训让我不断成长。今后，我将用自己的表现来报答母校的养育之恩! 最后感谢我的父母，是你们的鼓励和支持，给了我一个舒适安定的学习环境， 使我能够专心地在学校完成学业，衷心地感谢他们! 1 引言 高振荡微分方程是其解具有高振荡性的一类微分方程，在分子动力系统、天 体力学、量子化学以及原子物理等方面有着广泛的应用。因此，系统地研究其数 值解法具有重要意义。 对于高振荡微分方程而言，一般的数值解法难以得到好的计算效果。近年来， i s e r l e s 研究了利用m a g n u s 展开构造线性高振荡微分方程的数值解法【l 卅； s a n z ，s e r n a 等从计算精度及稳定性方面出发研究了高振荡微分方程组数值解法构 造问题1 6 ；e h a i r e r 等研究了高振荡微分方程组对称的数值解法【_ 7 8 】。 当代计算方法的一条不成文的基本法则是，数值离散应尽可能地保持原问题 的基本特征。因此对于具有一定结构的系统，计算方法也应保持相应结构。h a m i l t o n 系统具有相流保持面积的性质，数值解法也应具有类似的性质。正是从这种想法 出发，冯康教授从1 9 8 4 年开始系统地研究计算h a m i l t o n 系统的辛方法。经过多年 的研究，冯康和他的研究小组获得了许多重要理论结果并从计算实验方面论证了 辛格式具有独特的计算稳健性与长时间跟踪能力【9 1 3 1 。类似地，基于系统结构保守 性的数值解法还包括多辛方法【1 4 1 、李群方法【4 , 7 , 1 5 1 等等。对于这些保持原系统结构 的数值解法，我们可以统称为保结构算法。 厂, - q r ，、 在本文中，我们主要研究形如戈+ q 2 x = g ( x ) lg ( x ) = 一半l 的一类高振荡微分 麟 方程组。这类方程组可以写成h a m i l t o n 方程组形式。除了相流保持面积，h a m i l t o n 方程组的另一个重要性质是h a m i l t o n 函数为守恒量。考虑数值解法时，我们也希 望数值解法具有较好的能量保守性质。在文【8 】中，e h a i r e r 等从保持能量守恒性入 手研究了一些具有对称性质的数值解法。在e h a i r e r 等研究工作的基础上，我们构 造了两个新的对称数值解法，并以f p u 问题为例进行了数值实验，实验结果显示， 这两个格式具有较好的能量保守性。在文【8 】和【1 6 】中，e h a i r e r 等利用m o d u l a t e d f o u r i e r 展开研究数值解法的能量保守性。在本文我们对m o d u l a t e df o u r i e r 展开方 法作一简单介绍。 在论文的第二部分，我们首先介绍h a m i l t o n 方程组，并给出其具有的一些性 ，, - - q lr 、 质。然后给出线性高振荡常微分方程和形如孑+ q 2 石= g ( z ) ig ( 石) = 一半l 的一类高 缳 振荡微分方程组。最后介绍了f p u 问题，并给出其h a m i l t o n 形式。 第三部分，我们先介绍对称、组合、分裂三种方法的基本概念和基础知识， 之后介绍了h a m i l t o n 系统辛几何算法及其相关性质，最后介绍脉冲法的基本知识。 第四部分，以刀= 3 的f p u 问题介绍脉冲法并分析其缺陷，之后给出两个对称 的数值解法。同时，给出这两个方法与脉冲法及其它一些方法的数值比较结果。 最后一部分，我们对m o d u l a t e df o u r i e r 展开方法作简单介绍。 2 2h a miit o n 系统与高振荡微分方程 2 1h a mi lt o n 系统 经典力学问题有n e w t o n 力学、l a g r a n g e 力学及h a m i l t o n 力学三种表示形式。 这些不同的数学形式具有等价的关系，可以陈述同一物理规律，但是由于形式上 的差异，在实践中是不等效的。 首先，我们来看经典力学的这三种形式。考虑有疗个自由度的运动，位置向量 q = ( 9 1 ，一，q 。) r ，势能函数矿= 矿( g ) ，则n e w t o n 形式有 m 象：一毒y ，m ：fm 它是以维位形空i 刚r ”中的2 阶微分方程组，是经典力学的标准形式。 e u l e r 及l a g r a n g e 通过引进动能与势能函数的作用量 l ( q ，口) = r ( 口) 一矿( g ) = 寺( 圣，m ? i ) - v ( q ) ， 利用变分原理将运动方程写成如下等价形式 d o l 0 l ：o d ta q a q 上式被称为经典力学的变分形式，即l a g r a n g e 形式。 h a m i l t o n 利用动量p 与总能量h 将运动方程描述为 户：一i o h ，圣：i o h ， ( 2 1 1 )p2 一i 一，g 2 i 一， u l l j a g印 其中，日= 妻( 口，岣) + y ( g ) ，变量p = ( a ，p 。) r ，g = ( g l 一，吼) r 。方程组( 2 1 1 ) 称为h 锄i l t o nj e , 贝。方程组，日为( 2 ) 的h 锄；l t o n 函数。取z = ； ，有 警母阱1 以 i。j 那么h a m i l t o n 方程组( 2 1 1 ) 可写成 生=lhdt 0玎1 警， )一= i j i i = ，一 z 1 z - 昆 、 这里为以阶单位矩阵= 匕台 为反对称矩阵，厂1 r “。由方程 组( 2 1 2 ) ，可以看出h a m i l t o n 正则方程组具有简单而对称的形式。 h a m i l t o n 形式具有一些优点。首先，h a m i l t o n 形式可以最明显地体现运动的 规律性。其次，h a m i l t o n 形式具有远比n e w t o n 形式更广的遍在性和普适性，它覆 盖了经典的、相对论的、量子性的、有限或无限自由度的一切真实的、耗散效应 可忽略的物理过程。因此，研究h a m i l t o n 系统具有极为重要的现实意义和广阔的 发展前景。 h a m i l t o n 系统存在于一些我们所熟悉的经典物理模型中，同时也广泛应用于 结构生物学、药理学、半导体、等离子体、材料和偏微分方程等问题。例如，物 理学中经典的单摆模型具有下面形式的h a m i l t o n 函数： 日( p ，q ) = 去p 2 - c o s q ， 其运动方程可以表示成h a m i l t o n 正则方程( 2 1 1 ) 的形式： p = - s i n q ，矗= p 又例如，天体物理中双体运动模型具有以下运动方程形式： 蕴= 一i i 裔虢= 一i i 蕾。 该方程等价于具有以下h a m i l t o n 函数形式的h a m i l t o n 系统： 日( a ，岛，g l ，吼) = 三( 彳+ 雳) 一丽1 ，只2 幺， 其中p = ( p l ，岛) r ，g = ( g 。，q 2 ) r 。 h a m i l t o n 系统的数学框架是辛几何，h a m i l t o n 为研究牛顿力学，引进了广义 坐标和广义动量来表示系统的能量，即现在通称的h a m i l t o n 函数。对于自由度为n 的系统，刀个广义坐标和厅个广义动量可张成2 n 维相空间。于是，牛顿力学转化 成了辛几何学。 下面我们将辛几何与欧氏几何相比较，以说明辛几何的性质。r ”的欧几里德 结构取决于双线性对称的、非退化的内积： ( x , y ) = ( x ，i y ) ，= l 4 由于非退化，当x o 时( 毛x ) 恒正，从而可以确定长度= ( x ，石) o 。保持 内积( 即长度不变) ，即满足a r a = i 的线性算子a 组成正交群d ( 刀) 。这是一个典 型的李群，其李代数d ( 甩) 由反对称变换组成，即由无穷小正交变换所组成。 辛几何则是相空间中的几何学，辛空间是指相空间具有特定的辛结构，它取 决于满足双线性、反对称、非退化性质的辛内积： i x , y 如础m 。= 匕讣 当n = 1 时，有 k 小臣卦 这就是以向量x ，y 为边的平行四边形面积( 严格讲为有向面积) 。一般地，辛内 积是面积度量，由于其具有反对称性，对于任意向量石恒有f x ，石1 - 0 ，因此辛结构 中没有长度的概念。保持辛内积不变的线性变换满足彳r j a = j ，它们构成一个群 s p ( 2 n ) ，称为辛群，其李代数由无穷小辛变换b ( 满足曰7 - ，+ 捆= 0 ) 组成，并以 s p ( 2 n 1 表示。欧氏几何与辛几何根本的差别在于，欧氏几何是研究长度的几何学， 而辛几何是研究面积的几何学。 h a m i l t o n 系统具有相流保持面积和h a m i l t o n 函数日为守恒量的性质，下面对 其给出证明，我们首先证明相流保持面积。 对于相空间上的任意点z 0 ，方程组( 2 1 2 ) 的解z ( f ) 满足初值条件z ( o ) = z o ，存 映射识，使得 仍( z o ) = z 【f j ， 称映射仍为h a m i l t o n 相流。 给定初值( p ，g ) 1 ，方程组( 2 1 2 ) 的相流仍的导数具有如下形式： 慨) = 笺铲 h a m i l t o n 系统具有相流保持面积的性质，即有 ( ( p ，g ) ) 7 ，( 科( p ，口) ) = ， ( 2 1 3 ) 我们对此性质作一简单证明。当t = 0 时，上式显然成立，我们只需证明 烈d 仍i ( p ，g ) ) r ，( 科( 刖) ) ) = o ( 2 1 4 ) 注恿剑 丢帆) _ 薏薏卜g ) ， 易证( 2 1 4 ) 成立。 下面利用p o i s s o n 括号的相关知识证明h a m i l t o n 函数h 的保守性。 定义2 1 1设f 表示r 2 “上的光滑函数的集合。一个p o i s s o n 括号 ， 是一个 算子f x f 专f ，这个算子满足性质： ( i ) ， 是双线性反对称的， ( i i ) ，) 满足j a c o b i 等式： 尸，q ，尺 + r ，砖，q + q ，r ) ，p ) = o ， ( i i i ) ，) 满足l e i b n i t z 法则 p q ，r ) = p q ，r ) + 尸，r o ， 其中p ，q ，r f 。 r 2 ”上的两个光滑函数缈( z ) 和y ( z ) 的p o i s s o n 括号 伊( z ) ，y ( z ) 可定义为： 伊，杪) ( z ) = ( 鲠) 2 厂1 以 匕j 捌 叫， = 一( 形一硝) 可以验证( 2 1 5 ) 所定义的f ，一f 映射满足定义2 1 1 的要求9 1 。 定义2 1 2h a m i l t o n 相流的首次积分是这样的函数p ：r 2 4 专r ，p 满足 p ( z ( f ) ) = c ，v f r ， 其中c 是常数。 定理2 1 3p 是h 锄i l t o n 方程组( 2 1 2 ) 的一个首次积分当且仅当 只日) = o 。 证明：对p ( 仍( z ) ) 求导可有： 丢p ( 仍( z ) ) = 尸，日 ， 因此由首次积分的定义，p 是首次积分当且仅当 p ，h ) = o 。 定理2 1 4h 是h a m i l t o n 方程组( 2 1 2 ) 的首次积分。 证明：由p o i s s o n 括号定义， h ，h ) = 0 。因此，由定理2 1 3 ，h 是h a m i l t o n 方 程组( 2 1 2 ) 的首次积分。 6 在这一节中，我们简要地介绍了h a m i l t o n 方程的一些性质。在【1 7 】，【1 8 】中， 有关于这一部分内容的详细讨论。 2 2 高振荡微分方程 高振荡微分方程是指其解具有高振荡性的一类微分方程，它在各个领域具有 广泛的应用。i s e r l e s 在文【1 】中研究了线性高振荡微分方程 y 。+ g ( f ) y = o ，t 0 ，y ( o ) = y o ，y ( o ) = “ ( 2 2 1 ) 其中g ( f ) 0 ，f 0 ，熙g ( f ) = 佃，o f ( 1 z l ，g ( x ) = c ( o ) ，c 为常数。设初值而= o ，磊= 1 ，取乃彩：2 刀，有 量： = 妻 + 苫 = 1 + 以j 1 j l l c ，刀= 。，2 ，c 4 t 4 ， 而方程( 4 1 3 ) 的准确解为 在整数点o - 满足戈( 厶+ t ) = 南= 1 ，石以+ ) = 而= o 。因此，当刀较大时，毫+ 。与j 纯+ 。) 之间误差很大。 下面我们以f p u 问题为例，利用脉冲法进行数值实验，通过能量误差分析， 说明脉冲法在计算中存在的缺陷。对于栉= 3 的f p u 问题，选取初值“( o ) = ( 1 ，0 ，o ) ， 缈 甜 彻 埘 咖 一 一，一国!缈 “ 斗 卜 澎 研 泐 鼽 垮 薹 啪 五( o ) = ( 1 ，0 ，0 ) ，v ( o ) = ( 彩，o ，0 ) ，t ( o ) = ( 1 ，0 ，0 ) 。在不同的j i l 缈取值条件下计算脉 冲法的能量误差，则有以下一些实验结果： ( a ) h = 0 0 1 ，c o = 1 0 0 0 0 0 0 1 万h ( c ) h = o 0 1 ，c o = 2 0 0 0 0 0 0 1 万h( d ) h - - - - 0 ( 1 2 5 ，c o = 2 0 0 0 0 0 0 1 刀lh ( c ) h = o 0 1 ，c o = 2 5 ，r h( 0 h = 0 ( 1 2 5 ，c o = 2 5 n h 图4 1 脉冲法的能量误差图 图4 1 是h 分别取o 0 1 、0 0 2 5 时脉冲法的能量误差。由图4 1 可以看出，当h c o 接近7 和2 n 时，脉冲法的能量误差很大。 为了进一步说明脉冲法的缺陷，我们在取相同步长、不同h t o 的条件下，做出 脉冲法的最大能量误差图。保持原先f p u 实验的初值条件不变，取步长h = o 1 ， 2 l 实验结果如下： 图4 2 脉冲法的最大能量误差图 图4 2 是h = 0 1 时脉冲法的最大能量误差图，其横坐标表示h c o 的取值，纵坐 标表示最大能量误差。从图4 2 可以看出，脉冲法在远离万的整数倍时有较小的能 量误差，但在接近万的整数倍时误差较大。因此，我们需要研究能量保守性更好的 改进方法。 4 2 一些修正的脉冲方法 基于对脉冲法的修正，e h a i r e r 等研究了如下的数值解法： 酬= 篇固c o s s 警h 甜i i qti 尝( 甲。邑+ 甲。邑+ 。) 嬖 2 。4 ， ( 4 2 1 ) 这里岛= g ( 魄) ，= ( 删，、壬，= 以，v o = ( 刎，一= 惭( 目，其中g 、 、奶为光滑函数。选取不同的，甲，甲。，、壬，。，相应地可构造出不同的数值解 法。例如选取【8 】 矽( 0 = s i i l c ( 0 ，吵= 血c ( 回，= s ( 旬，= 1 ， ( 4 2 2 ) 其中s i n = s i n 孝善。利用方法( 4 2 2 ) ，对，l = 3 的f p u 问题进行数值实验，有以下 数值结果： ( a ) h = o 0 1 ，0 9 = 1 i x ) o o ( o l ：r h ( c ) h = o 0 1 ，缈= 2 0 0 0 0 0 0 1 n h ( e ) h = 0 0 1 ，c o = 2 5 7 r h ( d ) h - - q ( 1 2 5 ，0 9 = 2 0 0 0 0 0 0 1 r c h ( 0 h = 0 ( 1 2 5 ，e o = 2 5 万h 图4 3 方法( 4 2 2 ) 的能量误差图 类似于图4