差分法の基礎¶

著者: 松元亮治（千葉大学）

流体・磁気流体方程式を差分法を用いて数値的に解く際に必要になる基礎的事項について解説する。波の伝播をあらわす線形移流方程式や非線形のBurgers方程式をとりあげ、差分解法の数値的な安定性や数値振動について論じる。特に、数値的な安定性に優れ、非物理的な数値振動を起こさない差分法として風上差分法を紹介する。

差分近似¶

変数 $u$ が空間座標 $x,y$ に依存するという2次元問題を考える。

2次元メッシュの図¶

2次元空間を図のような格子に区切り、各格子点の座標を $(x_i, y_j)$ とする。格子間隔は $x$ 方向が $\Delta x$ 、 $y$ 方向が $\Delta y$ とする。 $x_{i \pm 1} = x_i \pm \Delta x$ 、 $y_{j \pm 1} = y_j \pm \Delta y$ である。以下、格子点番号 $(i,j)$ を用いて $u_{i,j} = u(x_i, y_j)$ のように略記する。

着目している点 $(x_i, y_j)$ のまわりでテイラー展開すると、

(1)¶ $u_{i+1,j} = u(x_i+\Delta x, y_j) = u_{i,j} +\Delta x \left (\displaystyle{{\partial u} \over {\partial x}} \right)_i +\displaystyle{{{\Delta x}^2} \over {2!}} \left (\displaystyle{{\partial^2 u} \over {\partial x^2}} \right)_i +\displaystyle{{{\Delta x}^3} \over {3!}} \left (\displaystyle{{\partial^3 u} \over {\partial x^3}} \right)_i + ...$

(2)¶ $u_{i-1,j} = u(x_i-\Delta x, y_j) = u_{i,j} -\Delta x \left (\displaystyle{{\partial u} \over {\partial x}} \right)_i +\displaystyle{{{\Delta x}^2} \over {2!}} \left (\displaystyle{{\partial^2 u} \over {\partial x^2}} \right)_i -\displaystyle{{{\Delta x}^3} \over {3!}} \left (\displaystyle{{\partial^3 u} \over {\partial x^3}} \right)_i + ...$

式 (1) から式 (2) を引くと

(3)¶ $u_{i+1,j} - u_{i-1,j} = 2\Delta x \left (\displaystyle{{\partial u} \over {\partial x}}\right )_i + O(\Delta x^3)$

したがって、

(4)¶ $\left (\displaystyle{{\partial u} \over {\partial x}} \right )_i = \displaystyle{{u_{i+1,j} - u_{i-1,j}} \over {2\Delta x}} + O(\Delta x^2)$

すなわち、 $(i,j)$ 点における $u$ の $x$ 方向の微分係数 $(\partial u/\partial x)_i$ が $\Delta x^2$ の誤差を含む近似のもとで( $\Delta x$ について2次の精度で)以下のように求まる。

(5)¶ $\left (\displaystyle{{\partial u} \over {\partial x}} \right)_i = \displaystyle{{u_{i+1,j} - u_{i-1,j}} \over {2\Delta x}}$

これを中心差分の式と言う。

同様にして、 $\Delta x$ について1次の精度で以下の差分近似式が得られる。

前進差分：

(6)¶ $\left (\displaystyle{{\partial u} \over {\partial x}}\right)_i = \displaystyle{{u_{i+1,j} - u_{i,j}} \over {\Delta x}}$

後退差分：

(7)¶ $\left(\displaystyle{{\partial u} \over {\partial x}}\right)_i = \displaystyle{{u_{i,j} - u_{i-1,j}} \over {\Delta x}}$

式 (1) と式 (2) を加えると

(8)¶ $u_{i+1,j} + u_{i-1,j} = 2u_{i,j} + \Delta x^2 \left(\displaystyle{{\partial^2 u} \over {\partial x^2}}\right)_i + O(\Delta x^4)$

したがって、

(9)¶ $\left(\displaystyle{{\partial^2 u} \over {\partial x^2}}\right)_i = \displaystyle{{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}} \over {\Delta x^2}} + O(\Delta x^2)$

これより、 $u$ の $x$ に関する2階微分の係数 $(\partial^2 u/\partial x^2)_i$ を $\Delta x$ について2次の精度で以下のように近似することができる。

(10)¶ $\left(\displaystyle{{\partial^2 u} \over {\partial x^2}}\right)_i = \displaystyle{{u_{i+1,j} - 2u_{i,j} + u_{i-1,j}} \over {\Delta x^2}}$

同様に、

(11)¶ $\left(\displaystyle{{\partial^2 u} \over {\partial y^2}}\right)_j = \displaystyle{{u_{i,j+1} - 2u_{i,j} + u_{i,j-1}} \over {\Delta y^2}}.$

線形スカラー移流方程式の差分解法¶

１次元線形スカラー移流方程式¶

流体・磁気流体方程式の本質は波の伝播にある。この部分だけを取り出して次のような方程式を考える。

(12)¶ $\displaystyle{{\partial u} \over {\partial t}} + c \displaystyle{{\partial u} \over {\partial x}} = 0$

ただし、 $c$ は定数で $c > 0$ とする。この方程式は、スカラー量 $u$ の空間分布が、一定の速度 $c$ で伝播することをあらわす波動方程式である。

方程式 (12) の厳密解は

(13)¶ $u(x,t) = u(x-ct,0)$

である。これは、時刻 $t > 0$ におけるスカラー量 $u$ のプロフィールは $t=0$ のスカラー量 $u$ のプロフィールが形を保って $ct$ だけ平行移動した形になることをあらわす。

1次元スカラー移流問題の初期条件と時間発展¶

いま、図 1次元スカラー移流問題の初期条件と時間発展のように初期に $x \ge 0$ で $u=u_1$ 、 $x < 0$ で $u=u_2$ のように $x=0$ で不連続な分布を考えてみると $t > 0$ での厳密解は右図のような形になる。

FTCSスキーム¶

1次元線形スカラー移流方程式 (12) を時間について現在の時刻 $t_n$ と $\Delta t$ 後の時刻 $t_{n+1}=t_n+\Delta t$ の間で前進差分、空間については中心差分をとって差分化すると次式を得る。ここで、 $u_j^n = u(x_j,t_n)$ である。

(14)¶ $\displaystyle{{u_{j}^{n+1}-u_{j}^{n}} \over {\Delta t}} + c \displaystyle{{u_{j+1}^{n}-u_{j-1}^n} \over {2\Delta x}} = 0$

このような差分のとり方をFTCSスキーム(Forward in Time and Centered Difference in Space)と言う。これを整理すると、

(15)¶ $u_j^{n+1} = u_j^n - \displaystyle{1 \over 2} \nu (u_{j+1}^n - u_{j-1}^n)$

ここで、 $\nu$ は次式で定義される数であり、 クーラン数 と呼ばれる。

(16)¶ $\nu \equiv c \displaystyle{{\Delta t} \over {\Delta x}}$

式 (15) の右辺は時刻 $t_n$ での値、左辺は時刻 $t_{n+1}=t_n+\Delta t$ での値だけで書けている。したがって、時刻 $t_n$ での各格子点での値がわかっていれば直ちに1タイムステップ後( $t_{n+1}$ )の各格子点での値を計算することができる。このような解法のことを 陽解法 と言う。

FTCSスキームにおける変数の依存関係¶

FTCSスキームにおける変数の依存関係を図示すると図 FTCSスキームにおける変数の依存関係のようになる。矢印は時刻 $t_{n+1}$ の白丸の点の値を計算するのに時刻 $t_n$ の黒丸の格子点の値を使うことを示す。

1次元波動伝播のシミュレーションを行うアルゴリズムは一般に次のようになる。

各メッシュ点の座標値 $x_j$ をセットする(メッシュ生成)
各メッシュ点の初期値 $u_j(t=0)$ をセットする(初期条件)
時刻 $t$ が、あらかじめ決められた終了時刻 $t_{\rm end}$ に達するまで、あるいは決められた回数だけ、以下を繰り返す

左右の境界を除く各格子点について $\Delta t$ 後の値を差分式にもとづいて計算する(時間積分)。たとえばFTCSスキームの場合には計算式 (15) を用いる。

左右の境界の値を境界条件から決める。たとえば隣接点と同じ値を入れる (境界条件の適用)

時刻を $\Delta t$ だけ進める

左図：FTCSスキームで、初期値として、 $j=1,...,50$ に対して $u=1$ 、 $j=51,...,100$ に対して $u=0$ とし、クーラン数 $\nu=c\Delta t/\Delta x=0.25$ で1ステップ、2ステップ、3ステップ計算したときの $u$ をプロットした図。右図：50ステップ、100ステップ計算したときの $u$ をプロットした図。¶

FTCSスキームを用いて１次元線形スカラー移流方程式を解いた結果を上図に示す。波は形を保って伝わらずに振動が発生してしまっている。この振動は物理的な理由で発生しているのではなく、数値的不安定性によるものである。なぜこのような数値振動が発生してしまうのか、次節で説明する。

FTCSスキームの数値的安定性¶

Von Neumannの安定性解析¶

前節のFTCSスキームによって1次元波動伝播のシミュレーションを行ってみると解が激しく振動して数値的に不安定になってしまうことがわかった（ FTCSの結果の図）。この不安定性の原因を調べるために

(17)¶ $u_j^n={\rm{cos}}(j\theta)$

を差分式 (15) に代入してみる。ここで $\theta$ は、波の波数を $k$ として $\theta=k\Delta x$ であらわされる量である。たとえば $\theta=\pi$ のとき $u_j^n$ は安定性解析結果左図のように2メッシュで1波長の波、 $\theta=\pi/3$ のときは右図のように6メッシュで1波長の波をあらわす。

安定性解析結果。メッシュ番号を $j$ としたときの $u_j^n={\rm cos}(j\theta)$ のプロフィール。左図: $\theta=\pi$ の場合。右図： $\theta=\pi/3$ の場合。¶

その結果は

(18)¶ $u_j^{n+1}={\rm cos}(j\theta)+\nu {\rm sin}\theta {\rm sin}(j\theta) =Re \left[(1-i\nu{\rm sin}\theta) e^{ij\theta}\right]$

これをもう一度差分式に代入すると

(19)¶ $u_j^{n+2}=(1-\nu^2{\rm sin}^2\theta){\rm cos}(j\theta)+2\nu {\rm sin}\theta {\rm sin}(j\theta)\\ = Re \left[(1-i\nu{\rm sin}\theta)^2 e^{ij\theta}\right]$

である。ここで、 $i$ は虚数単位、 $Re$ は実部をとることをあらわす。以上からわかるように、

(20)¶ $u_j^{n+k}=Re \left[(1-i\nu{\rm sin}\theta)^k e^{ij\theta}\right]$

が成り立つ。

以上と同様に、式 (17) を複素数に拡張した $u_j^{n}=e^{ij\theta}$ を差分式 (15) に代入すると

(21)¶ $u_j^{n+k}=(1-i\nu{\rm sin}\theta)^k u_j^{n}$

が成り立つ。

差分法(差分スキーム)の数値的安定性を導くひとつの方法として、

(22)¶ $u_j^{n} = g^n e^{ij\theta}$

を差分式に代入して複素増幅率 $g$ を求め、1タイムステップ間の振幅の増幅率 $|g| \le 1$ となる条件を求める方法がある。これを Von Neumann の安定性解析と言う。

$u_j^n = g^n {\rm exp} (i j \theta)$ をFTCS差分式に代入すると

(23)¶ $g = 1 - \displaystyle{1 \over 2} \nu (e^{i \theta} - e^{-i \theta})\\ = 1 - i \nu {\rm sin} \theta$

したがって

(24)¶ $|g|^2 = 1 + \nu^2 {\rm sin}^2 \theta \ge 1$

以上の結果より、 $\theta=0$ の場合を除いてFTCSスキームは常に不安定になる。

テイラー展開による方法¶

差分化した式にテイラー展開を適用して差分式が満たす偏微分方程式を導くことによってもFTCSスキームが数値的に不安定であることを示すことができる。 $t_{n+1} = t_n + \Delta t$ 、 $x_{j \pm 1} = x_j \pm \Delta x$ を用いると、

(25)¶ $u_j^{n+1} = u_j^n + \displaystyle{{\partial u} \over {\partial t}} \Delta t + \displaystyle{1 \over 2} \displaystyle{{\partial^2 u} \over {\partial t^2}} \Delta t^2 + ...$

(26)¶ $u_{j+1}^n = u_j^n + \displaystyle{{\partial u} \over {\partial x}} \Delta x + \displaystyle{1 \over 2} \displaystyle{{\partial^2 u} \over {\partial x^2}} \Delta x^2 + \displaystyle{1 \over 6} \displaystyle{{\partial^3 u} \over {\partial x^3}}\Delta x^3 + ...$

(27)¶ $u_{j-1}^n = u_j^n - \displaystyle{{\partial u} \over {\partial x}} \Delta x + \displaystyle{1 \over 2} \displaystyle{{\partial^2 u} \over {\partial x^2}} \Delta x^2 - \displaystyle{1 \over 6} \displaystyle{{\partial^3 u} \over {\partial x^3}}\Delta x^3 + ...$

FTCSスキームの差分式

(28)¶ $u_j^{n+1} - u_j^n = -\displaystyle{1 \over 2} \nu (u_{j+1}^n - u_{j-1}^n)$

の左辺に (25) 、右辺に (26) 、 (27) を代入すると、

(29)¶ $\displaystyle{{\partial u} \over {\partial t}}\Delta t + \displaystyle{1 \over 2} {{\partial ^2 u} \over {\partial t}^2} \Delta t^2 = -\nu \left (\displaystyle{{\partial u} \over {\partial x}}\Delta x + \displaystyle{1 \over 6} \displaystyle{{\partial^3 u} \over {\partial x^3}}\Delta x^3 + ... \right)$

これを整理すると

(30)¶ $\displaystyle{{\partial u} \over {\partial t}} + c \displaystyle{{\partial u} \over {\partial x}} = -\displaystyle{1 \over 2} \displaystyle{ {\partial^2 u} \over {\partial t^2}}\Delta t - \displaystyle{c \over 6} {{\partial^3 u} \over {\partial x^3}} \Delta x^2 + ...$

ここで、解くべき偏微分方程式

(31)¶ $\displaystyle{{\partial u} \over {\partial t}} = -c \displaystyle{{\partial u} \over {\partial x}}$

より

(32)¶ $\displaystyle{{\partial^2 u} \over {\partial t^2}} = c^2 \displaystyle{{\partial^2 u} \over {\partial x^2}}$

であることを用いると

(33)¶ $\displaystyle{{\partial u} \over {\partial t}} + c \displaystyle{{\partial u} \over {\partial x}} = -\displaystyle{{c^2} \over 2} \displaystyle{ {\partial^2 u} \over {\partial x^2}}\Delta t - \displaystyle{c \over 6} {{\partial^3 u} \over {\partial x^3}} \Delta x^2 + ...$

右辺が差分化によって新たに加わった項である。右辺第1項は負の拡散係数を持つ拡散項になっている。「正の拡散」は物理量の値のピークをなまらせる働きがあるが、「負の拡散」では物理量が周囲よりもわずかに高い値を持つ部分があるとこのピークがどんどん大きくなるという不安定性を生ずる。

よって、テイラー展開法からもスカラー移流方程式のFTCSスキームは数値的に不安定であることがわかる。

Lax-Friedrich のスキーム¶

この方法ではFTCSスキームの右辺の $u_j^n$ を $(u_{j+1}^n+u_{j-1}^n)/2$ で置き換え、以下のように差分化する。

(34)¶ $u_j^{n+1}=\displaystyle{1 \over 2}(u_{j+1}^n+u_{j-1}^n) -\displaystyle{\nu \over 2}(u_{j+1}^n-u_{j-1}^n)$

$u_j^n=g^n {\rm exp}(i j \theta)$ を代入して増幅率 $g$ を求めると

(35)¶ $g = \displaystyle{1 \over 2}(e^{i\theta}+e^{-i\theta})- \displaystyle{1 \over 2}(e^{i\theta}-e^{-i\theta})\\ = {\rm cos} \theta - i \nu {\rm sin} \theta$

したがって

(36)¶ $|g|^2={\rm cos}^2{\theta}+\nu^2 {\rm sin}^2 \theta$

Lax-Friedrichスキームの場合の増幅率¶

Lax-Friedrichスキームの場合の増幅率に増幅率 $|g|$ を $\theta$ の関数として極座標 $(g,\theta)$ で示す。Lax-Friedrich のスキームでは、クーラン数 $\nu=c\Delta t/\Delta x$ が $|\nu| \le 1$ を満たす場合、安定に計算を進めることができる。この条件のことを Courant, Friedrich, Lewy 条件 ( CFL条件 あるいはクーラン条件)と言う。

Lax-Friedrichスキームにおける依存関係。破線と点線はそれぞれ $\nu < 1$ 、 $\nu > 1$ の場合の波の伝播を示す。¶

クーラン条件の意味を考えてみよう。差分式 (34) を書き換えると以下の式を得る。

(37)¶ $u_j^{n+1}={{1-\nu} \over 2}u_{j+1}^n+{{1+\nu} \over 2}u_{j-1}^n$

クーラン条件 $|\nu|\le 1$ が満たされている場合、時刻 $t=t_{n+1}$ の値 $u_j^{n+1}$ は $t=t_n$ の $j-1$ 点の値 $u_{j-1}^n$ と $j+1$ 点の値 $u_{j+1}^n$ の内挿値になっている。上図にLax-Friedrichスキームにおける変数の依存関係を示す。実線は時刻 $t=t_{n+1}$ の白丸の格子点の値を計算する際に用いられる時刻 $t=t_n$ の格子点、破線は $\nu=c\Delta t/\Delta x < 1$ の場合、点線は $\nu > 1$ の場合の波の伝播を示す。クーラン条件は $|c|\Delta t < \Delta x$ 、すなわち時間間隔 $\Delta t$ の間に波が１メッシュ以上伝わってはいけないことを意味する。 $u_j^{n+1}$ は $u_{j-1}^n$ と $u_{j+1}^n$ だけから計算されるが、時間間隔が $\Delta t > \Delta x/|c|$ となると $x_{j-1} \le x \le x_{j+1}$ より外側からも情報が伝わってくるため計算を安定に進めることができなくなるのである。

Lax-Friedrichスキームを用いた１次元線形スカラー移流問題のシミュレーション結果。クーラン数 $\nu=0.25$ で50ステップ、100ステップ計算した結果を示す。¶

上図にLax-Friedrichスキームを用いて１次元線形スカラー移流方程式の解を求めた結果を示す。数値振動のない解が得られている。Lax-Friedrichスキームの欠点は数値散逸が大きく、不連続面が時間とともになまってしまうことである。

1次精度風上差分法¶

図のように波が正の方向に伝わっている場合を考える。このとき、 $j$ 点での空間微分を、 $j$ 点と風上にあたる $j-1$ 点の間の差分で近似する方法が風上差分である。

左図：右方向に伝わる波、右図：1次精度風上差分における依存関係。破線は波による情報の伝達を示す。¶

1次元スカラー移流方程式を時間については前進差分、空間については風上差分として差分化すると、 $c > 0$ の場合、以下の差分式を得る。

(38)¶ $\displaystyle{{u_j^{n+1}-u_j^n} \over {\Delta t}} + c \displaystyle{{u_j^n-u_{j-1}^n} \over \Delta x} = 0$

したがって

(39)¶ $u_j^{n+1} = u_j^n - c \displaystyle{{\Delta t} \over {\Delta x}} (u_j^n - u_{j-1}^n)$

右図に1次精度風上差分における変数の依存関係を示す。黒丸は時刻 $t_{n+1}$ の白丸の格子点の値を計算する際に用いられる時刻 $t_n$ の格子点、破線は時刻 $t_{n+1}$ に白丸の格子点に到達する波の伝播を示す。

増幅率は

(40)¶ $g = 1 - \nu (1-e^{-i\theta})\\ =(1-\nu+\nu {\rm cos} \theta)-i\nu {\rm sin} \theta$

したがって

(41)¶ $|g|^2 = (1-\nu+\nu {\rm cos} \theta)^2+\nu^2 {\rm sin}^2 \theta\\ = 1-2\nu(1-\nu)(1-{\rm cos} \theta)$

これより、 $0 \le \nu \le 1$ の場合、任意の $\theta$ について $|g| \le 1$ であり、安定であることがわかる。

風上差分の差分式 (39) は次のようにも書ける。

(42)¶ $u_j^{n+1}=(1-\nu)u_j^n+\nu u_{j-1}^n .$

クーラン条件 $0 \le \nu \le 1$ が満たされている場合、 $u_j^{n+1}$ は時刻 $t=t_{n+1}$ に $j$ 番目の格子点に到達する波の $t=t_n$ での位置（右図の破線矢印の出発点）における値 $u(x_j-c\Delta t,t_n)$ を $u_{j-1}^n$ と $u_j^n$ から線形内挿した値になっている。

1次精度風上差分法を用いた１次元線形スカラー移流問題のシミュレーション結果。左：クーラン数 $\nu=0.25$ で50ステップ、100ステップ計算した結果。右： $\nu=0.80$ で16ステップ、32ステップ計算した結果。¶

上図に１次精度風上差分法を用いたシミュレーション結果を示す。

練習問題¶

クーラン数 $\nu$ が1,0.75,0.5 の場合について1次精度風上差分の増幅率 $|g|$ を位相 $\theta$ の関数として求め、極座標 $(|g|,\theta)$ でプロットせよ。
1次精度風上差分法の差分式 (39) にテイラー展開を適用することによって、以下の偏微分方程式が得られることを示せ。右辺第1項が拡散項であることに注意して、クーラン条件を導け。

$\displaystyle{{\partial u} \over {\partial t}} +c\displaystyle{{\partial u} \over {\partial x}} = \displaystyle{1 \over 2} c\Delta x (1-\nu) \displaystyle{{\partial^2 u} \over {\partial x^2}} -\displaystyle{1 \over 6} c (\Delta x)^2(2\nu^2-3\nu+1) \displaystyle{{\partial^3 u} \over {\partial x^3}} + ...$

Lax-Wendroffのスキーム¶

Lax-Wendroffスキームはテイラー展開にもとづく差分法であり、以下のようにして導かれる。

(43)¶ $u_j^{n+1} = u_j^n+\Delta t \displaystyle{{\partial u} \over {\partial t}} + \displaystyle{1 \over 2}\Delta t^2 \displaystyle {{\partial^2 u} \over {\partial t^2}} + O(\Delta t^3)$

右辺第2項、第3項に $\partial u/\partial t = -c \partial u/\partial x$ 、 $\partial^2 u/\partial t^2=c^2 \partial^2 u/\partial x^2$ を代入すると

(44)¶ $u_j^{n+1}=u_j^n-c\Delta t \displaystyle{{\partial u} \over {\partial x}}+\displaystyle{1 \over 2}c^2 \Delta t^2 \displaystyle{{\partial^2 u} \over {\partial x^2}} + O(\Delta t^3)$

空間微分 $\partial u/\partial x$ 、 $\partial^2 u/\partial x^2$ をそれぞれ中心差分で近似すると

(45)¶ $u_j^{n+1} = u_j^n - \displaystyle{1 \over 2}c\Delta t \displaystyle{{u_{j+1}^n-u_{j-1}^n} \over {\Delta x}} + \displaystyle{1 \over 2}c^2 \left(\displaystyle{{\Delta t} \over {\Delta x}} \right )^2 (u_{j+1}^n-2u_j^n+u_{j-1}^n)$

これがLax-Wendroffスキームである。以上の導出過程からわかるように、Lax-Wendroffスキームは空間、時間についていずれも2次精度の解法になっている。

Lax-Wendroffスキームの安定性をvon Neumannの方法で調べてみる。増幅率は

(46)¶ $g = 1 - \displaystyle{\nu \over 2}(e^{i\theta}-e^{-i\theta}) +\displaystyle{\nu^2 \over 2}(e^{i\theta}-2+e^{-i\theta})\\ = 1 -i\nu {\rm sin} \theta + \nu^2 {\rm cos} \theta -\nu^2$

したがって、

(47)¶ $|g|^2 = [1-\nu^2(1-{\rm cos}\theta)]^2 + \nu^2 {\rm sin}^2 \theta\\ = 1 - 2\nu^2(1-\nu^2)(1-{\rm cos}\theta)$

これより、 $|\nu| \le 1$ であれば任意の $\theta$ について $|g| \le 1$ であり、安定であることがわかる。

2段階Lax-Wendroffスキームにおける依存関係。第一段階で時刻 $t_n$ における格子点 $j-1$ 、 $j$ 、 $j+1$ の値から時刻 $t_{n+1/2}$ における格子点 $j-1/2$ 、 $j+1/2$ の値が計算される。第二段階ではこれらの点の値を用いて時刻 $t_{n+1}$ の白丸の格子点の値が求まる。¶

1次元スカラー方程式の場合、Lax-Wendroffスキームは以下のように2段階に分けたスキームと同等である。この方法を 2段階Lax-Wendroff法 と呼ぶ。

(48)¶ $u_{j+1/2}^{n+1/2} = \displaystyle{{u_{j+1}^n+u_j^n} \over 2} - {1 \over 2} c \displaystyle{{\Delta t} \over {\Delta x}} (u_{j+1}^n - u_j^n)$

(49)¶ $u_j^{n+1} = u_j^n - c \displaystyle{{\Delta t} \over {\Delta x}} (u_{j+1/2}^{n+1/2}-u_{j-1/2}^{n+1/2})$

これを図示すると上図のようになる。

2段階Lax-Wendroff法による１次元線形スカラー移流方程式のシミュレーション結果。左図：クーラン数 $\nu=0.25$ で50ステップ、100ステップ計算した結果。右図：クーラン数 $\nu=0.80$ で16ステップ、32ステップ計算した結果。¶

図にLax-Wendroff法を用いて1次元線形スカラー移流方程式の数値解を求めた結果を示す。Lax-Wendroff法は空間、時間についていずれも2次精度の方法であるが、不連続面近傍で数値振動を生じるという欠点を持つ。これに関して、以下の定理が知られている。

Godunov の定理¶

1次元スカラー移流方程式 $\partial u/\partial t + c \partial u /\partial x = 0$ に対して、

(50)¶ $u_j^{n+1} = \sum_k a_k u_{j+k}^n$

の形の2次精度以上の精度を持つどのようなスキームも解の単調性を維持することはできない。

ここで、「解の単調性を維持する」とは、時刻 $t_n$ におけるプロフィール $u(x,t_n)$ が $x$ に関して単調増加または単調減少する関数であるならば時刻 $t_{n+1}$ における関数 $u(x,t_{n+1})$ も単調増加または単調減少関数でなければならないことを意味する。たとえば１次精度風上差分の場合、 $0 \le \nu \le 1$ なら $u_j^{n+1}$ は必ず $u_{j-1}^n$ と $u_j^n$ の間の値をとるため、もしも $u_{j-1}^n \le u_j^n \le u_{j+1}^n$ なら $u_j^{n+1} \le u_{j+1}^{n+1}$ となり、単調性が維持される。Godunovの定理の証明については、たとえば藤井 (1994) を参照されたい。

数値振動を抑える方法には以下のものがある。

人工粘性を加える

粘性係数を $\kappa$ として、

(51)¶ $\tilde{u}_j^{n+1} = u_j^{n+1} + \kappa \displaystyle {{\Delta t} \over {{\Delta x}^2}} (u_{j+1}^n - 2u_j^n + u_{j-1}^n)$

(52)¶ $\kappa_{j+1/2}= Q_v \Delta x |u_{j+1}^n-u_j^n|$

とする。拡散係数 $\kappa$ は、たとえば以下のように不連続面付近で大きな値をとるように決める。 $Q_v$ はパラメータである。

流束制限関数を用いる

これについては後述する。

保存形表示と数値流束¶

1次元スカラー移流方程式

(53)¶ $\displaystyle{{\partial u} \over {\partial t}} + c \displaystyle{{\partial u} \over {\partial x}} =0$

を以下の形に変形する。

(54)¶ $\displaystyle{{\partial u} \over {\partial t}} + \displaystyle{{\partial f} \over {\partial x}} = 0$

ここで、

(55)¶ $f = cu$

は流束をあらわす。式 (54) の形を保存形と呼ぶ。

メッシュ点とメッシュ境界を通って出入りする流束の関係¶

保存形式の物理的意味を考えるために、図に四角で囲って示した領域 ( $x_{j-1/2} < x < x_{j+1/2}$ )における保存量 $u$ の時間変化を求めてみよう。方程式 (54) を $x=x_{j-1/2}$ から $x=x_{j+1/2}$ まで積分すると次式を得る。

(56)¶ $\displaystyle{\partial \over {\partial t}} \int_{x_{j-1/2}}^{x_{j+1/2}} u dx + f(x_{j+1/2})-f(x_{j-1/2}) = 0 .$

したがって、保存量 $u$ の積分量

(57)¶ $u_j^n = \int_{x_{j-1/2}}^{x_{j+1/2}} u(x,t_n) dx$

の時間変化は、この時間の間に左右の境界 $x_{j \pm 1/2}$ を通って出入りする流束 $f_{j \pm 1/2}$ の差に等しい。これより次式を得る。

(58)¶ $u_j^{n+1} = u_j^n - \displaystyle{{\Delta t} \over {\Delta x}} (f_{j+1/2}^n - f_{j-1/2}^n)$

差分式 (58) は保存則を厳密に満たす。これが保存形式を用いる利点である。

メッシュ境界の流束 $f_{j \pm 1/2}^n$ は各メッシュ点での流束から近似的に計算することができる。これを 数値流束 と言い、 $\tilde{f}_{j \pm 1/2}^n$ であらわす。各種差分スキームの差分式から数値流束を求めると以下のようになる。

FTCS

(59)¶ $\tilde{f}_{j+1/2}^n = \displaystyle{1 \over 2}(f_{j+1}^n+f_j^n)$

Lax-Friedrich

(60)¶ $\tilde{f}_{j+1/2}^n = \displaystyle{1 \over 2} \left [(1-\displaystyle{1 \over \nu})f_{j+1}^n+(1+\displaystyle{1 \over \nu}) f_j^n \right ]$

Upwind (風上差分)

(61)¶ $\tilde{f}_{j+1/2}^n = \displaystyle{1 \over 2} \left [(f_{j+1}^n + f_j^n)-|c|(u_{j+1}^n-u_j^n) \right ]$

: この式は、 $c > 0$ の場合は $\tilde{f}_{j+1/2}^n= f_j^n$ 、 $c < 0$ の場合は $\tilde{f}_{j+1/2}=f_{j+1}^n$ と一致する。

Lax-Wendroff

(62)¶ $\tilde{f}_{j+1/2}^n = \displaystyle{1 \over 2} \left [ (1-\nu)f_{j+1}^n+(1+\nu)f_j^n \right ]$

Burgers方程式の数値解法¶

ここまでは、移流の速さ $c$ が一定の場合の１次元線形スカラー移流方程式を扱ってきた。本節では、以下のような 非線形 波動方程式を差分近似によって解くことを考える。これは、非粘性の場合のBurgers方程式である。

(63)¶ $\displaystyle{{\partial u} \over {\partial t}} + u \displaystyle{{\partial u} \over {\partial x}} = 0$

この方程式は流線に沿うラグランジュ微分 $d/dt=\partial /\partial t+ u \partial /\partial x$ を用いると次のように表現できる。

(64)¶ $\displaystyle{{du} \over {dt}} = 0$

Burgers方程式の解の様子。ある有限の時刻で後ろからきた粒子が前の粒子に追いつく。¶

粒子的な描像に立てば、この方程式は力を受けていない粒子の運動を記述しており、その解はもちろん $u=const.$ である。初期速度分布が正弦波的な場合、上図に示すように、振幅が正の領域は $+x$ 方向に、負の領域は $-x$ 方向に移動してしだいに波が突っ立ち、有限の時刻で後からきた粒子が前の粒子に追いついてしまう。連続系では空間の１点で速度が多価になることはできないため、このような場合に解に不連続が生じる。

Burgers方程式を１次精度風上差分法で解いた結果の例。左：初期に $u>0$ の場合。右：初期に $u<0$ の場合。図中の数字は時間ステップ数。時間きざみは $\Delta t/\Delta x=0.8$ とした。¶

Burgers方程式を差分法によって解くため、方程式をまず、以下のような保存系に変形する。

(65)¶ $\displaystyle{{\partial u} \over {\partial t}} + \displaystyle{{\partial} \over {\partial x}} \left(\displaystyle {{u^2} \over 2} \right)=0$

これは、流束 $f(x)$ が $f(x)=u^2/2$ の場合に相当し、各種差分スキームを適用することができる。たとえば１次精度の風上差分法を適用する場合、線形スカラー方程式で移流の速さ $c$ が一定である場合には $c > 0$ のとき $f_{j+1/2}^n=f_j^n$ であったことに注意し、メッシュ境界の $j+1/2$ 点での速さを $(u_j(t)+u_{j+1}(t))/2$ で近似すると、

$u_{j+1}(t)+u_j(t) > 0$ のとき、 $f_{j+1/2}^n = f_j^n=|u_j(t)|^2/2$
$u_{j+1}(t)+u_j(t) \le 0$ のとき、 $f_{j+1/2}^n = f_{j+1}^n=|u_{j+1}(t)|^2/2$

図に初期に $u(x)=1+\epsilon {\rm sin}(kx)$ ( $\epsilon=0.01$ , $0 \le kx \le 2\pi$ ) のような速度分布を与えた場合のBurgers方程式の解を1次精度の風上差分法で計算した結果を示す。左図の場合、初期に $u \sim 1$ であることから予想できるように非線形効果が小さい間の解は波の速さが $c=1$ の場合の線形スカラー移流方程式の解とほぼ一致し、波はほぼその形を保ちながら右側に伝わっていく。右図では初期に $u \sim -1$ であり、波は左に伝わる。

初期に $u(x)=1+0.1 {\rm sin}(kx)$ の速度分布から始めた場合のBurgers方程式の数値解。図中の数字は時間ステップ数。１次精度風上差分法で時間きざみは $\Delta t/\Delta x=0.8$ とした。¶

非線形性が強くなる場合のBurgers方程式の数値解の例を図に示す。この例では初期に $u(x)=1+0.1 {\rm sin}(kx)$ のような速度分布を与え、その後の時間発展を１次精度の風上差分法によって解いた。Burgers方程式の非線形項 $u\partial u/\partial x$ の効果により波がしだいに突っ立ち、不連続（衝撃波）が形成されることがわかる。

流束制限関数¶

以上、線形スカラー移流方程式とBurgers方程式を例にして差分解法について解説してきた。１次精度の風上差分法を用いるとこれらの方程式の解にあらわれる不連続面を数値振動を起こすことなくとらえることができる。しかしながら、Godunovの定理が示すように、空間２次精度以上の解法では数値振動があらわれてしまうことがわかった。

Lax-Wendroff法の数値流束を補正することによって、不連続面近傍での振動を抑えることができないかどうか考えてみよう。Lax-Wendroff法の数値流束は次のようにも書ける。

(66)¶ $\tilde{f}_{j+1/2}^n = c[u_j^n + \displaystyle{1 \over 2} (1-\nu)(u_{j+1}^n - u_j^n)]$

数値振動が生じない1次精度の風上差分の数値流束は $c > 0$ のとき $\tilde{f}_{j+1/2}^n=cu_j$ であり、Lax-Wendroff法の数値流束の右辺第1項と一致している。そこで、Lax-Wendroff法の数値流束の右辺第2項を次のように補正してみる。

(67)¶ $\tilde{f}_{j+1/2}^n = c[u_j^n + \displaystyle{1 \over 2} (1-\nu)B_{j+1/2}(u_{j+1}^n - u_j^n)]$

ここで導入した $B_{j+1/2}$ のことを 流束制限関数 と呼ぶ。数値流束 (67) を差分式 (58) に代入して変形すると次式を得る。

(68)¶ $\displaystyle{{u_j^{n+1}-u_j^n} \over {u_{j-1}^n-u_j^n}} = \nu [1-\displaystyle{1 \over 2}(1-\nu)B_{j-1/2}] +\displaystyle{1 \over 2}\nu(1-\nu) \displaystyle {{B_{j+1/2} \over r_j}}$

ここで、

(69)¶ $r_j \equiv \displaystyle{{u_j^n-u_{j-1}^n} \over {u_{j+1}^n-u_j^n}}$

である。

数値振動が生じないようにするために $u_j^{n+1}$ の値を制限する範囲¶

数値振動が生じないようにするために、図に示したように $u_{j}^{n+1}$ が $u_j^n$ と $u_{j-1}^{n}$ の間の値をとるように制限を加えることにしよう。これには、式 (68) の左辺の値を以下のように制限すればよい。

(70)¶ $0 \le \displaystyle{{u_j^{n+1}-u_j^n} \over {u_{j-1}^n-u_j^n}} \le 1$

式 (68) の右辺を代入すると以下の条件を得る。

(71)¶ $-\displaystyle{2 \over \nu} \le B_{j-1/2} - \displaystyle {{B_{j+1/2}} \over {r_j}} \le \displaystyle{2 \over {1-\nu}}$

CFL条件が満たされている場合 $0 \le \nu \le 1$ なので、

(72)¶ $-2 \le B_{j-1/2}-\displaystyle{{B_{j+1/2}} \over {r_j}} \le 2$

この関係式は、以下のふたつの条件がともに満たされれば成立する。

(73)¶ $0 \le B_{j+1/2} \le 2$

かつ

(74)¶ $0 \le \displaystyle{{B_{j+1/2}} \over {r_j}} \le 2$

この範囲を図示すると流束制限関数の許容範囲図の斜線のない領域になる。

流束制限関数 $B_{j+1/2}(r)$ の許容範囲。LWはLax-Wendroffスキームの数値流束に対応する制限関数。minmodは minmod関数。¶

$r < 0$ の場合は $B_{j+1/2}=0$ のみが許される。Lax-Wendroff法の数値流束では $B_{j+1/2}=1$ (図のLW)であるため、 $r < 1/2$ の領域で許容範囲外となり、数値振動が生じる。図の許容範囲内にある流束制限関数を用いることにより、数値振動が起こらないようにすることができる。その一例は以下のminmod関数(図のminmod)である。

(75)¶ ${\rm minmod}(r) = \left\{ \begin{array}{ll} 0 & (r < 0)\\ r & (0 \le r \le 1)\\ 1 & (r > 1) \end{array} \right.$

TVDスキーム¶

前節では、数値振動をおさえる方法として流束制限関数を導入した。ここでは、数値振動の発生を定量化する方法について考える。

このために、1次元線形スカラー移流方程式において以下の量を定義する。

(76)¶ $U = \int \left |\displaystyle{{du} \over {dx}}\right | dx .$

この量は波の振幅の総和に等しく、移流方程式の厳密解では波のプロフィールが保たれるため、 $dU/dt = 0$ である。

以上との類推により、メッシュ点ごとの物理量の変化量の総和を次式のように定義し、これをTotal Variation (TV)と言う。

(77)¶ $TV(u^n) \equiv \sum_j |u_{j+1}^n-u_j^n|$

Total Variation が時間とともに増大しないという条件

(78)¶ $TV(u^{n+1}) \le TV(u^n)$

のことをTotal Variation Diminishing (TVD)条件と呼ぶ。

流束制限関数を導入することによって、差分スキームがTVD条件を満たすようにすることができる。

放物型方程式の差分解法¶

天体シミュレーションにあらわれる放物型方程式¶

熱伝導方程式

$\centerline{$\displaystyle{{\partial T} \over {\partial t}}= \kappa {\bf \nabla}^2 T$}$

磁気拡散方程式

$\centerline{$\displaystyle{{\partial \boldmath{B}} \over {\partial t}}= \eta {\bf \nabla}^2 \boldmath{B}$}$

以下のような1次元拡散方程式を差分近似によって初期値問題として解くことを考えてみよう。すなわち、時刻 $t=0$ における $u(x,t)$ の値 $u(x,0)$ を与えて、任意の時刻 $t$ ( $> 0$ )における $u(x,t)$ を求める。

(79)¶ $\displaystyle{{\partial u} \over {\partial t}} = \kappa \displaystyle{{\partial^2 u} \over {\partial x^2}}$

拡散係数 $\kappa$ は x に依らないとする。よく知られているように、この方程式の解は初期条件をフーリエ変換することによって解析的に求めることができる。

拡散方程式の解の時間発展の様子¶

解のおおまかな様子を図に示す。

拡散方程式の陽解法 (explicit法)¶

1次元拡散方程式 (79) を時間について現在の時刻 $t_n$ と $\Delta t$ 後の時刻 $t_{n+1}=t_n+\Delta t$ の間で前進差分、空間については中心差分（FTCS差分）をとって差分化すると次式を得る。ここで、 $u_j^n = u(x_j,t_n)$ である。

(80)¶ $\displaystyle{{u_{j}^{n+1}-u_{j}^{n}} \over {\Delta t}} = \kappa \displaystyle{{u_{j+1}^{n}-2u_j^n+u_{j-1}^n} \over {\Delta x^2}}$

式 (80) を変形して次式を得る

(81)¶ $u_j^{n+1} = u_j^n + \displaystyle{{\kappa \Delta t} \over {\Delta x^2}} (u_{j+1}^n - 2u_j^n + u_{j-1}^n)$

右辺は時刻 $t_n$ での値、左辺は時刻 $t_{n+1}=t_n+\Delta t$ での値だけで書けている。したがって、時刻 $t_n$ での各格子点での値がわかっていれば直ちに1タイムステップ後( $t_{n+1}$ )の各格子点での値を計算することができる（陽解法）。

Von Neumannの安定性解析¶

FTCSスキームの数値的安定性を調べるために $u_j^n = g^n {\rm exp} (i j \theta)$ をFTCS差分式に代入すると

(82)¶ $g^{n+1}e^{i j \theta} = g^n e^{i j \theta} + \displaystyle{{\kappa \Delta t} \over {\Delta x^2}} g^n [e^{i (j+1) \theta} - 2 e^{i j \theta} + e^{i (j-1) \theta}].$

よって、

(83)¶ $g = 1 - 2 \displaystyle{{\kappa \Delta t} \over {\Delta x^2}} (1-{\rm cos} \theta).$

増幅率が $|g| \le 1$ であるためには

(84)¶ $-1 \le 1 - 2 \displaystyle{{\kappa \Delta t} \over {\Delta x^2}} (1-{\rm cos} \theta) \le 1 .$

したがって、

(85)¶ $0 \le \displaystyle{{\kappa \Delta t} \over {\Delta x^2}} (1-{\rm cos} \theta) = \displaystyle{{\kappa \Delta t} \over {\Delta x^2}} 2{\rm sin}^2 {\theta \over 2} \le 1 .$

任意の $\theta$ (任意の波長の波)について安定であるためには

(86)¶ $0 \le \displaystyle{{\kappa \Delta t} \over {\Delta x^2}} \le \displaystyle{1 \over 2} .$

以上により、FTCSスキームにより1次元拡散方程式のシミュレーションを行う場合、時間ステップ $\Delta t$ が上式を満たすようにコントロールする必要があることがわかる。たとえばメッシュサイズを半分にした場合、 $\Delta t$ は1/4にしなければならない。

拡散方程式の陰解法 (implicit法)¶

拡散方程式を差分化する際に右辺の空間差分の部分に、求めるべき $t_{n+1}$ での $u$ の値を含めて差分化する方法がある。このような方法を 陰解法（implicit法） と呼び、explicit法とは安定性条件が異ってくる。代表的な陰解法である Crank-Nicolson 法では、パラメータ $\lambda$ を導入して、以下のように差分化する。

(87)¶ $\displaystyle{{u_{j}^{n+1}-u_{j}^{n}} \over {\Delta t}} = \kappa \left [ \lambda \displaystyle{{u_{j+1}^{n+1}-2u_j^{n+1}+ u_{j-1}^{n+1}} \over {\Delta x^2}} + (1-\lambda) \displaystyle{{u_{j+1}^{n}-2u_j^{n}+ u_{j-1}^{n}} \over {\Delta x^2}} \right ]$

これを整理すると次のような行列を含む式になる。

(88)¶ $\boldmath{A} \boldmath{u}^{n+1} = \boldmath{b}(\boldmath{u}^n) .$

これを解いて $\boldmath{u}^{n+1}$ を求めればよい。

練習問題¶

行列 $\boldmath{A}$ とベクトル $\boldmath{b}$ の各要素を求めなさい。
Von Neumann の安定性解析により、 $\lambda > 1/2$ ならば $\kappa \Delta t/\Delta x^2 > 0$ を満たす任意の $\Delta t$ についてCrank-Nicolsonスキームは数値的に安定であることを示しなさい。

参考文献¶

[1]は、数値流体力学全般についてまとめられたテキストであり、必読文献である。