圧縮性流体の運動方程式(運動量保存則)

（一般の）圧縮性流体の運動方程式は、テンソルが出てきて分かりにくいので、レイノルズの輸送定理を用いた具体的な導出をメモします。

レイノルズの輸送定理
$$ \int_{B_t}f({\bm{x}},t) dv $$
空間座標でtの変化とともに積分領域が変化する時に上の物質微分をしたいとすると、積分領域が変化するのでそのままでは上手くいかない。
一度物質座標に変換（変形勾配テンソルを使う）すれば、物質座標では積分領域が変化しない。なぜなら、tに合わせて空間座標における物質座標も変化するからである。（ゴムをのばすとゴム上の座標ものびるイメージ）
なので、物質座標で積分してから、また空間座標に戻せば良い。すると以下の結果が得られる。
$$
\frac{D}{Dt}\int_{B_t}f({\bm{x},t}) dv= \int_{B_t}(\dot{f}({\bm{x},t})+f({\bm{x},t})\mbox{div}{\bm{v}})dv
$$
$ f $の物質微分を考えて、
$$
\frac{D}{Dt}\int_{B_t}f({\bm{x},t}) dv= \int_{B_t}(\frac{\partial f}{\partial t}+{\bm{v}}\cdot\mbox{grad}f+f({\bm{x},t})\mbox{div}{\bm{v}})dv $$
よって
$$ \frac{D}{Dt}\int_{B_t}f({\bm{x},t}) dv=\int_{B_t}(\frac{\partial f}{\partial t}+\mbox{div}(f{\bm{v}})) dv
$$

運動量変化は以下のように表される
$$ \frac{D}{Dt}\int_{V}\rho{\bm{v}}dv $$
速さがベクトルだと扱いづらいので$ {\bm{e_k}} $成分について考える。レイノルズの輸送定理をつかうと
\begin{equation} \frac{D}{Dt}\int_{V}\rho{v_k}dv=\int_V(\frac{\partial \rho v_k}{\partial t}+\frac{\partial \rho v_k v_j}{\partial x_j})dv\end{equation}

また、$$ \nabla\cdot(\rho {\bm{v}}\otimes{\bm{v}})= ({\bm{e_i}}\frac{\partial }{\partial x_i})\cdot (\rho v_j v_k{\bm{e_j}}\otimes{\bm{e_k}})$$
$$=\frac{\partial \rho v_j v_k}{\partial x_j}{\bm{e}}_k$$
は(1)の右辺第二項と一致する
よって圧縮性流体の運動方程式の右辺は
$$
\frac{\partial\rho{\bm{u}}}{\partial t}+\nabla\cdot(\rho{\bm{v}}\otimes{\bm{v}})=\mbox{外部からの力積}
$$

共役勾配法(CG法)の収束について

工学の様々な場面で登場する前処理付き共役勾配法。
共役勾配法は前処理をするのが重要(PCG法)とよく言われますが、その理由は省略されていることが多いのでメモしておきます。

$ A $は正定値対称行列
$$ {\bm{r}}_0={\bm{b}}-A{\bm{x}}_0 $$
$$ {\bm{p}}_{i+1}={\bm{r}}_{i+1}+\beta_i{\bm{p}}_i $$
$$ {\bm{x}}_{i+1}={\bm{x}}_i+\alpha_i{\bm{p}}_i $$
$$ {\bm{r}}_{i+1}={\bm{b}}-A{\bm{x}}_{i+1} ={\bm{b}}-A({\bm{x}}_i+\alpha_i{\bm{p}}_i )={\bm{r}}_i-\alpha_iA{\bm{p}}_i $$
$$ A{\bm{x^*}}={\bm{b}} $$
$$ {\bm{e}}_i={\bm{x^{*}}}-{\bm{x}}_i $$
$$ {\bm{e}}_{i+1}={\bm{{e}}}_i-\alpha_i{\bm{p}}_i $$
ベクトルの太字は面倒なので以下省略。
クリロフ部分空間
$$ \mathcal{K}_i=\mbox{span}\{ r_0,Ar_0,\cdots A^{i-1}r_0\}= \mbox{span}\{ Ae_0,A^2e_0,\cdots A^{i}e_0\}$$
$ i $次多項式の多項式$ P_i $を用いて(polynomial)
$$ e_i=P_i(A)e_0 $$
と書ける。アルゴリズムを書いてみて動かすとわかるが、$ P_i(0)=1 $
ここからが重要
$ A $の固有値、固有ベクトルを$ \lambda_i,u_i $とする。
$ e_0 $は以下のように表せる。$ a_i $は係数。
$$ e_0=\sum_{j=1}^n a_i u_i $$
$$ e_i=P_i(A)\sum_{j=1}^n a_j u_j =\sum_{j=1}^n a_jP_i(\lambda_j)u_j$$
$$ Ae_i=P_i(A)A\sum_{j=1}^n a_j u_j =\sum_{j=1}^n a_j\lambda_jP_i(\lambda_j)u_j$$
対称行列の固有ベクトルは直交するので
$$ (e_i,Ae_i)=\sum_{j=1}^n {a_j}^2 P_i(\lambda_j)^2\lambda_j $$
$ A $は正定値対称なので、$ (x-x^*,A(x-x^*)) $を最小にするのは$ x=x^* $のみである。
共役勾配法は各ステップで$ (e_{i+1},Ae_{i+1}) $が最小になるように$ \alpha_i $を選んでいるので、
$ \mathcal{K}_i $上で$ (e_{i},Ae_{i}) $が最小になるから、
$$ (e_i,Ae_i)=\min_{P_i}\sum_{j=1}^n {a_j}^2 P_i(\lambda_j)^2\lambda_j $$
ここで、$ 0\leq X_i,Y_i $とし、$ X_{max}=\max_{i }X_i $とすると、一般に$ x^x $
$$　\sum_i X_i Y_i\leq\sum_i X_{max} Y_i=X_{max}\sum_i Y_i $$
が成り立つ。Aは正定値なので$ X_i=P_i(\lambda_i)^2,Y_i=a_j^2\lambda_j $とすると、
\begin{eqnarray}
\sum_{j=1}^n {a_j}^2 P_i(\lambda_j)^2\lambda_j&\leq  &
\max_{\lambda_k\in{\mbox{Aの固有値の集合}}}P_i(\lambda_k)^2\sum_{j=1}^n {a_j}^2 \lambda_j\\
&\leq&\max_{\lambda\in[\lambda_{min},\lambda_{max}]}P_i(\lambda)^2\sum_{j=1}^n {a_j}^2 \lambda_j
\end{eqnarray}
よって、
\begin{eqnarray}
(e_i,Ae_i)&=&\min_{P_i}\sum_{j=1}^n {a_j}^2 P_i(\lambda)^2\lambda_j\\
&\leq&\min_{P_i}\max_{\lambda\in[\lambda_{min},\lambda_{max}]}P_i(\lambda_j)^2\sum_{j=1}^n {a_j}^2 \lambda_j\\
&=&\min_{P_i}\max_{\lambda\in[\lambda_{min},\lambda_{max}]}P_i(\lambda)^2(e_0,Ae_0)
\end{eqnarray}
上式は$ \lambda $がmin~maxまで動く時の関数の最大値を最小にするような$ P_i $を選んだ時の値を表している
最小偏差多項式を使えば収束の速さが評価できる。
$ k $次のチェビシェフ多項式を以下のように書ける
$$ T_k(x)=\cos(l\arccos x) $$
$$ \cfrac{\lambda_{min}}{\lambda_{max}}
$$
が大きいほど早いことがわかる