スペクトル法のソースを表示


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

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

大域的アプローチの性質から、スペクトル法は解が[[滑らかな関数]]である場合に誤差が指数関数に従い収束するという特性（「指数収束」）を持ち、有限要素法よりも遥かに高速に収束することが知られている。 ただし、指数収束は解が[[滑らかな関数|滑らか]]'''でない'''場合には保証されないため、たとえば単連結な三次元定義域における{{ill|衝撃捕捉|en|shock-capturing methods}}<ref name="CHQZ">[https://books.google.com/books?id=7COgEw5_EBQC pp 235, Spectral Methods]: evolution to complex geometries and applications to fluid dynamics, By Canuto, Hussaini, Quarteroni and Zang, Springer, 2007.</ref>といった課題に対しては一般に成立しない。これはインパルス波の微分不可能性に起因する。なお、有限要素法の分野でも、要素の次数がグリッド幅''h''と反比例して大きくなるような手法を「{{ill|スペクトル要素法|en|Spectral element method}}」と呼ぶことがあるが、これはスペクトル法とは厳密には異なる手法である。（「#スペクトル要素法との関係」で後述） 

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

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

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

これらのスペクトル法は、通常、[[選点法]]や{{ill|ガラーキン法|en|Galerkin method}}、および[[タウ法]]のいずれかを用いることで実装される。 

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

== スペクトル法の例 ==

=== 具体例（線形の場合） ===
ここでは、基本的な多変量[[微分積分学|計算]]と[[フーリエ級数|フーリエ級数の]]理解を前提としている。 もし<math>g(x,y)</math>が2つの実変数を取る既知の複素関数であり、''g''が''x'', ''y''に関して周期的であるとき（つまり、 <math>g(x,y)=g(x+2\pi,y)=g(x,y+2\pi)</math>である場合 ）、以下を満たす関数 ''f''(''x'', ''y'') を見つけることを考える。 

: <math>\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</math>

ただし、左辺は''x'', ''y''における''f''の2次偏微分係数をそれぞれ示している。これは[[ポアソン方程式]]であり、物理的には熱伝導の問題、またはポテンシャル理論の問題として解釈できる。 

フーリエ級数で''f''と''g''を書くと、 

: <math>f=:\sum a_{j,k}e^{i(jx+ky)}</math>
: <math>g=:\sum b_{j,k}e^{i(jx+ky)}</math>

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

: <math>\sum -a_{j,k}(j^2+k^2)e^{i(jx+ky)}=\sum b_{j,k}e^{i(jx+ky)}</math>

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

: （*） <math>a_{j,k}=-\frac{b_{j,k}}{j^2+k^2}</math> 

これは、フーリエ係数''a<sub>j</sub>''<sub>,''k''</sub>の陽な表現である。 

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

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

==== アルゴリズム ====

# ''g''のフーリエ変換(''b<sub>j,k</sub>'') を計算 
# 式(*)を用いて''f''のフーリエ変換(''a<sub>j,k</sub>'')を計算 
# (''a<sub>j,k</sub>'')の逆フーリエ変換を実行して''f''を計算 

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

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

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

: <math>\partial_{t} u + u \partial_{x} u = \rho \partial_{xx} u + f \quad \forall x\in\left[0,2\pi\right), \forall t>0</math>

ここで、 ''&#x3C1;''は[[粘度]]係数である。 弱保存形では、これは次式のようになる。 

: <math>\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</math>

ただし、<math>\langle f, g \rangle := \int_{0}^{2\pi} f(x)
 \overline{g(x)}\,dx</math>は[[計量ベクトル空間|内積]]である。 [[部分積分]]と周期性により、 

: <math>\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.</math>

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

: <math>\mathcal{U}^N := \left\{ u : u(x,t)=\sum_{k=-N/2}^{N/2-1} \hat{u}_{k}(t) e^{i k x}\right\}</math>  

: <math>\mathcal{V}^N :=\operatorname{span}\left\{ e^{i k x} : k\in -N/2,\dots,N/2-1\right\}</math>

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

: <math>\langle \partial_{t} u, e^{i k x} \rangle = \left\langle \frac{1}{2} u^2 - \rho \partial_{x} u , \partial_x e^{i k x} \right\rangle + \left\langle f, e^{i k x} \right\rangle \quad \forall k\in \left\{ -N/2,\dots,N/2-1 \right\}, \forall t>0.</math>

ここで、[[直交]]関係<math>\langle e^{i l x}, e^{i k x} \rangle = 2 \pi \delta_{lk}</math>を利用している。ただし<math>\delta_{lk}</math>は[[クロネッカーのデルタ|クロネッカーデルタ]]である。上記の3つの項を<math>k</math>について整理すると次のようになる。 

: <math>
\begin{align}
\left\langle \partial_{t} u, e^{i k x}\right\rangle &= \left\langle   \partial_{t} \sum_{l} \hat{u}_{l} e^{i l x}  ,   e^{i k x} \right\rangle = \left\langle  \sum_{l} \partial_{t} \hat{u}_{l} e^{i l x}  ,   e^{i k x} \right\rangle = 2 \pi \partial_t \hat{u}_k,
\\
\left\langle f, e^{i k x} \right\rangle &= \left\langle  \sum_{l} \hat{f}_{l} e^{i l x}  ,   e^{i k x}\right\rangle= 2 \pi \hat{f}_k, \text{ and}
\\
\left\langle
 \frac{1}{2} u^2 - \rho \partial_{x} u
 ,
 \partial_x e^{i k x}
\right\rangle
&=
\left\langle
  \frac{1}{2}
  \left(\sum_{p} \hat{u}_p e^{i p x}\right)
  \left(\sum_{q} \hat{u}_q e^{i q x}\right)
  - \rho \partial_x \sum_{l} \hat{u}_l e^{i l x}
  ,
  \partial_x e^{i k x}
\right\rangle
\\
&=
\left\langle
  \frac{1}{2}
  \sum_{p} \sum_{q} \hat{u}_p \hat{u}_q e^{i \left(p+q\right) x}
  ,
  i k e^{i k x}
\right\rangle
-
\left\langle
  \rho i \sum_{l} l \hat{u}_l e^{i l x}
  ,
  i k e^{i k x}
\right\rangle
\\
&=
-\frac{i k}{2}
\left\langle
  \sum_{p} \sum_{q} \hat{u}_p \hat{u}_q e^{i \left(p+q\right) x}
  ,
  e^{i k x}
\right\rangle
- \rho k
\left\langle
  \sum_{l} l \hat{u}_l e^{i l x}
  ,
  e^{i k x}
\right\rangle
\\
&=
- i \pi k \sum_{p+q=k} \hat{u}_p \hat{u}_q - 2\pi\rho{}k^2\hat{u}_k.
\end{align}
</math>

これらの3つの項を<math>k</math>についてまとめることで次式を得る。 

: <math>
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.
</math>

両辺を<math>2\pi</math>で除することで、最終的に次式を得る。 

: <math>
\partial_t \hat{u}_k
=
- \frac{i k}{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.
</math>

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

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

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

== 関連記事 ==

* [[有限要素法]] 
* {{ill|ガウス格子|en|Gaussian grid}} 
* {{ill|疑似スペクトル法|en|Pseudo-spectral method}}
* {{ill|スペクトル要素法|en|Spectral element method}}
* {{ill|ガラーキン法|en|Galerkin method}}
* [[選点法]] 

== 参考文献 ==
{{Reflist}}

* Bengt Fornberg (1996): ''A Practical Guide to Pseudospectral Methods'', Cambridge University Press, Cambridge, UK
* [https://www.researchgate.net/publication/2405798_Chebyshev_and_Fourier_Spectral_Methods John P. Boyd (2000): ''Chebyshev and Fourier Spectral Methods'', 2nd Ed.]
* Canuto C., [[:en:M._Yousuff_Hussaini|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: [http://epubs.siam.org/sam-bin/dbq/article/35198 ''A Spectral Element Method for the Navier–Stokes Equations with Improved Accuracy'']
* [http://cdm.unimo.it/home/matematica/funaro.daniele/bube.htm ''Polynomial Approximation of Differential Equations''], by Daniele Funaro, (Lecture Notes in Physics, Volume 8), Spinger-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
* {{Cite book|last1=Press|first1=WH|last2=Teukolsky|first2=SA|last3=Vetterling|first3=WT|last4=Flannery|first4=BP|year=2007|title=Numerical Recipes: The Art of Scientific Computing|edition=3rd|publisher=Cambridge University Press|location=New York|isbn=978-0-521-88068-8|chapter=Section 20.7. ''Spectral Methods'' |chapter-url=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), {{ISBN2|3-540-71040-X}}
* Lloyd N. Trefethen (2000): ''Spectral Methods in MATLAB'', SIAM, Philadelphia, PA
{{デフォルトソート:すへくとるほう}}
[[Category:数値微分方程式]]
[[Category:数値解析]]
[[Category:数学に関する記事]]