高次精度移流方程式解法 CIP 法

はじめに

流体力学分野において，東京工業大学名誉教授である矢部孝らによって提案された高次精度移流方程式解法である Constrained Interpolation Profile Method（CIP 法）*1を試してみた。

CIP 法は多次元の問題にも適用できるが，とりあえず，一次元の移流方程式に適用してみた。

方向分離解法を使えば，多次元化を一次元問題に帰着できる（M 型 CIP 法や C 型 CIP 法と呼ばれている）。

今後は，多次元の解法にも挑戦したいと思う。

CIP 法の特徴

CIP 法は，関数の値とその空間微分値を用いて，空間格子間を補間し，更新する手法である。

微分値を用いることで，波形を崩さずに更新できるという特長がある。

具体的には，関数値と微分値が，同じ速度で移流することを利用して，格子点間で 3 次関数補間（Hermite 補間）曲線を作成し，値を移流させる。

CIP法―原子から宇宙までを解くマルチスケール解法

作者:孝, 矢部,陽一, 尾形,隆行, 内海
発売日: 2003/11/01
メディア: 単行本

CIP法とJavaによるCGシミュレーション

作者:矢部孝,尾形陽一,滝沢研二
発売日: 2007/02/17
メディア: 単行本（ソフトカバー）

MATLAB ブログラム

下図に示すガウシアンパルスが速度 1 で伝搬していく様子を CIP 法によりシミュレーションする。

%===Parameter===%
NX=101;NT=100;
c=1;
dx=1;
CFL=0.5;
dt=CFL*dx/c;
x=0:dx:dx*(NX-1);
N=2;
sigma=5;
x0=50;
%===Memory Occupation===%
F=zeros(N,NX);
G=zeros(N,NX);
FN=zeros(N,NX);
GN=zeros(N,NX);
%===Initial Condition===%
F(1,:)=0.5*exp(-(x-x0).^2/2/sigma^2);
F(2,:)=0.5*exp(-(x-x0).^2/2/sigma^2);
G(1,:)=-(x-x0).*F(1,:)/sigma^2;
G(2,:)=-(x-x0).*F(2,:)/sigma^2;
%===Time Step===%
disp('Star Time Step !') 
for T=0:NT
    %===figure===%
    plot(x,F(1,:)+F(2,:),'b-');
    axis([0,max(x),-0.1,1.1]);
    xlabel('x');ylabel('F');
    grid on; drawnow;
    %===Step===%
    for n=1:N
        if n==1;
            D=-dx;xx=-c*dt;iup0=-1;
            for m=2:NX iup=m+iup0;
                a=(G(n,m)+G(n,iup))./D^2+2*(F(n,m)-F(n,iup))./D^3;
                b=3*(F(n,iup)-F(n,m))./D^2-(2*G(n,m)+G(n,iup))./D;
                FN(n,m)=a*xx^3+b*xx^2+G(n,m)*xx+F(n,m);
                GN(n,m)=3*a*xx^2+2*b*xx+G(n,m);
            end
        end
        if n==2;
            D= dx;xx= c*dt;iup0= 1;
            for m=NX-1:-1:1 iup=m+iup0;
                a=(G(n,m)+G(n,iup))./D^2+2*(F(n,m)-F(n,iup))./D^3;
                b=3*(F(n,iup)-F(n,m))./D^2-(2*G(n,m)+G(n,iup))./D;
                FN(n,m)=a*xx^3+b*xx^2+G(n,m)*xx+F(n,m);
                GN(n,m)=3*a*xx^2+2*b*xx+G(n,m);
            end
        end
    end
    %===Reset===%
    F=FN;G=GN;
end