12-16.rst

問題12.16 - 拡散が支配する1次元化学反応系のモンテカルロ・シミュレーション

シミュレーションの目的

単一の種類Aの粒子からなる系を考える、全ての粒子は拡散し、2個の粒子が衝突すると反応が起こり、1個の粒子が消えるか2個の粒子が結合して、もはや反応には関わらない不活性な種類の粒子になる。 後者の場合の化学反応は

A + A → 0   (不活性)

と表すことができる。A粒子の密度の空間的ゆらぎを無視すると、その時間変化は簡単な反応速度方程式

$$\frac{dA(t)}{dt} = -kA^{2}(t)$$

で与えることができる。ここでAはA粒子の時刻tの密度であり、kは反応速度定数である。簡単のために、全ての反応物質がt = 0に加えられ、その後には反応物質は何も加えられないとする(閉じた系)。1階微分方程式12-16-e2の解は

$$A(t) = \frac{1}{kt + 1/A(0)}$$

であり、長時間の極限では

A(t) ∼ t^− 1

となる。

上に述べた1種類の粒子が消滅する過程におけるAの時間依存性は簡単であるが、空間ゆらぎが無視されている。ここでは、この過程の時間発展のシミュレーションを行い、仮定が正しいか調べることにする。

作成したプログラム

本シミュレーションで作成したプログラムを以下に示す。

1次元格子上で反応して不活性化する粒子の拡散

• 12-16_diffusion-dominant_chemical_reaction.py(download <12-16_diffusion-dominant_chemical_reaction.py>)

12-16_diffusion-dominant_chemical_reaction.py

このプログラムは3つの関数からなっている。関数main_rw_d1は、1次元格子上を複数の粒子がランダムウォークし、同じ格子点上に2個の粒子が来たときその粒子を取り除くような挙動を実現する。関数plot_graphは、問題12−12 <12-12-label>でも用いたグラフ描画のための関数である。関数fittingは、1/A(t) − 1のtに関する両対数グラフを直線でフィッティングしたときのその傾きを求めるものである。何回か試行を繰り返し、その平均値と偏差を返すようにしてある。 関数main_rw_d1の動きについての解説を加えると、まずlatticeは行方向が粒子を表し、列方向は座標を表している。単位行列によってA(t = 0) = 1を表現し、時間発展ごとに右向きか左向きに等確率でシフトする。その後、列方向に和を計算し、その値が2であるときはその列の要素をすべて0にしている(2以上の値となることはない)。この操作を繰り返し、時刻t( = 2^p(p = 1, 2, 3, ⋯, pmax))のときの1/A(t) − 1の値をAに追加している。ここでtは対数プロットしたときのデータが等間隔になるように、指数関数的にとった。

実習課題

• N個の粒子が、長さLの1次元格子上を周期的境界条件の下でランダムウォークを行う場合を考え、全ての格子点が占有されている状態A(t = 0) = 1からシミュレーションを始めよ。量1/A(t) − 1を時間tに対して両対数でプロットせよ。対数プロットのデータが等間隔になるように、時間の区間を指数関数的にとれ。長時間の極限では、得られた両対数プロットは直線的か、直線ならばその傾きを求めよ。平均場近似は1次元で正しいといえるか。L = 100程度の小さい格子と、t = 10²程度の時間で大雑把な結果を得ることができる。10パーセント以内の精度で結果を得るためには、L = 10⁴程度の大きさの格子とt = 2¹³程度の時間が必要である。

上記のプログラムを用いて、N個の粒子が1次元周期境界条件のもとでランダムウォークを行い、反応によって不活性化する問題をシミュレーションした。このとき、条件として、格子の大きさL = 1000で、時刻t = 512( = 2⁹)までの計算を行った。A(t)を時刻tでの粒子の密度として、量1/A(t) − 1を時間tに対して両対数プロットし、これを図#fig-12-16-f1に示した。図から分かるように、この両対数プロットは直線的であり、A(t)とtの間にベキ乗則が成り立つことが分かる。実際にグラフにおける直線の傾きaを算出してみると、10回の試行の平均値としてa = 0.512988987355が得られた。このときの分散σはσ = 0.0521364760422であった。したがって、時刻tとその時刻での粒子の存在密度A(t)の間には、およそ以下のような関係が成り立つと分かる:

$$A(t) \sim t^{-\frac{1}{2}}$$

この結果は、空間ゆらぎを考慮せず解析的に求めた、長時間後のA(t)とtの関係12-16-e3と異なっている。したがって、式12-16-e2のように平均場近似をすることは妥当でなく、粒子密度の空間的なゆらぎの効果は無視できないことが分かった。

1/A(t) − 1の時間tに対する両対数プロット

まとめ

1次元の化学反応系の簡単な例について考え、粒子の存在密度A(t)と時刻tの間にベキ乗則が成り立ち、その指数の値が平均場近似より求められた値と異なることを示した。

参考文献

• ハーベイ・ゴールド,ジャン・トボチニク,石川正勝・宮島佐介訳『計算物理学入門』,ピアソン・エデュケーション, 2000.