スペクトル法 は、主に高速フーリエ変換 を用いた微分方程式 の数値解法の総称であり、応用数学 や科学計算 で使用されている。 微分方程式の解をある「基底関数 」の和によって近似し、方程式を充足する和の係数を求める。基底関数の選び方としては、例えば正弦波 を用いる方法があり、この場合の解の近似表現はフーリエ級数 になる。
スペクトル法は有限要素法 と密接に関連しており、基本的にはどちらも同じアイデアに基づいている。これらの主な違いは、近似に用いる基底関数の定義域にある。スペクトル法は近似対象とする関数の定義域全体に渡って非零になるような基底関数を使用するため全体をカバーできるのに対し、有限要素法はある点の近傍など限られた範囲にのみ基底関数を用い、残りはゼロであると仮定する。こうした理由から、スペクトル法と有限要素法はそれぞれ、大域的アプローチ、局所的アプローチと呼ばれ区別される。
大域的アプローチの性質から、スペクトル法は解が滑らかな関数 である場合に誤差が指数関数に従い収束するという特性(「指数収束」)を持ち、有限要素法よりも遥かに高速に収束することが知られている。 ただし、指数収束は解が滑らか でない 場合には保証されないため、たとえば単連結な三次元定義域における衝撃捕捉 (英語版 ) [1] といった課題に対しては一般に成立しない。これはインパルス波の微分不可能性に起因する。なお、有限要素法の分野でも、要素の次数がグリッド幅h と反比例して大きくなるような手法を「スペクトル要素法 (英語版 ) 」と呼ぶことがあるが、これはスペクトル法とは厳密には異なる手法である。(「#スペクトル要素法との関係」で後述)
スペクトル法を使用すると、 常微分方程式 (ODE)や偏微分方程式 (PDE)などの微分方程式を含む固有値 問題を解くことができる。 時間依存のPDEにスペクトル法を適用した場合、解は時間依存の係数を持つ基底関数の合計として記述される。これをPDEに代入すると、ODEの任意の数値法 を使用して解くことができる係数のODEシステムが生成される。 ODEの固有値問題も同様に行列固有値問題に変換される[要出典 ] 。
スペクトル法は、1969年 より数学者のスティーブン・オルザグ によって出版された複数編に渡る論文により確立されたものであるが、一連の論文は今日で多く実装されている周期幾何問題を対象にしたフーリエ級数を用いたもの以外にも、以下のような手法を含んでいる。
有限幾何・非有界幾何のための多項式スペクトル法 高次非線形問題のための擬球スペクトル法 定常問題の高速解法のためのスペクトル反復法 これらのスペクトル法は、通常、選点法 やガラーキン法 (英語版 ) 、およびタウ法 のいずれかを用いることで実装される。
スペクトル法は有限要素法よりも計算コストが低くなるが、複素幾何や不連続係数の問題では精度が低下する。 この誤差の増加は、 ギブス現象 によるものである。
スペクトル法の例 [ 編集 ] 具体例(線形の場合) [ 編集 ] ここでは、基本的な多変量計算 とフーリエ級数の 理解を前提としている。 もし g ( x , y ) {\displaystyle g(x,y)} が2つの実変数を取る既知の複素関数であり、g がx , y に関して周期的であるとき(つまり、 g ( x , y ) = g ( x + 2 π , y ) = g ( x , y + 2 π ) {\displaystyle g(x,y)=g(x+2\pi ,y)=g(x,y+2\pi )} である場合 )、以下を満たす関数 f (x , y ) を見つけることを考える。
( ∂ 2 ∂ x 2 + ∂ 2 ∂ y 2 ) f ( x , y ) = g ( x , y ) for all x , y {\displaystyle \left({\frac {\partial ^{2}}{\partial x^{2}}}+{\frac {\partial ^{2}}{\partial y^{2}}}\right)f(x,y)=g(x,y)\quad {\text{for all }}x,y} ただし、左辺はx , y におけるf の2次偏微分係数をそれぞれ示している。これはポアソン方程式 であり、物理的には熱伝導の問題、またはポテンシャル理論の問題として解釈できる。
フーリエ級数でf とg を書くと、
f =: ∑ a j , k e i ( j x + k y ) {\displaystyle f=:\sum a_{j,k}e^{i(jx+ky)}} g =: ∑ b j , k e i ( j x + k y ) {\displaystyle g=:\sum b_{j,k}e^{i(jx+ky)}} であり、これを微分方程式に代入すると、次の方程式が得られる。
∑ − a j , k ( j 2 + k 2 ) e i ( j x + k y ) = ∑ b j , k e i ( j x + k y ) {\displaystyle \sum -a_{j,k}(j^{2}+k^{2})e^{i(jx+ky)}=\sum b_{j,k}e^{i(jx+ky)}} ここで偏微分を無限和と交換している。これは、たとえばf に連続的な2次導関数があると仮定した場合に正当である。フーリエ展開の一意性定理により、フーリエ係数を項ごとに等しくする必要がある。
(*) a j , k = − b j , k j 2 + k 2 {\displaystyle a_{j,k}=-{\frac {b_{j,k}}{j^{2}+k^{2}}}} これは、フーリエ係数aj ,k の陽な表現である。
周期的境界条件から、ポアソン方程式 はb 0,0 = 0の場合に限り解を持つ。したがって、我々は自由に解の平均値a 0,0 を選択することができる。これは、積分定数の選択に対応する。
ここからアルゴリズムを構成するため、有限数の周波数のみを解く。 これにより、 h n {\displaystyle h^{n}} に比例する誤差が発生する。ただし h := 1 / n {\displaystyle h:=1/n} であり、 n {\displaystyle n} は処理対象の最大周波数である。
アルゴリズム [ 編集 ] g のフーリエ変換(bj,k ) を計算 式(*)を用いてf のフーリエ変換(aj,k )を計算 (aj,k )の逆フーリエ変換を実行してf を計算 ここでは幅n の周波数の有限窓のみに関心があるため、このアルゴリズムは高速フーリエ変換 を使用して実行できる。したがって、アルゴリズムはグローバルにO (n log n ) 時間で実行できる。
非線形の場合 [ 編集 ] スペクトルアプローチを使用し、強制的な非定常非線形バーガース方程式 を解く。
与えられた u ( x , 0 ) {\displaystyle u(x,0)} 周期領域で x ∈ [ 0 , 2 π ) {\displaystyle x\in \left[0,2\pi \right)} 、次式を満たす u ∈ U {\displaystyle u\in {\mathcal {U}}} を見つけることを考える。
∂ t u + u ∂ x u = ρ ∂ x x u + f ∀ x ∈ [ 0 , 2 π ) , ∀ t > 0 {\displaystyle \partial _{t}u+u\partial _{x}u=\rho \partial _{xx}u+f\quad \forall x\in \left[0,2\pi \right),\forall t>0} ここで、 ρ は粘度 係数である。 弱保存形では、これは次式のようになる。
⟨ ∂ t u , v ⟩ = ⟨ ∂ x ( − 1 2 u 2 + ρ ∂ x u ) , v ⟩ + ⟨ f , v ⟩ ∀ v ∈ V , ∀ t > 0 {\displaystyle \left\langle \partial _{t}u,v\right\rangle =\left\langle \partial _{x}\left(-{\frac {1}{2}}u^{2}+\rho \partial _{x}u\right),v\right\rangle +\left\langle f,v\right\rangle \quad \forall v\in {\mathcal {V}},\forall t>0} ただし、 ⟨ f , g ⟩ := ∫ 0 2 π f ( x ) g ( x ) ¯ d x {\displaystyle \langle f,g\rangle :=\int _{0}^{2\pi }f(x){\overline {g(x)}}\,dx} は内積 である。 部分積分 と周期性により、
⟨ ∂ t u , v ⟩ = ⟨ 1 2 u 2 − ρ ∂ x u , ∂ x v ⟩ + ⟨ f , v ⟩ ∀ v ∈ V , ∀ t > 0. {\displaystyle \langle \partial _{t}u,v\rangle =\left\langle {\frac {1}{2}}u^{2}-\rho \partial _{x}u,\partial _{x}v\right\rangle +\left\langle f,v\right\rangle \quad \forall v\in {\mathcal {V}},\forall t>0.} フーリエ- ガラーキン法を適用するには、以下の両方を選択する。
U N := { u : u ( x , t ) = ∑ k = − N / 2 N / 2 − 1 u ^ k ( t ) e i k x } {\displaystyle {\mathcal {U}}^{N}:=\left\{u:u(x,t)=\sum _{k=-N/2}^{N/2-1}{\hat {u}}_{k}(t)e^{ikx}\right\}} V N := span { e i k x : k ∈ − N / 2 , … , N / 2 − 1 } {\displaystyle {\mathcal {V}}^{N}:=\operatorname {span} \left\{e^{ikx}:k\in -N/2,\dots ,N/2-1\right\}} ただし、 u ^ k ( t ) := 1 2 π ⟨ u ( x , t ) , e i k x ⟩ {\displaystyle {\hat {u}}_{k}(t):={\frac {1}{2\pi }}\langle u(x,t),e^{ikx}\rangle } 。 これにより、 u ∈ U N {\displaystyle u\in {\mathcal {U}}^{N}} の探索は以下の問題に帰着される。
⟨ ∂ t u , e i k x ⟩ = ⟨ 1 2 u 2 − ρ ∂ x u , ∂ x e i k x ⟩ + ⟨ f , e i k x ⟩ ∀ k ∈ { − N / 2 , … , N / 2 − 1 } , ∀ t > 0. {\displaystyle \langle \partial _{t}u,e^{ikx}\rangle =\left\langle {\frac {1}{2}}u^{2}-\rho \partial _{x}u,\partial _{x}e^{ikx}\right\rangle +\left\langle f,e^{ikx}\right\rangle \quad \forall k\in \left\{-N/2,\dots ,N/2-1\right\},\forall t>0.} ここで、直交 関係 ⟨ e i l x , e i k x ⟩ = 2 π δ l k {\displaystyle \langle e^{ilx},e^{ikx}\rangle =2\pi \delta _{lk}} を利用している。ただし δ l k {\displaystyle \delta _{lk}} はクロネッカーデルタ である。上記の3つの項を k {\displaystyle k} について整理すると次のようになる。
⟨ ∂ t u , e i k x ⟩ = ⟨ ∂ t ∑ l u ^ l e i l x , e i k x ⟩ = ⟨ ∑ l ∂ t u ^ l e i l x , e i k x ⟩ = 2 π ∂ t u ^ k , ⟨ f , e i k x ⟩ = ⟨ ∑ l f ^ l e i l x , e i k x ⟩ = 2 π f ^ k , and ⟨ 1 2 u 2 − ρ ∂ x u , ∂ x e i k x ⟩ = ⟨ 1 2 ( ∑ p u ^ p e i p x ) ( ∑ q u ^ q e i q x ) − ρ ∂ x ∑ l u ^ l e i l x , ∂ x e i k x ⟩ = ⟨ 1 2 ∑ p ∑ q u ^ p u ^ q e i ( p + q ) x , i k e i k x ⟩ − ⟨ ρ i ∑ l l u ^ l e i l x , i k e i k x ⟩ = − i k 2 ⟨ ∑ p ∑ q u ^ p u ^ q e i ( p + q ) x , e i k x ⟩ − ρ k ⟨ ∑ l l u ^ l e i l x , e i k x ⟩ = − i π k ∑ p + q = k u ^ p u ^ q − 2 π ρ k 2 u ^ k . {\displaystyle {\begin{aligned}\left\langle \partial _{t}u,e^{ikx}\right\rangle &=\left\langle \partial _{t}\sum _{l}{\hat {u}}_{l}e^{ilx},e^{ikx}\right\rangle =\left\langle \sum _{l}\partial _{t}{\hat {u}}_{l}e^{ilx},e^{ikx}\right\rangle =2\pi \partial _{t}{\hat {u}}_{k},\\\left\langle f,e^{ikx}\right\rangle &=\left\langle \sum _{l}{\hat {f}}_{l}e^{ilx},e^{ikx}\right\rangle =2\pi {\hat {f}}_{k},{\text{ and}}\\\left\langle {\frac {1}{2}}u^{2}-\rho \partial _{x}u,\partial _{x}e^{ikx}\right\rangle &=\left\langle {\frac {1}{2}}\left(\sum _{p}{\hat {u}}_{p}e^{ipx}\right)\left(\sum _{q}{\hat {u}}_{q}e^{iqx}\right)-\rho \partial _{x}\sum _{l}{\hat {u}}_{l}e^{ilx},\partial _{x}e^{ikx}\right\rangle \\&=\left\langle {\frac {1}{2}}\sum _{p}\sum _{q}{\hat {u}}_{p}{\hat {u}}_{q}e^{i\left(p+q\right)x},ike^{ikx}\right\rangle -\left\langle \rho i\sum _{l}l{\hat {u}}_{l}e^{ilx},ike^{ikx}\right\rangle \\&=-{\frac {ik}{2}}\left\langle \sum _{p}\sum _{q}{\hat {u}}_{p}{\hat {u}}_{q}e^{i\left(p+q\right)x},e^{ikx}\right\rangle -\rho k\left\langle \sum _{l}l{\hat {u}}_{l}e^{ilx},e^{ikx}\right\rangle \\&=-i\pi k\sum _{p+q=k}{\hat {u}}_{p}{\hat {u}}_{q}-2\pi \rho {}k^{2}{\hat {u}}_{k}.\end{aligned}}} これらの3つの項を k {\displaystyle k} についてまとめることで次式を得る。
2 π ∂ t u ^ k = − i π k ∑ p + q = k u ^ p u ^ q − 2 π ρ k 2 u ^ k + 2 π f ^ k k ∈ { − N / 2 , … , N / 2 − 1 } , ∀ t > 0. {\displaystyle 2\pi \partial _{t}{\hat {u}}_{k}=-i\pi k\sum _{p+q=k}{\hat {u}}_{p}{\hat {u}}_{q}-2\pi \rho {}k^{2}{\hat {u}}_{k}+2\pi {\hat {f}}_{k}\quad k\in \left\{-N/2,\dots ,N/2-1\right\},\forall t>0.} 両辺を 2 π {\displaystyle 2\pi } で除することで、最終的に次式を得る。
∂ t u ^ k = − i k 2 ∑ p + q = k u ^ p u ^ q − ρ k 2 u ^ k + f ^ k k ∈ { − N / 2 , … , N / 2 − 1 } , ∀ t > 0. {\displaystyle \partial _{t}{\hat {u}}_{k}=-{\frac {ik}{2}}\sum _{p+q=k}{\hat {u}}_{p}{\hat {u}}_{q}-\rho {}k^{2}{\hat {u}}_{k}+{\hat {f}}_{k}\quad k\in \left\{-N/2,\dots ,N/2-1\right\},\forall t>0.} フーリエ変換後の初期条件 u ^ k ( 0 ) {\displaystyle {\hat {u}}_{k}(0)} と外力 f ^ k ( t ) {\displaystyle {\hat {f}}_{k}(t)} を入力として与えることで、この常微分方程式の結合系の時間発展は、ルンゲ=クッタ法 などを使った数値積分によって解くことができる。 第一項(非線形項)は畳み込み 演算であるため、これも効率的に評価するための変換がいくつか存在する。 BoydおよびCanutoらの参考文献を参照してください。詳細については。
スペクトル要素法との関係 [ 編集 ] g {\displaystyle g} が無限回微分可能な関数であるとき、高速フーリエ変換を使用する数値アルゴリズムは、グリッドサイズhのどの多項式よりも速く収束することを示すことができる。 つまり、n が正であるとき、任意の十分小さな値 h {\displaystyle h} に対し、誤差が C n h n {\displaystyle C_{n}h^{n}} 以下になるような C n < ∞ {\displaystyle C_{n}<\infty } が存在する。n >0に対し、 C n {\displaystyle C_{n}} が適切に選ばれることでこの誤差条件が満たされる手法は n {\displaystyle n} 次スペクトル法と呼ばれる。
スペクトル要素法 もまた非常に高次の有限要素法 であるため、収束特性には類似点がある。 ただし、スペクトル法は特定の境界値問題の固有分解を用いるためそれだけ適用範囲が狭くなるが、有限要素法はこうした固有分解に依存しないため、任意の楕円境界値問題 に対して適用することができる。
関連記事 [ 編集 ] 参考文献 [ 編集 ] ^ pp 235, Spectral Methods : evolution to complex geometries and applications to fluid dynamics, By Canuto, Hussaini, Quarteroni and Zang, Springer, 2007. Bengt Fornberg (1996) A Practical Guide to Pseudospectral Methods. Cambridge University Press, Cambridge, UK Chebyshev and Fourier Spectral Methods by John P. Boyd. Canuto C., Hussaini M. Y. , Quarteroni A., and Zang T.A. (2006) Spectral Methods. Fundamentals in Single Domains. Springer-Verlag, Berlin Heidelberg Javier de Frutos, Julia Novo: A Spectral Element Method for the Navier–Stokes Equations with Improved Accuracy Polynomial Approximation of Differential Equations , by Daniele Funaro, Lecture Notes in Physics, Volume 8, Springer-Verlag, Heidelberg 1992 D. Gottlieb and S. Orzag (1977) "Numerical Analysis of Spectral Methods : Theory and Applications", SIAM, Philadelphia, PA J. Hesthaven, S. Gottlieb and D. Gottlieb (2007) "Spectral methods for time-dependent problems", Cambridge UP, Cambridge, UK Steven A. Orszag (1969) Numerical Methods for the Simulation of Turbulence , Phys. Fluids Supp. II, 12, 250–257 Press, WH; Teukolsky, SA; Vetterling, WT; Flannery, BP (2007). “Section 20.7. Spectral Methods” . Numerical Recipes: The Art of Scientific Computing (3rd ed.). New York: Cambridge University Press. ISBN 978-0-521-88068-8 . http://apps.nrbook.com/empanel/index.html#pg=1083 Jie Shen, Tao Tang and Li-Lian Wang (2011) "Spectral Methods: Algorithms, Analysis and Applications" (Springer Series in Computational Mathematics, V. 41, Springer), ISBN 354071040X Lloyd N. Trefethen (2000) Spectral Methods in MATLAB. SIAM, Philadelphia, PA