第9-5回　内部反復[python, julia]

\begin{eqnarray}
\left[ -\frac{1}{\Delta x}\left(A^+ \right)^n_{i-1} \ \ \ \frac{1}{\Delta t}+\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}

前回はLU-SGSを用いて陰解法を進めてきました。しかし、この手法は近似をしすぎており、このままでは上手くいきません。

そこで今回は、陰解法で幅広く使用されている内部反復という手法を紹介します。


いくつか文字を定義していきますが、最終的にはなぜそのように文字を置くのかわかるようにしますので、読んでみてください。

まず、下記のように\(Q^m\)を定義します。ここで、\(m\)は疑似時間ステップ、\(n\)は実際の物理時間ステップです。

\begin{eqnarray}
\lim_{m \to \infty} Q^m = Q^{n+1}
\end{eqnarray}

この\(Q^m\)を用いて、今までの式に新しく第一項目を足します。この項は疑似時間項、\(\Delta \tau\)は疑似時間幅になります。

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

この式に対し\(m\)の極限をとった場合は、一項目がキャンセルされ、第二項が元の式の形と一致します。

では、これのうまみを理解するため陰解法を適用して式変形してみます。

気を付けてほしいのですが、最終的に\(Q^m = Q^{n+1}\)となるため、上記の式の\(m\)について時間を進めていきます。

\begin{eqnarray}
\frac{\Delta Q^m}{\Delta \tau} +\frac{\partial}{\partial x} \left( \frac{\partial E}{\partial Q} \Delta Q^m \right) = R^m \ \ \ \ \ \ \ \ \ ・・・(1)
\end{eqnarray}

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

今までの陰解法と比べて、右辺に第一項の分が増えていますが、それ以外は変わっていません。この左辺にLU-SGSを用いて近似していました。

では\(m\)についての時間ステップを進めるとどうなるでしょうか？

十分に時間を進めると、\(Q^m = Q^{n+1}, \ Q^{m+1} = Q^{n+1}\)となるため、\(\Delta Q^m = 0\)となります。

すなわち

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

以上の結果から、陰解法を適用した(1)式の左辺に対しどのような近似を用いたとしても、疑似時間を進め収束した場合は元の方程式を満たします。

これが内部反復のメリットになります。

デメリットとしては、近似が雑過ぎれば収束しないことや内部反復ステップを進める必要があるため、物理時間の1ステップの計算量が増えてしまうなどがあげられます。

実際のコードに落とし込む際は、初めに\(Q^m = Q^n\)として、\(Q^m\)の疑似時間ステップを進めます。そして、十分に\(Q^m\)が収束したとき、\(Q^m\)は\(Q^{n+1}\)に近づいたということなので、\(Q^m = Q^{n+1}\)として物理時間ステップを進めます。

また、疑似時間刻み幅\(\Delta \tau\)として、今回は物理時間刻み幅\(\Delta t\)と同じものを使いたいと思います。

そして、疑似時間ステップ数は経験測として5~50程度でとっているイメージがあります。しっかりと行うのであれば、どのくらいのステップを進めれば収束するか評価する必要があります。

第9-6回ではLUSGSと内部反復をコードに落とし込みたいと思います。

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

by hide