4.4. 離散力学系(ハミルトン系)

連続時間の時間発展は計算時間のコストがかかるため,次のハミルトニアン

\[H(q,p,t) = T(p) + V(q)\sum_n \delta(t-n)\]

の時間に周期的な劇力が加わるkicked rotatorと呼ばれるモデルがよく用いられる. kickから次のkickまでの運動方程式の解は求める(積分を実行する)ことでき, その解離散化された

\[\begin{split}p_{n+1} &= p_n - V'(q_n) \\ q_{n+1} &= q_n + T'(p_{n+1})\end{split}\]

漸化式で与えられる. ここで,プライムは引数に関する微分を意味する.

\[T(p) = \frac{p^2}{2},\qquad V(q) = k\cos(q)\]

の時、時間発展方程式はChirikov標準写像 [1] と呼ばれる.

[1]http://www.scholarpedia.org/article/Chirikov_standard_map
# map.py

import numpy as np
import matplotlib.pyplot as plt
twopi = 2.0*np.pi

def evolv(q,p,k):
	pp = p - k*func0(q)
	qq = q + func1(pp)
	return qq, pp
	
def func0(x):
	return -np.sin(x)
	
def func1(x):
	return x

sample = 100	
iter = 1000

q = np.random.random(sample)*twopi
p = np.random.random(sample)*twopi

trajectory = [np.array([]),np.array([])]
k=0.971635
for i in range(iter):
	q,p = evolv(q,p,k)
	q = q - np.floor(q/twopi)*twopi
	p = p - np.floor(p/twopi)*twopi
	trajectory[0] = np.append(trajectory[0],q)
	trajectory[1] = np.append(trajectory[1],p)	

plt.plot(trajectory[0],trajectory[1], ',')
plt.xlim(0.0,twopi)
plt.ylim(0.0,twopi)
plt.xlabel("q")
plt.ylabel("p")
plt.savefig("standardmap_k%.2f.png" % k)
plt.show()
Map1 Map2 Map3
左から順に \(k=0.5\)\(k=0.971635\)\(k=5.0\)