4.2. ハミルトン系(1自由度系)

ハミルトン方程式を数値的に解く際の手順も全く同様である. ハミルトニアンが

\[H(q,p) = T(p) + V(q)\]

の形で与えられた状況を考える.ハミルトン方程式は

\[\dot{q} = \frac{dH}{dp},\quad \dot{p} = -\frac{dH}{dq}\]

である.カオス系に移行することを前提としているので, 導入として剛体振り子の運動を考える剛体振り子

\[T(p)=p^2/2, \quad V(q)=k \cos q\]

を考えよう. \(q \in (0,2\pi]\) の周期境界条件を課して 初期条件を変化させた際の軌道を位相空間に示す.

# pendulum.py

import numpy as np
import pylab
import scipy.integrate

name = 'Pendulum'
def Pendulum(x,t=0,k=1.0):
    x0 = x[1]
    x1 = -k*np.sin(x[0])
    return np.array([x0,x1])

f=Pendulum
""" parameter setting """
k=1.0
args = (k,) # tuple

""" initial condition """
x0 = np.array([ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
x1 = np.array([-2.5, 2.5,-2.0, 2.0, 0.5, 1.0, 1.5])
initsets = [np.array(x) for x in zip(x0,x1)]

""" integration """
dt=0.01
t = np.arange(0, 20, dt)
trajectories = [scipy.integrate.odeint(f,x0, t, args, Dfun=None) for x0 in initsets ]

for i,x in enumerate(trajectories):
    x = x.T
    q = x[0] - np.floor(x[0]/(2.0*np.pi))*2.0*np.pi
    pylab.plot(q[::2],x[1][::2],'.k')    

pylab.xlabel("q",fontsize=20)
pylab.ylabel("p",fontsize=20,rotation='horizontal')
#pylab.savefig(name+'.eps')
pylab.show()

本質的ではないが,ベクトル場は

nb_points=10
xmin,xmax = pylab.xlim(0.0,2.0*np.pi)
ymin,ymax = pylab.ylim(-2.5,2.5)
x = np.linspace(xmin, xmax, nb_points)
y = np.linspace(ymin, ymax, nb_points)

X1 , Y1  = np.meshgrid(x, y)
DX1, DY1 = f([X1, Y1],1,k) # compute growth rate with f
M = (np.hypot(DX1, DY1))   # Norm of the growth rate 
M[ M == 0] = 1.            # Avoid zero division errors 
DX1 /= M                   # Normalize each arrows
DY1 /= M
pylab.quiver(X1, Y1, DX1, DY1, M, width=3e-2,pivot='mid', cmap=pylab.cm.jet)        
pylab.xlabel("q",fontsize=20)
pylab.ylabel("p",fontsize=20,rotation='horizontal')
#pylab.savefig(name+'.eps')
pylab.show()

で求めることができる.

_images/pendulum.png

4.3. 自励系(1.5自由度ハミルトン系)

さて,ポテンシャルに周期摂動を加えてみよう.ポテンシャルを

\[V(q,t)=k\cos q + \epsilon\cos(q)\sin(\omega t)\]

として揺する. 解は摂動の周期で振動するので, \(\omega t=2\pi n, (n=1,2,3,\cdots)\) を通過した際に \(q(t),p(t)\) を記録する(Poincare 断面). 計算時間が比較的かかるので注意.

# forced_pendulum.py

import numpy as np
import pylab
import scipy.integrate

twopi=2.0*np.pi
name="ForcedPendulum"
def ForcedPendulum(x,t=0,k=1.0,eps=1.0,omega=twopi):
    def f(t):
        return eps*np.sin(omega*t)        
    x0 = x[1]
    x1 = -k*np.sin(x[0])+eps*np.sin(x[0])*f(t)
    return np.array([x0,x1])

f=ForcedPendulum
""" parameter setting """
k, eps,omega = 1.0, 1.5,twopi
args = (k,eps,omega)

sample=50	
x0_0 = np.linspace(0.0,twopi,sample/2)
x1_0 = np.array([0.0]*int(sample/2))
x0_1 = np.array([np.pi]*int(sample/2))
x1_1 = np.linspace(-np.pi, np.pi,int(sample/2))
x0 = np.append(x0_0,x0_1)
x1 = np.append(x1_0,x1_1)
initsets = [np.array(x) for x in zip(x0,x1)]

dt=0.05
tmax=500
time = np.arange(0, tmax, dt)
trajectories = [scipy.integrate.odeint(f,x, time, args, Dfun=None) for x in initsets ]

T=2.0*np.pi/omega
index= [(a-T)>0 for a in time]
stb = len(time)-np.sum(index) - 1

with open("trajectories.dat","wb") as of: 
	for i,x in enumerate(trajectories):
		x = x.T
		q = x[0] - np.floor(x[0]/twopi)*twopi
		np.savetxt(of,np.array([ q[::stb], x[1][::stb]]).transpose())
		pylab.plot(q[::stb],x[1][::stb], '.',markersize=1)

pylab.xlim(0.0,twopi)
pylab.ylim(-np.pi,np.pi)
pylab.xlabel("q")
pylab.ylabel("p")
pylab.savefig("eps%.2f.png" % eps)
#pylab.show()
Poin1 Poin2 Poin3
左から順に \(\epsilon=0\)\(\epsilon=0.5\)\(\epsilon=1.5\)