第9-6回　LU-SGSと内部反復(Euler方程式)[python,julia]

\begin{eqnarray}
\left[ -\frac{1}{\Delta x}\left(A^+ \right)^n_{i-1} \ \ \ \frac{1}{\Delta \tau}+\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^m_{i-1} \\ \Delta Q^m_i \\ \Delta Q^m_{i+1} \end{array}\right]= R^n_i
\end{eqnarray}

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

第9-6回では、LU-SGSという手法を用いてsodの衝撃波管を解きたいと思います。

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

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

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

https://github.com/hide-dog/Euler_implicit_LU-SGS

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

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

それでは、解析結果になります。（青：厳密解、赤：解析解）

これまでの解析結果と同じく、しっかりと解けていることがわかります。

また、時間刻み幅の長いケース\(\Delta t =0.01 \)も計算してみました。


第9-3回で行った中心差分の陰解法では上手くいきませんでしたが、今回は発散することなく下記の通りになりました。（青：厳密解、赤：解析解）

こちらも解析解と一致しており、内部反復とLU-SGSを用いることで時間刻み幅が大きい場合も解くことができました。

考察として、時間刻み幅を大きくとった場合は解析解とのずれが大きくなっています。

これは、時間項に対しては1次精度前進差分をとっていること、LDU分解の際に時間刻み幅の二乗の項を無視したことがあげられます。

安定的に解析できるとはいえ、刻み幅が大きいと精度が落ちるので注意してください。

ここまで読んでくださり、ありがとうございます。

実用的に使用されているコードに、徐々に近づいておりどんどん楽しくなっています。


by hide

・補足１ \(\beta\)について

\(\frac{\partial}{\partial x}\frac{\partial E}{\partial Q}\)に対し中心差分を用いた場合は、対角優位になりにくいため、反復法が適用しにくいことが問題です。

そのため、LU-SGSでは対角優位にするため、定数\(\beta\)をかけています。こちらは経験的に\(1.1～1.01\)程度の値にするようです。


内部反復を用いているため、このように定数をかけても収束さえすれば問題ありません。