CIP法のソースを表示

'''CIP法'''（CIPほう、{{lang-en-short|Constrained Interpolation Profile Scheme}}）とは、[[矢部孝]]らによって提案された、[[双曲型偏微分方程式]]を解く'''高次精度差分法'''の一つである。CIP法は高精度[[差分法|差分]]スキームであるので、[[機械工学]]、[[流体]]、[[電磁気]]、[[弾性体力学]]などの分野で広く[[数値解析]]に使用されている。

== 概要==

[[File:Cipfig.jpg|thumb|right|400px|数値拡散が起こる様子]]

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

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

「CIP」とはもともと{{Lang|en|Cubic Interpolated Pseudo-Particle Scheme}}の略であったが、研究がすすむにつれて3次多項式以外の補間関数を用いる手法へと発展した。
これにより、開発者の[[矢部孝]]は「CIP法」という名称の意味を考え直し、CIP法の本質が3次多項式を用いることにあるのではなく、元の方程式から導かれるいろいろな拘束条件をプロファイル(波の形状)に反映させることこそが本質であるとして現在の名称に変更した。よって{{Lang|en|Cubic Interpolated Pseudo-Particle Scheme}}と{{Lang|en|Constrained Interpolation Profile Scheme}}のどちらも正式名称ということになる
<ref>CFD最前線 p.209参照</ref>
<ref>https://w.atwiki.jp/mynote/?page=CIP%20Memo</ref>。

== 1次元移流方程式でのCIP解法 ==
[[File:Cp2lns.jpg|thumb|right|500px|CIP法計算スキームのイメージ]]

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

:<math> \frac{\partial f}{\partial t} +c\frac{\partial f}{\partial x}= 0</math>

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

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

:<math> \frac{\partial g}{\partial t} +c\frac{\partial g}{\partial x}= 0</math>

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

:<math> 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</math>


このようにプロファイルの補間関数を3次関数で表現することが{{Lang|en|CubicInterpolated Pseudo-Particle Scheme}}たる所以である。
ここで、係数<math>a_i</math>、<math>b_i</math>は、

:<math> a_i=\frac{g^n_i+g^n_{iup}}{D^2}+\frac{2(f^n_i-f^n_{iup})}{D^3}</math>
:<math> b_i=\frac{3(f^n_{iup}-f^n_{i})}{D^2}-\frac{2g^n_i+g^n_{iup}}{D}</math>

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


[[適合条件]]式により、<math>F^n_i(0)=f^n_i</math>、<math>\frac{\partial F^n_i(0)}{\partial x}=g^n_i</math>、<math>F^n_i(D)=f^n_{iup}</math>、<math>\frac{\partial F^n_i(D)}{\partial x}=g^n_{iup}</math>が成り立つので、上式において係数<math>a_i</math>、<math>b_i</math>が求まる。

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

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

:<math> f^{n+1}_i=F^n_i(\xi)=a_i\xi^3+b_i\xi^2+g^n_i\xi+f^n_i</math>
:<math> g^{n+1}_i=\frac{\partial F^n_i(\xi)}{\partial x}=3a_i\xi^2+2b_i\xi+g^n_i</math>

== 多次元問題でのCIP解法 ==

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


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

:<math> \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}</math>

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

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

=== A型CIP法 ===

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

:<math> 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}</math>
::<math> +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</math>

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


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

=== 方向分離解法 ===

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

=== M型CIP法 ===

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

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

:<math> \frac{\partial \boldsymbol{f}}{\partial t}+\boldsymbol{C}_{1}\frac{\partial \boldsymbol{f}}{\partial x_{1}} =\boldsymbol{0}</math>
:<math> \frac{\partial \boldsymbol{f}}{\partial t}+\boldsymbol{C}_{2}\frac{\partial \boldsymbol{f}}{\partial x_{2}} =\boldsymbol{0}</math>

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


M型CIP法で<math>x_{\alpha}</math>方向への移流計算を行うとき、ベクトル<math>\boldsymbol{f}</math>と<math>x_{\alpha}</math>方向の空間微分値<math>\frac{\partial \boldsymbol{f}}{\partial x_{\alpha}}</math>については1節の1次元CIPスキームを使って解くことが出来るが、<math>x_{\beta}(\alpha \neq \beta)</math>方向の空間微分値<math>\frac{\partial \boldsymbol{f}}{\partial x_{\beta}}</math>については更なる空間微分値<math>\frac{\partial^{2} \boldsymbol{f}}{\partial x_{\alpha}\partial x_{\beta}}</math>を計算していないのでCIP法を使って求めることは出来ない。
よって<math>x_{\beta}</math>方向の空間微分値<math>\frac{\partial \boldsymbol{f}}{\partial x_{\beta}}</math>を求めるためには線形補間を行うことで、移流計算を行う。
<div style="font-size:larger;">
:<math title="M型CIP法計算の流れ">
\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}{|c|}
			\hline
			x_\alpha \ \text{方} \ \text{向} \ \text{へ} \ \text{の} \ \text{移} \ \text{流} \ \text{計} \ \text{算} \ \\
			\hline
		\end{array}
		&
		\begin{array}{|c|}
			\hline
			x_\beta \ \text{方} \ \text{向} \ \text{へ} \ \text{の} \ \text{移} \ \text{流} \ \text{計} \ \text{算} \ \\
			\hline
		\end{array}
	\end{array}
\end{align}
</math>
</div>

=== C型CIP法 ===

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


そこで、格子点上の2階空間微分値<math>\frac{\partial^{2} \boldsymbol{f}}{\partial x_{\alpha}\partial x_{\beta}}</math>を覚え、<math>x_{\beta}</math>方向にもCIP計算を行おうというのがC型CIP法である。
<div style="font-size:larger;">
:<math title="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}{|c|}
			\hline
			x_\alpha \ \text{方} \ \text{向} \ \text{へ} \ \text{の} \ \text{移} \ \text{流} \ \text{計} \ \text{算} \ \\
			\hline
		\end{array}
		&
		\begin{array}{|c|}
			\hline
			x_\beta \ \text{方} \ \text{向} \ \text{へ} \ \text{の} \ \text{移} \ \text{流} \ \text{計} \ \text{算} \ \\
			\hline
		\end{array}
	\end{array}
\end{align}
</math>
</div>

==脚注==
{{Reflist}}

== 参考文献 ==
*{{Cite book|和書
|author=矢部孝、内海隆行、尾形陽一
|year=2003-11-07
|title=CIP法：原子から宇宙までを解くマルチスケール解法
|publisher=森北出版
|isbn=978-4-627-91831-3
}}
* 矢部孝、尾形陽一、滝沢研二：「CIP法とJavaによるCGシミュレーション」、森北出版、ISBN 978-4-627-91911-2、(2007年2月17日）。
*{{Cite book|和書
|author=蔦原道久、渡利實、棚橋隆彦、矢部孝
|year=2007-04-25
|title=CFD最前線
|publisher=共立出版
|isbn=978-4-320-08166-6
}}
* 日本応用数理学会「応用数理」、Vol.18, No.2 (2008), (CIP法特集号）。 
* 肖鋒、小野寺直幸、伊井仁志：「計算流体力学：CIPマルチモーメント法による手法」、コロナ社、ISBN 978-4-339-04597-0（2009年3月26日）。
* 日本建築学会 (編)：「はじめての音響数値シミュレーションプログラミングガイド」（第6章：CIP(constrained interpolation profile)法）、コロナ社、ISBN 978-4-339-00838-8（2012年11月30日）。

== 関連項目 ==
* [[移流]]
* [[差分法]]
* [[移流方程式]]
* [[方向分離解法]]

{{DEFAULTSORT:しいあいひいほう}}
[[Category:応用力学]]