工学系大学院生のブログ

2人の大学院生による雑記ブログ

第9-4回 陰解法(LU-SGS)[python, julia]




\begin{eqnarray}
\frac{\Delta Q^n}{\Delta t} + \frac{\partial}{\partial x} \left(\frac{\partial E}{\partial Q} {\Delta Q}\right)^n = R^n, \ R^n = -\left(\frac{\partial E}{\partial x} \right)^n
\end{eqnarray}


第9-1回では陰解法の場合の式変形を行いました.


第9-4回以降ではLU-SGSという手法を用いて、Euler方程式を解きたいと思います.



こちらの手法は現在でも使用されている方法であり,少し複雑ですが噛み砕いて紹介したいと思います.




LU-SGS(Lower-Upper Symmetric Gauss Seidel)のポイントは行列の反転の際に計算量が少なくなるよう工夫が凝らされている点になります。




一番最初に示した式の左辺第二項について離散化を行うのですが,LU-SGSではヤコビアン行列を近似します.




まず,ヤコビアン行列を\(A\)と置き,中心差分を用いてセル\(i\)の離散化を行います.
\begin{eqnarray}
\frac{\partial}{\partial x} \left(A {\Delta Q}\right)^n = \frac{1}{\Delta x} \left((A\Delta Q)^n_{i+1/2}-(A\Delta Q)^n_{i-1/2})\right)
\end{eqnarray}




この第一項はFVSのようにヤコビアン行列の風向きを考えて近似します.

FVSについては過去の記事をご覧ください。



概念的には次のようになります.



\begin{eqnarray}
(A\Delta Q)^n_{i+1/2} = (A^+ \Delta Q)^n_{i} + (A^- \Delta Q)^n_{i+1}
\end{eqnarray}


(+)は固有値が正のもの,(-)は固有値が負のものを示しています.




まとめます.
\begin{eqnarray}
\left[ -\frac{1}{\Delta x}\left(A^+ \right)^n_{i-1} \ \ \ \frac{1}{\Delta t}\textbf{I}+\frac{1}{\Delta x}\left(A^+ \right)^n_{i}-\frac{1}{\Delta x}\left(A^- \right)^n_{i} \ \ \ \frac{1}{\Delta x}\left(A^- \right)^n_{i+1} \right] \left[ \begin{array}{c} \Delta Q^n_{i-1} \\ \Delta Q^n_i \\ \Delta Q^n_{i+1} \end{array}\right]= R^n_i
\end{eqnarray}




FVSで使用した固有値が正負のヤコビアン行列の計算手法を簡単に紹介します.
\begin{eqnarray}
A^+ = R\left(\frac{\Lambda+|\Lambda|}{2}\right) R^{-1} = \frac{1}{2}(A+R|\Lambda|R^{-1}) \\
A^- = R\left(\frac{\Lambda-|\Lambda|}{2}\right) R^{-1} = \frac{1}{2}(A-R|\Lambda|R^{-1})
\end{eqnarray}




LU-SGSでは,それぞれの第二項一つの波で近似します.後で理由を述べますが一番大きい値で近似します.


その値を\(\sigma\)とします.

\begin{eqnarray}
\textbf{A}^+ &=&\frac{1}{2}(\textbf{A}+\textbf{R}|\Lambda |\textbf{R}^{-1}) \simeq \frac{1}{2}(\textbf{A}+\textbf{R}\sigma \textbf{R}^{-1}) \\
&=& \frac{1}{2}(\textbf{A}+\sigma \textbf{I}) \\
\textbf{A}^- &=& \frac{1}{2}(\textbf{A}-\sigma \textbf{I})
\end{eqnarray}




\(I\)は対角成分が1でそれ以外が0の対角行列です.
\begin{eqnarray}
\textbf{I} = \left[\begin{array}{ccc} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \ \end{array}\right]
\end{eqnarray}




この結果を代入すると
\begin{eqnarray}
\left[ -\frac{1}{\Delta x}\left(A^+ \right)^n_{i-1} \ \ \ \frac{1}{\Delta t}\textbf{I}+\frac{1}{\Delta x} \sigma^n_{i}\textbf{I} \ \ \ \frac{1}{\Delta x}\left(A^- \right)^n_{i+1} \right] \left[ \begin{array}{c} \Delta Q^n_{i-1} \\ \Delta Q^n_i \\ \Delta Q^n_{i+1} \end{array}\right]= R^n_i
\end{eqnarray}




これで,対角行列の逆行列は下記の様に逆数になるため,簡単に求めることができます.
\begin{eqnarray}
\textbf{D} = \left[\begin{array}{ccc} a & 0 & 0 \\ 0 & a & 0 \\ 0 & 0 & a \end{array}\right]
\Leftrightarrow \textbf{D}^{-1} = \left[\begin{array}{ccc} 1/a & 0 & 0 \\ 0 & 1/a & 0 \\ 0 & 0 & 1/a \ \end{array}\right]
\end{eqnarray}



また,第9-2回で紹介したLDU分解を用いるガウスザイデル法の形に上式を適用すると反転の際に逆行列を求める必要がなく,計算速度が向上します.



少し大雑把ではありますが,このような手法がLU-SGSになります.




ところで,近似についての議論を飛ばしましたが,3つの固有値を1つに近似しても大丈夫なのでしょうか?


そんな大胆な近似がまかり通るなら,そもそも解く方程式は1つで十分ということになります.



もちろん、そんなわけありません.つまり,この近似のまま計算しても結果は上手く行きません.




なので,次回の第9-5回はこの問題を解決すべく内部反復法に触れたいと思います.

この内部反復法は,LU-SGSに限らず陰解法には幅広く用いられている手法になります.



第9-6回にLU-SGSのコードがあります。


by hide


・参考論文
A. Jameson, S. Yoon, “Lower-Upper Implicit Schemes with Multiple Grids for the Euler Equations,” AIAA Jounal,25(7):929-935. (1987).

このエントリーをはてなブックマークに追加

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です