4. 双曲型偏微分方程式の数値解法入門 双曲型の偏微分方程式(partial differential equation, PDE)の最も簡単なの例として1変数の線形PDE を考える; この方程式の意味は大雑把に言って、Δx の セル内に流入流出する f の量がフラックス その結果セル内で f が時間的に変化する割合が Δx この方程式は f (t-x/c) を解に持つ。 後の議論のために、放物型の偏微分方程式簡単な例として1変数の拡散方程式 も思い出しておこう。 αは定数である。
有限差分法 ( Finite difference method,PDEの数値解法の一つ) x t
線形波動方程式 の陽解法の幾つかの例 (1) FTCS (Forward in Time and Central difference Scheme,陽的オイラー法). に代入して、 x t
t x 影響領域 Domain of dependence 陽的 Euler 法では、 の値を計算するのに、点線のピラミッドの形の下の この数値スキームの影響領域は物理的な 影響領域を含んでいる必要がある。 また,数値スキームの影響領域は差分法が異なると違ったものになる。
線形波動方程式 の陽解法の幾つかの例 (1) FTCS (Forward in Time and Central difference Scheme,陽的オイラー法). (2) Lax (- Friedrich) scheme (ラックス-フリードリッヒ). (3) Leap-Flog scheme (リープフロッグ 蛙飛び) (4) Lax-Wendroff scheme (ラックス-ウェンドロフ) (5) 1st order upwind scheme,1次の風上差分
線形波動方程式 の陽解法の例 - 別の見方 Lax, Lax-Wendroff, 1次風上の各スキームは、FTCS scheme +.拡散項 の形で解釈することが出来る。 cf) 拡散方程式 の差分スキーム(陽的オイラー法) ここで、タイムステップ t と格子点間隔 h を適当に調整して、 となるように取ってみる。すると差分スキームは となる。これは、 の値を計算するのに 前の時間ステップでの左右の格子点での 値を平均していることになっている。 → 平均すると関数の値がだんだん鈍っていって 拡散していくように見える。
線形波動方程式 の陽解法の幾つかの例 (続き) Lax, Lax-Wendroff, 1次風上の各スキームは、FTCS scheme +.拡散項 の形で解釈することが出来る (1) 陽的オイラー法 (FTCS : 時間前進差分、空間中心差分). (2) Lax (- Friedrich) 法. (3) Lax-Wendroff 法. (4) 1次風上差分法 1st order upwind scheme (拡散項が一番小さい) (負の特性速度cの場合にもそのまま使える書き方。)
精度から要求されるステップサイズt とhに対する制限 の差分法の記法のひとつとして、
PDEの数値解法の基本概念の要約 (ODEの場合と類似している) 局所打切り誤差(Local truncation error) : PDEの厳密解が有限差分方程式を満たしていない度合い。 ex.) One-step method. 大域差分誤差(Global discretization error) Definitions: (適合性 Consistency) Definitions: (収束性 Convergence)
Definitions: (ゼロ安定性 Zero-Stability, h ! 0 の時の安定性条件) Definitions: (絶対安定性 - h が有限の値に固定された時の安定性条件) Remark: 安定性解析の目的は、摂動が増大しないことを保証する t0(h) の値を決めることである。摂動が増大しないためには、 ゼロ安定の場合には,行列 C0(h) のサイズは n ! 1 の極限で 増大することもあることに注意せよ。 Theorem: (Lax’s equivalence theorem) Convergence ) Zero (or Absolute) - stability. Zero (or Absolute) - stability and Consistency ) Convergence.
有限体積差分法 保存法則(保存形をした PDE)の有限差分近似を導くのに適した考え方。 j–1/2 j+1/2 j–1 j j+1 格子点間の中点を境界面にした仮想的な単位体積を考える。 境界面、 j–1/2, j+1/2 , で流束(flux)を定義する。そして、PDEを時間方向に陽的オイラー法で差分化する ただし、境界面 j-1/2, j+1/2… には格子点がない(データがない)! 格子点での値 を使って流束 を計算する。 陽的Euler, Lax, Lax-Wendroff, 1st order upwind 法を を用いて書き直す。 (1) 陽的 Euler :
(1) 陽的 Euler : (2) Lax : (3) Lax-Wendroff : (4) 1st order upwind :
ノイマンの安定性解析 von Neumann stability analysis. 等間隔格子点と周期的境界条件を空間方向に仮定した下で、線形PDE のを解くための数値スキームの安定性を解析する方法。 有限差分方程式が次の様な解を持っている場合を考える。 この場合、摂動 も同じく の様に書ける。 PDEの有限差分法の一段階公式(one-step formula)は、h とτに依存する 行列 C0 を用いて一般に ,従って摂動 も 同様に を満たす。 この式に上のフーリエ級数展開を代入すると、 となる。 係数 を増幅因子と呼ぶ。初期値 t = 0 まで戻れば、この式は とも書ける。
ノイマンの安定性解析 von Neumann stability analysis. 続き 増幅因子 は一般に複素数なので、振幅と位相にわけてやる。 増幅率 位相差 PDE の解析解は f (t-x/c) なので、厳密な増幅率と 位相差は、 これより、 と が一般に各モードの増幅率と位相差 を表す。 θが0に近いほうが、長波長、低周波モードに対応し、 πに近いほうが、短波長、高周波モードに対応する。
ノイマンのゼロ安定性( zero stability )条件: ノイマンの絶対安定性( absolute stability )条件: 線形のPDEに線形の差分法を適用した場合には、各モードは独立になる。従って、 が1以下ならばモードは増大しない。
(参考)離散フーリエ解析の公式 差分スキームのノイマンの安定性解析では離散フーリエ変換を利用した。 このために必要な公式を下にまとめる。 格子点を ととり、 のように周期境界条件を課す。独立な点は従って、 のN点になる。これに対応して、N個の独立な基底 を取ることが出来る。 これより、離散フーリエ変換を次の様に定義できる; 公式 より、逆変換は 実空間と波数空間での l2 ノルムの値の間には Perceval の公式が成り立つ;
課題4-1) Lax-Wendroff 法が2次精度の有限差分スキームになっていることを示せ。 課題4-2) FTCS, Lax, Lax-Wendroff, 1st order upwind, の各有限差分法で、 線形PDE を次の初期条件の場合に解け。 各数値スキームにより、シミュレーションの結果の様子が異なる。 この理由を考察せよ。また、初期条件を変えてシミュレーション を実行してみよ。(課題の解説を参照せよ。) 課題4-3) FTCS, Lax, Lax-Wendroff, 1st order upwind, の差分法について ノイマン解析を行い、シミュレーションの結果と比べてみよ。