キーワード
おさらい
『偏微分方程式の数値計算入門（2/3）』の続きです．テイラー展開とか使うと格子座標\( i,j \)についての微分係数を簡単な多項式に置き換えることができた．
1階偏微分
\[ \small{ \frac{\partial \psi(x_{i,j}, y_{i,j})}{\partial x} \approx \frac{\psi_{i+1, j} - \psi_{i, j}}{h}, \ \frac{\partial \psi(x_{i,j}, y_{i,j})}{\partial y} \approx \frac{\psi_{i, j+1} - \psi_{i, j}}{h} （前進差分近似）} \]
\[ \small{ \frac{\partial \psi(x_{i,j}, y_{i,j})}{\partial x} \approx \frac{\psi_{i, j} - \psi_{i-1, j}}{h}, \ \frac{\partial \psi(x_{i,j}, y_{i,j})}{\partial y} \approx \frac{\psi_{i, j} - \psi_{i, j-1}}{h} （後進差分近似）} \]
\[ \small{ \frac{\partial \psi(x_{i,j}, y_{i,j})}{\partial x} \approx \frac{\psi_{i+1, j} - \psi_{i-1, j}}{2h}, \ \frac{\partial \psi(x_{i,j}, y_{i,j})}{\partial y} \approx \frac{\psi_{i, j+1} - \psi_{i, j-1}}{2h} （中心差分近似）} \]
2階偏微分
\[
\small{
\frac{\partial^2 \psi(x_{i,j}, y_{i,j})}{\partial x^2} \approx \frac{\psi_{i+1, j} - 2\psi_{i, j} + \psi_{i-1, j}}{h^2}, \frac{\partial^2 \psi(x_{i,j}, y_{i,j})}{\partial y^2} \approx \frac{\psi_{i, j+1} - 2\psi_{i, j} + \psi_{i, j-1}}{h^2}
}
\]
拡散方程式を数値計算してみる
前回登場したラプラス方程式では，\( \nabla^2 \psi \)に"0"がイコールで繋がっていました．\( \nabla^2 \psi \)は"拡散項"と呼ばれることが多いです．今回は拡散項に非定常項を繋げてみたいと思います． また，単位時間（微小時間）における拡散の度合いを示す係数\( k \)を拡散項に乗じると以下のような式になります．\[
\frac{\partial \psi}{\partial t} = k \nabla^2 \psi
\]
この式には名前がついており，一般に"拡散方程式"と呼ばれています．ラプラス方程式で得られる解は「収束した最終的な状態」でしたが，非定常項がつくことで収束に至る経過を時間\( t \)の関数として得ることができます． つまり，「何秒後にはこんな状態になっている」という計算ができるということになります．また，境界条件を時間的に変化させたり，計算の自由度（次元）が増えることになります．
今回も2次元を考えているので\( \nabla^2 \)を展開すると下式のようになります．
\[
\frac{\partial \psi}{\partial t} = k \left( \frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} \right)
\]
このとき，\( \psi \)は時間と空間座標の関数になっていることに注意が必要です．
\[
\psi = \psi (x, y, t)
\]
時間微分を前進差分近似，空間2階微分を中心差分近似で多項式に置き換えます．
\[
\frac{\psi_{i,j}^{(n+1)} - \psi_{i,j}^{(n)}}{\Delta t} = k \left( \frac{\psi_{i+1, j}^{(n)} - 2\psi_{i, j}^{(n)} + \psi_{i-1, j}^{(n)}}{h^2} + \frac{\psi_{i, j+1}^{(n)} - 2\psi_{i, j}^{(n)} + \psi_{i, j-1}^{(n)}}{h^2} \right)
\]
\( \psi_{i,j}^{(n+1)} \)について解くと，
\[
\psi_{i,j}^{(n+1)} = \psi_{i,j}^{(n)} + \frac{k \Delta t}{h^2}(\psi_{i+1, j}^{(n)} + \psi_{i-1, j}^{(n)} + \psi_{i, j+1}^{(n)} + \psi_{i, j-1}^{(n)} - 4\psi_{i, j}^{(n)})
\]
この漸化式を繰り返し計算することで，時間\( t \)における\( \psi \)の状態を知ることができます．それでは実際に，初期条件と境界条件を定義して計算してみます．
まずは計算範囲を定義します．
\[
\frac{\partial \psi}{\partial t} = k \left( \frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} \right) ; 0 \leq x \leq 1, 0 \leq y \leq 1
\]
初期条件はこんな正規分布の形にしてみたいと思います．
\[
\psi(x, y, 0) = \exp \left\{ -100 \left( x - \frac{1}{2} \right)^2 \right\}
\]
境界条件も同様に，
\[
\psi(0, y, t) = \psi(1, y, t) = \psi(x, 0, t) = \psi(x, 1, t) = \exp \left\{ -100 \left( x - \frac{1}{2} \right)^2 \right\}
\]
その他，パラメータを以下のように定義します．実はこのパラメータの設定が結構シビアで，あまり大きい値にすると計算誤差による発散を起こすので注意が必要です．
\[
\Delta t = 10^{-2} \\
k = 10^{-4}
\]
こちらに計算プログラムを用意してみました．（ブラウザで動きます）
頑張って計算するとこんな感じになります．
波動方程式を数値計算してみる
拡散方程式では非定常項が「時間の1階微分」でしたが，非定常項を「時間の2階微分」としたらどうなるのでしょうか．つまり，以下のような式です．\[
\frac{\partial^2 \psi}{\partial t^2} = k \nabla^2 \psi
\]
実はこの式にも名前がついていて，"波動方程式"と呼ばれるものになります．中二病感溢れるネーミングでモチベーションが上がってきますよね． 名前の通り，この式は水面や膜の振動をモデル化しているものらしい．
このとき\( k \)は以下のように表されます．
\[
k = s^2
\]
\( s \)は振動の位相速度と解釈されることが多いです．要するに「波が伝わる速度」のことですね．
ここでも同じように2次元を考えていきたいので以下のように\( \nabla^2 \)を展開します．
\[
\frac{\partial^2 \psi}{\partial t^2} = k \left( \frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} \right)
\]
時間微分，空間微分を中心差分近似で多項式に置き換えます．
\[
\small{
\frac{\psi_{i,j}^{(n+1)} - 2 \psi_{i,j}^{(n)} + \psi_{i,j}^{(n-1)}}{\Delta t^2} = k \left( \frac{\psi_{i+1, j}^{(n)} - 2\psi_{i, j}^{(n)} + \psi_{i-1, j}^{(n)}}{h^2} + \frac{\psi_{i, j+1}^{(n)} - 2\psi_{i, j}^{(n)} + \psi_{i, j-1}^{(n)}}{h^2} \right)
}
\]
\( \psi_{i,j}^{(n+1)} \)について解くと，
\[
\small{
\psi_{i,j}^{(n+1)} = 2 \psi_{i,j}^{(n)} - \psi_{i,j}^{(n-1)} + \frac{k \Delta t^2}{h^2}(\psi_{i+1, j}^{(n)} + \psi_{i-1, j}^{(n)} + \psi_{i, j+1}^{(n)} + \psi_{i, j-1}^{(n)} - 4\psi_{i, j}^{(n)})
}
\]
この漸化式を繰り返し計算することで，時間\( t \)における\( \psi \)の状態を知ることができます．今回は時間を2階微分しているので，\( \psi \)の状態を格納する配列について，\( \psi^{(n-1)}, \psi^{(n)}, \psi^{(n+1)} \)の3つを用意する必要があります．
今回は格子の大きさ\( h \)を指定します．
\[
h = 0.1
\]
格子数400×400のときの計算範囲は，
\[
\frac{\partial^2 \psi}{\partial t^2} = k \left( \frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} \right) ; 0 \leq x \leq 400h, 0 \leq y \leq 400h
\]
初期条件はこんな感じに定義してみました．
\[
\small{
\psi(x, y, 0) = 20 \cdot \exp \left\{ -10\left( x - \frac{400h}{3} \right)^2 \right\} \cdot \exp \left\{ -10\left( y - \frac{400h}{4} \right)^2 \right\}
}
\]
端っこ（境界条件）は0の固定値としました．こうすることでいい感じに波紋が反射します．ちなみに，境界を特定の値に固定するような条件をディリクレ境界条件と呼ぶんだそうだ．
\[
\psi(0, y, t) = \psi(400h, y, t) = \psi(x, 0, t) = \psi(x, 400h, t) = 0
\]
その他，パラメータを以下のように定義します．
\[
\Delta t = 1 / 60 \\
k = 1
\]
こちらに計算プログラムを用意してみました．（ブラウザで動きます）
きれいだねえ．
最後に
今回は思ったより長くなってしまったのですが，基本的な部分はこれで一通り触れられているんじゃないかと思います．全体を通して考え方は非常にシンプルなので，そこまで理解が難しいものではなかったと思います．ここまで理解できていれば，外力の作用，障害物の設置みたいなこともどんどん考えていけると思います．何かの役に立てて頂ければ幸いです．