数値積分法

$\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\for[1]{\quad\mathrm{for}\quad #1}
\newcommand\set[1]{\left\{#1\right\}}
$　__定積分の結果が欲しくて悩んでいるほぼ全ての人達__は数値計算を専門としてはいない．
従ってそれが必要なら，「既知の卓越した積分ストラテジを適切に実装し，極限まで磨かれた最強の既存ルーチン」をcallしろというのは尤もな言い分である (例えばquadpack)．

ところが，__定積分の結果が欲しくて悩んでいるほぼ全ての人達__は，
定積分の実行が最終目的ではなく，全体の計算の途中に現れるほんの一部に過ぎない．
そのため，既存の卓越したルーチンを用いるためだけにプログラムの構造を書き換えようものなら，
__プログラム全体としてパフォーマンスがガタ落ちする__ということがしばしばある．

その結果，学生時代に使った数値解析の教科書や古代より伝わる Numerical Recipe を紐解き，
__こそこそと禁じ手とされる自作積分ルーチン__を組み始めることになるのである．

この記事は，

(1) あまり凝ったことをする気も余裕もない
(2) でも台形公式では心許ない
(3) 手っ取り早くそこそこの精度と速度で数値積分できれば満足
(4) 被積分関数の振る舞いについて，ある程度感覚を持っている
(5) 宗教上の理由で自分で書かなければならない

といった境遇の人を対象としたものである．

以下では1次元の積分に限る．
2次元，3次元は1次元積分の逐次積分として処理する (Fubiniの定理)．
(これは__最も__愚劣な方法である．しかし，我々は数値計算の専門家ではないから許される）

この記事で考慮していない内容で特に注意を要するもの：

・振動積分 (Fourier変換の類)
・Cauchy主値積分

これらの積分にここで紹介する方法をナイーブに適用すると__100%不幸になるので注意．これらは特別な配慮が必要な積分である__．
また，モンテカルロ法による求積法は扱っていない．

定積分の数値計算の一般論

実数$x$に対する複素数値関数$f(x)$の定積分の数値積分は全て

\int_a^b \D x f(x)=\sum_{j=1}^N w_j f(x_j)+R_N,

の形に表される ($a<b$)．
複素積分を実行したいときは適当な実パラメータを導入して左辺の形に帰着させる．
$w_j,x_j$は各積分法によって与えられる．$R_N$は誤差であり，上式を等号で結ぶために必要なものである．
このように書かれるとき，$N$は$f(x)$を評価する回数を表す．
我々はこの回数をなるべく減らし，しかも高精度を得たいというわけである．
$\set{x_j}$と$\set{w_j}$の表があれば，あとは$f(x_j)$を計算し，$w_j$と掛け算して和を取るだけである．
英語では前者をabscissa(s), 後者をweight(s)という．
　2次元の積分は単純に

\int_{a_1}^{b_1}\D x_1\int_{a_1}^{b_2} \D x_2 f(x_1,x_2)\simeq\sum_{j=1}^{N_1}\sum_{k=1}^{N_2} w_1^j w_2^k f(x_1^j,x_2^k),

とする．$f(x_1,x_2)$ 1点の評価に時間がかかる場合は並列計算することも検討する．
上式で$x_1$積分を並列化する場合，

Fortran

    !$OMP parallel do schedule(dynamic) default(none)&
    !$OMP &shared(integrands,weights,x1,x2,w1,w2) private(j,k)
    do j=1,size(x1)
        do k=1,size(x2)
            integrands(j+(k-1)*size(x1)) = function_f(x1(j),x2(k))
            weights   (j+(k-1)*size(x1)) = w1(j)*w2(k)
        end do
    end do
    !$OMP end parallel do

    result = dot(integrands,weights) ! BLAS95

と書いてdoループを並列化すればよい．

定積分の数値計算法としては，ひとまず台形公式とGauss型積分公式だけ知っておけばよい．
それ以外の数値積分法は無視してよい (必要が生じたときに再検討する)，
ただし，特殊関数の数値計算法は非常によく研究されているので，それで表示できる場合はそちらを優先する¹．
以下によく使う数値積分法をリストする (後々追加していくかもしれない)．

台形公式 (trapezoidal formula)

基本公式

数値積分の最も基本となる公式である．

\int_a^b \D x f(x)\simeq \frac{b-a}{N}\left[\frac{f(a)+f(b)}{2}+\sum_{j=1}^{N-1} f\left(a+\frac{b-a}{N}j\right)\right]\equiv\sum_{j=0}^{N} w_j f(x_j).

ここで

\begin{align}
x_j&=a+\frac{b-a}{N}j,\\
w_j&=
\begin{cases}
\frac{b-a}{2N},\for{j=0\ \mathrm{or}\ N}\\
\frac{b-a}{N},\quad\mathrm{for\ others}
\end{cases},
\end{align}

となる．$f(a)=f(b)$の場合はより簡単化され，

I_N=\frac{b-a}{N}\sum_{j=0}^{N-1} f\left(a+\frac{b-a}{N}j\right)\for{f(a)=f(b)},

となる．特に，周期関数をその周期で積分する場合は特別な精度向上が期待できる^{2 3 4}．
表式から明らかなように，台形公式は公式というよりRiemann和そのものである．
よって，Riemann積分可能な問題であれば必ず収束する．
収束しない場合は，式が間違っているか，被積分関数が積分範囲において精度の限界を超えた変化をしている可能性を疑おう．
例えば，変数の無次元化やエネルギの原点の移動が適切に行われているか確認しよう．
これらは数値積分に取り組む以前の問題で，数値計算の常識である．

この手法の特徴は，

(1) 公式が最も単純である
(2) $\set{w_j}$と$\set{x_j}$が最も単純である
(3) 刻み幅$(b-a)/N$を半分にしてゆくイタレーションを作ることにより，$f(x_j)$を__完全に__再利用できる

という点にある．これらにより，以下の特長につながる．
(1) バグの入る余地が少なくなる．
そのため，未知の数値積分に挑むときはまずこの方法で評価してみるべきである．
より複雑で高精度な積分法に進むときのリファレンスにもなる．

(2) どんなに多くの点が必要になっても，$\set{w_j}$と$\set{x_j}$は一瞬で用意できる．
全体的に被積分関数の変化が激しい or 激しい領域が一部あるがその場所が先験的に分からないといった事情から，
積分領域全体の刻み幅を小さくとらねばならないときに有利である．

(3) の特徴は$f(x)$1点の評価に高いコストを有する問題で極めて重要なことである．
(そしてそのような場合は応用上多い)
また，誤差評価が簡単という利点もある (凝った積分法ではこれが面倒である)．

最も単純な公式とは言え，特殊な条件下では異常とも言える高精度を叩き出す公式であるため侮れない⁵．

###誤差評価法
現在の刻み幅を半分にしたものと比較する．
すなわち，現在の$N$の値を$n$，数値積分値を$I$とするとき，次回の数値積分値$I'$は$N=2n$の値から評価する：

I'=\frac{I}{2}+\frac{b-a}{2n}\sum_{j=1}^n f\left(a+\frac{b-a}{n}\left(j-\frac{1}{2}\right)\right).

もし$\abs{I-I'}$が所望の精度に達していなければ，
$n\Leftarrow 2n$，$I\Leftarrow I'$と更新した後，上式を再度評価する．
これを繰り返す．

注意として，比較対象を$N=n+1$のものに選んではいけない．
なぜなら，$\set{x_j=a+j(b-a)/n}$と$\set{x_j'=a+j(b-a)/(n+1)}$を比べるとわかるように，
これらの間には端点以外に共通の座標がなく，せっかく評価した$\set{f(x_j)|j=0,1,\ldots,n}$がほぼ全て無駄になるからである．これにより，まるまる$\set{f(x_j')|j=1,\ldots,n}$の計算負荷が生じる．
これは$I'$の計算負荷と全く同じである割に，精度はほとんど上がらない．
$I'$の表式では前回の$I$を再利用しているため，$N=2n$で達成される精度でありながら，
追加の計算負荷は$n$個の$f(x)$の評価で済んでいるのである．
要するに$N$を制御変数としてdoループを回してはいけない．

二重指数関数型積分公式 (double exponential formula)

先述した通り，特殊条件下では台形公式が異常な高精度を叩き出すことがある．
DE (double exponential) 積分公式は，この特殊条件を人工的に作り出す台形公式の最終進化形態である⁶．

以下の公式は考案者本人の著書から引用^{7 8}したものである．

###有限区間

\int_a^b\D x f(x)=\sum_{j=-\infty}^{\infty} w_j^\mathrm{DE} f(x_j^\mathrm{DE}),

と書いたとき，座標$x_j^\mathrm{DE}$と重み$w_j^\mathrm{DE}$はDE型変換後の台形公式の刻み幅$h$の関数であって，

\begin{align}
x_j^\mathrm{DE}&=\frac{b-a}{2}\tanh\left(\frac{\pi}{2}\sinh(jh)\right)+\frac{b+a}{2},\\
w_j^\mathrm{DE}&=\pi h\frac{b-a}{4}\frac{\cosh(jh)}{\cosh^2\left(\frac{\pi}{2}\sinh(jh)\right)},
\end{align}

である．

###半無限区間

\int_0^\infty\D x f(x)=\sum_{j=-\infty}^{\infty} w_j^\mathrm{DE} f(x_j^\mathrm{DE}),

パターン1:

\begin{align}
x_j^\mathrm{DE}&=\E^{2\sinh(jh)},\\
w_j^\mathrm{DE}&=2h\cosh(jh)x_j^\mathrm{DE}.
\end{align}

パターン2:

\begin{align}
x_j^\mathrm{DE}&=\E^{jh-\E^{-jh}},\\
w_j^\mathrm{DE}&=h(1+\E^{-jh})x_j^\mathrm{DE}.
\end{align}

パターン3:

\begin{align}
x_j^\mathrm{DE}&=\E^{\frac{jh}{2}-\E^{-jh}},\\
w_j^\mathrm{DE}&=h\left(\frac{1}{2}+\E^{-jh}\right)x_j^\mathrm{DE}.
\end{align}

###全無限区間

\int_{-\infty}^\infty\D x f(x)=\sum_{j=-\infty}^{\infty} w_j^\mathrm{DE} f(x_j^\mathrm{DE}),

において

\begin{align}
x_j^\mathrm{DE}&=\sinh\left(\frac{\pi}{2}\sinh(jh)\right),\\
w_j^\mathrm{DE}&=h\frac{\pi}{2}\cosh(jh)\cosh\left(\frac{\pi}{2}\sinh(jh)\right).
\end{align}

いずれも和の範囲が無限大になっているが，
もちろん収束判定によって有限で打ち切る^{8 9}．

台形公式なので基本的には台形公式の利点をそのまま引き継いでいる．
変数変換 (DE変換) には幾つもバリエーションがあるため
完全に使いこなすにはノウハウが必要であり，奥が深い．
また，過激な関数を使うので，
オーバーフローとアンダーフローに気をつけなければならない^{10 11}．
とりあえず上記のパターンで試してみるのが良いと思う．

周期関数をその周期で積分する場合にも最適らしい¹²．

Gauss求積法 (Gaussian quadrature)

ある与えられた重み関数 $w(x)>0$ に関する直交多項式のゼロ点を$\set{x_j}$に選択し，$f(x)$をLagrange補間する方法であり，

\int_a^b\D x w(x)f(x)=\sum_{k=1}^N w_k f(x_k)+ R_N,

の形を持つ．ここで$\set{x_k}$と$\set{w_k}$は

\int_a^b \D x w(x) p_i(x)p_j(x)\propto \delta_{ij}\quad (i,j=1,2,\ldots),

を満たす直交多項式$\set{p_i(x)}$によって

\begin{align}
p_N(x_k)&=0,\\
w_k&=\int_a^b \D x w(x)\frac{p_N(x)}{(x-x_k)p_N'(x_k)},
\end{align}

で与えられる．誤差は

R_N\propto \frac{\D^{2N} f(\xi)}{\D\xi^{2N}} \quad \mathrm{where}\quad a<\xi<b,

であって，$f(x)$が$x$の$2N-1$次以下の多項式のときは厳密値を与えることがわかる．
この公式は，

(1) 精度が高い．
(2) $\set{x_j}$と$\set{x_j}$が簡単に与えられない場合が多い．
(3) 台形公式のように$f(x)$の値を再利用できない．

という特徴を持つ．
上に書いた通り，$N$点の$f(x)$の評価から誤差が$f^{(2N)}(\xi)$のオーダの結果を得ることができ，コスパが高い．
一方，$\set{x_j}$は多項式のゼロ点で与えられるため，それ自体を求めるために数値計算が必要なことが多い
(ただし，予め$\set{x_j}$と$\set{w_j}$を用意しておけば良い話ではある)．
$f(x)$の変化が激しく，要求される$N$に対して$\set{x_j}$と$\set{w_j}$の表を用意することが困難な場合，
計算パラメータを導入して積分範囲を分割する必要があり，プログラムの複雑化を招く．
何より，このような余分なパラメータの出現はノウハウの要求を予感させ，精神衛生上よくない．
また，$N$を変えると基本的には関数値の再利用ができない．
なお，後述のGauss--Legendre求積法において「小さい$N$＋領域分割」よりも「大きい$N$＋領域分割なし」の方が高効率であるとの報告もある¹³．

ネガティブなことばかり書いたが，この精度の高さは魅力的であり，非常によく用いられる方法である．

様々な$w(x)$に対して公式が存在するが，
この記事ではAbramowitz and Stegunから2例引用した (Gauss--Legendre, Gauss--Chebyshev)．
先に誤差評価法について述べておく．

##誤差評価法
Gauss求積法の枠組みでは，$N$を変えると$f(0)$以外の$f(x)$の値は再利用できない．
この欠点を改善するために，Kronrod¹⁴は$N$次のGauss点：$\set{x_j|j=1,2,\ldots,N}$に新たな点：$\set{x_j'|j=1,2,\ldots,N+1}$を加えた$2N+1$個の点列$\set{x_j^\mathrm{GK}}=\set{x_j}\cup\set{x'_j}$ (Kronrod点と呼ぶことにする) を作り，
$2N+1$次の求積法を考案した (参考)．この求積法をGauss--Kronrod求積法と呼ぶ．

ここで，$\set{x_j'}$は$w(x)$に関するStieltjes多項式$E_{N+1}(x)$のゼロ点で与えられる：

E_{N+1}(x_j')=0\quad(j=1,2,\ldots,N+1),

$E_{N+1}(x)$は具体的には

\int_a^b\D x w(x) p_N(x) E_{N+1}(x) x^k =0\quad(k=0,1,\ldots,N),

によって定義される ($E_{N+1}(x)=\sum_{i=0}^{N+1} c_i p_i(x)$と展開して$c_i$を求める．ただし，$c_{N+1}$は確定せず，振幅の自由度として残る)．
このように，Kronrod点は$w(x)$に関する$N$次の直交多項式のゼロ点 (Gauss点) と上式から成る．
したがって$2N+1$次のKronrod点と$2N+1$次のGauss点は一致せず，$N$次のGauss求積法よりは精度で優るが，同次数のGauss求積法より精度で劣る．
このため，$2N+1$次のGauss--Kronrod求積法は専ら$N$次のGauss求積法の二番手として誤差評価に使われるにとどまっている ($2N+1$次のGauss--Kronrod求積法を使うために$2N+1$回関数を評価するなら初めから$2N+1$次のGauss求積法を使えばよく，全く利点がない)．

なお，Gauss--Kronrod求積法の重み $\set{w_j^\mathrm{GK}}$ は$f(x)$をKronrod点でLagrange補間し，それに $w(x)$ をかけて積分すれば求まる：

\begin{align}
f(x)&\simeq \sum_{j=1}^{2N+1} g_j(x) f(x_j^\mathrm{GK}),\\
\int_a^b \D x w(x) f(x)&=
\sum_{j=1}^{2N+1} w_j^\mathrm{GK} f(x_j^\mathrm{GK})\quad\text{(Gauss--Kronrod)},\\
w_j^\mathrm{GK}&:=\int_a^b\D x w(x) g_j(x).
\end{align}

書いていたら長くなったので別記事にまとめた．
原理が気になる人は参照されたい．

後にPattersonはこの手法を一般化し，$f(x)$を再利用し続けることが可能な公式を考案している (Gauss--Patterson求積法)．
原論文¹⁵には，1例として3個のGauss--Legendre点から出発して$f(x)$の値を再利用し続ける表が掲載されている．この点列はスタートとなるGauss--Legendre点の次数によって異なり(7次のGauss--Legendre--Patterson点と7次のGauss--Legendre点が異なるのだから当然である)，それぞれFamilyを作る．
また，Fortranサブルーチンも作成されている．
ただし，そこにも記述されているように，
完全な再利用性を備えることと引き換えに精度が低下しているため手放しで喜べるような方法とは言えない．
精度低下については同次数のGauss--Legendre--Patterson点とGauss--Legendre点が異なっていることからも納得されよう．

Gauss--Legendre--Kronrod求積法の表はこちら．
なお，Gauss--Chebyshev積分公式タイプ ($\set{w_j}$が一定) のKronrod拡張は存在しないが¹⁶，
ほぼChebyshevタイプのKronrod拡張は存在し，よく知られている^{16 17}．

##Gauss--Legendre求積法

$w(x)=1$，$a=-1$，$b=1$に対するGauss型求積法である：

\int_{-1}^1\D x f(x)=\sum_{j=1}^N w_j f(x_j)+\frac{2^{2N+1}(N!)^4}{(2N+1)[(2N)!]^3}f^{(2N)}(\xi)\quad(-1<\xi<1),

ここで$x_j$は

P_N(x_j)=0,\quad (j=1,2,\ldots,N),

を解いて求める．$P_N(x)$は$N$次のLegendre多項式である．$w_j$は

w_j=\frac{2}{(1-x_j^2)P'_N(x_j)^2}

で与えられる．
具体的な数値はここで得られる．
自分で計算する場合は，例えば_Mathematica_なら

20次Gauss求積法

Clear[n, x, xs, wp, sol, xs, ws];
n = 20;(*ノード数*)
wp = 100;(*作業精度*)
xs = x /. 
  NSolve[LegendreP[n, x] == 0, x, Reals, WorkingPrecision -> wp];
ws = 2/LegendreP[n, 1, xs]^2;
DATAPOINTS = Table[{xs[[i]], ws[[i]]}, {i, Length@xs, 1, -1}];
DATAPOINTS // MatrixForm (*左列が{x_j}，右列が{w_j}*)

とすれば表示される．
ここではLegendre多項式の導関数はLegendre陪多項式で表現している．

Gauss--Chebyshev求積法

$w(x)=1/\sqrt{1-x^2}$，$a=-1$，$b=1$に対するGauss型求積法である：

\int_{-1}^1\D x \frac{f(x)}{\sqrt{1-x^2}}=\sum_{j=1}^N w_j f(x_j)+\frac{\pi}{(2N)!2^{2N-1}}f^{(2N)}(\xi)\quad(-1<\xi<1),

ここで

\begin{align}
x_j&=\cos\left(\frac{2j-1}{2N}\pi\right),\\
w_j&=\frac{\pi}{N},
\end{align}

である．$\set{x_j}$と$\set{w_j}$が簡単なため，使いやすい．
これは$1/\sqrt{1-x^2}$に関する直交多項式である第1種Chebyshev多項式のゼロ点が簡単に求められるためである．

Gauss--Chebyshev--Kronrod求積法

表式が簡単なのでおいておく¹⁸：

\int_{-1}^1\D x \frac{f(x)}{\sqrt{1-x^2}}\simeq\frac{\pi}{2N}
\left[
\frac{f(-1)}{2}
+
\sum_{j=1}^{2N-1}
f\left(
\cos\frac{j\pi}{2N}
\right)
+
\frac{f(1)}{2}
\right]
=
\frac{1}{2}\sum_{j=1}^N w_jf(x_j)
+\frac{\pi}{2}
\left[
\frac{f(-1)}{2}
+
\sum_{j=1}^{N-1}
f\left(
\cos \frac{j\pi}{N}
\right)
+
\frac{f(1)}{2}
\right]
\for{N\ge 2}.

上式第3辺の第2項が追加項である．

1. 例えば，第1種完全楕円積分は数値積分を実行するのでなく，Landen変換に基づく漸化式から計算すべきである．他にも，Bessel関数のような特殊関数は安定な漸化式がよく知られており，それを使うべきである．もっとも，この記事では特殊関数で表現できるような簡単な問題は対象としていない． ↩

2. 伊里正夫，藤野和建，「数値計算の常識」，共立 (2005)，p. 47. ↩

3. 藤原宏志，「球面上の高精度積分則とその非適切問題への応用」，第63回理論応用力学講演会，東工大 (2014年9月). ↩

4. DLMF3.5.4 ↩

5. 森正武，「FORTRAN77 数値計算プログラミング(増補版)」，岩波（2007），p. 169. ↩

6. 森正武，数学 50, 248 (1998). ↩

7. 森正武，「FORTRAN77 数値計算プログラミング(増補版)」，岩波（2007），pp. 170--171. ↩

8. 森正武，「FORTRAN77 数値計算プログラミング(増補版)」，岩波（2007），pp. 375--377. ↩ ↩²

9. 森正武，応用数理 8, 323 (1998). ↩

10. 渡辺二太，核融合研究 63, 397 (1990). ↩

11. 森正武，「FORTRAN77 数値計算プログラミング(増補版)」，岩波（2007），p. 173. ↩

12. 森正武，数理解析研究所講究録 1040, 143 (1998). ↩

13. 渡辺二太，大内昭司，核融合研究 64, 397 (1990). ↩

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

15. T. N. L. Patterson, Math. Comput. 22, 847 (1968). ↩

16. S. E. Notaris, Math. Comput. 58, 745 (1992). ↩ ↩²

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

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