第9-2回　保存形の陰解法(中心差分)[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-2回では左辺第二項に中心差分を適用した場合について進めます．

はじめに断っておきますが，こちらの手法は解を導けますが，少し問題がありますのでご注意ください．問題点については第9-3回で述べます。

セル\(i\)における式を考えて中心差分を用いて離散化します．
\begin{eqnarray}
\frac{\Delta Q^n_i}{\Delta t} + \frac{1}{2\Delta x}\left( \left(\frac{\partial E}{\partial Q} {\Delta Q}\right)^n_{i+1}-\left(\frac{\partial E}{\partial Q} {\Delta Q}\right)^n_{i-1}\right) = R^n_i
\end{eqnarray}

\begin{eqnarray}
\Leftrightarrow \left[ -\frac{1}{2\Delta x}\left(\frac{\partial E}{\partial Q} \right)^n_{i-1} \ \ \ \frac{1}{\Delta t} \ \ \ \frac{1}{2\Delta x}\left(\frac{\partial E}{\partial Q} \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}

気を付けてほしいのですが，この式は\(i\)セルにおける式なので，同じような式がセル分存在しているわけです．それをまとめるとこのようになります．
\begin{eqnarray}
\left[\begin{array}{ccccc} \ddots& \ddots & \ddots && \\ & -\frac{1}{2\Delta x}\left(\frac{\partial E}{\partial Q} \right)^n_{i-1} & \frac{1}{\Delta t} & \frac{1}{2\Delta x}\left(\frac{\partial E}{\partial Q} \right)^n_{i+1} & \\ && \ddots & \ddots & \ddots \\ \end{array}\right] \left[\begin{array}{c} \vdots \\ \Delta Q_{i-1}^n \\ \Delta Q_{i}^n \\ \Delta Q_{i+1}^n \\ \vdots \end{array}\right] =\left[\begin{array}{c} \vdots \\ \vdots \\ R_i^n \\ \vdots \\ \vdots \end{array}\right]
\end{eqnarray}

左の行列を反転させることで解\(\Delta Q\)を求めることができます．

この反転にはガウスザイデル法を使用します．
ガウスザイデル法は逐次的に更新しながら反復していく方法になります．

拡散方程式に対するガウスザイデル法の紹介は過去にしていたので，よかったら見てみてください．

ここでは，上記の内容と少し違う「LDU分解を用いたガウスザイデル法」を使用します．

まず行列の左辺を次のように変形，近似を行います．この近似については，最後に詳しく述べます．
\begin{eqnarray}
(\textbf{L}+\textbf{D}+\textbf{U})\Delta\textbf{Q} &=& (\textbf{L}+\textbf{D})\textbf{D}^{-1}(\textbf{D}+\textbf{U})\Delta\textbf{Q}-\textbf{L}\textbf{D}^{-1}\textbf{U}\Delta\textbf{Q} \\
&\simeq& (\textbf{L}+\textbf{D})\textbf{D}^{-1}(\textbf{D}+\textbf{U})\Delta\textbf{Q}
\end{eqnarray}

ここで，\(\textbf{L},\textbf{D},\textbf{U}\)はそれぞれ下三角行列成分，対角行列成分，上対角行列成分になります．

この近似を用いて下記のように二段階で方程式を解きます．
\begin{eqnarray}
(\textbf{L}+\textbf{D})\Delta \textbf{Q}^* = \textbf{R}
\end{eqnarray}
\begin{eqnarray}
(\textbf{D}+\textbf{U})\Delta \textbf{Q} = \textbf{D}\Delta\textbf{Q}^*
\end{eqnarray}

二つの式を繋ぐために\(\Delta \textbf{Q}^*\)を用いていますが，上式から\(\Delta \textbf{Q}^*\)を消去すれば元の式に戻ることが確認できると思います．

これがLDU分解を用いたガウスザイデル法の式になります．気になる点は要素ごとに計算して頂ければ，ガウスザイデル法と一致することがわかるかと思います．

第9-3回では、これらを用いてsodの衝撃波管の問題を解いて見たいと思います．

＊第9-3回に中心差分の場合のコードがあります。

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

by hide

・LDU分解のオーダーチェック
\begin{eqnarray}
\left[\begin{array}{ccccc} \ddots& \ddots & \ddots && \\ & -\frac{1}{2\Delta x}\left(\frac{\partial E}{\partial Q} \right)^n_{i-1} & \frac{1}{\Delta t} & \frac{1}{2\Delta x}\left(\frac{\partial E}{\partial Q} \right)^n_{i+1} & \\ && \ddots & \ddots & \ddots \\ \end{array}\right] \left[\begin{array}{c} \vdots \\ \Delta Q_{i-1}^n \\ \Delta Q_{i}^n \\ \Delta Q_{i+1}^n \\ \vdots \end{array}\right] =\left[\begin{array}{c} \vdots \\ \vdots \\ R_i^n \\ \vdots \\ \vdots \end{array}\right]
\end{eqnarray}

このような式になっているので，時間項に関するオーダーはそれぞれ次のようになります．

\begin{eqnarray}
\textbf{L}=O(1) , \ \textbf{D}=O(1/\Delta t) , \ \textbf{U}=O(1)
\end{eqnarray}

LDU分解で生じる二つの項のオーダーは次のようになります．
\begin{eqnarray}
(\textbf{L}+\textbf{D}) \textbf{D}^{-1} (\textbf{D}+\textbf{U}) &=& O(1/ \Delta t)O(\Delta t)O(1/ \Delta t) = O(1/ \Delta t) \\
\textbf{L}\textbf{D}^{-1} \textbf{U} &=& O(\Delta t)
\end{eqnarray}

よって第二項は，第一項に対して\(\Delta t^2\)のオーダーであるため，無視できます．