ラプラス方程式とは, \[\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}=0\]のように, ラプラシアン\[\Delta=\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\]を用いて\[\Delta u=0\]と書ける微分方程式です．なお, 既知関数\(f\)を用いて\[\Delta u=f\]となると, ポアソン(お魚ではなく, Siméon Denis Poisson)方程式と呼ばれます. 

• 流体力学では, Poiseuille流とかCouetteの時に出てきました．
• 調和関数とか言って複素関数論で出てきました．
• 完全流体の流れでも出てきました．
• 温度分布の方程式でもあります．
• 数値解析とかデザインとかの格子作成でも使います．
• まあ，ようするに，至る所で使います．

多次元配列

差分法でこれを解いてみましょう. 変数が2つあるので, 離散化は\[u(\Delta i,\Delta j)=u_{i,j}\]のように, 添字が複数ある多次元のデータになるはずです. 多次元の配列を取り扱う方法を考えてみましょう.

2次元ならEigen

2次元までしか使えませんが, EigenのMatrixを利用すると良いですね. やり方はここと同じです.

多次元配列は神への冒涜である派

そもそもコンピュータのメモリーは1次元的に「アドレス」に並んでいるものである. だから, u[i,j]のような多次元配列はマヤカシでありSyntaxSugarであり, 原理に対する許し難い冒涜と感じる原理主義者もいるでしょう．長所は計算がBoost並みに早いこと，短所は原理主義者は通常嫌われることです．とはいえ, 実戦でこれを理解していないと何かと不便なので, ここで試してみましょう．

Boost/Array

出来合いのBoost/Arrayを用いると，自分で考えなくても多次元配列が利用できます．長所は自分で書かなくても良いしBoostだから計算が早い，短所はマニュアルを見なければ使い方がわからない，ところです．ここで試してみましょう．

配列？ってな〜に派

そもそもデータを記録するのに配列が必要であることが理解できない異世界人もいるでしょう．まあN次元メモリーが実現されたら, 意外と正解かもしれませんけどね.

差分スキーム

多次元データを格納できる見通しがついたら, 差分スキームを考えます. ここの境界値問題で述べた通り\[\frac{\partial^2 u}{\partial x^2}\left(\Delta i, \Delta j\right)\sim\frac{u_{i-1,j}-2u_{i,j}+u_{i+1,j}}{\Delta^2}\]\[\frac{\partial^2 u}{\partial y^2}\left(\Delta i, \Delta j\right)\sim\frac{u_{i,j-1}-2u_{i,j}+u_{i,j+1}}{\Delta^2}\]であるから, 要は\[u_{i,j}=\frac{u_{i-1,j}+u_{i+1,j}+u_{i,j-1}+u_{i,j+1}}{4}\]を解けば良いことになります.

差分スキームを絵で示すと

こんな感じですね. 赤色格子点が未知数で, 同じ数だけの差分式があります．ですが，これを解くには, 緑色格子点の境界条件が既知でなければなりません. こういうのをDirichlet問題というのでしたね. 

定理：滑らかな有限な領域\(\Omega\)におけるLaplace方程式\[\Delta u=0\]は, なめらかな境界\(\partial \Omega\)における境界条件\[u=\phi\]を指定すれば, 滑らかな解は一意に持つ.

例：\(\cfrac{\text{d}^2 T}{\text{d}x^2}=0\) は, 有限領域\(0<x<1\)ならば, 例えば \(T(0)=0\), \(T(1)=1\)で解を求めることができる. てゆうか境界値がなんであっても解を求めることができる. ただし領域が\(0<x<\infty\)とかだと, ハマる可能性がある. \(0<x<0.5\)では\(T=4x\), \(0.5<x<1\)では\(T=3-2x)\)も, 1点で滑らかではないが, 解ではある. だが, そんなこと許したら, なんぼでも解が見つかってしまう.

ついでにNeumann問題ってのは

定理：滑らかな有限な領域\(\Omega\)におけるLaplace方程式\[\Delta u=0\]は, なめらかな境界\(\partial \Omega\)における境界条件\[{\boldsymbol n}\cdot\text{grad } u=\psi\]を指定すれば, 未定定数を除けば, 滑らかな解は一つだけ持つ. ただし, \[\int_{\partial\Omega} {\boldsymbol n}\cdot\text{grad } u \text{d}S=0\]が成立しないとやばい.

例：\(\cfrac{\text{d}^2 T}{\text{d}x^2}=0\) は, \(\cfrac{\text{d}T}{\text{d}x}(0)=\cfrac{\text{d}T}{\text{d}x}(1)=1\)で解を求めることができるのじゃが, \(T=x+\text{Const.}\)のように未定定数が残る.

\(\cfrac{\text{d}T}{\text{d}x}(0)=1, \cfrac{\text{d}T}{\text{d}x}(1)=0\)で解くのは無理ゲーだろ

Neumann問題の差分式を考えてみよ.

逐次近似解法

ここで学んだJacobi法で解いてみましょう. まあ要するに, \(n=0\)の初期値\(u^{(0)}_{i,j}\)をテキトーに選んでおき,\[u^{(n+1)}_{i,j}=\frac{u^{(n)}_{i-1,j}+u^{(n)}_{i+1,j}+u^{(n)}_{i,j-1}+u^{(n)}_{i,j+1}}{4}\]を繰り返していけば, 神のご加護があれば解に収束するはずである.

解答例は, こちら. こいつは, 100イタレーション毎に, \[\text{max}|u^{(n+1)}_{i,j}-u^{(n)}_{i,j}|\]を出力する設定なのだ. 実行結果はこんな感じである:

Step:Diff= 100 0.000132403
Step:Diff= 200 8.75579e-07
Step:Diff= 300 5.79332e-09
Step:Diff= 400 3.83319e-11
Step:Diff= 500 2.53686e-13
Step:Diff= 600 1.72085e-15
Finish Step:Diff= 610 9.99201e-16

ふむ. なかなか良さそうである. だが，結果をどうやって眺めたものか・・・\(u(x,y)\)は2次元関数なので, graphソフトでは対応できないな．

可視化

やむを得まい. 多次元関数\(u(x,y)\)を可視化する方法を考える.

オイラはMATLAB/OCTAVE

2次元くらいまでなら，MATLABとかOCTAVEでも描けそうな気がする. ここで試してみる.

出力はこんな感じ:

最終兵器ParaView

やはり複雑化すると立派な描画ソフトが必要であるね. ここで試してみる.

出力はこんな感じ：

一発解法

いやー，逐次近似は嫌いだ．メモリーは最近安いんで，一発で解いてほしい．という願いは当然あるでしょう．

差分式は連立方程式なので, 必ず\[\text{A}{\boldsymbol x}={\boldsymbol b}\]のように書ける. ここで, \(\boldsymbol{x}\)に該当するものは, 求めるべき\(u_{i,j}\)を一列化したベクトルである.

なんのことはない, \(u\)に関しては多次元配列を利用するのは無駄であった．原理主義に基づき, 1列ベクトルで作成する.

だが，先に進む前に考えるべきことがある．例えば 500x500 格子で計算する場合・・・・未知数の数は500x500=250000あるので・・・連立方程式を格納する行列サイズは 250000x250000 ってわけで・・・必要な計算機メモリーは, \[\frac{8\times 500^4}{10^9}\ \text{GB}\]となります．要するに．496GBです．あなたのノートPCのメモリーが512GB以上なら，進んでよし．実際には

　

ま，こんなものですね．というわけで，ここでやらかしたような，疎行列を面倒だからといって密行列でコードを書いて俺は楽ちん，コンピュータ大変，では動作しないわけです．

で，どうしましょうか．

• 数値計算業界において，力は，正義である．
• 力とは，搭載メモリーサイズのことである．
• こないだAmazonで見たご家庭用マザーボードの最大搭載メモリーは2TBであった．時代は変わるのだ．

という方針では，でっかい計算機を使う，というのが確実な手です．が，まあ，まばらにしか要素が存在しない疎行列であるので，まあなんとかなる場合もあります．この問題については, 次の節で議論しましょう.