はじめに
一般的な流体シミュレーションでは空間をグリッドに分割して流体を表現します．これを「格子法」と呼びます．
一方，流体を粒子で表現しようとするのが「粒子法」です．粒子法にはいろんなメリットがあります（こちらのサイトで詳しく説明されています）．何よりも，個人的には粒子法の方が物理的に自然なモデリングかなぁと思っています．
この記事は「MPS法解説シリーズ」の一部です．このシリーズでは粒子法の一種である「MPS法」を解説しています．初めましての方は，ぜひシリーズの全体像からご覧ください．
#目次
・やりたいこと
・シンプルなモデル
・妥当性の検討
・λの導出
・ベクトルのラプラシアンモデル
・まとめ
#やりたいこと
勾配モデルの記事で述べたように，流体シミュレーションにはスカラーの勾配モデルとラプラシアンモデルが必要です．この記事では，MPS法のラプラシアンモデルを導出します．係数 $\lambda$ の導出についていろいろ調べてみたのですが，「解析解と統計的に一致させるための係数」といった説明しか見つけられませんでした．「統計的に一致」ってなんだよ！！！
この記事では，単にラプラシアンモデルを紹介するだけでなく、 $\lambda$ の中身を数学的に導出することを目指します．
スカラー $\phi$ のラプラシアンモデル
\langle \boldsymbol{\nabla}^2 \phi \rangle _i
= \frac{2d}{\lambda}
\sum_{j\ne i} \phi_{ij} w_{ij}\\
\lambda = \frac{1}{n^0}
\sum_{j\ne i} |\boldsymbol{r}_{ij}|^2 w_{ij}
#シンプルなモデル
ラプラシアンの物理的な意味は「拡散」，すなわち物理量の受け渡しです．そこで，物理量 $\phi$ に重み関数を掛けた値が粒子間で受け渡されると考えてみましょう．
この時，粒子 $i$ は粒子 $j$ に $\phi_i w_{ij}$ だけの物理量を受け渡し，粒子 $j$ から $\phi_j w_{ji}$ だけの物理量を受け取ると考えられそうです．すると，ラプラシアンは次のように表せます．
\begin{align}
\langle \boldsymbol{\nabla}^2 \phi \rangle _i
&= \frac{d}{n^0} \sum_{j\ne i} \{
\phi_j w_{ji}
- \phi_i w_{ij} \} \\
&= \frac{d}{n^0} \sum_{j\ne i} \phi_{ij} w_{ij}
\quad \because w_{ji}=w_{ij}
\end{align}
なお，勾配モデル・発散モデルと同様に，次元数 $d$ をかけて粒子数密度の基準値 $n^0$ で割っています．こちらは物理的に納得がいくシンプルなモデルですね．
#妥当性の検討
上記のラプラシアンモデルの妥当性を検討するため，ラプラシアンを用いた微分方程式の解析解と，その方程式にラプラシアンモデルを適用した場合の解を比較します．用いる微分方程式は，スカラー $\phi$ に関する非定常拡散方程式
\frac{\partial \phi}{\partial t}
= D \boldsymbol{\nabla}^2 \phi
です．$D$ は拡散係数を表します．
粒子法では，各粒子が物理量を持つ状態はデルタ関数の重ね合わせと考えられます．そこで，初期分布を原点におけるデルタ関数とした場合の非定常拡散方程式の解析解を求めます．この解はガウス分布になることが知られています．
\phi(r)=\Biggl( \frac{1}{\sqrt{4\pi Dt}} \Biggr)^d
\exp{\Biggl( -\frac{r^2}{4Dt}} \Biggr)
つまり，デルタ関数が拡散するとガウス分布になるんです．美しいですね．したがって，物理量をガウス分布にしたがって配分するモデルを用いればモデルの解は解析解と一致するはずです．
しかし，ガウス分布を計算に用いるのには問題があります．ガウス分布は非有界の影響範囲を持つ（compuct supportでない）のですべての粒子に関して計算する必要があり，膨大な計算時間がかかってしまうのです．
そこで登場するのが「中心極限定理」です．中心極限定理とは，
「ガウス分布とは異なる分布であっても，その分布による分配操作を多数回反復すればその分布形はガウス分布に収束する」
というものです．例えば，サイコロを $n$ 回振った際に出た目の和が分布する形状は， $n=1$ でこそ直線的ですが， $n$ が増加するにつれてガウス分布に近づいていく様子がイメージできますね．本来は一部の例外があるそうですが，我々は（積分内の極限操作を外に取り出すことさえ気にしない）工学屋です．細かい話は数学屋さんに任せて先に進みましょう．
中心極限定理を考慮すると，上述のように重み関数を用いて分配するモデルであっても問題はなさそうです（もっと言うと，重み関数でなくてもたいていは問題ない）．しかし，一口にガウス分布といってもとんがったものから平べったいものまであります．それらを特徴づけるのが「統計的な分散 $\sigma^2$ 」です．MPS法のラプラシアンモデルに登場する $\lambda$ は，この $\sigma^2$ を一致させるために必要な係数だったんですね．分散について知りたい方やその定義についてはそういったサイトをご覧いただくとして，以降では $\lambda$ の中身を導出していきます．
λ の導出
$\lambda$ の具体的な形を求めるために，分散を解析解と一致させるための係数 $\alpha$ をかけたラプラシアンモデル
\langle \boldsymbol{\nabla}^2 \phi \rangle _i
= \alpha \frac{d}{n^0} \sum_{j\ne i} \phi_{ij} w_{ij}
を考えて，解析解の分散との比較から係数 $\alpha$ を求めます．
ラプラシアンモデルから分散を計算
時刻 $k$ で粒子 $i$ のみが大きさ1の物理量を持つ状態を考えましょう．
\phi_j^k=\Biggl\{
\begin{array}{ll}
1&(j=i)\\
0&(j\ne i)
\end{array}
この時，粒子 $i$ を中心とする分散 $(\sigma^2)_i^k$ は，
\begin{align}
(\sigma^2)_i^k
&=\sum_{j\ne i}|\boldsymbol{r}_{ij}|^2\phi_j^k\\
&=0
\end{align}
と求まります．この時，ラプラシアンモデルを用いて非定常拡散方程式を離散化して次の時刻の物理量 $\phi^{k+1}$ を求めると，
\begin{align}
\langle \phi \rangle _j^{k+1}
&=\alpha \frac{d}{n^0} D \Delta t
\sum_{l\ne j} \phi_{jl}^k w_{jl}\\
&=\alpha \frac{d}{n^0} D \Delta t w_{ij}
\end{align}
となるので，この時刻の分散 $(\sigma^2)_i^{k+1}$ は，
\begin{align}
(\sigma^2)_i^{k+1}
&= \sum_{j\ne i} |\boldsymbol{r}_{ij}|^2 \langle \phi \rangle _j^{k+1} \\
&=\alpha \frac{d}{n^0} D \Delta t
\sum_{j\ne i} |\boldsymbol{r}_{ij}|^2 w_{ij}
\end{align}
と求まります．つまり，時間 $\Delta t$ での分散の増加量 $\Delta(\sigma^2)_i$ は次のように求まります．
\begin{align}
\Delta (\sigma^2)_i
&=(\sigma^2)_i^{k+1}-(\sigma^2)_i^k \\
&=\alpha \frac{d}{n^0} D \Delta t
\sum_{j\ne i} |\boldsymbol{r}_{ij}|^2 w_{ij}
\end{align}
解析解から分散を計算
次に，物理量 $\phi$ が先ほどと同様に
\phi_j^k=\Biggl\{
\begin{array}{ll}
1&(j=i)\\
0&(j\ne i)
\end{array}
として分布している状態を初期配置とする非定常拡散方程式の解析解を求めましょう．このような分布は，粒子 $i$ の位置にデルタ関数が存在する状態と考えられます．したがって，粒子 $i$ を中心とする分散 $(\sigma^2)_i$ の計算に，先ほどの解析解
\phi(r)=\Biggl( \frac{1}{\sqrt{4\pi Dt}} \Biggr)^d
\exp{\Biggl( -\frac{r^2}{4Dt}} \Biggr)
を使うことができます．
\begin{align}
(\sigma^2)_i
&=\int r^2 \phi(r)dv\\
&=2dDt
\end{align}
つまり，時間 $\Delta t$ での分散の増加量 $\Delta(\sigma^2)_i$ は次のように求まります．
\Delta(\sigma^2)_i=2dD \Delta t
両者の比較
最後のステップです．モデルの解と解析解を比較します．
\begin{align}
&\alpha \frac{d}{n^0} D \Delta t
\sum_{j\ne i} |\boldsymbol{r}_{ij}|^2 w_{ij}
=2dD \Delta t\\
&\therefore \alpha=2\frac{n^0}{\sum_{j\ne i} |\boldsymbol{r}_{ij}|^2 w_{ij}}
\end{align}
ようやく $\alpha$ が求まりました．これを最初の式に代入すれば，MPS法のラプラシアンモデルが導出されると同時に，係数 $\lambda$ の中身も示されます．
#ベクトルのラプラシアンモデル
勾配モデルの時と同じように，ベクトルにラプラシアンを適用した場合の表式も求めておきます．
ベクトル $\boldsymbol{u}$ のラプラシアンは次のように表されます．
\boldsymbol{\nabla}^2\boldsymbol{u}
= \begin{pmatrix}
\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} \\
\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} + \frac{\partial^2 v}{\partial z^2} \\
\frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2}
\end{pmatrix}
= \begin{pmatrix}
\boldsymbol{\nabla}^2u \\
\boldsymbol{\nabla}^2v \\
\boldsymbol{\nabla}^2w \\
\end{pmatrix}
つまり，それぞれの成分ごとにラプラシアンを取ってやればいいわけですね．上で求めたラプラシアンモデルを適用してみましょう．
math \begin{align} \langle\boldsymbol{\nabla}^2\boldsymbol{u}\rangle_i &= \frac{2d}{\lambda}\sum_{j\ne i} \begin{pmatrix} u_{ij} \\ v_{ij} \\ w_{ij} \end{pmatrix} w_{ij} \\ &= \frac{2d}{\lambda}\sum_{j\ne i}\boldsymbol{u}_{ij}w_{ij} \end{align}
ちょー簡単ですね．そりゃそうだって感じです．
#まとめ
この記事ではMPS法のラプラシアンモデルを導出しました．
まずはシンプルで分かりやすいモデルを考えました．このモデルの妥当性を検討すると，「中心極限定理」によって大まかには問題ないことが分かりました．ただし，分布形状を解析解と一致させる必要があります．これを考慮して，ようやく $\lambda$ の導出ができました．
$\lambda$ の詳細な導出は探しても出てこなかったので，ずっともやもやしてました．（私にしては）深めの数学で大変でしたが，答えが求まったときは最高に幸せでした．皆さんもぜひ同じ気持ちを味わって頂ければと思います．