はじめに
数値シミュレーションをする時、平衡状態や定常状態を調べることがよくあります。この時、十分長いシミュレーションステップ後に定常状態になると期待されますが、その途中は、初期状態から定常状態への緩和過程が観測されます。多くの場合、この緩和過程は指数関数でよくフィッティングできます。
例えば以下は二次元強磁性イジング模型における自発磁化の緩和です。L=64の正方格子で、メトロポリス法を用いて時間発展させています。温度を臨界点より少し低温にしているため、自発磁化は有限の値に収束しています。シンボルが自発磁化、実線が指数関数です。よく一致していることがわかります。
また、以下は分子動力学法における温度の緩和過程です。32000個のLennard-Jones相互作用する原子を面心立方格子状に配置し、周期境界条件、密度0.5で、Langevin熱浴にて目的温度0.8として制御しています。こちらもシンボルが系の温度、実線が指数関数で、よくフィッティングできていることがわかります。
このように、数値シミュレーションでは指数関数緩和がよく現れます。これが何故かを簡単に説明するのが本稿の目的です。
緩和過程
なにか時間に依存する物理量A(t)を考えましょう。先ほどの例でいえば、スピン系の秩序変数や、分子動力学法における系の温度などの物理量に対応します。
十分に時間が経過した時、系が平衡状態に達するとします。すると、物理量も平衡状態における値A_\inftyになるでしょう。
\lim_{t \rightarrow \infty} A(t) = A_\infty
さて、系がまだ平衡状態に達していない時、物理量Aが平衡状態の値よりも大きいなら小さく、小さいなら大きくするような「力」が働くと考えられます。この抽象的な力をFとすると、こんな微分方程式を考えることができます。
さて、平衡状態においては、物理量Aは時間変化しませんから、A(t) = A_\inftyなら力はゼロです。従って、
であるはずです。いま、物理量Aが平衡状態での値に十分近い場合、
A(t) = A_\infty + \varepsilon
と書けるでしょう。この時、
\begin{aligned}
\frac{dA}{dt} &= F(A) \\
&= F(A_\infty + \varepsilon) \\
&= F(A_\infty) + F'(A_\infty)\varepsilon + O(\varepsilon^2)\\
&= F'(A_\infty)\varepsilon + (\varepsilon^2)
\end{aligned}
さて、Fは物理量AがA_\inftyよりも大きければ小さくしようとするので負になり、小さければ大きくしようとするから正になるため、F'(A_\infty) < 0です。
そこで、
とすると、先ほどの式は
\frac{dA}{dt} = \gamma(A_\infty - A(t))
と書けます。これは、F(A)=0となる点の周りで一次までのテイラー展開をしたことになっています。この微分方程式は簡単に解けて、A(t)の時間依存性は以下のように表すことができます。
A(t) = \left(A(0) - A_\infty \right) \exp(-\gamma t) + A_\infty
さて、ここで
という量を定義します。すると、先程の式は
A(t) = \left(A(0) - A_\infty \right) \exp(-t/\tau) + A_\infty
と書けます。この\tauは緩和時間と呼ばれ、系が初期状態から定常状態にどのくらいの時間をかけて緩和するかを表す時定数となります。
以上から、定常状態を持つような系の時間発展では、定常状態に十分近ければ、指数関数的な緩和を持つことが期待されます。
緩和時間について
数値シミュレーションにおいて、最終的に興味のあるのは平衡状態や定常状態での値ですが、その値を調べる前に緩和時間を調べるのは重要です。先ほどの式をもう一度見てみましょう。
A(t) = \left(A(0) - A_\infty \right) \exp(-t/\tau) + A_\infty
t=0を代入するとA(t) = A(0)に、t=\inftyを代入するとA(t) = A_\inftyになっていることがわかります。これは、初期状態A(0)を、どのくらいの速度で忘れて定常状態での値であるA_\inftyに近づいていくかを表しています。\tauは、初期状態の影響が1/e倍になるのに必要な時間です。e \sim 2.72ですから、例えば緩和時間の5倍、すなわちt=5\tauだけのシミュレーションを行うと、初期条件の影響はe^{-5}倍、つまり0.7\%程度となり、初期条件の影響が1\%より小さくなります。もう少し余裕を見て、緩和時間の10倍のシミュレーションをすれば初期状態の影響は2万分の1程度になるため、ほぼ間違いなく緩和していると言えます。
モンテカルロ法の例を再掲します。
こちらは、緩和時間が概ね50です。なので、250ステップでほぼ緩和しており、500ステップとれば十分であることがわかります。
こちらは分子動力学法の例です。
この系では、緩和時間が100程度です。従って、500まで計算すれば緩和しており、1000まで計算すれば安全でしょう、という感じです。
まとめ
数値シミュレーションでよく現れる指数緩和についてまとめました。我々が知りたいのは定常状態での値ですが、「本当に定常状態に達しているか」を判断するのに緩和時間の情報が必要になります。系は定常状態に指数関数的に緩和し、その時の緩和の時定数を緩和時間と呼びます。だいたいの目安として、緩和時間の5倍から10倍の長さをシミュレーションすれば安全です。
ただし、臨界点の近傍であったり、物理量が振動したりする場合は指数関数的でない緩和が現れる場合もあります。その時にはグラフを見て何が起きているかを注意深く調べる必要があります。
Discussion