ハミルトン系の数値計算 (PRML11章)

  • ????????????????????

PRML (→amazon;パターン認識と機械学習 下) の11章に出てくるハミルトン系の話。ここではハイブリッドモンテカルロのためにハミルトン系を導入して、リープフロッグスキームと呼ばれる数値積分(差分スキーム)を導入している。リープフロッグとかは昔ひと通り勉強したのだけど、久しぶりに見て懐かしかったので思い出しながら/遊びがてら計算してみた。

ハミルトン系

\(\boldsymbol p\) を一般化運動量、\(\boldsymbol q\) を一般化座標として(この2つのベクトルの次元は同じ)、\(\boldsymbol p, \boldsymbol q\) の関数 \(H(\boldsymbol p,\boldsymbol q)\) (=ハミルトニアン!!)に関する以下の微分方程式をハミルトン力学系という:
\begin{align}
\dot{\boldsymbol q}&=\frac{\partial H}{\partial \boldsymbol p}\\
\dot{\boldsymbol p}&=-\frac{\partial H}{\partial \boldsymbol q}
\end{align}
最も簡単なハミルトン系。\(p, q\) はそれぞれ1次元として、
\begin{align}
H(p,q)&=p^2+q^2
\end{align}
とする。このときハミルトン系は
\begin{align}
\dot{q}&=p\\
\dot{p}&=-q
\end{align}
となる。これは調和振動子と呼ばれる、単なる線型常微分方程式なので、初期値を \((p(0), q(0))=(1,0)\) とか置けば、解析的に解けて
\begin{align}
p(t)&=\cos t\\
q(t)&=\sin t
\end{align}
となる。つまり、単位円周上を反時計回りにぐーるぐーる回る。PRMLにはハミルトン系は

  • ハミルトニアン \(H\) の値
  • 相空間の体積

の2つの量が時間変化に対して不変である、と言っている。ハミルトニアンが不変であることは pp.265 に書いてあるとおり(ただしこのケースに限って言えば \(\sin^2 t+\cos^2 t=1\) より明らか)。ここでは相空間の体積について考える。

ハミルトン系を数値的に解く

上の例は解析解がわかっているけども、あえて数値的に解いてみる。2つのアルゴリズムを比較してみる。

オイラー法

1つ目は最も単純なオイラー法と呼ばれる離散化。時間を刻み幅 \(h\) でシンプルに離散化して、
\begin{align}
q_{t+1}&=q_t+h\cdot p_t\\
p_{t+1}&=p_t-h\cdot q_t
\end{align}
となる。\(h\to 0\) とすればハミルトン系そのものになる。

リープフロッグ法

2つ目は PRML でも登場するリープフロッグ法。これは
\begin{align}
p_{t+1/2}&=p_{t-1/2}-h\cdot q_t\\
q_{t+1}&=q_{t}+h\cdot p_{t+1/2}
\end{align}
と表される(この表記はPRMLと異なるけれども本質的には同じ)。\(1/2\) ステップというわかりづらい表記があるけど、これは刻み幅の半分の \(h/2\) だけ進んだ場所のイメージ。初期値 \(p_0\) が与えられた時、\(p_{0+1/2}\) は
\begin{align}
p_{0+1/2}&=p_{0}-\frac{h}2\cdot q_0
\end{align}
で与えることにする。イメージ的には次のような感じ(赤いのが運動量ガエルで緑が座標ガエル)

結果

数値解の相空間での挙動を見てみる

実際にオイラー法とリープフロッグ法で解いてみて、相空間上にプロットした図。

  • 解析解は図の t=0 から始まって単位円周を反時計回りに回る。
  • オイラー法は渦巻状に発散してしまう。
  • 一方、リープフロッグ法はほぼ真の解を再現している。

数値解のハミルトニアンの挙動を見てみる

時間発展に関して保存するはずの \(H(p,q)\) の挙動を数値的に計算してプロットしてみる。

  • 上で述べたように、解析解では変化しない。
  • オイラー法では指数的に発散していく。つまり、保存則が成立していない。
  • リープフロッグ法も保存しているように見える。
  • しかし、下に示すように、詳細に見てみると厳密な保存ではなくて、振動している

数値解の小さな領域の時間変化を見てみる

最後に、PRML pp.265~266 の記述

ハミルトン力学系の2つ目の重要な性質は,位相空間の体積が保存されるものであり,リュービルの定理として知られる.(中略)実際,リュービルの定理が依然として厳密に成り立つように積分方式を工夫することができる.(中略)これを実現する方式の一つはリープフロッグ離散化と呼ばれる

について、つまり、リープフロッグの体積保存性について。

これを数値的に確かめるために、初期値近辺に小さな三角形をつくり、その各頂点がどのように時間発展するかをプロットしてみる。水色の▲はオイラー法で、紫色の▲はリープフロッグ。

上図が示すように、オイラー法では三角形の面積は次第に大きくなっていくが、リープフロッグでは面積が保存されていることがわかる。

独断と偏見によるまとめ

調和振動子と呼ばれる単純なハミルトン系の数値実験をした。

  • オイラー法は単純だけど、かなり悪い。
  • 一方、リープフロッグは陽的スキーム(ステップ毎に連立方程式などを解く必要がない)であるにもかかわらずかなりよさげな性能を示した。
  • (以下、昔の記憶をたどりながら書きます)リープフロッグ公式はテイラー展開した時の近似オーダーは低いけど、体積保存則やハミルトニアンの保存則が成立する。前者は実はある意味で厳密で、後者は厳密ではない(実は別の量=離散版ハミルトニアンのようなもの、が保存している)
  • 陽的スキームにもかかわらずこの手の性能を示すような手法はあまりない。
  • 有名なルンゲ-クッタ法に関して言えば、体積保存性や離散的なハミルトニアンがないため、テイラー展開の意味での近似精度はいいにもかかわらず、長時間の計算では失敗するはず

ソースコード

こういうのはホントは Matlab か Python でやったほうがいい気もするけど、R。

# Hamiltonian
H <- function(pq){
  sum(pq^2)
}

# Euler method
step.euler <- function(p,q,h){
  c(p-h*q,q+h*p)
}

# Leapfrog method
step.leapfrog <- function(p,q,h){
  p.half <- p-(h/2)*q
  q.1 <- q+h*p.half
  p.1 <- p.half-(h/2)*q.1
  c(p.1,q.1)
}

# True solution
true.solution <- function(h,N,simulate=T){
  if(simulate)
    return(t(sapply(1:N,function(n) c(cos(h*n),sin(h*n)))))
  tt=seq(0,2*pi,2*pi/1000)
  t(sapply(tt,function(x) c(cos(x),sin(x))))
}



# Setting
h <- 0.05
N <- 1001 # Total number of the steps

# Solving hamilton system by Euler method
pq <- c(1,0)
pq.t <- matrix(0,nc=2,nr=N)
pq.t[1,] <- pq

for(tt in 2:N){
  pq <- step.euler(pq[1],pq[2],h)
  pq.t[tt,] <- pq
}

# Solving hamilton system by leapfrog method
pq <- c(1,0)
pq.t.lf <- matrix(0,nc=2,nr=N)
pq.t.lf[1,] <- pq

for(tt in 2:N){
  pq <- step.leapfrog(pq[1],pq[2],h)
  pq.t.lf[tt,] <- pq
}

# Plotting orbit
plot(pq.t,type="n",xlab="p(t)",ylab="q(t)",xlim=c(-4,4),ylim=c(-4,4))
lines(pq.t,lw=3,col="steelblue")
lines(pq.t.lf,type="l",col=3,lw=3)
lines(true.solution(h,N,F),col=2)
title("Phase space plot of numerical solution of Hamilton system")
legend("topleft",c("Euler","Leapfrog","Analytic"),lt=1,col=c("steelblue",3,2),lw=c(2,2,1))
text(pq.t[c(1,N),],c("t=0","t=50"),adj=1.2)
points(pq.t[c(1,N),],pch=20)
dev2bitmap(file="evo_euler.jpg",taa=4,gaa=4)

# Plotting Hamiltonian
bitmap(file="H_evo_euler.jpg",taa=4,gaa=4,height=5,width=7)
plot(h*(0:(N-1)),apply(pq.t,1,H),xlab="t",ylab="H(p,q)",type="n",ylim=c(0,13))
lines(h*(0:(N-1)),apply(pq.t,1,H),col="steelblue",lw=5)
lines(h*(0:(N-1)),apply(pq.t.lf,1,H),xlab="t",ylab="H(p,q)",type="l",ylim=c(0,13),col=3,lw=5)
lines(h*(0:(N-1)),apply(true.solution(h,N),1,H),col=2)
title("Time evolution of H(p,q)")
legend("topleft",c("Euler","Leapfrog","Analytic"),lt=1,col=c("steelblue",3,2),lw=c(5,5,1))
dev.off()

# Plotting Hamiltonian (2) -- only leapfrog
bitmap(file="H_evo_leapfrog.jpg",taa=4,gaa=4,height=5,width=7)
plot(h*(0:(N-1)),apply(pq.t.lf,1,H),xlab="t",ylab="H(p,q)",type="n",ylim=c(0.99,1.01))
lines(h*(0:(N-1)),apply(pq.t.lf,1,H),xlab="t",ylab="H(p,q)",type="l",ylim=c(0,13),col=3,lw=5)
lines(h*(0:(N-1)),apply(true.solution(h,N),1,H),col=2)
title("Time evolution of H(p,q) (Leapfrog)")
legend("topleft",c("Leapfrog","Analytic"),lt=1,col=c(3,2),lw=c(5,1))
dev.off()

# Small triangle evolution
pq <- c(1.2,0)  # Modifying initial condition!!
pq.lf <- pq

pq.t.2 <- matrix(0,nc=2,nr=N)
pq.t.2[1,] <- pq

pq.t.2.lf <- matrix(0,nc=2,nr=N)
pq.t.2.lf[1,] <- pq

for(tt in 2:N){
  pq <- step.euler(pq[1],pq[2],h)
  pq.lf <- step.leapfrog(pq.lf[1],pq.lf[2],h)
  pq.t.2[tt,] <- pq
  pq.t.2.lf[tt,] <- pq.lf
}

pq <- c(1.1,0.2)  # Modifying initial condition!!
pq.lf <- pq

pq.t.3 <- matrix(0,nc=2,nr=N)
pq.t.3[1,] <- pq

pq.t.3.lf <- matrix(0,nc=2,nr=N)
pq.t.3.lf[1,] <- pq

for(tt in 2:N){
  pq <- step.euler(pq[1],pq[2],h)
  pq.lf <- step.leapfrog(pq.lf[1],pq.lf[2],h)
  pq.t.3[tt,] <- pq
  pq.t.3.lf[tt,] <- pq.lf
}

# Plotting invariance of area
nn <- seq(1,N,200)
plot(pq.t,type="n",xlab="p(t)",ylab="q(t)",xlim=c(-4,4),ylim=c(-4,4))
lines(pq.t,lw=3,col="steelblue")
lines(pq.t.lf,type="l",col=3,lw=3)
lines(true.solution(h,N,F),col=2)
title("Time evolution of small polygon")
legend("topleft",c("Euler","Leapfrog","Analytic"),lt=1,col=c("steelblue",3,2),lw=c(3,3,1))

triangles.lf <- triangles.euler <- cbind(pq.t.lf,pq.t.2.lf,pq.t.3.lf)[nn,]
apply(triangles.lf,1,function(row) polygon(matrix(row,nc=2,byrow=T),border=NA,col="purple"))
points(pq.t.lf[nn,],pch=20)
text(pq.t.lf[nn,],paste("t=",(nn-1)*h),adj=1.4,col="purple")

triangles.euler <- cbind(pq.t,pq.t.2,pq.t.3)[nn,]
apply(triangles.euler,1,function(row) polygon(matrix(row,nc=2,byrow=T),border=NA,col="turquoise"))
points(pq.t[nn,],pch=20)
text(pq.t[nn,],paste("t=",(nn-1)*h),adj=1.4)

dev2bitmap(file="poly_evo_euler.jpg",taa=4,gaa=4)
はてなブックマーク - ハミルトン系の数値計算 (PRML11章)
Pocket