第9-3回　陰解法(Euler方程式、中心差分コード)[python, julia]

\begin{eqnarray}
\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}

\begin{eqnarray}
R^n_i = -\left(\frac{E^n_{i+1/2}-E^n_{i-1/2}}{\Delta x} \right)^n
\end{eqnarray}

第9-3回では、オイラー方程式に対し陰解法を用いてsodの衝撃波管に対する解析を行います。

第6回では、muscl法を用いてsodの衝撃波管の計算を行っており、そのコードを元に書き直したので、muscl法等は第6回をご覧ください。

また、反復法の収束条件として、残差の2-ノルムが\(10^{-4}\)以下のとき収束したとみなします。
残差や2-ノルムについては、下記の記事の最後の方に書いてありますのでご参照ください。

解析条件は下記の通りです。

時間ステップ数 : \(t=300\)
時間刻み幅 : \(\Delta t=0.002\)
空間分割数 : \(nx=100\)
空間刻み幅 : \(\Delta x=0.01\)

コードはgithubにアップしております。（pythonとjulia)

https://github.com/hide-dog/Euler_implicit_central

解析結果に移ります。（青：厳密解、赤：解析解）

見ての通り、厳密解と一致しています。

陰解法によって生じるヤコビアン行列に対し、中心差分で評価しましたが、上記のように計算することができました。

しかし、陰解法の強みは、陽解法と比べ時間刻み幅を長くとれることです。

それでは5倍長い\(\Delta t = 0.01\)として、計算してみましょう。

結果は「発散」でした。

陰解法は安定なはずなのに、どうしてでしょうか？

\begin{eqnarray}
\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}

改めて今回解いた行列式を書きました。

ガウスザイデル法の収束条件は対角優位です。つまり、対角成分がほかの成分と比べて大きくなければいけません。

仮に\(\Delta t\)を大きくすると、対角成分は徐々に小さくなってしまい、対角優位を満たさなくなります。このことから、陰解法にも関わらず、\(\Delta t\)を自由にとることができなかったのです。

では、現在でも使用されているLU-SGSという手法はどのようになっているのでしょうか。

第9-4回からLU-SGSについて掘り下げていきたいと思います。

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

今回使用したコードは下記にアップしております。

https://github.com/hide-dog/Euler_implicit_central

by hide