選点法

数学 > 数値解析 > 偏微分方程式の数値解法 > 選点法

数学 > 数値解析 > 常微分方程式の数値解法 > 選点法

選点法（せんてんほう、英: Collocation method） とは、数値解析において常微分方程式、偏微分方程式と積分方程式に対して数値解を与える方法である。この方法のアイディアは、解候補（通常はある次数以下の多項式）からなる有限次元のベクトル空間と定義域から幾つかの点を先に選び、それらの点で与えられた方程式を満足する解を解候補の空間から選択することである。そのように選ばれた点は、選点（collocation points）と呼ぶ。

常微分方程式

区間 [t₀, t₀ + h] 上の常微分方程式の初期値問題を以下のように定める。

$y'=f(t,y),\quad y(t_{0})=y_{0}.$

解候補の空間を n 次以下の多項式からなるベクトル空間とし、点 0 ≤ c₁ < c₂ < ... < c_n ≤ 1 を選ぶ（時々、c_k も選点と呼ばれるが、方程式の定義域が [0,1] とは限らないので最初に述べた選点とは少し違う定義となる）。

それらの点に対応する選点法による近似（多項式）解 p(t) は、初期値条件 p(t₀) = y₀ と選点 t_k = t₀ + c_kh (1 ≤ k ≤ n) での微分方程式 p'(t_k) = f (t_k, p(t_k)) を満足する。そのように与えられた方程式の数は n+1 であり、p の未知係数の数と一致する。そのため、近似解は一意に定められる。したがって、最後の時刻 t₀ + h での近似解は p(t₀ + h) となる。

以上のような方法は実際にすべて陰的ルンゲ＝クッタ法であり、以下のブッチャー配列で与えられる。

${\begin{array}{c|ccccc}c_{1}&a_{11}&a_{12}&\cdots &a_{1n}\\c_{2}&a_{11}&a_{12}&\cdots &a_{2n}\\\vdots &\vdots &\vdots &\ddots &\vdots \\c_{n}&a_{11}&a_{12}&\cdots &a_{nn}\\\hline &b_{1}&b_{2}&\cdots &b_{n}\end{array}}$

この中で、c_k は選ばれた点に対応し、b_k と a_kj は以下の公式で与えられる。

$b_{k}=\int _{0}^{1}{\frac {q_{k}(\tau )}{q_{k}(c_{k})}}d\tau$

$a_{kj}=\int _{0}^{c_{k}}{\frac {q_{j}(\tau )}{q_{j}(c_{j})}}d\tau$

ここで、$q_{j}(\tau )=(\tau -c_{j})^{-1}\prod _{i=1}^{n}(\tau -c_{i})$ である。

しかし、陰的ルンゲ＝クッタ法がすべて選点法であるわけではない。

計算例：台形公式

例として、上記の常微分方程式に対し点 c₁ = 0 と c₂ = 1 を選ぶとしよう（よって n = 2）。近似解 p(t) は2次多項式であり、選点法による方程式

$p(t_{0})=y_{0},\,$

$p'(t_{0})=f(t_{0},p(t_{0})),\,$

$p'(t_{0}+h)=f(t_{0}+h,p(t_{0}+h))\,$

を満足する。計算を簡単にするために、p を以下の形で書く。

$p(t)=\alpha (t-t_{0})^{2}+\beta (t-t_{0})+\gamma \,$

上記の方程式系を用いて未知の係数を解くと

${\begin{aligned}\alpha &={\frac {1}{2h}}\left(f(t_{0}+h,p(t_{0}+h))-f(t_{0},p(t_{0}))\right),\\\beta &=f(t_{0},p(t_{0})),\\\gamma &=y_{0}.\end{aligned}}$

であることがわかる。したがって対応する選点法は次の公式で（陰的に）与えられる。

$y_{1}=p(t_{0}+h)=y_{0}+{\frac {1}{2}}h\left(f(t_{0}+h,y_{1})+f(t_{0},y_{0})\right),\,$

ここで、y₁ = p(t₀ + h) は t = t₀ + h での近似解である。

この方法は微分方程式における台形公式として知られている。確かに、方程式を以下のように書き換えって（数値積分における）台形公式で近似することから上記の公式を導出することも可能である。

$y(t)=y(t_{0})+\int _{t_{0}}^{t}f(\tau ,y(\tau ))\,{\textrm {d}}\tau ,\,$

ガウス・ルジャンドル法

上述の点 c_i は n 次のルジャンドル多項式の根として選べる。それらの点に対応する選点法は、n 段ガウス・ルジャンドル法（Gauss-Legendre method）として知られている。n 段ガウス・ルジャンドルは次数 2n を持ち、n 段ルンゲ＝クッタ法の中にでも一番精度の高い方法である。さらに、ガウス・ルジャンドル法はすべてA-安定であり、硬い方程式にも適用できる。しかし n が4以上の時、ルジャンドル多項式の根を効率良く計算することが困難であり、加えて対応する方法の係数も極めて複雑であるので、4段以上のガウス・ルジャンドル法はあまり使われない。

2段ガウス・ルジャンドル法は以下のブッチャー配列で与えられる。

${\begin{array}{c|cc}{\frac {1}{2}}-{\frac {\sqrt {3}}{6}}&{\frac {1}{4}}&{\frac {1}{4}}-{\frac {\sqrt {3}}{6}}\\{\frac {1}{2}}+{\frac {\sqrt {3}}{6}}&{\frac {1}{4}}+{\frac {\sqrt {3}}{6}}&{\frac {1}{4}}\\\hline &{\frac {1}{2}}&{\frac {1}{2}}\\\end{array}}$

そして3段ガウス・ルジャンドル法は以下のブッチャー配列で与えられる。

${\begin{array}{c|ccc}{\frac {1}{2}}-{\frac {\sqrt {15}}{10}}&{\frac {5}{36}}&{\frac {2}{9}}-{\frac {\sqrt {15}}{15}}&{\frac {5}{36}}-{\frac {\sqrt {15}}{30}}\\{\frac {1}{2}}&{\frac {5}{36}}+{\frac {\sqrt {15}}{24}}&{\frac {2}{9}}&{\frac {5}{36}}-{\frac {\sqrt {15}}{24}}\\{\frac {1}{2}}+{\frac {\sqrt {15}}{10}}&{\frac {5}{36}}+{\frac {\sqrt {15}}{30}}&{\frac {2}{9}}+{\frac {\sqrt {15}}{15}}&{\frac {5}{36}}\\\hline &{\frac {5}{18}}&{\frac {4}{9}}&{\frac {5}{18}}\end{array}}$

偏微分方程式

偏微分方程式における選点法は、単独の方法より、あるスキームを実装するために必要とされる手法の一つに近い。例えば、放物型偏微分方程式の解を計算する場合、有限差分法より方程式の空間変数を離散化すると時間変数の常微分方程式となるので上述の選点法が適用できる。

現在よく使われているほとんどのスキームに選点法（必ずしも常微分方程式における選点法ではない）が使える。特に重要と見なされるスキームの中には、有限要素法（重み付き残差法も参照）、スペクトル法（選点法に基づいたスペクトル法は時々擬スペクトル法として知られている）がある（具体的な方法は、それぞれの記事を参照）。ここでは例として、フーリエ選点法と呼ばれるスペクトル法を簡単に紹介する。

フーリエ選点法は（理想的に）指数的収束速度を持ち、周期的境界条件を持つ方程式に対し特に効果的である。空間領域 [0, 2π] 上の（周期的境界条件付き）移流方程式

$u_{t}+a(x)u_{x}=0,\;u(x,0)=f(t)$

を考える。偶数の N に対し、等距離の点 x_j = jh （h = 2π/N）を選点として選ぶ。簡単のために、まず方程の時間変数を（前進）差分法より離散化し、次のようにする。

${\frac {U^{n+1}-U^{n}}{\Delta t}}=-a(x)u_{x}$

ここで、Uⁿ は時間 t_n での近似解である。したがって、正しく u_x を近似することで、次の時刻での近似解がわかるようになる。

既知の数値 v_j = Uⁿ(x_j) から空間上のグリッド関数 v = {v_j} を定義し、周期的に拡張する。そして a(x)u_x を次のように各選点で近似する。

$a(x_{j})u_{x}(x_{j})\approx a(x_{j})(Dv)_{j}$

ここで、D は スペクトル微分作用素 (spectral differentiation operator) という線型作用素であり、以下のように定義される。

$Dv={\mathcal {F}}_{h}^{-1}(i\xi {\mathcal {F}}_{h}(v))$

ここで、${\mathcal {F}}_{h}$ は半離散フーリエ変換 (semi-discrete Fourier transform) といい、以下のように定義される。

${\mathcal {F}}_{h}(v)=h\sum _{j=-\infty }^{\infty }e^{i\xi x_{j}}v_{j},\;\xi \in [-\pi /h,\pi /h]$

そして対応する逆変換は

$({\mathcal {F}}_{h}^{-1}({\hat {v}}))_{j}={\frac {1}{2\pi }}\int _{-\pi /h}^{\pi /h}e^{i\xi x_{j}}{\hat {v}}(\xi )d\xi ,\;j\in \mathbb {Z}$

である。上述の近似から、次の時刻 t_n+1 での近似解がわかる。

また、グリッド関数を空間全体に周期的に拡張する代わりに、そのまま離散フーリエ変換を使って近似することも可能である。この方法は高速フーリエ変換を活用できるため、計算速度が相応に上がる。離散フーリエ変換は以下のように定義される。

$({\mathcal {F}}(v))_{j}={\frac {1}{N}}\sum _{k=0}^{N}v_{k}e^{ijx_{k}},\;j=0,\pm 1,\ldots ,\pm N/2$

そして対応する逆変換は、

$({\mathcal {F}}^{-1}({\hat {v}}))_{j}=\sum _{k=-N/2}^{N/2}{\hat {v}}_{k}e^{ikx_{j}},\;j=0,1,\ldots ,N$

である。スペクトル微分作用素 D を使わずに、まず周波数領域の導関数を以下のように設定する。

${\hat {v}}_{j}={\begin{cases}ij({\mathcal {F}}(v))_{j}\;&j\neq \pm N/2\\0\;&j=\pm N/2\end{cases}}$

それから、u_x を以下のように近似できる。

$a(x_{j})u_{x}(x_{j})\approx a(x_{j})({\mathcal {F}}^{-1}({\hat {v}}))_{j}$

最後に、次の時刻での近似解を同じように計算する。

参考文献

• Ascher, Uri M.; Petzold, Linda R. (1998), Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, Philadelphia: Society for Industrial and Applied Mathematics, ISBN 978-0-89871-412-8 .
• Hairer, Ernst; Nørsett, Syvert Paul; Wanner, Gerhard (1993), Solving ordinary differential equations I: Nonstiff problems, Berlin, New York: Springer-Verlag, ISBN 978-3-540-56670-0 .
• Iserles, Arieh (2008), A First Course in the Numerical Analysis of Differential Equations (Second Edition), Cambridge University Press, ISBN 978-0-521-73490-5 .
• Trefethen, Lloyd N. (1996年). “Finite Difference and Spectral Methods”. 2017年1月2日閲覧。