我编写了一个测试程序，在上向前积分，然后从哈密顿系统的正向积分端点向后积分 和的表达式（在下面的程序中给出并不重要）。$[0,T_f]$$[T_f,0]$

$$\dot q(t) = \frac{\partial H}{\partial p}(q(t),p(t)), \qquad \dot p(t) = -\frac{\partial H}{\partial q}(q(t),p(t))$$ $H$

我使用SageMath象征性地计算动力学的哈密顿量和雅可比的导数。我使用scipy python 库的“dopri5”例程进行数值积分。

该程序做了我们想要的，但是有一些数字错误，因为我们没有得到我们开始的确切初始值，我认为这些错误来自积分器？我们看到中间点 (zmid) 的顺序是 +1e12，所以它会导致错误，但其他原因可能是什么？

为了更加准确和准确，我修改了以下设置：

• 减少“dopri5”例程的atol、rtol和增加nsteps

您对提高方法的精密度和准确度有更多建议吗？从数学的角度来看，这些错误的原因是什么？例如，如果我们在动力学中有不同的时间尺度......

from scipy.integrate import ode
from sage.all import *

#Dimension of the state variables
N=2
#Final time
Tf = 10.
#Symbolic variables of the Hamiltonian
var('q1 q2 p1 p2 v')
qs = [q1,q2]
ps = [p1,p2]
zs = qs + ps

#Hamiltonian function and jacobian of the dynamics
def hfun():
    H = p1*(cos(q1^2)+1/3*sin(q2)^2) + p2*(cos(q1)-2*sin(q2*q1)^2)
    dHdp = jacobian(H,ps)
    dHdq = jacobian(H,qs)
    matzdot=block_matrix(SR,[[dHdp.T],[-dHdq.T]])
    jaczdot = jacobian(matzdot,zs)
    return H, jaczdot

#Hamiltonian dynamics
def dynz(t,z):
    H,aux = hfun()
    jacHp = jacobian(H,ps)
    jacHp = list(jacHp[0].subs({zs[i]:z[i] for i in range(0,2*N)}))
    dzdt = jacHp
    jacHq = jacobian(H,qs)
    dzdt.extend(list(-jacHq[0].subs({zs[i]:z[i] for i in range(0,2*N)})))
    return dzdt

#integration
def main():
    z0 = [1.,-0.5,-2.,1.]
    print("z0: ")
    print(z0)

    solverz = ode(dynz).set_integrator('dopri5')
    solverz.set_initial_value(z0,0.)
    solverz.integrate(Tf)

    zf = solverz.y
    print("zmid: ")
    print(zf)

    solverz.set_initial_value(zf,Tf)
    solverz.integrate(0.)
    print("z0: ")
    print(solverz.y)

    return solverz.y


main()

这是输出：

z0: 
[1.00000000000000, -0.500000000000000, -2.00000000000000, 1.00000000000000]
zmid: 
[ 1.35249010e+00 -2.07489135e+00 -3.88286701e+12  8.51614447e+11]
z0: 
[ 0.99956526 -0.50018124 -2.00172766  1.00161313]