您需要正确计算哈密顿量的导数. |y-x|^n对x的导数是

You need to compute the derivatives of the Hamiltonian correctly. The derivative of |y-x|^n for x is

n*(x-y)*|x-y|^(n-2)=n*sign(x-y)*|x-y|^(n-1)

和y的导数几乎相同，但不完全相同(如您的代码中一样)

and the derivative for y is almost, but not exactly (as in your code), the same,

n*(y-x)*|x-y|^(n-2)=n*sign(y-x)*|x-y|^(n-1)，

请注意符号差异.通过这种校正，您可以采取更大的时间步长，并且采用正确的线性插值甚至可能更大的步长来获取图像

note the sign difference. With this correction you can take larger time steps, with correct linear interpolation probably even larger ones, to obtain the images

我将ODE的集成更改为

I changed the integration of the ODE to

t = linspace(0, 1000.0, 2000+1)
...
E_kin = E-total_energy([x0,y0,0,0])
init_cons = [[x0, y0, (2*E_kin-py**2)**0.5, py] for py in np.linspace(-10,10,8)]

outs = [ odeint(pqdot, con, t, atol=1e-9, rtol=1e-8) ) for con in init_cons[:8] ]

显然，初始条件的数量和参数化可能会发生变化.

Obviously the number and parametrization of initial conditions may change.

过零的计算和显示更改为

The computation and display of the zero-crossings was changed to

def refine_crossing(a,b):
    tf = -a[0]/a[2]
    while abs(b[0])>1e-6:
        b = odeint(pqdot, a, [0,tf], atol=1e-8, rtol=1e-6)[-1];
        # Newton step using that b[0]=x(tf) and b[2]=x'(tf)
        tf -= b[0]/b[2]
    return [ b[1], b[3] ]

# Plot Poincare sections at x=0 and px>0
fig2 = figure(2)
for ii in xrange(8):
    #subplot(4, 2, ii+1)
    xcrossings = findcrossings(outs[ii][:,0], outs[ii][:,3])
    ycrossings = [ refine_crossing(outs[ii][cross], outs[ii][cross+1]) for cross in xcrossings]
    yints, pyints = array(ycrossings).T
    plot(yints, pyints,'.')
    ylabel("py")
    xlabel("y")
    title("Poincare section x = 0")

并评估更长的积分间隔的结果

and evaluating the result of a longer integration interval