直感的に理解する拡散方程式 - プログラムdeタマゴ

　研究室の友人が「反応拡散系」というのがわからない、とか言っていました。反応拡散系ないし反応拡散方程式というのはその名前の通り反応方程式+拡散方程式というだけだったりします。んじゃ、反応方程式って何よ、拡散方程式って何よ、ということになりますね。今日は拡散方程式について見ていきます。

　さて、「拡散」方程式というぐらいなのですから、「拡散」を表します。すなわち、なんか広がっていく様を表します。
　本当はフィックの法則とか、流体力学の分野から導くんだけど、そんな難しい知識は使いません。ここでは、最も単純なモデルを考えてみます。
　下の図の様に、ある位置にたくさんのボールがある状態を考えます。この状態は不自然なので、いつか勝手に山が崩れてボールが広がっていきますね。
　ボールが崩れて広がっていく様子は刻一刻と変わっていくのですが、ここでは非常に単純に考えて、ある時間たった後に玉は一定の割合に従って隣の位置に移動するということにしましょう。

　ここで、時刻tにおいて位置x上にある玉の数をf(x,t)と表記します。ある時間Δtだけ時間がたった後、nの割合だけ隣の位置に移動するとします。隣の位置との間隔はΔxとしておきましょう。（ここで、nは0以上0.5以下の数値になります。）
　そうすると、時刻t+Δtにおいて、x上の玉の数は時刻tにあった玉が移動せずに残った数と隣の位置から移動してきた玉の数の合計になります。
$f(x,t+\Delta t) = (1-2n)f(x,t) + nf(x-\Delta x,t)+nf(x+\Delta x,t)$

　これを変形すると
$f(x,t+\Delta t) - f(x,t) = n\{ f(x+\Delta x)+f(x-\Delta x)-2f(x)\}$
　右辺はf(x,t)のtは省略しました。これを両辺Δtで割ると

$\frac{f(x,t+\Delta t)-f(x,t)}{\Delta t} = \frac{n}{\Delta t}\{f(x+\Delta x)+f(x-\Delta x)-2f(x)\}$

$\frac{f(x,t+\Delta t)-f(x,t)}{\Delta t} = \frac{n\Delta x^2}{\Delta t}\frac{f(x+\Delta x)+f(x-\Delta x)-2f(x)}{\Delta x^2}$

　さて、ここでΔt→0を考えてみます。当然、時間間隔が短くなれば、この時間の間に移動できるボールの距離Δxは短くなる為、Δx→0になります。
　ということは、偏微分の定義より

$\frac{\partial f}{\partial t} = \frac{n\Delta x^2}{\Delta t} \frac{\partial^2 f}{\partial x^2}$

　さて、なんか微妙に $\frac{n\Delta x^2}{\Delta t}$ という部分が残っちゃいました。この値はΔtが大体Δxの二乗オーダーだったら、何らかの値に収束するわけです。というわけで、根拠はここでは特にないですが、勝手にこれをDという係数でおいてしまいましょう。一応、ざっくり説明をすると、時刻0にある点にあったボール達が広がった時、ボールが存在する面積と時間が比例すると考えると、ΔtはΔxの二乗オーダーになります。これはそんなに変な仮定ではありませんよね？

　nはただの定数の割合なので無次元ですから、Dの次元は[ $m^2/sec$ ]。というわけで、拡散を表す方程式は
$\frac{\partial f}{\partial t} = D \frac{\partial^2 f}{\partial x^2}$
　となります。

　ここで少しだけ、Dについて考えてみましょう。
$D= \frac{n\Delta x^2}{\Delta t}$
とおいたわけですから、Dの値が大きければ、Δtを定数としてみると、Δtの時間の間に玉が移動する距離Δxが長いことになります。逆にDの値が小さければ、Δtの時間の間に玉が移動する距離Δxが短いことになります。つまり、Dの値が大きいほど玉の移動距離が長い＝拡散が速いという事になります。

　さて、これで1次元の場合は終わりました。では、2次元、3次元ならどうなのでしょうか。2次元で見てみますと
$f(x,y,t+\Delta t) = (1-4n)f(x,y,t) + nf(x-\Delta x,y,t)+nf(x+\Delta x,y,t)+nf(x,y+\Delta y,t)+nf(x,y-\Delta y,t)$

とあらわせます。これって、上の議論を当てはめると
$\frac{\partial f}{\partial t} = D \frac{\partial^2 f}{\partial x^2}+ D \frac{\partial^2 f}{\partial y^2}$

という風に単純にxとyの偏微分の足し算って言うだけです。これをナブラを使って短く表現すると
$\frac{\partial f}{\partial t} = D \nabla^2 f$
と書き表すことができます。3次元になっても同じです。

　ここでは説明は完全に省きましたが、定数nが位置やfの値（玉の数）よって変化することがあります。その場合も考慮すると、拡散方程式は一般に
$\frac{\partial f}{\partial t} = \nabla\cdot(D(f,x)\nabla f)$
となります。この形式を数値解析などで解くのは中々難しいです。