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

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

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

f:id:StatModeling:20201106184256p:plain

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

1次元の場合

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

f:id:StatModeling:20201106184300p:plain

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

f:id:StatModeling:20201106184303p:plain

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

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

f:id:StatModeling:20201106184307p:plain

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

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

f:id:StatModeling:20201106184311p:plain

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

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

f:id:StatModeling:20201106184314p:plain

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

f:id:StatModeling:20201106184244p:plain

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

2次元円対称の場合

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

f:id:StatModeling:20201106184318p:plain

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

f:id:StatModeling:20201106184322p:plain

f:id:StatModeling:20201106184326p:plain

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

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

f:id:StatModeling:20201106184249p:plain

3次元球対称の場合

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

f:id:StatModeling:20201106184330p:plain

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

f:id:StatModeling:20201106184334p:plain

f:id:StatModeling:20201106184338p:plain

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

f:id:StatModeling:20201106184341p:plain

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

f:id:StatModeling:20201106184345p:plain

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

f:id:StatModeling:20201106184252p:plain

参考文献

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

Random Walks in Biology

作者:Berg, Howard C.
発売日: 1993/09/27
メディア: ペーパーバック

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