0

我正在尝试使用三阶辛积分器 Ruth 绘制 Henon-Heiles 哈密顿量的庞加莱曲面部分，但我遇到了一些问题，因为经过几个步骤后，变量开始取巨大值导致错误。我不明白为什么会发生这种情况，因为我使用的积分器应该是辛的。

我曾多次尝试编写此代码，但我总是遇到某种错误，所以如果有人能解释我做错了什么，我将不胜感激。

我用python写的代码是这样的：

import math as m
import numpy as np
from matplotlib import pyplot as plt

def xpunto(vx0): #time derivative of x
    return vx0
def pxpunto(x0,y0,lambda0): #time derivative of px
    return -x0-(2.*lambda0*x0*y0)
def ypunto(vy0): #time derivative of y
    return vy0
def pypunto(x0,y0,lambda0): #time derivative of py
    return -y0-(lambda0*((x0*x0)-(y0*y0)))
def energia_poinc(y0,py0, H0): #initial condition for px using the constraints
    return m.sqrt((2.*H0)-(py0*py0)-(y0*y0)+((2./3.)*(y0*y0*y0)))


iteraciones=100000
dt=0.01
t=[0.]
x=[]
y=[]
px=[]
py=[]

def ruth3rd(xx,yy,ppx,ppy,lambda00):#define function with symplectic integrator

    aa=[1,-1./24.,3./4.,7./24.] #constants for position variable of the integrator
    bb=[1,1.,-2./3., 2./3.] #constants of velocity of the integrator

    x.append(xx)
    y.append(yy)
    px.append(ppx)
    py.append(ppy)



    xi , yi , pxi , pyi= xx, yy, ppx , ppy
    for i in range (iteraciones):

        vvx = pxi+bb[1]*pxpunto(xi,yi,lambda00)*dt
        VV1y = pyi+bb[1]*pypunto(xi,yi,lambda00)*dt
        xx = xi + aa[1]*vvx*dt
        yy =yi + aa[1]*VV1y*dt
        for j in range(2,len(aa)):
            vvx += bb[j]*pxpunto(xx,yy,lambda00)*dt
        VV1y += bb[j]*pxpunto(xx,yy,lambda00)*dt
            xx += aa[j]*vvx*dt
            yy += aa[j]*VV1y*dt

        xi = xx
        pxi = vvx
        yi = yy
        pyi = VV1y

        x.append(xi)
        px.append(pxi) 
        y.append(yy)
        py.append(pyi) 
        t.append(t[-1] + dt)

setsrange=10
sets=np.array(np.linspace(-0.4, 0.4, 10),dtype=np.float64)


for I in range(1):

    xpoinc=[]
    ppoinc=[]

    for ii in range (1):
        for jj in range (setsrange):

            Hinicial=1./12.
            x=[0.] 
            y=[0.4]
            py=[sets[jj]] 
            print (py[0],y[0],x[0],Hinicial)
            px=[energia_poinc(y[0],py[0],Hinicial)] 

            ruth3rd(x[0],y[0],px[0],py[0],2.)       


        crossings=[]
            crossingsindex=[]    
            for idx, item in enumerate(x[:-1]): #intersections
                if x[idx] < 0 and x[idx+1] > 0:
                    crossings.append(x[idx])
                    crossingsindex.append(idx)


            for index in crossingsindex:
                xpoinc.append(y[index])
                ppoinc.append(py[index])
4

0 回答 0