拡散方程式
拡散方程式は、その名の通り、物質の拡散を表現した偏微分方程式です。
例えば二次元なら、以下のようになります。
\frac { { \partial }u(x,y,t) }{ { \partial }{ t } } =D(\frac { { \partial }^{ 2 }u(x,y,t) }{ { \partial }{ x }^{ 2 } } +\frac { { \partial }^{ 2 }u(x,y,t) }{ { \partial }y^{ 2 } } )
Dについては比例定数です
物質の流れる速度は濃度勾配に比例します。
わりとなんでもそうですが、濃度というのは濃い方から薄い方へ流れていきます。
なんでかというと、例えば空間を細切れにして小さな部屋がいっぱいある様子をイメージします。
そして隣り合う２つの部屋を考えます。
原子はあらゆる方向へ無作為に移動を繰り返しているため、たくさん原子がある部屋からはたくさん原子が飛んできますが、何もない部屋からは出ていくものがありませんから、やはり入ってくるばかりです。逆に２つの部屋の原子の数が同じなら、出て行く量と入ってくる量は均衡していることが想像できます。
式のイメージはこれくらいにして、
導出のほうは、拡散方程式の手短な導出あたりを参照するとして、早速離散化して動かしてみます。
例のごとく、差分法でやってみましょう。
できました。
\frac { { \partial }u(x,y,t) }{ { \partial }{ t } } =D\left(\frac { { \partial }^{ 2 }u(x,y,t) }{ { \partial }{ x }^{ 2 } } +\frac { { \partial }^{ 2 }u(x,y,t) }{ { \partial }y^{ 2 } } \right)\\ \frac { { \partial }u(x,y,t) }{ { \partial }{ t } } =D(\frac { u(x-\Delta x,y,t)-2u(x,y,t)+u(x+\Delta x,y,t) }{ { \Delta x }^{ 2 } } +\frac { u(x,y-\Delta y,t)-2u(x,y,t)+u(x,y+\Delta y,t) }{ { \Delta y }^{ 2 } } )\\ \Delta x=\Delta y\\ \frac { { \partial }u(x,y,t) }{ { \partial }{ t } } =\frac { D }{ { \Delta x }^{ 2 } } \left\{ u(x-\Delta x,y,t)+u(x+\Delta x,y,t)+u(x,y-\Delta x,t)+u(x,y+\Delta x,t)-4u(x,y,t) \right\} \\ \frac { u(x,y,t+\Delta t)\quad -\quad u(x,y,t) }{ \Delta t } =\frac { D }{ { \Delta x }^{ 2 } } \left\{ u(x-\Delta x,y,t)+u(x+\Delta x,y,t)+u(x,y-\Delta x,t)+u(x,y+\Delta x,t)-4u(x,y,t) \right\} \\ u(x,y,t+\Delta t)=u(x,y,t)\quad +\frac { D\Delta t }{ { \Delta x }^{ 2 } } \left\{ u(x-\Delta x,y,t)+u(x+\Delta x,y,t)+u(x,y-\Delta x,t)+u(x,y+\Delta x,t)-4u(x,y,t) \right\} \\
あとはコードに落とし込みましょう。
float deltaX = float(kWaveWidth) / float(kWaveGrid);
float deltaT = 0.5f;
float c = D * deltaT / (deltaX * deltaX);
std::fill(_values.begin(), _values.end(), 0.0f);
for (int x = 1; x < kWaveGrid - 1; ++x) {
for (int y = 1; y < kWaveGrid - 1; ++y) {
int index = this->valueIndex(x, y);
float uL = _u[this->valueIndex(x - 1, y)];
float uR = _u[this->valueIndex(x + 1, y)];
float uT = _u[this->valueIndex(x, y - 1)];
float uB = _u[this->valueIndex(x, y + 1)];
float u = _u[index];
_values[index] = u + c * (uL + uR + uT + uB - 4 * u);
}
}
std::swap(_values, _u);
本質的なコードはこれだけです。
さくっとビジュアライズすると、
いい感じです。
※マウスの動きは単調なので、コクのある乱数でペンの位置を振動させています
※思考停止で前方差分で解いていますが、大きなdeltaTなどをとると不安定になります。これを嫌う場合は陰解法をつかうのがわりと定石？のようです。実際Real-Time Fluid Dynamics for Gamesによると拡散の計算にはガウス＝ザイデル法を使うと良いよ、とあります。
謝辞
@SatoshiTerasaki さん、数式の修正ありがとうございました！