CIP法
CIP法(CIPほう、Constrained Interpolation Profile Scheme)とは、矢部孝らによって提案された、双曲型偏微分方程式を解く高次精度差分法の一つである。CIP法は高精度差分スキームであるので、機械工学、流体、電磁気、弾性体力学などの分野で広く数値解析に使用されている。
CIP法では3次関数を補間関数として使用することで、双曲型問題に対して分散エラーが少ない、数値拡散が小さい、局所的な高精度補間ができる、等間隔格子を使う必要がないなどの利点が得られる。
右図で、左上の絵は移流の様子を表している。
これを格子点上の値として計算機に記憶させると右上の絵のようになる。
ここで、線形補間を行うと左下の絵のようになり、本当の波形であるピンク色の破線とはかなり開きが出てくる。
これが数値拡散であり、次の段階ではこの数値拡散によるなまりがさらに数値拡散を進展させることになる。
対して、右下の絵は傾きを考慮して補間を行っており、数値拡散が少ないことが分かる。
「CIP」とはもともとCubic Interpolated Pseudo-Particle Schemeの略であったが、研究がすすむにつれて3次多項式以外の補間関数を用いる手法へと発展した。
これにより、開発者の矢部孝は「CIP法」という名称の意味を考え直し、CIP法の本質が3次多項式を用いることにあるのではなく、元の方程式から導かれるいろいろな拘束条件をプロファイル(波の形状)に反映させることこそが本質であるとして現在の名称に変更した。よってCubic Interpolated Pseudo-Particle SchemeとConstrained Interpolation Profile Schemeのどちらも正式名称ということになる
。
CIP法は1次元移流方程式を高精度に解く解法である。1次元移流方程式は次式で与えられる。
\frac{\partial f}{\partial t} +c\frac{\partial f}{\partial x}= 0
ここで、cは移流速度である。
CIP法では、格子点の値g(=\frac{\partial f}{\partial x})についても同時に移流計算を行うことが特徴である。
空間微分値gに対する移流方程式は、上の移流方程式を空間に関して微分することで得られ、以下のようになる。
\frac{\partial g}{\partial t} +c\frac{\partial g}{\partial x}= 0
時刻nにおける値fとその微分値gが格子点上の点i、iup(点iupは点iの上流点、つまり移流速度
c>0ならiup=i-1である)において既知とすると、この2点間のプロファイル(つまり形状)は以下のように3次多項式で表される。
ここで、上付き添字は時刻を、下付き添字は格子点番号をあらわす。
F^n_i(x)=a_i(x-x_{i})^3+b_i(x-x_{i})^2+g^n_i(x-x_{i})+f^n_i
このようにプロファイルの補間関数を3次関数で表現することがCubicInterpolated Pseudo-Particle Schemeたる所以である。
ここで、係数a_i、b_iは、
a_i=\frac{g^n_i+g^n_{iup}{D^2}+\frac{2(f^n_i-f^n_{iup})}{D^3}
b_i=\frac{3(f^n_{iup}-f^n_{i})}{D^2}-\frac{2g^n_i+g^n_{iup}{D}
のようになる。
ただし、移流速度c>0のときD=-\Delta x、iup=i-1であり、移流速度c<0のときD=\Delta x、iup=i+1である。
適合条件式により、F^n_i(0)=f^n_i、\frac{\partial F^n_i(0)}{\partial x}=g^n_i、F^n_i(D)=f^n_{iup}、\frac{\partial F^n_i(D)}{\partial x}=g^n_{iup}が成り立つので、上式において係数a_i、b_iが求まる。
このように、格子点上の点において微分値gも与えられるので、格子間のプロファイルを3次多項式で補間することができ、精度の高い計算が可能となる。
対象とする問題は移流方程式であるので、次の時刻n+1での値f^{n+1}_iと微分値g^{n+1}_iは、この2点間のプロファイルをc\Delta tだけ移動することで得られる。
つまり、\xi=-c\Delta tとして次式のようになる。
f^{n+1}_i=F^n_i(\xi)=a_i\xi^3+b_i\xi^2+g^n_i\xi+f^n_i
g^{n+1}_i=\frac{\partial F^n_i(\xi)}{\partial x}=3a_i\xi^2+2b_i\xi+g^n_i
CIP法は多次元問題での移流方程式についても適用可能である。例として、2次元での移流方程式を考えてみるが、一般にn次元での移流方程式にも適用できることを断っておく。
さて、2次元移流方程式は以下のように表せる。
\frac{\partial \boldsymbol{f}{\partial t}+\boldsymbol{C}_{1}\frac{\partial \boldsymbol{f}{\partial x_{1} + \boldsymbol{C}_{2}\frac{\partial \boldsymbol{f}{\partial x_{2}=\boldsymbol{0}
ここで、\boldsymbol{f}は変数ベクトル、\boldsymbol{C}_{\alpha}(\alpha=1,2)は係数マトリクスである。
2次元移流方程式にCIP法を適用する方法として、2元3次の多項式を補間関数として使う方法(A型CIP法)や方向分離解法により1次元移流方程式に落とし込んで計算する方法(M型CIP法)などが考えられる。
A型CIP法では、2次元移流方程式を解くにあたり、2元3次の多項式を補間関数として用いる。
つまり、x_{1}軸、x_{2}軸からなる平面空間において、補間関数は以下のようになる。
F^n_{i,j}(x_{1},x_{2})=C_{3,0}(x_{1}-x_{1iup})^3+C_{2,0}(x_{1}-x_{1iup})^2+g^n_{xi,j}(x_{1}-x_{1iup})+f^n_{i,j}
+C_{0,3}(x_{2}-x_{2jup})^3+C_{0,2}(x_{2}-x_{2jup})^2+g^n_{yi,j}(x_{2}-x_{2jup})+C_{2,1}(x_{1}-x_{1iup})^2(x_{2}-x_{2jup})+C_{1,1}(x_{1}-x_{1iup})(x_{2}-x_{2jup})+C_{1,2}(x_{1}-x_{1iup})(x_{2}-x_{2jup})^2
ここでも、点iupと点jupはそれぞれ、点iと点jの上流点である。
また、Cは係数であり、1次元での場合と同様に適合条件式より
格子点(i,j)、(iup,j)、(i,jup)の値f^nと微分値g^n、格子点(iup,jup)での値f^n_{iup,jup}を用いて求められる。
A型CIP法では、点(iup,jup)において値f^nの連続性しか要求していない。
しかし、他の3点(i,j)、(iup,j)、(i,jup)では
値f^nと微分値g^nの連続性も保証している。
このため、求めようとしている点(i,j)に対して対角線上にあり最も遠い点(iup,jup)のプロファイルが
不正確であるために、このプロファイルを持ってくるような大きな時間ステップをとってはならない。
方向分離解法は一般に多次元問題を1次元問題に帰着させるために行われる。
A型CIP法では補間関数の係数の数が多く、これを解析に適用しようとすると格子点上で覚えさせる値の数が多くなり、計算を行う上で合理的でない。
M型CIP法では、多次元の移流方程式に方向分離を行うことで幾つかの1次元移流方程式に帰着させ、1次元のCIPスキームで計算を行う。
方向分離解法を適用することで、上の2次元移流方程式をつぎのように分解できる。
\frac{\partial \boldsymbol{f}{\partial t}+\boldsymbol{C}_{1}\frac{\partial \boldsymbol{f}{\partial x_{1} =\boldsymbol{0}
\frac{\partial \boldsymbol{f}{\partial t}+\boldsymbol{C}_{2}\frac{\partial \boldsymbol{f}{\partial x_{2} =\boldsymbol{0}
このように方向分離を行うと、x_{1}方向へ分離した式を解くことによって時刻nの値\boldsymbol{f}^{n}から中間の値\boldsymbol{f}^{*}が得られ、x_{2}方向へ分離した式を解くことによって中間の値\boldsymbol{f}^{*}から時刻n+1の値\boldsymbol{f}^{n+1}が得られる。
M型CIP法でx_{\alpha}方向への移流計算を行うとき、ベクトル\boldsymbol{f}とx_{\alpha}方向の空間微分値\frac{\partial \boldsymbol{f}{\partial x_{\alpha}については1節の1次元CIPスキームを使って解くことが出来るが、x_{\beta}(\alpha \neq \beta)方向の空間微分値\frac{\partial \boldsymbol{f}{\partial x_{\beta}については更なる空間微分値\frac{\partial^{2} \boldsymbol{f}{\partial x_{\alpha}\partial x_{\beta}を計算していないのでCIP法を使って求めることは出来ない。
よってx_{\beta}方向の空間微分値\frac{\partial \boldsymbol{f}{\partial x_{\beta}を求めるためには線形補間を行うことで、移流計算を行う。
\begin{align}
& \qquad \mathsf{M} \ \text{型} \ \ \mathsf{CIP} \ \text{法} \ \text{の} \ \text{計} \ \text{算} \ \text{順} \ \text{序} \ \ \left(n \rightarrow n + 1 \right) \\
\\
& \begin{array}{cc}
\begin{array}{c}
\sigma_{ij}^n \\
\\
\frac{\partial\sigma_{ij}^n}{\partial x_\alpha} \\
\\
\frac{\partial\sigma_{ij}^n}{\partial x_\beta}
\end{array}
\begin{array}{c}
\xrightarrow{\mathsf{CIP} \ \text{法} \\ \\ \\ \\
\xrightarrow{\mathsf{FDM} \\
\end{array}
\begin{array}{c}
\sigma_{ij}^* \\
\\
\frac{\partial\sigma_{ij}^*}{\partial x_\alpha} \\
\\
\frac{\partial\sigma_{ij}^*}{\partial x_\beta}
\end{array}
&
\begin{array}{c}
\sigma_{ij}^* \\
\\
\frac{\partial\sigma_{ij}^*}{\partial x_\beta} \\
\\
\frac{\partial\sigma_{ij}^*}{\partial x_\alpha}
\end{array}
\begin{array}{c}
\xrightarrow{\mathsf{CIP} \ \text{法} \\ \\ \\ \\
\xrightarrow{\mathsf{FDM} \\
\end{array}
\begin{array}{c}
\sigma_{ij}^{n+1} \\
\\
\frac{\partial\sigma_{ij}^{n+1}{\partial x_\beta} \\
\\
\frac{\partial\sigma_{ij}^{n+1}{\partial x_\alpha}
\end{array}
\\
\begin{array}
\hline
x_\alpha \ \text{方} \ \text{向} \ \text{へ} \ \text{の} \ \text{移} \ \text{流} \ \text{計} \ \text{算} \ \\
\hline
\end{array}
&
\begin{array}
\hline
x_\beta \ \text{方} \ \text{向} \ \text{へ} \ \text{の} \ \text{移} \ \text{流} \ \text{計} \ \text{算} \ \\
\hline
\end{array}
\end{array}
\end{align}
M型CIP法では2階の空間微分値\frac{\partial^{2} \boldsymbol{f}{\partial x_{\alpha}\partial x_{\beta}(\alpha \neq \beta)を計算していなかったので、空間微分値\frac{\partial \boldsymbol{f}{\partial x_{\beta}を計算する際は線形補間を行っていた。
この方法ではある程度の精度は保証されるが、格子間隔を広くとった場合などにはx_{\beta}方向の線形補間の影響が大きく出て、CIP法によるうまみを生かしきれなくなってしまう。
そこで、格子点上の2階空間微分値\frac{\partial^{2} \boldsymbol{f}{\partial x_{\alpha}\partial x_{\beta}を覚え、x_{\beta}方向にもCIP計算を行おうというのがC型CIP法である。
\begin{align}
& \qquad \mathsf{C} \ \text{型} \ \ \mathsf{CIP} \ \text{法} \ \text{の} \ \text{計} \ \text{算} \ \text{順} \ \text{序} \ \ \left(n \rightarrow n + 1 \right) \\
\\
& \begin{array}{cc}
\begin{array}{c}
\sigma_{ij}^n \\
\\
\frac{\partial\sigma_{ij}^n}{\partial x_\alpha} \\
\\
\frac{\partial\sigma_{ij}^n}{\partial x_\beta} \\
\frac{\partial^2\sigma_{ij}^n}{\partial x_\alpha \partial x_\beta}
\end{array}
\begin{array}{c}
\xrightarrow{\mathsf{CIP} \ \text{法} \\ \\ \\ \\ \\
\xrightarrow{\mathsf{CIP} \ \text{法} \\ \\
\end{array}
\begin{array}{c}
\sigma_{ij}^* \\
\\
\frac{\partial\sigma_{ij}^*}{\partial x_\alpha} \\
\\
\frac{\partial\sigma_{ij}^*}{\partial x_\beta} \\
\frac{\partial^2\sigma_{ij}^*}{\partial x_\alpha \partial x_\beta}
\end{array}
&
\begin{array}{c}
\sigma_{ij}^* \\
\\
\frac{\partial\sigma_{ij}^*}{\partial x_\beta} \\
\\
\frac{\partial\sigma_{ij}^*}{\partial x_\alpha} \\
\frac{\partial^2\sigma_{ij}^*}{\partial x_\alpha \partial x_\beta}
\end{array}
\begin{array}{c}
\xrightarrow{\mathsf{CIP} \ \text{法} \\ \\ \\ \\ \\
\xrightarrow{\mathsf{CIP} \ \text{法} \\ \\
\end{array}
\begin{array}{c}
\sigma_{ij}^{n+1} \\
\\
\frac{\partial\sigma_{ij}^{n+1}{\partial x_\beta} \\
\\
\frac{\partial\sigma_{ij}^{n+1}{\partial x_\alpha} \\
\frac{\partial^2\sigma_{ij}^{n+1}{\partial x_\alpha \partial x_\beta}
\end{array}
\\
\begin{array}
\hline
x_\alpha \ \text{方} \ \text{向} \ \text{へ} \ \text{の} \ \text{移} \ \text{流} \ \text{計} \ \text{算} \ \\
\hline
\end{array}
&
\begin{array}
\hline
x_\beta \ \text{方} \ \text{向} \ \text{へ} \ \text{の} \ \text{移} \ \text{流} \ \text{計} \ \text{算} \ \\
\hline
\end{array}
\end{array}
\end{align}