スペクトル法の一部の基礎の初歩への はじめの一歩 高橋芳幸
はじめに 偏微分方程式の数値解法としては, これまでに様々なものが考案されている. ここではスペクトル法に注目し, その「一部の基礎の初歩の概要を大雑把に」説明したい. なお, ここでは非常に簡単な例を「多少誤魔化しつつ」述べているので, 正確な記述と詳細は様々な文献をあたっていただきたい.
参考文献 石岡 (2004), スペクトル法による数値計算入門, 東京大学出版会. Durran (1999), Numerical Methods for Wave Equations in Geophysical Fluid Dynamics, Springer-Verlag. Haltiner and Williams (1980), Numerical prediction and dynamic meteorology - 2nd ed., Wiley.
題材とする方程式 ここでは例として以下に示す 1 次元移流方程式を離散化を扱う. (1) ただし, ここでは c は定数で, 0 x 2 とし, 周期境界条件とする.
有限差分法 (参考までに) (1)式は, 有限差分法を用いると以下のように離散化される. 例えば, 2 次精度の中心差分で離散化すると以下のようになる.
スペクトル法 (1) スペクトル法では, 解を有限個の直交関数で展開し, その係数を求める. ここでは, 展開関数としてフーリエ級数を用いる. 注:以下では, 非常に簡単な系の話を, 詳細を飛ばして説明しています. 本当はもっとちゃんとした議論がありますので詳細は文献をあたってください. 解を以下のように展開することにする. ここで M は切断波数. (1) 式に代入すると, となり, 常微分方程式になった. (2)
スペクトル法 (2) 時間微分を有限差分法で近似し, 整理すると となる. 実空間での解の構造を求めるには, 離散フーリエ逆変換すればよい. ( (2) 式に代入すればよい.)
スペクトル法 (3) 実空間データ 離散フーリエ変換 時間積分 離散フーリエ逆変換 実空間データ
スペクトル法 (4) スペクトル法の利点 スペクトル法の欠点 微分値の見積もり精度が差分法に比べて非常に高い. 格子点配置に任意性がない. 滑らかな関数を展開関数に用いると, 物理量の急激な変化を表すことができない. 境界条件によっては適した展開関数がない.
スペクトル法 (5) 線型方程式 非線型方程式 非線型方程式を考える際には工夫が必要. 波数空間においても方程式がとても簡単. 計算量が少ない. 非線型方程式 波数空間での計算が”若干のややこしい”. (何も考えずに)計算すると計算量が多い. 非線型方程式を考える際には工夫が必要.
スペクトル法 (6) 以下の非線型方程式を考える. 解を以下のように展開することにして上の式に代入すると, より,
スペクトル法 (7) となる.
スペクトル法 (8) 非線形項の計算量に関する考察. 右辺の計算量は O(M2). 計算量が多すぎるので, 減らしたい.
スペクトル法 (9) 変換法 非線形項(の掛け算)を実空間で計算(微分は波数空間で評価). 変換法(非線形項を実空間で求め逆変換)を用いると計算量は O(M log M)
スペクトル法 (10) 実空間データ 波数空間で微分の評価 離散フーリエ変換 離散フーリエ逆変換 時間積分 実空間での非線型項の評価(掛け算) 離散フーリエ逆変換 離散フーリエ変換 時間積分 実空間データ
スペクトル法 (11) 変換法を用いた場合の計算量 FFT を用いると, 何も考えない場合の計算量(O(M2)) よりも少なくてすむ. 実質的には離散フーリエ(逆)変換の計算量で決まる 離散フーリエ変換の計算量は O(M2) 高速フーリエ変換(Fast Fourier Transform: FFT)を用いると, 計算量は O(MlogM) FFT を用いると, 何も考えない場合の計算量(O(M2)) よりも少なくてすむ.
スペクトル法 (12) spmodel の特徴(の一部) 離散フーリエ変換(を含む各種変換)が簡単. 微分も簡単 離散フーリエ変換 :関数 e_g e_data = e_g( g_data ) 離散フーリエ逆変換 :関数 g_e g_data = g_e( e_data ) 微分も簡単 x 微分 g_dfdx = g_Dx_e( e_g( g_data ) )
スペクトル法 (13) 実空間データ 波数空間で微分の評価 離散フーリエ変換 離散フーリエ逆変換 時間積分 実空間での非線型項の評価(掛け算) 離散フーリエ逆変換 離散フーリエ逆変換 g_psi * g_Dx_e( e_psi ) 時間積分 実空間データ
スペクトル法 (14) 格子点数 N の時, 最初に考える切断波数の選択は M=N/2 変換法を用いる際にはエイリアシングに注意が必要. 格子点数 N の時, 最初に考える切断波数の選択は M=N/2 実空間でのデータ数(自由度)と波数空間でのデータ数が等しい. 波数 M までの成分を含む 2 項の積は, 波数 2M までの成分を含む. 実際には折り返されて“偽の”成分が生じる(エイリアシング). 畳み込み積分の結果は, 波数 M までしか含まない. そこで, エイリアシングが起こらない程度まで切断波数を小さくする. M = (N-1)/3 とすると, 2 項までの積を含む系ではエイリアシングは起こらない. M 2M 波数 N/2 k1+k2 k’
参考文献 石岡 (2004), スペクトル法による数値計算入門, 東京大学出版会. Durran (1999), Numerical Methods for Wave Equations in Geophysical Fluid Dynamics, Springer-Verlag. Haltiner and Williams (1980), Numerical prediction and dynamic meteorology - 2nd ed., Wiley.