移流項の差分近似

初めに

最近，有限差分法の基礎的な勉強をしていた．

これまで，「数値粘性」という言葉は耳にしながらも，其の実，何者か良く分かっていなかった．

少しだけ，自分なりに解釈できた気がするので纏めてみる（理解を深化して補足すべき点が未だ多いが，，，）．

移流項の差分近似

移流方程式を考える（非保存型）．

\frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x} = 0

ここで，$\phi \left( = \phi(t, x); (t, x) \in (0, T) \times \Omega \right)$は解（例えば，物質の濃度），$u$は移流の速度である．このような式を考えるのは，例えば，流体の基礎方程式であるNavier-Stokes方程式にも移流項が登場し，その取り扱いを注意深く考えなければならないため，まずは簡略化した式を考えよう，というマインドである．

1次精度風上差分スキーム

移流項を差分近似する．最も単純な方法として，1次精度風上差分スキームがある．

u \frac{\partial \phi}{\partial x} \approx 
    \begin{cases}
      u \frac{\phi_{i} - \phi_{i-1}}{\Delta x} & \text{if $u \ge 0$} \\
      u \frac{\phi_{i+1} - \phi_{i}}{\Delta x} & \text{if $u \lt 0$}
    \end{cases}

場合分けは，移流が正の方向 / 負の方向どちら向きかによって定まる（本当は図で説明したいが，今は図を作るよりも文字に残しておきたいので，後々アップデートすることにする）．

さて，移流の向きによって場合分けが，，，と記したが，

u = \frac{u + |u|}{2} + \frac{u - |u|}{2}

という関係から，上式を以下のように纏めることができる．

\begin{align}
u \frac{\partial \phi}{\partial x} 
    &\approx
    \frac{u + |u|}{2} \frac{\phi_{i} - \phi_{i-1}}{\Delta x}
     + \frac{u - |u|}{2} \frac{\phi_{i+1} - \phi_{i}}{\Delta x} \\
    &= \frac{u}{2} \frac{(\phi_{i} - \phi_{i-1}) + (\phi_{i+1} - \phi_{i})}{\Delta x} 
    + \frac{|u|}{2} \frac{(\phi_{i} - \phi_{i-1}) + (- \phi_{i+1} + \phi_{i})}{\Delta x} \\
    &= u \frac{\phi_{i+1} - \phi_{i-1}}{2 \Delta x} 
    - \frac{|u| \Delta x}{2} \frac{\phi_{i+1} - 2 \phi_{i} + \phi_{i-1}}{(\Delta x)^2}
\end{align}

ここで，

• 第1項は移流項を2次精度中心差分近似したもの
• 第2項は拡散項を2次精度中心差分近似し，何らかの係数（$|u| \Delta x / 2$）を乗じたもの

になっている．本来，移流方程式を離散化したつもりだったのに，移流「拡散」方程式を離散化したんだっけ？と思ってしまうような形になった．

移流項を差分近似した結果，新たに登場してしまった第2項が，「数値粘性」と呼ばれるものである．これは分子粘性（教科書で$\nu$などの文字で登場する物理的な粘性）とは異なり，本来興味のある現象の，粘性の影響を変えてしまう存在である（今回，分子粘性はゼロだが，このスキームでは，粘性が$|u| \Delta x / 2$だけ過大評価されてしまう）．

実際，分子粘性$\nu$の移流拡散方程式：

\frac{\partial \phi}{\partial t} + u \frac{\partial \phi}{\partial x} = \nu \frac{\partial^2 \phi}{\partial x^2}

の移流項に1次精度風上差分，拡散項に2次精度中心差分を適用すると，

\frac{\partial \phi}{\partial t} + u \frac{\phi_{i+1} - \phi_{i-1}}{2 \Delta x} = \left( \nu + \frac{|u| \Delta x}{2} \right) \frac{\phi_{i+1} - 2 \phi_{i} + \phi_{i-1}}{(\Delta x)^2}

となってしまい，分子粘性$\nu$に数値粘性$\bar{\nu} = |u| \Delta x / 2$が加えられている．

既に十分かもしれないが，再度，シンプルな移流方程式に立ち戻り，これを1次精度風上差分でシミュレートする．初期条件に，

\phi =
    \begin{cases}
      1 & \text{if $0.1 \le x \le 0.3$} \\
      0 & \text{elsewhere}
    \end{cases}

を与え，移流速度$u = 0.5$とする．時間方向には前進Euler法を適用する．空間$\Omega:=[0, 1]$を$\Delta x = 0.01$で離散化し，時間増分$\Delta t = 0.1 \Delta x / u$で$T = 1$までシミュレートする（Courant数: $0.1$）．

シミュレーション結果から，初期条件として与えられた矩形波が，高さを失い，裾を広げながら移流していることが確認できる（起きるはずのない拡散が発生している）．本来興味のある現象とは異なる現象をシミュレートしてしまったことになるので，これは大きな問題である（本来追跡したい現象を追跡できていない）．

数値粘性が悪いのならば，

u \frac{\partial \phi}{\partial x} 
    \approx
    u \frac{\phi_{i+1} - \phi_{i-1}}{2 \Delta x} 

としてしまい，2次精度中心差分だけで近似すれば良いのでは？と感じてしまう．実際に計算してみる．

今回は数値解が激しく振動してしまい，非物理的な結果を得た．ここで，Godunovの定理というものが知られており，その主張は，極めて簡略化して言ってしまえば「数値振動を起こさない離散化スキームは，高々1次精度である」というものである（と私は解釈している）¹．上に示したシミュレーションは2次精度であり，振動してしまうスキームであることが分かる（かつ，振動の程度がとても激しい）．

（2024/06追記）
「上記は誤りである」とは未だ言い辛いが（今の知識では自信を持って判断できない．．．），より明解な説明を得た．von Neumannの安定性解析を実行すれば，これが不安定なスキームであることが分かる．

数値粘性が入り込んでしまうのは困る，かと思いきや，数値振動が発生しても困る．実際，数値粘性が入っているお陰で振動が抑えられ，安定にシミュレーションが実行できているのである．ただし，やはり数値粘性が強すぎると，追従したい現象でないものをシミュレートしてしまうので，少々の数値粘性には目を瞑るものの，これはなるべく小さくしたいというのが本心である．ここから，数値粘性を小さく留めつつ，数値振動も抑えて安定にシミュレートする方法が提案されてきた．

Lax-Wendroffの方法

Lax-Wendroffスキーム（Lax&Wendroff1960）は，1次精度風上差分スキームの次に名前が挙がる，最も有名な方法の一つ．Taylor展開からシンプルに導出することができる（導出には複数の方法があるらしいが，初めに勉強したこれがしっくり来たので他の方法は良く知らない）．

解$\phi$を時間方向に離散化し，$\phi (t + \Delta t)$を$\phi (t)$まわりでTaylor展開する．

\phi^{(n+1)}_{i} = \phi^{(n)}_{i} 
    + \Delta t \left. \frac{\partial \phi}{\partial t} \right|^{(n)}_{i}
    + \frac{(\Delta t)^2}{2} \left. \frac{\partial^2 \phi}{\partial t^2} \right|^{(n)}_{i}
    + \mathcal{O} \left[ (\Delta t)^3 \right]

ここで，元の移流方程式から，

\begin{align}
    \frac{\partial \phi}{\partial t} &= - u \frac{\partial \phi}{\partial x} \\
    \frac{\partial^2 \phi}{\partial t^2}
    &= - u \frac{\partial}{\partial t} \left(\frac{\partial \phi}{\partial x}\right)
    = u^2 \frac{\partial^2 \phi}{\partial x^2}
\end{align}

が得られる．これより，上式を以下のように改める．

\phi^{(n+1)}_{i} = \phi^{(n)}_{i} 
    - \Delta t \left. u \frac{\partial \phi}{\partial x} \right|^{(n)}_{i}
    + \frac{(\Delta t)^2}{2} \left. u^2 \frac{\partial^2 \phi}{\partial x^2} \right|^{(n)}_{i}

既に陽解法になっている．シミュレーションにはこのまま使えば良い．

他のスキームとの比較のため，移流項の近似に着目したいなら，元の移流方程式に戻せば良い．すなわち，

\frac{\phi^{(n+1)}_{i} - \phi^{(n)}_{i}}{\Delta t}
+ u \left. \frac{\partial \phi}{\partial x} \right|^{(n)}_{i}
- \frac{|u|^2 \Delta t}{2} \left. \frac{\partial^2 \phi}{\partial x^2} \right|^{(n)}_{i}
= 0 

と整理すれば，Lax-Wendroffスキームにおける移流項の近似は

\begin{align}
    u \frac{\partial \phi}{\partial x}
    &\approx u \left. \frac{\partial \phi}{\partial x} \right|_{i}
    - \frac{|u|^2 \Delta t}{2} \left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i} \\
    &\approx u \frac{\phi_{i+1} - \phi_{i-1}}{2 \Delta x}
    - \frac{|u|^2 \Delta t}{2} \frac{\phi_{i+1} - 2 \phi_{i} + \phi_{i-1}}{(\Delta x)^2}
\end{align}

となっていることが分かる．

これは，先に記した1次精度風上差分スキームに良く似ている．具体的には，数値粘性の大きさが異なるのみであり，

• 1次精度風上差分では，$|u| \Delta x / 2 \left( = \bar{\nu}_{\mathrm{FU}} \right)$
• Lax-Wendroffでは，$|u|^2 \Delta t / 2 \left( = \bar{\nu}_{\mathrm{LW}} \right)$

だけ，粘性を過大評価する（$\tilde{\nu} = \nu + \bar{\nu}$，$\nu$は分子粘性，$\bar{\nu}$は数値粘性）．ここで，時間増分$\Delta t$には空間解像度$\Delta x$より小さい値を採用することが一般的である．なぜなら，（陽解法において）移流のシミュレーションにはCFL条件が課されるためである²．

C := \sum_{j}^{d} \frac{|u_{j}| \Delta t}{\Delta x_{j}} \le C_{\mathrm{max}}

ここで，$j$は空間の軸方向，$d$は問題の空間次元（$\Omega \subset \mathbb{R}^d$）．$C_{\mathrm{max}} = 1$として整理すると，以下を得る．

\Delta t \le \left( \sum_{j}^{d} \frac{|u_{j}|}{\Delta x_{j}} \right)^{-1}

さらに，1次元の問題であるとして簡略化すると，

\begin{align}
    \Delta t &\le \frac{\Delta x}{|u|} \\
    |u|^2 \Delta t &\le |u|^2 \frac{\Delta x}{|u|} \\
    \therefore \bar{\nu}_{\mathrm{LW}} &\le \bar{\nu}_{\mathrm{FU}} \\
\end{align}

なる関係を得るから，Lax-Wendroffは1次精度風上差分より小さい数値粘性に留まると期待できる．上記の矩形波の移流をシミュレートすると，以下の結果を得る．

1次精度風上差分と比較してLax-Wendroffは矩形波前方の立ち上がりをよりシャープに捉えている．また，高さを失い裾を広げる傾向は改善されているようである．しかしながら，矩形波後方では振動が残されており，さらに高度な手法を模索したい気持ちになる．

Leaonardの方法

Leaonard（Leaonard1979）は，

• QUICKスキーム
• QUICKESTスキーム

の2つの手法を提案している．

QUICK (Quadratic Upstream Interpolation for Convective Kinematics)

1次精度風上差分，およびLax-Wendroffは3点ステンシルを用いたスキームであった．i.e., 注目している節点$\phi_{i}$の前後1点ずつ（$\phi_{i+1}$，$\phi_{i-1}$）の物理量で以て移流を表現する手法であった．

LeaonardのQUICKスキームは，もう少しだけ腕を伸ばして5点ステンシルの計算を行う．まずは，$\phi_{i+1/2}$を$\phi_{i}$まわりでTaylor級数展開するところから始める．

\begin{split}
    \phi_{i+1/2}
    &= \phi_{i}
    + \frac{\Delta x}{2} \left. \frac{\partial \phi}{\partial x} \right|_{i}
    + \frac{1}{2} \left( \frac{\Delta x}{2} \right)^{2} \left. \frac{\partial^{2} \phi}{\partial x^{2}} \right|_{i}
    + \mathcal{O} \left[ (\Delta x)^{3} \right] \\
    &= \phi_{i}
    + \frac{\Delta x}{2} \frac{\phi_{i+1} - \phi_{i-1}}{2 \Delta x}
    + \frac{1}{2} \left( \frac{\Delta x}{2} \right)^{2} \frac{\phi_{i+1} - 2 \phi_{i} + \phi_{i-1}}{(\Delta x)^2}
    + \mathcal{O} \left[ (\Delta x)^{3} \right] \\
    &= \frac{3 \phi_{i+1} + 6 \phi_{i} - \phi_{i-1}}{8}
    + \mathcal{O} \left[ (\Delta x)^{3} \right]
\end{split}

続いて，$\phi_{i-1/2}$を$\phi_{i-1}$まわりでTaylor展開する．

\begin{align}
    \phi_{i-1/2}
    &= \phi_{i-1}
    + \frac{\Delta x}{2} \left. \frac{\partial \phi}{\partial x} \right|_{i-1}
    + \frac{1}{2} \left( \frac{\Delta x}{2} \right)^{2} \left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i-1}
    + \mathcal{O} \left[ (\Delta x)^{3} \right] \\
    &= \phi_{i-1}
    + \frac{\Delta x}{2} \frac{\phi_{i} - \phi_{i-2}}{2 \Delta x}
    + \frac{1}{2} \left( \frac{\Delta x}{2} \right)^{2} \frac{\phi_{i} - 2 \phi_{i-1} + \phi_{i-2}}{(\Delta x)^2}
    + \mathcal{O} \left[ (\Delta x)^{3} \right] \\
    &= \frac{3 \phi_{i} + 6 \phi_{i-1} - \phi_{i-2}}{8}
    + \mathcal{O} \left[ (\Delta x)^{3} \right]
\end{align}

以上より，

\begin{align}
    \phi_{i\pm1/2}
    &= \phi_{i} 
    \pm \frac{1}{2} \Delta x \left. \frac{\partial \phi}{\partial x} \right|_{i}
    + \frac{1}{2} \left( \frac{1}{2} \Delta x \right)^2 \left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i}
    + \mathcal{O} \left[ (\Delta x)^3 \right] \\
    \phi_{i+1/2} - \phi_{i-1/2}
    &= \Delta x \left. \frac{\partial \phi}{\partial x} \right|_{i}
    + \mathcal{O} \left[ (\Delta x)^3 \right] \\
    \left. \frac{\partial \phi}{\partial x} \right|_{i}
    &= \frac{3 \phi_{i+1} + 3 \phi_{i} - 7 \phi_{i-1} + \phi_{i-2}}{8 \Delta x} + \mathcal{O} \left[ (\Delta x)^2 \right]
\end{align}

したがって，移流の方向に応じて以下のように近似を与える．

u \frac{\partial \phi}{\partial x} \approx 
    \begin{cases}
      u \frac{3 \phi_{i+1} + 3 \phi_{i} - 7 \phi_{i-1} + \phi_{i-2}}{8 \Delta x} & \text{if $u \ge 0$} \\
      u \frac{- \phi_{i+2} + 7 \phi_{i+1} - 3 \phi_{i} - 3 \phi_{i-1}}{8 \Delta x} & \text{if $u \lt 0$}
    \end{cases}

また，再度$u = (u + |u|) / 2 + (u - |u|) / 2$の関係を用いて式を整理する．

\begin{align}
u \frac{\partial \phi}{\partial x} 
    &\approx \frac{u + |u|}{2} \frac{3 \phi_{i+1} + 3 \phi_{i} - 7 \phi_{i-1} + \phi_{i-2}}{8 \Delta x}
    + \frac{u - |u|}{2} \frac{- \phi_{i+2} + 7 \phi_{i+1} - 3 \phi_{i} - 3 \phi_{i-1}}{8 \Delta x} \\
    &= u \frac{- \phi_{i+2} + 10 \phi_{i+1} - 10 \phi_{i-1} + \phi_{i-2}}{16 \Delta x}
    + \frac{|u| (\Delta x)^3}{16} \frac{\phi_{i+2} - 4 \phi_{i+1} + 6 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{(\Delta x)^4}
\end{align}

ここで，

• 第1項は移流項を4次精度中心差分近似したもの
• 第2項は4次の空間微分を中心差分近似し，何らかの係数を乗じたもの

である．上記の2手法との相違点は，これまでは数値粘性が2次の空間微分であり分子粘性（物理的な拡散）と同じオーダーであったのに対し，QUICKでは4次の空間微分となっており，より局所化された拡散効果を持つという点にある³．

QUICKEST (QUICK with Estimated Stream Terms)

QUICKESTスキームは，QUICKと同じ論文の中でLeonardによって導入されている．Lax-Wendroffと同様に，Taylor展開からスタートする：

\phi^{(n+1)}_{i} = \phi^{(n)}_{i} 
    + \Delta t \left. \frac{\partial \phi}{\partial t} \right|^{(n)}_{i}
    + \frac{(\Delta t)^2}{2} \left. \frac{\partial^2 \phi}{\partial t^2} \right|^{(n)}_{i}
    + \frac{(\Delta t)^3}{6} \left. \frac{\partial^3 \phi}{\partial t^3} \right|^{(n)}_{i}
    + \mathcal{O} \left[ (\Delta t)^4 \right]

また，元の移流方程式の時間方向への1階，および2階の微分を求めておく．

\begin{align}
    \frac{\partial \phi}{\partial t} &= - u \frac{\partial \phi}{\partial x} \\
    \frac{\partial^2 \phi}{\partial t^2} &= \frac{1}{\partial t} \left(  - u \frac{\partial \phi}{\partial x} \right) = u^2 \frac{\partial^2 \phi}{\partial x^2} \\
    \frac{\partial^3 \phi}{\partial t^3} &= \frac{1}{\partial t} \left( u^2 \frac{\partial^2 \phi}{\partial x^2} \right) = - u^3 \frac{\partial^3 \phi}{\partial x^3}
\end{align}

先のTaylor展開に代入する．

\phi^{(n+1)}_{i} = \phi^{(n)}_{i} 
    - \Delta t \left. u \frac{\partial \phi}{\partial x} \right|^{(n)}_{i}
    + \frac{(\Delta t)^2}{2} \left. u^2 \frac{\partial^2 \phi}{\partial x^2} \right|^{(n)}_{i}
    - \frac{(\Delta t)^3}{6} \left. u^3 \frac{\partial^3 \phi}{\partial x^3} \right|^{(n)}_{i}

以前と同じように既に陽解法になっているので，適切に空間離散化できれば直ちに解ける．

元の移流方程式の形に戻せば，以下が移流項の離散化であると分かる．

u \frac{\partial \phi}{\partial x} \approx
    \left. u \frac{\partial \phi}{\partial x} \right|_{i}
    - \frac{\Delta t}{2} \left. u^2 \frac{\partial^2 \phi}{\partial x^2} \right|_{i}
    + \frac{(\Delta t)^2}{6} \left. u^3 \frac{\partial^3 \phi}{\partial x^3} \right|_{i}

ここで，2階・3階の微分には中心差分を適用する．

\begin{align}
\left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i}
&\approx \frac{\phi_{i+1} - 2 \phi_{i} + \phi_{i-1}}{(\Delta x)^2} \\
\left. \frac{\partial^3 \phi}{\partial x^3} \right|_{i}
&\approx \frac{\phi_{i+2} - 2 \phi_{i+1} + 2 \phi_{i-1} - \phi_{i-2}}{2 (\Delta x)^2}
\end{align}

また，1階の微分は以下のように離散化する．

\left. u \frac{\partial \phi}{\partial x} \right|_{i} \approx 
    \begin{cases}
      u \frac{2 \phi_{i+1} + 3 \phi_{i} - 6 \phi_{i-1} + \phi_{i-2}}{6 \Delta x} & \text{if $u \ge 0$} \\
      u \frac{- \phi_{i+2} + 6 \phi_{i+1} - 3 \phi_{i} - 2 \phi_{i-1}}{6 \Delta x} & \text{if $u \lt 0$}
    \end{cases}

より，

\begin{align}
u \frac{\partial \phi}{\partial x} 
  &\approx u \frac{- \phi_{i+2} + 8 \phi_{i+1} - 8 \phi_{i-1} + \phi_{i-2}}{12 \Delta x}
  + \frac{|u| (\Delta x)^3}{12} \frac{\phi_{i+2} - 4 \phi_{i+1} + 6 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{(\Delta x)^4}
\end{align}

である．以上より，

\begin{align}
u \frac{\partial \phi}{\partial x} 
& \approx
    \left. u \frac{\partial \phi}{\partial x} \right|^{(n)}_{i}
    - \frac{\Delta t}{2} \left. u^2 \frac{\partial^2 \phi}{\partial x^2} \right|^{(n)}_{i}
    + \frac{(\Delta t)^2}{6} \left. u^3 \frac{\partial^3 \phi}{\partial x^3} \right|^{(n)}_{i} \\
& \approx
    u \frac{- \phi_{i+2} + 8 \phi_{i+1} - 8 \phi_{i-1} + \phi_{i-2}}{12 \Delta x} \\
    &+ \frac{|u| (\Delta x)^3}{12} \frac{\phi_{i+2} - 4 \phi_{i+1} + 6 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{(\Delta x)^4} \\
    &- \frac{\Delta t u^2}{2} \frac{\phi_{i+1} - 2 \phi_{i} + \phi_{i-1}}{(\Delta x)^2} \\
    &+ \frac{(\Delta t)^2 u^3}{6} \frac{\phi_{i+2} - 2 \phi_{i+1} + 2 \phi_{i-1} - \phi_{i-2}}{2 (\Delta x)^2}
\end{align}

を得る．

これまで通り，QUICK，およびQUICKESTによる矩形波の移流シミュレーションを実施する．

これまでの手法と比較して，かなり改善したように感じる．やはり完璧ではないものの，1次精度風上差分で見られた強い（強すぎる）拡散の効果や，Lax-Wendroffで現れた後方の振動も緩和している．特に，QUICKESTスキームは，元の矩形波の高さを良く保ちながら移流することができており，オーバーシュート，アンダーシュートも比較的小さい．

Kawamura-Kuwaharaの方法

Kawamura-Kuwaharaスキーム（Kawamura&Kuwahara1984⁴），あるいはKKスキームと呼ばれる本手法は，Leonardの2手法と同様に5点のステンシルを伸ばし，4次精度中心差分+4階微分の数値粘性によって移流項を近似する．

KKスキームは2次精度風上差分からスタートするので，その導出から始める．

\begin{align}
    \phi_{i-1} 
        &= \phi_{i} 
        - \Delta x \left. \frac{\partial \phi}{\partial x} \right|_{i}
        + \frac{(\Delta x)^{2}}{2} \left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i}
        + \mathcal{O} \left[ (\Delta x)^{3} \right] \\
    \phi_{i-2} 
        &= \phi_{i} 
        - 2 \Delta x \left. \frac{\partial \phi}{\partial x} \right|_{i}
        + \frac{(2 \Delta x)^{2}}{2} \left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i}
        + \mathcal{O} \left[ (\Delta x)^{3} \right] \\
\end{align}

両者の差を取り，2階微分を除去する．

\begin{align}
    4 \phi_{i-1} - \phi_{i-2}
        &= 3 \phi_{i} 
        - 2 \Delta x \left. \frac{\partial \phi}{\partial x} \right|_{i}
        + \mathcal{O} \left[ (\Delta x)^{3} \right] \\
    \therefore \left. \frac{\partial \phi}{\partial x} \right|_{i}
        &= \frac{3 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{2 \Delta x}
        + \mathcal{O} \left[ (\Delta x)^{2} \right]
\end{align}

より，次の2次精度風上差分を得る．

u \frac{\partial \phi}{\partial x} \approx 
\begin{cases}
  u \frac{3 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{2 \Delta x} & \text{if $u \ge 0$} \\
  u \frac{- \phi_{i+2} + 4 \phi_{i+1} -3 \phi_{i}}{2 \Delta x} & \text{if $u \lt 0$}
\end{cases}

これまでと同様に，$u = (u + |u|) / 2 + (u - |u|) / 2$の関係を利用して整理する．

\begin{align}
u \frac{\partial \phi}{\partial x} 
  &\approx \frac{u + |u|}{2} \frac{3 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{2 \Delta x} 
  + \frac{u - |u|}{2} \frac{- \phi_{i+2} + 4 \phi_{i+1} -3 \phi_{i}}{2 \Delta x} \\
  &= u \frac{- \phi_{i+2} + 4 \phi_{i+1} - 4 \phi_{i-1} + \phi_{i-2}}{4 \Delta x} 
  + \frac{|u| (\Delta x)^3}{4} \frac{\phi_{i+2} - 4 \phi_{i+1} + 6 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{(\Delta x)^4}
\end{align}

ここで，第1項に現れる$\phi_{i+1}$などのTaylor級数展開を考えると，以下を得る．

\begin{align}
    \phi_{i+1} - \phi_{i-1}
    &= 2 \Delta x \left. \frac{\partial \phi}{\partial x} \right|_{i}
    + 2 \frac{(\Delta x)^{3}}{6} \left. \frac{\partial^3 \phi}{\partial x^3} \right|_{i}
    + \mathcal{O} \left[ (\Delta x)^{5} \right] \\
    \phi_{i+2} - \phi_{i-2} 
    &= 4 \Delta x \left. \frac{\partial \phi}{\partial x} \right|_{i}
    + 2 \frac{(2 \Delta x)^{3}}{6} \left. \frac{\partial^3 \phi}{\partial x^3} \right|_{i}
    + \mathcal{O} \left[ (\Delta x)^{5} \right]
\end{align}

したがって，

\begin{align}
u \frac{- \phi_{i+2} + 4 \phi_{i+1} - 4 \phi_{i-1} + \phi_{i-2}}{4 \Delta x} 
&= u \frac{1}{4 \Delta x}
\left(
    4 \Delta x \left. \frac{\partial \phi}{\partial x} \right|_{i}
    + 8 \frac{(\Delta x)^{3}}{6} \left. \frac{\partial^3 \phi}{\partial x^3} \right|_{i}
    + \mathcal{O} \left[ (\Delta x)^{5} \right]
\right) \\
&= u \left. \frac{\partial \phi}{\partial x} \right|_{i}
+ \frac{u (\Delta x)^{2}}{3} \left. \frac{\partial^3 \phi}{\partial x^3} \right|_{i}
+ \mathcal{O} \left[ (\Delta x)^{4} \right]
\end{align}

これより，上式の移流項の近似では，誤差の主要項が$(\Delta x)^{2}$のオーダーにあることが分かる．

ここで，再度Taylor展開で得る関係式から，3階の微分を取り除くと，

\begin{align}
u \frac{- \phi_{i+2} + 8 \phi_{i+1} - 8 \phi_{i-1} + \phi_{i-2}}{12 \Delta x} 
&= u \left. \frac{\partial \phi}{\partial x} \right|_{i}
+ \mathcal{O} \left[ (\Delta x)^{4} \right]
\end{align}

を得る．したがって，先ほどの第1項を上式で置き換え，4次精度中心差分と4次の数値粘性で移流項を近似できると結論付けるのがKKスキームである．

u \frac{\partial \phi}{\partial x} 
\approx u \frac{- \phi_{i+2} + 8 \phi_{i+1} - 8 \phi_{i-1} + \phi_{i-2}}{12 \Delta x}
+ \frac{|u| (\Delta x)^3}{4} \frac{\phi_{i+2} - 4 \phi_{i+1} + 6 \phi_{i} - 4 \phi_{i-1} + \phi_{i-2}}{(\Delta x)^4}

さて，QUICKの紹介のときにも，4次の数値粘性が局所化された拡散効果を持つと述べた．KKにも4次の数値粘性が現れるが，これについてKawamura+1986は，

Physically, the fourth-order derivative term stands for short-range diffusion and plays an important role in stabilizing the computation.

と言葉で記すに留まっており，何故，4次の数値粘性が局所化された拡散効果に対応するのかは未だ良く分かっていない（理解を深めて補足したい）．

取り敢えず，導出はできたので，矩形波の移流シミュレーションを示す．KKに加え，QUICK, QUICKESTの結果を併記する．

KKは，QUICKよりも振動が小さく抑えられ，より良い近似となっていそうである．一方，KKとQUICKESTを比較すると，KKの方が不連続点周辺でオーバーシュート気味であり，この問題の範囲内では，QUICKESTの方が高さを維持する近似能力と振動を抑える数値粘性のバランスが良いように感じる．

最後に，三角波の移流シミュレーションを実施する．見易さのため，幾つかのスキームだけ示す．

やはり，1次精度風上差分では強い数値粘性が導入され，解のなまりが顕著である．高次スキームのQUICKESTとKKは三角波を凡そ保持している．QUICKESTはピークを減衰させているものの，裾の振動は小さい．KKはピークをより正確にキャプチャしているものの，裾付近での振動が顕著である（より長時間ではピークをオーバーシュートしそう？）．QUICKESTとKKは甲乙付け難い気がするが，実装の容易さと精度のバランスから選ぶのが良いだろう．

各手法の比較

(2023/12追記)
各手法での誤差をきちんと数値で測るようにした．初期条件に三角波を与えた問題とガウス関数を与えた問題を新たに作り，厳密解との誤差を評価した．

$t=0.2$	$t=0.6$	$t=1.0$

ここで，max errは領域内で最も大きいズレ，mean errは領域内での平均的な大きいズレを相対的に数値化したものである．三角波の問題では，max errが三角波のピーク部分で測られたり裾の振動している部分で測られたりしており，若干解釈し辛い．だが，全体の傾向から，誤差の小さいものから順に並べると，

• QUICKEST < KK < QUICK < LW < Upwind

となっていそうだ．式の複雑さからも，QUICKESTが（本稿に記した中では）最も高級な手法らしい．

終わりに

移流項の差分近似について，最も原始的な1次精度風上差分をレビューし，分子粘性とは異なる数値粘性の存在を紹介した．

また，その後提案されてきた高次スキームを幾つか紹介した．

Appendix A: 1次精度風上差分スキームの数値粘性

本文では，1次精度風上差分スキームの議論の際，移流の方向に応じて式を整理し直して数値粘性を導出したが，Taylor展開を用いても同様の導出ができる．

簡単のため，$u \ge 0$の場合を考えると，移流項は以下のように離散化できるのであった．

u \frac{\partial \phi}{\partial x} \approx u \frac{\phi_{i} - \phi_{i-1}}{\Delta x}

ここで，右辺に登場する$\phi_{i-1}$を，$\phi_{i}$まわりでTaylor展開する．

\begin{align}
u \frac{\phi_{i} - \phi_{i-1}}{\Delta x} 
&= u \frac{1}{\Delta x} \left( \phi_{i} - \phi_{i-1} \right) \\
&= u \frac{1}{\Delta x} \left( \phi_{i} - \left( \phi_{i} - \Delta x \left. \frac{\partial \phi}{\partial x} \right|_{i} + \frac{(\Delta x)^2}{2} \left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i} + \mathcal{O} \left[ (\Delta x)^3 \right] \right) \right) \\
&= u \frac{1}{\Delta x} \left( \Delta x \left. \frac{\partial \phi}{\partial x} \right|_{i} - \frac{(\Delta x)^2}{2} \left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i} + \mathcal{O} \left[ (\Delta x)^3 \right] \right) \\
&= u \left. \frac{\partial \phi}{\partial x} \right|_{i} - \frac{u \Delta x}{2} \left. \frac{\partial^2 \phi}{\partial x^2} \right|_{i} + \mathcal{O} \left[ (\Delta x)^2 \right]
\end{align}

より，やはり風上差分近似した結果に拡散の効果が現れることが分かる（それぞれの項を中心差分近似すると，本文で登場した式と同じものになる）．

Appendix B: 各スキームの収束性

各スキームはどれくらいの速度で収束するか（例えば，1次精度風上差分スキームは数値的にも本当に空間について1次精度か）を確認する．

解析解が知られている，Taylor-Green流れをベンチマークとする．まず，非圧縮性Navier-Stokes方程式は以下の通り．

\begin{alignat}{2}
    \nabla \cdot \mathbf{u}
    &= 0 
    \quad &&\text{in} \quad \Omega \times (0, T) \\
    \partial_t \mathbf{u} 
    + (\mathbf{u} \cdot \nabla) \mathbf{u}
    &= - \rho^{-1} \nabla p 
    + \nu \nabla^2 \mathbf{u}
    + \mathbf{f}
    \quad &&\text{in} \quad \Omega \times (0, T) \\
\end{alignat}

ここで，$\Omega := (0, L)^2$とする．解析解は以下の通りである（例えば，Lind&Stansby2016）．

\begin{align}
    u &=   \sin(c x) \cos(c y) \exp(- 2 c^2 \nu t) \\
    v &= - \cos(c x) \sin(c y) \exp(- 2 c^2 \nu t) \\
    p &=   \frac{\rho}{4} ( \cos(2 c x) + \cos(2 c y) ) \exp(- 4 c^2 \nu t) \\
\end{align}

$L = 1 \ [m], \ T = 0.1 \ [\mathrm{s}], \ \mathbf{f} = \mathbf{0} \ [\mathrm{kg \cdot m / s^2}], \ \rho = 1 \ [\mathrm{kg / m^3}], \ \nu = 0.01 \ [\mathrm{m^2 / s}], \ c = 2 \pi$としてシミュレーションを行い，$t = T \ [\mathrm{s}]$での数値解と解析解との間の相対$L^2$誤差を測る．特に，水平速度$u$に関する誤差を尺度とする．

$\Delta t = 1 \times 10^{-5} \ [\mathrm{s}]$に固定してChorinの射影法（Chorin1968）を用いて時間積分を実行し，拡散項と圧力勾配項は2次精度中心差分を用いる．移流項の近似スキームを1次精度風上差分と3次精度風上差分（Kawamura&Kuwahara1984）に切り替えて比較する．結果は以下の通りである．

解の例	収束性

左図は，$\Delta x = 1 \times 10^{-2} \ [\mathrm{m}], \Delta t = 1 \times 10^{-3} \ [\mathrm{s}]$でのKKスキームによる数値解であり，コードのチェックのためである．

右図は，$\Delta t = 1 \times 10^{-5} \ [\mathrm{s}]$に固定して空間解像度を変化させたものである．両者には若干だけ違いがあるが，ともに凡そ1のスロープに乗っており，違いが殆ど無い．恐らく時間進行スキームが1次精度なので，これに空間スキームの精度が埋もれてしまったのだろうと考える（cf. 非圧縮性粘性流れの数値計算プログラムにおける時間積分精度を調査する，非圧縮性 Navier-Stokes 方程式の数値解法４：ソルバー精度検証）．Crank-Nikolson法，Adams-Bashforth/Moulton法，4段4次Runge-Kutta法など，時間スキームを高精度化して再度検討する必要があるだろう．また各手法を実装できたときに帰ってこようと思う．

(2024/06追記)
過去の自分は，何故わざわざTaylor-Green渦という難しい問題を選んだのだろう．．．本文で議論していたシンプルな1次元の移流を対象に，収束性を見る．

$\Delta t = 1 \times 10^{-5}$に固定して陽的Euler法で時間積分を実行した（Courant数は常に$5 \times 10^{-3}$を下回った）．$T = 1$での，厳密解との誤差の平均値と最大値をプロットする．

	矩形波	Gauss波
$L^2$
$L^{\infty}$

切り立った矩形波では，どうしても，不連続面付近で大きい誤差が発生してしまう様子である．解像度を上昇させても，誤差の最大値は低減しない．平均値を見ても，期待ほどの収束の速度は出ないようだ．いずれの方法も1次の速度すら出ていない．

滑らかなGauss波では，誤差の平均値，最大値のいずれを見ても，ある程度直感に合う結果となった．最も原始的な風上差分が~1次精度ほど，Lax-Wendroffが凡そ2次精度を達成している．QUICKは殆どLax-Wendroffと同じくらいに見える．QUICKESTとKKは2次精度より早く，凡そ3次の速度が出ていると言って良いだろうか．

1. 私は厳密に数学を勉強してきた身ではないので，この日本語は誤解を招く，ないし，間違いを含む可能性があります．私の解釈を日本語に起こしたに過ぎません． ↩

2. 移流・非線形性に関する時間増分の制約はCourant+1967が，拡散に関する時間増分の制約はNeumann&Richtmyer1950が与えています．例えば，Navier-Stokes方程式には移流・拡散両者の効果が入っていますから，そのシミュレーションでは，両者を満たす数値の中から，なるべく大きい数値を時間増分として採用することになります． ↩

3. ここは本当は良く理解できていません．講義を受けた先生の言葉をそのまま記しています．背景の数学が理解できたら補足したいと思っています． ↩

4. 「Direct and Large Eddy Simulation of Turbulence」のpp 227–244:「Direct Simulation of High-Reynolds-Number Flows by Finite-Difference Methods」のチャプターも参考に． ↩