読者です 読者をやめる 読者になる 読者になる

StatModeling Memorandum

StanとRでベイズ統計モデリングします. たまに書評.

拡散方程式の解と粒子数の時間変化

理論

拡散現象は幅広く観測されます。この記事では拡散方程式の解と粒子数の時間変化のグラフを備忘録として残します。

拡散方程式は以下です。導出方法はWebにあふれているので検索してください。

image001.png

ここで t は経過時間、 p は位置を表す座標と経過時間で決まる確率分布、 D は拡散係数です。

1次元の場合

ラプラシアン x の2階微分になるので以下になります。

image002.png

初期値を \delta 関数(≒位置 x = 0 にちょぴっとある溶液をスポイトで入れた)とすると、解は以下になります。

image003.png

標準偏差 \sqrt{2 D t} 正規分布になります。経過時間 t のルートでしか幅が広がりません。ここで、 t = 0 に溶液を入れた瞬間に x = 10 などにおいて密度関数の値が非ゼロになってしまっている(粒子が瞬間移動している)のは、拡散係数の定義(仮定)に由来します。モデルは近似ですのでこういうことあります。

 t = 0  n 個(単位はmolとか)の粒子をスポイトで入れたとすると、位置 x , 経過時間 t における粒子数は、分布に n を掛けて以下になります(★式とします)。

image004.png

次に、 t = 0 から t = T までスポイトでぶちゅーと入れることを考えます。微小時間あたり n 個の粒子が入るとします。この場合に位置 x における合計の粒子数の時間変化 N(x,t) はどうなるでしょうか?

まず t \lt T の場合を考えます。 t 時間経つと粒子が拡散して位置 x において(★)式となるので、合計の粒子数は「入れたてホヤホヤ( t = 0 )で x に来た粒子」と「ちょっと前に入れて x に来た粒子」と、…、「はじめにぶちゅーと入れて t 時間経って x に来た粒子」の和(積分)になります。数式で表すと以下の積分になります(WolframAlpha使いました)。

image005.png

ここで \Phi は標準正規分布の累積分布関数です。

次に t \gt T の場合を考えます。この場合は t - T 時間前にスポイトを抜いてますので、入れたてホヤホヤなどは積分の範囲に入らず、「スポイト抜く直前に入れて t - T 時間経って x に来た粒子」と、…、「はじめにぶちゅーと入れて t 時間経って x に来た粒子」の和(積分)になります。数式で表すと以下になります。

image006.png

仮に n = 1, x = 1, T = 1 と固定し、 D の値をいろいろ変えて、 t \lt T  t \gt T の場合をあわせた関数をプロットすると以下になります。

fig1D.png

こういう時系列を僕は幾度となく見てきた覚えがあります。今度から背後にあるメカニズムとして拡散現象も考えてみようと思います。同様に2次元円対称と3次元球対称の場合の拡散方程式とその解も求めておきます。

2次元円対称の場合

拡散方程式は2次元極座標のラプラシアンを知っていれば簡単で以下になります。

image007.png

初期値を \delta 関数とすると、解は以下になります。

image008.png

image009.png

残念ながら t に関する積分が指数積分exponential integral)となっていて解析的に N(r,t) を表すことができません。そこで数値積分します。

1次元の場合と同じように、 n = 1, r = 1, T = 1 と固定し、 D の値をいろいろ変えてプロットすると以下になります。

fig2D.png

3次元球対称の場合

拡散方程式は3次元極座標のラプラシアンを知っていれば簡単で以下になります。

image010.png

初期値を \delta 関数とすると、解は以下になります。

image011.png

image012.png

 N(r,t)  t \lt T の場合に以下になり、

image013.png

 t \gt T の場合に以下になります。

image014.png

1次元の場合と同じように、 n = 1, r = 1, T = 1 と固定し、 D の値をいろいろ変えてプロットすると以下になります。

fig3D.png

参考文献

今回参考にした書籍は以下です。特に1~3章がオススメです。

Random Walks in Biology

Random Walks in Biology

2章に拡散方程式の導出と3次元球対称の場合の解が載っています。その他にも、例えば3章では3次元球対称の場合に r = b からランダムウォークをスタートして、 r = a \ (a \lt b) にたどり着くか、 r = c \ (c \gt b) にたどり着くかを計算しています。到達時間の期待値や到達確率の期待値などを簡潔に導いています。石島先生の補足ページ『「生物学におけるランダムウォーク」を理解する』で補いながら読むとよいでしょう。