【数値計算】風上差分の中心差分+数値拡散表示

風上差分の中心差分+数値拡散表示の導出をやります。ポイントは風向きの正負によらない表示です。まず、 $u_i$ を節点 $i$ における濃度ないしは熱、 $V_i$ を節点 $i$ における流速、 $x$ を空間座標としたとき、風上差分は

$\displaystyle \left.V\frac{\partial u}{\partial x}\right|_i = \begin{cases} \displaystyle V_i\frac{u_i - u_{i-1}}{\Delta x} & (V_i \geq 0)\\ \displaystyle V_i\frac{u_{i+1} - u_i}{\Delta x} & (V_i < 0) \end{cases}$
のように書けます。風上側（情報が伝わってくる側）の値を使うことで安定的に計算を行うことができます。ここで

$\displaystyle \frac{V_i+|V_i|}{2}= \begin{cases} \displaystyle V_i & (V_i \geq 0)\\ \displaystyle 0 & (V_i < 0) \end{cases}$ 　 $\displaystyle \frac{V_i-|V_i|}{2}= \begin{cases} \displaystyle 0 & (V_i \geq 0)\\ \displaystyle V_i & (V_i < 0) \end{cases}$

というちょっとテクニカルな関係を使います。落ち着いて絶対値をはずしてみればたいしたことはありません。上の関係をよく見ると、 $V_i \geq 0$ のとき $V_i=\frac{V_i+|V_i|}{2}$ に、 $V_i < 0$ のとき $V_i=\frac{V_i-|V_i|}{2}$ になっていることに気づきます（よく見て下さい）。なので、これらを風上差分の式の $V_i$ に代入してみましょう。すると

$\displaystyle \left.V\frac{\partial u}{\partial x}\right|_i = \begin{cases} \displaystyle \frac{V_i+|V_i|}{2} \frac{u_i - u_{i-1}}{\Delta x} & (V_i \geq 0)\\ \displaystyle \frac{V_i-|V_i|}{2} \frac{u_{i+1} - u_i}{\Delta x} & (V_i < 0) \end{cases}$
となります。さらに、 $\boldsymbol{\frac{V_i+|V_i|}{2} \frac{u_i - u_{i-1}}{\Delta x}}$ と $\boldsymbol{\frac{V_i-|V_i|}{2} \frac{u_{i+1} - u_i}{\Delta x}}$ を足してしまいます！するとそれは風上差分と等しくなっているのです！すなわち

$\displaystyle \left.V\frac{\partial u}{\partial x}\right|_i = \frac{V_i+|V_i|}{2} \frac{u_i - u_{i-1}}{\Delta x} + \frac{V_i-|V_i|}{2} \frac{u_{i+1} - u_i}{\Delta x}$
です。何故そんなことができるか？それは場合分けして考えるとわかります。まず $V_i \geq 0$ のときは

$\displaystyle \frac{V_i+|V_i|}{2}=V_i$ 　かつ　 $\displaystyle \frac{V_i-|V_i|}{2}=0$

なので和の二項目が消えて風上差分と等しくなります。また $V_i < 0$ のときは

$\displaystyle \frac{V_i+|V_i|}{2}=0$ 　かつ　 $\displaystyle \frac{V_i-|V_i|}{2}=V_i$

なので和の一項目が消えて風上差分と等しくなります。

さらに式を $V_i$ と $|V_i|$ で整理すると結局

$\displaystyle \left.V\frac{\partial u}{\partial x}\right|_i = V_i \frac{u_{i+1} - u_{i-1}}{2 \Delta x} - \frac{|u_i| \Delta x}{2} \frac{u_{i+1} - 2 u_i + u_{i-1}}{(\Delta x)^2}$
となります。この式の一項目は、一階微分に対する中心差分に、二項目は二階微分に対する中心差分となっているので、風上差分の中心差分+数値拡散表示を得ることができました。これは風向きの正負によらない表示でもあります。