Next: 再度数値流束 Up: 数値流体力学の基礎的事項 Previous: 空間１次精度のスキーム

安定性

前節で見たように、空間２次精度のスキームでは、余分な振動が発生した。空間１次精度のスキームは余分な振動は発生しなかったが、密度の飛びは時間がすすむとともに鈍ってゆくことがわかった。これらの差はどこからくるのだろうか。差分方程式に含まれる元の微分方程式からのずれ（誤差）に起因することは前に述べた。

ここでは、von Neumannの安定性解析を使って余分な振動や解の鈍りと差分スキームの間の関係を定量的に見てゆく。解をフーリエ分解して、

と書く。ここでjは格子点の番号で、は格子間隔を単位とした波数つまりで。

この式をスキームに代入しに比例する項を取り出す。 FTCSスキーム

に対しては、

となる。これから複素振幅gが

となっていることがわかる。

この問題では厳密解は、を満足していなければならないので、

が成り立たねばならず、これから厳密解に対する複素振幅gは

で与えられる。これは当然で、厳密解では位相だけが動いてゆくことになる。これからどれだけ異なっているかを調べれば、差分化にともなって、微分方程式の厳密解からずれた部分がわかることになる。

複素振幅gの絶対値はある波の１ステップ後との振幅の比を与えるが、式(1.45)に対しては

で与えられる。に対しては波数の波は増幅されることがわかる。

つぎに複素振幅gの位相は、

で、厳密解の位相に対する比は

のように計算できる。これを図に書くと図1.6のようになる。 |g|とを図示するMathematicaのプログラムは

G[nu_,theta_]=Sqrt[1+nu^2 Sin[theta]^2]
Plot[G[0.5,theta],{theta,0,Pi},Frame -> True]

および

Phi[nu_,theta_]=ArcSin[nu Sin[theta]/Sqrt[1+nu^2 Sin[theta]^2]]/(nu theta)
Plot[Phi[0.5,theta],{theta,0.0001,Pi},Frame->True]

である。

Figure 1.6: FTCSスキームで、クーラン数でにする複素振幅gの絶対値(左)と複素振幅gの位相の厳密解の位相に対する比（右）。

ここで、複素振幅gの絶対値が（0と以外）に対して 1を越えているので、初期値に含む周期的な波はかならず増幅されることになることがわかる。また、の大きな（短波長の）波の方で厳密解に比べて、位相が遅れることもわかる。図1.2で短波長の波が後方に取り残されているように見えることもこのように説明される。

同じようにして、Laxのスキームでは、

となり、絶対値と位相の比は、

で与えられる。これを図示するMathematicaのプログラムは

G[nu_,theta_]=Sqrt[Cos[theta]^2+nu^2 Sin[theta]^2]
Plot[G[0.5,theta],{theta,0,Pi},Frame -> True]

Phi[nu_,theta_]=ArcCos[Cos[theta]/G[nu,theta]]/(nu theta)
Plot[Phi[0.5,theta],{theta,0.0001,Pi},Frame->True]

のようになる。

Figure 1.7: Laxスキームで、クーラン数に対する複素振幅gの絶対値(左)と複素振幅gの位相の厳密解の位相に対する比（右）。

図1.7（左）を見ると、すべてので|g|<1であり、周期的な振動はその振幅が減少してゆくことがわかる。特に、つまり4グリッドで1波長となる波が大きく減衰することがわかる。図1.7（右）からは、の大きい波に対して、正の位相誤差とを生じており、Laxスキームで飛びが鈍ってゆくにつれて前方に崩れてゆく構造がこれが原因となっていると見ることができる。