レッスン 特2

分離解法

1次元移流方程式の右辺の値をとりあえずひとまとめに外力項（source項）として G と置きます。

・・・（Ex2.1）

これが新たな基礎式となります。このままの形ではCIP法を適用出来ません。右辺はゼロじゃなきゃ駄目なんです。そこで、(Ex2.1)式を便宜的に以下の様に分離出来ると仮定します。

・・・（Ex2.2 a, b）

そして、これをどうするかと言うと、f ( t ), f ( t + Δt )をそれぞれ、fⁿ, fⁿ⁺¹、と置くと、(Ex2.1)式でfⁿ→fⁿ⁺¹を求めるのを、いったん(Ex2.2a)式でfⁿ→と中間量を求め、(extra.2-2b)式で→fⁿ⁺¹を求める。このように2段階で計算します。これを分離解法と言い、移流を計算する段階をAdvection Phase、それ以外をNon-Advection Phaseと呼びます。イメージとしては、後で図にも示しますが、仮に G が拡散項だったとしたら、まず拡散する部分(Non-Advection)だけ計算して、その後、拡散したものを移流(Advection)させるといった感じです。整理すると次のようになります。


Advection PhaseにはCIP法がそのまま使えます。で、CIP法を使うとなると、 f の空間微分量 fx が必要なので、同様に更新します。


こんな文章と式だけ見てもよく分からないと思うので、図で見てみましょう。イメージが湧きやすいと思います。


[ 図-Ex2.1 ]
ちなみに、この分離解法という手法は計算上のテクニックです。理論的に分離してよいという根拠はありません。

Non-Advection Phase（非移流項）の計算

では、実際に解いてみましょう。(Extra.2-2a)式より

・・・（Ex2.3）

よって、

・・・（Ex2.4）

となります。計算点を i とすれば、

・・・（Ex2.5）

です。この式でNon-Advection Phaseの f が更新されます。
続いて微分量の更新です。 （Ex2.4）式の両辺を x で微分すると、

・・・（Ex2.6）

G の空間微分を中央差分で表すと、

・・・（Ex2.7）

です。一方、G は(Ex2.3)式で表されるように、

・・・（Ex2.8）

なので、これを(Ex2.7)式に代入して、

・・・（Ex2.9）

となり、この式でNon-Advection Phaseの微分量が更新出来ます。

Advection Phase（移流項）の計算

Advection Phaseの(Ex2.2b)にはこれまでにやってきたCIP法が使えるので、使います。

・・・（Ex2.10）

係数 a, b は

・・・（Ex2.11）

続いて微分量です。(Ex2.2b)式の両辺を x で微分すると、

・・・（Ex2.12）

これは、どっかで見た形。そう、(Ex2.1)式と同じ形です。と言うわけで、分離します。

・・・（Ex2.13a, b）

じゃあ、またCIP法？それはかったるい。CIP法では F の空間微分量を必ず出すので、それを使います。すなわち、(Ex2.13a)式は、

・・・（Ex2.14）

で求め、(Ex2.13b)式は、右辺に中央差分を適用して差分化すると、

・・・（Ex2.15）

これをまとめて、

・・・（Ex2.16）

これで、Advection Phaseの微分量も更新出来ました。

まとめ

以上をまとめると、


[ 図-Ex2.2 ]

-- f の更新 --
	・・・（Ex2.5）
-- f の微分量の更新 --
	・・・（Ex2.9）

[ 図-Ex2.3 ]

-- ｆ の更新 --
	・・・（Ex2.10）
ただし、
	・・・（Ex2.11）
-- f の微分量の更新 --
	・・・（Ex2.14）
	・・・（Ex2.16）

問題Ex2-1.

それでは、計算行ってみよう。 source項はとりあえず中央差分のところでやった拡散項とします。これだとただの移流拡散方程式で面白くないので、流れが時間的、空間的に変化するように設定してみましょう。

解説

[ DL ソースコード ]

[ イメージ（クリックすると大きくなります。） ]


[ 動-Ex2.1 ]

解説

いまいち空間的な流れの変化が分かりにくいですかね．．．でも、こんなのどうやって計算結果を検証したら良いんでしょうね？だれかこれを理論的に解いて。