平林孝太＠東京大学研究のはなし、個人的なメモ書きなどのページになる予定です。。。 |
5th-order WENO interpolation
概要
まず３次精度 ENO (Essentially Non-Oscillatory) スキームとは、 ５つのステンシル \(x_{i-2} \sim x_{i+2}\) を用いて３つの多項式 \(p_1(x_{i-2},x_{i-1},x_i), \ \ p_2(x_{i-1}, x_i, x_{i+1}), \ \ p_3(x_i, x_{i+1}, x_{i+2})\) を構成し、 一番滑らかなものを補間関数（あるいは再構築）として採用するもの。 衝撃波のような不連続面があるときには、不連続を含まないステンシルから作られた多項式を自動的に選べば、 衝撃波近傍でも滑らかな補間ができる。 ただし厳密にTVD条件を満たしているわけではないので、Essentially Non-Oscillatory である。
この ENO スキームには、ステンシルを５つ使っているにもかかわらず３次精度しか出ないという欠点がある。 WENO (Weighted ENO) スキームは、滑らかな領域で \(p_1 \sim p_3\) の重み付き平均で補間関数を構成することにより ５次精度を達成する。つまり、５つのステンシルを使った補間関数に一致するように、３次の打ち切り誤差を相殺する（ような重みで足し合わせる）。 不連続近傍では ENO と同じように、もっとも滑らかな補間に漸近するように重みを調節し数値振動を回避する。
Godunov の定理によれば、２次精度以上のすべての線型スキームでは数値振動が発生するが、 WENO は非線型の重み付けを用いることでこれを克服したひとつの例と言える。
実装
変数 \(u\) を \(x=x_{i+1/2}\)に補間することを考えよう。 ここで扱うのは単なる interpolation であって、reconstruction とは係数が異なるので注意。
まず３つの補間関数 \(p_1 \sim p_3\) から \(x_{i+1/2}\) での値を求める。
\begin{align*} u_1 = p_1(x_{i+1/2}) &= \frac{3}{8}u_{i-2} - \frac{5}{4}u_{i-1} + \frac{15}{8}u_{i } \\ u_2 = p_2(x_{i+1/2}) &= - \frac{1}{8}u_{i-1} + \frac{3}{4}u_{i } + \frac{ 3}{8}u_{i+1} \\ u_3 = p_3(x_{i+1/2}) &= \frac{3}{8}u_{i } + \frac{3}{4}u_{i+1} - \frac{ 1}{8}u_{i+2} \end{align*}
これらは次の重みで足し合わせるとちょうど５次精度の補間に一致する。
\[ u_{i+1/2} = \frac{1}{16} u_1 + \frac{5}{8} u_2 + \frac{5}{16} u_3 \]
WENO の真髄は、補間関数の滑らかさに応じてこの重み付けを非線型に調整することであった。 すべてのステンシルで滑らかなら上の式に一致し、 滑らかでないステンシルを検出したら他の補間関数を使うように重みを決定する。 滑らかさの指標は Smoothness Indicator と呼ばれる次の式で評価される。
\begin{align*} \beta_1 &= \frac{1}{3} \left( 4 u_{i-2} u_{i-2} - 19 u_{i-2} u_{i-1} + 25 u_{i-1} u_{i-1} + 11 u_{i-2} u_{i } - 31 u_{i-1} u_{i } + 10 u_{i } u_{i } \right) \\ \beta_2 &= \frac{1}{3} \left( 4 u_{i-1} u_{i-1} - 13 u_{i-1} u_{i } + 13 u_{i } u_{i } + 5 u_{i-1} u_{i+1} - 13 u_{i } u_{i+1} + 4 u_{i+1} u_{i+1} \right) \\ \beta_3 &= \frac{1}{3} \left( 10 u_{i} u_{i} - 31 u_{i+1} u_{i+1} + 25 u_{i+1} u_{i+1} + 11 u_{i } u_{i+2} - 19 u_{i+1} u_{i+2} + 4 u_{i+2} u_{i+2} \right) \end{align*}
\(\beta\) は（パッと見た感じよくわからないけど）２階と４階の差分から構成されていて、 滑らかな領域でゼロに近づく関数になっているので、 逆数をとって適当に規格化したものを重みとすればよさそうだと予想できる。 実際に用いられるのは、
\begin{align*} u_{i+1/2} &= \tilde{w}_1 u_1 + \tilde{w}_2 u_2 + \tilde{w}_3 u_3 \\ \tilde{w}_k &= \frac{w_k}{\sum_k w_k} \\ w_k &= \frac{d_k}{(\varepsilon + \beta)^2} \end{align*}
のような表式。 \(d_k\)は上で述べた線型補間の重みで、\(\varepsilon\)はゼロ割を防ぐための微小定数。 標準的には \(\varepsilon=10^{-6}\) あたりが使われる。
以上から必要な手順は、
- 線型補間値 \(u_1, u_2, u_3\) を求める
- Smoothness Indicator \(\beta_1, \beta_2, \beta_3\) を評価
- 非線型重み \(\tilde{w}_1, \tilde{w}_2, \tilde{w}_2\) を求める
- 重み付き平均を最終的な補間値 \(u_{i+1/2}\) として採用
任意座標への補間
何かの都合で、半整数座標 \(x_{i+1/2}\) のみでなく任意の座標
\(x_{\alpha}=x_{i}+\alpha \Delta x\) での補間値が必要になることもある。
この時にも WENO スキームをそのまま応用することができる。
補間関数 \(p_1 \sim p_3\) を用いて \(x_{\alpha}\) での値を求めると、
\begin{align*}
u_1 = p_1(x_{\alpha}) &= \frac{u_{i } - 2u_{i-1} + u_{i-2}}{2} \alpha^2
+ \frac{3u_{i } - 4u_{i-1} + u_{i-2}}{2} \alpha + u_{i} \\
u_2 = p_2(x_{\alpha}) &= \frac{u_{i+1} - 2u_{i } + u_{i-1}}{2} \alpha^2
+ \left( u_{i+1} - u_{i-1} \right) \alpha + u_{i} \\
u_3 = p_3(x_{\alpha}) &= \frac{u_{i+2} - 2u_{i+1} + u_{i }}{2} \alpha^2
+ \frac{3u_{i+2} - 4u_{i+1} + u_{i }}{2} \alpha + u_{i}
\end{align*}
これらの線型結合で５次精度を達成するには、
上と同様の手続きで \(d_k\) を修正するだけでよくて
\begin{align*}
d_1 &= \frac{1}{12} (1-\alpha)(2-\alpha) \\
d_2 &= \frac{1}{6} (2-\alpha)(2+\alpha) \\
d_3 &= \frac{1}{12} (1+\alpha)(2+\alpha)
\end{align*}
とすると \(x_{\alpha}\)
での Non-oscillatory な補間値が求まる。
システム方程式への拡張
MHDのようなシステム方程式に拡張したとき、
各物理量に独立に WENO 補間を施すと、
バラバラなステンシルを用いることによる不整合が原因で数値振動が発生する。
これを回避するためには、Global Smoothness Indicator を用いてすべての物理量で共通の重み関数を採用するとよい。
Global Smoothness Indicator は、各物理量 \(u\) に対して評価した Smoothness Indicator \(\beta\)
に物理量 \(u\) の空間平均を重み付けして足し合わせたもの。
別の流派としては、密度 \(\rho\) の Smoothness Indicator をすべての物理量で共通して用いることもある。
いずれにしても、各物理量で異なる重みを用いることは数値振動の温床となるので避けるべきであろう。