LoginSignup
2
3

More than 3 years have passed since last update.

@Baruim

Gauss--Kronrod求積法

Last updated at Posted at 2020-06-04

$\phantom{}$$
\newcommand\nr{\notag\\}%タグなし改行用
\newcommand\ret{\notag\\&\qquad}%数式改行用
\newcommand\I{\mathrm{i}}
\newcommand\D{\mathrm{d}}
\newcommand\E{\mathrm{e}}
\newcommand\hc{\mathrm{h.c.}}
\newcommand\cc{\mathrm{c.c.}}
\newcommand\O[1]{\mathscr{O}\left(#1\right)}
\newcommand\abs[1]{{\left\rvert #1 \right\lvert}}
\newcommand\Res{\mathop{\mathrm{Res}}}
\newcommand\bra[1]{\mathinner{\langle{#1}|}}
\newcommand\ket[1]{\mathinner{|{#1}\rangle}}
\newcommand\braket[1]{\mathinner{\langle{#1}\rangle}}
\newcommand\Bra[1]{\left\langle#1\right|}
\newcommand\Ket[1]{\left|#1\right\rangle}
\newcommand\Braket[1]{\left\langle #1 \right\rangle}
\newcommand|{\middle|}%\Braket用
\DeclareMathOperator{\Log}{Log}
\DeclareMathOperator{\Arg}{Arg}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\Tr}{Tr}
\newcommand\dashint{\mathchoice
{\int\kern-10pt-}
{\int\kern-8.5pt-}
{\int\kern-6.1pt-}
{\int\kern-4.58pt-}}
\newcommand\set[1]{\left\{#1\right\}}
\newcommand\for[1]{\quad\mathrm{for}\quad #1}
\newcommand\LHS{\mathrm{(LHS)}}
\newcommand\RHS{\mathrm{(RHS)}}
$

Gauss求積法の解説記事の多さに比べて，
そのKronrod拡張1である表題の方法のそれはあまり見当たらないため，詳述することにする．

Gauss求積法 (Gauss quadrature) とは$N$点の関数評価によって誤差が$2N$次となる高精度数値積分法であって，次のようなものであった2 3 4 5 6．
すなわち，実数$x$を引数に持つ複素数値関数$f(x)$と，ある重み関数$w(x) > 0$の積の
区間$[a,b]$に関して積分は

\begin{align}
\int_a^b\D x w(x) f(x)
 = 
\sum_{j = 1}^N
w_j^\mathrm{G} f(x_j^\mathrm{G})
+
\frac{f^{(2N)}(\xi)}{(2N)!}h_N\quad (a < \xi < b),
 \tag{1}\label{eq:G}
\end{align}

と表示できる．ここで$\int_a^b\D x w(x)$は有限に定まるとする．
また，$w_j^\mathrm{G}$はChristoffel数である [($\ref{eq:wg}$)参照]．
$\set{x_j^\mathrm{G}}_{j = 1}^{N}$は

\begin{align}
\pi_N(x_j^\mathrm{G}) = 0\quad (j = 1,2,\ldots,N),
\end{align}

を満たす$x_j^\mathrm{G}\in[a,b]$なる数であり，Gauss点と呼ぶことにする．
$\pi_j(x)$は$w(x)$に関する$j$次の直交多項式であり，次式で定義される：

\begin{align}
\int_a^b\D x w (x)\pi_i(x)\pi_j(x)& = \delta_{ij} h_i, \tag{2}\label{eq:pi}
\end{align}

ここで

\begin{align}
h_i&: = \int_a^b\D x w (x)\pi_i(x)^2 > 0. \tag{3}\label{eq:h}
\end{align}

とおいた．直交性($\ref{eq:pi}$)だけでは振幅が定まらないため，本稿では付加条件として

\begin{align}
\pi_j^{(j)}(0) = j!,
\end{align}

を加えているが，表式を簡単にするための条件であって必須ではない．
これは多項式の最大ベキの係数を1とする意味を持っているだけである (monic polynomials)：

\begin{align}
\pi_j(x) = x^j+ ? x^{j-1}+? x^{j-2}+\cdots+?,
 \tag{4}\label{eq:aexp}
\end{align}

ここで「$?$」の部分にはある数が入る (記号の節約のため以後時々使う)．
ここまでに述べたGauss求積法は一般論であるが，
単にGauss求積法と言及した場合は$w(x) = 1$の特殊例を指すことが多い．
本稿では$w(x) = 1$の場合をGauss--Legendre求積法と呼んで区別することにする．

さて，次数が異なる直交多項式のゼロ点は0を除いて決して重ならないという特徴があるため，
Gauss求積法で次数$N$を上げると0以外のGauss点は互いに全て異なる．
これは$\set{f(x_j^\mathrm{G})}_{j = 1}^N$の再利用が出来ないことを意味し，
計算コストの点で問題があり，Gauss求積法が自動数値積分法として適さない理由となっていた．
これに対し，KronrodはGauss求積法の拡張版，すなわちGauss--Kronrod求積法 (以下，GK公式と呼ぶ) を考案した：

\begin{align}
\int_a^b\D x w(x) f(x)
 = 
\sum_{j = 1}^N
v_j^\mathrm{G} f(x_j^\mathrm{G})
+
\sum_{j = 1}^{N+1}
v_j^\mathrm{K} f(x_j^\mathrm{K})
+
R (f^{(3N+2)}),
 \tag{5}\label{eq:GK}
\end{align}

ここで$v_j^\mathrm{G}$はGauss点における重み，
$v_j^\mathrm{K}$は追加点$\set{x_j^\mathrm{K}}_{j = 1}^{N+1}$における重みである．
このように，追加の計算負荷$\set{f(x_j^\mathrm{K})}_{j = 1}^{N+1}$によって
計算精度が$R(f^{(2N)})\to R(f^{(3N+2)})$に増強される7．
もし同じ計算負荷でGuass求積法を用いたとすると，増強は$R(f^{(2N)})\to R(f^{(2N+2)})$に留まることから，
その効果は明白である．
ただし，これは予め$\set{ f(x_j^\mathrm{G})}_{j = 1}^N$が評価済みだった場合の話であり，
もしこれがなかった場合は当然$2N+1$個の関数評価の負荷が生じる．
この場合ははじめから$2N+1$次のGauss求積法を使えばよいのであって ($R(f^{(4N+2)})$の精度が出る)，
GK公式の出番はない．
そのため，GK公式は専らGauss求積法のお供 (誤差評価法) として利用されている．

GK公式を使うためには，予め

\begin{align}
\set{x_j^\mathrm{GK}}&: = \set{x_j^\mathrm{G}}\cup\set{x_j^\mathrm{K}},\\
w_j^\mathrm{G}& = \frac{h_{N-1}}{\pi_{N-1}(x_j^\mathrm{G})\pi_{N}'(x_j^\mathrm{G})}, \tag{6}\label{eq:wg}\\
v_j^\mathrm{G}
& = 
w_j^\mathrm{G}
\left(
1
+
\frac{h_N}{h_{N-1}}
\frac{\pi_{N-1}(x_j^\mathrm{G})}{E_{N+1}(x_j^\mathrm{G})}
\right)
,
 \tag{7}\label{eq:vg}
\\
v_j^\mathrm{K}& = \frac{h_N}{\pi_N(x_j^\mathrm{K})
E_{N+1}'(x_j^\mathrm{K})}
 \tag{8}\label{eq:vk}
,
\end{align}

を計算して表を作っておく (例えばここで入手可能)．
$2N+1$個の点集合$\set{x_j^\mathrm{GK}}$はKronrod点と呼ばれる．
$E_{N+1}(x)$は$w(x)$に関する$N+1$次のStieltjes (スティルチェス) 多項式である．

以上がGauss--Kronrod求積法の一般論である．
単にGauss--Kronrod求積法と呼んだ場合は$w(x) = 1$の場合を指すことが多い．
本稿では$w(x) = 1$の場合をGauss--Legendre-Kronrod求積法と呼んで区別することにする．

以下では上式の導出を行うが，
直交多項式の基本的性質は既知のものとしたので，
詳細は文献8を当たられたい．
(普段から量子論を使っている人ならば問題ない)

導出

文献を検索してみたところ，Monegatoの計算ノート9が大変簡潔にまとまっていると感じたのでこれに従うことにする．
以下はそれに色々と補足を追加あるいは部分的にカットしたものである．

最初に追加点$\set{x_j^\mathrm{K}}_{j = 1}^{N+1}$を定義しよう．
そのために，$N+1$次の多項式$E_{N+1}(x)$を次式で定義する：

\begin{align}
\int_a^b\D x w(x) \pi_N(x) E_{N+1}(x) x^k: = 0\quad (k = 0,1,2,\ldots,N),
 \tag{9}\label{eq:sp}
\end{align}

$E_{N+1}(x)$を$w(x)$に関する$N+1$次のStieltjes多項式と呼ぶ1011．
この表式だけでは振幅が定まらないため，本稿ではmonicであるとする：

\begin{align}
E_{N+1}^{(N+1)}(0) = (N+1)!.
\end{align}

このようにすると以下のようにして$E_{N+1}(x)$が完全に定まる．
まず，$E_{N+1}(x)$は$N+1$次の多項式だから$\set{\pi_j(x)}_{j = 1}^{N+1}$によって展開できる：

\begin{align}
E_{N+1}(x) = \sum_{j = 0}^{N+1} C_j \pi_j(x). \tag{10}\label{eq:exp}
\end{align}

これを($\ref{eq:sp}$)に代入すると

\begin{align}
0 = \sum_{j = 0}^{N+1} C_j \int_a^b\D x w(x) \pi_N(x)\pi_j(x) x^k\quad(k = 0,1,2\ldots,N),
 \tag{11}\label{eq:sp2}
\end{align}

となる．上式に$k = 0$を代入してみると，直交性($\ref{eq:pi}$)によって

\begin{align}
0 = C_N h_N,
\end{align}

となるので$C_N = 0$である．また，我々は多項式をmonicなものに選んでいるから，
$\pi_{N+1}(x) = x^{N+1}+\cdots$および$E_{N+1}(x) = x^{N+1}+\cdots$であって，
($\ref{eq:exp}$)より直ちに$C_{N+1} = 1$である．
結果，($\ref{eq:exp}$)は

\begin{align}
E_{N+1}(x) = \pi_{N+1}(x)+\sum_{j = 0}^{N-1} C_j \pi_j(x),
 \tag{12}\label{eq:s1}
\end{align}

である．残る未知数は$\set{C_j}_{j = 0}^{N-1}$の$N$個である．
一方，$C_N = 0$と$C_{N+1} = 1$を($\ref{eq:sp2}$)に代入すると，$N$本の連立1次方程式：

\begin{align}
\sum_{j = 0}^{N-1} C_j \int_a^b\D x w(x) \pi_N(x)\pi_j(x) x^k
 = 
-\int_a^b\D x w(x) \pi_N(x)\pi_{N+1}(x) x^k
\quad(k = 1,2\ldots,N)
,
 \tag{13}\label{eq:s2}
\end{align}

が得られるから，$\set{C_j}_{j = 0}^{N-1}$は全て確定し，すなわち$E_{N+1}(x)$が確定することがわかる．
このように導入したStieltjes多項式のゼロ点によって$\set{x_j^\mathrm{K}}_{j = 1}^{N+1}$を定義する：

\begin{align}
E_{N+1}(x_j^\mathrm{K}) = 0\quad (j = 1,2,\ldots,N+1).
\end{align}

ここで問題となるのは，$\set{x_j^\mathrm{K}}_{j = 1}^{N+1}$が全て実数になるかわからないということであるが，
本稿ではこの問題に立ち入らず，これらが全て実数となるような$w(x)$だけを考えているものとする．
また，$E_{N+1}(x)$の具体的な計算例は次節に回した．

ゼロ点がわかれば次のように因数分解できる：

\begin{align}
E_{N+1}(x) = \prod_{j = 1}^{N+1}(x-x_j^\mathrm{K}).
\end{align}

同様に

\begin{align}
\pi_N(x) = \prod_{j = 1}^{N+1}(x-x_j^\mathrm{G}),
\end{align}

である．最初に$v_j^\mathrm{K} $を求めよう．
そのために次の$2N$次の多項式を考える：

\begin{align}
f_k(x): = \frac{\pi_N(x)E_{N+1}(x)}{x-x_k^\mathrm{K}}
 = 
\pi_N(x)
\prod_{j = 1,j\neq k}^{N+1}(x-x_j^\mathrm{K})
 = 
\pi_N(x)
(x^N+?x^{N-1}+?x^{N-2}+\cdots+?x+?)
,
\end{align}

この場合，GK公式($\ref{eq:GK}$)の誤差項はゼロであるから，

\begin{align}
\int_a^b\D x w(x) f_k(x)
 = 
\sum_{j = 1}^N
v_j^\mathrm{G} f_k(x_j^\mathrm{G})
+
\sum_{j = 1}^{N+1}
v_j^\mathrm{K} f_k(x_j^\mathrm{K}),
 \tag{14}\label{eq:GK2}
\end{align}

は厳密に成立する．この両辺を比較して$v_j^\mathrm{K} $を求めよう．
まず左辺は

\begin{align}
\LHS
 = 
\int_a^b\D x w(x)
\pi_N(x)
(x^N+? x^{N-1}+\cdots)
,
\end{align}

である．上式で$A x^{N-1}$以下の次数は$\pi_N(x)$の直交性($\ref{eq:pi}$)により全て消える．
このことは，これらの項が$\set{\pi_j(x)}_{j = 0}^{N-1}$で展開可能であることを考えれば自明であろう．
よって

\begin{align}
\int_a^b\D x w(x) f_k(x)
 = 
\int_a^b\D x w(x)
\pi_N(x)
x^N
,
\end{align}

である．更に，($\ref{eq:aexp}$)から$\pi_N(x) = x^N+?x^{N-1}+\cdots$であるから，上式の$x^N$を
$\pi_N(x)$で置き換えても結果は変わらない ($?x^{N-1}+\cdots$を加えても再び直交性($\ref{eq:pi}$)で消える)．
よって($\ref{eq:h}$)を用いて

\begin{align}
\int_a^b\D x w(x) f_k(x)
 = 
\int_a^b\D x w(x)
\pi_N(x)^2
\equiv
h_N
,
 \tag{15}\label{eq:temp1}
\end{align}

を得る．今度は($\ref{eq:GK2}$)の右辺から，

\begin{align}
\RHS = 
\sum_{j = 1}^{N+1}
v_j^\mathrm{K} \pi_N(x_j^\mathrm{K})
\prod_{i = 1,i\neq k}^{N+1}(x_j^\mathrm{K}-x_i^\mathrm{K})
 = 
v_k^\mathrm{K} \pi_N(x_k^\mathrm{K})
\prod_{i = 1,i\neq k}^{N+1}(x_k^\mathrm{K}-x_i^\mathrm{K}),
\end{align}

を得る (Gauss項は$\pi_N(x_j^\mathrm{G})\equiv 0$から消える)．上式の制限付総乗は

\begin{align}
\prod_{i = 1,i\neq k}^{N+1}(x_k^\mathrm{K}-x_i^\mathrm{K})
 = 
E_{N+1}'(x_k^\mathrm{K}),
\end{align}

のように書けることに注意して，

\begin{align}
\sum_{j = 1}^N
v_j^\mathrm{G} f_k(x_j^\mathrm{G})
+
\sum_{j = 1}^{N+1}
v_j^\mathrm{K} f_k(x_j^\mathrm{K})
 = 
v_k^\mathrm{K} \pi_N(x_k^\mathrm{K})
E_{N+1}'(x_k^\mathrm{K})
,
\end{align}

が得られる．上式と($\ref{eq:temp1}$)より

\begin{align}
v_j^\mathrm{K} = \frac{h_N}{\pi_N(x_j^\mathrm{K})
E_{N+1}'(x_j^\mathrm{K})},
\end{align}

であることがわかった．

次にGauss点での重みの表式：$\set{v_j^\mathrm{G}}_{j = 1}^N$を求めよう．
再び$2N$次の多項式を改めて導入する：

\begin{align}
f_k(x): = \frac{\pi_N(x)E_{N+1}(x)}{x-x_k^\mathrm{G}}
 = 
E_{N+1}(x)
\prod_{j = 1,j\neq k}^{N}(x-x_j^\mathrm{G})
.
 \tag{16}\label{eq:f2}
\end{align}

これをGauss求積法($\ref{eq:G}$)に代入すると厳密な表式が得られる：

\begin{align}
\int_a^b\D x w(x) f_k(x)
 = 
\sum_{j = 1}^N w_j^\mathrm{G} f_k(x_j^\mathrm{G})
+
h_N
 = 
w_k^\mathrm{G}
E_{N+1}(x_k^\mathrm{G})
\prod_{i = 1,i\neq k}^{N}(x_k^\mathrm{G}-x_i^\mathrm{G})
+
h_N
 = 
w_k^\mathrm{G}
E_{N+1}(x_k^\mathrm{G})
\pi_N'(x_k^\mathrm{G})
+
h_N.
 \tag{17}\label{eq:G2}
\end{align}

一方，GK公式($\ref{eq:GK}$)も厳密であり，同様に

\begin{align}
\int_a^b\D x w(x) f_k(x)
 = 
v_k^\mathrm{G} 
E_{N+1}(x_k^\mathrm{G})
\pi_N'(x_k^\mathrm{G})
,
\end{align}

を得る．両者の比較から

\begin{align}
v_j^\mathrm{G}
 = 
w_j^\mathrm{G}
+
\frac{h_N}{E_{N+1}(x_j^\mathrm{G})
\pi_N'(x_j^\mathrm{G})},
\end{align}

が得られる．

Christoffel数$w_j^\mathrm{G}$も求めておこう．
次の$N-1$次の多項式：

\begin{align}
f_k(x): = \frac{\pi_N(x)}{x-x_k^\mathrm{G}},
\end{align}

をGauss求積法($\ref{eq:G}$)に代入し，これまでと同様に右辺を計算すると

\begin{align}
\int_a^b\D x w(x)\frac{\pi_N(x)}{x-x_k^\mathrm{G}}
 = 
w_k^\mathrm{G}\pi'_N(x_k^\mathrm{G}),
\end{align}

が得られる (厳密)．
よって

\begin{align}
w_j^\mathrm{G}
 = \frac{1}{\pi'_N(x_k^\mathrm{G})}
\int_a^b\D x w(x)\frac{\pi_N(x)}{x-x_j^\mathrm{G}},
\end{align}

となる．
分母に$x$が入っているのが厄介に見えるが，
直交多項式が一般に満たす関係式 (Christoffel--Darbouxの公式) を用いれば処理できる：

\begin{align}
\frac{\pi_N(x)}{x-x_k^\mathrm{G}}
 = 
\frac{h_{N-1}}{\pi_{N-1}(x_k^\mathrm{G})}
\sum_{i = 0}^{N-1}\frac{\pi_i(x)\pi_i(x_k^\mathrm{G})}{h_i}
. \tag{18}\label{eq:pd}
\end{align}

これにより，

\begin{align}
\int_a^b\D x w(x) f_j(x)
& = 
\frac{h_{N-1}}{\pi_{N-1}(x_j^\mathrm{G})}
\sum_{i = 0}^{N-1}\frac{\pi_i(x_j^\mathrm{G})}{h_i}
\int_a^b\D x w(x)
\pi_i(x)
 = 
\frac{h_{N-1}}{\pi_{N-1}(x_j^\mathrm{G})}
,
\end{align}

である．ここで，上式第2辺で現れた積分は$\pi_0(x) = 1$と直交性($\ref{eq:pi}$)を用いて

\begin{align}
\int_a^b\D x w(x)
\pi_i(x)
 = 
\int_a^b\D x w(x)
\pi_i(x)
\pi_0(x)
 = 
\delta_{i0}
h_0,
\end{align}

のように計算している．それゆえ

\begin{align}
w_j^\mathrm{G}
 = \frac{h_{N-1}}{\pi_{N-1}(x_j^\mathrm{G})\pi_{N}'(x_j^\mathrm{G})},
\end{align}

を得る．

Stieltjes多項式の計算

Gauss求積法では与えられた$w(x)$に関する直交多項式のセット$\set{\pi_j(x)}$が必要であった．
GK公式を使うためには，更に$E_{N+1}(x)$を求める必要がある．
ここでは定義通りの方法として，($\ref{eq:s2}$)を解いて係数列$\set{C_j}_{j = 0}^{N-1}$を求め，
($\ref{eq:s1}$)に代入することを考える．
再掲すると，

\begin{align}
E_{N+1}(x) = \pi_{N+1}(x)+\sum_{j = 0}^{N-1} C_j \pi_j(x),
 \tag{19}\label{eq:saikei}
\end{align}

の$N$個の展開係数：$\set{C_j}_{j = 0}^{N-1}$は
$N$本の連立1次方程式($\ref{eq:s2}$)：

\begin{align}
\sum_{j = 0}^{N-1} C_j \int_a^b\D x w(x) \pi_N(x)\pi_j(x) x^i
 = 
-\int_a^b\D x w(x) \pi_N(x)\pi_{N+1}(x) x^i
\quad(i = 1,2\ldots,N)
,
 \tag{20}\label{eq:s3}
\end{align}

を解いて得られるのであった．
左辺の添字の動く範囲を次のように調整しよう：

\begin{align}
\sum_{j = 1}^{N} C_{N-j} \int_a^b\D x w(x) \pi_N(x)\pi_{N-j}(x) x^i
 = 
-\int_a^b\D x w(x) \pi_N(x)\pi_{N+1}(x) x^i
\quad(i = 1,2\ldots,N)
,
\end{align}

ここで行列要素が

\begin{align}
A_{ij}: = \frac{1}{h_N}\int_a^b\D x w(x) \pi_N(x)\pi_{N-j}(x) x^i\quad (i,j = 1,2,\ldots,N),
 \tag{21}\label{eq:Aij}
\end{align}

で与えられる$N\times N$の行列を導入する．
同様に

\begin{align}
\tilde C& = 
\begin{pmatrix}
C_{N-1}\\
C_{N-2}\\
\vdots\\
C_{0}
\end{pmatrix},
\quad
\mu = 
\begin{pmatrix}
\mu_{1}\\
\mu_{2}\\
\vdots\\
\mu_{N}
\end{pmatrix}
,
\quad
\mu_i: = 
\frac{1}{h_N}
\int_a^b\D x w(x) \pi_N(x)\pi_{N+1}(x) x^{i}
\end{align}

と書こう ($h_N$で割った理由はすぐにわかる)．そうすると($\ref{eq:s3}$)は

\begin{align}
A\tilde C = -\mu,
 \tag{22}\label{eq:linear}
\end{align}

と表される．行列$A$は左三角行列であることに注意しよう．
なぜなら，$A_{ij}$の$\pi_{N-j}(x) x^i$の部分を$x$のベキで展開したときの最大次数は$x^{N+i-j}$であるが，
$i < j$を仮定すると$N+i-j < N$であるため，直交性($\ref{eq:pi}$)によって消えて

\begin{align}
A_{ij} = 0\for{i < j},
\end{align}

となるからである．また，対角成分は1である ($h_N$で割った理由)：

\begin{align}
A_{jj}
 = \frac{1}{h_N}\int_a^b\D x w(x) \pi_N(x)\pi_{N-j}(x) x^j
 = \frac{1}{h_N}\int_a^b\D x w(x) \pi_N(x)\pi_{N}(x)
 = 1
,
\end{align}

故に

\begin{align}
A = 
\begin{pmatrix}
1 & & & 0 \\
A_{21} & 1 & & \\
\vdots&\vdots &\ddots\\
A_{N1} & A_{N2}&\cdots& 1
\end{pmatrix}
,
\end{align}

という構造 (左三角行列) をしている．
($\ref{eq:linear}$)を陽に書けば

\begin{align}
\begin{pmatrix}
1 & & & 0 \\
A_{21} & 1 & & \\
\vdots&\vdots &\ddots\\
A_{N1} & A_{N2}&\cdots& 1 
\end{pmatrix}
\begin{pmatrix}
C_{N-1}\\
C_{N-2}\\
\vdots\\
C_{0}
\end{pmatrix}
 = 
-
\begin{pmatrix}
\mu_{1}\\
\mu_{2}\\
\vdots\\
\mu_{N}
\end{pmatrix},
 \tag{23}\label{eq:explicit}
\end{align}

となる．$\det A = 1$なので必ずStieltjes多項式が存在することもはっきりした．

更に楽することを考えよう12．
一般にmonicな直交多項式は次の形の3項漸化式を満たす：

\begin{align}
\pi_j(x) = (a_j+x)\pi_{j-1}(x)-\frac{h_{j-1}}{h_{j-2}}\pi_{j-2}(x)\quad(j = 2,3,\ldots),
 \tag{24}\label{eq:rec}
\end{align}

ここで$a_j$はある係数である．
上式で$j\to N-j$とおくと

\begin{align}
\pi_{N-j}(x) = (a_{N-j}+x)\pi_{N-j-1}(x)-\frac{h_{N-j-1}}{h_{N-j-2}}\pi_{N-j-2}(x),
\end{align}

となるが，これを($\ref{eq:Aij}$)の$\pi_{N-j}(x)$に代入して適当に添字を調整すると

\begin{align}
A_{i,j}
 = 
A_{i-1,j-1}
-
a_{N-j+1}A_{i-1,j}
+
\frac{h_{N-j}}{h_{N-j-1}}A_{i-1,j+1},
 \tag{25}\label{eq:recA}
\end{align}

を得る．上式から$A_{ij}$は上列にある3つのものから計算できることがわかる (図1)．
これにより，行列$A$は$\mu_i$の値と漸化式から全て決定されることがわかる (図2)．

image.png
図1. $A_{ij}$の依存関係．

image.png
図2. 漸化式によって$\mu_i$の値から行列$A$が決定されていく様子を表した模式図 ($N=5$)．
赤の数字は要素が計算されていく順序の一例を表している．

応用例 (Gauss--Legendre--Kronrod求積法)

Kronrodが示したもので，単にGauss--Kronrod積分公式と言えばこれを指す．
$w(x) = 1, a = -1, b = 1$に選んだものである．すなわち，

\begin{align}
\int_{-1}^1\D x f(x)
 = 
\sum_{j = 1}^N
v_j^\mathrm{G} f(x_j^\mathrm{G})
+
\sum_{j = 1}^{N+1}
v_j^\mathrm{K} f(x_j^\mathrm{K})
+
R (f^{(3N+2)}),
 \tag{28}\label{eq:GLK}
\end{align}

である．
$w(x) = 1$に関する直交多項式はLegendre多項式$P_i(x)$であり，通常の記法では直交性は

\begin{align}
\int_{-1}^1
\D x
P_i(x) P_j(x)
 = 
\frac{2}{2i+1}\delta_{ij}.
\end{align}

と表される．
しかし，我々は一貫してmonic多項式を仮定していたので，それに合わせて，

\begin{align}
\pi_i(x): = \frac{i!}{(2i-1)!!}P_i(x) = x^i+\cdots,
\end{align}

とおこう(ここの$k_n$で$P_i(x)を割る)．これよりまず，

\begin{align}
h_i
 = \int_{-1}^1\D x \pi_i(x)^2
 = 
\left(\frac{i!}{(2i-1)!!}\right)^2
\int_{-1}^1\D x P_i(x)^2
 = 
\left(\frac{i!}{(2i-1)!!}\right)^2
\frac{2}{2i+1}
\end{align}

が得られる．次にChristoffel数$w_j^\mathrm{G}$を求めよう．
Legendre多項式は次の漸化式を満たす：

\begin{align}
(1-x^2)P_N'(x) = -Nx P_N(x)+N P_{N-1}(x).
\end{align}

上式に$x = x_j^\mathrm{G}$を代入すると，

\begin{align}
(1-(x_j^\mathrm{G})^2)P_N'(x_j^\mathrm{G}) = N P_{N-1}(x_j^\mathrm{G}),
\end{align}

である．ここで$x_j^\mathrm{G}$の定義より$P_N(x_j^\mathrm{G}) = 0$となることに注意．
それゆえ

\begin{align}
\pi_N'(x_j^\mathrm{G})
 = 
\frac{N!}{(2N-1)!!}P_N'(x_j^\mathrm{G})
 = 
\frac{N!}{(2N-1)!!}
\frac{N}{1-(x_j^\mathrm{G})^2}
P_{N-1}(x_j^\mathrm{G}),
\end{align}

となる．よって

\begin{align}
w_j^\mathrm{G}
 = 
\frac{h_{N-1}}{\pi_{N-1}(x_j^\mathrm{G})\pi'_N(x_j^\mathrm{G})}
 = 
\frac{
\left(\frac{(N-1)!}{(2(N-1)-1)!!}\right)^2
\frac{2}{2(N-1)+1}
}{
\frac{(N-1)!}{[2(N-1)-1]!!}P_{N-1}(x_j^\mathrm{G})
\cdot
\frac{N!}{(2N-1)!!}
\frac{N}{1-(x_j^\mathrm{G})^2}
P_{N-1}(x_j^\mathrm{G})
}
 = 
\frac{
2[1-(x_j^\mathrm{G})^2]
}{
N^2
P_{N-1}(x_j^\mathrm{G})^2
}
 \tag{29}\label{eq:wjg2}
\end{align}

が得られる．従ってまずGauss--Legendre求積法が得られた：

\begin{align}
\int_{-1}^1\D x f(x)\simeq \sum_{j = 1}^N
\frac{
2[1-(x_j^\mathrm{G})^2]
}{
N^2
P_{N-1}(x_j^\mathrm{G})^2
}
f(x_j^\mathrm{G}).
\end{align}

$x_j^\mathrm{G}$は$P_N(x_j^\mathrm{G}) = 0$を数値的に解いて求める．

更にGauss--Legendre--Kronrod積分公式に進むには，
$w(x) = 1$に関するStieltjes多項式$E_{N+1}(x)$を知る必要がある (Mathematicaで得る方法)．
$A_{ij}$は($\ref{eq:Aij}$)より

\begin{align}
A_{ij}& = 
\frac{1}{h_N}
\frac{N!}{(2N-1)!!}
\frac{(N-j)!}{(2N-2j-1)!!}
\int_{-1}^1\D x P_N(x)x^iP_{N-j}(x)
\nr
& = 
\left(\frac{(2N-1)!!}{N!}\right)^2
\frac{2N+1}{2}
\frac{N!}{(2N-1)!!}
\frac{(N-j)!}{(2N-2j-1)!!}
\int_{-1}^1\D x P_N(x)x^iP_{N-j}(x)
\nr
& = 
\alpha_j
\int_{-1}^1\D x P_N(x)x^iP_{N-j}(x)
,
 \tag{30}\label{eq:alij}
\end{align}

から求められる．ここで

\begin{align}
\alpha_j: = \frac{(N-j)!(2N+1)!!}{2N!(2N-2j-1)!!},
\end{align}

とおいた．最初に$\mu_i = A_{i,-1}$を計算してみよう：

\begin{align}
\mu_i
& = 
-
\frac{2N+1}{2}
\int_{-1}^1\D x P_N(x)x^iP_{N+1}(x)
. \tag{31}\label{eq:mu}
\end{align}

Legendre多項式は偶関数または奇関数であり，かつ次数が1つ変化すると偶関数と奇関数が入れ替わる．
すなわち，$P_N(x)P_{N-1}(x)$は常に奇関数である．
従って上式で$i$が偶数のものは全てゼロである：

\begin{align}
\mu_i = 0\for{\mathrm{even}\ i}.
\end{align}

そこで我々は$i$が奇数のものだけを考えればよい．
次に漸化式によって三角行列$A$を作ろう．
Legendre多項式が満たす3項漸化式は

\begin{align}
nP_n(x) = (2n-1)xP_{n-1}(x)-(n-1)P_{n-2}(x),
\end{align}

であった．これより

\begin{align}
\pi_j(x) = x\pi_{j-1}-\frac{(j-1)^2}{(2j-1)(2j-3)}\pi_{j-2},
\end{align}

であって，($\ref{eq:rec}$)で$a_j = 0$とおいたものであることがわかる．
従って($\ref{eq:recA}$)から

\begin{align}
A_{i,j}
 = 
A_{i-1,j-1}
+
\frac{
(N-j)^2
}{(2N-2j-1)(2N-2j+1)}
A_{i-1,j+1},
\end{align}

を得る．

簡単のため，$N = 5$を考えよう．
$\mu$は($\ref{eq:mu}$)から直接計算し

\begin{align}
\mu = 
\begin{pmatrix}
A_{1,-1}\\
A_{2,-1}\\
A_{3,-1}\\
A_{4,-1}\\
A_{5,-1}
\end{pmatrix}
 = 
\begin{pmatrix}
\frac{36}{143}\\
0\\
\frac{136}{715}\\
0\\
\frac{8156}{51051}
\end{pmatrix}
,
\end{align}

となる．これを出発点として後は漸化式を繰り返し使用することで

\begin{align}
A = 
\begin{pmatrix}
1 & 0 & 0 & 0 & 0 \\
0 &1 & 0 & 0 & 0 \\
\frac{69}{91} &0 & 1 & 0 & 0 \\
0 & \frac{66}{65} & 0 & 1 & 0 \\
\frac{291}{453} &0 & \frac{50}{39} & 0 & 1 
\end{pmatrix}
,
\end{align}

を得ることができる．これより($\ref{eq:explicit}$)は

\begin{align}
\begin{pmatrix}
1 & 0 & 0 & 0 & 0 \\
0 &1 & 0 & 0 & 0 \\
\frac{69}{91} &0 & 1 & 0 & 0 \\
0 & \frac{66}{65} & 0 & 1 & 0 \\
\frac{291}{453} &0 & \frac{50}{39} & 0 & 1 
\end{pmatrix}
\begin{pmatrix}
C_4\\
C_3\\
C_2\\
C_1\\
C_0
\end{pmatrix}
 = 
-
\begin{pmatrix}
\frac{36}{143}\\
0\\
\frac{136}{715}\\
0\\
\frac{8156}{51051}
\end{pmatrix},
\end{align}

となって連立1次方程式が完成した．
上式を解くと

\begin{align}
\begin{pmatrix}
C_4\\
C_3\\
C_2\\
C_1\\
C_0
\end{pmatrix}
 = 
\begin{pmatrix}
-\frac{36}{143}\\
0\\
\frac{4}{5915}\\
0\\
\frac{496}{1307215}
\end{pmatrix},
\end{align}

となる．この結果を($\ref{eq:saikei}$)に代入すると$w(x) = 1$に関する6次のStieltjes多項式：

\begin{align}
E_{6}(x) = x^6-\frac{21}{13}x^4+\frac{567}{845}x^2-\frac{8043}{186745},
\end{align}

が得られる．$x_j^\mathrm{K}$は$E_6(x_j^\mathrm{K}) = 0$を数値的に解いて求める．

これでようやくGK公式の重みを計算する準備が整った．
($\ref{eq:wjg2}$)と($\ref{eq:vg}$)，($\ref{eq:vk}$)より

\begin{align}
w_j^\mathrm{G}
& = 
\left.
\frac{
128(1-x^2)
}{
25(35x^4-30x^2+3)^2
}\right|_{x = x_j^\mathrm{G}}
,
\\
v_j^\mathrm{G}
& = 
w_j^\mathrm{G}
\left(
1
+
\frac{25}{99}
\frac{x^4-\frac{6}{7}x^2+\frac{3}{35}}{
x^6-\frac{21}{13}x^4+\frac{567}{845}x^2-\frac{8043}{186745}
}
\right)_{x = x_j^\mathrm{G}}
,\\
v_j^\mathrm{K}& = 
\left.
\frac{128}{43659}
\frac{1}{\left(x^5-\frac{10}{9}x^3+\frac{5}{21}x\right)
\left(
6x^5-\frac{84}{13}x^3+\frac{1134}{845}x
\right)}
\right|_{x = x_j^\mathrm{K}}
,
\end{align}

となる．

まずGauss点$\set{x_j^\mathrm{G}}_{j = 1}^5$は

\begin{align}
x_1^\mathrm{G}& = -0.9061798459386640,\\
x_2^\mathrm{G}& = -0.5384693101056831,\\
x_3^\mathrm{G}& = 0,\\
x_4^\mathrm{G}& = +0.5384693101056831,\\
x_5^\mathrm{G}& = +0.9061798459386640,
\end{align}

となった．これは$P_5(x) = 0$を数値計算によって解いたものである．
Mathematicaなら1行である：

Mathematica
x /. NSolve[LegendreP[5, x], x, WorkingPrecision -> 15] // MatrixForm

実行結果：

image.png

次に追加点$\set{x_j^\mathrm{K}}_{j = 1}^5$は

\begin{align}
x_1^\mathrm{K}& = -0.9840853600948425,\\
x_2^\mathrm{K}& = -0.7541667265708492,\\
x_3^\mathrm{K}& = -0.2796304131617832,\\
x_4^\mathrm{K}& = +0.2796304131617832,\\
x_5^\mathrm{K}& = +0.7541667265708492,\\
x_6^\mathrm{K}& = +0.9840853600948425,
\end{align}

となった．これは$E_6(x) = 0$を数値計算によって解いたものである．
Mathematicaなら1行である：

Mathematica
x /. NSolve[x^6 - 21/13 x^4 + 567/845 x^2 - 8043/186745, x, 
   WorkingPrecision -> 16] // MatrixForm

実行結果：

image.png

これらを合わせた11点をKronrod点と呼ぶ．
Kronrod点は次のように入れ子構造になっていることがわかる (一般論であるが，証明は省いた9)：

\begin{align}
x_1^\mathrm{K}& = -0.9840853600948425,\\
x_1^\mathrm{G}& = -0.9061798459386640,\\
x_2^\mathrm{K}& = -0.7541667265708492,\\
x_2^\mathrm{G}& = -0.5384693101056831,\\
x_3^\mathrm{K}& = -0.2796304131617832,\\
x_3^\mathrm{G}& = 0,\\
x_4^\mathrm{K}& = +0.2796304131617832,\\
x_4^\mathrm{G}& = +0.5384693101056831,\\
x_5^\mathrm{K}& = +0.7541667265708492,\\
x_5^\mathrm{G}& = +0.9061798459386640,\\
x_6^\mathrm{K}& = +0.9840853600948425,
\end{align}

Kronrod点に対応する重みを評価すると

\begin{align}
v_1^\mathrm{K}& = 0.042582036751082,\\
v_1^\mathrm{G}& = 0.11523331662247,\\
v_2^\mathrm{K}& = 0.18680079655649,\\
v_2^\mathrm{G}& = 0.24104033922865,\\
v_3^\mathrm{K}& = 0.272849801912559,\\
v_3^\mathrm{G}& = 0.2829874178574912,\\
v_4^\mathrm{K}& = 0.272849801912559,\\
v_4^\mathrm{G}& = 0.24104033922865,\\
v_5^\mathrm{K}& = 0.18680079655649,\\
v_5^\mathrm{G}& = 0.11523331662247,\\
v_6^\mathrm{K}& = 0.042582036751082,
\end{align}

となる．以上の結果はたしかに11点GK公式の表と一致している．

Gauss--Legendre求積法の重みの計算結果も載せておく：

\begin{align}
w_1^\mathrm{G}& = 0.23692688505619,\\
w_2^\mathrm{G}& = 0.47862867049937,\\
w_3^\mathrm{G}& = 0.5688888888888889,\\
w_4^\mathrm{G}& = 0.47862867049937,\\
w_5^\mathrm{G}& = 0.23692688505619,
\end{align}

この結果も5点Gauss求積法の表と一致している．

  1. A. S. Kronrod, Nodes and Weights of Quadrature Formulas, English transl. from Russian, Consultants Bureau (New York, 1965). ↩

  2. DLMF. ↩

  3. F. B. Hildebarnd, Introduction to Numerical Analysis: Second Edition, McGraw-Hill (New York, 1794). ↩

  4. 青本和彦，「直交多項式入門」，数学書房 (2013)，p. 153. ↩

  5. 加古富志雄，「数値解析」，第六回， ↩

  6. 千葉豪，Gauss--Legendre求積公式の導出． ↩

  7. S. E. Notaris, Electron. Trans. Numer. Anal. 45, 371 (2016). ↩

  8. 青本和彦，「直交多項式入門」，数学書房 (2013)，p. 33; 伏見康治，赤井逸，「直交関数系　増補版」，共立 (2011)，p. 58; T. H. Koornwinder, arXiv:1303.2825 等． ↩

  9. G. Monegato, Math. Comp. 30, 812 (1976). ↩

  10. G. Monegato, SIAM Review 24, 137 (1982). ↩

  11. Encyclopedia of Mathematics. ↩

  12. C. Blumstein and J. C. Wheeler, Phys. Rev. B 8 , 1764 (1973). ↩

2
3
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
3