CIP法

CIP法では3次関数を補間関数として使用することで、双曲型問題に対して分散エラーが少ない、数値拡散が小さい、局所的な高精度補間ができる、等間隔格子を使う必要がないなどの利点が得られる。

右図で、左上の絵は移流の様子を表している。 これを格子点上の値として計算機に記憶させると右上の絵のようになる。 ここで、線形補間を行うと左下の絵のようになり、本当の波形であるピンク色の破線とはかなり開きが出てくる。 これが数値拡散であり、次の段階ではこの数値拡散によるなまりがさらに数値拡散を進展させることになる。 対して、右下の絵は傾きを考慮して補間を行っており、数値拡散が少ないことが分かる。

「CIP」とはもともとCubic Interpolated Pseudo-Particle Schemeの略であったが、研究がすすむにつれて3次多項式以外の補間関数を用いる手法へと発展した。 これにより、開発者の矢部孝は「CIP法」という名称の意味を考え直し、CIP法の本質が3次多項式を用いることにあるのではなく、元の方程式から導かれるいろいろな拘束条件をプロファイル(波の形状)に反映させることこそが本質であるとして現在の名称に変更した。よってCubic Interpolated Pseudo-Particle SchemeとConstrained Interpolation Profile Schemeのどちらも正式名称ということになる ^{[1] [2]}。

Summarize Timeline Top Qs Fact Check

CIP法は1次元移流方程式を高精度に解く解法である。1次元移流方程式は次式で与えられる。

${\displaystyle {\frac {\partial f}{\partial t}}+c{\frac {\partial f}{\partial x}}=0}$

ここで、cは移流速度である。

CIP法では、格子点の値${\displaystyle g(={\frac {\partial f}{\partial x}})}$についても同時に移流計算を行うことが特徴である。 空間微分値${\displaystyle g}$に対する移流方程式は、上の移流方程式を空間に関して微分することで得られ、以下のようになる。

${\displaystyle {\frac {\partial g}{\partial t}}+c{\frac {\partial g}{\partial x}}=0}$

時刻${\displaystyle n}$における値${\displaystyle f}$とその微分値${\displaystyle g}$が格子点上の点${\displaystyle i}$、${\displaystyle iup}$(点${\displaystyle iup}$は点${\displaystyle i}$の上流点、つまり移流速度 ${\displaystyle c>0}$なら${\displaystyle iup=i-1}$である)において既知とすると、この2点間のプロファイル(つまり形状)は以下のように3次多項式で表される。 ここで、上付き添字は時刻を、下付き添字は格子点番号をあらわす。

${\displaystyle F_{i}^{n}(x)=a_{i}(x-x_{i})^{3}+b_{i}(x-x_{i})^{2}+g_{i}^{n}(x-x_{i})+f_{i}^{n}}$

このようにプロファイルの補間関数を3次関数で表現することがCubicInterpolated Pseudo-Particle Schemeたる所以である。 ここで、係数${\displaystyle a_{i}}$、${\displaystyle b_{i}}$は、

${\displaystyle a_{i}={\frac {g_{i}^{n}+g_{iup}^{n}}{D^{2}}}+{\frac {2(f_{i}^{n}-f_{iup}^{n})}{D^{3}}}}$

${\displaystyle b_{i}={\frac {3(f_{iup}^{n}-f_{i}^{n})}{D^{2}}}-{\frac {2g_{i}^{n}+g_{iup}^{n}}{D}}}$

のようになる。 ただし、移流速度${\displaystyle c>0}$のとき${\displaystyle D=-\Delta x}$、${\displaystyle iup=i-1}$であり、移流速度${\displaystyle c<0}$のとき${\displaystyle D=\Delta x}$、${\displaystyle iup=i+1}$である。

適合条件式により、${\displaystyle F_{i}^{n}(0)=f_{i}^{n}}$、${\displaystyle {\frac {\partial F_{i}^{n}(0)}{\partial x}}=g_{i}^{n}}$、${\displaystyle F_{i}^{n}(D)=f_{iup}^{n}}$、${\displaystyle {\frac {\partial F_{i}^{n}(D)}{\partial x}}=g_{iup}^{n}}$が成り立つので、上式において係数${\displaystyle a_{i}}$、${\displaystyle b_{i}}$が求まる。

このように、格子点上の点において微分値${\displaystyle g}$も与えられるので、格子間のプロファイルを3次多項式で補間することができ、精度の高い計算が可能となる。

対象とする問題は移流方程式であるので、次の時刻${\displaystyle n+1}$での値${\displaystyle f_{i}^{n+1}}$と微分値${\displaystyle g_{i}^{n+1}}$は、この2点間のプロファイルを${\displaystyle c\Delta t}$だけ移動することで得られる。 つまり、${\displaystyle \xi =-c\Delta t}$として次式のようになる。

${\displaystyle f_{i}^{n+1}=F_{i}^{n}(\xi )=a_{i}\xi ^{3}+b_{i}\xi ^{2}+g_{i}^{n}\xi +f_{i}^{n}}$

${\displaystyle g_{i}^{n+1}={\frac {\partial F_{i}^{n}(\xi )}{\partial x}}=3a_{i}\xi ^{2}+2b_{i}\xi +g_{i}^{n}}$

Summarize Timeline Top Qs Fact Check

CIP法は多次元問題での移流方程式についても適用可能である。例として、2次元での移流方程式を考えてみるが、一般に${\displaystyle n}$次元での移流方程式にも適用できることを断っておく。

さて、2次元移流方程式は以下のように表せる。

${\displaystyle {\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}}}$

ここで、${\displaystyle {\boldsymbol {f}}}$は変数ベクトル、${\displaystyle {\boldsymbol {C}}_{\alpha }(\alpha =1,2)}$は係数マトリクスである。

2次元移流方程式にCIP法を適用する方法として、2元3次の多項式を補間関数として使う方法(A型CIP法)や方向分離解法により1次元移流方程式に落とし込んで計算する方法(M型CIP法)などが考えられる。

A型CIP法では、2次元移流方程式を解くにあたり、2元3次の多項式を補間関数として用いる。 つまり、${\displaystyle x_{1}}$軸、${\displaystyle x_{2}}$軸からなる平面空間において、補間関数は以下のようになる。

${\displaystyle F_{i,j}^{n}(x_{1},x_{2})=C_{3,0}(x_{1}-x_{1iup})^{3}+C_{2,0}(x_{1}-x_{1iup})^{2}+g_{xi,j}^{n}(x_{1}-x_{1iup})+f_{i,j}^{n}}$

${\displaystyle +C_{0,3}(x_{2}-x_{2jup})^{3}+C_{0,2}(x_{2}-x_{2jup})^{2}+g_{yi,j}^{n}(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}}$

ここでも、点${\displaystyle iup}$と点${\displaystyle jup}$はそれぞれ、点${\displaystyle i}$と点${\displaystyle j}$の上流点である。 また、${\displaystyle C}$は係数であり、1次元での場合と同様に適合条件式より 格子点${\displaystyle (i,j)}$、${\displaystyle (iup,j)}$、${\displaystyle (i,jup)}$の値${\displaystyle f^{n}}$と微分値${\displaystyle g^{n}}$、格子点${\displaystyle (iup,jup)}$での値${\displaystyle f_{iup,jup}^{n}}$を用いて求められる。

A型CIP法では、点${\displaystyle (iup,jup)}$において値${\displaystyle f^{n}}$の連続性しか要求していない。 しかし、他の3点${\displaystyle (i,j)}$、${\displaystyle (iup,j)}$、${\displaystyle (i,jup)}$では 値${\displaystyle f^{n}}$と微分値${\displaystyle g^{n}}$の連続性も保証している。 このため、求めようとしている点${\displaystyle (i,j)}$に対して対角線上にあり最も遠い点${\displaystyle (iup,jup)}$のプロファイルが 不正確であるために、このプロファイルを持ってくるような大きな時間ステップをとってはならない。

方向分離解法は一般に多次元問題を1次元問題に帰着させるために行われる。

A型CIP法では補間関数の係数の数が多く、これを解析に適用しようとすると格子点上で覚えさせる値の数が多くなり、計算を行う上で合理的でない。

M型CIP法では、多次元の移流方程式に方向分離を行うことで幾つかの1次元移流方程式に帰着させ、1次元のCIPスキームで計算を行う。 方向分離解法を適用することで、上の2次元移流方程式をつぎのように分解できる。

${\displaystyle {\frac {\partial {\boldsymbol {f}}}{\partial t}}+{\boldsymbol {C}}_{1}{\frac {\partial {\boldsymbol {f}}}{\partial x_{1}}}={\boldsymbol {0}}}$

${\displaystyle {\frac {\partial {\boldsymbol {f}}}{\partial t}}+{\boldsymbol {C}}_{2}{\frac {\partial {\boldsymbol {f}}}{\partial x_{2}}}={\boldsymbol {0}}}$

このように方向分離を行うと、${\displaystyle x_{1}}$方向へ分離した式を解くことによって時刻${\displaystyle n}$の値${\displaystyle {\boldsymbol {f}}^{n}}$から中間の値${\displaystyle {\boldsymbol {f}}^{*}}$が得られ、${\displaystyle x_{2}}$方向へ分離した式を解くことによって中間の値${\displaystyle {\boldsymbol {f}}^{*}}$から時刻${\displaystyle n+1}$の値${\displaystyle {\boldsymbol {f}}^{n+1}}$が得られる。

M型CIP法で${\displaystyle x_{\alpha }}$方向への移流計算を行うとき、ベクトル${\displaystyle {\boldsymbol {f}}}$と${\displaystyle x_{\alpha }}$方向の空間微分値${\displaystyle {\frac {\partial {\boldsymbol {f}}}{\partial x_{\alpha }}}}$については1節の1次元CIPスキームを使って解くことが出来るが、${\displaystyle x_{\beta }(\alpha \neq \beta )}$方向の空間微分値${\displaystyle {\frac {\partial {\boldsymbol {f}}}{\partial x_{\beta }}}}$については更なる空間微分値${\displaystyle {\frac {\partial ^{2}{\boldsymbol {f}}}{\partial x_{\alpha }\partial x_{\beta }}}}$を計算していないのでCIP法を使って求めることは出来ない。 よって${\displaystyle x_{\beta }}$方向の空間微分値${\displaystyle {\frac {\partial {\boldsymbol {f}}}{\partial x_{\beta }}}}$を求めるためには線形補間を行うことで、移流計算を行う。

M型CIP法では2階の空間微分値${\displaystyle {\frac {\partial ^{2}{\boldsymbol {f}}}{\partial x_{\alpha }\partial x_{\beta }}}(\alpha \neq \beta )}$を計算していなかったので、空間微分値${\displaystyle {\frac {\partial {\boldsymbol {f}}}{\partial x_{\beta }}}}$を計算する際は線形補間を行っていた。 この方法ではある程度の精度は保証されるが、格子間隔を広くとった場合などには${\displaystyle x_{\beta }}$方向の線形補間の影響が大きく出て、CIP法によるうまみを生かしきれなくなってしまう。

そこで、格子点上の2階空間微分値${\displaystyle {\frac {\partial ^{2}{\boldsymbol {f}}}{\partial x_{\alpha }\partial x_{\beta }}}}$を覚え、${\displaystyle x_{\beta }}$方向にもCIP計算を行おうというのがC型CIP法である。