数学とか語学とか楽しいよね

ドイツ語、ロシア語、オランダ語、英語、フランス語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。

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

風上差分の中心差分+数値拡散表示の導出をやります。ポイントは風向きの正負によらない表示です。まず、 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}

となります。この式の一項目は、一階微分に対する中心差分に、二項目は二階微分に対する中心差分となっているので、風上差分の中心差分+数値拡散表示を得ることができました。これは風向きの正負によらない表示でもあります。