# 数値積分基礎のキソ
###### tags: `数値解析`

$%2020/04/24
\newcommand\setN[0]{\mathbb{N}}
\newcommand\setZ[0]{\mathbb{Z}}
\newcommand\setQ[0]{\mathbb{Q}}
\newcommand\setR[0]{\mathbb{R}}
\newcommand\setC[0]{\mathbb{C}}
\newcommand\setH[0]{\mathbb{H}}
\newcommand\setE[0]{\mathbb{E}}
\newcommand\pare[1]{{(#1)}}
\newcommand\Pare[1]{\left(#1\right)}
\newcommand\curl[1]{\{#1\}}
\newcommand\Curl[1]{\left\{#1\right\}}
\newcommand\squa[1]{[#1]}
\newcommand\Squa[1]{\left[#1\right]}
\newcommand\abs[1]{\lvert#1\rvert}
\newcommand\Abs[1]{\left\lvert#1\right\rvert}
\newcommand\norm[1]{\|#1\|}
\newcommand\Norm[1]{\left\|#1\right\|}
\newcommand\floor[1]{\lfloor#1\rfloor}
\newcommand\Floor[1]{\left\lfloor#1\right\rfloor}
\newcommand\ceil[1]{\lceil#1\rceil}
\newcommand\Ceil[1]{\left\lceil#1\right\rceil}
\newcommand\angl[1]{\langle#1\rangle}
\newcommand\Angl[1]{\left\langle#1\right\rangle}
\newcommand\conj[1]{\overline{#1}}
\newcommand\transpose[1]{\,{\vphantom{#1}}^t\!#1}
\newcommand\sfrac[2]{#1/#2}
\newcommand\od[2]{\frac{d#1}{d#2}}
\newcommand\pd[2]{\frac{\partial#1}{\partial#2}}
\newcommand\sod[2]{\sfrac{d#1}{d#2}}
\newcommand\spd[2]{\sfrac{\partial#1}{\partial#2}}
\newcommand\tensor[2]{\boldsymbol{\mathrm{T}}^#1_#2}
\newcommand\bm[1]{\boldsymbol{#1}}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator{\Re}{Re}
\DeclareMathOperator{\Im}{Im}
\DeclareMathOperator{\supp}{supp}
\DeclareMathOperator{\sgn}{sgn}
\DeclareMathOperator{\Diff}{Diff}
\DeclareMathOperator{\Sym}{Sym}
\DeclareMathOperator{\spn}{span}
\DeclareMathOperator{\grad}{grad}
\DeclareMathOperator{\rot}{rot}
\DeclareMathOperator{\div}{div}
\newcommand\tensor[2]{\mathbf{T}^{#1}_{#2}}$

## 数値積分とは
一般に, 任意の関数の積分を厳密に(解析的に)得ることは困難である.
例えば[楕円積分](https://ja.wikipedia.org/wiki/%E6%A5%95%E5%86%86%E7%A9%8D%E5%88%86)はその代表例であるし

\begin{align}
\int_{1/2}^{1}\frac{\sin(\sin(x))}{x^3+\tan(x)}dx
\end{align}

とか絶対できねーだろって思う訳である.

しかしの目的が数値解析であれば、解析的な積分結果よりも「おおよその値」さえ知れたら十分という状況は十分にある.
これが数値積分である.

## 数値積分の方法
1次元の閉区間上の数値積分を考えよう.
$f$が連続関数であれば

\begin{align}
I(f)
&=\int_a^b f(x)dx \\
&= \lim_{n\to\infty} \sum_{i=1}^{n} \frac{b-a}{n}f\Pare{a+\frac{b-a}{n}i} \\
&= \lim_{n\to\infty} \sum_{i=0}^{n-1} \frac{b-a}{n}f\Pare{a+\frac{b-a}{n}i} \\
\end{align}

が成り立つ(短冊状に区切るRiemann積分を思い出せば良い)から, これを適当な$n$を選んで有限和で打ち切ればおおよその値が得られそうである.

ここでは短冊の左端や右端の値を参照して区分定数関数で近似したが, 短冊の中点を参照した

\begin{align}
I(f)
&=\int_a^b f(x)dx \\
&= \lim_{n\to\infty} \sum_{i=1}^{n} \frac{b-a}{n}f\Pare{a+\frac{b-a}{n}(i-1/2)}
\end{align}

も考えられる.
さらに後述のような放物線で近似する方法(Simpson則)も考えられる.

いずれにしても, 数値積分は有限回の計算しか使えないから, 有限個の点列$\curl{x_i}$と重みと呼ばれる実数列$\curl{w_i}$を用意して

\begin{align}
I(f)=
\int_a^b f(x)dx
\approx\sum_{i=1}^{n} w_i f(x_i)
=I_n(f)
\end{align}

のような形で計算する必要がある.

## 数値積分の評価
数値積分は前述のような区分定数関数による近似から多項式近似, 直交多項式を使う方法など様々な手法があり, それらによって$x_i$や$w_i$の選び方が異なる.

では, これらの手法の優劣をどのようにして調べれば良いだろうか.
その評価の一つが多項式に関する精度評価である.

$\tilde{I}(f)$で積分$I(f)$を近似した数値積分を表すこととする.
このとき

\begin{align}
I(x^q)-\tilde{I}(x^q) &= 0 & (0 \le q \le p) \\
I(x^{p+1})-\tilde{I}(x^{p+1}) &\ne 0
\end{align}

が成立するとき, 数値積分$\tilde{I}$は$p$次の精度を持つという.
要するに, $p$次までの多項式なら厳密に計算できるよ！という訳だ.

## Newton-Cotes法
異なる$p+1$点$\curl{x_i}$を$\setR$上に取り, その各点に任意に実数$\curl{\alpha_i}$を対応させれば, $f(x_i)=\alpha_i$を充たすような$p$次多項式が一意的に存在する.
詳しくは[Lagrange多項式](https://en.wikipedia.org/wiki/Lagrange_polynomial)を参照のこと.

積分区間$[a,b]$を$p$等分すれば$p+1$点$\curl{x_i}$が得られるから,
これによって重み$\curl{w_i}$も決定できる.
これが[Newton-Cotes法](https://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas)である.
$n=2$が[台形則](https://en.wikipedia.org/wiki/Trapezoidal_rule), $n=3$が[Simpson則](https://en.wikipedia.org/wiki/Simpson%27s_rule)であり, Newton-Cotes法はこれらを一般化したものとなっている.

これまでの議論から, Newton-Cotes法は$n$点で**少なくとも**$n-1$次の精度を持つ数値積分と言える.
(実は$n$が奇数であれば$n$次の精度, 偶数であれば$n-1$次の精度である.)

## 複合Newton-Cotes法
* 被積分関数が多項式で近似しにくい
* 被積分関数が区分多項式

このような場合は, 積分区間を分割してからそれぞれの区間に対してNewton-Cotes法を使う方が都合が良い場合がある.

## 直交多項式を使う方法
Newton-Cotes法を使う方法では, $n$点に対して高々$n$次の精度しか得られなかったが, 別の方法を使えば$n$点に対して$2n-1$次の精度を得ることができる.
これが直交多項式を使う方法であり, [Gauss求積](https://en.wikipedia.org/wiki/Gaussian_quadrature)とよばれる.

さらに$n$点公式のうち, $2n-1$次を超える精度を持つものは存在しないことが知られている.
よって直交多項式を使う方法は最高の精度を持つと言える.
理論の詳細は参考文献を参照されたい.

Juliaで数値計算する際には[FastGaussQuadrature.jl](https://github.com/JuliaApproximation/FastGaussQuadrature.jl)という便利なパッケージがあるのでそちらを利用すると良い.

## 複合Gauss求積
直交多項式に関しても, 区間を分割してから積分する方が良い場合が多い.
B-spline基底関数を扱う場合などは特にそうである.

## 発展
### 多次元
矩形の領域$[a_1,b_1]\times\cdots\times[a_d,b_d]$であれば, 単に各次元の方向に1次元の積分を適用すれば良い.

### もっと多い次元
上記の多次元の手法では$n^d$個の点の参照が必要になるが, $d$がべらぼうに大きい場合は参照点が指数関数的に増えるので良い方法とは言えない.
この場合はMonte-Carlo法が使われる
https://en.wikipedia.org/wiki/Monte_Carlo_method#Integration

## 参考文献
[数値解析入門](https://www.amazon.co.jp/dp/4339060712)
やすい(2000えん！)なので買いましょう