スペクトル法

スペクトル法は、主に高速フーリエ変換を用いた微分方程式の数値解法の総称であり、応用数学や科学計算で使用されている。 微分方程式の解をある「基底関数」の和によって近似し、方程式を充足する和の係数を求める。基底関数の選び方としては、例えば正弦波を用いる方法があり、この場合の解の近似表現はフーリエ級数になる。

スペクトル法は有限要素法と密接に関連しており、基本的にはどちらも同じアイデアに基づいている。これらの主な違いは、近似に用いる基底関数の定義域にある。スペクトル法は近似対象とする関数の定義域全体に渡って非零になるような基底関数を使用するため全体をカバーできるのに対し、有限要素法はある点の近傍など限られた範囲にのみ基底関数を用い、残りはゼロであると仮定する。こうした理由から、スペクトル法と有限要素法はそれぞれ、大域的アプローチ、局所的アプローチと呼ばれ区別される。

大域的アプローチの性質から、スペクトル法は解が滑らかな関数である場合に誤差が指数関数に従い収束するという特性（「指数収束」）を持ち、有限要素法よりも遥かに高速に収束することが知られている。 ただし、指数収束は解が滑らかでない場合には保証されないため、たとえば単連結な三次元定義域における衝撃捕捉（英語版）^[1]といった課題に対しては一般に成立しない。これはインパルス波の微分不可能性に起因する。なお、有限要素法の分野でも、要素の次数がグリッド幅hと反比例して大きくなるような手法を「スペクトル要素法（英語版）」と呼ぶことがあるが、これはスペクトル法とは厳密には異なる手法である。（「#スペクトル要素法との関係」で後述）

スペクトル法を使用すると、 常微分方程式 （ODE）や偏微分方程式 （PDE）などの微分方程式を含む固有値問題を解くことができる。 時間依存のPDEにスペクトル法を適用した場合、解は時間依存の係数を持つ基底関数の合計として記述される。これをPDEに代入すると、ODEの任意の数値法を使用して解くことができる係数のODEシステムが生成される。 ODEの固有値問題も同様に行列固有値問題に変換される^[要出典] 。

スペクトル法は、1969年より数学者のスティーブン・オルザグによって出版された複数編に渡る論文により確立されたものであるが、一連の論文は今日で多く実装されている周期幾何問題を対象にしたフーリエ級数を用いたもの以外にも、以下のような手法を含んでいる。

• 有限幾何・非有界幾何のための多項式スペクトル法
• 高次非線形問題のための擬球スペクトル法
• 定常問題の高速解法のためのスペクトル反復法

これらのスペクトル法は、通常、選点法やガラーキン法（英語版）、およびタウ法のいずれかを用いることで実装される。

スペクトル法は有限要素法よりも計算コストが低くなるが、複素幾何や不連続係数の問題では精度が低下する。 この誤差の増加は、 ギブス現象によるものである。

ここでは、基本的な多変量計算とフーリエ級数の理解を前提としている。 もし${\displaystyle g(x,y)}$が2つの実変数を取る既知の複素関数であり、gがx, yに関して周期的であるとき（つまり、 ${\displaystyle g(x,y)=g(x+2\pi ,y)=g(x,y+2\pi )}$である場合 ）、以下を満たす関数 f(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を書くと、

${\displaystyle f=:\sum a_{j,k}e^{i(jx+ky)}}$

${\displaystyle g=:\sum b_{j,k}e^{i(jx+ky)}}$

であり、これを微分方程式に代入すると、次の方程式が得られる。

${\displaystyle \sum -a_{j,k}(j^{2}+k^{2})e^{i(jx+ky)}=\sum b_{j,k}e^{i(jx+ky)}}$

ここで偏微分を無限和と交換している。これは、たとえばfに連続的な2次導関数があると仮定した場合に正当である。フーリエ展開の一意性定理により、フーリエ係数を項ごとに等しくする必要がある。

（*） ${\displaystyle a_{j,k}=-{\frac {b_{j,k}}{j^{2}+k^{2}}}}$

これは、フーリエ係数a_j,kの陽な表現である。

周期的境界条件から、ポアソン方程式はb_0,0 = 0の場合に限り解を持つ。したがって、我々は自由に解の平均値a_0,0を選択することができる。これは、積分定数の選択に対応する。

ここからアルゴリズムを構成するため、有限数の周波数のみを解く。 これにより、${\displaystyle h^{n}}$に比例する誤差が発生する。ただし${\displaystyle h:=1/n}$であり、${\displaystyle n}$は処理対象の最大周波数である。

1. gのフーリエ変換(b_j,k) を計算
2. 式(*)を用いてfのフーリエ変換(a_j,k)を計算
3. (a_j,k)の逆フーリエ変換を実行してfを計算

ここでは幅nの周波数の有限窓のみに関心があるため、このアルゴリズムは高速フーリエ変換を使用して実行できる。したがって、アルゴリズムはグローバルにO(n log n)時間で実行できる。

スペクトルアプローチを使用し、強制的な非定常非線形バーガース方程式を解く。

与えられた${\displaystyle u(x,0)}$周期領域で${\displaystyle x\in \left[0,2\pi \right)}$ 、次式を満たす${\displaystyle u\in {\mathcal {U}}}$を見つけることを考える。

${\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}$

ここで、 ρは粘度係数である。 弱保存形では、これは次式のようになる。

${\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}$

ただし、${\displaystyle \langle f,g\rangle :=\int _{0}^{2\pi }f(x){\overline {g(x)}}\,dx}$は内積である。 部分積分と周期性により、

${\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.}$

フーリエ- ガラーキン法を適用するには、以下の両方を選択する。

${\displaystyle {\mathcal {U}}^{N}:=\left\{u:u(x,t)=\sum _{k=-N/2}^{N/2-1}{\hat {u}}_{k}(t)e^{ikx}\right\}}$

${\displaystyle {\mathcal {V}}^{N}:=\operatorname {span} \left\{e^{ikx}:k\in -N/2,\dots ,N/2-1\right\}}$

ただし、${\displaystyle {\hat {u}}_{k}(t):={\frac {1}{2\pi }}\langle u(x,t),e^{ikx}\rangle }$ 。 これにより、${\displaystyle u\in {\mathcal {U}}^{N}}$の探索は以下の問題に帰着される。

${\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.}$

ここで、直交関係${\displaystyle \langle e^{ilx},e^{ikx}\rangle =2\pi \delta _{lk}}$を利用している。ただし${\displaystyle \delta _{lk}}$はクロネッカーデルタである。上記の3つの項を${\displaystyle 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つの項を${\displaystyle k}$についてまとめることで次式を得る。

${\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.}$

両辺を${\displaystyle 2\pi }$で除することで、最終的に次式を得る。

${\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.}$

フーリエ変換後の初期条件${\displaystyle {\hat {u}}_{k}(0)}$と外力${\displaystyle {\hat {f}}_{k}(t)}$ を入力として与えることで、この常微分方程式の結合系の時間発展は、ルンゲ＝クッタ法などを使った数値積分によって解くことができる。 第一項（非線形項）は畳み込み演算であるため、これも効率的に評価するための変換がいくつか存在する。 BoydおよびCanutoらの参考文献を参照してください。詳細については。

${\displaystyle g}$が無限回微分可能な関数であるとき、高速フーリエ変換を使用する数値アルゴリズムは、グリッドサイズhのどの多項式よりも速く収束することを示すことができる。 つまり、nが正であるとき、任意の十分小さな値${\displaystyle h}$に対し、誤差が${\displaystyle C_{n}h^{n}}$以下になるような ${\displaystyle C_{n}<\infty }$が存在する。n>0に対し、${\displaystyle C_{n}}$が適切に選ばれることでこの誤差条件が満たされる手法は${\displaystyle n}$次スペクトル法と呼ばれる。

スペクトル要素法もまた非常に高次の有限要素法であるため、収束特性には類似点がある。 ただし、スペクトル法は特定の境界値問題の固有分解を用いるためそれだけ適用範囲が狭くなるが、有限要素法はこうした固有分解に依存しないため、任意の楕円境界値問題に対して適用することができる。

• 有限要素法
• ガウス格子（英語版）
• 疑似スペクトル法（英語版）
• スペクトル要素法（英語版）
• ガラーキン法（英語版）
• 選点法

• 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