引 言
RK法是工程中最常应用的数值方法，它格式简单、几何意义明确、运算效率高，但是存在耗散机制，误差会累积，不适宜做长时间的计算.为此，Feng(冯康)［1］在1984年国际微分几何与微分方程会议上发表的《差分格式与辛几何》，首次系统地提出了Hamilton方程及其算法——辛几何算法(symplectic geometric algorithms)，这是根据Hamilton系统的解能保持辛几何结构不变的性质而设计的辛算法.
辛算法虽然克服了RK法的精度耗散机制，然而在计算过程中如RK法一样会产生相位误差，从而影响数值精度.Grtz研究了单步辛算法的相位误差问题［2］，邢誉峰等针对经典的无阻尼Hamilton系统阐述了辛算法导致动态响应“面目全非”的原因——相位误差累积［3］.Monovasilis和Simos等设计了具有最高阶相位精度的某些辛算法的计算格式(譬如SPRK法)等［4-13］.但是这些方法只能有限地提高辛算法在时域上的精度，长时间计算后相位误差的精度还是会大幅度降低.
采用小步长是提高计算精度的另一种途径，但单纯地将步长细化会大大增加计算量.为此，钟万勰设计了精细算法［14］，能够在步长允许的范围内多次计算，避免了大数吃小数，克服了计算机的舍入误差，使得计算精度几乎可以达到计算机的精度，并且精细矩阵生成后可以做到“一次生成，终身使用”，从而计算量并没有增加，但计算精度却大大提高了，因此得到了广泛的应用［15-16］.Huang(黄永安)等［17］、曾进等［18］利用 Padé逼近，将辛算法与精细积分法相结合，设计了精细辛算法，使得解在保持计算长期稳定性的同时也大大提高了计算精度.但是计算过程中涉及矩阵求逆，只能近似保辛［19-20］.
因此，本文从相位误差的角度出发，将精细算法引入到一阶辛差分格式中，设计了基于相位误差的精细辛算法，证明了此算法每一小步是辛算法，合成一大步后仍然是辛算法，当精细参数N趋于无穷大时，相位误差趋于零.因此，所设计的精细辛算法不仅可以保持原Hamilton系统的性质，还使相位精度得到了大幅度的提高，可以有效逼近精确的相位，长时间计算后仍可逼近精确解.
1 Hamilton系统的辛算法
由流形的理论可知，Hamilton方程的解能保持正则变换.即若Hamilton方程的解在t0时刻为z，到达t0+t时刻的解为z，则其Jacobi矩阵M=z/z是辛矩阵(M T JM=J).这种特性使得Hamilton系统的解能保持相空间的辛结构.
在数值计算中，若从z k=z(tk)到z k+1的差分格式的Jacobi矩阵能满足正则变换(即差分格式的Jacobi矩阵z k+1/z k是辛矩阵)，则称此差分格式为辛差分格式.
因此，根据辛矩阵的定义可知
所以其耗散误差为
这就说明了辛算法在计算过程中不会发生耗散现象，能保持系统的辛结构，保证整体性规律的守恒.
实际中，往往更关心时域上数值解精度的高低.Brusa和Nigro［21］给出了谐振子系统
q=－ω2q
的差分格式的色散量(即phase lag，相位误差):
通过上式即可计算出算法的相位误差.也可证明，对于n阶辛差分格式，其相位精度为2n阶.
2 基于相位误差分析的精细辛算法
考虑线性Hamilton正则方程:
初始条件为
其精确解为
2.1 基于相位误差的精细辛算法格式构造
系统(2)的一阶辛差分格式为
式中τ为时间步长，qk为tk时刻q的计算值，pk为tk时刻p的计算值.
定理1 从(qk，pk)到(qk+1，pk+1)的格式变为
该格式是辛格式.并称为基于相位误差的精细辛算法.
证明 首先，将时间步长 τ划分为 2N(N为自然数)段，根据辛格式(4)在小段
在每一小段上的传递矩阵均为
因其是辛差分格式，故有
格式(3)是由格式(5)组合而成的，即
其传递矩阵为
因此
即Mτ是辛矩阵，故格式(4)为辛格式.
此格式的耗散误差为零(det(Mτ)=1).因此，此精细辛算法在长时间计算过程中能保持原有Hamilton系统的辛结构.
2.2 程序实现
传递矩阵(6)又可记为
这里
在传递矩阵Mτ的计算过程中采用精细算法:
而(E+B)(E+B)=E+2B+B·B.从而计算机上的算法实现为
for i=1:N
B=B*B+2*B
end
M=E+B.
2.3 相位误差分析
定理2 格式(5)具有2阶相位精度.
证明 根据格式(5)可以得到
再根据相位误差的公式，得到
因此，格式(5)具有2阶相位精度.□
定理3 设N为精细参数，当N趋近于+∞时，此精细辛算法的相位误差趋于零.
证明 精细辛算法的传递矩阵又可转化为
这里
若时间步长τ的划分越来越细，即N→+∞，则
因此，精细辛算法的相位精度就取决于e C中C的特征值.
根据特征方程，可以得出特征值为
即
因此，根据相位误差的公式，可得
φ=ωτ－ωτ=0.
即当时间的划分越来越细时，其相位精度(色散误差)几乎可以趋于零. □
因此，本文所设计的精细辛算法既能满足Hamilton系统辛结构的要求，又能大幅度提高辛算法的相位精度.这样，时域上解的精度就能大幅度提高，并且在长时间计算过程中结果不失真.
3 数值算例
例1 考虑简谐振子系统
的精确解:
1)相图
取时间步长τ=0.1，精细常数N=1，经过100步的计算后，各种算法的相图如图1所示.
图1 各种算法关于谐振子系统的相图
Fig.1 The phase diagram
图2 Hamilton函数的计算图
Fig.2 The Hamiltonian function diagram
从图1中可以看出，精细辛算法、三阶辛算法(SPRK3算法)的相图与精确解的相图完全吻合，说明了这些算法能保持Hamilton系统的辛结构，同时也从数值算例上证实了精细辛算法是辛算法.另外，精细算法(HPD)并不是辛算法，当精细参数比较小时，它不能保持系统的辛结构.
2)Hamilton函数
H(q，p)=q2+p2=1.
分别采用精细辛算法(时间步长τ=0.1)、精细积分法(时间步长τ=0.1)和三阶辛算法(时间步长τ=0.01)计算Hamilton函数，其结果如图2所示，纵坐标ΔH= H－H compute.从图2中可以看出，精细辛算法比精细算法至少高3个精度，比三阶辛算法高6个精度.
3)时域上解的精度
取时间步长τ=0.1，常用精细常数N=20，观察p在时域上的绝对误差(图3).从图3中可以看出，精细辛算法比三阶辛算法至少高5个精度.精细辛算法有了辛结构的保证，比一般的精细算法高9个精度，数值精度在长期的计算过程中得到了保证.
当精细常数N=40时，计算q的误差(图4)，精细辛算法的计算精度非常高，已经基本上达到了计算机的计算精度.
图3 关于p的绝对误差
Fig.3 The absolute error of p
图4 关于q的绝对误差
Fig.4 The absolute error of q
例2 对于高低混频系统，文献［22］指出数值方法的步长和系统频率之间必须满足CFL条件，即hw≈1.这意味着对于高频信号的仿真，数值上必须采用极小的步长.考虑Hamilton函数:
正则方程组为
其初值为
解析解为
数值结果分别从高频系统、低频系统、Hamilton函数、计算时间4个方面去探讨分析.
1)高频系统
从图5中可看出，精细辛算法在大步长下，对高频信号仿真的精度达到10－11，比常规精细积分法高4～5个精度，比三阶辛格式高6个精度.图6计算了q1的误差，精细辛算法的计算精度非常高.图7展示了在计算35 s后，由于相位误差的存在，三阶辛算法在时域上已经完全颠倒，图像严重失真.而精细辛算法数值解的结果与精确解基本上吻合.
图5 关于p1的绝对误差
Fig.5 The absolute error of p1
图6 关于q1的绝对误差Fig.6 The absolute error of q1
图7 q1的数值结果
Fig.7 The result of q1
2)低频系统
从图8、图9中可看出，对于低频信号，几种算法的精度都非常高，几乎达到了计算机的精度 10－14.
图8 关于q2的绝对误差
Fig.8 The absolute error of q2
图9 关于p2的绝对误差
Fig.9 The absolute error of p2
3)Hamilton函数
从表1中可以看出，随着精细参数N的不断增大，Hamilton函数的误差不断减小，尤其是精细辛算法，由于辛算法的保证，即使在步长比较大(τ=0.1)，精细参数比较小(N=10)的情况下，也能保证Hamilton函数的守恒，达到10－3.当精细参数N=40时，Hamilton函数的精度达到 10－12.
从表2可以看出，三阶辛算法在步长比较大(τ=0.1)的情况下是无法计算的，只有步长非常小(τ=0.001)的情况下，才能计算得比较好，Hamilton函数的精度只达到了10－5.这说明精细辛算法优于常规辛算法和精细算法.
表1 不同精细参数下Hamilton函数的相对误差(τ=0.1)
Table 1 The relative errors of the Hamiltonian function for different parameters(τ =0.1)
t/s 5 10 50 100 HPD(N=10) 1.168 5E+2 1.578 9E+4 1.452 8E+21 2.3279E+42 HPD(N=20) 0.004 7 0.009 5 0.048 7 0.099 9 HPD(N=40) 4.45E－9 9.00E－9 4.537E－8 9.082E－8 HPD-symplectic(N=10) 8E－4 3.5E－3 3.5E－3 2.1E－3 HPD-symplectic(N=20) 8.39E－7 3.408E－6 3.549E－6 2.407E－6 HPD-symplectic(N=40) 8.35E－13 3.179E－12 3.027E－12 3.007E－12
表2 三阶辛算法在不同步长下Hamilton函数的相对误差
Table 2 The relative errors of the Hamiltonian function for different steps by SPRK3
t/s 5 10 50 100 SPRK3(τ =0.1) inf － － －SPRK3(τ=0.01) 7.5E－3 2.41E－2 9.9E－3 2.19E－2 SPRK3(τ=0.001) 1.928E－5 2.678E－5 4.07E－6 1.135E－5
4)计算时间
从表3中可以看出，随着精细参数N的增加，计算时间的增加非常少，而随着三阶辛算法步长变小，计算时间却增加非常多.这说明精细辛算法在基本不增加工作量的前提下，能大幅度提高数值精度.
表3 各种算法计算100 s所需的CPU时间
Table 3 The CPU time for computing
method t/s HPD-symplectic(N=20)0.228 HPD-symplectic(N=40)0.264 SPRK3(τ=0.01) 1.897 SPRK3(τ=0.001) 146.924
例3 上述两个系统均为线性Hamilton系统，下面以Kepler二体问题为例阐述精细辛算法计算非线性Hamilton系统的效能:
若初始值选椭圆轨道的偏心率e=0，则此问题的精确解为
计算过程中采用局部线性化处理的方法，使得每一步长内的传递矩阵为常数矩阵，并且这一矩阵为辛矩阵.其数值结果分别从相图、时域上解的精度、精细参数的选择3个方面去探讨.
1)相图
取精细参数N=10，时间步长τ=0.1，经过1 000 s的计算以后，各种算法的相图如图10所示.
图10 关于(p2，q2)和(p1，q1)的轨道图
Fig.10 The phase diagrams of(p2，q2)and(p1，q1)
从图10中可以看出，虽然此系统是非线性的Hamilton系统，但是精细辛算法的相图与精确解的相图完全吻合，说明了此精细辛算法仍能保持Hamilton系统的辛结构.但是，精细算法并不是辛算法，不能保持系统的辛结构，且随时间的推移，误差具有放大的趋势，使得结果失真.
图11 各种算法关于(p1，q1)的绝对误差
Fig.11 The absolute errors of(p1，q1)
2)时域上解的精度
图11中，纵坐标lg(Δp1)=lg( p1－p1)，lg(Δq1)=lg( q1－ q1)，其中 p1，q1是精确解，p1，q1是数值解.从图11和表4、表5中可以看出，大步长的精细辛算法比小步长的三阶辛算法高3～5个精度，而常规的精细算法与三阶辛算法计算结果相差不大.
表4 各种算法关于p2的相对误差
Table 4 The relative errors of p2 for different parameters
t/s 50 100 500 1 000 SPRK3(τ=0.01) 3.293E－8 1.330E－7 5.479E－7 3.169E－6 HPD(N=40，τ =0.1) 9.650E－11 8.069E－10 1.799E－8 2.004E－7 HPD-symplectic(N=20，τ =0.1) 1.232E－8 2.784E－8 2.510E－8 6.971E－8 HPD-symplectic(N=40，τ =0.1) 2.474E－14 2.024E－13 5.302E－13 3.151E－11
表5 各种算法关于q2的相对误差
Table 5 The relative errors of q2 for different parameters
t/s 50 100 500 1 000 SPRK3(τ=0.01) 3.986E－7 3.697E－7 1.930E－6 1.467E－6 HPD(N=40，τ =0.1) 1.253E－9 2.332E－9 6.431E－8 9.258E－8 HPD-symplectic(N=20，τ =0.1) 3.419E－7 1.619E－7 1.797E－7 6.468E－8 HPD-symplectic(N=40，τ =0.1) 4.591E－13 6.569E－13 1.822E－12 1.453E－11
3)精细参数的选择
图 12 中，纵坐标 lg(εmax)=lg(max1≤i≤m p1，i － p1，i)，其中 p1，i是 p1在 i·τ处的精确解，p1，i是p1在i·τ处的数值解，m=100为最大步数.从图中可以看出，随着精细参数N的增大和步长τ的减小，数值精度得到了极大的提高.当精细参数N从1增加到20时，数值结果提高了6个精度，从20增加到40时，再一次提高了8个精度.N=40后，数值精度并没有得到明显的提高.另外，在相同的精细参数N下，步长增大为50倍，数值结果只减少了1个精度，说明步长对精度影响并不大，当N=40时，其精度也可达到O(10－13).因此，在数值计算过程中可选取大步长，精细参数N=40.
图12 数值精度随精细参数和时间步长的变化
Fig.12 The changes of numerical precision with parameters and time steps
3 结 论
精细辛算法继承了辛算法在求解Hamilton系统中保辛的特性，同时精细化时间步长后，在不增加原有计算量的前提下，极大地减少了相位误差，使其时域上解的精度几乎达到了计算机的精度.数值结果令人满意，说明该方法的有效性，具有广泛的工程实践价值.
致谢 本文作者衷心感谢上海第二工业大学应用数学重点学科(XXKPY1604)对本文的资助.
［1］ FENG K.Difference schemes for Hamiltonian formalism and symplectic geometry［J］.Journal of Computational Mathematics，1986，4(3):279-289.
［2］ GRTZ P.Backward error analysis of symplectic integrators for linear separable Hamiltonian systems［J］.Journal of Computational Mathematics，2002，20(5):449-460.
［3］ 邢誉峰，杨蓉.单步辛算法的相位误差分析及修正［J］.力学学报，2007，39(5):668-671.(XING Yufeng，YANG Rong.Phase errors and their correction in symplectic implicit singlestep algorithm［J］.Chinese Journal of Theoretical and Applied Mechanics，2007，39(5):668-671.(in Chinese))
［4］ MONOVASILIS T，KALOGIRATOU Z，SIMOS T E.Symplectic partitioned Runge-Kutta methods with minimal phase-lag［J］.Computer Physics Communications，2010，181(7):1251-1254.
［5］ SIMOS T E.A two-step method with vanished phase-lag and its first two derivatives for the numerical solution of the Schrdinger equation［J］.Journal of Mathematical Chemistry，2011，49(10):2486-2518.
［6］ MONOVASILIS T，KALOGIRATOU Z，SIMOS T E.Two new phase-fitted symplectic partitioned Runge-Kutta methods［J］.International Journal of Modern Physics C，2011，22(12):1343-1355.
［7］ 陈璐，王雨顺.保结构算法的相位误差分析及其修正［J］.计算数学，2014，36(3):271-290.(CHEN Lu，WANG Yushun.Phase error analysis and correction of structure preserving algorithms［J］.Mathematica Numerica Sinica，2014，36(3):271-290.(in Chinese))
［8］ VYVER H V D.A symplectic Runge-Kutta-Nystrom method with minimal phase-lag［J］.Physics Letters A，2007，367(1/2):16-24.
［9］ MONOVASILIS Th，KALOGIRATOU Z，SIMOS T E.Exponentially fitted symplectic Runge-Kutta-Nystrm methods［J］.Applied Mathematics ＆ Information Sciences，2013，7(1):81-85.
［10］ 刘晓梅，周钢，王永泓，等.辛算法的纠飘研究［J］.北京航空航天大学学报，2013，39(1):22-26.(LIU Xiaomei，ZHOU Gang，WANG Yonghong，et al.Rectifying drifts of symplectic algorithm［J］.Journal of Beijing University of Aeronautics and Astronautics，2013，39(1):22-26.(in Chinese))
［11］ MONOVASILIS Th.Symplectic partitioned Runge-Kutta methods with the phase-lag property［J］.Applied Mathematics and Computation，2012，218(18):9075-9084.
［12］ XI X P，SIMOS T E.A new high algebraic order four stages symmetric two-step method with vanished phase-lag and its first and second derivatives for the numerical solution of the Schrdinger equation and related problems［J］.Journal of Mathematical Chemistry，2016，54(7):1417-1439.
［13］ 朱帅，周钢，刘晓梅，等.精细辛有限元方法及其相位误差研究［J］.力学学报，2016，48(2):399-405.(ZHU Shuai，ZHOU Gang，LIU Xiaomei，et al.Precise symplectic time finite element method and the study of phase error［J］.Chinese Journal of Theoretical and Applied Mechanics，2016，48(2):399-405.(in Chinese))
［14］ 钟万勰.结构动力方程的精细时程积分法［J］.大连理工大学学报，1994，37(2):131-136.(ZHONG Wanxie.On precise time-integration method for structural dynamics［J］.Journal of Dalian University of Technology，1994，37(2):131-136.(in Chinese))
［15］ 钟万勰，吴峰，孙雁，等.保辛水波动力学［J］.应用数学和力学，2018，39(8):855-874.(ZHONG Wanxie，WU Feng，SUN Yan，et al.Symplectic water wave dynamics［J］.Applied Mathematics and Mechanics，2018，39(8):855-874.(in Chinese))
［16］ 富明慧，李勇息.求解病态线性方程组的预处理精细积分法［J］.应用数学和力学，2018，39(4):462-469.(FU Minghui，LI Yongxi.A preconditioned precise integration method for solving ill-conditioned linear equations［J］.Applied Mathematics and Mechanics，2018，39(4):462-469.(in Chinese))
［17］ HUANG Y A，DENG Z C，YAO L X.An improved symplectic precise integration method for analysis of the rotating rigid-flexible coupled system［J］.Journal of Sound and Vibration，2007，299(1/2):229-246.
［18］ 曾进，周钢.精细辛算法［J］.上海交通大学学报，1997，31(9):31-33.(ZENG Jin，ZHOU Gang.Precise symplectic algorithm［J］.Journal of Shanghai Jiaotong University，1997，31(9):31-33.(in Chinese))
［19］ 黄永安，尹周平，邓子辰，等.多体动力学的几何积分方法研究进展［J］.力学进展，2009，39(1):44-57.(HUANG Yongan，YIN Zhouping，DENG Zicheng，et al.Progress in geometric integration method for multibody dynamics［J］.Advances in Mechanics，2009，39(1):44-57.(in Chinese))
［20］ 徐明毅，张勇传.精细辛算法的高效格式和简化计算［J］.力学与实践，2005，27(1):55-57.(XU Mingyi，ZHANG Yongchuan.Efficient format and simple computation of precise symplectic integration method［J］.Mechanics in Engineering，2005，27(1):55-57.(in Chinese))
［21］ BRUSA L，NIGRO L.A one-step method for direct integration of structural dynamic equations［J］.International Journal for Numerical Methods in Engineering，1980，15(5):685-699.
［22］ DAVID C，ERNST H，LUBICH C.Numerical energy conservation for multi-frequency oscillatory differential equations［J］.BIT Numerical Mathematics，2005，45(2):287-305.