1次元の移流方程式をMUSCL使って解くことについて

それなりに前から仕事の合間に圧縮性のコードを解こうとしてるんだけど、まあMUSCLが多いですから、それ使おうとしてます。で、1次元の移流方程式を解こうとしてるんですが、そもそも衝撃波の問題を解くには質量保存と、運動量保存と、エネルギー保存の3つの保存則を解くんだけど、それに至るのに1変数の問題を解くのになんか意義があるんだろうか。
いや、あるんだけど、なんかもう俺みたいな休日だけコードを組むようなのがこんな悠長なやり方でいけるのか謎。どっかにコードが転がってればそれ使いたいんだけど、やっぱ自分で書いた方がすべてを分って使える訳で悩みはつきません。
結局は、

$\frac{\partial Q}{\partial t} + \frac{\partial F}{\partial x} = 0$ ,
を解くのに、流束FをF=UQにして、流束と差分を分離して、格子の間の不連続を2次関数で補完しますみたいなことなんだけど、

$Q_{n+1} = Q_n - \frac{\Delta t}{\Delta x}\left(F_{i+1/2}-F_{i-1/2}\right)$ ,
という差分の流束を、離散化した格子の境界上で補完する。

$F_{i+1/2}=\frac{1}{2}\left\{F(Q)^R_{i+1/2}+F(Q)^L_{i-1/2}-|U_{i+1/2}\left(Q^R_{i+1/2}-Q^L_{i-1/2}\right)\right\}$ ,
φを任意の変数とすると、色々とあって空間的な補完は、

$\emptyset^R_{i+1/2}=\emptyset_{i+1}-\frac{1}{4}\left\{(1-\kappa)\bar{\Delta}_+|_{i+1}+(1+\kappa)\bar{\Delta}_-|_{i+1}\right\}$ ,
$\emptyset^L_{i+1/2}=\emptyset_{i}+\frac{1}{4}\left\{(1-\kappa)\bar{\Delta}_-|_{i}+(1+\kappa)\bar{\Delta}_+|_{i}\right\}$ ,
$\bar{\Delta}_+|_i=\text{minmod}(\Delta_+, b\Delta_-)$ ,
$\bar{\Delta}_-|_i=\text{minmod}(\Delta_-, b\Delta_+)$ ,
$\Delta_+=\emptyset_{i+1}-\emptyset_i$
$\Delta_-=\emptyset_i-\emptyset_{i-1}$

なんてかける*1。

これが1変数のときで、3変数になると、そもそもが、支配方程式が、

$\frac{\partial Q^j}{\partial t}+\frac{\partial F^j}{\partial x}=0$ ,
になって、これを解こうとすると、

$Q^j_{n+1}=Q^j_{n}-\frac{\Delta t}{\Delta x}\left\{F^j_{i+1/2}-F^j_{i-1/2}\right\}$ ,
になって、まあ色々とあって、流束差分離なんてすると、

$F^j_{i+1/2}=\frac{1}{2}\left\{F^j(Q)^R_{i+1/2}+F^j(Q)^L_{i+1/2}-|U^{jk}_{i+1/2}|\left(\left.Q^k\right.^R_{i+1/2}-\left.Q^k\right.^L_{i+1/2}\right)\right\}$ ,
になるわけだけど、やっぱここで2次関数で補完するのとか、制限関数のサブルーチンとか、その辺のものを作らないといけない訳だから、これはこれでやっぱり意味のあることなのだとやっと気がついた。因に上の式では添字kについてEinsteinの縮訳記号を使ってます。
それにしてもまだまだ先が長い。

*1:minmod関数を制限関数に使ったとき。slopeは核の面倒なのでスルー。東大の藤井先生の本なんかに詳しく載ってます。まとまりが悪いし、誤記があったりするけど、日本語で書いてある数少ない本のうちの一つではある。