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

中文摘要 摘要：高振荡微分方程是其解具有高振荡性的一类微分方程，它广泛应用于诸如 分子动力学、天体力学、量子化学以及原子物理等方面。因此，研究其数值解法 具有重要意义。 对于高振荡微分方程，一般的数值解法难以给出好的计算结果。例如，对于 形如y 。+ g ( f ) j ，= 0 的线性高振荡方程，经典的方法在处理该类问题时均会产生较大 的误差。近来，i s e r l e s 利用m a g n u s 展开方法详细研究了该类方程数值解法问题， 给出了计算结果较好的数值算法。e h a i r e r 等研究了高振荡微分方程对称的数值解 法。 本文介绍了h a m i l t o n 方程的性质、辛几何算法、对称方法、m a g n u s 展开和 n e u m a n n 展开方法。在m a g n u s 展开和n e u m a n n 展开方法中，迭代都起着重要作 用，我们考虑了利用迭代法构造高振荡微分方程的数值解法。对于基于迭代的数 值解法，以f p u 问题为例进行了数值实验，实验显示，该方法可给出较好的数值 结果。 关键词：高振荡微分方程；h a m i l t o n 方程；辛几何算法；m a g n u s 展开方法；n e u m a n n 展开方法；迭代法 分类号：0 2 4 1 8 1 a bs t r a c t a b s t r a c t ：h 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 s a reak i n do fe q u a t i o n sw h o s e s o l u t i o n sa r eh i g h l y o s c i l l a t o r y , i ti se x t e n s i v e l ya p p l i e di na s p e c t ss u c ha sm o l e c u l a r d 姗i c s ，c e l e s t i a lm e c h a n 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 i g n i f i c a n tt os t u d 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 yd i f f e r e n t i a le q u a t i o n s ，i ti sh a r dt og i v eg o o dc o m p u t a t i o n a l r e s u l t sw i t hg e n e r a ln u m e r i c a lm e t h o d s f o re x a m p l e ，d e a l l i n g w i t ht h el i n e a r h i g h l y - o s c i l l a t o r ys y s t e m sy 。+ g ( t ) y = 0 ，c l a s s i c a l n u m e r i c a lm e t h o d sw i l lp r o d u c e b i g g e re 仃0 r r e c e n t l y , u s i n gm a g n u se x p a n s i o ni s e r l e sh a ss t u d i e dn u m e r i c a lm e t h o d s f o rt h i sk i n do fe q u a t i o n si nd e t a i la n dg i v e ng o o dn u m e r i c a lm e t h o d s e h a i r e r e ta l h a v es t u d i e ds y m m e t r i cn u m e r i c a lm e t h o d sf o rh 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 s i nt h i sp a p e r , w ei n t r o d u c et h ep r o p e r t i e so fh a m i l t o n i a ne q u a t i o n s ，s y m p l e c t l c g e o m e t r i ca l g o r i t h m s a n ds y m m e t r i c ，m a g n u se x p a n s i o na n dn e u m a n ne x p 锄s l o n m e t h o d s i t e r a t i o n p l a y s a ni m p o r t a n tr o l e i nm a g n u se x p a n s i o na n dn e u m a l l n e x p a n s i o nm e t h o d s ，u s i n gi t e r a t i v e m e t h o dw ec o n s i n d e rn u m e r i c a lm e t h o d s 士o r l l i 龇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 t h en u m e r i c a le x p e r i m e n tr e s u l t sf o rf p u p r o b l 锄ss h o wt h a tt h em e t h o dw h i c hb a s e d o ni t e r a t i o nc a ng i v eb e t t e rn u m e n c a l r e s u l t s k e y - w o r d s ：h i g h l y o s c i l l a t o r y s y m p l e c t i cg e o m e t r i ca l g o r i t h m s ；m a g n u s m e t h o d ；i t e r a t i v em e t h o d c l a s s n 0 ：0 2 4 1 8 1 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得的研 究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经发表或 撰写过的研究成果，也不包含为获得北京交通大学或其他教育机构的学位或证书 而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中作 了明确的说明并表示了谢意。 学位论文作者签名：终 氛签字日期加7 年7 月幺日 3 3 学位论文版权使用授权书 本学位论文作者完全了解北京交通大学有关保留、使用学位论文的规定。特 授权北京交通大学可以将学位论文的全部或部分内容编入有关数据库进行检索， 提供阅览服务，并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。 同意学校向国家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 廖构撬 7a 分日 导师签名： 赵哥福 签字日期：训7 年7 月2 日 致谢 本论文的工作是在我的导师赵平福老师的悉心指导下完成的，从课题的选择 到最终完成，赵老师都始终给予我指导和支持，在此特向赵老师表示深深的敬意 和感激! 两年来，赵老师严谨的治学态度和科学的工作方法给了我极大的帮助和影响， 从赵老师身上，我不仅学到了许多专业知识，而且懂得了做人的道理。通过赵老 师细心地讲解，使我对知识有了更深入的理解，为今后学习和工作打下坚实的基 础，在此再次向赵老师表示由衷的感谢! 真诚地感谢北京交通大学理学院的老师和同学们，两年的时光里，你们教会 我太多太多。在我遇到困难的时候指导我，在我感到挫折的时候鼓励我，正是你 们无私的帮助，让我体会到了师生情谊、同窗情谊的美好! 也正是由于你们的帮 助和支持，我才能克服一个一个的困难和疑惑，直至本文的顺利完成。 感谢我的家人，给了我一个舒适安定的学习环境，正是你们的无私奉献和默 默支持，让我能够专心学习，顺利完成我的学业。 最后感谢北京交通大学理学院给我这样一个学习深造的机会! 1引言 高振荡微分方程是其解具有高振荡性的一类微分方程，它广泛应用于诸如分 子动力学、天体力学、量子化学以及原子物理等方面。因此，系统地研究其数值 解法具有重要的现实意义。 高振荡微分方程的数值求解历来是微分方程的数值求解问题中，最富有挑战 性难点问题之一。经典的数值算法效果甚微，这也激励人们探讨更富有挑战性的 数值方法【l 。2 1 。对于线性高振荡微分方程的数值求解问题，近年来，i s e r l e s 利用 m a g n u s 展开方法对其进行了详尽的研究【3 。7 】，并给出了好的数值解法；对于非线性 高振荡微分方程的数值求解问题，e h a i r e r 等研究了对称的数值解法【2 ，8 】。 m a g n u s 展开式已经有效地用于数值求解线性矩阵微分方程，甚至用于求解李 群上的非线性微分方程，h a m i l t o n 系统的微分方程。用m a g n u s 展开式求解上述微 分方程，其优点在于适当的截断展开式仍能保持方程精确解的部分几何性质【9 】， i s e r l e s 等在【6 】中用m a g n u s 展开式求解线性高振荡微分方程，得到了比传统数值方 法( 如r u n g e k u t t a 法、线性多步法) 更好的结果，使计算更加精确、有效。在文 1 0 中，i s e r l e s 介绍了n e u m a n n 展开方法。 h a m i l t o n 体系是动力系统的重要体系，一切真实的，耗散可忽略不计的物理 过程都可表示为h a m i l t o n 系统，这类系统的本质特征是相流保持相空间的辛结构 不变( 即保持面积) ，针对这一特性导致了辛算法( 保持辛结构的数值方法) 的产 生，冯康教授于1 9 8 4 年开始系统的研究计算h a m i l t o n 系统的辛算法，深入的理论 分析和大量的数值实验令人信服的表明辛算法在稳定性和长期跟踪能力上具有独 特的优越性【l l - 1 5 】。 我们主要研究形如 m 2 一( 咖) = 一罢) 的一类高振荡微分方程组，可将这类方程组改写成h a m i l t o n 方程组的形式。由于 h a m i l t o n 函数( 能量) 为守恒量，故希望数值解法也具有较好的能量保守性。在 m a g n u s 和n e u m a n n 展开中，迭代起着重要作用，我们利用迭代法对上述非线性高 振荡微分方程进行求解，进行适当截断后得到一个数值解法，基于此解法给出一 个修正方法，并以f p u 问题为例进行了数值实验，实验结果显示修正的数值解法 具有较好的能量保守性。 论文第二部分，我们分别介绍了线性高振荡微分方程组和形如 2 一( 工) k ) = 一豢) 的一类高振荡微分方程组，并各自给出了一种求解上述微分方程组的传统方法， 分析了数值解法存在的缺陷。 第三部分，我们先介绍了h a m i l t o n 系统辛几何算法，之后介绍了对称方法的 基本概念，最后介绍了三个对称的数值解法，同时给出了这三个数值解法与脉冲 法的数值比较结果。 第四部分，介绍了m a g n u s 、修正的m a g n u s 和n e u m a n n 展开方法，并分别给 出用修正的m a g n u s 展开方法和n e u m a n n 展开方法构造线性高振荡微分方程的数 值解法。 第五部分，我们首先介绍了迭代法求解非线性高振荡微分方程的数值解法， 之后给出了修正的数值解法，也给出了关于f p u 问题的数值结果。 最后一部分，通过对文中数值实验结果的分析，我们给出最后的结论。 2 扯毫奎丑盘堂亟堂僮监塞 2直拯蓬邀盐直捏 2 高振荡微分方程 2 1 线性高振荡微分方程 高振荡微分方程是其解具有高振荡性的一类微分方程，它在各个领域具有广 泛的应用。 我们称方程 y 4 + 妒= o ，y ( o ) = 托，( o ) = “ ( 2 1 1 ) 为a i r y 方程，若取a i r y 方程的初始倩为y ( o ) = 1 ，( o ) = 0 ，则它的精确解为 y ( t ) = ，r a i ( - t ) b i ( o ) 一a i ( o ) 目( 一f ) 其中a i ( t ) ，b i ( t ) 分别是第一类和第二类a i r y 函数。 图2 l ：a 打方程的精确解 由图2 1 可见，a i r y 方程是线性高振荡微分方程。较一般的线性高振荡微分方 程可写作 ，+ g ( f ) ，= o ，t - 0 ，y ( o ) = y o ，( o ) = “， ( 2 12 ) g e p g ( f ) ) 0 ，f ，0 ，熙g ( f ) = 一，o - = ( q ，也，吃) r ， 1 ，= ( m ，v 2 ，v 3 ) r ， = 匕 “= ( ”l ，1 2 ，吩) r ， 计 ( 2 2 6 ) ( 2 2 7 ) ( 2 2 8 ) 此时刀= 3 的f p u 问题可表示成方程( 2 2 2 ) 的形式，光滑函数g ( “，1 ，) 为： g ( u ，) = ( g 。，g ：，g ，g 。，g ，g 。) r ， 6 蜀= 一( u t - v , ) 3 + ( 一屹一“一m ) 3 ，g ：= 一u - - v 三一“。一m ) 3 + ( u 3 - v 3 - “：一心) 3 ， 岛：一u + b ) 3 一( 吩一屹一“：一v j ) 3 ，9 4 = ( q h ) 3 + ( “：一吃一l l i m ) 3 ， ：u ：一1 ，2 - - u ! - - m ) 3 + ( u 3 - - v j 一“：一v ! ) 3 ，9 6 = 一( “，+ b ) 3 + ( “，一b 一“：一吃) 3 方程组( 2 2 2 ) 的准确解形式如下： ) _ ( 盎一f 一s i n t q 舳f f ：c o + fh t , n - s m i no 吖- s ) n 卜yt 小肛，( 2 2 聊 其中x o ，矗是方程组的初值。 我们对h a m i l t o nf f 撇h ( z ( t ) ) 求导，可得 掣：( 阳) r - ，一阳 由 ，= e o l 讣 计算可得 广= 巴讨 厂1 为一反对称矩阵，由反对称矩阵性质，可以得到 必：( v 日) r 厂1 跗：0 d t 、7 因此，日( z ( f ”= c ( c 为常数) ，所以h a m i l t o n 函数h 为守恒量。 由于h a m i l t o n 方程组中h 为守恒量，所以考虑到数值解法我们也希望数值解 法能较好的保持h 的守恒 1 5 。 对( 2 2 9 ) 右端的积分使用梯形公式可导出脉冲法，用脉冲法对方程( 2 2 2 ) 进行 求解得： _ ( = 一嚣 其中岛= g ( 吒) ，s i n c 孝= s i n e ( 7 ( 2 2 1 0 ) 一、i， 岛 岛 _ b 0 是 2 c 西 m 刚 两 矿一2 磊一2 ，。_ 卜 韭夏至遭五生班芏位监塞2 商握蓝丝盐直捏 初值取为“( o ) = ( 1 ，0 , 0 ) ，( o ) = ( 1 ，0 ，0 ) v ( o ) = ( m ，o , 0 ) ，i ( o ) = ( 1 ，0 ，0 ) ，利 用格式( 2 2 1o ) 对月= 3 的f p u 问题进行数值实验在不同的 甜取值条件下计算 方法( 22l0 1 的能量误差，有以下实验结果： _ i m _ m i f a ) h = 0 0 1 ，o = 1 0 0 0 0 0 0 1 z _ m _ _ ” o ) = 伽西= 10 0 0 0 0 0 1 z ：阿t 1 ：r 一 一 _ i c ) h ：0 0 1 ，= 2 0 0 0 0 0 0 1 # h( d ) = q 晒，m = 2 0 0 0 0 0 0 1 z h 删懈懒丌 p 1 讯即哪即 ( e ) h = c t 0 1 to = 2 新h ( n = 0 ( 1 2 5 ，甜= 2 5 口h 图2 3 冉法f 22l o ) 的能量误差图 图2 3 是h 分别取oo l ，0 0 2 5 时方法f 22 10 ) 的能量误差图。由图2 3 可以看 出，当妇，接近f 和2 口时，方法( 22 1 0 ) 的能量误差很大，因此，我们需要研究能 量保守性更好的改进方法。 另外，文献1 1 7 中给出了关于高振荡h a m i l t o n 偏微分方程的一些研究结果。 3h a m i l t o n 方程及辛算法、对称方法 3 1h a m i l t o n 方程及辛几何算法 经典力学有n e w t o n 力学、l a g r a n g e 力学、h a m i l t o n 力学- - 矛t 表示形式，这些 不同的数学形式陈述同一物理规律，但它们在实践中对解决问题会提供不同的途 径，因此，在实践中是不等效的。下面我们以单摆运动( 假设质点质量为1 ，绳长 为1 ，重力加速度为1 ) 为例来说明经典力学的这三种形式。 按照牛顿第二定律，质点在外力，的作用下所获得的加速度与所受的力f 有 下列关系： f = 砌， ( 3 1 1 ) 其中m 为质量，a 为质点某一时刻的瞬时加速度 按此定律，单摆运动方程为 牙= - s i n q ， ( 3 1 2 ) 其中q 为绳偏离竖直位置角度。此方程为单摆运动的牛顿方程形式。 引入作用量 l ( q ，口) = t ( o ) - v ( q ) ， 其中丁( g ) = 去口2 为动能，y ( g ) = 1 一c o s q 为势能，利用变分原理，单摆运动方程写 成l a g r a n g e 形式为： j di o l 一罢：m i n q ：0 - i - 一= ，7 r1= 出8 q l 由于动量p = 香，在h a m i l t o n 系统中，单摆运动具有下面形式的h a m i l t o n 函 数 h ( p , q ) ：要p z + ( 1 - c o s q ) ， 其运动方程可以表示为 o ho h p2 一- = ；- 2 - s i n q ，q2 _ 2p - o ao p 对于有以个自由度的运动，h a m i l t o n 正则方程为 9 p ：一i o h ，香：i o h ， ( 3 1 3 )p2 一一，g2i 一，【j l j j 其中变量p = ( p l ，一，以) r ，q = ( 9 1 一，吼) r ，h 为h a m i l t o n 函数。 取z - 阶有 警母阱 则h a m i l t o n 方程( 3 1 3 ) 可写成 妄= 0 俐j l - 。j 玎1 警， 即为方程( 2 2 1 ) 的形式 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 形式有远比牛顿形 式更广的遍在性和普适性，它覆盖了经典的、量子性的、有限或无限自由度的一 切真实的、耗散效应可忽略的物理过程。 h a m i l t o n 系统具有相流保持面积和h a m i l t o n 函数h 为守恒量的两个重要性 质前面我们已经证明了h 为守恒量，下面我们来看相流保持面积的性质。 对于相空间上的任意点z o ，方程组( 2 2 1 ) 的解z ( f ) 满足初值条件z ( o ) = z 0 ，存 在映射g 刍，使得 磊( 气) = z ( f ) ， 称映射g 二为h a m i l t o n 相流。 给定初值( p ，g ) 7 ，方程组( 2 2 1 ) 的相流g 二的导数具有如下形式： “沪笺铲 h a m i l t o n 系统具有相流保持面积的性质，即有 ( 舀( 舢) ) 2 ( 西( 删) ) = ， ( 3 1 4 ) 我们对此性质作一简单证明。当t = 0 时，上式显然成立，我们只需证明 丢( ( “( 圳r ，( 州舢) ) ) - o ( 3 ) 1 0 注意到 翱舢，= 嚣薏b 劫， 易证( 3 1 5 ) 成立。 在尺2 “空间中，一个映射的切映射就是该映射的j a c o b i 矩阵。令p ( 乙) ，q ( t 。) 和 p ( t o + 。) ，g ( 乙+ 。) 分别表示乙和o 。时p ，q 的值，由( 3 1 4 ) 可有 ( 渊h 渊卜 对于2 玎阶矩阵么，若满足崩= 则称彳为辛矩阵。因此，矩阵渊 为辛矩阵。当用数值方法求解h a m i l t o n 方程组( 3 1 3 ) 时，称使矩阵剿为辛 矩阵的数值方法为辛几何算法( 简称辛算法) 。 显式和隐式的差分格式都可看成是从这一时刻到下一时刻的映射。如果这个 映射是辛的( 正则的) ，则称该差分格式是辛几何格式( 简称辛格式) 。 下面，我们介绍辛的r u n g e - k u t t a 方法，首先简单介绍r u n g e k u r a 方法( 详见 【18 一 2 0 1 ) 。对于方程组夕= f ( y ) ，r u n g e k u a a 方法的一般形式是 以+ l = 咒+ 伊( 吒，儿，h ) ， ( 3 1 6 ) 其中够有如下形式： 伊( ，儿，h ) = 岛k ， ( 3 1 7 ) k i = ( 毛，y 。) ，k = f l 吒+ q h ，+ 五嘞巧l ，i = 2 ，s ， 其中q ，岛，均为常数，整数s 1 ，称( 3 1 6 ) 和( 3 1 7 ) 为s 一级显式r u n g e k u t t a 方法 ( 简称r - k 方法) 。选取适当的c f ，岛，r - k 方法可给出解的高阶逼近。 h a m i l t o n 方程组( 2 2 1 ) 可以写成 d d z t j _ 1 h z = 八南 此时，一个一般的单步s 级r k 方法为： z “1 = z 七+ 6 ：( e ) ， i = l 鬈= z + j l l 口 ( ) ，l i i ( 1 i s ) 时，a ，= 0 ，且在对角线上存在一些 a 。0 ( 1 i s 1 ，则称格式( 3 1 8 ) 为对角隐式r k 格式；不满足上述条件的格式 ( 3 1 8 ) n 称为隐式r k 格式。 b u t c h e r 向量法是一种表达格式( 3 1 8 ) 的简便方法，其形式如下： c么 6 丁 c la 1 1 。 a l s c 2a 2 1 。 a 2 z o c 5a 5 1 气5 龟 6 ； ( 3 1 9 ) 其中c f = 嘞( f = 1 ，2 ，s ) 。这样，一个r - k 差分格式完全由( 3 1 9 ) 中的数来确定， 这种表示称为b u t c h e r 格式 1 1 , 1 9 。 定义3 2 2 将单步的差分格式看成是a k t , 至j jt k + ，的变换，称为格式推进映射。 如果格式( 3 1 8 ) 的推进映射是辛的，即j a c o b i 矩阵警是辛的，则称此格式是辛 r k 方法。 r k 方法的辛条件由s a n z s e r n a 给出，此处我们不加证明给出以下定理( 详 细证明可参看 11 】) 定理3 2 3 设矩阵曰= d i a g ( b l ，吃，吃) ，m = b a + a r b - b b 7 ，如果m = o ， 则格式( 3 1 8 ) 是辛的。 3 2 对称方法 对于自治的微分方程 夕= 厂( y ) ，y ( t o ) = y o ， ( 3 2 1 ) 其相流仍满足等式伊j = 纪成立。在进行数值计算时，我们也希望数值解法能够保 持这种性质。具有类似性质的数值解法称为对称数值方法。下面，我们对该方法 1 2 作一简单介绍。 定义3 1 1 设方程( 3 2 1 ) 的一个数值方法为 z 斛1 = 吼( z ”) ， ( 3 2 2 ) 令：= 二：，称数值方法z n + i = ：( z ”) 为式( 3 2 2 ) 的伴随方法。此时若满足 ：= 哦，则称方法( 3 2 2 ) 是对称的。 对于单步方法。和甲。，伴随方法满足： ( ：) + = 。，( 。甲。) = 甲：。： 下面由上述对称方法的定义来考查e u l e r 方法、梯形方法的对称性。 e u l e r 方法表达式为 此+ l = y n + h f ( y ) ( 3 2 3 ) 将( 3 2 3 ) 中的h 改为一h ，且将 + 1 与，z 互换，计算可得 虬+ l = 以+ 7 矿( y 州) ， 上式与( 3 2 3 ) 不同，故e u l e r 方法不是对称方法 梯形方法形如 1 虼+ l 一只= i a 办( 厂( 儿) + 厂( 虼+ 1 ) ) ( 3 2 4 ) 将( 3 2 4 ) 巾的h 改为一h 得 1 以+ l 一只= 一i a ( 厂( y 。) + ( 以+ l ”， 将上式中n + l 与刀互换可得方法( 3 2 4 ) ，因此格式( 3 2 4 ) 是对称的计算方法 3 3 辛算法及对称方法在高振荡微分方程中的应用 对于方程( 2 2 2 ) ，h a i r e r 和l u n i e h 研究了数值解法 - ( 盎趣e o 豇s h 警心+ qj izj 兰( 甲。邑+ 甲。g 川) 姿 2 。“ ( 3 3 1 ) 其中晶= g ( 呶) ，中= 痧( j i l q ) ，甲= y ( j i l q ) ，甲o = ( q ) ，甲。= 奶( j j l q ) ，并且 有沙( 孝) = s i i l c f ( 孝) ，( 孝) = e o s 孝( 孝) ，其中s i n c 孝= s i n 善孝，痧，为偶函 数。 1 3 韭基至通厶芏鳕尘芏世垃卫丛! l ! ! 应攫且圣显选、盟整应然 数值解法( 33i ) 当满足 y ( f ) = s i n c ( f ) n ( 掌) ，( f ) = c o s 4 ( f )( 332 ) 条件时，该方法为对称的数值解法下面我们作一简单的证明，将( 331 ) 中的h 改 为一h ，有 = 一茹兰。絮等) + 然后将( 333 ) r 卜的n + 1 与h 互换，得 即 f _ 1 ：f c o s h i f 2 l 矗ji n 。s i n h s q m 汕0 1 一- 1 + c o s h n 八j 一罢( 甲。岛+ 甲岛+ 。) 嬖， 昙( 甲。岛+ + 甲岛) 篓。 = ( # s 幽i n q h 删c o s h q 一( 篇s i n n h n 地c o s 汕h f ，1l + l ji 口1八jl 盯。j ( 333 ) 一：( 甲。+ ，+ 甲禹) 篓崛。 将( 332 ) 代入上式计算可得方法( 331 ) ，即证明格式( 331 ) 当满足( 332 ) 时是对称的 数值解法 对( 33 1 ) 选取不同的中，、i ，甲o ，甲，相应的可给出不同的数值解法t 例如 选取【8 l 烈0 = ”c ( 毋“封= m c ( 旬，( 0 = 叫日，( 日= 1 烈0 = 1 ，州旬= m ，( 曰( 自= 盘1 c ( 日雠( 自，( 旬= 一c ( 旬 对n = 3 的f p u 问题，麻用方法( 3 34 1 ，有如下数值实验结果： ( 334 ) f 335 1 ： jh 忡! ；1 i ；i ；卞；r l r = r 言盲育中j r ；一i 茅肓_ ( 的h = 0 0 1 ，= 10 0 0 0 0 0 1 z h( 钟h = 0 0 2 5 ，= 1 0 0 0 0 0 0 l z h 韭匿窑丑盔堂亟堂位监窑! 旦e 女i ! 垃g 直捏区奎差皓：丑挂直挂 ( a ) h = o 0 1 ，= 1 0 0 0 0 0 0 1 ，r h( m = n 位5 ，出= 1 , 0 0 0 0 0 0 1 z c h j e 夏至堕厶芏熊主芏笪监塞3 旦! ! j i q b 止捏且皇簋浊、盟控直些 c 1h = 0 0 1 ，m = 2 0 0 0 0 0 0 1 t f h 4 。r i i 女mm “t ( d ) 如晒，o d - 2 0 0 0 0 0 0 1 r f h ( e ) h = o o l ，= 2 5 ，r h( o = o t e 5 ，f o = 2 5 z ：h 图3 2 方法f 335 1 的能量误差图 图3 1 、32 是h 分别取0o l ，00 2 5 对方法( 334 1 、( 3 35 ) 的能量误差，由图3 1 口t 以看出，当 m 接近2 口时，f 3 34 ) 的能量误差较大；由图3 2 可以看出，当枷接 近口和2 口时，方_ ：击( 3 3 5 ) 的能量误差都较小。 对于脉冲法( 辛方法) ，我们在第二部分图2 3 中已经给出了它的能量误差为 了进一步说明脉冲法的缺陷，可在取相同步长、不同h o j 的条件下，做出脉冲法的 最大能量误差图，以考察其能量保守性。取步长h = 0l ，实验结果如下： 图3 3 脉冲法的最大能量误差图 1 6 j 夏至通盘生甄土生擅监毫3 旦! ! ! 1 9 直握厘茔篡挂：盟琏直挂 图3 3 是 = 01 时脉冲法的最大能量误差图，其横坐标表示h o j 的取值，纵坐 标表示最大能量误差。从图3 3 可以看出，脉冲法在远离”的整数倍时有较小的能 量误差，但在接近f 的整数倍时误差较大。 类似于图3 3 的情形，方法( 334 ) 、( 335 ) 的最大能量误差图为： 抽) 方法( 334 ) 的最大能量误差罔( b ) 方法( 335 ) 的最太能量误差图 图3 a ( 334 ) 、( 335 ) 两种方法的最人能量误差图 比较图2 3 、3 1 、3 2 、3 3 、3 4 可以看出，当 m 接近口的奇数倍时，方法( 334 ) 的能量误差明显小于脉冲法，但是当 口接近口的偶数倍时，方法( 334 ) 的能量误 差仍然较大；而方法( 335 】当 m 接近z 的整数倍时，能量误差都明显小于脉冲法。 通过以上数值实验可以看出，这两个数值解法都比脉冲法有了较为明显的改 进，但当枷接近z 的偶数倍时，方法034 ) 的能量误差仍然较大，而方法( 3 3 5 ) 具有较好的能量保守性。 h a i r e r 等对线性和非线性高振荡问题( 2 2 2 ) 的分析及一些数值实验结果显示， 具有对称性的方法( 33 1 ) ，当( f ) = s i n c 2 ( f ) ( f ) 时，其能量误差较小口】。在此条 件的基础上，文 2 1 1 0 0 给出一个对称数值解法如下，即选取 爿旬= d n “自，( 封= 出，( 旬，蝎( 白= d n ? ( 旬慨 旬，嘶( 0 = 血，( 0 ， ( 3 36 ) 类似于脉冲法，对h = 3 的f p u 问题- 应用方法( 336 ) ，数值实验结果如下 j e 基窒盈去坐亟芏位监毫3h 自女i l 些直程屋圭差掂：盟鏊左选 im 齑寄；j i # 苜_ ” ( 曲h = o 【0 1 = 1 0 0 0 0 0 0 1 n h ( 0h = o 0 1 ，m = 2 0 0 0 0 0 0 1 r c h棚) h - - 0 ( 2 5 ，= 2 0 0 0 0 0 0 1 m h ( 0h = o o l ，a j = 2 5 n l h( oh - - o ( p 5 ，a j = 2 5 n h 图3 5 方法( 33 回的能量误差图 方法( 3 3 国的最大能量误差图： 图3 6 疗法( 336 ) 的最大能量误差罔 由图3 5 、3 6 可以看出，当 接近f 的整数倍时，方法( 33 6 ) 的能量误差结 果明显优于脉冲法，一致地具有较小的能量误差。 通过以上数值实验结粜，我们可以看到数值解法( 335 ) 、( 3 36 ) 具有较好的能 量保守性。 4 线性高振荡方程的m a g n u s 和n e u m a n n 展开方法 4 1m a g n u s 展开方法 4 1 1 指数导数和它的逆 下面简单介绍一下在推导m a g n u s 展开时需要用到的指数导数的知识f 2 2 2 1 。 指数导数和它的逆可以用矩阵换位子【q ，a 】_ d _ a - a q 来表示。我们假设q 是 固定的，这样就定义了一个线性算子彳h 【q ，a 】，贝, l ja d n ( a ) - - i n ，彳】称为伴随算子。 群的导数为 ( 嘉q ) 日= h q 扣l + q 日q 扣2 + + q 扣1 日 其中q 和h 一般是不可交换的。 当尼= 2 ，k ：3 时，我们有 d h + h f 2 = 2 h q + 峨( ) ， h f 2 2 + 鲫q + q 2 h = 3 h f 2 2 + 3 ( a d o ( h ) ) f 2 + a d 三( 日) 这里口以是线性算子口如的f 次复合。我们约定口以( h ) = 日。对q “1 = f 2 q 应 用l e i b n i z i 去贝l j ，并利用q ( 口以( 日) ) = ( 口以( 日) ) q + 口站1 ( 日) ，采用数学归纳法可 以证明 ( 杀群) 日= 乳冰k 以( h ) 矿1 ， x - ( 4 1 1 ) 两端乘( 尼! ) 卅相加，a 8 4 - j = k i 一1 ，交换求和号顺序，计算后可得： ( 杀唧qh = ( d e 冲。( h ) ) e 坤q ， ( 4 2 ) 其中 d e x p 。( 日) 2 篆巧击矿碰( h ) ( 4 1 3 ) 引理4 1 1 如果线性算子口如的特征值不同于2 ，万f ，其中z 1 ，控，) ，那 么d e x p q 是可逆的而且，对所有的l | q 0 o 证明 设a d n 的特征值是五，则de x p q 的特征值是= 名( 七+ 1 ) ! = t 加 ( 一1 ) 五。由假设的值非零， 此d e x p n 是可逆的。由伯努利数的定义，结合 ( 4 1 3 ) 可得等式( 4 1 4 ) 。r ol l 口如i i 20 q 0 和函数衫( 矿一1 ) 的幂级数的收敛半径是2 万 可知对所有的0 q 0 t n ，彳= 彳( 乙+ h 2 ) ， ( 4 2 4 ) 则x ( f 一乙) 满足方程 x ( t - t ) = b ( t ) x ( t - t ) ，x ( o ) = 】，( 乙) ， ( 4 2 5 ) 其中 曰( f ) = p 一吲j ( 彳( f ) 一j ) 训 计算可得 训j = ( 一：溉端w 刚- s i n ( w 心( t - t 料) w = 、一g ( t + h 2 ) 令t = f 一乞，( 4 2 5 ) 可写成 x 7 ( 丁) = 曰( 乙+ 丁) x ( r ) ，r 0 工( 丁) = y ( ) + r 口( r 。+ x 皿y ( ) + r r b + 王p ( r + 乇) 如凼y ( ) + 若t = 则t = h ，上式变为 x ( ) = y ( 。) + r 口( t ，+ x 皿r ( ) + r rb ( t + x t p “+ x 2 ) d x j x r ( 。) + 其中 眠+ x ) = v ( z ) l w ts i n ( 。2 ( w 。x ) ) 2 w 2 s 。i 。n ( 2 2 ( 驯w x ) 2 ，咄) = 地叫一 设= 2 w h ，则有 r 础。+ z 胁2 ；e 苫j f v c h 胁一；l ? j f v c 蛐c o s c x 冲 + 就轴c h ，s i n ( , u x ) 出 f 427 1 对于a i r y 方程，取h = o1 在r 426 冲取两项进行计算，整体截断误差图为 o 5 。i - o5 o 删l i 1 0 0 02 0 0 03 0 0 04 0 0 05 0 0 0 罔4 2 ：a i r y 冉程的整体截断蜞差 由图42 可见，基于n e u m a n n 展开的数值解法可长时间给出好的计算结果 5 迭代法求解非线性高振荡微分方程 方程( 2 2 2 ) 的准确解形式如下： ) _ ( 篇等心+ f ( 拦五韫如肛 t ， ； = ( 恶等心+ f ( 舞i 韫柚妫凼 2 ， c o s t l - i 一蚓倒咖， 3 ， 么= 百s i n t f l ，召= ( s i n 堡2 马2 2 f q 把而( f ) 代k ( 5 2 ) q h 第- 个方程，对积分使用梯形公式可得 ：乏o ) = c o s f q 蠕一q s i n f q +c o s f q g ( 而) + g ( x a ( f ) ) ) - ( 恶 考虑到数值解法，我们有 一qs i n tl)-xo，j+cost q x n八， = ( 盎 乏( c o s f a g ( x o ) + g ( 加) ) ) i t 2 船( 而) 趣c 0 8 i s h 剖卧 q 八兄j 2 6 互h ( c 。s h 嗡。+ 岛+ 。) 了h 2 s m 。c 2l 了h f 2 ) 岛 2、2 “ ( 5 5 ) 对于”= 3 的f p u 问题，选取初值u ( o ) = ( 1 ，0 , 0 ) ，( o ) = ( 1 ，0 ，0 ) v ( o ) = ( 甜，0 ，0 ) ，v ( o ) = ( 1 ，0 ，0 ) 。在不同的枷取值条件下计算( 55 ) 的能量误差 则有以下一些实验结果； m f a ) h = o 0 1 = 1 0 0 0 ( 0 ) 0 i r e h ( c ) h = 0 0 1 ，c o = 2o f l 0 0 0 0 1 # h 。j j 。o1 m ， m ) = q 晒，o j = 1 0 0 0 0 0 0 1 n h d 1 = q 晒m = 20 0 0 0 0 0 1 r h ：厂一 ：7 一1 。i j 喜m 差懈i ! l l i i i l l一1一i 1 ( 0h = o 0 1 ，= 2 斩，h( 0h = 0 ( 1 2 5 ，= 2 5 ：r h 图5 l 方法幅5 ) 的能量误差囤 图5 l 是h 分别取0 、0 0 2 5 时方法( 5 5 ) 的能量误差。通过数值实验可知当 由 接近口和2 口时，( 5 5 ) 能量误差也很大。 基于对迭代法的修正，我们考虑数值解法 一f c o s h o i 钿j _ l 盯1 s i n h q 其巾晶= g ( m h ) ，中= ( q ) ，甲= 妒( q ) ，甲o = t u o ( h n ) ，甲= ( n ) 由( 332 ) 可知，数值解法( 56 ) 当满足 ( f ) s l n c 2 ( 婴) ：s l n c ( f ) ，： 条件时( ，为偶函数) 方法