移流方程式

自然界には波に関連する現象が多くあります．水面に伝わる波や携帯電話から出る電磁波、人の声(音)、など様々な波の形態で伝搬が行われます．この伝搬（何かしらの物理量が伝わっていくこと）を表す以下のような方程式を移流方程式と言います．

${ \displaystyle \dfrac {\partial f}{\partial t}+u\dfrac {\partial f}{\partial x}=0 }$

(f : 物理量, u : 流速)

導出

今回は流速が一定の場合を考えていきます．以下のようにf(x, t)のときの物理量はΔt秒後(次の時間でのサンプリング点)にf(x+cΔt, t+Δt)となります．

（※関数fはt,xの2変数関数なので、図に示す時は以下のように時間を決めたスナップショットのようになる．） f:id:takaishi78:20191130010547p:plain

以上より

${ \displaystyle f\left( x,t\right) =f\left( x+u\Delta t,t+\Delta t\right) }$

次の時間での物理量を求める形に変形 ${ \displaystyle f\left( x,t+\Delta t\right) =f\left( x-u\Delta t,t\right) }$

一次の項までテイラー展開を行うと ${ \displaystyle f\left( x,t\right) +\Delta t\dfrac {\partial f\left( x,t\right) }{\partial t}\simeq f\left( x,t\right) -u\Delta t\dfrac {\partial f\left( x,t\right) }{\partial x} }$

よって ${ \displaystyle \dfrac {\partial f}{\partial t}+u\dfrac {\partial f}{\partial x}=0 }$

初期条件（時刻 t = 0 のときの値）が与えられれば、この方程式の解（位置と時間から求まる物理量）を求めることができる

数値解法

移流方程式を数値解析的に解くためのいくつかのアルゴリズムをあげるてみると

1次風上差分
Lax-Wendroff法
CIP法

などがあげられる．

数値拡散

以下の画像のように波が進んでいってそれぞれの離散点だけをみて物理量の値を線形補間すると実際のグラフとの誤差が生まれる．このことを数値拡散という． (a)の黒い波が進むと青い波に位置に進む．

f:id:takaishi78:20200131194550j:plain — 数値拡散説明

1次風上差分

導出1

移流方程式 ${ \frac{\partial f} {\partial t}+ u\frac{\partial f}{\partial x} = 0 \tag{1.1} }$

時間微分に関しては前進差分を用いて、位置微分に関しては後退差分を用いる.位置微分に後退差分を用いることで出現するi-1が上流を示すため風上と呼ばれる．

${ \frac{f{_i}{^{n+1}} - f{_i}{^n}} {{\Delta t}} + u \frac{f{_i}{^{n}} - f{_{i-1}}{^n}} {{\Delta x}} = 0 \tag{1.2} }$

(2)式を変形すると

${ f{_i}{^{n+1}} = f{_i}{^n} + c(f{_{i-1}}{^n} - f{_i}{^n}) \tag{1.3} }$

cをクーラン数という

${ c = u\frac{\Delta t}{\Delta x} \tag{1.4} }$

導出2

黒い波が1ブロック分移動した際の新しい波を青で示している．

f:id:takaishi78:20200131212846j:plain — 風上差分-解説01

地点iにおいて、新しい値を算出したい場合、前の時間での波（黒）の波を元に考えることが出来る．下の図でピンクの直線の方程式は以下のようになる ${ F{_i} = \frac{f{_i} - f{_{i-1}}} {{\Delta x}} {(x - x_i)} + f{_i}{^n} \tag{2.1} }$

f:id:takaishi78:20200131213200j:plain — 風上差分-解説02

(2-1)式を元に地点iでの現在の値は先ほど求めた関数Fに移動した分だけ(今回で言うと1ブロック分)戻った位置の値と等しくなるので

${ F{_i}(x{_i} - u{\Delta t}) = f{^i}{_{n+1}} \tag{2.2} }$

${ f{_i}{^{n+1}} = f{_i}{^n} - \frac{u{\Delta t}}{\Delta x} (f{_i}{^n} - f{_{i-1}}{^n} ) \tag{2.3} }$

導出2の補足

例では微少時間として、移動分は1ブロックとしているがもし移動量が100ブロックなどになった場合上の方程式で得られる値はとても大きな値となり振動してしまう．これを防ぐためクーラン数は1以下でないとう条件が存在する．

CIP法

上記に示したように風上差分だと数値拡散(なまり)が発生してしまうので誤差がつもり、本来の波の形状からどんどん形が変形してしまう．この問題の解消の方法として2点をただ線形補完するのではなく多項式で補完する手法が提案されている．イメージとしてはこのような感じ．一次風上差分では青のような曲線になるところを多項式補完することで赤の曲線になる．

f:id:takaishi78:20200131222036j:plain — 多項式による近似イメージ図

CIP法では今までの位置だけでなく微分の情報を加えて伝搬させる．以下のようにサンプリング点での粒子の動き（緑の矢印）の情報をク加えて伝搬させてあげることで元の波形に近づけようとしている．

${ \frac{\partial f} {\partial t}+ u\frac{\partial f}{\partial x} = 0 \tag{3.1} }$ 移流方程式(3.1)を微分してもわかるように、微分した値も流速uによって情報が伝搬されている． ${ \frac{\partial g} {\partial t}+ u\frac{\partial g}{\partial x} = 0 \tag{3.2} }$

f:id:takaishi78:20200131223750j:plain — CIP-method01

補間用の曲線生成

CIP法は3次多項式によるエルミート補間をおこなうので、i-1, i点を３次多項式で表すと ${ F{_i}(x) = a{_i} (x - x{_i}) {^3} + b{_i} (x - x{_i}) {^2} + c{_i} (x - x{_i}) + d{_i} \tag{3.3} }$

次にi, i-1での値を求める(F{_i}(x)はi, i-1を元に作った曲線のこと)

${ F{_i}(x{_i}) = d{_i} = f{_i} \tag{3.4} }$

${ \frac {dF{_i}(x{_i})} {dx} = c{_i} = g{_i} \tag{3.5} }$

${ F{_i}(x{_{i-1}}) = - a{_i}{\Delta x}{^3} + b{_i}{\Delta x}{^2} - c{_i}{\Delta x}+ d{_i} = f{_{i-1}} \tag{3.6} }$

${ \frac {dF{_i}(x{_{x-1}})} {dx} = - 3a{_i}{\Delta x}{^2} + 2b{_i}{\Delta x} + c{_i} = g{_{i-1}} \tag{3.7} }$

以上の4つの関係式から係数a, bを求める ${ a{_i} = \frac { g{_i} +g{_{i-1}} }{ {\Delta x}{^2} } - \frac { 2(g{_i} -g{_{i-1}}) }{ {\Delta x}{^3} } \tag{3.8} }$

${ b{_i} = \frac { 3(g{_{i-1}} -g{_i}) }{ {\Delta x}{^2} } + \frac { 2g{_i} +g{_{i-1}} }{ {\Delta x} } \tag{3.9} }$

次時間n+1の算出

ここまでで補間方程式を求めることが出来た．これを用いて次時刻n+1での値を求めるにはこれまでのように ${u{\Delta t}}$ 進めた、 ${ f{^{n+1}} = F({x-u{\Delta t}}) }$ 、 ${ g{^{n+1}} = dF({x-u{\Delta t}}) }$ を求めれば良い．

${ f{^{n+1}}{_i}(x) = a{_i} (u{\Delta x}) {^3} + b{_i} (u{\Delta x}) {^2} + g{_i} (u{\Delta x}) + f{_i} \tag{3.10} }$

${ g{^{n+1}}{_i}(x) = 3a{_i} (u{\Delta x}) {^2} + 2b{_i} (u{\Delta x}) + g{_i} (u{\Delta x}) \tag{3.11} }$

サンプルコード

void ofApp::methodCIP(){
    int ip;
    double Fp, Gp;
    double c3, c2, x, s;
    
    //流速の方向を元に符号を決定
    if(speed > 0.0) {
        s = 1.0;
    }else{
        s = -1.0;
    }
    for(int i = 0; i <= nMesh; i++)
    {
        ip = i-(int)s;//上流点
        if(ip < 0.0 || ip > nMesh){
            Fp = 0.0;
            Gp = 0.0;
        }else{
            Fp = f0[ip];
            Gp = g0[ip];
        }
        
        double dx = -s * deltaX;
        double dx2 = dx * dx;
        double dx3 = dx2 * dx;
        
        c3 = (g0[i] + Gp) / dx2 + 2.0 * (f0[i] - Fp) / dx3;
        c2 = 3.0 * (Fp - f0[i]) / dx2 - (2.0 * g0[i] + Gp) / dx;
        x = -speed * deltaT;//uΔt
        
//こっちの方がわかりやすいかも
//f1[i] = ((c3*x+c2)*x+g0[i])*x + f1[i];
//g1[i] = (3.0 * c3 * x + 2.0 * c2) * x + g1[i];
        f1[i] += ((c3*x+c2)*x+g0[i])*x;
        g1[i] += (3.0 * c3 * x + 2.0 * c2) * x;
        
    }
//値を更新
    for(int i = 0; i <= nMesh; i++) {
        f0[i] = f1[i];
        g0[i] = g1[i];
    }
}