Transcript スペクトル法の一部の基礎の初歩への はじめの一歩 高橋芳幸

スペクトル法の一部
の基礎の初歩への
はじめの一歩
高橋芳幸
はじめに
• 偏微分方程式の数値解法としては, これまで
に様々なものが考案されている. ここではス
ペクトル法に注目し, その「一部の基礎の初
歩の概要を大雑把に」説明したい.
• なお, ここでは非常に簡単な例を「多少誤魔
化しつつ」述べているので, 正確な記述と詳
細は様々な文献をあたっていただきたい.
参考文献
• 石岡 (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 次元移流方
程式を離散化を扱う.


 c
t
x
(1)
ただし, ここでは c は定数で, 0  x  2 とし, 周期境界条件とする.
有限差分法
（参考までに）
• (1)式は, 有限差分法を用いると以下のよう
に離散化される.
 nj  ( x j , tn ), x j  jx, tn  nt
例えば, 2 次精度の中心差分で離散化すると以下のようになる.
 nj1  nj1
2t

n1
j

n1
j
 c
 nj1  nj1
2x

t n
 c  j 1  nj1
x

スペクトル法 (1)
• スペクトル法では, 解を有限個の直交関数で展開し, その係数を求
める.
• ここでは, 展開関数としてフーリエ級数を用いる.
注：以下では, 非常に簡単な系の話を, 詳細を飛ばして説明していま
す. 本当はもっとちゃんとした議論がありますので詳細は文献をあ
たってください.
• 解を以下のように展開することにする.

M
ikx
a
(
t
)
e
k
k  M
ここで M は切断波数. (1) 式に代入すると,
dak (t )
 ick ak (t )
dt
となり, 常微分方程式になった.
(2)


 c
t
x
スペクトル法 (2)
• 時間微分を有限差分法で近似し, 整理すると
akn1  akn1
 ick akn
2t
akn1  akn1  2ick t akn
となる.
• 実空間での解の構造を求めるには, 離散フーリエ逆変換すれ
ばよい. ( (2) 式に代入すればよい.)
スペクトル法 (3)
 (x j , tn1), (x j , tn )
実空間データ
1 N n ikxj
ak (t )   j e
N j 1
離散フーリエ変換
akn1  akn1  2ick t akn
時間積分
離散フーリエ逆変換
 ( x j , tn1)
 
n
j
M
jkx
a
(
t
)
e
k
k  M
実空間データ
スペクトル法 (4)
• スペクトル法の利点
– 微分値の見積もり精度が差分法に比べて非常に高い.
– 格子点配置に任意性がない.
• スペクトル法の欠点
– 滑らかな関数を展開関数に用いると, 物理量の急激な変化を表
すことができない.
– 境界条件によっては適した展開関数がない.
スペクトル法 (5)
• 線型方程式
– 波数空間においても方程式がとても簡単.
– 計算量が少ない.
• 非線型方程式
– 波数空間での計算が”若干のややこしい”.
– (何も考えずに)計算すると計算量が多い.
• 非線型方程式を考える際には工夫が必要.
スペクトル法 (6)
• 以下の非線型方程式を考える.


 
t
x
• 解を以下のように展開することにして上の式に代入すると,
 
n
j
M
M
ikx
a
(
t
)
e
k
k  M

M
M
d
ikx
ilx
im x
a
(
t
)
e

a
(
t
)
e

ima
(
t
)
e



k
l
m
dt k  M
l  M
m  M
M
 i 
 ma (t)a
M
l   M m  M
より,
l

i ( l  m) x
(
t
)
e
m

スペクトル法 (7)
M
dak (t )

imal (t )am (t )

dt
l  m k ; l , m  M
となる.
スペクトル法 (8)
• 非線形項の計算量に関する考察.
M
dak (t )

imal (t )am (t )

dt
l  m k ; l , m  M
右辺の計算量は O(M2).
• 計算量が多すぎるので, 減らしたい.
スペクトル法 (9)
• 変換法
– 非線形項(の掛け算)を実空間で計算（微分は波
数空間で評価）.
dak    
  

dt   x k
   
akn1  akn1  2t  


x
k
 
n
 
n
j
M
a
m  M
m
(t )eim x
n
n
   
1 N  n     ikx
    j 
 e
 
 x  j 
  x k N j 1 

M
  

  imal (t )eim x
 x  j mM
n
n
j
スペクトル法 (10)
 (x j , tn1), (x j , tn ) 実空間データ
1 2M 1 n ikxj
ak (t ) 
 e
2M 1 j 1 j
  

  imamn eim x
 x m
n
波数空間で微分の評価
離散フーリエ変換
M
  

  imam (t )eim x
 x  j mM
n
離散フーリエ逆変換
時間積分
実空間での非線型項の評価（掛け算）
 nj 
離散フーリエ逆変換
M
 a (t)e
k  M
jkx
k
n
n
M 
   
1
 n     ikx






 
  j  x  j e
  x k 2M  1 k  M 

離散フーリエ変換
   
 t  

  x k
n
 ( x j , tn1)
  
 j 

 x  j
実空間データ
時間積分
n1
k
a
n1
k
a
スペクトル法 (11)
• 変換法を用いた場合の計算量
– 実質的には離散フーリエ（逆）変換の計算量で決
まる
• 離散フーリエ変換の計算量は 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)
 (x j , tn1), (x j , tn ) 実空間データ
1 2M 1 n ikxj
ak (t ) 
 e
2M 1 j 1 j
  

  imamn eim x
 x m
n
波数空間で微分の評価
離散フーリエ変換
M
  

  imam (t )eim x
 x  j mM
n
離散フーリエ逆変換
時間積分
実空間での非線型項の評価（掛け算）
 nj 
離散フーリエ逆変換
  
 j 

 x  j
M
 a (t)e
k  M
jkx
k
n
n
M 
   
1
 n     ikx






 
  j  x  j e
  x k 2M  1 k  M 

離散フーリエ逆変換
g_psi * g_Dx_e( e_psi )
   
 t  

  x k
n
 ( x j , tn1)
実空間データ
時間積分
n1
k
a
n1
k
a
スペクトル法 (14)
• 変換法を用いる際にはエイリアシングに注意が必要.
– 格子点数 N の時, 最初に考える切断波数の選択は
M=N/2
• 実空間でのデータ数（自由度）と波数空間でのデータ数が等しい.
• 波数 M までの成分を含む 2 項の積は, 波数 2M までの成分を含む. 実際
には折り返されて“偽の”成分が生じる（エイリアシング）.
• 畳み込み積分の結果は, 波数 M までしか含まない.
– そこで, エイリアシングが起こらない程度まで切断波数
を小さくする.
• M = (N-1)/3 とすると, 2 項までの積を含む系ではエイリアシングは起こら
ない.
k1+k2
k’
0
M
N/2
波数
2M
参考文献
• 石岡 (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.