Euler-Lagrange 方程式からSPH法の運動方程式を求める

[SPH法の理論と実装 ~理論編1~ - Qiita](https://qiita.com/mitchiinaga/items/1e88c22c1708b19290cf) ですっ飛ばしたところを少し詳しく書きます。自分用のメモも兼ねて。
より正確に知りたい人は、Rosswog (2009) や Springel (2010) を眺めたほうがいいです。

# 記号
|記号 |説明 |
|----|---|
|$\boldsymbol{r}$|位置|
|$\boldsymbol{v}$|速度|
|$m$|質量|
|$p$|圧力|
|$\rho$|密度|
|$u$|単位質量当たり内部エネルギー|
|$L$|ラグランジアン|
|$W(r,h)$|カーネル関数|
|$h$|smoothing length|
|$d$|次元|

# 導出
ラグランジアン
$$
L = \sum_j \left( m_j \frac{\boldsymbol{v}^2_j}{2} - m_j u_j \right)
$$
に対して、Euler-Lagrange 方程式

```math
\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\boldsymbol{r}}_i} \right) - \frac{\partial L}{\partial \boldsymbol{r}_i} = 0
```
を考えます。第1項は簡単で、

```math
\frac{d}{dt} \left( \frac{\partial L}{\partial \dot{\boldsymbol{r}}_i} \right) = m_i \frac{d \boldsymbol{v}_i}{dt}
```
です。
第2項は

```math
\frac{\partial L}{\partial \boldsymbol{r}_i} = - \sum_j m_j \frac{\partial u_j}{\partial \boldsymbol{r}_i} = - \sum_j m_j \frac{p_j}{\rho_j^2} \frac{\partial \rho_j}{\partial \boldsymbol{r}_i}
```
となります。等エントロピーを仮定し、
$$
du = -p dV = -p d\left( \frac{1}{\rho} \right) = \frac{p}{\rho^2} d\rho
$$
の関係式 (熱力学第一法則)を使用しています。SPH法の密度の計算式

```math
\rho_j = \sum_k m_k W_{jk}\\
W_{jk} \equiv W(|\boldsymbol{r}_j - \boldsymbol{r}_k|, h)
```
から

```math
\frac{\partial \rho_j}{\partial \boldsymbol{r}_i} =
\begin{cases}
\sum_k m_k \nabla_i W_{ik} & (i = j)\\
m_i \nabla_i W_{ij} & (i \neq j)
\end{cases}
```
となり、

```math
\frac{\partial L}{\partial \boldsymbol{r}_i} = - m_i \frac{p_i}{\rho_i^2} \sum_k m_k \nabla_i W_{ik} - \sum_{j \neq i} m_j \frac{p_j}{\rho_j^2} m_i \nabla_i W_{ij}
```
と変形できます。添字を付け替えると

```math
\frac{\partial L}{\partial \boldsymbol{r}_i} = - m_i \sum_j m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla_i W_{ij}
```
のように整理できます。ここで、 $ \nabla_i W_{ii} = 0 $であることを使っています。
第1項と第2項を合わせて、運動方程式

```math
m_i \frac{d \boldsymbol{v}_i}{dt} = - m_i \sum_j m_j \left( \frac{p_i}{\rho_i^2} + \frac{p_j}{\rho_j^2} \right) \nabla_i W_{ij}
```
が導出できます。

# 可変smoothing lengthの場合の導出
smoothing length が粒子ごとに与えられる場合を考えます。密度の計算式は

```math
\rho_j = \sum_k m_k W_{jk}(h_j)\\
W_{jk}(h_j) \equiv W(|\boldsymbol{r}_j - \boldsymbol{r}_k|, h_j)
```
を使います。このとき、密度の空間微分は、smoothing lengthが変数であることを考慮して

```math
\frac{\partial \rho_j}{\partial \boldsymbol{r}_i} = \sum_k m_k \left[ \nabla_i W_{jk}(h_j) + \frac{\partial W_{jk}(h_j)}{\partial h_j} \frac{\partial h_j}{\partial \rho_j} \frac{\partial \rho_j}{\partial \boldsymbol{r}_i} \right]
```
となります。smoothing length の条件として

$$
\rho_j h_j^d = {\rm const.}
$$
を与えると、以下のように整理できます。

```math
\frac{\partial \rho_j}{\partial \boldsymbol{r}_i} = f_j \times
\begin{cases}
\sum_k m_k \nabla_i W_{ik}(h_i) & (i = j)\\
m_i \nabla_i W_{ij}(h_j) & (i \neq j)
\end{cases}\\

f_j = \left(1 + \frac{h_j}{\rho_j d} \sum_k m_k \frac{\partial W_{jk}(h_j)}{\partial h_j} \right)^{-1}
```
これを使うと、Euler-Lagrange方程式の第2項は

```math
\frac{\partial L}{\partial \boldsymbol{r}_i} = - m_i \sum_j m_j \left( f_i \frac{p_i}{\rho_i^2} \nabla_i W_{ij}(h_i) + f_j \frac{p_j}{\rho_j^2} \nabla_i W_{ij}(h_j) \right)
```
となるので、運動方程式が

```math
m_i \frac{d \boldsymbol{v}_i}{dt} = - m_i \sum_j m_j \left( f_i \frac{p_i}{\rho_i^2} \nabla_i W_{ij}(h_i) + f_j \frac{p_j}{\rho_j^2} \nabla_i W_{ij}(h_j) \right)
```
に変わります。

# 参考文献
* Springel, V. (2010). Smoothed Particle Hydrodynamics in Astrophysics. Annual Review of Astronomy and Astrophysics, 48(1), 391–430. https://doi.org/10.1146/annurev-astro-081309-130914
* Rosswog, S. (2009). Astrophysical smooth particle hydrodynamics. New Astronomy Reviews, 53(4–6), 78–104. https://doi.org/10.1016/j.newar.2009.08.007