拡散方程式の数値解法

  • 5
    Like
  • 0
    Comment

拡散方程式

拡散方程式は、その名の通り、物質の拡散を表現した偏微分方程式です。
例えば二次元なら、以下のようになります。

\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については比例定数です

物質の流れる速度は濃度勾配に比例します。
わりとなんでもそうですが、濃度というのは濃い方から薄い方へ流れていきます。
なんでかというと、例えば空間を細切れにして小さな部屋がいっぱいある様子をイメージします。
そして隣り合う2つの部屋を考えます。
原子はあらゆる方向へ無作為に移動を繰り返しているため、たくさん原子がある部屋からはたくさん原子が飛んできますが、何もない部屋からは出ていくものがありませんから、やはり入ってくるばかりです。逆に2つの部屋の原子の数が同じなら、出て行く量と入ってくる量は均衡していることが想像できます。

式のイメージはこれくらいにして、
導出のほうは、拡散方程式の手短な導出あたりを参照するとして、早速離散化して動かしてみます。
例のごとく、差分法でやってみましょう。
できました。

\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 } } )\\ \frac { { \partial  }u(x,y,t) }{ { \partial  }{ t } } =D(\frac { u(x-⊿x,y,t)-2u(x,y,t)+u(x+⊿x,y,t) }{ { ⊿x }^{ 2 } } +\frac { u(x,y-⊿y,t)-2u(x,y,t)+u(x,y+⊿y,t) }{ { ⊿y }^{ 2 } } )\\ ⊿x=⊿y\\ \frac { { \partial  }u(x,y,t) }{ { \partial  }{ t } } =\frac { D }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} \\ \frac { u(x,y,t+⊿t)\quad -\quad u(x,y,t) }{ ⊿t } =\frac { D }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿x,t)-4u(x,y,t) \right\} \\ u(x,y,t+⊿t)=u(x,y,t)\quad +\frac { D⊿t }{ { ⊿x }^{ 2 } } \left\{ u(x-⊿x,y,t)+u(x+⊿x,y,t)+u(x,y-⊿x,t)+u(x,y+⊿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);

本質的なコードはこれだけです。
さくっとビジュアライズすると、

diffuse.gif

いい感じです。

※マウスの動きは単調なので、コクのある乱数でペンの位置を振動させています
※思考停止で前方差分で解いていますが、大きなdeltaTなどをとると不安定になります。これを嫌う場合は陰解法をつかうのがわりと定石?のようです。実際Real-Time Fluid Dynamics for Gamesによると拡散の計算にはガウス=ザイデル法を使うと良いよ、とあります。