Hamilton力学的辛算法课件

Hamilton力学的辛算法

和分子动力学模拟陈敏伯中国科学院上海有机化学有机所计算化学课题组2006年10月1Hamilton力学的辛算法

和分子动力学模拟陈敏伯1内容冯康对世界科学的重大贡献Euclid空间辛空间Hamilton力学的辛结构正则变换的辛结构辛算法应用实例2内容冯康对世界科学的重大贡献2Schrödinger：“Hamilton原理已经成为现代物理学的基石。”Hamilton原理将不同的物理规律纳入了统一的数学形式。现在问题就归结到：怎样才能对Hamilton力学的运动方程作正确的数值计算。一切Hamilton体系的动力学演化都使辛度量保持不变，即都是辛（正则）变换。一切解Hamilton方程“正确”的离散算法都应当是辛变换的。（冯康，1997年国家自然科学一等奖“哈密尔顿系统辛几何算法”）Lax：“他的声望是国际性的。”丘成桐：“中国…在数学历史上很出名的有三个：一个是陈省身教授在示性类方面的工作，一个是华罗庚在多复变函数方面的工作，一个是冯康在有限元计算方面的工作。”(1998年3月11日《中国科学报》)3Schrödinger：“Hamilton原理已经成为现代物“冯氏大定理”同一物理定律的不同的数学表述，尽管在物理上是等价的；但在计算上是不等价的。冯康：如果在算法中能够保持辛几何的对称性，将可避免人为耗散性这类算法的缺陷，成为具有高保真性的算法。在天体力学的轨道计算，粒子加速器中的轨道计算和分子动力学计算中得到广泛的应用。4“冯氏大定理”同一物理定律的不同的数学表述，尽管在物理上是等冯康（1920－1993）的学术成就1965年发表论文“基于变分原理的差分格式”。国际学术界承认冯康独立发展了有限元方法。（仅获1982年国家自然科学二等奖。冯康得悉非常难过，曾打算将申请撤回。）前国际数学会理事长J.–L.Lions教授1981年说：“中国学者在对外隔绝的环境下独立创造了有限元，在世界上是最早之列。今天这一贡献已为全人类所共享。”1984年以后创建的“哈密尔顿系统的辛几何算法”。（1991年评为国家自然科学奖二等奖。冯康获悉后撤回申请。直到1997年底，在冯康去世四年之后，终于授予了国家自然科学一等奖。）石钟慈：“国际上最早系统地研究并建立辛几何算法的。”5冯康（1920－1993）的学术成就1965年发表论文“基于数学地位6数学地位6外微分辛几何

辛几何的基础是外微分形式。

外微分形式是如下概念推广到高维的产物：1、作功—在场中沿某一路径所作的功；2、流量—单位时间内流体穿过某曲面的量3、面积或体积—平行四边形面积或平行六面体体积。

外微分形式中有“1-形式”、“2-形式”等辛构造就是非简并的闭2-形式。7外微分辛几何辛几何的基础是外微分形式。7Euclid空间

对称性：

线性：（k为任意实数）（c是V中的任意向量）

非简并性：，当且仅当时才符合如下内积定义的线性空间V称为“Euclid空间”。然后就可以给出向量的长度、正交、单位向量等概念。8Euclid空间对称性：符辛空间（SimplecticSpace）

反对称性：

双线性：

非简并性：若向量a对于W中的任意向量b均有，则具有如下内积定义的线性空间W为“辛空间”。这种内积称为“辛内积”。9辛空间（SimplecticSpace）反对称性：具辛空间度量：作功、面积（或体积）、流量等辛内积：2维：a、b平行四边形面积2n维：单位辛矩阵：10辛空间度量：作功、面积（或体积）、流量等单位辛矩阵：10单位辛矩阵的性质

若A为对称阵，且，则证明：▌11单位辛矩阵的性质证明：▌11Euclid空间和辛空间的对应关系Euclid空间辛空间内积——长度内积——面积单位矩阵单位辛矩阵正交辛正交正交归一基共轭辛正交归一基正交矩阵辛正交矩阵对称变换Hamilton变换实对称矩阵的本征值均为实数若Hamilton矩阵的本征值为，则也是它的本征值实对称矩阵的不同本征值的本征向量必正交Hamilton矩阵的非辛共轭本征值的本征向量必辛正交实对称矩阵的所有本征向量组成一组正交归一基Hamilton矩阵的所有本征向量组成一组共轭辛正交归一基12Euclid空间和辛空间的对应关系Euclid空间辛空间内Hamilton力学的辛结构13Hamilton力学的辛结构13正则变换的辛结构正则变量从变换到记为：即：M辛变换：14正则变换的辛结构正则变量从变换到正则变换M的性质15正则变换M的性质15无穷小辛阵定义：若，则该2n阶矩阵称为“无穷小辛阵”设为对称阵，当且仅当时，为无穷小辛阵。（证明略）若为无穷小辛阵，则为辛阵。若为无穷小辛阵，又若非奇异，则为辛阵16无穷小辛阵定义：若辛阵2、当且仅当和，则、都为辛阵3、是辛阵4、当且仅当，则是辛阵5、当且仅当和，则是辛阵1、是辛阵的充要条件：17辛阵2、当且仅当和线性Hamilton体系的辛差分格式线性Hamilton体系——Hamilton函数是的二次型且其中为无穷小辛阵为辛阵积分18线性Hamilton体系的辛差分格式线性Hamilton体系1919中点Euler法的辛格式h为时间步长因为为无穷小辛阵，且非奇异即，故步进算符为辛阵，故为辛格式。20中点Euler法的辛格式h为时间步长因为为无穷小可分、线性Hamilton体系的中点Euler公式

——“可分、线性Hamilton体系”21可分、线性Hamilton体系的中点Euler公式21Euler中点法演绎见后页22Euler中点法演绎见后页22演绎细节：23演绎细节：23前面我们已经证明了是辛阵，所以上面算法是辛格式。24前面我们已经证明了基于Padé逼近的辛格式线性Hamilton体系相流有理Padé逼近：称为“l+m阶对ex的Padé逼近”即可分体系：25基于Padé逼近的辛格式线性Hamilton体系相流称为“用以下构造的差分格式都是辛格式*：“(1,1)逼近”，就是Euler中点格式。精度2阶4阶6阶8阶26用以下构造的差分格式都是辛格式*：“(可分线性Hamilton体系的交叉显式辛格式差分27可分线性Hamilton体系的交叉显式辛格式差分27当且仅当和时,和都为辛阵。当且仅当，则是辛阵。现在是，所以也是辛阵。故为辛格式。28当且仅当和时,演绎细节：29演绎细节：29可分线性Hamilton体系的交叉显式辛格式差分h：时间步长30可分线性Hamilton体系的交叉显式辛格式差分h：时间步验证：验证完毕31验证：验证完毕31实例1－谐振子的相空间轨迹

(a)Runge-Kutta法

3000步，步长0.4。人为耗散，轨道收缩(b)Adams法

步长0.2。人为反耗散，轨道发散(c)蛙跳法

*步，步长0.1。初、中、末各取三段1000步的结果完全吻合32实例1－谐振子的相空间轨迹(a)Runge-Kutta实例2－非谐振子的相空间轨迹（a）与（b）为同一个蛙跳法模拟的分段取样结果(c)二阶辛算法1000步.初、中、末三段结果完全吻合最初1000步轨道失真第9000－10000步轨道继续失真*：蛙跳法即二步中心差分法，它对于非线性方程不是辛算法33实例2－非谐振子的相空间轨迹（a）与（b）为同一个蛙跳法模实例3－Huygens振子

(a)Runge-Kutta法

步长0.10000005；9x105步。趋于左吸引子

(b)Runge-Kutta法

步长0.10000004；9x105步。趋于右吸引子

(c)二阶辛算法4条轨道，每条各108步；步长0.1每条轨道的初、中、末各取三段500步的结果完全吻合。具有超长期跟踪能力

位于双纽线之外的任意初始相点趋于左右两个假吸引子的几率相同。34实例3－Huygens振子(a)Runge-Kutta实例4－椭球面上的测地线

(a)Runge-Kutta法轨道不趋稠密步长0.05658，104步频率比：(b)辛算法轨道趋于趋稠密无理数35实例4－椭球面上的测地线(a)Runge-Kutta法实例5－椭球面上的测地线步长0.033427，105步周期：25频率比：11/16有理数(a)Runge-Kutta法轨道不封闭(b)辛算法轨道封闭36实例5－椭球面上的测地线步长0.033427，105步有实例6－Kepler轨道当频率比为有理数时，应当形成封闭轨道步长0.01605，2.5x105步频率比：11/20有理数(a)Runge-Kutta法轨道不封闭(b)辛算法轨道封闭37实例6－Kepler轨道当频率比为有理数时，应当形成封闭轨实例7－Li2分子的经典轨迹法设原子位置折合质量广义位置广义动量动能势能取Morse势Hamiltian量38实例7－Li2分子的经典轨迹法设原子位置Li2分子的经典轨迹的正则方程Li2分子态的参数：，Å，Å-1。设初态为：步长0.005。39Li2分子的经典轨迹的正则方程Li2分子1.振幅、周期(a)辛算法：长达106步时还保持振幅恒定，周期性恒定。(b)Runge-Kutta法：5000步之后振幅变小、周期变短401.振幅、周期(a)辛算法：(b)Runge-Kutt2.相空间轨迹(a)辛算法：长达106步时还保持总能量恒定、相空间轨迹稳定、。(b)Runge-Kutta法：104步之后总能量急剧下降；相空间轨迹沿q方向收缩，5104步时已经面目全非。412.相空间轨迹(a)辛算法：(b)Runge-Kutta实例说明：8种实例：简谐振子、Duffing振子（非线性振子）Huygens振子、Cassini振子、二维多晶格与准晶格定常流、Lissajous图形、椭球面测地线流、Kepler运动。说明了在整体性、结构性和长期跟踪能力上辛算法的优越性。一切传统非辛算法，无论精度高低均无例外地全然失效。一切辛算法无论精度高低均无例外地过关，均具有长期稳健的跟踪能力。显示了压倒性的优越性。42实例说明：8种实例：简谐振子、Duffing振子（非线性振子Hamilton体系的守恒律辛算法保持了Hamilton体系具有的两个守恒律：1、相空间体积的不变性——Liouville-Poincaré守恒律2、运动不变量：如能量、动量、角动量的守恒辛算法能够在数值计算中保持辛变换的结构，于是就会得到高的稳定性。辛算法的差分方法被认为是目前最稳定、高效的计算方法。最适合用于经典力学体系。辛算法不含人为耗散性，先天性地免于一切非哈污染，是“干净”的算法。43Hamilton体系的守恒律辛算法保持了Hamilton体系传统算法除了极个别例外，均为非辛算法。大都是为了渐近稳定体系设计的，都含有耗散机制以保证计算稳定性。Hamilton体系不具有渐近稳定性，所以传统算法都不可避免地带进人为耗散性、虚假吸引子及其它种种非哈体系本身具有的寄生效应。44传统算法除了极个别例外，均为非辛算法。大都是为了渐近稳定体系Refs[1]余扬政，冯承天，《物理学中的几何方法》，高等教育出版社，施普林格出版社，1998年。[2]Arnold,V.I.,MathematicalMethodsofClassicalMechanics,Springer-Verlag,Heidelberg,1978;中译本：齐民友译，《经典力学的数学方法》（第4版），高等教育出版社，北京，1992年（中译本第1版），2006年（中译本第2版）。[3]李德明，陈昌民，《经典力学》，高等教育出版社，北京，2006年，第4章。[4]冯康，秦孟兆，《哈密尔顿系统的辛几何算法》，浙江科学技术出版社，杭州，2003年。[5]余德浩，汤华中，《微分方程数值解法》（中国科学院研究生院教材），科学出版社，2003年。[6]姚伟岸，钟万勰，《辛弹性力学》，高等教育出版社，北京，2002年。45Refs[1]余扬政，冯承天，《物理学中的几何方法》，高等Hamilton力学的辛算法

和分子动力学模拟陈敏伯中国科学院上海有机化学有机所计算化学课题组2006年10月46Hamilton力学的辛算法

和分子动力学模拟陈敏伯1内容冯康对世界科学的重大贡献Euclid空间辛空间Hamilton力学的辛结构正则变换的辛结构辛算法应用实例47内容冯康对世界科学的重大贡献2Schrödinger：“Hamilton原理已经成为现代物理学的基石。”Hamilton原理将不同的物理规律纳入了统一的数学形式。现在问题就归结到：怎样才能对Hamilton力学的运动方程作正确的数值计算。一切Hamilton体系的动力学演化都使辛度量保持不变，即都是辛（正则）变换。一切解Hamilton方程“正确”的离散算法都应当是辛变换的。（冯康，1997年国家自然科学一等奖“哈密尔顿系统辛几何算法”）Lax：“他的声望是国际性的。”丘成桐：“中国…在数学历史上很出名的有三个：一个是陈省身教授在示性类方面的工作，一个是华罗庚在多复变函数方面的工作，一个是冯康在有限元计算方面的工作。”(1998年3月11日《中国科学报》)48Schrödinger：“Hamilton原理已经成为现代物“冯氏大定理”同一物理定律的不同的数学表述，尽管在物理上是等价的；但在计算上是不等价的。冯康：如果在算法中能够保持辛几何的对称性，将可避免人为耗散性这类算法的缺陷，成为具有高保真性的算法。在天体力学的轨道计算，粒子加速器中的轨道计算和分子动力学计算中得到广泛的应用。49“冯氏大定理”同一物理定律的不同的数学表述，尽管在物理上是等冯康（1920－1993）的学术成就1965年发表论文“基于变分原理的差分格式”。国际学术界承认冯康独立发展了有限元方法。（仅获1982年国家自然科学二等奖。冯康得悉非常难过，曾打算将申请撤回。）前国际数学会理事长J.–L.Lions教授1981年说：“中国学者在对外隔绝的环境下独立创造了有限元，在世界上是最早之列。今天这一贡献已为全人类所共享。”1984年以后创建的“哈密尔顿系统的辛几何算法”。（1991年评为国家自然科学奖二等奖。冯康获悉后撤回申请。直到1997年底，在冯康去世四年之后，终于授予了国家自然科学一等奖。）石钟慈：“国际上最早系统地研究并建立辛几何算法的。”50冯康（1920－1993）的学术成就1965年发表论文“基于数学地位51数学地位6外微分辛几何

辛几何的基础是外微分形式。

外微分形式是如下概念推广到高维的产物：1、作功—在场中沿某一路径所作的功；2、流量—单位时间内流体穿过某曲面的量3、面积或体积—平行四边形面积或平行六面体体积。

外微分形式中有“1-形式”、“2-形式”等辛构造就是非简并的闭2-形式。52外微分辛几何辛几何的基础是外微分形式。7Euclid空间

对称性：

线性：（k为任意实数）（c是V中的任意向量）

非简并性：，当且仅当时才符合如下内积定义的线性空间V称为“Euclid空间”。然后就可以给出向量的长度、正交、单位向量等概念。53Euclid空间对称性：符辛空间（SimplecticSpace）

反对称性：

双线性：

非简并性：若向量a对于W中的任意向量b均有，则具有如下内积定义的线性空间W为“辛空间”。这种内积称为“辛内积”。54辛空间（SimplecticSpace）反对称性：具辛空间度量：作功、面积（或体积）、流量等辛内积：2维：a、b平行四边形面积2n维：单位辛矩阵：55辛空间度量：作功、面积（或体积）、流量等单位辛矩阵：10单位辛矩阵的性质

若A为对称阵，且，则证明：▌56单位辛矩阵的性质证明：▌11Euclid空间和辛空间的对应关系Euclid空间辛空间内积——长度内积——面积单位矩阵单位辛矩阵正交辛正交正交归一基共轭辛正交归一基正交矩阵辛正交矩阵对称变换Hamilton变换实对称矩阵的本征值均为实数若Hamilton矩阵的本征值为，则也是它的本征值实对称矩阵的不同本征值的本征向量必正交Hamilton矩阵的非辛共轭本征值的本征向量必辛正交实对称矩阵的所有本征向量组成一组正交归一基Hamilton矩阵的所有本征向量组成一组共轭辛正交归一基57Euclid空间和辛空间的对应关系Euclid空间辛空间内Hamilton力学的辛结构58Hamilton力学的辛结构13正则变换的辛结构正则变量从变换到记为：即：M辛变换：59正则变换的辛结构正则变量从变换到正则变换M的性质60正则变换M的性质15无穷小辛阵定义：若，则该2n阶矩阵称为“无穷小辛阵”设为对称阵，当且仅当时，为无穷小辛阵。（证明略）若为无穷小辛阵，则为辛阵。若为无穷小辛阵，又若非奇异，则为辛阵61无穷小辛阵定义：若辛阵2、当且仅当和，则、都为辛阵3、是辛阵4、当且仅当，则是辛阵5、当且仅当和，则是辛阵1、是辛阵的充要条件：62辛阵2、当且仅当和线性Hamilton体系的辛差分格式线性Hamilton体系——Hamilton函数是的二次型且其中为无穷小辛阵为辛阵积分63线性Hamilton体系的辛差分格式线性Hamilton体系6419中点Euler法的辛格式h为时间步长因为为无穷小辛阵，且非奇异即，故步进算符为辛阵，故为辛格式。65中点Euler法的辛格式h为时间步长因为为无穷小可分、线性Hamilton体系的中点Euler公式

——“可分、线性Hamilton体系”66可分、线性Hamilton体系的中点Euler公式21Euler中点法演绎见后页67Euler中点法演绎见后页22演绎细节：68演绎细节：23前面我们已经证明了是辛阵，所以上面算法是辛格式。69前面我们已经证明了基于Padé逼近的辛格式线性Hamilton体系相流有理Padé逼近：称为“l+m阶对ex的Padé逼近”即可分体系：70基于Padé逼近的辛格式线性Hamilton体系相流称为“用以下构造的差分格式都是辛格式*：“(1,1)逼近”，就是Euler中点格式。精度2阶4阶6阶8阶71用以下构造的差分格式都是辛格式*：“(可分线性Hamilton体系的交叉显式辛格式差分72可分线性Hamilton体系的交叉显式辛格式差分27当且仅当和时,和都为辛阵。当且仅当，则是辛阵。现在是，所以也是辛阵。故为辛格式。73当且仅当和时,演绎细节：74演绎细节：29可分线性Hamilton体系的交叉显式辛格式差分h：时间步长75可分线性Hamilton体系的交叉显式辛格式差分h：时间步验证：验证完毕76验证：验证完毕31实例1－谐振子的相空间轨迹

(a)Runge-Kutta法

3000步，步长0.4。人为耗散，轨道收缩(b)Adams法

步长0.2。人为反耗散，轨道发散(c)蛙跳法

*步，步长0.1。初、中、末各取三段1000步的结果完全吻合77实例1－谐振子的相空间轨迹(a)Runge-Kutta实例2－非谐振子的相空间轨迹（a）与（b）为同一个蛙跳法模拟的分段取样结果(c)二阶辛算法1000步.初、中、末三段结果完全吻合最初1000步轨道失真第9000－10000步轨道继续失真*：蛙跳法即二步中心差分法，它对于非线性方程不是辛算法78实例2－非谐振子的相空间轨迹（a）与（b）为同一个蛙跳法模实例3－Huygens振子

(a)Runge-Kutta法

步长0.10000005；9x105步。趋于左吸引子

(b)Runge-Kutta法

步长0.10000004；9x105步。趋于右吸引子

(c)二阶辛算法4条轨道，每条各108步；步长0.1每条轨道的初、中、末各取三段500步的结果完全吻合。具有超长期跟踪能力

位于双纽线之外的任意初始相点趋于左右两个假吸引子的几率相同。79实例3－Huygens振子(a)Runge-Kutta实例4－椭球面上的测地线

(a)Runge-Kutta法轨道不趋稠密步长0.05658，104步频率比：(b)辛算法轨道趋于趋稠密无理数80实例4－椭球面上的测地线(a)Runge-Kutta法实例5－椭球面上的测地线步长0.033427，105步周期：25频率比：11/16有理数(a)Runge-Kutta法轨道不封闭(b)辛算法轨道封闭81实例5－椭球面上的测地线步长0.033427，105步有实例6－Kepler轨道当频率比为有理数时，应当形成封闭轨道步长0.01605，2.5x105步频率比：11/20有理数(a)Runge-Kutta法轨道不封闭(b)辛算法轨道封闭82实例6－Kepler轨道当频率比为有理数时，应当形成封闭轨实例7－Li2分子的经典轨迹法设原子位置折合质量广义位置广义动量动能势能取Morse势Hamiltian量83实例7－Li2分子的经典轨迹法设原子位置Li2分子的经典轨迹的正则方程Li2分子态的参数：，Å，Å-1。设初态为：步长0.005。84Li2分子的经典轨迹的正则方程Li2