1、1Hamilton力学的辛算法和分子动力学模拟 陈敏伯中国科学院上海有机化学有机所计算化学课题组2006年10月2内容 冯康对世界科学的重大贡献 Euclid空间 辛空间 Hamilton力学的辛结构 正则变换的辛结构 辛算法应用实例3 Schrdinger：“Hamilton原理已经成为现代物理学的基石。” Hamilton原理将不同的物理规律纳入了统一的数学形式。 现在问题就归结到：怎样才能对Hamilton力学的运动方程作正确的数值计算。 一切Hamilton体系的动力学演化都使辛度量保持不变，即都是辛（正则）变换。 一切解Hamilton方程 “正确”的离散算法都应当是辛变换的。 （冯
2、康，1997年国家自然科学一等奖“哈密尔顿系统辛几何算法” ） Lax：“他的声望是国际性的。” 丘成桐：“中国在数学历史上很出名的有三个：一个是陈省身教授在示性类方面的工作，一个是华罗庚在多复变函数方面的工作，一个是冯康在有限元计算方面的工作。” (1998年3月11日中国科学报) 4“冯氏大定理” 同一物理定律的不同的数学表述，尽管在物理上是等价的；但在计算上是不等价的。 冯康：如果在算法中能够保持辛几何的对称性，将可避免人为耗散性这类算法的缺陷，成为具有高保真性的算法。 在天体力学的轨道计算，粒子加速器中的轨道计算和分子动力学计算中得到广泛的应用。5冯康（19201993）的学术成就 1
3、965年发表论文“基于变分原理的差分格式”。国际学术界承认冯康独立发展了有限元方法。（仅获1982年国家自然科学二等奖。冯康得悉非常难过，曾打算将申请撤回。） 前国际数学会理事长J. L. Lions教授1981年说：“中国学者在对外隔绝的环境下独立创造了有限元，在世界上是最早之列。今天这一贡献已为全人类所共享。” 1984年以后创建的“哈密尔顿系统的辛几何算法”。 （1991年评为国家自然科学奖二等奖。冯康获悉后撤回申请。直到1997年底，在冯康去世四年之后，终于授予了国家自然科学一等奖。 ） 石钟慈：“国际上最早系统地研究并建立辛几何算法的。” 6数学地位现代微分几何规范场理论微分拓扑辛几
4、何.线性空间线性泛函对偶空间多重线性泛函张量空间张量分析流形理论7外微分 辛几何 辛几何的基础是外微分形式。 外微分形式是如下概念推广到高维的产物： 1、作功在场中沿某一路径所作的功； 2、流量单位时间内流体穿过某曲面的量 3、面积或体积平行四边形面积 或平行六面体体积。 外微分形式中有“1-形式”、“2-形式”等 辛构造就是非简并的闭2-形式。8Euclid空间 对称性： 线性： （ k 为任意实数） （ c是V中的任意向量） 非简并性： ，当且仅当 时才 ,0a aa0,0a a,a bb a ,ac ba bc b,kka ba b符合如下内积定义的线性空间V称为“Euclid空间”。然
5、后就可以给出向量的长度、正交、单位向量等概念。 9辛空间（Simplectic Space） 反对称性： 双线性： 非简并性：若向量a对于W中的任意向量b均有 ，则 具有如下内积定义的线性空间W为“辛空间” 。 这种内积称为“辛内积”。, a bb a1212,aa ba ba b1212,a bba ba b,0a ba010辛空间 度量：作功、面积（或体积）、流量等 辛内积： 2维： a、b平行四边形面积 2n维：1212,aabba b122.Tna aaa122.Tnb bbb21,nnin in iiiababa ba J b2nnn01J10单位辛矩阵 ：1221112122112
6、2,TTTTnnTTTnin in iiinababaa01aabba ba J ba ba b10bb11单位辛矩阵 的性质 若 A 为对称阵，且 ，则 2nJ2 J11T JJJ20TnR a Jaa2nBJ A22TnnB JJ B0222nnnnnnnn 000J0001111111证明：12Euclid空间和辛空间的对应关系 Euclid空间辛空间内积 长度内积 面积单位矩阵 单位辛矩阵 正交辛正交正交归一基共轭辛正交归一基正交矩阵辛正交矩阵对称变换 Hamilton变换 实对称矩阵的本征值均为实数若Hamilton矩阵的本征值为 ，则 也是它的本征值实对称矩阵的不同本征值的本征向量
7、必正交Hamilton矩阵的非辛共轭本征值的本征向量必辛正交实对称矩阵的所有本征向量组成一组正交归一基Hamilton矩阵的所有本征向量组成一组共轭辛正交归一基, a b, a bTTC CCC11TS JSJ,0TTa ba bab11J,0Ta ba Jb,abbaA AA A,abbaH HH H13Hamilton力学的辛结构 12.Tfq qqq12.Tfp ppp112 . .Tfffzz zz pzq12 . TfHHHzzz11 . . Tffpp qq pzq 1HHHHH 01ppqJz10qzqp, iiiiHHqpipq 1HJzz14正则变换的辛结构 pzq1HJzz
8、PyQ正则变量从 变换到 记为：yMzydydzMdzzTHHMzy即：yMz1HHHHH01PPQJy10yQQP111TTHHHHH JyMzMJMJMyzyyMJMJMyTMJMJM1HJyy1HJzzM辛变换： 1T JJJyHHHMzzyy15正则变换M的性质1TMJMJ、2 det1M2、16无穷小辛阵 定义：若 ， 则该2n阶矩阵 称为“无穷小辛阵” 设 为对称阵，当且仅当 时 ， 为无穷小辛阵。 （证明略） 若 为无穷小辛阵， 则 为辛阵 。 若 为无穷小辛阵，又若 非奇异， 则 为辛阵 22TnnB JJ B0BC2nBJ CBeBBBB 1BB1117辛阵2、当且仅当 和
9、，则 、 都为辛阵 3、 是辛阵 4、当且仅当 ，则 是辛阵 5、当且仅当 和 ，则 是辛阵 TBBTDDnnB011nn0D111TB00BTTMJMNJN1M N2QQTQQnnQ1Q1QQ1、 是辛阵的充要条件：ABCD,TTTTTTTTTTTTnnABBA0CDDC0A CC A0B DD B0ADBCA DC B1122TnnS J SJ18线性Hamilton体系的辛差分格式 线性Hamilton体系Hamilton函数是 的二次型 pzq 12THzz CzTCC且11111122TTTdHdt zzJJz CzJCzz CzJzzCddtzzBz其中1BJ C为无穷小辛阵 Be
10、B为辛阵 积分19 100lnln0ttddttzzzzBBzz 0tteBzz20中点Euler法的辛格式1112mmmmhzzB zzddtzzBz111122hmmmhhzBF zBzh 为时间步长111222hhhhFBBB 111 因为 为无穷小辛阵，且非奇异即 ，故步进算符 为辛阵，故为辛格式。 2hBdet02hB111122hhhFBB21可分、线性Hamilton体系的中点Euler公式 “可分、线性Hamilton体系” ,Hq pT pV q111,222TTTTH pqT0qq pp Tpq Vq0VpT0C0V1BJ C1nn0T00VBJ C00VT01122111
11、2222nnmmmmnnhhhh VVppqqTT1111Euler中点法1112mmmmhzzB zz1111112mmmmmmmmhqqqq0VppppT0演绎见后页23演绎细节：111122mmmmmmmmhhqqV ppppT qq1111112mmmmmmmmhqqqq0VppppT01111121 2mmmmmmmmhh qqV ppppT qq11112222mmmmmmmmhhhhqVpqVppTqpTq11112222mmmmmmmmhhhhqVpqVppTqpTq112222nnmmmmnnhhhhVVppqqTT11112411112222mmmmmmnnhnnhhhhp
12、pzVVFTTzqq111112222nnhnnhhhh VVFTT1111前面我们已经证明了 是辛阵，所以上面算法是辛格式。11122hhhFBB25基于Pad逼近的辛格式 线性Hamilton体系相流 有理Pad逼近： 0lltgtzB z 0tteBzz lmxlmlmnxegxdx 00! ! ! ! ! !mklmklklmklmkmnxxlm kmklmkldxxlm klk称为“l+m阶对ex的Pad逼近” 00lllltgttgtqB qpB p即2n0VBJ ST0可分体系：26用以下 构造的差分格式都是辛格式 llgx, l l0,01,12,23,34,423231210
13、1201210120 xxxxxx2342343122884168031228841680 xxxxxxxx2212121212xxxx1212xx11 llgx*： “(1,1) 逼近” ，就是Euler中点格式。 11gx精度2阶4阶6阶8阶27可分线性Hamilton体系的交叉显式辛格式 11,22TTHq pp Tpq Vq111222TTTH pp Tpq Vqq VqVqqqq111222TTTHqp Tpq Vqp TpTpppp11231122 11 mmmmmmhh pVqqTpppVqqqTp差分11231122 mmmmmmhhppVqqqTp2813122 mmhmmp
14、pFqq1nnhnnhh0VFT01111当且仅当 和 时, 和 都为辛阵。 TTTTVVnnh1V01nnh10T1当且仅当 ，则 是辛阵。现在是 ，所以也是辛阵。故为辛格式。 TTMJMNJN1M NTTMJMNJNJ1nnnnhh101VT10129演绎细节：11231122 mmmmmmhhppVqqqTp11213122mmmmmmhhpVqpTpqq13122 mmnnmmnnhhpp0VqqT0111111121232 mmnnmmnnmhmhhpp0qqqV0pFT111130可分线性Hamilton体系的交叉显式辛格式 3112211211 mmmmmmhh qqTpppVq
15、 qTppVq差分31122112mmmmmmhhqqTpppVq13122mmhmmppFqq1nnhnnhh101VFT101h：时间步长31验证：113122mmnnmmnnhhpp101VqqT10113122mmnnmmnnhhpp101VqqT10111213122mmmmmmhhpVqpTpqq11231122mmmmmmhhppVqqqTp验证完毕32实例1 谐振子的相空间轨迹 (a) Runge-Kutta法 3000步 ，步长0.4。人为耗散，轨道收缩 (b) Adams法 步长0.2 。人为反耗散，轨道发散 (c)蛙跳法 * 步，步长0.1 。初、中、末各取三段1000步
16、的结果完全吻合 71033实例2 非谐振子的相空间轨迹（a）与（b）为同一个蛙跳法模拟的分段取样结果 (c) 二阶辛算法1000步 .初、中、末三段结果完全吻合 最初1000步 轨道失真 第900010000步 轨道继续失真 *：蛙跳法即二步中心差分法，它对于非线性方程不是辛算法 34实例3 Huygens振子 (a) Runge-Kutta法 步长0.10000005 ；9x105步 。趋于左吸引子 (b) Runge-Kutta法 步长0.10000004 ；9x105步 。趋于右吸引子 (c) 二阶辛算法4条轨道，每条各108步；步长0.1 每条轨道的初、中、末各取三段500步的结果完全
17、吻合。具有超长期跟踪能力 位于双纽线之外的任意初始相点趋于左右两个假吸引子的几率相同。 35实例4 椭球面上的测地线 5 4(a) Runge-Kutta法轨道不趋稠密步长0.05658 ，104步频率比：(b)辛算法轨道趋于趋稠密无理数36实例5椭球面上的测地线步长0.033427 ，105步周期：25频率比：11/16有理数(a) Runge-Kutta法轨道不封闭(b)辛算法轨道封闭37实例6Kepler轨道 当频率比为有理数时，应当形成封闭轨道 步长0.01605 ，2.5x105步频率比：11/20有理数(a) Runge-Kutta法轨道不封闭(b)辛算法轨道封闭38实例7Li2
18、分子的经典轨迹法 设原子位置 折合质量 广义位置 广义动量 动能 势能取Morse势 Hamiltian量12xx、1212m mmm12qxxdqpdt 22pT p 22eea q qa q qV qD ee ,H q pT pV q39Li2分子的经典轨迹的正则方程 dpdVf qdtdqdqdTg pdtp Li2分子 态的参数： ， ， -1。设初态为：步长0.005。1gX18541 cmD2.67328 eq 0.867a 0eqq 020.0001pD401. 振幅、周期(a) 辛算法：长达106步时还保持振幅恒定，周期性恒定。(b) Runge-Kutta法：5000步之后振
19、幅变小、周期变短 412.相空间轨迹(a) 辛算法：长达106步时还保持总能量恒定、相空间轨迹稳定、。(b) Runge-Kutta法：104步之后总能量急剧下降 ；相空间轨迹沿q方向收缩，5104步时已经面目全非。 42实例说明： 8种实例：简谐振子、Duffing振子（非线性振子）Huygens振子、Cassini振子、二维多晶格与准晶格定常流、Lissajous图形、椭球面测地线流、Kepler运动。 说明了在整体性、结构性和长期跟踪能力上辛算法的优越性。一切传统非辛算法，无论精度高低均无例外地全然失效。一切辛算法无论精度高低均无例外地过关，均具有长期稳健的跟踪能力。显示了压倒性的优越性
20、。43Hamilton体系的守恒律辛算法保持了Hamilton体系具有的两个守恒律 ： 1、相空间体积的不变性 Liouville-Poincar守恒律 2、运动不变量：如能量、动量、角动量的守恒 辛算法能够在数值计算中保持辛变换的结构，于是就会得到高的稳定性。 辛算法的差分方法被认为是目前最稳定、高效的计算方法。最适合用于经典力学体系。 辛算法不含人为耗散性，先天性地免于一切非哈污染，是“干净”的算法。44 传统算法除了极个别例外，均为非辛算法。大都是为了渐近稳定体系设计的，都含有耗散机制以保证计算稳定性。 Hamilton体系不具有渐近稳定性，所以传统算法都不可避免地带进人为耗散性、虚假吸
21、引子及其它种种非哈体系本身具有的寄生效应。45Refs1 余扬政，冯承天，物理学中的几何方法，高等教育出版社，施普林格出版社，1998年。2 Arnold, V. I., Mathematical Methods of Classical Mechanics, Springer-Verlag, Heidelberg, 1978; 中译本：齐民友译，经典力学的数学方法（第4版），高等教育出版社，北京，1992年（中译本第1版），2006年（中译本第2版）。3 李德明，陈昌民，经典力学，高等教育出版社，北京，2006年，第4章。4 冯康，秦孟兆，哈密尔顿系统的辛几何算法，浙江科学技术出版社，杭州，2003年。5 余德浩，汤华中，微分方程数值解法（中国科学院研究生院教材），科学出版社，2003年。6 姚伟岸，钟万勰，辛弹性力学，高等教育出版社，北京，2002年。
